Test failures with recent numpy versions

[avalentino]

During a rebuild of the Debian archive we had 4 test failures. The tests were previously working and we thing that the issue could related to update of numpy (currently 1.23.5) or some of the other dependencies. Please find the relevant part of the log below.

I see that the test_reverse_terra_with_zero_iterations has been already fixed in git. For other test I see that the error only exceeds slightly the tolerance so I think that it should be fine to relax a little bit. Some guidance on how to solve/workaround the issue would be appreciated. If considered useful, I could open a PR with the proposed fix.

> ============================= test session starts ==============================
> platform linux -- Python 3.10.9, pytest-7.2.0, pluggy-1.0.0+repack
> rootdir: /<<PKGBUILDDIR>>
> plugins: astropy-header-0.2.2, doctestplus-0.12.1, mock-3.8.2, filter-subpackage-0.1.1, cov-4.0.0, arraydiff-0.5.0, hypothesis-6.60.0, astropy-0.10.0, openfiles-0.5.0, remotedata-0.3.3
> collected 521 items / 226 deselected / 295 selected
> 
> skyfield/tests/test_against_novas.py ................................... [ 11%]
> ........................................                                 [ 25%]
> skyfield/tests/test_api.py .                                             [ 25%]
> skyfield/tests/test_broadcastability.py .                                [ 26%]
> skyfield/tests/test_constellations.py .                                  [ 26%]
> skyfield/tests/test_curves.py .                                          [ 26%]
> skyfield/tests/test_data_iers.py ..                                      [ 27%]
> skyfield/tests/test_earth_satellites.py ......                           [ 29%]
> skyfield/tests/test_earthlib.py F                                        [ 29%]
> skyfield/tests/test_elementslib.py ......F.F....                         [ 34%]
> skyfield/tests/test_frames.py ....                                       [ 35%]
> skyfield/tests/test_functions.py .....                                   [ 37%]
> skyfield/tests/test_geometry.py .                                        [ 37%]
> skyfield/tests/test_io.py .......                                        [ 40%]
> skyfield/tests/test_io_parsing.py ........                               [ 42%]
> skyfield/tests/test_keplerlib.py ...................                     [ 49%]
> skyfield/tests/test_positions.py .........F..                            [ 53%]
> skyfield/tests/test_satellite_events.py ..                               [ 53%]
> skyfield/tests/test_searchlib.py ......                                  [ 55%]
> skyfield/tests/test_strs_and_reprs.py ....                               [ 57%]
> skyfield/tests/test_text_pck.py ..                                       [ 57%]
> skyfield/tests/test_timelib.py ...........s...................s......... [ 71%]
> ................................                                         [ 82%]
> skyfield/tests/test_topos.py ...................................         [ 94%]
> skyfield/tests/test_trigonometry.py .                                    [ 94%]
> skyfield/tests/test_units.py ...............                             [100%]
> 
> =================================== FAILURES ===================================
> ___________________ test_reverse_terra_with_zero_iterations ____________________
> 
>     def test_reverse_terra_with_zero_iterations():
>         # With zero iterations, should return "geocentric" rather than
>         # "geodetic" (="correct") longitude and latitude.
>         lat, lon, elevation = reverse_terra(array([1, 0, 1]), 0, iterations=0)
>         assert abs(lat - tau / 8) < 1e-16
>         assert lon == 0.0
> >       assert abs(elevation - (AU_M * sqrt(2) - ERAD)) < 1e-16
> E       assert 6.103515625e-05 < 1e-16
> E        +  where 6.103515625e-05 = abs((211556959509.4766 - ((149597870700 * 1.4142135623730951) - 6378136.6)))
> E        +    where 1.4142135623730951 = sqrt(2)
> 
> skyfield/tests/test_earthlib.py:10: AssertionError
> ________________________________ test_parabolic ________________________________
> 
> ts = <skyfield.timelib.Timescale object at 0x7fbaecba04f0>
> 
>     def test_parabolic(ts):
>         e = 1
>         check_orbit(300000, e, angles1[:4], 0, 4, 3, ts)
>         check_orbit(300000, e, 2, angles1, 4, 3, ts)
> >       check_orbit(300000, e, 2, 3, angles1, 3, ts)
> 
> skyfield/tests/test_elementslib.py:512: 
> _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
> skyfield/tests/test_elementslib.py:428: in check_orbit
>     compare(elements.semi_latus_rectum.km, p, 1e-9)
> _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
> 
> value = array([300000., 300000., 300000., 300000., 300000., 300000., 300000.,
>        300000.])
> expected_value = 300000, epsilon = 1e-09, mod = False
> 
>     def compare(value, expected_value, epsilon, mod=False):
>         """Compares value to expected value, and works if one or both are arrays.
>     
>         Also allows epsilon to be an array.
>     
>         If mod==True, then compare(0, tau, 0) is True.
>         """
>         if mod:
>             diff = normpi(value - expected_value)
>         else:
>             diff = value - expected_value
>     
>         if hasattr(value, '__len__') or hasattr(expected_value, '__len__'):
>             if hasattr(epsilon, '__len__'):
>                 assert (abs(diff) <= epsilon).all()
>             else:
> >               assert max(abs(diff)) <= epsilon
> E               assert 1.0477378964424133e-09 <= 1e-09
> E                +  where 1.0477378964424133e-09 = max(array([0.00000000e+00, 4.65661287e-10, 0.00000000e+00, 1.04773790e-09,\n       0.00000000e+00, 1.16415322e-10, 1.74622983e-10, 2.91038305e-10]))
> E                +    where array([0.00000000e+00, 4.65661287e-10, 0.00000000e+00, 1.04773790e-09,\n       0.00000000e+00, 1.16415322e-10, 1.74622983e-10, 2.91038305e-10]) = abs(array([ 0.00000000e+00, -4.65661287e-10,  0.00000000e+00, -1.04773790e-09,\n        0.00000000e+00, -1.16415322e-10,  1.74622983e-10,  2.91038305e-10]))
> 
> skyfield/tests/test_elementslib.py:51: AssertionError
> _____________________________ test_parabolic_polar _____________________________
> 
> ts = <skyfield.timelib.Timescale object at 0x7fbaecdfdde0>
> 
>     def test_parabolic_polar(ts):
>         e = 1
>         i = pi/2
>         check_orbit(300000, e, i, angles1, 2, 3, ts)
> >       check_orbit(300000, e, i, 1, angles1, 3, ts)
> 
> skyfield/tests/test_elementslib.py:539: 
> _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
> skyfield/tests/test_elementslib.py:428: in check_orbit
>     compare(elements.semi_latus_rectum.km, p, 1e-9)
> _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
> 
> value = array([300000., 300000., 300000., 300000., 300000., 300000., 300000.,
>        300000.])
> expected_value = 300000, epsilon = 1e-09, mod = False
> 
>     def compare(value, expected_value, epsilon, mod=False):
>         """Compares value to expected value, and works if one or both are arrays.
>     
>         Also allows epsilon to be an array.
>     
>         If mod==True, then compare(0, tau, 0) is True.
>         """
>         if mod:
>             diff = normpi(value - expected_value)
>         else:
>             diff = value - expected_value
>     
>         if hasattr(value, '__len__') or hasattr(expected_value, '__len__'):
>             if hasattr(epsilon, '__len__'):
>                 assert (abs(diff) <= epsilon).all()
>             else:
> >               assert max(abs(diff)) <= epsilon
> E               assert 1.57160684466362e-09 <= 1e-09
> E                +  where 1.57160684466362e-09 = max(array([1.16415322e-10, 0.00000000e+00, 1.74622983e-10, 1.57160684e-09,\n       0.00000000e+00, 0.00000000e+00, 1.74622983e-10, 1.16415322e-10]))
> E                +    where array([1.16415322e-10, 0.00000000e+00, 1.74622983e-10, 1.57160684e-09,\n       0.00000000e+00, 0.00000000e+00, 1.74622983e-10, 1.16415322e-10]) = abs(array([-1.16415322e-10,  0.00000000e+00,  1.74622983e-10,  1.57160684e-09,\n        0.00000000e+00,  0.00000000e+00,  1.74622983e-10, -1.16415322e-10]))
> 
> skyfield/tests/test_elementslib.py:51: AssertionError
> ________________________ test_velocity_in_ITRF_to_GCRS2 ________________________
> 
>     def test_velocity_in_ITRF_to_GCRS2():
>         # TODO: Get test working with these vectors too, showing it works
>         # with a non-zero velocity vector, but in that case the test will
>         # have to be fancier in how it corrects.
>         # r = np.array([(1, 0, 0), (1, 1 / DAY_S, 0)]).T
>         # v = np.array([(0, 1, 0), (0, 1, 0)]).T
>     
>         ts = api.load.timescale()
>         t = ts.utc(2020, 7, 17, 8, 51, [0, 1])
>         r = np.array([(1, 0, 0), (1, 0, 0)]).T
>         v = np.array([(0, 0, 0), (0, 0, 0)]).T
>     
>         r, v = ITRF_to_GCRS2(t, r, v, True)
>     
>         # Rotate back to equinox-of-date before applying correction.
>         r = mxv(t.M, r)
>         v = mxv(t.M, v)
>     
>         r0, r1 = r.T
>         v0 = v[:,0]
>     
>         # Apply a correction: the instantaneous velocity does not in fact
>         # carry the position in a straight line, but in an arc around the
>         # origin; so use trigonometry to move the destination point to where
>         # linear motion would have carried it.
>         angvel = (t.gast[1] - t.gast[0]) / 24.0 * tau
>         r1 = mxv(rot_z(np.arctan(angvel) - angvel), r1)
>         r1 *= np.sqrt(1 + angvel*angvel)
>     
>         actual_motion = r1 - r0
>         predicted_motion = v0 / DAY_S
>     
>         relative_error = (length_of(actual_motion - predicted_motion)
>                           / length_of(actual_motion))
>     
>         acceptable_error = 4e-12
> >       assert relative_error < acceptable_error
> E       assert 9.700550724172958e-12 < 4e-12
> 
> skyfield/tests/test_positions.py:188: AssertionError
> =========================== short test summary info ============================
> FAILED skyfield/tests/test_earthlib.py::test_reverse_terra_with_zero_iterations
> FAILED skyfield/tests/test_elementslib.py::test_parabolic - assert 1.04773789...
> FAILED skyfield/tests/test_elementslib.py::test_parabolic_polar - assert 1.57...
> FAILED skyfield/tests/test_positions.py::test_velocity_in_ITRF_to_GCRS2 - ass...
> =========== 4 failed, 289 passed, 2 skipped, 226 deselected in 5.08s ===========

