Rewrite core in C using Numpy ufuncs

[lpsinger]

As a result, broadcasting across nside and ipix is supported automatically.

Fixes #72, closes #75.

[lpsinger]

There are a few Python 2 issues with this patch. Is it OK to drop Python 2 support from this package?

[cdeil]

@lpsinger - Thanks!

Unless this is urgent for you, I'd suggest we decide about Python 2 support first in #111 and then continue here. (This discussion will probably take a few days, but hopefully within a week we should know.)

[lpsinger]

Unless this is urgent for you, I'd suggest we decide about Python 2 support first in #111 and then continue here. (This discussion will probably take a few days, but hopefully within a week we should know.)

I have a header six.h that I could contribute so that we can decouple the two issues.

[cdeil]

Sure, if you want this in quickly, and you can keep Python 2 support without much extra work, OK to add six.h.

[cdeil]

@lpsinger - Can you have a look at the CI fails?

• https://travis-ci.org/astropy/astropy-healpix/jobs/448633260#L1248
• https://ci.appveyor.com/project/Astropy/astropy-healpix/builds/19931294/job/oro13g9fou0hlm4a#L321

Apart from those, there's a few fails because Astropy isn't compatible with Matplotlib 3.0:

• https://travis-ci.org/astropy/astropy-healpix/jobs/448633256#L1129
• https://travis-ci.org/astropy/astropy-healpix/jobs/448633257#L1243

@astrofrog - What should be done about those? Probably tweak the versions in .travis.yml, e.g. pin MPL to 2.x for now? Can you do the appropriate edit in master?

[cdeil]

Oops, I didn't refresh the page, so was looking at an older status.

The only CI fail that remains here is the unrelated one for MPL 3.0.

@astrofrog - Could you please do a code review on this PR?

[cdeil]

@lpsinger - Two small inline comments from me.

Overall this looks good and I'm +1 to merge.

But I'm not really familiar with Python C interfacing much, so @astrofrog - if you could review this, or ping someone from the Astropy team that knows C / Numpy / ufunc well, that would be good. Although if no-one has time, I think we could also just merge this in for now, I don't really have a concern.

[astrofrog]

This looks good to me - the only thing perhaps missing (unless I didn't see it) would be a test to demonstrate the broadcasting functionality of the ufuncs. Otherwise, our tests are pretty extensive, so if they all pass then this is good to go.

[cdeil]

I see that you have this everywhere in Python when calling into C:

    if _validate_order(order) == 'ring':
        func = _core.lonlat_to_healpix_ring
    else:  # _validate_order(order) == 'nested'
        func = _core.lonlat_to_healpix_nested

and then there is the order_funcs in the C code with the IMO quite complex way to select the func.

I really don't know, just as a question: is there a simpler way to do this, avoiding the whole *data and order_funcs stuff? E.g. by having a wrapper function that selects or helper function that applies the conversion? Related: I also see strcmp in one place - would int there to select ordering be better that str?

One small suggestion: you use -1 in many places as "NA" placeholder in the C code, right? Maybe make that a global constant instead of the hard-coded -1 for better readability?

@lpsinger - Do as you think best, and go ahead and merge this in.

Should we make a release to see if this works OK in package distribution?

[cdeil]

Why do we have mostly int64_t for integers, but then int nside? I guess it doesn't matter because nside will fit in an int, but wouldn't it be simpler to have int64_t everywhere so that we never have to think about what to put as type? - just a thought. Do as you think best and 👍 to merge this in.

[lpsinger]

Such a change would affect the astrometry.net code as well, so that's a bigger question.

[lpsinger]

I see that you have this everywhere in Python when calling into C:

    if _validate_order(order) == 'ring':
        func = _core.lonlat_to_healpix_ring
    else:  # _validate_order(order) == 'nested'
        func = _core.lonlat_to_healpix_nested

and then there is the order_funcs in the C code with the IMO quite complex way to select the func.

I really don't know, just as a question: is there a simpler way to do this, avoiding the whole *data and order_funcs stuff? E.g. by having a wrapper function that selects or helper function that applies the conversion?

In Healpy, the ordering is treated as a flag (nested=False or nested=True) and it is broadcastable like any other parameter. Switching between the two indexing modes is done inside the inner loop in C. That would certainly be simpler.

I judged that the user would never want to pass the ordering as a vector, that it would never need to be broadcastable, and I tried to avoid the potential performance hit of a conditional in the inner loop. However, the function pointer stuff might be equally costly, and is a bit less readable.

We could use preprocessor macros to make two versions of each loop function, one for each ordering. That would be more performant, but also more code and less readable code.

Another option is that the _core module could wrap healpixl_nested_to_xy, healpixl_xy_to_ring, healpixl_to_radec, and radec_to_healpixlf as ufuncs; then healpix_to_lonlat and lonlat_to_healpix could be composed from those functions in Python. That might be the simplest solution of all, but the memory access pattern would be sub-optimal.

Related: I also see strcmp in one place - would int there to select ordering be better that str?

We could use PyUnicode_RichCompare to do the comparison with Python and take advantage of string interning, but that would take several more lines of code and would come with its own reference counting overhead. But it's only compared once, so none of that would make a measurable difference in performance.

One small suggestion: you use -1 in many places as "NA" placeholder in the C code, right? Maybe make that a global constant instead of the hard-coded -1 for better readability?

Done.

[lpsinger]

Should we make a release to see if this works OK in package distribution?

I tried some of my own code that uses astropy-healpix through reproject, and the results looked good.

[cdeil]

In Healpy, the ordering is treated as a flag (nested=False or nested=True) and it is broadcastable like any other parameter. Switching between the two indexing modes is done inside the inner loop in C. That would certainly be simpler.

I judged that the user would never want to pass the ordering as a vector, that it would never need to be broadcastable, and I tried to avoid the potential performance hit of a conditional in the inner loop. However, the function pointer stuff might be equally costly, and is a bit less readable.

That does sound a bit simpler and easier to read to me.

Agreed that the broadcasting feature isn't needed here, and performance is probably very similar, so the only improvement here would be code simplicity, not needing to know about the not well-documented data option of ufuncs and the function pointer complexity.

@lpsinger @astrofrog - What do you think? Just keep as-is or change to this?