WIP: Reimplement core_cython routines as ufuncs

[lpsinger]

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.

[coveralls]

Coverage increased (+0.6%) to 96.767% when pulling 1a12149 on lpsinger:core-cython-ufuncs into f99445a on astropy:master.

[cdeil]

@lpsinger - Thanks!

I've never seen this, so can't really comment on the proposal.

@astrofrog @dstndstn - Thoughts?

@lpsinger - So you think this is better because it allows to remove type casting and array shape handling code in our Python wrapper functions, right? Could you make the corresponding changes for one example so that we can see the difference here on a concrete example? This would also mean that this code is then actually executed by the tests, and we can see if it works, especially on Windows.

[astrofrog]

I agree with @cdeil that I'd like to see what code can be removed as a result to properly understand the benefits.

[lpsinger]

I agree with @cdeil that I'd like to see what code can be removed as a result to properly understand the benefits.

OK, done. Note that I moved the handling of order from core_cython to core because it's more natural to do the string comparison outside of the ufunc. I also combined the with_offset and no offset versions into a single ufunc, because the ufunc can broadcast the default (scalar) values as necessary.

[cdeil]

For me, as a typical Astropy coder with a little, but mostly superficial familiarity with Cython, and no familiarity with ufuncs, the ufunc version seems even a little more cryptic that what we have now.

I mean things like this:

cdef void healpix_to_lonlat_nested_loop(
        char** args, intp_t* dimensions, intp_t* steps, void* data) nogil:
dx = (<double*> &args[2][i * steps[2]])[0]

Probably it's mostly that I have never seen this ufunc stuff before, and once one gets used to it, then for experts like @lpsinger it is the nicer version?

Uniformity of Cython code within Astropy core and affiliated packages is one concern for maintainability, i.e. if this alternative way is better, than possibly larger fractions of code should be migrated long-term. (but the process could still start here with astropy-healpix)

Given that the goal of this package is to be merged as astropy.healpix into astropy core when in good shape, I think it should be up to @astrofrog or other Astropy maintainers to make this decision.

@astrofrog - Thoughts?

[lpsinger]

Probably it's mostly that I have never seen this ufunc stuff before, and once one gets used to it, then for experts like @lpsinger it is the nicer version?

I've written lots of ufuncs in C, but I'm definitely not a Cython expert. The strided pointer dereferencing (such as the <double*> &args[4][i * steps[5]] part) could probably be done more idiomatically in Cython using typed memoryviews, though I couldn't figure out how to make that work.

I agree that the ufunc approach is going down a bit closer to the C and Numpy C API level, but the payoff is greater robustness and simplicity for the user by letting Numpy deal with the broadcasting, casting, and type checking for you.

In Cython, it may be possible to reduce the boilerplate to create a ufunc by writing a little function decorator.

[lpsinger]

@astrofrog, what do you think?

[astrofrog]

@lpsinger - in general I am in favor of this approach, but have a few high-level comments:

• To address @cdeil's concerns, I think that the first time you use notation like <double*> &args[4][i * steps[5]] I think it would be good to write out a comment explaining what this does just so that contributors not familiar with C are not completely lost.

• Since they basically look like C code at this point, is there any reason to write the loop functions in Cython or would it make sense to write them in pure-C? If the OpenMP stuff harder to do? I'm not saying we should necessarily do this, but just curious about the pros and cons.

• I originally tried to avoid writing separate functions for order='ring' and order='nested' to avoid code duplication, which becomes important for routines like the interpolation. So I wonder if the loop functions should behind the scenes call a common function with an order argument to avoid this?

Thanks for your work on this!

[lpsinger]

Since they basically look like C code at this point, is there any reason to write the loop functions in Cython or would it make sense to write them in pure-C? If the OpenMP stuff harder to do? I'm not saying we should necessarily do this, but just curious about the pros and cons.

If that's an option, I'd prefer pure C. Personally, I find it easier to read C than to try to figure out how Cython syntax maps to C. And OpenMP is just as easy in C.