[brandon-rhodes]

Good morning, @avalentino! I have had a busy few months, but am now reading back through my open Skyfield issues, and I see that I never replied to this one—probably because I didn't have a way to reproduce the problem locally, so I wasn't sure how to get started answering.

I wish that I knew how to improve Skyfield's math to recover the precision it's customarily had with floating point, but for the moment you have my go-ahead to try submitting a pull request that simply relaxes these limits—as you note, these aren't hugely different numbers.

[avalentino]

Thanks @brandon-rhodes, just opened PR #844

[brandon-rhodes]

@avalentino — Two quick notes: first, I think we'll need to adjust the commit message a bit, because it's not a NumPy problem; I was just able to run the tests successfully against NumPy 1.24.2 here on my laptop. My guess would be that the actual problem is broken (by which I mean non-IEEE) floating point on the platform that Debian uses for its test cluster — do you have any idea whether the platform is an unusual processor type? Test failures like this always seem the result of running floating point tests on the cloud.

[brandon-rhodes]

Oh, and, my other note: could you revert the _IERS2 change? I should really do that as a separate diff.

[avalentino]

Dear @brandon-rhodes OK I can remove the commit regarding _IERS2. Regarding numpy I would like to clarify that right now I'm running tests on my ubuntu 22.10 laptop.

