test_position_of_radec failure on 32bit arm and s390x

[avalentino]

Recently I packaged skyfield for debian and I'm having some failures of the test_position_of_radec test on some architectures:

• ARMEL

=== python3.10 ===
============================= test session starts ==============================
platform linux -- Python 3.10.6, pytest-7.1.2, pluggy-1.0.0+repack -- /usr/bin/python3.10
cachedir: .pytest_cache
hypothesis profile 'default' -> database=DirectoryBasedExampleDatabase('/tmp/autopkgtest-lxc.w7ek_adm/downtmp/autopkgtest_tmp/.hypothesis/examples')
rootdir: /tmp/autopkgtest-lxc.w7ek_adm/downtmp/autopkgtest_tmp
plugins: arraydiff-0.5.0, astropy-0.10.0, remotedata-0.3.3, hypothesis-6.36.0, filter-subpackage-0.1.1, openfiles-0.5.0, cov-3.0.0, mock-3.8.2, doctestplus-0.12.0, astropy-header-0.2.1
collecting ... collected 518 items / 225 deselected / 293 selected

[...]

=================================== FAILURES ===================================
____________________________ test_position_of_radec ____________________________

    def test_position_of_radec():
        epsilon = _GIGAPARSEC_AU * 1e-16
    
        p = api.position_of_radec(0, 0)
        assert length_of(p.position.au - [_GIGAPARSEC_AU, 0, 0]) < epsilon
    
        p = api.position_of_radec(6, 0)
        assert length_of(p.position.au - [0, _GIGAPARSEC_AU, 0]) < epsilon
    
        epsilon = 2e-16
    
        p = api.position_of_radec(12, 90, 2)
        assert length_of(p.position.au - [0, 0, 2]) < epsilon
    
        p = api.position_of_radec(12, 90, distance_au=2)
        assert length_of(p.position.au - [0, 0, 2]) < epsilon
    
        ts = api.load.timescale()
        epoch = ts.tt_jd(api.B1950)
        p = api.position_of_radec(0, 0, 1, epoch=epoch)
        assert length_of(p.position.au - [1, 0, 0]) > 1e-16
        ra, dec, distance = p.radec(epoch=epoch)
        assert abs(ra.hours) < 1e-12
        assert abs(dec.degrees) < 1e-12
>       assert abs(distance.au - 1) < 1e-16
E       assert 1.1102230246251565e-16 < 1e-16
E        +  where 1.1102230246251565e-16 = abs((0.9999999999999999 - 1))
E        +    where 0.9999999999999999 = <Distance 1.0 au>.au

/usr/lib/python3/dist-packages/skyfield/tests/test_positions.py:142: AssertionError
=========================== short test summary info ============================
FAILED tests/test_positions.py::test_position_of_radec - assert 1.11022302462...
========== 1 failed, 290 passed, 2 skipped, 225 deselected in 19.62s ===========

• ARMHF

=== python3.10 ===
============================= test session starts ==============================
platform linux -- Python 3.10.6, pytest-7.1.2, pluggy-1.0.0+repack -- /usr/bin/python3.10
cachedir: .pytest_cache
hypothesis profile 'default' -> database=DirectoryBasedExampleDatabase('/tmp/autopkgtest-lxc.2z77qe2u/downtmp/autopkgtest_tmp/.hypothesis/examples')
rootdir: /tmp/autopkgtest-lxc.2z77qe2u/downtmp/autopkgtest_tmp
plugins: mock-3.8.2, astropy-header-0.2.1, doctestplus-0.12.0, filter-subpackage-0.1.1, cov-3.0.0, remotedata-0.3.3, openfiles-0.5.0, arraydiff-0.5.0, astropy-0.10.0, hypothesis-6.36.0
collecting ... collected 518 items / 225 deselected / 293 selected

[...]

=================================== FAILURES ===================================
____________________________ test_position_of_radec ____________________________

    def test_position_of_radec():
        epsilon = _GIGAPARSEC_AU * 1e-16
    
        p = api.position_of_radec(0, 0)
        assert length_of(p.position.au - [_GIGAPARSEC_AU, 0, 0]) < epsilon
    
        p = api.position_of_radec(6, 0)
        assert length_of(p.position.au - [0, _GIGAPARSEC_AU, 0]) < epsilon
    
        epsilon = 2e-16
    
        p = api.position_of_radec(12, 90, 2)
        assert length_of(p.position.au - [0, 0, 2]) < epsilon
    
        p = api.position_of_radec(12, 90, distance_au=2)
        assert length_of(p.position.au - [0, 0, 2]) < epsilon
    
        ts = api.load.timescale()
        epoch = ts.tt_jd(api.B1950)
        p = api.position_of_radec(0, 0, 1, epoch=epoch)
        assert length_of(p.position.au - [1, 0, 0]) > 1e-16
        ra, dec, distance = p.radec(epoch=epoch)
        assert abs(ra.hours) < 1e-12
        assert abs(dec.degrees) < 1e-12
