From 87bf0c1990825f51ac4fa5a289461b50608ddcd0 Mon Sep 17 00:00:00 2001 From: Leo Singer Date: Tue, 30 Oct 2018 17:36:38 -0400 Subject: [PATCH 1/6] Rewrite core in C using Numpy ufuncs As a result, broadcasting across nside and ipix is supported automatically. Fixes #72, closes #75. --- astropy_healpix/_core.c | 368 ++++++++++++++++ astropy_healpix/core.py | 156 +++---- astropy_healpix/core_cython.pxd | 44 -- astropy_healpix/core_cython.pyx | 492 ---------------------- astropy_healpix/setup_package.py | 11 +- astropy_healpix/tests/test_core.py | 18 +- astropy_healpix/tests/test_core_cython.py | 129 ------ astropy_healpix/tests/test_healpy.py | 6 +- cextern/astrometry.net/healpix-utils.h | 2 +- cextern/lalsuite/six.h | 86 ++++ cextern/numpy/ieee754.h | 101 +++++ 11 files changed, 640 insertions(+), 773 deletions(-) create mode 100644 astropy_healpix/_core.c delete mode 100644 astropy_healpix/core_cython.pxd delete mode 100644 astropy_healpix/core_cython.pyx delete mode 100644 astropy_healpix/tests/test_core_cython.py create mode 100644 cextern/lalsuite/six.h create mode 100644 cextern/numpy/ieee754.h diff --git a/astropy_healpix/_core.c b/astropy_healpix/_core.c new file mode 100644 index 0000000..ccb0090 --- /dev/null +++ b/astropy_healpix/_core.c @@ -0,0 +1,368 @@ +#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION + +#include +#include +#include +#include "healpix.h" +#include "healpix-utils.h" +#include "interpolation.h" + +/* FIXME: Remove this once we drop Python 2 support. */ +#include "six.h" + +/* FIXME: We need npy_set_floatstatus_invalid(), but unlike most of the Numpy + * C API it is only available on some platforms if you explicitly link against + * Numpy, which is not typically done for building C extensions. This bundled + * header file contains a static definition of _npy_set_floatstatus_invalid(). + */ +#include "ieee754.h" + + +typedef struct { + int64_t (*order_to_xy)(int64_t, int); + int64_t (*xy_to_order)(int64_t, int); +} order_funcs; + + +static int64_t nside2npix(int nside) +{ + return 12 * ((int64_t) nside) * ((int64_t) nside); +} + + +static int pixel_nside_valid(int64_t pixel, int nside) +{ + return pixel >= 0 && pixel < nside2npix(nside); +} + + +static void healpix_to_lonlat_loop( + char **args, npy_intp *dimensions, npy_intp *steps, void *data) +{ + order_funcs *funcs = data; + npy_intp i, n = dimensions[0]; + + for (i = 0; i < n; i ++) + { + int64_t pixel = *(int64_t *) &args[0][i * steps[0]]; + int nside = *(int *) &args[1][i * steps[1]]; + double dx = *(double *) &args[2][i * steps[2]]; + double dy = *(double *) &args[3][i * steps[3]]; + double *lon = (double *) &args[4][i * steps[4]]; + double *lat = (double *) &args[5][i * steps[5]]; + int64_t xy = -1; + + if (pixel_nside_valid(pixel, nside)) + xy = funcs->order_to_xy(pixel, nside); + + if (xy >= 0) + healpixl_to_radec(xy, nside, dx, dy, lon, lat); + else + { + *lon = *lat = NPY_NAN; + _npy_set_floatstatus_invalid(); + } + } +} + + +static void lonlat_to_healpix_loop( + char **args, npy_intp *dimensions, npy_intp *steps, void *data) +{ + order_funcs *funcs = data; + npy_intp i, n = dimensions[0]; + + for (i = 0; i < n; i ++) + { + double lon = *(double *) &args[0][i * steps[0]]; + double lat = *(double *) &args[1][i * steps[1]]; + int nside = *(int *) &args[2][i * steps[2]]; + int64_t *pixel = (int64_t *) &args[3][i * steps[3]]; + double *dx = (double *) &args[4][i * steps[4]]; + double *dy = (double *) &args[5][i * steps[5]]; + int64_t xy = -1; + + xy = radec_to_healpixlf(lon, lat, nside, dx, dy); + if (xy >= 0) + *pixel = funcs->xy_to_order(xy, nside); + else { + *pixel = -1; + *dx = *dy = NPY_NAN; + _npy_set_floatstatus_invalid(); + } + } +} + + +static void nested_to_ring_loop( + char **args, npy_intp *dimensions, npy_intp *steps, void *data) +{ + npy_intp i, n = dimensions[0]; + + for (i = 0; i < n; i ++) + { + int64_t nested = *(int64_t *) &args[0][i * steps[0]]; + int nside = *(int *) &args[1][i * steps[1]]; + int64_t *ring = (int64_t *) &args[2][i * steps[2]]; + int64_t xy = -1; + + if (pixel_nside_valid(nested, nside)) + xy = healpixl_nested_to_xy(nested, nside); + if (xy >= 0) + *ring = healpixl_xy_to_ring(xy, nside); + else { + *ring = -1; + _npy_set_floatstatus_invalid(); + } + } +} + + +static void ring_to_nested_loop( + char **args, npy_intp *dimensions, npy_intp *steps, void *data) +{ + npy_intp i, n = dimensions[0]; + + for (i = 0; i < n; i ++) + { + int64_t ring = *(int64_t *) &args[0][i * steps[0]]; + int nside = *(int *) &args[1][i * steps[1]]; + int64_t *nested = (int64_t *) &args[2][i * steps[2]]; + int64_t xy = -1; + + if (pixel_nside_valid(ring, nside)) + xy = healpixl_ring_to_xy(ring, nside); + if (xy >= 0) + *nested = healpixl_xy_to_nested(xy, nside); + else { + *nested = -1; + _npy_set_floatstatus_invalid(); + } + } +} + + +static void bilinear_interpolation_weights_loop( + char **args, npy_intp *dimensions, npy_intp *steps, void *data) +{ + npy_intp i, n = dimensions[0]; + + for (i = 0; i < n; i ++) + { + double lon = *(double *) &args[0][i * steps[0]]; + double lat = *(double *) &args[1][i * steps[1]]; + int nside = *(int *) &args[2][i * steps[2]]; + int64_t indices[4]; + double weights[4]; + int j; + + interpolate_weights(lon, lat, indices, weights, nside); + for (j = 0; j < 4; j ++) + { + *(int64_t *) &args[3 + j][i * steps[3 + j]] = indices[j]; + *(double *) &args[7 + j][i * steps[7 + j]] = weights[j]; + } + } +} + + +static void neighbours_loop( + char **args, npy_intp *dimensions, npy_intp *steps, void *data) +{ + order_funcs *funcs = data; + npy_intp i, n = dimensions[0]; + + for (i = 0; i < n; i ++) + { + int64_t pixel = *(int64_t *) &args[0][i * steps[0]]; + int nside = *(int *) &args[1][i * steps[1]]; + int64_t neighbours[] = {-1, -1, -1, -1, -1, -1, -1, -1}; + int j; + int64_t xy = -1; + + if (pixel_nside_valid(pixel, nside)) + xy = funcs->order_to_xy(pixel, nside); + if (xy >= 0) + healpixl_get_neighbours(xy, neighbours, nside); + + for (j = 0; j < 8; j ++) + { + int k = 4 - j; + if (k < 0) + k += 8; + xy = neighbours[k]; + if (xy >= 0) + pixel = funcs->xy_to_order(xy, nside); + else { + pixel = -1; + _npy_set_floatstatus_invalid(); + } + *(int64_t *) &args[2 + j][i * steps[2 + j]] = pixel; + } + } +} + + +static PyObject *healpix_cone_search( + PyObject *self, PyObject *args, PyObject *kwargs) +{ + PyObject *result = NULL; + static const char *kws[] = {"lon", "lat", "radius", "nside", "order", NULL}; + double lon, lat, radius; + int nside; + char *order; + int64_t *indices, n_indices; + int64_t *result_data; + + if (!PyArg_ParseTupleAndKeywords( + args, kwargs, "dddis", kws, &lon, &lat, &radius, &nside, &order)) + return NULL; + + n_indices = healpix_rangesearch_radec_simple( + lon, lat, radius, nside, 0, &indices); + if (!indices) + { + PyErr_SetString( + PyExc_RuntimeError, "healpix_rangesearch_radec_simple failed"); + goto done; + } + + { + npy_intp dims[] = {n_indices}; + result = PyArray_SimpleNew(1, dims, NPY_INT64); + } + if (!result) + goto done; + + result_data = PyArray_DATA((PyArrayObject *) result); + + if (strcmp(order, "nested") == 0) + { + int i; + for (i = 0; i < n_indices; i ++) + result_data[i] = healpixl_xy_to_nested(indices[i], nside); + } else { + int i; + for (i = 0; i < n_indices; i ++) + result_data[i] = healpixl_xy_to_ring(indices[i], nside); + } + +done: + free(indices); + return result; +} + + +static PyMethodDef methods[] = { + {"healpix_cone_search", healpix_cone_search, METH_VARARGS | METH_KEYWORDS, NULL}, + {NULL, NULL, 0, NULL} +}; + + +static PyModuleDef moduledef = { + PyModuleDef_HEAD_INIT, + "_core", NULL, -1, methods +}; + +static const PyUFuncGenericFunction + healpix_to_lonlat_loops [] = {healpix_to_lonlat_loop}, + lonlat_to_healpix_loops [] = {lonlat_to_healpix_loop}, + nested_to_ring_loops [] = {nested_to_ring_loop}, + ring_to_nested_loops [] = {ring_to_nested_loop}, + bilinear_interpolation_weights_loops[] = {bilinear_interpolation_weights_loop}, + neighbours_loops [] = {neighbours_loop}; + +static const order_funcs + nested_order_funcs = {healpixl_nested_to_xy, healpixl_xy_to_nested}, + ring_order_funcs = {healpixl_ring_to_xy, healpixl_xy_to_ring}; + +static void *no_ufunc_data[] = {NULL}, + *nested_ufunc_data[] = {&nested_order_funcs}, + *ring_ufunc_data[] = {&ring_order_funcs}; + +static const char + healpix_to_lonlat_types[] = { + NPY_INT64, NPY_INT, NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE}, + lonlat_to_healpix_types[] = { + NPY_DOUBLE, NPY_DOUBLE, NPY_INT, NPY_INT64, NPY_DOUBLE, NPY_DOUBLE}, + healpix_to_healpix_types[] = { + NPY_INT64, NPY_INT, NPY_INT64}, + bilinear_interpolation_weights_types[] = { + NPY_DOUBLE, NPY_DOUBLE, NPY_INT, + NPY_INT64, NPY_INT64, NPY_INT64, NPY_INT64, + NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE}, + neighbours_types[] = { + NPY_INT64, NPY_INT, + NPY_INT64, NPY_INT64, NPY_INT64, NPY_INT64, + NPY_INT64, NPY_INT64, NPY_INT64, NPY_INT64}; + + +PyMODINIT_FUNC PyInit__core(void) +{ + PyObject *module; + + import_array(); + import_umath(); + + module = PyModule_Create(&moduledef); + + PyModule_AddObject( + module, "healpix_nested_to_lonlat", PyUFunc_FromFuncAndData( + healpix_to_lonlat_loops, nested_ufunc_data, + healpix_to_lonlat_types, 1, 4, 2, PyUFunc_None, + "healpix_nested_to_lonlat", NULL, 0)); + + PyModule_AddObject( + module, "healpix_ring_to_lonlat", PyUFunc_FromFuncAndData( + healpix_to_lonlat_loops, ring_ufunc_data, + healpix_to_lonlat_types, 1, 4, 2, PyUFunc_None, + "healpix_ring_to_lonlat", NULL, 0)); + + PyModule_AddObject( + module, "lonlat_to_healpix_nested", PyUFunc_FromFuncAndData( + lonlat_to_healpix_loops, nested_ufunc_data, + lonlat_to_healpix_types, 1, 3, 3, PyUFunc_None, + "lonlat_to_healpix_nested", NULL, 0)); + + PyModule_AddObject( + module, "lonlat_to_healpix_ring", PyUFunc_FromFuncAndData( + lonlat_to_healpix_loops, ring_ufunc_data, + lonlat_to_healpix_types, 1, 3, 3, PyUFunc_None, + "lonlat_to_healpix_ring", NULL, 0)); + + PyModule_AddObject( + module, "nested_to_ring", PyUFunc_FromFuncAndData( + nested_to_ring_loops, no_ufunc_data, + healpix_to_healpix_types, 1, 2, 1, PyUFunc_None, + "nested_to_ring", NULL, 0)); + + PyModule_AddObject( + module, "ring_to_nested", PyUFunc_FromFuncAndData( + ring_to_nested_loops, no_ufunc_data, + healpix_to_healpix_types, 1, 2, 1, PyUFunc_None, + "ring_to_nested", NULL, 0)); + + PyModule_AddObject( + module, "bilinear_interpolation_weights", PyUFunc_FromFuncAndData( + bilinear_interpolation_weights_loops, no_ufunc_data, + bilinear_interpolation_weights_types, 1, 3, 8, PyUFunc_None, + "bilinear_interpolation_weights", NULL, 0)); + + PyModule_AddObject( + module, "neighbours_nested", PyUFunc_FromFuncAndData( + neighbours_loops, nested_ufunc_data, + neighbours_types, 1, 2, 8, PyUFunc_None, + "neighbours_nested", NULL, 0)); + + PyModule_AddObject( + module, "neighbours_ring", PyUFunc_FromFuncAndData( + neighbours_loops, ring_ufunc_data, + neighbours_types, 1, 2, 8, PyUFunc_None, + "neighbours_ring", NULL, 0)); + + return module; +} + + +/* FIXME: Remove this once we drop Python 2 support. */ +SIX_COMPAT_MODULE(_core) diff --git a/astropy_healpix/core.py b/astropy_healpix/core.py index 9cd0a42..39f8bf0 100644 --- a/astropy_healpix/core.py +++ b/astropy_healpix/core.py @@ -9,8 +9,7 @@ from astropy import units as u from astropy.coordinates import Longitude, Latitude -from . import core_cython -from .core_cython import _validate_order +from . import _core __all__ = [ 'nside_to_pixel_area', @@ -44,6 +43,18 @@ def _restore_shape(*args, **kwargs): return np.asscalar(args[0]) +def _validate_order(order): + # We also support upper-case, to support directly the values + # ORDERING = {'RING', 'NESTED'} in FITS headers + # This is currently undocumented in the docstrings. + if order == 'nested' or order == 'NESTED': + return 'nested' + elif order == 'ring' or order == 'RING': + return 'ring' + else: + raise ValueError("order must be 'nested' or 'ring'") + + def _validate_healpix_index(label, healpix_index, nside): npix = nside_to_npix(nside) if np.any((healpix_index < 0) | (healpix_index > npix - 1)): @@ -51,6 +62,7 @@ def _validate_healpix_index(label, healpix_index, nside): def _validate_offset(label, offset): + offset = np.asarray(offset) if np.any((offset < 0) | (offset > 1)): raise ValueError('d{0} must be in the range [0:1]'.format(label)) @@ -341,7 +353,7 @@ def healpix_to_lonlat(healpix_index, nside, dx=None, dy=None, order='ring'): ---------- healpix_index : int or `~numpy.ndarray` HEALPix indices (as a scalar or array) - nside : int + nside : int or `~numpy.ndarray` Number of pixels along the side of each of the 12 top-level HEALPix tiles dx, dy : float or `~numpy.ndarray`, optional Offsets inside the HEALPix pixel, which must be in the range [0:1], @@ -357,42 +369,30 @@ def healpix_to_lonlat(healpix_index, nside, dx=None, dy=None, order='ring'): The latitude values """ - if (dx is None and dy is not None) or (dx is not None and dy is None): - raise ValueError('Either both or neither dx and dy must be specified') + _validate_nside(nside) - healpix_index = np.asarray(healpix_index, dtype=np.int64) + if _validate_order(order) == 'ring': + func = _core.healpix_ring_to_lonlat + else: # _validate_order(order) == 'nested' + func = _core.healpix_nested_to_lonlat - if dx is None and dy is not None: + if dx is None: dx = 0.5 - elif dx is not None and dy is None: - dy = 0.5 - - if dx is not None: - dx = np.asarray(dx, dtype=np.float) - dy = np.asarray(dy, dtype=np.float) + else: _validate_offset('x', dx) + if dy is None: + dy = 0.5 + else: _validate_offset('y', dy) - healpix_index, dx, dy = np.broadcast_arrays(healpix_index, dx, dy) - dx = dx.ravel() - dy = dy.ravel() - shape = healpix_index.shape - healpix_index = healpix_index.ravel() - nside = int(nside) + nside = np.asarray(nside, dtype=np.intc) - _validate_healpix_index('healpix_index', healpix_index, nside) - _validate_nside(nside) - order = _validate_order(order) - - if dx is None: - lon, lat = core_cython.healpix_to_lonlat(healpix_index, nside, order) - else: - lon, lat = core_cython.healpix_with_offset_to_lonlat(healpix_index, dx, dy, nside, order) + lon, lat = func(healpix_index, nside, dx, dy) lon = Longitude(lon, unit=u.rad, copy=False) lat = Latitude(lat, unit=u.rad, copy=False) - return _restore_shape(lon, lat, shape=shape) + return lon, lat def lonlat_to_healpix(lon, lat, nside, return_offsets=False, order='ring'): @@ -404,7 +404,7 @@ def lonlat_to_healpix(lon, lat, nside, return_offsets=False, order='ring'): lon, lat : :class:`~astropy.units.Quantity` The longitude and latitude values as :class:`~astropy.units.Quantity` instances with angle units. - nside : int + nside : int or `~numpy.ndarray` Number of pixels along the side of each of the 12 top-level HEALPix tiles order : { 'nested' | 'ring' } Order of HEALPix pixels @@ -422,26 +422,22 @@ def lonlat_to_healpix(lon, lat, nside, return_offsets=False, order='ring'): center of the HEALPix pixels """ - lon = lon.to(u.rad).value - lat = lat.to(u.rad).value + if _validate_order(order) == 'ring': + func = _core.lonlat_to_healpix_ring + else: # _validate_order(order) == 'nested' + func = _core.lonlat_to_healpix_nested - lon, lat = np.broadcast_arrays(lon, lat) + nside = np.asarray(nside, dtype=np.intc) - shape = np.shape(lon) + lon = lon.to_value(u.rad) + lat = lat.to_value(u.rad) - lon = lon.astype(float).ravel() - lat = lat.astype(float).ravel() - - nside = int(nside) - _validate_nside(nside) - order = _validate_order(order) + healpix_index, dx, dy = func(lon, lat, nside) if return_offsets: - healpix_index, dx, dy = core_cython.lonlat_to_healpix_with_offset(lon, lat, nside, order) - return _restore_shape(healpix_index, dx, dy, shape=shape) + return healpix_index, dx, dy else: - healpix_index = core_cython.lonlat_to_healpix(lon, lat, nside, order) - return _restore_shape(healpix_index, shape=shape) + return healpix_index def nested_to_ring(nested_index, nside): @@ -452,7 +448,7 @@ def nested_to_ring(nested_index, nside): ---------- nested_index : int or `~numpy.ndarray` Healpix index using the 'nested' ordering - nside : int + nside : int or `~numpy.ndarray` Number of pixels along the side of each of the 12 top-level HEALPix tiles Returns @@ -461,16 +457,9 @@ def nested_to_ring(nested_index, nside): Healpix index using the 'ring' ordering """ - nested_index = np.asarray(nested_index, dtype=np.int64) - shape = nested_index.shape - nested_index = nested_index.ravel() - nside = int(nside) + nside = np.asarray(nside, dtype=np.intc) - _validate_healpix_index('nested_index', nested_index, nside) - _validate_nside(nside) - - ring_index = core_cython.nested_to_ring(nested_index, nside) - return _restore_shape(ring_index, shape=shape) + return _core.nested_to_ring(nested_index, nside) def ring_to_nested(ring_index, nside): @@ -481,7 +470,7 @@ def ring_to_nested(ring_index, nside): ---------- ring_index : int or `~numpy.ndarray` Healpix index using the 'ring' ordering - nside : int + nside : int or `~numpy.ndarray` Number of pixels along the side of each of the 12 top-level HEALPix tiles Returns @@ -490,16 +479,9 @@ def ring_to_nested(ring_index, nside): Healpix index using the 'nested' ordering """ - ring_index = np.asarray(ring_index, dtype=np.int64) - shape = ring_index.shape - ring_index = ring_index.ravel() - nside = int(nside) - - _validate_healpix_index('ring_index', ring_index, nside) - _validate_nside(nside) + nside = np.asarray(nside, dtype=np.intc) - nested_index = core_cython.ring_to_nested(ring_index, nside) - return _restore_shape(nested_index, shape=shape) + return _core.ring_to_nested(ring_index, nside) def bilinear_interpolation_weights(lon, lat, nside, order='ring'): @@ -527,23 +509,19 @@ def bilinear_interpolation_weights(lon, lat, nside, order='ring'): interpolation """ - lon = lon.to(u.rad).value - lat = lat.to(u.rad).value + lon = lon.to_value(u.rad) + lat = lat.to_value(u.rad) - lon, lat = np.broadcast_arrays(lon, lat) + _validate_nside(nside) - shape = np.shape(lon) + nside = np.asarray(nside, dtype=np.intc) - lon = lon.astype(float).ravel() - lat = lat.astype(float).ravel() - nside = int(nside) + result = _core.bilinear_interpolation_weights(lon, lat, nside) + indices = np.stack(result[:4]) + weights = np.stack(result[4:]) - order = _validate_order(order) - _validate_nside(nside) - - indices, weights = core_cython.bilinear_interpolation_weights(lon, lat, nside, order) - indices = _restore_shape(indices, shape=(4,) + shape) - weights = _restore_shape(weights, shape=(4,) + shape) + if _validate_order(order) == 'nested': + indices = ring_to_nested(indices, nside) return indices, weights @@ -602,19 +580,16 @@ def neighbours(healpix_index, nside, order='ring'): ``neigh`` has shape (8, 2, 3). """ - healpix_index = np.asarray(healpix_index, dtype=np.int64) - nside = int(nside) - - shape = np.shape(healpix_index) + _validate_nside(nside) - healpix_index = healpix_index.ravel() + nside = np.asarray(nside, dtype=np.intc) - _validate_healpix_index('healpix_index', healpix_index, nside) - _validate_nside(nside) - order = _validate_order(order) + if _validate_order(order) == 'ring': + func = _core.neighbours_ring + else: # _validate_order(order) == 'nested' + func = _core.neighbours_nested - neigh = core_cython.neighbours(healpix_index, nside, order) - return _restore_shape(neigh, shape=(8,) + shape) + return np.stack(func(healpix_index, nside)) def healpix_cone_search(lon, lat, radius, nside, order='ring'): @@ -643,15 +618,14 @@ def healpix_cone_search(lon, lat, radius, nside, order='ring'): 1-D array with all the matching HEALPix pixel indices. """ - lon = float(lon.to(u.deg).value) - lat = float(lat.to(u.deg).value) - radius = float(radius.to(u.deg).value) - nside = int(nside) + lon = lon.to_value(u.deg) + lat = lat.to_value(u.deg) + radius = radius.to_value(u.deg) _validate_nside(nside) order = _validate_order(order) - return core_cython.healpix_cone_search(lon, lat, radius, nside, order) + return _core.healpix_cone_search(lon, lat, radius, nside, order) def boundaries_lonlat(healpix_index, step, nside, order='ring'): diff --git a/astropy_healpix/core_cython.pxd b/astropy_healpix/core_cython.pxd deleted file mode 100644 index bea16e5..0000000 --- a/astropy_healpix/core_cython.pxd +++ /dev/null @@ -1,44 +0,0 @@ -# Licensed under a 3-clause BSD style license - see LICENSE.rst - -# This file is needed in order to be able to cimport functions into other Cython -# files - -from numpy cimport int64_t - -cdef extern from "healpix.h": - - # Converts a healpix index from the XY scheme to the RING scheme. - int64_t healpixl_xy_to_ring(int64_t hp, int Nside) nogil - - # Converts a healpix index from the RING scheme to the XY scheme. - int64_t healpixl_ring_to_xy(int64_t ring_index, int Nside) nogil - - # Converts a healpix index from the XY scheme to the NESTED scheme. - int64_t healpixl_xy_to_nested(int64_t hp, int Nside) nogil - - # Converts a healpix index from the NESTED scheme to the XY scheme. - int64_t healpixl_nested_to_xy(int64_t nested_index, int Nside) nogil - - # Converts a healpix index, plus fractional offsets (dx,dy), into (x,y,z) - # coordinates on the unit sphere. (dx,dy) must be in [0, 1]. (0.5, 0.5) - # is the center of the healpix. (0,0) is the southernmost corner, (1,1) is - # the northernmost corner, (1,0) is the easternmost, and (0,1) the - # westernmost. - void healpixl_to_radec(int64_t hp, int Nside, double dx, double dy, - double* ra, double* dec) nogil - - # Converts (RA, DEC) coordinates (in radians) to healpix XY index - int64_t radec_to_healpixlf(double ra, double dec, int Nside, double* dx, double* dy) nogil - - # Finds the healpixes neighbouring the given healpix, placing them in the - # array "neighbour". Returns the number of neighbours. You must ensure - # that "neighbour" has 8 elements. - # - # Healpixes in the interior of a large healpix will have eight neighbours; - # pixels near the edges can have fewer. - int healpixl_get_neighbours(int64_t hp, int64_t* neighbours, int Nside) nogil - - int healpix_rangesearch_radec_simple(double ra, double dec, double radius, int Nside, int approx, int64_t** indices) nogil - - void interpolate_weights(double lon, double lat, int64_t *ring_indices, - double *weights, int Nside) nogil diff --git a/astropy_healpix/core_cython.pyx b/astropy_healpix/core_cython.pyx deleted file mode 100644 index 5c77b5f..0000000 --- a/astropy_healpix/core_cython.pyx +++ /dev/null @@ -1,492 +0,0 @@ -# Licensed under a 3-clause BSD style license - see LICENSE.rst -""" -This module contains vectorized Cython functions for common HEALPix operations. -Since they are written in Cython rather than Python, their input types are -strict and the functions will fail if the incorrect types are passed in. -""" - -import numpy as np -cimport numpy as np -import cython -from libc.stdlib cimport abort, malloc, free -from libc.math cimport sin, cos, sqrt - -ctypedef np.intp_t intp_t -ctypedef np.double_t double_t - -npy_double = np.double -npy_intp = np.intp -npy_int64 = np.int64 - - -def _validate_order(str order): - # We also support upper-case, to support directly the values - # ORDERING = {'RING', 'NESTED'} in FITS headers - # This is currently undocumented in the docstrings. - if order == 'nested' or order == 'NESTED': - return 'nested' - elif order == 'ring' or order == 'RING': - return 'ring' - else: - raise ValueError("order must be 'nested' or 'ring'") - - -@cython.boundscheck(False) -def healpix_to_lonlat(np.ndarray[int64_t, ndim=1, mode="c"] healpix_index, - int nside, str order): - """ - Convert HEALPix indices to longitudes/latitudes. - - This returns the longitudes/latitudes of the center of the HEALPix pixels. - If you also want to provide relative offsets inside the pixels, see - :func:`healpix_with_offset_to_lonlat`. - - Parameters - ---------- - healpix_index : `~numpy.ndarray` - 1-D array of HEALPix indices - nside : int - Number of pixels along the side of each of the 12 top-level HEALPix tiles - order : { 'nested' | 'ring' } - Order of HEALPix pixels - - Returns - ------- - lon, lat : `~numpy.ndarray` - 1-D arrays of longitude and latitude in radians - """ - - cdef intp_t n = healpix_index.shape[0] - cdef intp_t i - cdef int64_t xy_index - cdef double dx, dy; - cdef np.ndarray[double_t, ndim=1, mode="c"] lon = np.zeros(n, dtype=npy_double) - cdef np.ndarray[double_t, ndim=1, mode="c"] lat = np.zeros(n, dtype=npy_double) - - dx = 0.5 - dy = 0.5 - - order = _validate_order(order) - - if order == 'nested': - for i in range(n): - xy_index = healpixl_nested_to_xy(healpix_index[i], nside) - healpixl_to_radec(xy_index, nside, dx, dy, &lon[i], &lat[i]) - elif order == 'ring': - for i in range(n): - xy_index = healpixl_ring_to_xy(healpix_index[i], nside) - healpixl_to_radec(xy_index, nside, dx, dy, &lon[i], &lat[i]) - - return lon, lat - - -@cython.boundscheck(False) -def healpix_with_offset_to_lonlat(np.ndarray[int64_t, ndim=1, mode="c"] healpix_index, - np.ndarray[double_t, ndim=1, mode="c"] dx, - np.ndarray[double_t, ndim=1, mode="c"] dy, - int nside, str order): - """ - Convert HEALPix indices to longitudes/latitudes - - This function takes relative offsets in x and y inside the HEALPix pixels. - If you are only interested in the centers of the pixels, see - `healpix_to_lonlat`. - - Parameters - ---------- - healpix_index : `~numpy.ndarray` - 1-D array of HEALPix indices - dx, dy : `~numpy.ndarray` - 1-D arrays of offsets inside the HEALPix pixel, which must be in the - range [0:1] (0.5 is the center of the HEALPix pixels) - nside : int - Number of pixels along the side of each of the 12 top-level HEALPix tiles - order : { 'nested' | 'ring' } - Order of HEALPix pixels - - - Returns - ------- - lon, lat : `~numpy.ndarray` - 1-D arrays of longitude and latitude in radians - """ - - cdef intp_t n = healpix_index.shape[0] - cdef intp_t i - cdef int64_t xy_index - cdef np.ndarray[double_t, ndim=1, mode="c"] lon = np.zeros(n, dtype=npy_double) - cdef np.ndarray[double_t, ndim=1, mode="c"] lat = np.zeros(n, dtype=npy_double) - - order = _validate_order(order) - - if order == 'nested': - for i in range(n): - xy_index = healpixl_nested_to_xy(healpix_index[i], nside) - healpixl_to_radec(xy_index, nside, dx[i], dy[i], &lon[i], &lat[i]) - elif order == 'ring': - for i in range(n): - xy_index = healpixl_ring_to_xy(healpix_index[i], nside) - healpixl_to_radec(xy_index, nside, dx[i], dy[i], &lon[i], &lat[i]) - - return lon, lat - - -@cython.boundscheck(False) -def lonlat_to_healpix(np.ndarray[double_t, ndim=1, mode="c"] lon, - np.ndarray[double_t, ndim=1, mode="c"] lat, - int nside, str order): - """ - Convert longitudes/latitudes to HEALPix indices - - This returns only the HEALPix indices. If you also want to get relative - offsets inside the pixels, see :func:`lonlat_to_healpix_with_offset`. - - Parameters - ---------- - lon, lat : `~numpy.ndarray` - 1-D arrays of longitude and latitude in radians - nside : int - Number of pixels along the side of each of the 12 top-level HEALPix tiles - order : { 'nested' | 'ring' } - Order of HEALPix pixels - - - Returns - ------- - healpix_index : `~numpy.ndarray` - 1-D array of HEALPix indices - """ - - cdef intp_t n = lon.shape[0] - cdef intp_t i - cdef int64_t xy_index - cdef double dx, dy; - cdef np.ndarray[int64_t, ndim=1, mode="c"] healpix_index = np.zeros(n, dtype=npy_int64) - - order = _validate_order(order) - - if order == 'nested': - for i in range(n): - xy_index = radec_to_healpixlf(lon[i], lat[i], nside, &dx, &dy) - healpix_index[i] = healpixl_xy_to_nested(xy_index, nside) - elif order == 'ring': - for i in range(n): - xy_index = radec_to_healpixlf(lon[i], lat[i], nside, &dx, &dy) - healpix_index[i] = healpixl_xy_to_ring(xy_index, nside) - - return healpix_index - - -@cython.boundscheck(False) -def lonlat_to_healpix_with_offset(np.ndarray[double_t, ndim=1, mode="c"] lon, - np.ndarray[double_t, ndim=1, mode="c"] lat, - int nside, str order): - """ - Convert longitudes/latitudes to healpix indices - - This returns the HEALPix indices and relative offsets inside the pixels. If - you want only the HEALPix indices, see :func:`lonlat_to_healpix`. - - Parameters - ---------- - lon, lat : `~numpy.ndarray` - 1-D arrays of longitude and latitude in radians - nside : int - Number of pixels along the side of each of the 12 top-level HEALPix tiles - order : { 'nested' | 'ring' } - Order of HEALPix pixels - - Returns - ------- - healpix_index : `~numpy.ndarray` - 1-D array of HEALPix indices - dx, dy : `~numpy.ndarray` - 1-D arrays of offsets inside the HEALPix pixel in the range [0:1] (0.5 - is the center of the HEALPix pixels) - """ - - cdef intp_t n = lon.shape[0] - cdef intp_t i - cdef int64_t xy_index - cdef np.ndarray[int64_t, ndim=1, mode="c"] healpix_index = np.zeros(n, dtype=npy_int64) - cdef np.ndarray[double_t, ndim=1, mode="c"] dx = np.zeros(n, dtype=npy_double) - cdef np.ndarray[double_t, ndim=1, mode="c"] dy = np.zeros(n, dtype=npy_double) - - order = _validate_order(order) - - if order == 'nested': - for i in range(n): - xy_index = radec_to_healpixlf(lon[i], lat[i], nside, &dx[i], &dy[i]) - healpix_index[i] = healpixl_xy_to_nested(xy_index, nside) - elif order == 'ring': - for i in range(n): - xy_index = radec_to_healpixlf(lon[i], lat[i], nside, &dx[i], &dy[i]) - healpix_index[i] = healpixl_xy_to_ring(xy_index, nside) - - return healpix_index, dx, dy - - -@cython.boundscheck(False) -def nested_to_ring(np.ndarray[int64_t, ndim=1, mode="c"] nested_index, int nside): - """ - Convert a HEALPix 'nested' index to a HEALPix 'ring' index - - Parameters - ---------- - nested_index : `~numpy.ndarray` - Healpix index using the 'nested' ordering - nside : int - Number of pixels along the side of each of the 12 top-level HEALPix tiles - - Returns - ------- - ring_index : `~numpy.ndarray` - Healpix index using the 'ring' ordering - """ - - cdef intp_t n = nested_index.shape[0] - cdef intp_t i - cdef np.ndarray[int64_t, ndim=1, mode="c"] ring_index = np.zeros(n, dtype=npy_int64) - - for i in range(n): - ring_index[i] = healpixl_xy_to_ring(healpixl_nested_to_xy(nested_index[i], nside), nside) - - return ring_index - - -@cython.boundscheck(False) -def ring_to_nested(np.ndarray[int64_t, ndim=1, mode="c"] ring_index, int nside): - """ - Convert a HEALPix 'ring' index to a HEALPix 'nested' index - - Parameters - ---------- - ring_index : `~numpy.ndarray` - Healpix index using the 'ring' ordering - nside : int - Number of pixels along the side of each of the 12 top-level HEALPix tiles - - Returns - ------- - nested_index : `~numpy.ndarray` - Healpix index using the 'nested' ordering - """ - - cdef intp_t n = ring_index.shape[0] - cdef intp_t i - cdef np.ndarray[int64_t, ndim=1, mode="c"] nested_index = np.zeros(n, dtype=npy_int64) - - for i in range(n): - nested_index[i] = healpixl_xy_to_nested(healpixl_ring_to_xy(ring_index[i], nside), nside) - - return nested_index - - - -@cython.boundscheck(False) -def bilinear_interpolation_weights(np.ndarray[double_t, ndim=1, mode="c"] lon, - np.ndarray[double_t, ndim=1, mode="c"] lat, - int nside, str order): - """ - Get the four neighbours for each (lon, lat) position and the weight - associated with each one for bilinear interpolation. - - Parameters - ---------- - lon, lat : `~numpy.ndarray` - 1-D arrays of longitude and latitude in radians - nside : int - Number of pixels along the side of each of the 12 top-level HEALPix tiles - order : { 'nested' | 'ring' } - Order of HEALPix pixels - - Returns - ------- - indices : `~numpy.ndarray` - 2-D array with shape (4, N) giving the four indices to use for the - interpolation - weights : `~numpy.ndarray` - 2-D array with shape (4, N) giving the four weights to use for the - interpolation - """ - - cdef intp_t n = lon.shape[0] - cdef intp_t i, j - cdef np.ndarray[int64_t, ndim=2, mode="c"] indices = np.zeros((4, n), dtype=npy_int64) - cdef np.ndarray[double_t, ndim=2, mode="c"] weights = np.zeros((4, n), dtype=npy_double) - cdef int order_int - - cdef double *weights_indiv - cdef int64_t *indices_indiv - - order = _validate_order(order) - - if order == 'nested': - order_int = 0 - elif order == 'ring': - order_int = 1 - - with nogil: - - indices_indiv = malloc(sizeof(int64_t) * 4) - if indices_indiv == NULL: - abort() - - weights_indiv = malloc(sizeof(double) * 4) - if weights_indiv == NULL: - abort() - - for i in range(n): - interpolate_weights(lon[i], lat[i], indices_indiv, weights_indiv, nside) - for j in range(4): - if order_int == 0: - indices[j, i] = healpixl_xy_to_nested(healpixl_ring_to_xy(indices_indiv[j], nside), nside) - else: - indices[j, i] = indices_indiv[j] - weights[j, i] = weights_indiv[j] - - return indices, weights - - -@cython.boundscheck(False) -def neighbours(np.ndarray[int64_t, ndim=1, mode="c"] healpix_index, - int nside, str order): - """ - Find all the HEALPix pixels that are the neighbours of a HEALPix pixel - - Parameters - ---------- - healpix_index : `~numpy.ndarray` - 1-D array of HEALPix indices - nside : int - Number of pixels along the side of each of the 12 top-level HEALPix tiles - order : { 'nested' | 'ring' } - Order of HEALPix pixels - - Returns - ------- - neighbours : `~numpy.ndarray` - 2-D array with shape (8, N) giving the neighbours starting SW and - rotating clockwise. - """ - - cdef intp_t n = healpix_index.shape[0] - cdef intp_t i - cdef int64_t xy_index - cdef int j, k - cdef np.ndarray[int64_t, ndim=2, mode="c"] neighbours = np.zeros((8, n), dtype=npy_int64) - cdef int64_t * neighbours_indiv - cdef int order_int - - order = _validate_order(order) - - if order == 'nested': - order_int = 0 - elif order == 'ring': - order_int = 1 - - with nogil: - - neighbours_indiv = malloc(sizeof(int64_t) * 8) - if neighbours_indiv == NULL: - abort() - - # The neighbours above are ordered as follows: - # - # 3 2 1 - # 4 X 0 - # 5 6 7 - # - # but we want: - # - # 2 3 4 - # 1 X 5 - # 0 7 6 - # - # so we reorder these on-the-fly - - if order_int == 0: - - for i in range(n): - - xy_index = healpixl_nested_to_xy(healpix_index[i], nside) - healpixl_get_neighbours(xy_index, neighbours_indiv, nside) - - for j in range(8): - k = 4 - j - if k < 0: - k = k + 8 - if neighbours_indiv[k] < 0: - neighbours[j, i] = -1 - else: - neighbours[j, i] = healpixl_xy_to_nested(neighbours_indiv[k], nside) - - elif order_int == 1: - - for i in range(n): - - xy_index = healpixl_ring_to_xy(healpix_index[i], nside) - - healpixl_get_neighbours(xy_index, neighbours_indiv, nside) - - for j in range(8): - k = 4 - j - if k < 0: - k = k + 8 - if neighbours_indiv[k] < 0: - neighbours[j, i] = -1 - else: - neighbours[j, i] = healpixl_xy_to_ring(neighbours_indiv[k], nside) - - free(neighbours_indiv) - - return neighbours - - -@cython.boundscheck(False) -def healpix_cone_search(double lon, double lat, double radius, int nside, str order): - """ - Find all the HEALPix pixels within a given radius of a longitude/latitude. - - Note that this returns all pixels that overlap, including partially, with - the search cone. This function can only be used for a single lon/lat pair at - a time, since different calls to the function may result in a different - number of matches. - - Parameters - ---------- - lon, lat : float - The longitude and latitude to search around, in degrees - radius : float - The search radius, in degrees - nside : int - Number of pixels along the side of each of the 12 top-level HEALPix tiles - order : { 'nested' | 'ring' } - Order of HEALPix pixels - - Returns - ------- - healpix_index : `~numpy.ndarray` - 1-D array with all the matching HEALPix pixel indices. - """ - cdef intp_t i - cdef int64_t *indices - cdef int64_t n_indices - cdef int64_t index - - n_indices = healpix_rangesearch_radec_simple(lon, lat, radius, nside, 0, &indices) - - cdef np.ndarray[int64_t, ndim=1, mode="c"] result = np.zeros(n_indices, dtype=npy_int64) - - order = _validate_order(order) - - if order == 'nested': - for i in range(n_indices): - index = indices[i] - result[i] = healpixl_xy_to_nested(index, nside) - elif order == 'ring': - for i in range(n_indices): - index = indices[i] - result[i] = healpixl_xy_to_ring(index, nside) - - return result diff --git a/astropy_healpix/setup_package.py b/astropy_healpix/setup_package.py index cb73fb3..d33ea1a 100644 --- a/astropy_healpix/setup_package.py +++ b/astropy_healpix/setup_package.py @@ -16,6 +16,9 @@ C_DIR = os.path.join('cextern', 'astrometry.net') +C_DIRS = ['numpy', C_DIR, HEALPIX_ROOT, + os.path.join('cextern', 'lalsuite'), + os.path.join('cextern', 'numpy')] def get_extensions(): @@ -24,14 +27,12 @@ def get_extensions(): sources = [os.path.join(C_DIR, filename) for filename in C_FILES] sources.append(os.path.join(HEALPIX_ROOT, 'interpolation.c')) - sources.append(os.path.join(HEALPIX_ROOT, 'core_cython.pyx')) - - include_dirs = ['numpy', C_DIR, HEALPIX_ROOT] + sources.append(os.path.join(HEALPIX_ROOT, '_core.c')) extension = Extension( - name="astropy_healpix.core_cython", + name="astropy_healpix._core", sources=sources, - include_dirs=include_dirs, + include_dirs=C_DIRS, libraries=libraries, language="c", extra_compile_args=['-O2']) diff --git a/astropy_healpix/tests/test_core.py b/astropy_healpix/tests/test_core.py index e390a4e..8af5b0a 100644 --- a/astropy_healpix/tests/test_core.py +++ b/astropy_healpix/tests/test_core.py @@ -159,9 +159,11 @@ def test_healpix_to_lonlat_invalid(): dx = [0.1, 0.4, 0.9] dy = [0.4, 0.3, 0.2] - with pytest.raises(ValueError) as exc: + with pytest.warns(RuntimeWarning, match='invalid value'): lon, lat = healpix_to_lonlat([-1, 2, 3], 4) - assert exc.value.args[0] == 'healpix_index must be in the range [0:192]' + + with pytest.warns(RuntimeWarning, match='invalid value'): + lon, lat = healpix_to_lonlat([192, 2, 3], 4) with pytest.raises(ValueError) as exc: lon, lat = healpix_to_lonlat([1, 2, 3], 5) @@ -193,14 +195,14 @@ def test_healpix_to_lonlat_shape(): def test_lonlat_to_healpix_shape(): healpix_index = lonlat_to_healpix(2 * u.deg, 3 * u.deg, 8) - assert isinstance(healpix_index, integer_types) + assert np.can_cast(healpix_index, np.int64) lon, lat = np.ones((2, 4)) * u.deg, np.zeros((2, 4)) * u.deg healpix_index = lonlat_to_healpix(lon, lat, 8) assert healpix_index.shape == (2, 4) healpix_index, dx, dy = lonlat_to_healpix(2 * u.deg, 3 * u.deg, 8, return_offsets=True) - assert isinstance(healpix_index, integer_types) + assert np.can_cast(healpix_index, np.int64) assert isinstance(dx, float) assert isinstance(dy, float) @@ -214,7 +216,7 @@ def test_lonlat_to_healpix_shape(): @pytest.mark.parametrize('function', [nested_to_ring, ring_to_nested]) def test_nested_ring_shape(function): index = function(1, 8) - assert isinstance(index, integer_types) + assert np.can_cast(index, np.int64) index = function([[1, 2, 3], [2, 3, 4]], 8) assert index.shape == (2, 3) @@ -330,9 +332,11 @@ def test_neighbours(order): def test_neighbours_invalid(): - with pytest.raises(ValueError) as exc: + with pytest.warns(RuntimeWarning, match='invalid value'): neighbours([-1, 2, 3], 4) - assert exc.value.args[0] == 'healpix_index must be in the range [0:192]' + + with pytest.warns(RuntimeWarning, match='invalid value'): + neighbours([192, 2, 3], 4) with pytest.raises(ValueError) as exc: neighbours([1, 2, 3], 5) diff --git a/astropy_healpix/tests/test_core_cython.py b/astropy_healpix/tests/test_core_cython.py deleted file mode 100644 index 9977794..0000000 --- a/astropy_healpix/tests/test_core_cython.py +++ /dev/null @@ -1,129 +0,0 @@ -# Licensed under a 3-clause BSD style license - see LICENSE.rst - -from __future__ import absolute_import, print_function, division - -from itertools import product - -import pytest - -import numpy as np -from numpy.testing import assert_equal, assert_allclose - -from astropy import units as u -from astropy.coordinates.angle_utilities import angular_separation - -from ..core import nside_to_pixel_resolution -from .. import core_cython - -NSIDE_POWERS = range(0, 21) -ORDERS = ('nested', 'ring') - - -def get_test_indices(nside): - # For large number of pixels, we only compute a random subset of points - if nside > 2 ** 8: - try: - return np.random.randint(0, 12 * nside ** 2, 12 * 8 ** 2, dtype=np.int64) - except TypeError: # Numpy 1.9 and 1.10 - return (np.random.random(12 * 8 ** 2) * (12 * float(nside) ** 2)).astype(np.int64, copy=False) - else: - return np.arange(12 * nside ** 2, dtype=np.int64) - - -# NOTE: we use capfd in all tests here to make sure no errors/warnings are being -# raised by the C code. - -@pytest.mark.parametrize(('order', 'nside_power'), product(ORDERS, NSIDE_POWERS)) -def test_roundtrip_healpix_no_offsets(order, nside_power, capfd): - nside = 2 ** nside_power - index = get_test_indices(nside) - lon, lat = core_cython.healpix_to_lonlat(index, nside, order) - index_new = core_cython.lonlat_to_healpix(lon, lat, nside, order) - assert_equal(index, index_new) - out, err = capfd.readouterr() - assert out == "" and err == "" - - -@pytest.mark.parametrize(('order', 'nside_power'), product(ORDERS, NSIDE_POWERS)) -def test_roundtrip_healpix_with_offsets(order, nside_power, capfd): - nside = 2 ** nside_power - index = get_test_indices(nside) - dx = np.random.random(index.shape) - dy = np.random.random(index.shape) - lon, lat = core_cython.healpix_with_offset_to_lonlat(index, dx, dy, nside, order) - index_new, dx_new, dy_new = core_cython.lonlat_to_healpix_with_offset(lon, lat, nside, order) - assert_equal(index, index_new) - assert_allclose(dx, dx_new, atol=1e-8) - assert_allclose(dy, dy_new, atol=1e-8) - out, err = capfd.readouterr() - assert out == "" and err == "" - - -@pytest.mark.parametrize('nside_power', NSIDE_POWERS) -def test_roundtrip_nested_ring(nside_power, capfd): - nside = 2 ** nside_power - nested_index = get_test_indices(nside) - ring_index = core_cython.nested_to_ring(nested_index, nside) - nested_index_new = core_cython.ring_to_nested(ring_index, nside) - assert_equal(nested_index, nested_index_new) - if nside == 1: - assert np.all(nested_index == ring_index) - else: - assert not np.all(nested_index == ring_index) - out, err = capfd.readouterr() - assert out == "" and err == "" - - -@pytest.mark.parametrize(('order', 'nside_power'), product(ORDERS, NSIDE_POWERS)) -def test_neighbours(order, nside_power, capfd): - # This just makes sure things run, but doesn't check the validity of result - nside = 2 ** nside_power - index = get_test_indices(nside) - neighbours = core_cython.neighbours(index, nside, order) - assert np.all(neighbours >= -1) and np.all(neighbours < 12 * nside ** 2) - out, err = capfd.readouterr() - assert out == "" and err == "" - - -@pytest.mark.parametrize(('order', 'nside_power'), product(ORDERS, NSIDE_POWERS)) -def test_healpix_cone_search(order, nside_power, capfd): - nside = 2 ** nside_power - lon0, lat0 = 12., 40. - radius = nside_to_pixel_resolution(nside).to(u.degree).value * 10 - index_inside = core_cython.healpix_cone_search(lon0, lat0, radius, nside, order) - n_inside = len(index_inside) - dx = np.array([[0.0, 0.0, 1.0, 1.0]]) - dy = np.array([[0.0, 1.0, 1.0, 0.0]]) - dx = np.repeat(dx, n_inside, axis=0).ravel() - dy = np.repeat(dy, n_inside, axis=0).ravel() - index_inside = np.repeat(index_inside, 4).ravel() - lon, lat = core_cython.healpix_with_offset_to_lonlat(index_inside, dx, dy, nside, order) - lon, lat = lon.reshape((n_inside, 4)), lat.reshape((n_inside, 4)) - sep = angular_separation(lon0 * u.deg, lat0 * u.deg, lon * u.rad, lat * u.rad) - sep = sep.min(axis=1) - assert np.all(sep.to(u.degree).value < radius * 1.05) - out, err = capfd.readouterr() - assert out == "" and err == "" - - -# Regression tests for fixed bugs - -def test_regression_healpix_to_lonlat_sqrt(): - - # Regression test for a bug that caused the ring index decomposition to fail - # and return a negative longitude. - - index = np.array([9007199120523263], dtype=np.int64) - lon, lat = core_cython.healpix_to_lonlat(index, 2**26, order='ring') - assert_allclose(lon, 6.283185295476241, rtol=1e-14) - assert_allclose(lat, 0.729727669554970, rtol=1e-14) - - index = np.array([720575940916150240], dtype=np.int64) - lon, lat = core_cython.healpix_to_lonlat(index, 2**28, order='ring') - assert_allclose(lon, 6.283185122851909, rtol=1e-14) - assert_allclose(lat, -0.729727656226966, rtol=1e-14) - - index = np.array([180143985363255292], dtype=np.int64) - lon, lat = core_cython.healpix_to_lonlat(index, 2**27, order='ring') - assert_allclose(lon, 6.283185266217880, rtol=1e-14) - assert_allclose(lat, -0.729727656226966, rtol=1e-14) diff --git a/astropy_healpix/tests/test_healpy.py b/astropy_healpix/tests/test_healpy.py index 8afd164..758e165 100644 --- a/astropy_healpix/tests/test_healpy.py +++ b/astropy_healpix/tests/test_healpy.py @@ -4,8 +4,6 @@ from itertools import product -from .six import integer_types - import pytest import numpy as np @@ -80,7 +78,7 @@ def test_ang2pix(nside_pow, lon, lat, nest, lonlat): def test_ang2pix_shape(): ipix = hp_compat.ang2pix(8, 1., 2.) - assert isinstance(ipix, integer_types) + assert np.can_cast(ipix, np.int64) ipix = hp_compat.ang2pix(8, [[1., 2.], [3., 4.]], [[1., 2.], [3., 4.]]) assert ipix.shape == (2, 2) @@ -142,7 +140,7 @@ def test_pix2vec(nside_pow, frac, nest): def test_vec2pix_shape(): ipix = hp_compat.vec2pix(8, 1., 2., 3.) - assert isinstance(ipix, integer_types) + assert np.can_cast(ipix, np.int64) ipix = hp_compat.vec2pix(8, [[1., 2.], [3., 4.]], [[5., 6.], [7., 8.]], [[9., 10.], [11., 12.]]) assert ipix.shape == (2, 2) diff --git a/cextern/astrometry.net/healpix-utils.h b/cextern/astrometry.net/healpix-utils.h index a59aa71..f1ad43e 100644 --- a/cextern/astrometry.net/healpix-utils.h +++ b/cextern/astrometry.net/healpix-utils.h @@ -15,7 +15,7 @@ il* healpix_rangesearch_xyz(const double* xyz, double radius, int Nside, il* hps il* healpix_rangesearch_xyz_approx(const double* xyz, double radius, int Nside, il* hps); il* healpix_rangesearch_radec_approx(double ra, double dec, double radius, int Nside, il* hps); il* healpix_rangesearch_radec(double ra, double dec, double radius, int Nside, il* hps); -int healpix_rangesearch_radec_simple(double ra, double dec, double radius, int Nside, int **indices); +int64_t healpix_rangesearch_radec_simple(double ra, double dec, double radius, int Nside, int approx, int64_t **indices); /** Starting from a "seed" or list of "seeds" healpixes, grows a region diff --git a/cextern/lalsuite/six.h b/cextern/lalsuite/six.h new file mode 100644 index 0000000..a4176fe --- /dev/null +++ b/cextern/lalsuite/six.h @@ -0,0 +1,86 @@ +/* + * Copyright (C) 2016 Leo Singer + * + * Fake implementations of some parts of the Python C API functions to aid + * 2 to 3 porting. Named after the 'six' package that serves a similar + * purpose for pure Python code (https://pypi.python.org/pypi/six). + * + * + * Copied from https://github.com/lscsoft/lalsuite/blob/master/gnuscripts/six.h + * but re-released under a 3-clause BSD style license - see LICENSE.rst + */ + +#include + +#ifndef _SIX_H +#define _SIX_H + +#if PY_MAJOR_VERSION < 3 +#ifdef PyMODINIT_FUNC +#undef PyMODINIT_FUNC +#define PyMODINIT_FUNC PyObject* +#endif + +#define SIX_COMPAT_MODULE(name) _SIX_COMPAT_MODULE(name) +#define _SIX_COMPAT_MODULE(name) \ +void init##name(void); /* Silence -Wmissing-prototypes */ \ +void init##name(void) { \ + PyInit_##name(); \ +} + +typedef struct PyModuleDef { + int m_base; + const char *m_name; + const char *m_doc; + Py_ssize_t m_size; + PyMethodDef *m_methods; + inquiry m_reload; + traverseproc m_traverse; + inquiry m_clear; + freefunc m_free; +} PyModuleDef; + +#define PyModuleDef_HEAD_INIT 0 + +static inline PyObject *PyModule_Create(PyModuleDef *def) +{ + (void)PyModule_Create; /* Suppress unused function warning */ + + if (def->m_size != -1) + { + PyErr_SetString(PyExc_NotImplementedError, + "Python 2 does not support per-module state."); + return NULL; + } + return Py_InitModule3(def->m_name, def->m_methods, def->m_doc); +} + +static inline PyObject *PyModule_Create2(PyModuleDef *def, int module_api_version) +{ + (void)PyModule_Create2; /* Suppress unused function warning */ + + if (def->m_size != -1) + { + PyErr_SetString(PyExc_NotImplementedError, + "Python 2 does not support per-module state."); + return NULL; + } + return Py_InitModule4(def->m_name, def->m_methods, def->m_doc, + NULL, module_api_version); +} + +#ifdef NUMPY_IMPORT_ARRAY_RETVAL +#undef NUMPY_IMPORT_ARRAY_RETVAL +#endif +#define NUMPY_IMPORT_ARRAY_RETVAL NULL + +#ifdef NUMPY_IMPORT_UMATH_RETVAL +#undef NUMPY_IMPORT_UMATH_RETVAL +#endif +#define NUMPY_IMPORT_UMATH_RETVAL NULL + +#else /* PY_MAJOR_VERSION >= 3 */ +#define SIX_COMPAT_MODULE(name) +#endif + +#endif /* _SIX_H */ diff --git a/cextern/numpy/ieee754.h b/cextern/numpy/ieee754.h new file mode 100644 index 0000000..99fc1c1 --- /dev/null +++ b/cextern/numpy/ieee754.h @@ -0,0 +1,101 @@ +/* -*- c -*- */ +/* + * vim:syntax=c + * + * Low-level routines related to IEEE-754 format + * + * Adapted from https://github.com/numpy/numpy/blob/master/numpy/core/src/npymath/ieee754.c.src, + * removing all functions except for npy_set_floatstatus_invalid, which is + * made static and has is renamed to _npy_set_floatstatus_invalid. + * + * + * Copyright (c) 2005-2017, NumPy Developers. + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are + * met: + * + * * Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * + * * 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. + * + * * Neither the name of the NumPy Developers nor the names of any + * 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 + * OWNER 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. + */ +#include "numpy/utils.h" + + +#if (defined(__unix__) || defined(unix)) && !defined(USG) +#include +#endif + + +/* Solaris --------------------------------------------------------*/ +/* --------ignoring SunOS ieee_flags approach, someone else can +** deal with that! */ +#if defined(sun) || defined(__BSD__) || defined(__OpenBSD__) || \ + (defined(__FreeBSD__) && (__FreeBSD_version < 502114)) || \ + defined(__NetBSD__) +#include + +static void _npy_set_floatstatus_invalid(void) +{ + fpsetsticky(FP_X_INV); +} + +#elif defined(_AIX) +#include +#include + +static void _npy_set_floatstatus_invalid(void) +{ + fp_raise_xcp(FP_INVALID); +} + +#elif defined(_MSC_VER) || (defined(__osf__) && defined(__alpha)) + +/* + * By using a volatile floating point value, + * the compiler is forced to actually do the requested + * operations because of potential concurrency. + * + * We shouldn't write multiple values to a single + * global here, because that would cause + * a race condition. + */ +static volatile double _npy_floatstatus_x, _npy_floatstatus_inf; + +static void _npy_set_floatstatus_invalid(void) +{ + _npy_floatstatus_inf = NPY_INFINITY; + _npy_floatstatus_x = _npy_floatstatus_inf - NPY_INFINITY; +} + +#else +/* General GCC code, should work on most platforms */ +# include + +static void _npy_set_floatstatus_invalid(void) +{ + feraiseexcept(FE_INVALID); +} + +#endif From b8e37628fb7a95ba0d84a4dcfdf57a4018ee1ef9 Mon Sep 17 00:00:00 2001 From: Leo Singer Date: Wed, 31 Oct 2018 15:14:09 -0400 Subject: [PATCH 2/6] MSVC does not support inline keyword --- cextern/lalsuite/six.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cextern/lalsuite/six.h b/cextern/lalsuite/six.h index a4176fe..7e62c4f 100644 --- a/cextern/lalsuite/six.h +++ b/cextern/lalsuite/six.h @@ -42,7 +42,7 @@ typedef struct PyModuleDef { #define PyModuleDef_HEAD_INIT 0 -static inline PyObject *PyModule_Create(PyModuleDef *def) +static PyObject *PyModule_Create(PyModuleDef *def) { (void)PyModule_Create; /* Suppress unused function warning */ @@ -55,7 +55,7 @@ static inline PyObject *PyModule_Create(PyModuleDef *def) return Py_InitModule3(def->m_name, def->m_methods, def->m_doc); } -static inline PyObject *PyModule_Create2(PyModuleDef *def, int module_api_version) +static PyObject *PyModule_Create2(PyModuleDef *def, int module_api_version) { (void)PyModule_Create2; /* Suppress unused function warning */ From ac245b79f221007b0d76f454ec9935a95186b168 Mon Sep 17 00:00:00 2001 From: Leo Singer Date: Thu, 1 Nov 2018 18:06:05 -0400 Subject: [PATCH 3/6] Add comment for ufunc data --- astropy_healpix/_core.c | 21 +++++++++++++-------- 1 file changed, 13 insertions(+), 8 deletions(-) diff --git a/astropy_healpix/_core.c b/astropy_healpix/_core.c index ccb0090..52a62b0 100644 --- a/astropy_healpix/_core.c +++ b/astropy_healpix/_core.c @@ -18,11 +18,24 @@ #include "ieee754.h" +/* Data structure for storing function pointers for routines that are specific + * to the HEALPix ordering scheme. When we create the ufuncs using + * PyUFunc_FromFuncAndData, we will set them up to pass a pointer to this + * data structure through as the void *data parameter for the loop functions + * defined below. */ typedef struct { int64_t (*order_to_xy)(int64_t, int); int64_t (*xy_to_order)(int64_t, int); } order_funcs; +static const order_funcs + nested_order_funcs = {healpixl_nested_to_xy, healpixl_xy_to_nested}, + ring_order_funcs = {healpixl_ring_to_xy, healpixl_xy_to_ring}; + +static void *no_ufunc_data[] = {NULL}, + *nested_ufunc_data[] = {&nested_order_funcs}, + *ring_ufunc_data[] = {&ring_order_funcs}; + static int64_t nside2npix(int nside) { @@ -272,14 +285,6 @@ static const PyUFuncGenericFunction bilinear_interpolation_weights_loops[] = {bilinear_interpolation_weights_loop}, neighbours_loops [] = {neighbours_loop}; -static const order_funcs - nested_order_funcs = {healpixl_nested_to_xy, healpixl_xy_to_nested}, - ring_order_funcs = {healpixl_ring_to_xy, healpixl_xy_to_ring}; - -static void *no_ufunc_data[] = {NULL}, - *nested_ufunc_data[] = {&nested_order_funcs}, - *ring_ufunc_data[] = {&ring_order_funcs}; - static const char healpix_to_lonlat_types[] = { NPY_INT64, NPY_INT, NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE}, From 015edd38833e86a25a32a34c216f8e41e27d89ff Mon Sep 17 00:00:00 2001 From: Leo Singer Date: Thu, 1 Nov 2018 18:08:51 -0400 Subject: [PATCH 4/6] Avoid goto --- astropy_healpix/_core.c | 36 +++++++++++++++++------------------- 1 file changed, 17 insertions(+), 19 deletions(-) diff --git a/astropy_healpix/_core.c b/astropy_healpix/_core.c index 52a62b0..45092af 100644 --- a/astropy_healpix/_core.c +++ b/astropy_healpix/_core.c @@ -219,13 +219,14 @@ static void neighbours_loop( static PyObject *healpix_cone_search( PyObject *self, PyObject *args, PyObject *kwargs) { - PyObject *result = NULL; + PyObject *result; static const char *kws[] = {"lon", "lat", "radius", "nside", "order", NULL}; double lon, lat, radius; int nside; char *order; int64_t *indices, n_indices; int64_t *result_data; + npy_intp dims[1]; if (!PyArg_ParseTupleAndKeywords( args, kwargs, "dddis", kws, &lon, &lat, &radius, &nside, &order)) @@ -237,30 +238,27 @@ static PyObject *healpix_cone_search( { PyErr_SetString( PyExc_RuntimeError, "healpix_rangesearch_radec_simple failed"); - goto done; + return NULL; } + dims[0] = n_indices; + result = PyArray_SimpleNew(1, dims, NPY_INT64); + if (result) { - npy_intp dims[] = {n_indices}; - result = PyArray_SimpleNew(1, dims, NPY_INT64); - } - if (!result) - goto done; - - result_data = PyArray_DATA((PyArrayObject *) result); + result_data = PyArray_DATA((PyArrayObject *) result); - if (strcmp(order, "nested") == 0) - { - int i; - for (i = 0; i < n_indices; i ++) - result_data[i] = healpixl_xy_to_nested(indices[i], nside); - } else { - int i; - for (i = 0; i < n_indices; i ++) - result_data[i] = healpixl_xy_to_ring(indices[i], nside); + if (strcmp(order, "nested") == 0) + { + int i; + for (i = 0; i < n_indices; i ++) + result_data[i] = healpixl_xy_to_nested(indices[i], nside); + } else { + int i; + for (i = 0; i < n_indices; i ++) + result_data[i] = healpixl_xy_to_ring(indices[i], nside); + } } -done: free(indices); return result; } From 70e756a6db0aebd8fda2ed7b51770be02d0f43df Mon Sep 17 00:00:00 2001 From: Leo Singer Date: Mon, 5 Nov 2018 10:25:47 -0500 Subject: [PATCH 5/6] Introduce preprocessor macro for invalid index (-1) --- astropy_healpix/_core.c | 24 ++++++++++++++---------- 1 file changed, 14 insertions(+), 10 deletions(-) diff --git a/astropy_healpix/_core.c b/astropy_healpix/_core.c index 45092af..e0a44f6 100644 --- a/astropy_healpix/_core.c +++ b/astropy_healpix/_core.c @@ -18,6 +18,8 @@ #include "ieee754.h" +#define INVALID_INDEX (-1) + /* Data structure for storing function pointers for routines that are specific * to the HEALPix ordering scheme. When we create the ufuncs using * PyUFunc_FromFuncAndData, we will set them up to pass a pointer to this @@ -63,7 +65,7 @@ static void healpix_to_lonlat_loop( double dy = *(double *) &args[3][i * steps[3]]; double *lon = (double *) &args[4][i * steps[4]]; double *lat = (double *) &args[5][i * steps[5]]; - int64_t xy = -1; + int64_t xy = INVALID_INDEX; if (pixel_nside_valid(pixel, nside)) xy = funcs->order_to_xy(pixel, nside); @@ -93,13 +95,13 @@ static void lonlat_to_healpix_loop( int64_t *pixel = (int64_t *) &args[3][i * steps[3]]; double *dx = (double *) &args[4][i * steps[4]]; double *dy = (double *) &args[5][i * steps[5]]; - int64_t xy = -1; + int64_t xy = INVALID_INDEX; xy = radec_to_healpixlf(lon, lat, nside, dx, dy); if (xy >= 0) *pixel = funcs->xy_to_order(xy, nside); else { - *pixel = -1; + *pixel = INVALID_INDEX; *dx = *dy = NPY_NAN; _npy_set_floatstatus_invalid(); } @@ -117,14 +119,14 @@ static void nested_to_ring_loop( int64_t nested = *(int64_t *) &args[0][i * steps[0]]; int nside = *(int *) &args[1][i * steps[1]]; int64_t *ring = (int64_t *) &args[2][i * steps[2]]; - int64_t xy = -1; + int64_t xy = INVALID_INDEX; if (pixel_nside_valid(nested, nside)) xy = healpixl_nested_to_xy(nested, nside); if (xy >= 0) *ring = healpixl_xy_to_ring(xy, nside); else { - *ring = -1; + *ring = INVALID_INDEX; _npy_set_floatstatus_invalid(); } } @@ -141,14 +143,14 @@ static void ring_to_nested_loop( int64_t ring = *(int64_t *) &args[0][i * steps[0]]; int nside = *(int *) &args[1][i * steps[1]]; int64_t *nested = (int64_t *) &args[2][i * steps[2]]; - int64_t xy = -1; + int64_t xy = INVALID_INDEX; if (pixel_nside_valid(ring, nside)) xy = healpixl_ring_to_xy(ring, nside); if (xy >= 0) *nested = healpixl_xy_to_nested(xy, nside); else { - *nested = -1; + *nested = INVALID_INDEX; _npy_set_floatstatus_invalid(); } } @@ -189,9 +191,11 @@ static void neighbours_loop( { int64_t pixel = *(int64_t *) &args[0][i * steps[0]]; int nside = *(int *) &args[1][i * steps[1]]; - int64_t neighbours[] = {-1, -1, -1, -1, -1, -1, -1, -1}; + int64_t neighbours[] = { + INVALID_INDEX, INVALID_INDEX, INVALID_INDEX, INVALID_INDEX, + INVALID_INDEX, INVALID_INDEX, INVALID_INDEX, INVALID_INDEX}; int j; - int64_t xy = -1; + int64_t xy = INVALID_INDEX; if (pixel_nside_valid(pixel, nside)) xy = funcs->order_to_xy(pixel, nside); @@ -207,7 +211,7 @@ static void neighbours_loop( if (xy >= 0) pixel = funcs->xy_to_order(xy, nside); else { - pixel = -1; + pixel = INVALID_INDEX; _npy_set_floatstatus_invalid(); } *(int64_t *) &args[2 + j][i * steps[2 + j]] = pixel; From 3b9fd4b87fb6b8926e33afdc681a902af3fcdf8a Mon Sep 17 00:00:00 2001 From: Leo Singer Date: Mon, 5 Nov 2018 11:02:53 -0500 Subject: [PATCH 6/6] Update changelog --- docs/changes.rst | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/docs/changes.rst b/docs/changes.rst index 8618eb1..506c8f6 100644 --- a/docs/changes.rst +++ b/docs/changes.rst @@ -7,7 +7,8 @@ Changes 0.4 (unreleased) ================ -- No changes yet. +- Rewrite core module in C to make ``healpix_to_lonlat`` and + ``lonlat_to_healpix`` broadcastable over both pixel index and nside. [#110] 0.3.1 (2018-10-24) ==================