I have performed tests with several configurations:

• builtin Python 3.10.7 + builtin numpy 1.21.5 - OK
• builtin Python 3.10.7 + numpy 1.23.5 from PyPI - ERROR
• builtin Python 3.10.7 + numpy 1.24.2 from PyPI - ERROR
• builtin Python 3.10.7 + numpy 1.21.5 from PyPI - OK
• conda Python 3.11.2 + numpy 1.23.5 from PyPI - ERROR
• conda Python 3.11.2 + numpy 1.24.2 from PyPI - ERROR

I agree that the core problem is in non-IEEE floating point. By the way it seems to me that the problem can be triggered changing the numpy version. From my tests the problem seems to be present only in newer versions of numpy but the problem could be rather in the build flags used for numpy or the libraries used (openblas vs mkl vs whatever)

[avalentino]

None of the commit messages refers to numpy. Please feel free to change the PR title in the way you think is more appropriate.

[brandon-rhodes]

@avalentino — Ah, yes, I mean to reference the PR title, not the commit messages, thanks!

So, thank you for providing such detailed test results! I had thought that the failures were only out on a build server somewhere; I had no idea you were able to reproduce them locally, as that's something I've never been able to do successfully! Since the tests all pass on my own laptop with numpy==1.24.2, I had assumed they would pass anywhere, so it's an unhappy result that the exact same NumPy version fails for you.

Someday I should take a couple of weeks and try to track down what's going on at the machine level with these tests. But I suspect that I won't be tackling that odyssey soon.

What processor does your machine use? Mine looks like a Intel(R) Core(TM) i5-3427U CPU @ 1.80GHz.

On a lighter note, did you see the following blog post from last year? 🙂

https://moyix.blogspot.com/2022/09/someones-been-messing-with-my-subnormals.html

In any event, I think I'll go ahead and merge this as soon as I can confirm it won't change the IERS URL. Thanks!

[avalentino]

My CPU is "Intel(R) Core(TM) i7-1065G7 CPU @ 1.30GHz"