I originally tried to avoid writing separate functions for order='ring' and order='nested' to avoid code duplication, which becomes important for routines like the interpolation. So I wonder if the loop functions should behind the scenes call a common function with an order argument to avoid this?

They could, but I think that it might be confusing if it was possible to broadcast over the order argument. And although it probably doesn't matter from a performance standpoint, I'd rather not have an if statement to select between orders inside the loop function.

[astrofrog]

If that's an option, I'd prefer pure C. Personally, I find it easier to read C than to try to figure out how Cython syntax maps to C. And OpenMP is just as easy in C.

That is an option as far as I'm concerned - the only question then is do we keep Cython just for the high-level declaration of the ufuncs? (since that might avoid some boilerplate code when writing solely a C extension?)

They could, but I think that it might be confusing if it was possible to broadcast over the order argument. And although it probably doesn't matter from a performance standpoint, I'd rather not have an if statement to select between orders inside the loop function.

Ok, in that case we should probably have one source file for each mode to make it easy to diff them and make sure they are in sync (and in future we could even make a jinja template to avoid duplication of code in the repo, but let's not worry for now).

[cdeil]

@lpsinger @astrofrog - Can we try to move this discussion along and make a decision within the next week or two? I want to spend some time before and after the meeting July 2/3 on astropy-healpix, but at the moment I'm not quite sure in what direction to move.

@mhvk as Numpy ufunc expert and proponent (see e.g. this discussion) - maybe you could chime in here as well and make a suggestion how you think this new astropy.healpix should be written?

As for me: I'm reading the HEALPix paper and the code in this repo at the moment and trying to understand the details. It's not clear to me yet why we have 4000 lines of code in cextern/astrometry.net (plus another 200 in interpolation.c and 500 lines of wrapper code in astropy_healpix/core_cython.pyx.

Compare this to the license-compatible https://github.com/michitaro/healpix/blob/master/src/index.ts or https://github.com/ziotom78/Healpix.jl/blob/master/src/Healpix.jl that have ~ 500 lines of code and have very similar functionality (except interpolation missing) and could also serve as basis for a rewrite. I'm not saying we should do that, I'm just saying that we could consider a rewrite in whatever form we consider the best long-term solution for astropy.healpix. Given that we have tests and an API, a port of 500 or 1000 lines from another similar language should only be a day or two of work.

Some concrete questions:

• I think a pure Python & Numpy solution clearly not possible, because some functions require loops and can't be done as Numpy extressions?
• Numba is not an option yet for Astropy core yet, because of the heavy runtime dependency, right?
• If Numpy only or Numba aren't an option, that would always leave us with a solution involving Cython, right?
• Would a rewrite in Cython only bring benefits, or would you prefer a C library with Cython wrapper like we have now?
• Ufuncs or no ufuncs?

All (also @dstndstn): please comment and make a recommendation here, so that we can move astropy.healpix forward.

[dstndstn]

To be fair, most of the code imported from astrometry.net is "bl" (block-list) stuff -- just a generic C linked-list/array hybrid container type, which is only used in the healpix-rangesearch code to return its results. When wrapped in python, one would be better to perhaps use callbacks and native python list appends.

On the licensing front, you probably only need stuff that is 3-clause BSD licensed. If that is a problem, I am willing to re-license.

[cdeil]

@dstndstn - Apologies, I didn't mean to criticise your cextern/astrometry.net, I'm still in the phase of learning what it does and trying to figure out what we want for astropy.healpix.

One thing that would help would be a README there with a few notes, e.g. that one should start reading in healpix.h and which of the other ~ 20 files are relevant, i.e. a brief overview of the C library.

The part of the code that I find hardest to understand is the query stuff.

• I see http://healpy.readthedocs.io/en/latest/healpy_query.html
• Is it correct that this C code here only does disk query? Or are other Healpix queries also supported? E.g. I see a point_in_polygon function, but is this enough that it would allow us to implement the equivalent of healpy.query_polygon?
• The only description / paper I could find that describes HEALPix query functions is here. I'll have a look at that now, but if there is a description of the algorithm we use in this package, a link would be very much appreciated.

When wrapped in python, one would be better to perhaps use callbacks and native python list appends.

That sounds very slow and memory inefficient. Should we really consider or try to implement that?

[astrofrog]

Having thought about it more and seen what @mhvk did with ERFA too, I'm in favor of ufuncs. I'm happy to leave it up to @lpsinger to decide whether to do a mix of Cython + C or just C.

In terms of re-writing the C code to shorten it: no objections, but at the same time not really urgent.

I think a pure Python & Numpy solution clearly not possible, because some functions require loops and can't be done as Numpy expressions?

I think that's correct.

Numba is not an option yet for Astropy core yet, because of the heavy runtime dependency, right?

Yes

If Numpy only or Numba aren't an option, that would always leave us with a solution involving Cython, right?

Or pure C, without Cython

Would a rewrite in Cython only bring benefits, or would you prefer a C library with Cython wrapper like we have now?

If we were doing a re-write, then actually writing it as a standalone C library that we wrap could be good in the sense that there would then also be a BSD C implementation. If we care about performance then C will always be faster than Cython.

Ufuncs or no ufuncs?

I think ufuncs

[cdeil]

I think ufuncs; Or pure C, without Cython

Can you point to an example, please?
So the ufunc is written in C and it's a hand-written Python C extension file?

[astrofrog]

So the ufunc is written in C and it's a hand-written Python C extension file?

Yes - Astropy does this in astropy._erfa but that is a bit more complicated because the C file is generated on the fly and is very long!

[mhvk]

@cdeil - there is a nice tutorial on writing your own ufuncs [1], though for generalized ufuncs (which act on arrays), I learned most from looking at the numpy test cases [2].

[1] https://docs.scipy.org/doc/numpy/user/c-info.ufunc-tutorial.html
[2] https://github.com/numpy/numpy/blob/master/numpy/core/src/umath/_umath_tests.c.src#L74
(typing will not be relevant for you)

[cdeil]

I think this discussion has concluded, and surprisingly all agree on how to implement this (right?): using ufuncs in C.

@lpsinger - Next Monday / Tuesday we'll have the coding sprint on astropy-healpix, and I think we should try to resolve all open PRs before. OK if we close this one, and then a new one is started using C? @lpsinger - Do you have time to write that, or should I give it a try?

@mreineck - To follow up on #90 (comment) : indeed, if you'd be willing to put BSD on the Python interface in https://gitlab.mpcdf.mpg.de/ift/pyHealpix , that would probably be the perfect example to follow for us here.

[mreineck]

@mreineck - To follow up on #90 (comment) : indeed, if you'd be willing to put BSD on the Python interface in https://gitlab.mpcdf.mpg.de/ift/pyHealpix , that would probably be the perfect example to follow for us here.

If you are only interested in the interface (i.e. the specification which functionality is provided, and in which form), you don't need to worry about licensing at all - interfaces are not copyrightable (even if Oracle apparently disagrees).

But the pyHealpix interface is designed around the concept that the Nside parameter and the ordering scheme are fixed within each call, and as far as I understand, this is not how you want to do it.

In any case, feel free to take ideas concerning the interface! If you still want a copy of the files with a BSD header instead of GPL, let me know!

[cdeil]

I wanted a good example how to do ufuncs in C; but now I remember that you're using pybind11, so never mind. I'll look at what @mhvk suggested above in #75 (comment) .

[mhvk]

@cdeil - feel free to ping when you really start implementing the ufuncs. Being not much of a C programmer myself, the main thing that tripped me up was that one has to ensure that items passed in to PyUFunc_FromFuncAndDataAndSignature are kept in memory, i.e., are static variables (or use memory explicitly allocated). Obviously, the examples do this right, but without comment.

The other part that may need some care is that for arrays your function gets strides passed in, i.e., the axes are not necessarily contiguous and you need to make a copy if the functions that do the actual work expect contiguous input (or do the relevant copy in a python wrapper).