#147 ring2xyf differs between C and cxx versions, resulting in different pixel boundaries

[Anonymous]

ring2xyf differs between C and cxx versions, resulting in different pixel boundaries.

specifically this C ring2xyf:
https://github.com/healpy/healpixmirror/blob/a44dd367e50a038ea38f1916ff78c20607fde315/src/C/subs/chealpix.c#L185

and this cxx ring2xyf:
https://github.com/healpy/healpixmirror/blob/a44dd367e50a038ea38f1916ff78c20607fde315/src/cxx/Healpix_cxx/healpix_base.cc#L675

I was transcribing ring2xyf to pure python, and was wondering why some pixels had boundaries that didn't match up with healpy.boundaries(). I realized this was due to me transcribing the C ring2xyf function, and not the cxx version. I think this only affect pixels that have boundaries touching the longitude=0,90,180,270, where the pixel's top corner gets shifted over the wrong direction.

I've attached a plot of how this affected the pixel boundaries. Blue is the correct (cxx) boundary, while the C ring2xyf gives the orange boundary. Neighboring cxx pixels are shown in grey.

I dont know enough about the healpix code organization to know if that C function is still used somewhere, but I thought I'd bring it to people's attention.

[Martin Reinecke]

Thanks for the report!

To help me understand wat's going on, could you please describe how you get from the information provided by ring2xyf to the shape of pixel boundaries? It's not immediately obvious to me, since ring2xyf only deals with pixel centers. Are you working with highly supersampled pixels to trace the boundary?
If you could provide an example program demonstrating the problem that would of course be even better!

[Martin Reinecke]

... or even simpler: could you give us values for nside and pix, for which the C and C++ versions provide different results?

[Anonymous]

describe how you get from the information provided by ring2xyf to the shape of pixel boundaries

I don't really fully understand the whole algorithm, but the function calls are roughly nested like this:
healpy.boundaries()
boundaries_single()
healpix_base::boundaries()
pix2xyf() < might be a ufunc, or a pywrapper kludge
pix2xyf_ring() < might be a ufunc, or a pywrapper kludge
ring2xyf()

Are you working with highly supersampled pixels to trace the boundary?

I found it by transcribing the code from C to python, then I wrote a testing program (to check my work) to pick 100k random nsides and pixel numbers, and for each of those nside/pixel combinations, compare the boundary points between healpy (truth) and my transcribed boundary function. I kept getting a handful of random mismatched pixel boundaries, which led me down this rabbit hole.

In this case, I found the bug with nside=8, pixel=136, scheme=RING .

If you could provide an example program demonstrating the problem that would of course be even better!

I've attached a complete python script that

• requires healpy, matplotlib, astropy, and numpy to run, (ran on python 3.12.1)
• includes the transcribed C function and all other necessary stuff
• recreates the plot of the first bugged pixel I mentioned, saved as glitch_pixel.png

It additionally:

• does 100k random pixel comparisons at random nsides
• does a full scan of nsides 2,4,8,16 to search for bugged pixels
• plots all bugged pixels from the full scan, saved as bugged_pixels.png

could you give us values for nside and pix

From the random (not-exhaustive) scan, it found the following bugged RING pixels:
nside 4 : pixels = [24, 28, 32, 36, 152, 156, 160, 164]
nside 8 : pixels = [112, 120, 128, 136, 624, 632, 640, 648]
nside 16 : pixels = [480, 496, 528, 2544, 2560, 2576]
nside 32 : pixels = [1984, 10240, 10272]
nside 128 : pixels = [32896]

From the full scan (all pixels checked at these nsides), I found the following bugged RING pixels:
nside=2 has bugged pixels [4, 6, 8, 10, 36, 38, 40, 42]
nside=4 has bugged pixels [24, 28, 32, 36, 152, 156, 160, 164]
nside=8 has bugged pixels [112, 120, 128, 136, 624, 632, 640, 648]
nside=16 has bugged pixels [480, 496, 512, 528, 2528, 2544, 2560, 2576]

From the bugged_pixels.png, you can see that it only seems to affect 8 pixels. From their positioning, I think this code divergence happens at the intersection of:

• the left edge (in that plot) of each base pixel,
• when the code transitions from equatorial to polar regions.

Huge gotcha: I only transcribed the RING scheme code for all of the above, so I have no idea if this bug also affects the NESTed scheme.

[Martin Reinecke]

OK, I think I found he problem ... it happens in the translation to Python.
In ring2xyf there are integer divisions, for example:

iring = (ip/(4*nside_)) + nside_; /* counted from North pole */

int ifm = (iphi - ire/2 + nside_ -1) / nside_;
int ifp = (iphi - irm/2 + nside_ -1) / nside_;

When translating this to Python, it is mandatory to use the integer division operator "//", otherwise the results wil be wrong:

iring  = ( ip // ( 4 * nside ) ) + nside

ifm = ( iphi - ire//2 + nside - 1 ) // nside
ifp = ( iphi - irm//2 + nside - 1 ) // nside

Just converting the result to "int" won't do the trick in this case.