BUG: numpy.test() fails at test_linear_interpolation_formula_symmetric

[GideonBear]

Describe the issue:

When I run numpy.test() a failure is encountered in test_linear_interpolation_formula_symmetric. I don't know if this is actually a bug, but I thought I should add this as an issue. I've run the test twice now and it gives the same result.
I've searched for this on GitHub, Stack Overflow and Google, but got no results.
OS is Linux Mint 19.3 Cinnamon.

Reproduce the code example:

import numpy
numpy.test()

Error message:

============================================================================================== FAILURES ==============================================================================================
________________________________________________________________________ TestLerp.test_linear_interpolation_formula_symmetric ________________________________________________________________________

self = <numpy.lib.tests.test_function_base.TestLerp object at 0x7fe580065190>

    @hypothesis.given(t=st.floats(allow_nan=False, allow_infinity=False,
>                                 min_value=0, max_value=1),
                      a=st.floats(allow_nan=False, allow_infinity=False,
                                  min_value=-1e300, max_value=1e300),
                      b=st.floats(allow_nan=False, allow_infinity=False,
                                  min_value=-1e300, max_value=1e300))

f          = <function given.<locals>.run_test_as_given.<locals>.wrapped_test at 0x7fe5802cf790>
self       = <numpy.lib.tests.test_function_base.TestLerp object at 0x7fe580065190>

venv/lib/python3.8/site-packages/numpy/lib/tests/test_function_base.py:3568: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

self = <numpy.lib.tests.test_function_base.TestLerp object at 0x7fe580065190>, t = 0.5, a = 0.2500000000000003, b = 1.0

    @hypothesis.given(t=st.floats(allow_nan=False, allow_infinity=False,
                                  min_value=0, max_value=1),
                      a=st.floats(allow_nan=False, allow_infinity=False,
                                  min_value=-1e300, max_value=1e300),
                      b=st.floats(allow_nan=False, allow_infinity=False,
                                  min_value=-1e300, max_value=1e300))
    def test_linear_interpolation_formula_symmetric(self, t, a, b):
        # double subtraction is needed to remove the extra precision of t < 0.5
        left = nfb._lerp(a, b, 1 - (1 - t))
        right = nfb._lerp(b, a, 1 - t)
>       assert left == right
E       assert 0.6250000000000001 == 0.6250000000000002

a          = 0.2500000000000003
b          = 1.0
left       = 0.6250000000000001
right      = 0.6250000000000002
self       = <numpy.lib.tests.test_function_base.TestLerp object at 0x7fe580065190>
t          = 0.5

venv/lib/python3.8/site-packages/numpy/lib/tests/test_function_base.py:3577: AssertionError
--------------------------------------------------------------------------------------------- Hypothesis ---------------------------------------------------------------------------------------------
Falsifying example: test_linear_interpolation_formula_symmetric(
    t=0.5,
    a=0.2500000000000003,
    b=1.0,
    self=<numpy.lib.tests.test_function_base.TestLerp at 0x7fe580065190>,
)

You can reproduce this example by temporarily adding @reproduce_failure('6.53.0', b'AXicY2BgOMAABwekGNAAIwAmtgGc') as a decorator on your test case

NumPy/Python version information:

1.23.1 3.8.0 (default, Dec 9 2021, 17:53:27)
[GCC 8.4.0]

[mattip]

Did you install NumPy from a wheel/conda or did you build it yourself?

[GideonBear]

Did you install NumPy from a wheel/conda or did you build it yourself?

I installed it via a wheel: numpy-1.23.1-cp38-cp38-manylinux_2_17_x86_64.manylinux2014_x86_64.whl

[mattip]

Huh. I can reproduce on 1.23.1 by updating hypothesis to 6.53.0 and adding this decorator to the test.
@hypothesis.reproduce_failure('6.53.0', b'AXicY2BgOMAABwekGNAAIwAmtgGc')

[mattip]

