Alternative implementations

[cdeil]

To learn about HEALPix and explore the possibility of a rewrite of what we use here in astropy-healpix, I've translated https://github.com/michitaro/healpix by @michitaro to Python.

IMO one option could be to translate this to C and evolve it, the other is to evolve cextern/astrometry.net.

@dstndstn - How do the two implementations compare? Are they very similar, or does yours differ in the method used, or have numerical tricks already that makes it more robust?

I think a pure-Numpy implementation is possible, but can't give great performance, so not even worth trying?

Using the one test case from #84 (comment) I see the same precision error near the south pole as our implementation has (see #34 (comment)). @dstndstn or @michitaro - Does one of you have time to attempt to improve this, or should I try to understand this and attempt a fix?

[astrofrog]

@cdeil - I'd rather not include playground in the master branch to keep the history clean. How about creating a new branch playground for that (forked from master)

[cdeil]

@astrofrog - We don't have to merge this, this is just for discussion.
Others can check out my branch and hack on it here if they like.
Or if you think it's better, I could make a separate repo under my account.

That said: one option could be that we just make a single-commit clean PR to astropy.healpix for the core repo. I don't think it's useful to drag this version history in, it includes the package template and other version history that's IMO not appropriate for the Astropy core repo. If you agree, then we could also merge this and use playground in this repo to play a bit.

[astrofrog]

