Skip to content

Commit

Permalink
WIP: Reimplement core_cython routines as ufuncs
Browse files Browse the repository at this point in the history
This eliminates the need for much of the explicit broadasting and
type checking.

At the moment, the patch just implements `healpix_to_lonlat_ring`.
I'd like to ask for feedback on the basic approach before
undertaking conversion of the other functions.

Fixes #72.
  • Loading branch information
lpsinger committed Apr 24, 2018
1 parent fcb4089 commit bb48842
Showing 1 changed file with 48 additions and 0 deletions.
48 changes: 48 additions & 0 deletions astropy_healpix/core_cython.pyx
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,9 @@ npy_double = np.double
npy_intp = np.intp
npy_int64 = np.int64

np.import_array()
np.import_ufunc()


def _validate_order(str order):
# We also support upper-case, to support directly the values
Expand All @@ -32,6 +35,51 @@ def _validate_order(str order):
raise ValueError("order must be 'nested' or 'ring'")


@cython.boundscheck(False)
cdef void healpix_to_lonlat_ring_loop(char** args, intp_t* dimensions,
intp_t* steps, void* data) nogil:
cdef intp_t n = dimensions[0]
cdef intp_t i
cdef int64_t xy_index, healpix_index, nside
cdef double dx, dy
cdef double* lon
cdef double* lat

dx = 0.5
dy = 0.5

for i in prange(n, nogil=True, schedule='static'):
healpix_index = (<int64_t*> &args[0][i * steps[0]])[0]
nside = (<int64_t*> &args[1][i * steps[1]])[0]
lon = <double*> &args[2][i * steps[2]]
lat = <double*> &args[3][i * steps[3]]
xy_index = healpixl_ring_to_xy(healpix_index, nside)
healpixl_to_radec(xy_index, nside, dx, dy, lon, lat)


cdef void *no_data[1]
cdef np.PyUFuncGenericFunction loop_funcs[1]
cdef char types[4]
no_data[:] = [NULL]
loop_funcs[:] = [healpix_to_lonlat_ring_loop]
types[:] = [np.NPY_INT64, np.NPY_INT64, np.NPY_DOUBLE, np.NPY_DOUBLE]


healpix_to_lonlat_ring = np.PyUFunc_FromFuncAndData(
loop_funcs,
no_data,
types,
1, # number of supported input types
2, # number of input args
2, # number of output args
0, # `identity` element, never mind this
"healpix_to_lonlat_ring", # function name
"TODO", # docstring
0 # unused
)



@cython.boundscheck(False)
def healpix_to_lonlat(np.ndarray[int64_t, ndim=1, mode="c"] healpix_index,
int nside, str order):
Expand Down

0 comments on commit bb48842

Please sign in to comment.