From 9ba3c28ae91b440bb5d5350bf4e63db470bb7b34 Mon Sep 17 00:00:00 2001 From: Nicolas Tessore Date: Tue, 1 Nov 2022 10:19:35 +0000 Subject: [PATCH] initial commit --- .github/dependabot.yml | 6 + .github/workflows/cibuildwheel.yml | 61 +++ .gitignore | 3 + LICENSE.md | 28 + MANIFEST.in | 1 + README.md | 63 +++ pyproject.toml | 4 + python/chealpix.c | 740 +++++++++++++++++++++++++++ python/healpix/__init__.py | 178 +++++++ python/healpix/test/__init__.py | 11 + python/healpix/test/conftest.py | 59 +++ python/healpix/test/test_chealpix.py | 105 ++++ python/healpix/test/test_healpix.py | 75 +++ setup.cfg | 31 ++ setup.py | 15 + src/healpix.c | 409 +++++++++++++++ src/healpix.h | 168 ++++++ tools/memory_overhead.py | 116 +++++ tools/nside_precision.py | 18 + tools/randang_distribution.py | 54 ++ 20 files changed, 2145 insertions(+) create mode 100644 .github/dependabot.yml create mode 100644 .github/workflows/cibuildwheel.yml create mode 100644 .gitignore create mode 100644 LICENSE.md create mode 100644 MANIFEST.in create mode 100644 README.md create mode 100644 pyproject.toml create mode 100644 python/chealpix.c create mode 100644 python/healpix/__init__.py create mode 100644 python/healpix/test/__init__.py create mode 100644 python/healpix/test/conftest.py create mode 100644 python/healpix/test/test_chealpix.py create mode 100644 python/healpix/test/test_healpix.py create mode 100644 setup.cfg create mode 100644 setup.py create mode 100644 src/healpix.c create mode 100644 src/healpix.h create mode 100644 tools/memory_overhead.py create mode 100644 tools/nside_precision.py create mode 100644 tools/randang_distribution.py diff --git a/.github/dependabot.yml b/.github/dependabot.yml new file mode 100644 index 0000000..5ace460 --- /dev/null +++ b/.github/dependabot.yml @@ -0,0 +1,6 @@ +version: 2 +updates: + - package-ecosystem: "github-actions" + directory: "/" + schedule: + interval: "weekly" diff --git a/.github/workflows/cibuildwheel.yml b/.github/workflows/cibuildwheel.yml new file mode 100644 index 0000000..7c42e48 --- /dev/null +++ b/.github/workflows/cibuildwheel.yml @@ -0,0 +1,61 @@ +name: Build and upload to PyPI + +on: + push: + branches: + - main + pull_request: + branches: + - main + release: + types: + - published + +jobs: + build_wheels: + name: Build wheels on ${{ matrix.os }} + runs-on: ${{ matrix.os }} + strategy: + matrix: + os: [ubuntu-20.04, windows-2019, macos-11] + + steps: + - uses: actions/checkout@v3 + + - name: Build wheels + uses: pypa/cibuildwheel@v2.11.2 + env: + CIBW_SKIP: 'pp{37,38,39}-{macosx_x86_64,win_amd64}' + + - uses: actions/upload-artifact@v3 + with: + path: ./wheelhouse/*.whl + + build_sdist: + name: Build source distribution + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v3 + + - name: Build sdist + run: pipx run build --sdist + + - uses: actions/upload-artifact@v3 + with: + path: dist/*.tar.gz + + upload_pypi: + name: Upload to PyPI + needs: [build_wheels, build_sdist] + runs-on: ubuntu-latest + if: github.event_name == 'release' && github.event.action == 'published' + steps: + - uses: actions/download-artifact@v3 + with: + name: artifact + path: dist + + - uses: pypa/gh-action-pypi-publish@v1.5.1 + with: + user: __token__ + password: ${{ secrets.pypi_password }} diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..5c187c2 --- /dev/null +++ b/.gitignore @@ -0,0 +1,3 @@ +.*.swp +python/healpix.c +__pycache__ diff --git a/LICENSE.md b/LICENSE.md new file mode 100644 index 0000000..2696da7 --- /dev/null +++ b/LICENSE.md @@ -0,0 +1,28 @@ +Copyright (c) 2022, Nicolas Tessore, based on the healpix_bare library + +Copyright (c) 1997-2019, Krzysztof M. Gorski, Eric Hivon, Martin Reinecke, Benjamin D. Wandelt, Anthony J. Banday, Matthias Bartelmann, Reza Ansari & Kenneth M. Ganga + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + +1. Redistributions of source code must retain the above copyright notice, this + list of conditions and the following disclaimer. + +2. Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + +3. Neither the name of the copyright holder nor the names of its + contributors may be used to endorse or promote products derived from + this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE +FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/MANIFEST.in b/MANIFEST.in new file mode 100644 index 0000000..57f820d --- /dev/null +++ b/MANIFEST.in @@ -0,0 +1 @@ +include src/healpix.h diff --git a/README.md b/README.md new file mode 100644 index 0000000..196fb71 --- /dev/null +++ b/README.md @@ -0,0 +1,63 @@ +healpix +======= + +**Python and C package for HEALPix discretisation of the sphere** + +This package implements a lean set of routines for working with the HEALPix +discretisation of the sphere. It supports NSIDE parameters up to 2^29. + +The C library is based on the *healpix_bare* library, which was released under +the 3-clause BSD license. + +If you are using this code in your research, please consider citing the +original paper in your publications: + +* [Gorski et al., 2005, ApJ, 622, 759][Gorski+2005] + +[Gorski+2005]: http://adsabs.harvard.edu/abs/2005ApJ...622..759G + + +Python +------ + +The Python package provides functions that deal with the discretisation of the +sphere itself. It is not meant to be a replacement of *healpy*, and does not +contain things such as functions for the visualisation of maps, or spherical +harmonic transforms. + +The Python package consists of two modules: + +* The low-level `chealpix` module, which is a native C extension, and + efficiently vectorises the C library functions over arbitrary numpy array + inputs (including scalars, of course). +* The high-level `healpix` module, which contains a more streamlined interface, + and additional functionality: + + * Random point picking in HEALPix pixels. + +The high-level functions in the `healpix` module can be used more or less +interchangeably with functions from the *healpy* package. However, in some +cases, compatibility is sacrificed for consistency. + +The Python package requires only *numpy*, and can be installed using pip: + + pip install healpix + +The vectorised C functions carefully avoid the creation of temporary arrays, +and therefore have minimal memory overhead: + +```py +>>> import numpy as np, healpix, tracemalloc +>>> +>>> # random vectors with 1G of memory per component (less than a NSIDE=4K map) +>>> x, y, z = np.random.randn(3, 125_000_000) +>>> +>>> tracemalloc.start() +>>> +>>> lon, lat = healpix.vec2ang(x, y, z, lonlat=True) +>>> +>>> tracemalloc.get_traced_memory() +(2000010342, 2000013889) # current, peak +>>> +>>> # no memory overhead: only the 2G output arrays were used +``` diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..c22c505 --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,4 @@ +[build-system] +requires = ["setuptools", "wheel", "oldest-supported-numpy"] + +build-backend = "setuptools.build_meta" diff --git a/python/chealpix.c b/python/chealpix.c new file mode 100644 index 0000000..d7d48ab --- /dev/null +++ b/python/chealpix.c @@ -0,0 +1,740 @@ +// author: Nicolas Tessore +// license: BSD-3-Clause + +#define PY_SSIZE_T_CLEAN +#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION +#include +#include +#include + +#include "healpix.h" + +// maximum number of operands for vectorized operations +#ifndef VEC_MAXOP +#define VEC_MAXOP 6 +#endif + + +static void setang(t_ang ang, double* theta, double* phi) { + *theta = ang.theta; + *phi = ang.phi; +} + + +static void setvec(t_vec vec, double* x, double* y, double* z) { + *x = vec.x; + *y = vec.y; + *z = vec.z; +} + + +typedef void (*vecfunc)(void*, npy_intp, void**); + + +PyObject* vectorize(vecfunc func, void* args, npy_intp nin, npy_intp nout, + PyObject** op, int* types) { + npy_intp i, nop; + NpyIter* iter; + PyArrayObject* op_[VEC_MAXOP] = {0}; + npy_uint32 flags; + npy_uint32 op_flags[VEC_MAXOP]; + PyArray_Descr* op_dtypes[VEC_MAXOP] = {0}; + PyArrayObject** out; + PyObject* ret = NULL; + + nop = nin + nout; + + if (nop > VEC_MAXOP) { + PyErr_Format(PyExc_RuntimeError, + "chealpix internal error: increase VEC_MAXOP to %llu", + nop); + return NULL; + } + + for (i = 0; i < nin; ++i) { + op_[i] = (PyArrayObject*)PyArray_FromAny(op[i], NULL, 0, 0, 0, NULL); + op_flags[i] = NPY_ITER_READONLY | NPY_ITER_NBO | NPY_ITER_ALIGNED | + NPY_ITER_CONTIG; + op_dtypes[i] = PyArray_DescrFromType(types[i]); + if (!op_[i] || !op_dtypes[i]) + goto fail; + } + + for (; i < nop; ++i) { + if (op[i] && op[i] != Py_None) { + op_[i] = (PyArrayObject*)PyArray_FromAny(op[i], NULL, 0, 0, 0, + NULL); + if (!op_[i]) + goto fail; + } + op_flags[i] = NPY_ITER_WRITEONLY | NPY_ITER_ALLOCATE | NPY_ITER_NBO | + NPY_ITER_ALIGNED | NPY_ITER_CONTIG; + op_dtypes[i] = PyArray_DescrFromType(types[i]); + if (!op_dtypes[i]) + goto fail; + } + + flags = NPY_ITER_EXTERNAL_LOOP | NPY_ITER_BUFFERED | NPY_ITER_GROWINNER | + NPY_ITER_ZEROSIZE_OK; + + iter = NpyIter_MultiNew(nop, op_, flags, NPY_KEEPORDER, NPY_SAFE_CASTING, + op_flags, op_dtypes); + if (!iter) + goto fail; + + if (NpyIter_GetIterSize(iter) != 0) { + NpyIter_IterNextFunc* iternext = NpyIter_GetIterNext(iter, NULL); + char** dataptr = NpyIter_GetDataPtrArray(iter); + npy_intp* sizeptr = NpyIter_GetInnerLoopSizePtr(iter); + if (!iternext) + goto iterfail; + do { + func(args, *sizeptr, (void**)dataptr); + } while (iternext(iter)); + } + + out = &NpyIter_GetOperandArray(iter)[nin]; + if (nout == 0) { + ret = Py_None; + Py_INCREF(ret); + } else if (nout == 1) { + Py_INCREF(out[0]); + ret = PyArray_Return(out[0]); + } else { + ret = PyTuple_New(nout); + if (!ret) + goto iterfail; + for (i = 0; i < nout; ++i) { + Py_INCREF(out[i]); + PyTuple_SET_ITEM(ret, i, PyArray_Return(out[i])); + } + } + + if (NpyIter_Deallocate(iter) != NPY_SUCCEED) + goto fail; + + for (i = 0; i < nop; ++i) { + Py_XDECREF(op_[i]); + Py_DECREF(op_dtypes[i]); + } + + return ret; + +iterfail: + NpyIter_Deallocate(iter); + +fail: + for (i = 0; i < nop; ++i) { + Py_XDECREF(op_[i]); + Py_XDECREF(op_dtypes[i]); + } + Py_XDECREF(ret); + return NULL; +} + + +static void vang2vec(void* args, npy_intp size, void** data) { + double* theta = data[0], *phi = data[1]; + double* x = data[2], *y = data[3], *z = data[4]; + for (npy_intp i = 0; i < size; ++i) + setvec(ang2vec((t_ang){theta[i], phi[i]}), &x[i], &y[i], &z[i]); +} + + +PyDoc_STRVAR(cang2vec_doc, +"ang2vec(nside, theta, phi, x=None, y=None, z=None, /)\n" +"--\n" +"\n"); + + +static PyObject* cang2vec(PyObject* self, PyObject* args) { + PyObject* op[] = {NULL, NULL, NULL, NULL, NULL}; + int types[] = {NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE}; + + if (!PyArg_ParseTuple(args, "OO|OOO:ang2vec", + &op[0], &op[1], &op[2], &op[3], &op[4])) { + return NULL; + } + + return vectorize(vang2vec, NULL, 2, 3, op, types); +} + + +static void vvec2ang(void* args, npy_intp size, void** data) { + double* x = data[0], *y = data[1], *z = data[2]; + double* theta = data[3], *phi = data[4]; + for (npy_intp i = 0; i < size; ++i) + setang(vec2ang((t_vec){x[i], y[i], z[i]}), &theta[i], &phi[i]); +} + + +PyDoc_STRVAR(cvec2ang_doc, +"vec2ang(nside, x, y, z, theta=None, phi=None, /)\n" +"--\n" +"\n"); + + +static PyObject* cvec2ang(PyObject* self, PyObject* args) { + PyObject* op[] = {NULL, NULL, NULL, NULL, NULL}; + int types[] = {NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE}; + + if (!PyArg_ParseTuple(args, "OOO|OO:vec2ang", + &op[0], &op[1], &op[2], &op[3], &op[4])) { + return NULL; + } + + return vectorize(vvec2ang, NULL, 3, 2, op, types); +} + + +static void vang2nest(void* args, npy_intp size, void** data) { + int64_t nside = *(int64_t*)args; + double* theta = data[0], *phi = data[1]; + int64_t* ipix = data[2]; + for (npy_intp i = 0; i < size; ++i) + ipix[i] = ang2nest(nside, (t_ang){theta[i], phi[i]}); +} + + +PyDoc_STRVAR(cang2nest_doc, +"ang2nest(nside, theta, phi, ipix=None, /)\n" +"--\n" +"\n"); + + +static PyObject* cang2nest(PyObject* self, PyObject* args) { + Py_ssize_t nside; + PyObject* op[] = {NULL, NULL, NULL}; + int types[] = {NPY_DOUBLE, NPY_DOUBLE, NPY_INT64}; + + if (!PyArg_ParseTuple(args, "nOO|O:ang2nest", &nside, + &op[0], &op[1], &op[2])) { + return NULL; + } + + return vectorize(vang2nest, &nside, 2, 1, op, types); +} + + +static void vang2ring(void* args, npy_intp size, void** data) { + int64_t nside = *(int64_t*)args; + double* theta = data[0]; + double* phi = data[1]; + int64_t* ipix = data[2]; + for (npy_intp i = 0; i < size; ++i) + ipix[i] = ang2ring(nside, (t_ang){theta[i], phi[i]}); +} + + +PyDoc_STRVAR(cang2ring_doc, +"ang2ring(nside, theta, phi, ipix=None, /)\n" +"--\n" +"\n"); + + +static PyObject* cang2ring(PyObject* self, PyObject* args) { + Py_ssize_t nside; + PyObject* op[] = {NULL, NULL, NULL}; + int types[] = {NPY_DOUBLE, NPY_DOUBLE, NPY_INT64}; + + if (!PyArg_ParseTuple(args, "nOO|O:ang2ring", &nside, + &op[0], &op[1], &op[2])) { + return NULL; + } + + return vectorize(vang2ring, &nside, 2, 1, op, types); +} + + +static void vnest2ang(void* args, npy_intp size, void** data) { + int64_t nside = *(int64_t*)args; + int64_t* ipix = data[0]; + double* theta = data[1], *phi = data[2]; + for (npy_intp i = 0; i < size; ++i) + setang(nest2ang(nside, ipix[i]), &theta[i], &phi[i]); +} + + +PyDoc_STRVAR(cnest2ang_doc, +"nest2ang(nside, ipix, theta=None, phi=None, /)\n" +"--\n" +"\n"); + + +static PyObject* cnest2ang(PyObject* self, PyObject* args) { + Py_ssize_t nside; + PyObject* op[] = {NULL, NULL, NULL}; + int types[] = {NPY_INT64, NPY_DOUBLE, NPY_DOUBLE}; + + if (!PyArg_ParseTuple(args, "nO|OO:nest2ang", &nside, + &op[0], &op[1], &op[2])) { + return NULL; + } + + return vectorize(vnest2ang, &nside, 1, 2, op, types); +} + + +static void vring2ang(void* args, npy_intp size, void** data) { + int64_t nside = *(int64_t*)args; + int64_t* ipix = data[0]; + double* theta = data[1], *phi = data[2]; + for (npy_intp i = 0; i < size; ++i) + setang(ring2ang(nside, ipix[i]), &theta[i], &phi[i]); +} + + +PyDoc_STRVAR(cring2ang_doc, +"ring2ang(nside, ipix, theta=None, phi=None, /)\n" +"--\n" +"\n"); + + +static PyObject* cring2ang(PyObject* self, PyObject* args) { + Py_ssize_t nside; + PyObject* op[] = {NULL, NULL, NULL}; + int types[] = {NPY_INT64, NPY_DOUBLE, NPY_DOUBLE}; + + if (!PyArg_ParseTuple(args, "nO|OO:ring2ang", &nside, + &op[0], &op[1], &op[2])) + return NULL; + + return vectorize(vring2ang, &nside, 1, 2, op, types); +} + + +static void vvec2nest(void* args, npy_intp size, void** data) { + int64_t nside = *(int64_t*)args; + double* x = data[0], *y = data[1], *z = data[2]; + int64_t* ipix = data[3]; + for (npy_intp i = 0; i < size; ++i) + ipix[i] = vec2nest(nside, (t_vec){x[i], y[i], z[i]}); +} + + +PyDoc_STRVAR(cvec2nest_doc, +"vec2nest(nside, x, y, z, ipix=None, /)\n" +"--\n" +"\n"); + + +static PyObject* cvec2nest(PyObject* self, PyObject* args) { + Py_ssize_t nside; + PyObject* op[] = {NULL, NULL, NULL, NULL}; + int types[] = {NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE, NPY_INT64}; + + if (!PyArg_ParseTuple(args, "nOOO|O:vec2nest", &nside, + &op[0], &op[1], &op[2], &op[3])) { + return NULL; + } + + return vectorize(vvec2nest, &nside, 3, 1, op, types); +} + + +static void vvec2ring(void* args, npy_intp size, void** data) { + int64_t nside = *(int64_t*)args; + double* x = data[0], *y = data[1], *z = data[2]; + int64_t* ipix = data[3]; + for (npy_intp i = 0; i < size; ++i) + ipix[i] = vec2ring(nside, (t_vec){x[i], y[i], z[i]}); +} + + +PyDoc_STRVAR(cvec2ring_doc, +"vec2ring(nside, x, y, z, ipix=None, /)\n" +"--\n" +"\n"); + + +static PyObject* cvec2ring(PyObject* self, PyObject* args) { + Py_ssize_t nside; + PyObject* op[] = {NULL, NULL, NULL, NULL}; + int types[] = {NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE, NPY_INT64}; + + if (!PyArg_ParseTuple(args, "nOOO|O:vec2ring", &nside, + &op[0], &op[1], &op[2], &op[3])) { + return NULL; + } + + return vectorize(vvec2ring, &nside, 3, 1, op, types); +} + + +static void vnest2vec(void* args, npy_intp size, void** data) { + int64_t nside = *(int64_t*)args; + int64_t* ipix = data[0]; + double* x = data[1], *y = data[2], *z = data[3]; + for (npy_intp i = 0; i < size; ++i) + setvec(nest2vec(nside, ipix[i]), &x[i], &y[i], &z[i]); +} + + +PyDoc_STRVAR(cnest2vec_doc, +"nest2vec(nside, ipix, x=None, y=None, z=None, /)\n" +"--\n" +"\n"); + + +static PyObject* cnest2vec(PyObject* self, PyObject* args) { + Py_ssize_t nside; + PyObject* op[] = {NULL, NULL, NULL, NULL}; + int types[] = {NPY_INT64, NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE}; + + if (!PyArg_ParseTuple(args, "nO|OOO:nest2vec", &nside, + &op[0], &op[1], &op[2], &op[3])) { + return NULL; + } + + return vectorize(vnest2vec, &nside, 1, 3, op, types); +} + + +static void vring2vec(void* args, npy_intp size, void** data) { + int64_t nside = *(int64_t*)args; + int64_t* ipix = data[0]; + double* x = data[1], *y = data[2], *z = data[3]; + for (npy_intp i = 0; i < size; ++i) + setvec(ring2vec(nside, ipix[i]), &x[i], &y[i], &z[i]); +} + + +PyDoc_STRVAR(cring2vec_doc, +"ring2vec(nside, ipix, x=None, y=None, z=None, /)\n" +"--\n" +"\n"); + + +static PyObject* cring2vec(PyObject* self, PyObject* args) { + Py_ssize_t nside; + PyObject* op[] = {NULL, NULL, NULL, NULL}; + int types[] = {NPY_INT64, NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE}; + + if (!PyArg_ParseTuple(args, "nO|OOO:ring2vec", &nside, + &op[0], &op[1], &op[2], &op[3])) { + return NULL; + } + + return vectorize(vring2vec, &nside, 1, 3, op, types); +} + + +static void vang2nest_uv(void* args, npy_intp size, void** data) { + int64_t nside = *(int64_t*)args; + double* theta = data[0], *phi = data[1]; + int64_t* ipix = data[2]; + double* u = data[3], *v = data[4]; + for (npy_intp i = 0; i < size; ++i) + ipix[i] = ang2nest_uv(nside, (t_ang){theta[i], phi[i]}, &u[i], &v[i]); +} + + +PyDoc_STRVAR(cang2nest_uv_doc, +"ang2nest_uv(nside, theta, phi, ipix=None, u=None, v=None, /)\n" +"--\n" +"\n"); + + +static PyObject* cang2nest_uv(PyObject* self, PyObject* args) { + Py_ssize_t nside; + PyObject* op[] = {NULL, NULL, NULL, NULL, NULL}; + int types[] = {NPY_DOUBLE, NPY_DOUBLE, NPY_INT64, NPY_DOUBLE, NPY_DOUBLE}; + + if (!PyArg_ParseTuple(args, "nOO|OOO:ang2nest_uv", &nside, + &op[0], &op[1], &op[2], &op[3], &op[4])) { + return NULL; + } + + return vectorize(vang2nest_uv, &nside, 2, 3, op, types); +} + + +static void vang2ring_uv(void* args, npy_intp size, void** data) { + int64_t nside = *(int64_t*)args; + double* theta = data[0], *phi = data[1]; + int64_t* ipix = data[2]; + double* u = data[3], *v = data[4]; + for (npy_intp i = 0; i < size; ++i) + ipix[i] = ang2ring_uv(nside, (t_ang){theta[i], phi[i]}, &u[i], &v[i]); +} + + +PyDoc_STRVAR(cang2ring_uv_doc, +"ang2ring_uv(nside, theta, phi, ipix=None, u=None, v=None, /)\n" +"--\n" +"\n"); + + +static PyObject* cang2ring_uv(PyObject* self, PyObject* args) { + Py_ssize_t nside; + PyObject* op[] = {NULL, NULL, NULL, NULL, NULL}; + int types[] = {NPY_DOUBLE, NPY_DOUBLE, NPY_INT64, NPY_DOUBLE, NPY_DOUBLE}; + + if (!PyArg_ParseTuple(args, "nOO|OOO:ang2ring_uv", &nside, + &op[0], &op[1], &op[2], &op[3], &op[4])) { + return NULL; + } + + return vectorize(vang2ring_uv, &nside, 2, 3, op, types); +} + + +static void vnest2ang_uv(void* args, npy_intp size, void** data) { + int64_t nside = *(int64_t*)args; + int64_t* ipix = data[0]; + double* u = data[1], *v = data[2]; + double* theta = data[3], *phi = data[4]; + for (npy_intp i = 0; i < size; ++i) + setang(nest2ang_uv(nside, ipix[i], u[i], v[i]), &theta[i], &phi[i]); +} + + +PyDoc_STRVAR(cnest2ang_uv_doc, +"nest2ang_uv(nside, ipix, u, v, theta=None, phi=None, /)\n" +"--\n" +"\n"); + + +static PyObject* cnest2ang_uv(PyObject* self, PyObject* args) { + Py_ssize_t nside; + PyObject* op[] = {NULL, NULL, NULL, NULL, NULL}; + int types[] = {NPY_INT64, NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE}; + + if (!PyArg_ParseTuple(args, "nOOO|OO:nest2ang_uv", &nside, + &op[0], &op[1], &op[2], &op[3], &op[4])) { + return NULL; + } + + return vectorize(vnest2ang_uv, &nside, 3, 2, op, types); +} + + +static void vring2ang_uv(void* args, npy_intp size, void** data) { + int64_t nside = *(int64_t*)args; + int64_t* ipix = data[0]; + double* u = data[1], *v = data[2]; + double* theta = data[3], *phi = data[4]; + for (npy_intp i = 0; i < size; ++i) + setang(ring2ang_uv(nside, ipix[i], u[i], v[i]), &theta[i], &phi[i]); +} + + +PyDoc_STRVAR(cring2ang_uv_doc, +"ring2ang_uv(nside, ipix, u, v, theta=None, phi=None, /)\n" +"--\n" +"\n"); + + +static PyObject* cring2ang_uv(PyObject* self, PyObject* args) { + Py_ssize_t nside; + PyObject* op[] = {NULL, NULL, NULL, NULL, NULL}; + int types[] = {NPY_INT64, NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE}; + + if (!PyArg_ParseTuple(args, "nOOO|OO:ring2ang_uv", &nside, + &op[0], &op[1], &op[2], &op[3], &op[4])) + return NULL; + + return vectorize(vring2ang_uv, &nside, 3, 2, op, types); +} + + +static void vvec2nest_uv(void* args, npy_intp size, void** data) { + int64_t nside = *(int64_t*)args; + double* x = data[0], *y = data[1], *z = data[2]; + int64_t* ipix = data[3]; + double* u = data[4], *v = data[5]; + for (npy_intp i = 0; i < size; ++i) + ipix[i] = vec2nest_uv(nside, (t_vec){x[i], y[i], z[i]}, &u[i], &v[i]); +} + + +PyDoc_STRVAR(cvec2nest_uv_doc, +"vec2nest_uv(nside, x, y, z, ipix=None, u=None, v=None, /)\n" +"--\n" +"\n"); + + +static PyObject* cvec2nest_uv(PyObject* self, PyObject* args) { + Py_ssize_t nside; + PyObject* op[] = {NULL, NULL, NULL, NULL, NULL, NULL}; + int types[] = {NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE, NPY_INT64, NPY_DOUBLE, + NPY_DOUBLE}; + + if (!PyArg_ParseTuple(args, "nOOO|OOO:vec2nest_uv", &nside, + &op[0], &op[1], &op[2], &op[3], &op[4], &op[5])) { + return NULL; + } + + return vectorize(vvec2nest_uv, &nside, 3, 3, op, types); +} + + +static void vvec2ring_uv(void* args, npy_intp size, void** data) { + int64_t nside = *(int64_t*)args; + double* x = data[0], *y = data[1], *z = data[2]; + int64_t* ipix = data[3]; + double* u = data[4], *v = data[5]; + for (npy_intp i = 0; i < size; ++i) + ipix[i] = vec2ring_uv(nside, (t_vec){x[i], y[i], z[i]}, &u[i], &v[i]); +} + + +PyDoc_STRVAR(cvec2ring_uv_doc, +"vec2ring_uv(nside, x, y, z, ipix=None, u=None, v=None, /)\n" +"--\n" +"\n"); + + +static PyObject* cvec2ring_uv(PyObject* self, PyObject* args) { + Py_ssize_t nside; + PyObject* op[] = {NULL, NULL, NULL, NULL, NULL, NULL}; + int types[] = {NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE, NPY_INT64, NPY_DOUBLE, + NPY_DOUBLE}; + + if (!PyArg_ParseTuple(args, "nOOO|OOO:vec2ring_uv", &nside, + &op[0], &op[1], &op[2], &op[3], &op[4], &op[5])) { + return NULL; + } + + return vectorize(vvec2ring_uv, &nside, 3, 3, op, types); +} + + +static void vnest2vec_uv(void* args, npy_intp size, void** data) { + int64_t nside = *(int64_t*)args; + int64_t* ipix = data[0]; + double* u = data[1], *v = data[2]; + double* x = data[3], *y = data[4], *z = data[5]; + for (npy_intp i = 0; i < size; ++i) + setvec(nest2vec_uv(nside, ipix[i], u[i], v[i]), &x[i], &y[i], &z[i]); +} + + +PyDoc_STRVAR(cnest2vec_uv_doc, +"nest2vec_uv(nside, ipix, u, v, x=None, y=None, z=None, /)\n" +"--\n" +"\n"); + + +static PyObject* cnest2vec_uv(PyObject* self, PyObject* args) { + Py_ssize_t nside; + PyObject* op[] = {NULL, NULL, NULL, NULL, NULL, NULL}; + int types[] = {NPY_INT64, NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE, + NPY_DOUBLE}; + + if (!PyArg_ParseTuple(args, "nOOO|OOO:nest2vec_uv", &nside, + &op[0], &op[1], &op[2], &op[3], &op[4], &op[5])) { + return NULL; + } + + return vectorize(vnest2vec_uv, &nside, 3, 3, op, types); +} + + +static void vring2vec_uv(void* args, npy_intp size, void** data) { + int64_t nside = *(int64_t*)args; + int64_t* ipix = data[0]; + double* u = data[1], *v = data[2]; + double* x = data[3], *y = data[4], *z = data[5]; + for (npy_intp i = 0; i < size; ++i) + setvec(ring2vec_uv(nside, ipix[i], u[i], v[i]), &x[i], &y[i], &z[i]); +} + + +PyDoc_STRVAR(cring2vec_uv_doc, +"ring2vec_uv(nside, ipix, u, v, x=None, y=None, z=None, /)\n" +"--\n" +"\n"); + + +static PyObject* cring2vec_uv(PyObject* self, PyObject* args) { + Py_ssize_t nside; + PyObject* op[] = {NULL, NULL, NULL, NULL, NULL, NULL}; + int types[] = {NPY_INT64, NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE, + NPY_DOUBLE}; + + if (!PyArg_ParseTuple(args, "nOOO|OOO:ring2vec_uv", &nside, + &op[0], &op[1], &op[2], &op[3], &op[4], &op[5])) { + return NULL; + } + + return vectorize(vring2vec_uv, &nside, 3, 3, op, types); +} + + +PyDoc_STRVAR(cnside2npix_doc, "nside2npix(nside, /)\n--\n\n"); + + +static PyObject* cnside2npix(PyObject* self, PyObject* args) { + Py_ssize_t nside; + int64_t npix; + + if (!PyArg_ParseTuple(args, "n:nside2npix", &nside)) + return NULL; + + npix = nside2npix(nside); + + return Py_BuildValue("n", (Py_ssize_t)npix); +} + + +PyDoc_STRVAR(cnpix2nside_doc, "npix2nside(npix, /)\n--\n\n"); + + +static PyObject* cnpix2nside(PyObject* self, PyObject* args) { + Py_ssize_t npix; + int64_t nside; + + if (!PyArg_ParseTuple(args, "n:npix2nside", &npix)) + return NULL; + + nside = npix2nside(npix); + + return Py_BuildValue("n", (Py_ssize_t)nside); +} + + +static PyMethodDef methods[] = { + {"ang2vec", cang2vec, METH_VARARGS, cang2vec_doc}, + {"vec2ang", cvec2ang, METH_VARARGS, cvec2ang_doc}, + {"ang2nest", cang2nest, METH_VARARGS, cang2nest_doc}, + {"ang2ring", cang2ring, METH_VARARGS, cang2ring_doc}, + {"nest2ang", cnest2ang, METH_VARARGS, cnest2ang_doc}, + {"ring2ang", cring2ang, METH_VARARGS, cring2ang_doc}, + {"vec2nest", cvec2nest, METH_VARARGS, cvec2nest_doc}, + {"vec2ring", cvec2ring, METH_VARARGS, cvec2ring_doc}, + {"nest2vec", cnest2vec, METH_VARARGS, cnest2vec_doc}, + {"ring2vec", cring2vec, METH_VARARGS, cring2vec_doc}, + {"ang2nest_uv", cang2nest_uv, METH_VARARGS, cang2nest_uv_doc}, + {"ang2ring_uv", cang2ring_uv, METH_VARARGS, cang2ring_uv_doc}, + {"nest2ang_uv", cnest2ang_uv, METH_VARARGS, cnest2ang_uv_doc}, + {"ring2ang_uv", cring2ang_uv, METH_VARARGS, cring2ang_uv_doc}, + {"vec2nest_uv", cvec2nest_uv, METH_VARARGS, cvec2nest_uv_doc}, + {"vec2ring_uv", cvec2ring_uv, METH_VARARGS, cvec2ring_uv_doc}, + {"nest2vec_uv", cnest2vec_uv, METH_VARARGS, cnest2vec_uv_doc}, + {"ring2vec_uv", cring2vec_uv, METH_VARARGS, cring2vec_uv_doc}, + {"nside2npix", cnside2npix, METH_VARARGS, cnside2npix_doc}, + {"npix2nside", cnpix2nside, METH_VARARGS, cnpix2nside_doc}, + {NULL, NULL} +}; + + +static struct PyModuleDef module_def = { + PyModuleDef_HEAD_INIT, + "chealpix", + PyDoc_STR("healpix C library interface"), + -1, + methods +}; + + +PyMODINIT_FUNC PyInit_chealpix(void) { + PyObject* module = PyModule_Create(&module_def); + if (!module) + return NULL; + PyModule_AddIntConstant(module, "NSIDE_MAX", NSIDE_MAX); + import_array(); + if (PyErr_Occurred()) + return NULL; + return module; +} diff --git a/python/healpix/__init__.py b/python/healpix/__init__.py new file mode 100644 index 0000000..716d006 --- /dev/null +++ b/python/healpix/__init__.py @@ -0,0 +1,178 @@ +# author: Nicolas Tessore +# license: BSD-3-Clause +'''Python package for HEALPix discretisation of the sphere''' + +import numpy as np +import chealpix as _chp + + +NSIDE_MAX = _chp.NSIDE_MAX +'''Maximum admissible value for the NSIDE parameter. + +If nside > NSIDE_MAX is used, the resulting pixel indices can overflow their +underlying integer type in the C library functions. + +''' + +def is_power_of_two(n): + '''Return True if n is a power of two, or False otherwise.''' + return (n & (n-1)) == 0 + + +def check_nside(nside, nest=False): + '''Check whether nside is a valid resolution parameter.''' + if nside < 1 or nside > NSIDE_MAX: + raise ValueError('nside must be a positive integer 1 <= nside <= ' + f'2^{np.log2(NSIDE_MAX):.0f}') + if nest and (nside & (nside-1)) != 0: + raise ValueError('nside must be power of two in the NEST scheme'); + + +def thetaphi_from_lonlat(lon, lat, theta=None, phi=None): + '''convert longitude and latitude in degree to theta and phi in radian''' + theta, phi = np.deg2rad(lat, out=theta), np.deg2rad(lon, out=phi) + theta *= -1 + theta += np.pi/2 + return theta, phi + + +def lonlat_from_thetaphi(theta, phi, lon=None, lat=None): + '''convert theta and phi in radian to longitude and latitude in degree''' + lon, lat = np.rad2deg(phi, out=lon), np.rad2deg(theta, out=lat) + lat *= -1 + lat += 90 + return lon, lat + + +def ang2vec(theta, phi, lonlat=False): + out1, out2 = None, None + if lonlat: + s = np.broadcast(theta, phi).shape + out1, out2 = np.empty(s), np.empty(s) + theta, phi = thetaphi_from_lonlat(theta, phi, theta=out1, phi=out2) + return _chp.ang2vec(theta, phi, out1, out2) + + +def vec2ang(x, y, z, lonlat=False): + theta, phi = _chp.vec2ang(x, y, z) + if lonlat: + out1, out2 = (phi, theta) if np.ndim(theta) else (None, None) + theta, phi = lonlat_from_thetaphi(theta, phi, lon=out1, lat=out2) + return theta, phi + + +def ang2pix(nside, theta, phi, nest=False, lonlat=False): + check_nside(nside, nest) + if not lonlat: + if nest: + return _chp.ang2nest(nside, theta, phi) + else: + return _chp.ang2ring(nside, theta, phi) + with np.nditer([theta, phi, None], + ['buffered', 'external_loop', 'zerosize_ok'], + [['readonly']]*2 + [['writeonly', 'allocate']]*1, + [None, None, int]) as it: + theta, phi = None, None + for lon, lat, ipix in it: + theta, phi = thetaphi_from_lonlat(lon, lat, theta, phi) + if nest: + _chp.ang2nest(nside, theta, phi, ipix) + else: + _chp.ang2ring(nside, theta, phi, ipix) + return it.operands[2] + + +def pix2ang(nside, ipix, nest=False, lonlat=False): + if nest: + theta, phi = _chp.nest2ang(nside, ipix) + else: + theta, phi = _chp.ring2ang(nside, ipix) + if lonlat: + out1, out2 = (phi, theta) if np.ndim(theta) else (None, None) + theta, phi = lonlat_from_thetaphi(theta, phi, lon=out1, lat=out2) + return theta, phi + + +def vec2pix(nside, x, y, z, nest=False): + check_nside(nside, nest) + if nest: + return _chp.vec2nest(nside, x, y, z) + else: + return _chp.vec2ring(nside, x, y, z) + + +def pix2vec(nside, ipix, nest=False): + if nest: + return _chp.nest2vec(nside, ipix) + else: + return _chp.ring2vec(nside, ipix) + + +def nside2npix(nside): + return _chp.nside2npix(nside) + + +def npix2nside(npix): + return _chp.npix2nside(npix) + + +def randang(nside, ipix, nest=False, lonlat=False, rng=None): + '''Sample random spherical coordinates from the given HEALPix pixels. + + This function produces one pair of random spherical coordinates from each + HEALPix pixel listed in `ipix`, which can be scalar or array-like. The + indices use the `nside` resolution parameter and either the RING scheme + (if `nest` is false, the default) or the NEST scheme (if `nest` is true). + + Returns either a tuple `theta, phi` of mathematical coordinates in radians + if `lonlat` is False (the default), or a tuple `lon, lat` of longitude and + latitude in degrees if `lonlat` is True. The output is of the same shape + as the input. + + An optional numpy random number generator can be provided using `rng`; + otherwise, a new numpy.random.default_rng() is used. + ''' + + if rng is None: + rng = np.random.default_rng() + + u = rng.random(ipix.shape, np.double) + v = rng.random(ipix.shape, np.double) + + if nest: + theta, phi = _chp.nest2ang_uv(nside, ipix, u, v, u, v) + else: + theta, phi = _chp.ring2ang_uv(nside, ipix, u, v, u, v) + + if lonlat: + out1, out2 = (phi, theta) if np.ndim(theta) else (None, None) + theta, phi = lonlat_from_thetaphi(theta, phi, lon=out1, lat=out2) + + return theta, phi + + +def randvec(nside, ipix, nest=False, rng=None): + '''Sample random unit vectors from the given HEALPix pixels. + + This function produces one random unit vector from each HEALPix pixel + listed in `ipix`, which can be scalar or array-like. The pixel indices use + the `nside` resolution parameter and either the RING scheme (if `nest` is + false, the default) or the NEST scheme (if `nest` is true). + + Returns a tuple `x, y, z` of normalised vector components with the same + shape as the input. + + An optional numpy random number generator can be provided using `rng`; + otherwise, a new numpy.random.default_rng() is used. + ''' + + if rng is None: + rng = np.random.default_rng() + + u = rng.random(ipix.shape, np.double) + v = rng.random(ipix.shape, np.double) + + if nest: + return _chp.nest2vec_uv(nside, ipix, u, v, u, v) + else: + return _chp.ring2vec_uv(nside, ipix, u, v, u, v) diff --git a/python/healpix/test/__init__.py b/python/healpix/test/__init__.py new file mode 100644 index 0000000..9df55c4 --- /dev/null +++ b/python/healpix/test/__init__.py @@ -0,0 +1,11 @@ +import numpy as np + + +def assert_shape(a, shape): + s = np.shape(a) + if shape is None: + assert s == (), f'expected scalar, got shape {s}' + else: + if not isinstance(shape, tuple): + shape = (shape,) + assert s == shape, f'expected shape {shape}, got shape {s}' diff --git a/python/healpix/test/conftest.py b/python/healpix/test/conftest.py new file mode 100644 index 0000000..76528b0 --- /dev/null +++ b/python/healpix/test/conftest.py @@ -0,0 +1,59 @@ +import pytest +import numpy as np + +np.random.seed(12345) + + +def is_power_of_two(n): + return (n & (n-1)) == 0 + + +def fixture_get(request, name, default): + if name in request.fixturenames: + return request.getfixturevalue(name) + else: + return default + + +@pytest.fixture(params=[1, 15, 16, 128]) +def nside(request): + nside = request.param + if (fixture_get(request, 'nest', False) or request.node.get_closest_marker('nest')) \ + and not is_power_of_two(nside): + pytest.skip(reason='nside is not power of two') + return nside + + +@pytest.fixture(params=[True, False]) +def nest(request): + return request.param + + +@pytest.fixture(params=[True, False]) +def lonlat(request): + return request.param + + +@pytest.fixture(params=[None, (), 1, 1000, (7, 11)]) +def size(request): + return request.param + + +@pytest.fixture +def base_pixel_ang(request): + theta = np.repeat(np.arccos([2/3, 0, -2/3]), 4) + phi = np.pi/4*np.array([1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7]) + if fixture_get(request, 'lonlat', False): + lon = np.rad2deg(phi) + lat = np.rad2deg(np.pi/2 - theta) + theta, phi = lon, lat + return theta, phi + + +@pytest.fixture +def base_pixel_vec(request): + a, b = (5/18)**0.5, 2/3 + x = np.array([a, -a, -a, a, 1, 0, -1, 0, a, -a, -a, a]) + y = np.array([a, a, -a, -a, 0, 1, 0, -1, a, a, -a, -a]) + z = np.array([b, b, b, b, 0, 0, 0, 0, -b, -b, -b, -b]) + return x, y, z diff --git a/python/healpix/test/test_chealpix.py b/python/healpix/test/test_chealpix.py new file mode 100644 index 0000000..967979d --- /dev/null +++ b/python/healpix/test/test_chealpix.py @@ -0,0 +1,105 @@ +import pytest +import numpy as np +import numpy.testing as npt + + +def test_ang2vec(): + from chealpix import ang2vec + theta = np.arccos(np.random.uniform(-1, 1, size=100)) + phi = np.random.uniform(0, 2*np.pi, size=100) + x, y, z = ang2vec(theta, phi) + npt.assert_allclose(x, np.sin(theta)*np.cos(phi)) + npt.assert_allclose(y, np.sin(theta)*np.sin(phi)) + npt.assert_allclose(z, np.cos(theta)) + + +def test_vec2ang(): + from chealpix import vec2ang + x, y, z = np.random.randn(3, 100) + theta, phi = vec2ang(x, y, z) + npt.assert_allclose(theta, np.arctan2(np.hypot(x, y), z)) + npt.assert_allclose(phi, np.arctan2(y, x)) + + +@pytest.mark.nest +def test_ang_nest(nside): + from chealpix import ang2nest, nest2ang + ipix = np.arange(12*nside**2) + npt.assert_array_equal(ang2nest(nside, *nest2ang(nside, ipix)), ipix) + + +def test_ang_ring(nside): + from chealpix import ang2ring, ring2ang + ipix = np.arange(12*nside**2) + npt.assert_array_equal(ang2ring(nside, *ring2ang(nside, ipix)), ipix) + + +@pytest.mark.nest +def test_vec_nest(nside): + from chealpix import vec2nest, nest2vec + ipix = np.arange(12*nside**2) + npt.assert_array_equal(vec2nest(nside, *nest2vec(nside, ipix)), ipix) + + +def test_vec_ring(nside): + from chealpix import vec2ring, ring2vec + ipix = np.arange(12*nside**2) + npt.assert_array_equal(vec2ring(nside, *ring2vec(nside, ipix)), ipix) + + +@pytest.mark.nest +def test_ang_nest_uv(nside): + from chealpix import ang2nest_uv, nest2ang_uv + ipix = np.arange(12*nside**2) + u, v = np.random.rand(2, ipix.size) + result = ang2nest_uv(nside, *nest2ang_uv(nside, ipix, u, v)) + npt.assert_array_equal(result[0], ipix) + npt.assert_allclose(result[1], u) + npt.assert_allclose(result[2], v) + + +def test_ang_ring_uv(nside): + from chealpix import ang2ring_uv, ring2ang_uv + ipix = np.arange(12*nside**2) + u, v = np.random.rand(2, ipix.size) + result = ang2ring_uv(nside, *ring2ang_uv(nside, ipix, u, v)) + npt.assert_array_equal(result[0], ipix) + npt.assert_allclose(result[1], u) + npt.assert_allclose(result[2], v) + + +@pytest.mark.nest +def test_vec_nest_uv(nside): + from chealpix import vec2nest_uv, nest2vec_uv + ipix = np.arange(12*nside**2) + u, v = np.random.rand(2, ipix.size) + result = vec2nest_uv(nside, *nest2vec_uv(nside, ipix, u, v)) + npt.assert_array_equal(result[0], ipix) + npt.assert_allclose(result[1], u) + npt.assert_allclose(result[2], v) + + +def test_vec_ring_uv(nside): + from chealpix import vec2ring_uv, ring2vec_uv + ipix = np.arange(12*nside**2) + u, v = np.random.rand(2, ipix.size) + result = vec2ring_uv(nside, *ring2vec_uv(nside, ipix, u, v)) + npt.assert_array_equal(result[0], ipix) + npt.assert_allclose(result[1], u) + npt.assert_allclose(result[2], v) + + +def test_nside2npix(nside): + from chealpix import nside2npix + assert nside2npix(nside) == 12*nside**2 + + +def test_npix2nside(nside): + from chealpix import npix2nside + assert npix2nside(12*nside**2) == nside + + +def test_npix2nside_invalid(): + from chealpix import npix2nside + assert npix2nside(7) == -1 + assert npix2nside(49) == -1 diff --git a/python/healpix/test/test_healpix.py b/python/healpix/test/test_healpix.py new file mode 100644 index 0000000..c75ac6a --- /dev/null +++ b/python/healpix/test/test_healpix.py @@ -0,0 +1,75 @@ +import pytest +import numpy as np +import numpy.testing as npt + +from . import assert_shape + + +def test_ang2vec(lonlat, size): + from healpix import ang2vec + if lonlat: + theta_lon = np.random.uniform(0, 360, size=size) + phi_lat = np.random.uniform(-90, 90, size=size) + x_ = np.cos(np.deg2rad(phi_lat))*np.cos(np.deg2rad(theta_lon)) + y_ = np.cos(np.deg2rad(phi_lat))*np.sin(np.deg2rad(theta_lon)) + z_ = np.sin(np.deg2rad(phi_lat)) + else: + theta_lon = np.random.uniform(0, np.pi, size=size) + phi_lat = np.random.uniform(0, 2*np.pi, size=size) + x_ = np.sin(theta_lon)*np.cos(phi_lat) + y_ = np.sin(theta_lon)*np.sin(phi_lat) + z_ = np.cos(theta_lon) + x, y, z = ang2vec(theta_lon, phi_lat, lonlat=lonlat) + npt.assert_allclose(x, x_) + npt.assert_allclose(y, y_) + npt.assert_allclose(z, z_) + assert_shape(x, size) + assert_shape(y, size) + assert_shape(z, size) + + +def test_vec2ang(lonlat, size): + from healpix import vec2ang + x = np.random.standard_normal(size) + y = np.random.standard_normal(size) + z = np.random.standard_normal(size) + theta_, phi_ = np.arctan2(np.hypot(x, y), z), np.arctan2(y, x) + if lonlat: + theta_, phi_ = np.rad2deg(phi_), np.rad2deg(np.pi/2-theta_) + theta, phi = vec2ang(x, y, z, lonlat=lonlat) + npt.assert_allclose(theta, theta_) + npt.assert_allclose(phi, phi_) + assert_shape(theta, size) + assert_shape(phi, size) + + +def test_ang2pix_base(base_pixel_ang, nest, lonlat): + from healpix import ang2pix + npt.assert_array_equal(ang2pix(1, *base_pixel_ang, nest, lonlat), np.arange(12)) + + +def test_pix2ang_base(base_pixel_ang, nest, lonlat): + from healpix import pix2ang + npt.assert_allclose(pix2ang(1, np.arange(12), nest, lonlat), base_pixel_ang) + + +def test_ang_pix(nside, nest, lonlat): + from healpix import ang2pix, pix2ang + ipix = np.arange(12*nside**2) + npt.assert_array_equal(ang2pix(nside, *pix2ang(nside, ipix, nest, lonlat), nest, lonlat), ipix) + + +def test_vec2pix_base(base_pixel_vec, nest): + from healpix import vec2pix + npt.assert_array_equal(vec2pix(1, *base_pixel_vec, nest), np.arange(12)) + + +def test_pix2vec_base(base_pixel_vec, nest): + from healpix import pix2vec + npt.assert_allclose(pix2vec(1, np.arange(12), nest), base_pixel_vec, atol=1e-15) + + +def test_vec_pix(nside, nest): + from healpix import vec2pix, pix2vec + ipix = np.arange(12*nside**2) + npt.assert_array_equal(vec2pix(nside, *pix2vec(nside, ipix, nest), nest), ipix) diff --git a/setup.cfg b/setup.cfg new file mode 100644 index 0000000..8fa4346 --- /dev/null +++ b/setup.cfg @@ -0,0 +1,31 @@ +[metadata] +name = healpix +version = 2022.10.31 +maintainer = Nicolas Tessore +maintainer_email = n.tessore@ucl.ac.uk +description = Python package for HEALPix discretisation of the sphere +long_description = file: README.md +long_description_content_type = text/markdown +license = BSD-3-Clause +license_files = + LICENSE.md +url = https://github.com/ntessore/healpix +classifiers = + Programming Language :: Python :: 3 + License :: OSI Approved :: BSD License + Operating System :: OS Independent + +[options] +package_dir = + =python +packages = + healpix + healpix.test +python_requires = >=3.6 +install_requires = + numpy + + +[tool:pytest] +markers = + nest: test uses NEST scheme diff --git a/setup.py b/setup.py new file mode 100644 index 0000000..edb35c7 --- /dev/null +++ b/setup.py @@ -0,0 +1,15 @@ +from setuptools import setup, Extension +import os.path +import numpy as np + + +setup( + ext_modules=[ + Extension('chealpix', + sources=['python/chealpix.c', 'src/healpix.c'], + include_dirs=[np.get_include(), 'src'], + define_macros=[('NPY_NO_DEPRECATED_API', 'NPY_1_7_API_VERSION')], + extra_compile_args=['-std=c99'], + ), + ], +) diff --git a/src/healpix.c b/src/healpix.c new file mode 100644 index 0000000..68dc943 --- /dev/null +++ b/src/healpix.c @@ -0,0 +1,409 @@ +// based on the healpix_bare library +// license: BSD-3-Clause + +#include +#include "healpix.h" + +static const double PI = 3.141592653589793238462643383279502884197; + +static const int jrll[] = { 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4 }; +static const int jpll[] = { 1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7 }; + + +// maximum nside for 64bit signed integer types +const int64_t NSIDE_MAX = 1<<29; + + +/* conversions between continuous coordinate systems */ + +// A structure describing a location in cylindrical coordinates. +// z = cos(theta), s = sin(theta), phi +typedef struct { + double z, s, phi; +} t_loc; + + +// A structure describing the continuous Healpix coordinate system. +// f takes values in [0, 11], u and v lie in [0.0, 1.0]. +typedef struct { + double u, v; + int32_t f; +} t_hpc; + + +static t_hpc loc2hpc(t_loc loc) { + double za = fabs(loc.z); + double x = loc.phi*(1./(2.*PI)); + if (x < 0.) { + x += (int64_t)x + 1.; + } else if (x >= 1.) { + x -= (int64_t)x; + } + double tt = 4.*x; + + if (za <= 2./3.) { + // Equatorial region + double temp1 = 0.5+tt; // [0.5; 4.5) + double temp2 = loc.z*0.75; // [-0.5; +0.5] + double jp = temp1-temp2; // index of ascending edge line, [0; 5) + double jm = temp1+temp2; // index of descending edge line, [0; 5) + int ifp = (int)jp; // in {0,4} + int ifm = (int)jm; + return (t_hpc){ + jm-ifm, + 1+ifp - jp, + (ifp==ifm) ? (ifp|4) : ((ifp= 4) + ntt = 3; + double tp = tt - ntt; // [0;1) + double tmp = sqrt(3.*(1.-za)); + + double jp = tp*tmp; // increasing edge line index + double jm = (1.0-tp)*tmp; // decreasing edge line index + if (jp > 1.) + jp = 1.; // for points too close to the boundary + if (jm > 1.) + jm = 1.; + return (loc.z >= 0) ? (t_hpc){1.-jm, 1.-jp, ntt} + : (t_hpc){jp, jm, ntt+8}; +} + + +static t_loc hpc2loc(t_hpc hpc) { + double z, s, phi; + const double x = hpc.u - hpc.v; + const double y = (hpc.u - 0.5) + (hpc.v - 0.5); + const int r = 1 - hpc.f/4; + const double ry = r*y; + if (ry > 0.) { + // polar cap + const double tmp = (1.-ry)*(1.-ry)*(1./3.); + z = r*(1. - tmp); + s = sqrt(tmp*(2.-tmp)); + phi = (PI/4.)*(jpll[hpc.f] + x/(1.-ry)); + } else { + // equatorial region + z = (y+r)*(2./3.); + s = sqrt((1.+z)*(1.-z)); + phi = (PI/4.)*(jpll[hpc.f] + x); + } + return (t_loc){z, s, phi}; +} + + +static t_loc ang2loc(t_ang ang) { + double z = cos(ang.theta), s = sin(ang.theta); + if (s < 0) { + s = -s; + ang.phi += PI; + } + return (t_loc){z, s, ang.phi}; +} + + +static t_ang loc2ang(t_loc loc) { + return (t_ang){atan2(loc.s, loc.z), loc.phi}; +} + + +static t_loc vec2loc(t_vec vec) { + double s = hypot(vec.x, vec.y); + double r = hypot(s, vec.z); + return (t_loc){vec.z/r, s/r, atan2(vec.y, vec.x)}; +} + + +static t_vec loc2vec(t_loc loc) { + return (t_vec){loc.s*cos(loc.phi), loc.s*sin(loc.phi), loc.z}; +} + + +t_vec ang2vec(t_ang ang) { + return loc2vec(ang2loc(ang)); +} + + +t_ang vec2ang(t_vec vec) { + return (t_ang){atan2(hypot(vec.x, vec.y), vec.z), atan2(vec.y, vec.x)}; +} + + +// conversions between discrete coordinate systems + + +static int64_t isqrt(int64_t v) { + int64_t res = sqrt(v+0.5); + if (v < ((int64_t)(1)<<50)) + return res; + if (res*res > v) { + --res; + } else if ((res+1)*(res+1) <= v) { + ++res; + } + return res; +} + + +static int64_t spread_bits(int64_t v) { + int64_t res = v & 0xffffffff; + res = (res^(res<<16)) & 0x0000ffff0000ffff; + res = (res^(res<< 8)) & 0x00ff00ff00ff00ff; + res = (res^(res<< 4)) & 0x0f0f0f0f0f0f0f0f; + res = (res^(res<< 2)) & 0x3333333333333333; + res = (res^(res<< 1)) & 0x5555555555555555; + return res; +} + + +static int64_t compress_bits(int64_t v) { + int64_t res = v & 0x5555555555555555; + res = (res^(res>> 1)) & 0x3333333333333333; + res = (res^(res>> 2)) & 0x0f0f0f0f0f0f0f0f; + res = (res^(res>> 4)) & 0x00ff00ff00ff00ff; + res = (res^(res>> 8)) & 0x0000ffff0000ffff; + res = (res^(res>>16)) & 0x00000000ffffffff; + return res; +} + + +// A structure describing the discrete Healpix coordinate system. +// f takes values in [0, 11], x and y lie in [0, nside). +typedef struct { + int64_t x, y; + int32_t f; +} t_hpd; + + +static int64_t hpd2nest(int64_t nside, t_hpd hpd) { + return (hpd.f*nside*nside) + spread_bits(hpd.x) + (spread_bits(hpd.y)<<1); +} + + +static t_hpd nest2hpd(int64_t nside, int64_t pix) { + int64_t npface_=nside*nside, p2=pix&(npface_-1); + return (t_hpd){compress_bits(p2), compress_bits(p2>>1), pix/npface_}; +} + + +static int64_t hpd2ring (int64_t nside_, t_hpd hpd) + { + int64_t nl4 = 4*nside_; + int64_t jr = (jrll[hpd.f]*nside_) - hpd.x - hpd.y - 1; + + if (jrnl4) ? jp-nl4 : ((jp<1) ? jp+nl4 : jp); + return 2*jr*(jr-1) + jp - 1; + } + else if (jr > 3*nside_) + { + jr = nl4-jr; + int64_t jp = (jpll[hpd.f]*jr + hpd.x - hpd.y + 1) / 2; + jp = (jp>nl4) ? jp-nl4 : ((jp<1) ? jp+nl4 : jp); + return 12*nside_*nside_ - 2*(jr+1)*jr + jp - 1; + } + else + { + int64_t jp = (jpll[hpd.f]*nside_ + hpd.x - hpd.y + 1 + ((jr-nside_)&1)) / 2; + jp = (jp>nl4) ? jp-nl4 : ((jp<1) ? jp+nl4 : jp); + return 2*nside_*(nside_-1) + (jr-nside_)*nl4 + jp - 1; + } + } + +static t_hpd ring2hpd (int64_t nside_, int64_t pix) + { + int64_t ncap_=2*nside_*(nside_-1); + int64_t npix_=12*nside_*nside_; + + if (pix>1; /* counted from North pole */ + int64_t iphi = (pix+1) - 2*iring*(iring-1); + int64_t face = (iphi-1)/iring; + int64_t irt = iring - (jrll[face]*nside_) + 1; + int64_t ipt = 2*iphi- jpll[face]*iring -1; + if (ipt>=2*nside_) ipt-=8*nside_; + return (t_hpd) {(ipt-irt)>>1, (-(ipt+irt))>>1, face}; + } + else if (pix<(npix_-ncap_)) /* Equatorial region */ + { + int64_t ip = pix - ncap_; + int64_t iring = (ip/(4*nside_)) + nside_; /* counted from North pole */ + int64_t iphi = (ip%(4*nside_)) + 1; + int64_t kshift = (iring+nside_)&1; + int64_t ire = iring-nside_+1; + int64_t irm = 2*nside_+2-ire; + int64_t ifm = (iphi - ire/2 + nside_ -1) / nside_; + int64_t ifp = (iphi - irm/2 + nside_ -1) / nside_; + int64_t face = (ifp==ifm) ? (ifp|4) : ((ifp=2*nside_) ipt-=8*nside_; + return (t_hpd) {(ipt-irt)>>1, (-(ipt+irt))>>1, face}; + } + else /* South Polar cap */ + { + int64_t ip = npix_ - pix; + int64_t iring = (1+isqrt(2*ip-1))>>1; /* counted from South pole */ + int64_t iphi = 4*iring + 1 - (ip - 2*iring*(iring-1)); + int64_t face=8+(iphi-1)/iring; + int64_t irt = 4*nside_ - iring - (jrll[face]*nside_) + 1; + int64_t ipt = 2*iphi- jpll[face]*iring -1; + if (ipt>=2*nside_) ipt-=8*nside_; + return (t_hpd) {(ipt-irt)>>1, (-(ipt+irt))>>1, face}; + } + } + + +int64_t nest2ring(int64_t nside, int64_t ipnest) { + // power of two check + if ((nside & (nside-1)) != 0) + return -1; + return hpd2ring(nside, nest2hpd(nside, ipnest)); +} + + +int64_t ring2nest(int64_t nside, int64_t ipring) { + // power of two check + if ((nside&(nside-1)) != 0) + return -1; + return hpd2nest(nside, ring2hpd(nside, ipring)); +} + + +// mixed conversions + + +static t_hpd loc2hpd(int64_t nside, t_loc loc) { + t_hpc hpc = loc2hpc(loc); + return (t_hpd){hpc.u*nside, hpc.v*nside, hpc.f}; +} + + +static t_loc hpd2loc(int64_t nside, t_hpd hpd) { + return hpc2loc((t_hpc){(hpd.x + 0.5)/nside, (hpd.y + 0.5)/nside, hpd.f}); +} + + +int64_t npix2nside(int64_t npix) { + int64_t res = isqrt(npix/12); + return (res*res*12 == npix) ? res : -1; +} + + +int64_t nside2npix(int64_t nside) { + return 12*nside*nside; +} + + +double vec_angle(t_vec v1, t_vec v2) { + t_vec cross = { + v1.y*v2.z - v1.z*v2.y, + v1.z*v2.x - v1.x*v2.z, + v1.x*v2.y - v1.y*v2.x + }; + double len = sqrt(cross.x*cross.x + cross.y*cross.y + cross.z*cross.z); + double dot = v1.x*v2.x + v1.y*v2.y + v1.z*v2.z; + return atan2(len, dot); +} + + +int64_t ang2ring(int64_t nside, t_ang ang) { + return hpd2ring(nside, loc2hpd(nside, ang2loc(ang))); +} + + +int64_t ang2nest(int64_t nside, t_ang ang) { + return hpd2nest(nside, loc2hpd(nside, ang2loc(ang))); +} + + +int64_t vec2ring(int64_t nside, t_vec vec) { + return hpd2ring(nside, loc2hpd(nside, vec2loc(vec))); +} + + +int64_t vec2nest(int64_t nside, t_vec vec) { + return hpd2nest(nside, loc2hpd(nside, vec2loc(vec))); +} + + +t_ang ring2ang(int64_t nside, int64_t ipix) { + return loc2ang(hpd2loc(nside, ring2hpd(nside, ipix))); +} + + +t_ang nest2ang(int64_t nside, int64_t ipix) { + return loc2ang(hpd2loc(nside, nest2hpd(nside, ipix))); +} + + +t_vec ring2vec(int64_t nside, int64_t ipix) { + return loc2vec(hpd2loc(nside, ring2hpd(nside, ipix))); +} + + +t_vec nest2vec(int64_t nside, int64_t ipix) { + return loc2vec(hpd2loc(nside, nest2hpd(nside, ipix))); +} + + +// conversions with sub-pixel positions + + +static t_hpd loc2hpd_uv(int64_t nside, t_loc loc, double* u, double* v) { + t_hpc hpc = loc2hpc(loc); + *u = modf(hpc.u*nside, &hpc.u); + *v = modf(hpc.v*nside, &hpc.v); + return (t_hpd){hpc.u+0.1, hpc.v+0.1, hpc.f}; +} + + +static t_loc hpd2loc_uv(int64_t nside, t_hpd hpd, double u, double v) { + return hpc2loc((t_hpc){(hpd.x + u)/nside, (hpd.y + v)/nside, hpd.f}); +} + + +int64_t ang2ring_uv(int64_t nside, t_ang ang, double* u, double* v) { + return hpd2ring(nside, loc2hpd_uv(nside, ang2loc(ang), u, v)); +} + + +int64_t ang2nest_uv(int64_t nside, t_ang ang, double* u, double* v) { + return hpd2nest(nside, loc2hpd_uv(nside, ang2loc(ang), u, v)); +} + + +int64_t vec2ring_uv(int64_t nside, t_vec vec, double* u, double* v) { + return hpd2ring(nside, loc2hpd_uv(nside, vec2loc(vec), u, v)); +} + + +int64_t vec2nest_uv(int64_t nside, t_vec vec, double* u, double* v) { + return hpd2nest(nside, loc2hpd_uv(nside, vec2loc(vec), u, v)); +} + + +t_ang ring2ang_uv(int64_t nside, int64_t ipix, double u, double v) { + return loc2ang(hpd2loc_uv(nside, ring2hpd(nside, ipix), u, v)); +} + + +t_ang nest2ang_uv(int64_t nside, int64_t ipix, double u, double v) { + return loc2ang(hpd2loc_uv(nside, nest2hpd(nside, ipix), u, v)); +} + + +t_vec ring2vec_uv(int64_t nside, int64_t ipix, double u, double v) { + return loc2vec(hpd2loc_uv(nside, ring2hpd(nside, ipix), u, v)); +} + + +t_vec nest2vec_uv(int64_t nside, int64_t ipix, double u, double v) { + return loc2vec(hpd2loc_uv(nside, nest2hpd(nside, ipix), u, v)); +} diff --git a/src/healpix.h b/src/healpix.h new file mode 100644 index 0000000..42d7266 --- /dev/null +++ b/src/healpix.h @@ -0,0 +1,168 @@ +// license: BSD-3-Clause + +#ifndef HEALPIX_H_ +#define HEALPIX_H_ + +#include + +#ifdef __cplusplus +extern "C" { +#endif + +/* This code is written in C99; you may have to adjust your compiler flags. */ + + +// Continuous coordinate systems +// +// theta is the co-latitude in radians (0 at the north Pole, increasing to pi +// at the south pole) +// +// phi is the azimuth in radians +// +// admissible values for theta: 0 <= theta <= pi +// +// admissible values for phi: in principle unconstrained, but best accuracy is +// obtained for -2pi <= phi <= 2pi + + +// A structure describing a location on the sphere. +typedef struct { + double theta; + double phi; +} t_ang; + + +// A structure describing a 3-vector with coordinates x, y and z. +typedef struct { + double x; + double y; + double z; +} t_vec; + + +// Returns a normalized 3-vector pointing in the same direction as ang. +t_vec ang2vec(t_ang ang); + + +// Returns a t_ang describing the same direction as the 3-vector vec. +// vec need not be normalized +t_ang vec2ang(t_vec vec); + + +/* Discrete coordinate systems */ + +/* Admissible values for nside parameters: + any integer power of 2 with 1 <= nside <= 1<<29 + + Admissible values for pixel indices: + 0 <= idx < 12*nside*nside */ + +extern const int64_t NSIDE_MAX; // 1<<29 + +/*! Returns the RING pixel index of pixel \a ipnest at resolution \a nside. + On error, returns -1. */ +int64_t nest2ring(int64_t nside, int64_t ipnest); +/*! Returns the NEST pixel index of pixel \a ipring at resolution \a nside. + On error, returns -1. */ +int64_t ring2nest(int64_t nside, int64_t ipring); + + +// Conversions between continuous and discrete coordinate systems + + +// Returns the pixel number in NEST scheme at resolution nside, which contains +// the position ang. +int64_t ang2nest(int64_t nside, t_ang ang); + + +// Returns the pixel number in RING scheme at resolution nside, which contains +// the position ang. +int64_t ang2ring(int64_t nside, t_ang ang); + + +// Returns a t_ang corresponding to the angular position of the center of pixel +// ipix in NEST scheme at resolution nside. +t_ang nest2ang(int64_t nside, int64_t ipix); + + +// Returns a t_ang corresponding to the angular position of the center of pixel +// ipix in RING scheme at resolution nside. +t_ang ring2ang(int64_t nside, int64_t ipix); + + +// Returns the pixel number in NEST scheme at resolution nside, which contains +// the direction described by the 3-vector vec. +int64_t vec2nest(int64_t nside, t_vec vec); + + +// Returns the pixel number in RING scheme at resolution nside, which contains +// the direction described by the 3-vector vec. +int64_t vec2ring(int64_t nside, t_vec vec); + + +// Returns a normalized 3-vector pointing in the direction of the center of +// pixel ipix in NEST scheme at resolution nside. +t_vec nest2vec(int64_t nside, int64_t ipix); + + +// Returns a normalized 3-vector pointing in the direction of the center of +// pixel ipix in RING scheme at resolution nside. +t_vec ring2vec(int64_t nside, int64_t ipix); + + +// Miscellaneous utility routines + + +// Returns 12*nside*nside. +int64_t nside2npix(int64_t nside); + + +// Returns sqrt(npix/12) if this is an integer number, otherwise -1. +int64_t npix2nside(int64_t npix); + + +// Returns the angle (in radians) between the vectors v1 and v2. +// The result is accurate even for angles close to 0 and pi. +double vec_angle(t_vec v1, t_vec v2); + + +// conversions with sub-pixel positions + + +// Variant of ang2nest that also returns the sub-pixel position in u and v. +int64_t ang2nest_uv(int64_t nside, t_ang ang, double* u, double* v); + + +// Variant of ang2ring that also returns the sub-pixel position in u and v. +int64_t ang2ring_uv(int64_t nside, t_ang ang, double* u, double* v); + + +// Variant of nest2ang that also takes the sub-pixel position in u and v. +t_ang nest2ang_uv(int64_t nside, int64_t ipix, double u, double v); + + +// Variant of ring2ang that also takes the sub-pixel position in u and v. +t_ang ring2ang_uv(int64_t nside, int64_t ipix, double u, double v); + + +// Variant of vec2nest that also returns the sub-pixel position in u and v. +int64_t vec2nest_uv(int64_t nside, t_vec vec, double* u, double* v); + + +// Variant of vec2ring that also returns the sub-pixel position in u and v. +int64_t vec2ring_uv(int64_t nside, t_vec vec, double* u, double* v); + + +// Variant of nest2vec that also takes the sub-pixel position in u and v. +t_vec nest2vec_uv(int64_t nside, int64_t ipix, double u, double v); + + +// Variant of ring2vec that also takes the sub-pixel position in u and v. +t_vec ring2vec_uv(int64_t nside, int64_t ipix, double u, double v); + + +#ifdef __cplusplus +} // extern "C" +#endif + +#endif // HEALPIX_H_ diff --git a/tools/memory_overhead.py b/tools/memory_overhead.py new file mode 100644 index 0000000..eafe1ae --- /dev/null +++ b/tools/memory_overhead.py @@ -0,0 +1,116 @@ +import argparse +import tracemalloc +import numpy as np +import healpix + + +TEST_FUNCTIONS = [ + 'ang2vec', + 'vec2ang', + 'ang2pix', + 'pix2ang', + 'vec2pix', + 'pix2vec', +] + + +def get_nbytes(obj): + if isinstance(obj, tuple): + return sum(map(get_nbytes, obj)) + else: + return obj.nbytes + + +def params_add(params, *args, **kwargs): + return [((*a, *args), {**kw, **kwargs}) for a, kw in params] + + +def params_mul(params, *args, **kwargs): + if args: + params = [((*a, *a_), kw) for a, kw in params for a_ in args] + if kwargs: + kwargs = [{k: v} for k in kwargs for v in kwargs[k]] + params = [(a, {**kw, **kw_}) for a, kw in params for kw_ in kwargs] + return params + + +def printargs(*args, **kwargs): + def fmt(x): + return '[...]' if np.ndim(x) > 0 else repr(x) + + args_ = ', '.join(fmt(a) for a in args) + kwargs_ = ', '.join(f'{k}={fmt(v)}' for k, v in kwargs.items()) + return ', '.join(filter(None, [args_, kwargs_])) + + +if __name__ == '__main__': + parser = argparse.ArgumentParser() + parser.add_argument('-N', '--nside', type=int, default=1024, + help='NSIDE parameter, sets array sizes') + parser.add_argument('--healpy', action='store_true', + help='check memory overhead of healpy functions') + args = parser.parse_args() + + nside = args.nside + use_healpy = args.healpy + + size = 12*nside**2 + + if use_healpy: + healpix = __import__('healpy') + + print('') + print('using %s functions with NSIDE=%d' % ('healpix' if not use_healpy else 'healpy', nside)) + print('') + print(f'{"Function":60} {"Overhead":8}') + print(f'{"-"*60:}--{"-"*8:}') + + for name in TEST_FUNCTIONS: + + base_params = [((), {})] + + if 'ang' in name: + base_params = params_mul(base_params, lonlat=[False, True]) + + for base_args, base_kwargs in base_params: + params = [(base_args, base_kwargs)] + if 'pix' in name: + params = params_add(params, nside) + params = params_mul(params, nest=[False, True]) + if name.startswith('pix2'): + ipix = np.arange(size) + params = params_add(params, ipix) + elif name.startswith('ang2'): + if base_kwargs['lonlat']: + lon = np.random.uniform(0, 360, size=size) + lat = np.random.uniform(-90, 90, size=size) + params = params_add(params, lon, lat) + else: + theta = np.random.uniform(0, np.pi, size=size) + phi = np.random.uniform(0, 2*np.pi, size=size) + params = params_add(params, theta, phi) + elif name.startswith('vec2'): + if use_healpy and name == 'vec2ang': + xyz = np.random.randn(size, 3) + params = params_add(params, xyz) + else: + x, y, z = np.random.randn(3, size) + params = params_add(params, x, y, z) + + func = getattr(healpix, name) + + for args, kwargs in params: + label = f'{name}({printargs(*args, **kwargs)})' + try: + tracemalloc.start() + result = func(*args, **kwargs) + current, peak = tracemalloc.get_traced_memory() + tracemalloc.stop() + except Exception as e: + raise RuntimeError(label) from e + output = get_nbytes(result) + assert current - output < 1_000_000 + overhead = peak/current - 1 + print(f'{label:60} {overhead:-8.0%}') + + print('') diff --git a/tools/nside_precision.py b/tools/nside_precision.py new file mode 100644 index 0000000..4d00710 --- /dev/null +++ b/tools/nside_precision.py @@ -0,0 +1,18 @@ +# test that numerical operations retain precision for high nside + +import healpix + +for k in range(30): + nside = 1<