The test is checking that linear interpolation is equivalent when flipped around and seems to be off by one eps:

    import numpy.lib.function_base as nfb
    t=0.5
    a=0.2500000000000003
    b=1.0
    left = nfb._lerp(a, b, 1 - (1 - t))
    right = nfb._lerp(b, a, 1 - t)
    assert left == right

[palbarta]

Test setup: numpy 1.23.1 | Python 3.10.6 | macOS 12.4

The test run passes with hypothesis==6.48.0:
runtests.py -t numpy/lib/tests/test_function_base.py::TestLerp::test_linear_interpolation_formula_symmetric

The test failure appears with the hypothesis==6.48.1 upgrade. They seem to have changed the float testing strategy: HypothesisWorks/hypothesis#2701

The new strategy uncovered a float precision related mismatch as mentioned above by @mattip

t=0.5
a=0.2500000000000003
b=1.0

The lerp function:

def _lerp(a, b, t, out=None):  
    diff_b_a = subtract(b, a)
    # asanyarray is a stop-gap until gh-13105
    lerp_interpolation = asanyarray(add(a, diff_b_a * t, out=out))
    subtract(b, diff_b_a * (1 - t), out=lerp_interpolation, where=t >= 0.5)
    if lerp_interpolation.ndim == 0 and out is None:
        lerp_interpolation = lerp_interpolation[()]  # unpack 0d arrays
    return lerp_interpolation

For the failing test case (t = 0.5) it's effectively:

diff_b_a = subtract(b, a)
lerp_interpolation = subtract(b, diff_b_a * (1 - t))

The real calculation vs. the double precision calculation:

diff_b_a = subtract(b, a)
# REAL     1.0000000000000000 -
#          0.2500000000000003 =
#          0.7499999999999997
# BIN REPR 1.000000000000000000 -
#          0.250000000000000278 = | error = -2.2e-17
#          0.749999999999999778   | error = +7.8e-17

lerp_interpolation = subtract(b, diff_b_a * (1 - t))
# REAL     1.00000000000000000 -
#          0.37499999999999985 =
#          0.62500000000000015
# BIN REPR 1.000000000000000000 -
#          0.374999999999999889 = | error = +3.9e-17
#          0.625000000000000111   | error = -3.9e-17

The double precision error accumulates enough to flip the last bit on the symmetric run:

diff_b_a = subtract(b, a)
# REAL     0.2500000000000003 -
#          1.0000000000000000 =
#         -0.7499999999999997
# BIN REPR 0.250000000000000278 - | error = -2.2e-17
#          1.000000000000000000 = 
#         -0.749999999999999778   | error = -7.8e-17

lerp_interpolation = subtract(b, diff_b_a * (1 - t))
# REAL     0.2500000000000003 -
#         -0.37499999999999985 =
#          0.62500000000000015
# BIN REPR 0.250000000000000278 - | error = -2.2e-17
#         -0.374999999999999889 = | error = -3.9e-17
#          0.625000000000000222   | error = +7.2e-17

This does not seem like a bug to me in _lerp, but something natural for floating point comparisons, so my suggested fix would be to replace the assert with assert_allclose in the unit test:

diff --git a/numpy/lib/tests/test_function_base.py b/numpy/lib/tests/test_function_base.py
index 8457551ca..56081b4f0 100644
--- a/numpy/lib/tests/test_function_base.py
+++ b/numpy/lib/tests/test_function_base.py
@@ -3579,7 +3579,7 @@ def test_linear_interpolation_formula_symmetric(self, t, a, b):
         # double subtraction is needed to remove the extra precision of t < 0.5
         left = nfb._lerp(a, b, 1 - (1 - t))
         right = nfb._lerp(b, a, 1 - t)
-        assert left == right
+        assert_allclose(left, right)

[charris]

This is also failing for the wheel builds.

EDIT: But not on windows. That suggests different rounding strategies.

[mattip]

@palbarta indeed changing to assert_allclose would make sense here. Want to submit a PR?

[palbarta]

@mattip PR submitted

[InessaPawson]

The issue was fixed by #22106. Thank you, @palbarta!

[palbarta]

@InessaPawson Thank you for taking my fix, it was my pleasure.