>       assert abs(distance.au - 1) < 1e-16
E       assert 2.220446049250313e-16 < 1e-16
E        +  where 2.220446049250313e-16 = abs((1.0000000000000002 - 1))
E        +    where 1.0000000000000002 = <Distance 1.0 au>.au

/usr/lib/python3/dist-packages/skyfield/tests/test_positions.py:142: AssertionError
=========================== short test summary info ============================
FAILED tests/test_positions.py::test_position_of_radec - assert 2.22044604925...
========== 1 failed, 290 passed, 2 skipped, 225 deselected in 11.31s ===========

• S390X

=== python3.10 ===
============================= test session starts ==============================
platform linux -- Python 3.10.6, pytest-7.1.2, pluggy-1.0.0+repack -- /usr/bin/python3.10
cachedir: .pytest_cache
hypothesis profile 'default' -> database=DirectoryBasedExampleDatabase('/tmp/autopkgtest-lxc.2gzq6va4/downtmp/autopkgtest_tmp/.hypothesis/examples')
rootdir: /tmp/autopkgtest-lxc.2gzq6va4/downtmp/autopkgtest_tmp
plugins: doctestplus-0.12.0, arraydiff-0.5.0, cov-3.0.0, astropy-header-0.2.1, openfiles-0.5.0, remotedata-0.3.3, astropy-0.10.0, mock-3.8.2, hypothesis-6.36.0, filter-subpackage-0.1.1
collecting ... collected 518 items / 225 deselected / 293 selected

[...]

=================================== FAILURES ===================================
____________________________ test_position_of_radec ____________________________

    def test_position_of_radec():
        epsilon = _GIGAPARSEC_AU * 1e-16
    
        p = api.position_of_radec(0, 0)
        assert length_of(p.position.au - [_GIGAPARSEC_AU, 0, 0]) < epsilon
    
        p = api.position_of_radec(6, 0)
        assert length_of(p.position.au - [0, _GIGAPARSEC_AU, 0]) < epsilon
    
        epsilon = 2e-16
    
        p = api.position_of_radec(12, 90, 2)
        assert length_of(p.position.au - [0, 0, 2]) < epsilon
    
        p = api.position_of_radec(12, 90, distance_au=2)
        assert length_of(p.position.au - [0, 0, 2]) < epsilon
    
        ts = api.load.timescale()
        epoch = ts.tt_jd(api.B1950)
        p = api.position_of_radec(0, 0, 1, epoch=epoch)
        assert length_of(p.position.au - [1, 0, 0]) > 1e-16
        ra, dec, distance = p.radec(epoch=epoch)
        assert abs(ra.hours) < 1e-12
        assert abs(dec.degrees) < 1e-12
>       assert abs(distance.au - 1) < 1e-16
E       assert 2.220446049250313e-16 < 1e-16
E        +  where 2.220446049250313e-16 = abs((1.0000000000000002 - 1))
E        +    where 1.0000000000000002 = <Distance 1.0 au>.au

/usr/lib/python3/dist-packages/skyfield/tests/test_positions.py:142: AssertionError
=========================== short test summary info ============================
FAILED tests/test_positions.py::test_position_of_radec - assert 2.22044604925...
=========== 1 failed, 290 passed, 2 skipped, 225 deselected in 4.12s ===========

It seems that it is mostly a matter of floating point accuracy.
Does it make sense to relax the error threshold form 1e-16 to e.g. 3e-16?

[rmathar]

The native IEEE double precision is nominally 2^-53, which is 1.1e-16. If one is lucky and all information is kept in the registers (which may have e.g. 80 bits like it had in the the M68k series), and one does not use compiler options that enforce IEEE compatibility, this sort of relative accuracy can be achieved. If the processor is moving all the floating point numbers in and out of memory, because it does not have many registers, a precision of 1.e-16 is almost impossible to keep. Backports that emulate 32-bit systems are also likely to get exactly a 64-bit lookalike output (I haven't looked into the armel specs, though). So to keep tests compatible with all sorts of host systems, I'm usually taking 1.e-15 as reasonable error margins in my coding.

[avalentino]

Thanks for the feedback @rmathar.
The patch currently applied in debian is at https://salsa.debian.org/debian-astro-team/skyfield/-/blob/master/debian/patches/0003-Fix-unittests-on-32bit-ARM.patch.
It uses 3e-16 as threshold and it seems to work at the moment.
Anyway I totally agree with you and I think that using 1e15 would be safer.

If all agree on it I could provide a small PR.

[brandon-rhodes]

I think I have a branch somewhere that reduces several test thresholds—the Anaconda packaging folks, I think, may have asked for it so the tests could pass on some other platforms they test? But I got bogged down this summer with other responsibilities, and haven't gotten back to it yet.

Based on branch names, my guess is that this is my most recent work:

https://github.com/skyfielders/python-skyfield/tree/fix-float

And this might have been an even earlier attempt?

https://github.com/skyfielders/python-skyfield/tree/precisions

I wish Python and NumPy would use IEEE compatibility options when compiling, it would make it easier to write programs that print the same output on different platforms!

[brandon-rhodes]

As we have recently addressed this problem in fac3cff, I am going to go ahead and close this issue now.