@cdeil - ok no problem - even if we end up merging in one go into astropy core I think we should still strive to keep the repo here clean for now (since it's not clear on what timescale the merging would happen). So I'd still advocate either just keeping this PR open or, if you need to merge it, merging it into a playground branch.

[dstndstn]

I don't really have any numerical tricks; it's a fairly naive implementation.

I think a pure numpy implementation could work, as long as one was careful with vectorization. What is the main operation you want to optimize? radec_to_healpix and healpix_to_radec?

I can spend little time working on the accuracy near the poles.

[cdeil]

I think a pure numpy implementation could work, as long as one was careful with vectorization.

But it would always use a factor of 2-3 more memory and CPU than a C implementation that avoids temp arrays, no?

What is the main operation you want to optimize? radec_to_healpix and healpix_to_radec?

Yes, pix2ang and ang2pix are the main performance focus:
http://astropy-healpix.readthedocs.io/en/latest/performance.html
Although I have to say: for me, performance is not critical; I'm not even sure we should use OpenMP at all. To be discussed next Monday at the HEALPix meeting.

I can spend little time working on the accuracy near the poles.

That would be wonderful, I would probably spend a ton of time to understand and fix that issue.

[lpsinger]

I prefer the official CHEALPix library: it's written in very plain C, is short and sweet, and I maintain the package for Debian and MacPorts.

[cdeil]

@lpsinger - A major motivation for the work here is the GPL vs BSD license issue.
See http://astropy-healpix.readthedocs.io/en/latest/about.html and https://docs.google.com/presentation/d/1ewbKWRneSWkMVyio-WL3GA9gir1oUuF2OfaWObT20G4/edit?usp=sharing

The original CHEALPix is not possible given that we want to use BSD in Astropy and Astropy-affiliated packages.

[lpsinger]

Have you asked the CHEALPix maintainers if they would be willing to dual-license it?

[cdeil]

Yes, there was a request: https://sourceforge.net/p/healpix/mailman/message/34929929/
If you know them personally, you could ask again, but otherwise I would suggest we proceed with a fresh implementation here. Preferring GPL is a valid point, there's good reasons to pick it over a liberal license.

[mreineck]

Just a quick technical comment (sorry, I can't really help you with the licensing terms of chealpix, even though I rewrote most of it...):

If you are concerned about performance and applying Healpix functions in ufunc style, you should consider using a pybind11 interface (https://github.com/pybind/pybind11). Compared to cython, this is much, much cleaner, and, as far as I can tell, way more efficient.

If you want to play around with a chealpix-like Python Healpix library, you can have a look at https://gitlab.mpcdf.mpg.de/ift/pyHealpix.

Yes, it wraps the (GPL'd) Healpix C++ library ... but if you want to use the glue code (i.e. all sources sitting directly in the root directory of the package) for anything, I can give it to you under BSD terms.

[cdeil]

@mreineck - Thank you!

@astrofrog - Can you please comment whether C++ / pybind11 would be accepted for astropy.healpix in Astropy core?
The question and suggestion came up recently also here:
https://groups.google.com/d/msg/astropy-dev/0xYsvGl2AGs/iwW_LUuiAAAJ

I think maybe we should avoid it for now, always cast to a single type when calling into C?

And then in a few years, if C++ / pybind11 is used in Astropy, it could be rewritten to be a little more efficient for float32 or small int types, while keeping backward-compatibility.

The performance hit IMO isn't large enough to warrant introducing the C++ / pybind11 dependency in Astropy now or to complicate our code here via some C or Cython templating.

[cdeil]

@mreineck - Is there a way to find https://gitlab.mpcdf.mpg.de/ift/pyHealpix except via Google?

Would significantly extending https://en.wikipedia.org/wiki/HEALPix and trying to have that as the overview / up to date summary on HEALPix could be an option?

Probably you're maintaining http://healpix.sourceforge.net/ but there e.g. https://gitlab.mpcdf.mpg.de/ift/pyHealpix or astropy-healpix and others aren't mentioned, and it already has several dead links and probably only very few people have edit previliges (and sourceforge shows annoying ads on many pages).

[mreineck]

@mreineck - Is there a way to find https://gitlab.mpcdf.mpg.de/ift/pyHealpix except via Google?

Not really - and that's deliberate to some extent. This is a local toy project, and should in no way give the impression that we want to compete with healpy. Also, time for support and to implement feature requests is quite limited. So I don't want to draw too much attention to this.

I agree that SourceForge is becoming increasingly obnoxious - I'd be happy to move, but now that Github has been taken over by Microsoft I'm not feeling at home there either. My personal favorite would be to move all Healpix ports maintained by myself to the MPCDF, but the other Healpix maintainers prefer having everything under one roof.

[mreineck]

BTW, I did a quick performance comparison running pix2ang for all pixels in an Nside=1024 map. Using the pure Python implementation above, this took roughly 2 minutes on my laptop. Using pyHealpix with one call per pixel took 12.5 seconds; doing all work in one vectorized call took 0.52 seconds.

[cdeil]

I didn't necessarily mean that you should move code repos. For users, I think http://healpix.sourceforge.net/ and https://en.wikipedia.org/wiki/HEALPix are the entry points to info about HEALPix. Sourceforge is annoying because of the ads, so you could buy a domain like healpix.org and have a static page at a stable location and host anywhere. But that's a lot of work, thus my suggestion to use the free wiki at https://en.wikipedia.org/wiki/HEALPix and extend the information there to have an up-to-date summary.

Using the pure Python implementation above, this took roughly 2 minutes on my laptop

Sure, Python is slow. That's just for me to learn how HEALpix works, we don't plan to use that. Probably well-written C can't be beaten, right? Although it could be interesting to try to add @jit from http://numba.pydata.org/ and see if that gives reasonable performance. (just for fun, we don't want the dependency on Numba in astropy.healpix)

Note that we have a little benchmark script here in astropy-healpix that shows that we're a few times slower than healpy: http://astropy-healpix.readthedocs.io/en/latest/performance.html
If https://gitlab.mpcdf.mpg.de/ift/pyHealpix is even faster than healpy or interesting to add to the comparison, we could add that.

[mreineck]

Ah, sorry, I had misunderstood your above comment

The performance hit IMO isn't large enough to warrant introducing the C++ / pybind11 dependency in Astropy now or to complicate our code here via some C or Cython templating.

I somehow assumed you were talking about using compiled code vs. using pure Python, which of course was totally wrong.

(In my opinion, you don't need to worry about float32 vs float64 performance at all, if you are calling from Python. This feels very much like a micro-optimization.)

Concerning the performance of pyHealpix vs. healpy, they should be very similar for large arrays, and pyHealpix should be somewhat faster for small ones; nothing really exciting there.

[lpsinger]

Yes, there was a request: https://sourceforge.net/p/healpix/mailman/message/34929929/
If you know them personally, you could ask again, but otherwise I would suggest we proceed with a fresh implementation here. Preferring GPL is a valid point, there's good reasons to pick it over a liberal license.

OK, unfortunately, re-licensing chealpix is not an option. Why, again, is a GPL license unacceptable?

[cdeil]

Why, again, is a GPL license unacceptable?

Because Astropy (and Python, Numpy, ...) has a liberal license, and (I presume) the Astropy devs don't want to put GPL restrictions on their use now, so can't include GPL code or the whole would effectively be under GPL terms.

[mreineck]

I really don't understand why this is repeated over and over again ...
No, the whole would not be under GPL terms; BSD code can freely call GPL code and does not magically become tainted by it.

One exception: if you distribute a binary package that contains both GPL and BSD code, you have to make the BSD code available under GPL as well ... but definitely not exclusively.

As you probably remember, this is the answer I got from the FSF ... and I'm pretty sure they should know.

I'm happy to post the email exchange here if requested.

[cdeil]

I don't think it makes sense to continue license discussion. Apologies if the statement I made wasn't fully correct! The bottom line is that Astropy (and Numpy, Scipy, scikit-image, scikit-learn, ...) chose BSD and doesn't want to include or link against GPL code, so including the CHEALPix GPL code isn't an option.

[coveralls]

Coverage remained the same at 96.214% when pulling 8991e23 on cdeil:alternatives into b7a787c on astropy:master.

[fxpineau]

@mreineck The best way to find things is probably to compile the Javadoc and to look at it before going through the code.
For ang2pix(), you should look at the class 'HealpixNestedHashComputer', implementing the interface 'HashComputer'. I call the method 'hash()' since we can see ang2pix() as computing the hash() value of a position on the sphere.
For pix2ang(), you should look at the class 'HealpixNestedVerticesAndPathComputer' implementing the interface 'VerticesAndPathComputer'. I called the method 'center()' since it returns the coordinates on the unit-sphere of the center (in the Euclidean projection plane) of a cell.

[mreineck]

Thanks a lot (especially the terminology explanation is very helpful)!

[mreineck]

I know that the timing looks strange, but a few days ago I have in fact re-started the Healpix-internal discussion to release a reference implementation of Healpix core functionality under 3-clause BSD. The responses seem slightly more favorable this time, so I have some hope...
The code in question can be inspected at https://gitlab.mpcdf.mpg.de/mtr/healpix_reference; note that I'm still using GPL for the moment.

Of course this is by no means a complete package with cone/polygon queries etc., but it supports accurate conversion between a set of continuous coordinate systems:

• 3-vectors
• colatitude/longitude
• continuous Healpix coordinates (face number plus x, y in [0;1]) (only internal at the moment)

and of course the discrete Healpix RING and NESTED schemes.
Performance is not the best that can be done (although it's still pretty good), but since this is supposed to be a reference implementation, I considered compactness and clarity more important.

[cdeil]

@mreineck @fxpineau - From the Astropy perspective, this is wonderful, that these license-compatible HEALPix libs are emerging, thank you!

Is one of you willing to schedule a call before Christmas, try to find a date and time that works for all involved / interested? Maybe only offer dates after Nov 26 to give @mreineck and the others some time to do their license discussion?

---

The best way to find things is probably to compile the Javadoc and to look at it before going through the code.

If you could put that up online, that would be very nice.
It's easy to do with Github pages:
https://vaadin.com/blog/host-your-javadoc-s-online-in-github
Although it's tedious to have to update it manually.

I don't know about Java, there might also be a similar service as ReadTheDocs for Python (which we use for https://astropy-regions.readthedocs.io ) that has integration with Github, and builds and host updated docs versions automatically after each commit to master.

Generally for CI, if you're willing to set it up, I'd recommend https://github.com/marketplace/azure-pipelines. It's very new and I just started using it a bit, but it's very nice, many open-source projects are migrating there, and it'll give you a full and free and integrated CI solution for anything on https://github.com/cds-astro - of course there's the investment to need one person to learn it and set it up once, after that the maintenance effort should be very small to keep it running for the next years.

[cdeil]

One thing many people need is interpolation. @astrofrog implemented this here specifically for this package: https://github.com/astropy/astropy-healpix/blob/master/astropy_healpix/interpolation.c
and it's used e.g. in https://reproject.readthedocs.io/en/stable/healpix.html .

As far as I can see that's a feature that's still missing in https://github.com/cds-astro/cds-healpix-java , and would be something that we'd want to part of whatever C lib we use for astropy-healpix. At least I think it would be nicer to have / maintain that interpolation code there, and have a pretty clear separation into three parts: a C lib that we can then bundle in https://github.com/astropy/astropy/tree/master/cextern like Astropy does for WCSLib and ERFA already, the C wrapper code (recently rewritten in C as ufuncs by @lpsinger ), and then the Python code which would be in a new astropy/healpix subfolder (for now developed here, then migrated to Astropy when ready, hopefully landing next year).

[mreineck]

As far as I can see that's a feature that's still missing in https://github.com/cds-astro/cds-healpix-java ,

Most likely because this package doesn't have support for the RING numbering scheme (and RING indices are needed to perform physically meaningful interpolation).

I probably won't be able to provide that either; the internal agreement is that the reference implementation will only contain the essentials, especially those bits where special care is needed to achieve good accuracy. This means coordinate transforms, but not neighbor search, cone/polygon queries, SHTs or interpolation.

[cdeil]

But ring to nested index conversion is available or possible to add with limited effort, no? Then with an extra index conversion or a data copy, interpolation can be implemented without much effort, it just wouldn't be efficient memory and CPU wise?

[fxpineau]

I'am not familiar with the "interpolation" problem. Do you have docs about it?
Also, which function would you need in the RING scheme? (I have not implemented RING so far since it is not part of the internal CDS needs).
For information, I am working on a RUST version of the library. As it is a compiled language, RUST libraries can be called from Python just like C one's (and generate WebAssembly files). I will probably put it on github in a near future.

[mreineck]

I'am not familiar with the "interpolation" problem. Do you have docs about it?

It's basically bilinear interpolation in RING scheme.
The function takes a pointing on the sphere, and then

• determines the ring above and below the pointing (special treatment for the poles is necessary)
• on those rings, finds the closest pixels to the left and right of the pointing
• computes bilinear interpolation weights for these pixels with respect to the pointing
• returns pixel indices and weights.
  This is somewhat more consistent than to do bilinear interpolation in NESTED scheme, since this leads to problems at the seams between the base pixels.

[fxpineau]

Thank you Martin.

If I understand correctly, you have a value for each healpix cell and you want to perform a bilinear interpolation for a given position (pointing) on the sky.
I guess that the pointing is, for example, the center of a 2D image pixel.

First, I think that the problem is exaclty the same in RING and in NESTED schemes, the only difference comes from the way it is implemented, right?

If you consider a position very close from (ra, dec) = (90, 41) in degrees, I think that the interpolation algorithm should consider only 3 cells, not 4.

From a quick thought, I have the following algorithm in mind (maybe it is slow in practice, but I don't think so):

• find the cell (C) the position belong to
• look if the the position is in the S, E, N or W quarter of the cell (from shifted, rotated, scaled positions used to compute the cell value it should be fast)
  • if S: look at C, SW, S, SE cells
  • if E: look at C, SE, E, N cells
  • if N: look at C, NE, N, NW cells
  • if W: look at C, NW, W, SW cells
• compute the weights from the 4 cells, except for the 24 cells having only 3 neighbours (no S or N or E or W), for witch the weights will be based on 3 cells.

I tend to think that the weights are computed from the Euclidean (dx, dy) in the projection plane (in witch (dx, dy) are the coordinates differences between the pointing and the cells centres), right?

[mreineck]

First, I think that the problem is exaclty the same in RING and in NESTED schemes, the only difference comes from the way it is implemented, right?

If I understood you correctly, no.

The selected "neighbors" that take part in the interpolation are different for the scheme you propose and the one currently implemented in the Healpix library. Our approach always uses four neighbors, which are located on two different rings (let's ignore the poles for the moment); yours selects typically four, but sometimes three, and they lie on up to three different rings.

If we assume for a moment a map with perfect azimuthal symmetry (i.e. all pixels on a ring have the same value), then the results of our interpolation will also have this symmetry, but yours most likely won't. (That's just to illustrate the differences, not to make a point which algorithm is better.)

It's hard to say which approach is "better"; and exactly because we can't answer this question objectively, we never advertised the interpolation functionality very much.

Edit: I forgot to mention that we compute interpolation weights directly from theta and phi, no projection is made.

[cdeil]

Just to say: the use case for interpolation is mainly to transform HEALPix <-> WCS images as shown here: https://reproject.readthedocs.io/en/stable/healpix.html

This is very common for visualisation, and it can be useful for analysis, e.g. in gamma-ray astronomy we have all-sky diffuse models in HEALPix, but then analyses are done for certain fields of view of ~ 10 deg across in a local TAN projection at the pointing position, and for that one reprojects the diffuse model onto that WCS image.

If we assume for a moment a map with perfect azimuthal symmetry (i.e. all pixels on a ring have the same value), then the results of our interpolation will also have this symmetry, but yours most likely won't. (That's just to illustrate the differences, not to make a point which algorithm is better.)

@mreineck - is there a real advantage / use case where this matters?

At least for the use cases I'm aware of, it doesn't matter.

Although, from a testing / compatibility / documentation standpoint, if possible, the benefit of agreeing on how to do it would be very advantageous. But if it's not possible or not easily, then something is much better than nothing.

This has worked well to migrate some scripts or codes to astropy-healpix and obtain the same results (https://astropy-healpix.readthedocs.io/en/latest/healpy_compat.html), and the main issues are where results differ (the precision issues we had at the poles, different behaviour for interpolation and disk query).

Here's what we currently have:
https://github.com/astropy/astropy-healpix/blob/master/astropy_healpix/interpolation.c
I'm not sure what exactly it does, @astrofrog wrote it.

[mreineck]

@mreineck - is there a real advantage / use case where this matters?

I don't know ... so far I've only seen interpolation used in visualization, where it shouldn't matter.

For serious reprojection applications I'm not sure whether interpolation should be used at all; shouldn't such algorithms conserve flux? Interpolation won't do the trick then...

The rationale for choosing the ring-based approach was simply that theta/phi coordinates are more "physical" than Healpix-specific coordinates and that we avoid problems which would otherwise occur at the edges/corners of the base pixels.

[cdeil]

For serious reprojection applications I'm not sure whether interpolation should be used at all; shouldn't such algorithms conserve flux? Interpolation won't do the trick then...

Yes, agreed. I guess very few people do analysis and measurements on images on a mix of HEALPix and WCS. If an algorithm exists to reproject HEALPix that's flux-conserving, it would be great to have. That could be discussed at https://github.com/astrofrog/reproject/issues/70

But I think the interpolation is really important to have, because pretty much anyone that has HEALPix data wants to visualise it and compare with other data and for that very often wants to go to WCS.
See http://docs.astropy.org/en/stable/visualization/wcsaxes/index.html in Astropy.

[fxpineau]

I consider the HEALPix projection plane, rotated of 45 degrees such that each cell is a square (pixel).
If you add a border of width = 1 pixel in which you copy the value of the neighbours on the sphere, you
see that you have a problem at the north and south of the 4 equatorial base cells.
(See right part of Fig. 1 in my ADASS poster) )

yours selects typically four, but sometimes three, and they lie on up to three different rings.

The four of three cells are always on 2 different rings (2 strait lines of slope = -1 in the rotated frame), except if the pointing is abose the first ring, or under the last one (poles).

but yours most likely won't.

It depends on how you choose to compute the weights in the 3 cells case.

For serious reprojection applications I'm not sure whether interpolation should be used at all

I agree, you should probably compute the surface of each cell contributing to a given image pixel (which is very tricky if you want to do it exactly)

[mreineck]

Thanks for linking the poster, there's a lot of very helpful information! I like your X/Y coordinate system a lot :)

The four of three cells are always on 2 different rings (2 strait lines of slope = -1 in the rotated frame), except if the pointing is abose the first ring, or under the last one (poles).

I must be missing something then. Imagine a location somewhere in between base pixels 4 and 5; I would have expected that you use pixels 0, 4, 5 and 8 for the interpolation, which lie on three different rings.

[fxpineau]

Oups! You are right, I am using 3 rings!
Is there any interest in using only 2?
Talking to Pierre Fernique, Aladin also uses 3 rings.

[mreineck]

I guess the consensus is that for visualization purposes it shouldn't really matter, although you might see strange features at the seams between base pixels at low Nside.

Personally I'm perfectly happy if we all agree that interpolation on the Healpix grid should not be used for "serious" work.

[fxpineau]

For information, I am working on a RUST version of the library.

Here the very first version of the Rust implementation: https://github.com/cds-astro/cds-healpix-rust

[bmatthieu3]

I worked during the past days to interface the Rust code from fx's healpix implementation https://github.com/cds-astro/cds-healpix-rust with Python. I am using the CFFI python package for this purpose which allows python to call Rust functions as it would call C ones.

The code is available here : https://github.com/cds-astro/cds-healpix-python

For the moment I have deployed the library as a pip cdshealpix package for different architectures:

• windows 32/64 bit for python 2.7/3.4->3.7
• linux i638/x86_64 for python 2.7/3.4->3.7
• I plan on creating the wheels for MacOS too in the future.

This way, people can use pip to download the good package wheel without the need of having the Rust compiler except if they want to compile the library (python setup.py install) or if they want to contribute.

You can check the README for code examples and contributing.

[mreineck]

FYI: there is now an "official" C package containing the basic Healpix functions under the 3-clause BSD license; see
https://sourceforge.net/projects/healpix/files/healpix_bare_1.0/

It's not much, but I tried to make the available code as robust as possible, while still maintaining good performance and readability - so that it can serve as a starting point for ports.

[astrofrog]

Oh that's great about healpix_bare! It would be good to see if we can just bundle it and use that instead.

[mreineck]

It probably won't be sufficient. There is no neighbor and no disk query functionality.

[mreineck]

... but still, it could be a good start.

[cdeil]

I'm closing this old PR. It was just me learning a bit, I don't have time to continue and develop something from scratch.

It's still not decided what we do for astropy-healpix, but there's a concrete proposal now in #128, so let's continue discussion there. That proposal to use the Rust lib was already discussed in March 2019 at PyGamma19, see presentations here and here which contain a lot of information.

Thank you everyone for your work and the feedback here!