[0.7] test_vec2pix and test_vec2ang fail

[olebole]

The Debian package of astropy-healpix got bug report Debian#1026012, that the build time test now fails. The package was version 0.6, but0.7 shows the same behavior. Specifically, the following two tests fail:

_________________________________ test_vec2pix _________________________________

    @given(nside_pow=integers(0, 29), nest=booleans(),
>          x=floats(-1, 1, allow_nan=False, allow_infinity=False).filter(lambda x: abs(x) > 1e-10),
           y=floats(-1, 1, allow_nan=False, allow_infinity=False).filter(lambda y: abs(y) > 1e-10),
           z=floats(-1, 1, allow_nan=False, allow_infinity=False).filter(lambda z: abs(z) > 1e-10))

astropy_healpix/tests/test_healpy.py:116: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

nside_pow = 2, x = 1.0, y = 1.0, z = 0.5, nest = False

    @given(nside_pow=integers(0, 29), nest=booleans(),
           x=floats(-1, 1, allow_nan=False, allow_infinity=False).filter(lambda x: abs(x) > 1e-10),
           y=floats(-1, 1, allow_nan=False, allow_infinity=False).filter(lambda y: abs(y) > 1e-10),
           z=floats(-1, 1, allow_nan=False, allow_infinity=False).filter(lambda z: abs(z) > 1e-10))
    @settings(max_examples=2000, derandomize=True, deadline=None)
    def test_vec2pix(nside_pow, x, y, z, nest):
        nside = 2 ** nside_pow
        ipix1 = hp_compat.vec2pix(nside, x, y, z, nest=nest)
        ipix2 = hp.vec2pix(nside, x, y, z, nest=nest)
>       assert ipix1 == ipix2
E       assert 42 == 58
E       Falsifying example: test_vec2pix(
E           nside_pow=2, nest=False, x=1.0, y=1.0, z=0.5,
E       )

astropy_healpix/tests/test_healpy.py:124: AssertionError
_________________________________ test_vec2ang _________________________________

    @given(vectors=arrays(float, (3,), elements=floats(-1, 1)).filter(not_at_origin),
>          lonlat=booleans(), ndim=integers(0, 4))

astropy_healpix/tests/test_healpy.py:210: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

vectors = array([  1.00000000e+00,  -2.22044605e-16,   1.00000000e+00])
lonlat = False, ndim = 0

    @given(vectors=arrays(float, (3,), elements=floats(-1, 1)).filter(not_at_origin),
           lonlat=booleans(), ndim=integers(0, 4))
    @settings(max_examples=500, derandomize=True, deadline=None)
    def test_vec2ang(vectors, lonlat, ndim):
        vectors = np.broadcast_to(vectors, (2,) * ndim + (3,))
        theta1, phi1 = hp_compat.vec2ang(vectors, lonlat=lonlat)
        theta2, phi2 = hp.vec2ang(vectors, lonlat=lonlat)
        # Healpy sometimes returns NaNs for phi (somewhat incorrectly)
        phi2 = np.nan_to_num(phi2)
        assert_allclose(theta1, theta1, atol=1e-10)
>       assert_allclose(phi1, phi2, atol=1e-10)
E       AssertionError: 
E       Not equal to tolerance rtol=1e-07, atol=1e-10
E       
E       Mismatched elements: 1 / 1 (100%)
E       Max absolute difference: 6.283185307179586
E       Max relative difference: 1.0
E        x: array([ 0.])
E        y: array([ 6.283185])
E       Falsifying example: test_vec2ang(
E           vectors=array([  1.00000000e+00,  -2.22044605e-16,   1.00000000e+00]),
E           lonlat=False,
E           ndim=0,
E       )

astropy_healpix/tests/test_healpy.py:219: AssertionError

This is with

• Python 3.10.9, 3.11.1
• Astropy 5.1.1 (bug report), 5.2 (updated test)
• Numpy: 1.23.5
• Scipy: 1.8.1
• Matplotlib: 3.5.2

The package was reported to pass the integration testing for Astropy 5.2.

[roehling]

Appendix D of the HEALPix Introduction states:

Note that, due to finite precision of floating-point arithmetics and differences between
numerical libraries, the result of HEALPix functions like ang2pix (which converts the
angular coordinates of a point into the index of the pixel to which it belongs) may depend on
the underlying hardware, compilers, compiler options and linked libraries, if the requested
position is very close to (roughly 10^−15 radians or less) a pixel border. [...] This may be surprising when
testing apparently "simple" locations like the poles.

Maybe the vec2pix test would be more robust if it did not compare two different vec2pix implementations but checked if hp.pix2vec(hp_compat.vec2pix(...)) (and vice versa) ended up close to the original vector again?

[olebole]

@roehling astropy-healpix was originally introduced as a BSD-licensed alternative to healpy (which is GPL); so it IMO makes sense to make sure they both produce the same results. I'd however agree that this does not need to extend beyond FP arithmetic limits.

[roehling]

That's why I proposed the forward-backward round trip through both implementations, which effectively shows that they behave equivalently (up to FP precision).

[mreineck]

@roehling, if I understand it correctly, your proposed hp.pix2vec(hp_compat.vec2pix(...)) test will return a vector that's within half a pixel size of the input vector, and "pixel size" is not very nicely defined in Healpix. I'm not sure whether this kind of test would be very informative.

One could try something like hp.vec2pix(hp_compat.pix2vec(pix))==pix, but that only tests whether the functions are equivalent if provided with pixel centers...

[roehling]

@mreineck Given that the functions are numerically unstable if provided with vectors which are on (or very close to) a pixel border, I see two options:

1. Use your proposed test (which is probably better than mine, for the reasons you describe)
2. Make sure that all tested inputs are sufficiently far away from pixel borders