Allow for numerical precision in test

[llimeht]

On x86 hardware (Python 2.7.14 on Debian's i386 port) the test against a known problem fails due to numerical precision. The attached patch uses numpy's assert_almost_equal to test more robustly; one could reinvent an assert abs(a - b) / a < delta instead but since numpy is already used throughout, there's no harm in using it here too.

-> assert py == poly(px), "direct call to poly function fails"
(Pdb) py
array([ 13.18486645])
(Pdb) px
array([ 1.5])
(Pdb) poly(px)
array([ 13.18486645])
(Pdb) py - poly(px)
array([  1.77635684e-15])

[pkienzle]

Yes the results could very well differ. polyval uses Horner's method

(...(a[n] x + a[n-1]) x + ... + a[1])x + a[0] 

while wpolyfit.ci uses the more direct evaluation

a[n] x^n + ... + a[0]

From the numpy docs:

It is recommended to use one of assert_allclose, assert_array_almost_equal_nulp or assert_array_max_ulp instead of this function for more consistent floating point comparisons.

In this case assert_array_almost_equal_nulp would be the better choice, with nulp=8 allowing the last three bits to differ.

>>> from numpy.testing import assert_array_almost_equal_nulp as isclose_nulp
>>> isclose_nulp(1, 1+8*2**-52, 8)
>>> isclose_nulp(1, 1+9*2**-52, 8)
AssertionError: X and Y are not equal to 8 ULP (max is 9)
>>> isclose_nulp(2, 2+8*2**-51, 8)
>>> isclose_nulp(2, 2+9*2**-51, 8)
AssertionError: X and Y are not equal to 8 ULP (max is 9)

[llimeht]

Do I understand your explanation of why there might be differences as saying that it's not a matter of limited numerical precision in floats but of fundamentally different algorithms giving different values? My reading of assert_array_almost_equal_nulp is that it is for comparisons where the implementation details of the float storage is a problem, while relative tolerance tests (such as assert_allclose or assert_almost_equal) are for comparisons where it is not expected that the results would be identical if it were not for finite precision floats being used.

(Happy either way, just want to pick the right function for the job.)

[pkienzle]

Short answer: in this case assert_array_almost_equal_nulp would be the better choice, with nulp=8 allowing the last three bits to differ. Or even 4.

Looking a bit closer, this is a linear system, so effectively it is the same algorithm, though one code path goes through BLAS and the other is presumably cython (haven't looked at the polyval implementation). Depending on exact register usage and processor (Intel processors keep a guard bit on numbers retained in registers that can be discarded if the value is written to memory). Regardless, since the test has been passing on my machines using floating point equality (linux, windows and mac but all x86), that's a pretty good hint that the values should be the same except for the last bit or so.

The all_close test is using 7 digits, which is suitable for single precision, but is way too loose for this case. There are some cases where that much precision is lost compared to an analytic solution (e.g., optimizers may stop when df is less than 1e-7), but not here. In general, I will adjust a test to use as many bits as I can get away with, sometimes adjusting the algorithm to make the results more precise.

[llimeht]

Updated to use assert_array_almost_equal_nulp with nulp=8.