Repository: numpy/numpy
Pull request: MAINT: adjustments to test_ufunc_noncontigous
Fixes #14087.

For ufunc where results are floating point arrays allow for some fuss (currently 4*machine epsilon, up for debates) in comparing results evaluated on contiguous and non-contiguous arrays of the same values. 

@charris 

Problem statement:
BUG: `test_ufunc_noncontiguous` too restrictive
Numpy 1.17 has added `numpy.core.tests.test_ufunc.py::test_ufunc_noncontiguous`, which insists that result of ufunc applied to an array be bit-wise identical, using `assert_equal` for strided and non-strided loops.

These tests might fail (as they do in Intel patched NumPy), if different implementations of universal function loop are used depending on the stride.

For example:

```
>>> svml_res = np.arccosh(np.arange(1,7, dtype='f'))
>>> libm_res = np.array([np.arccosh(x) for x in np.arange(1,7, dtype='f')], 'f')
>>> svml_res == libm_res
array([ True,  True,  True,  True,  True, False])
>>> svml_res - libm_res
array([0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,
       0.0000000e+00, 2.3841858e-07], dtype=float32)
>>> (svml_res - libm_res)[-1]
2.3841858e-07
>>> np.finfo('f').eps
1.1920929e-07
``` 

The failure is at 2 epsilons of the type.

So this test essentially precludes using different implementations for different strides.

How do you feel about relaxing `np.testing.assert_equal` to `np.testing.assert_almost_equal` with absolute/relative tolerance set to a multiple of type's epsilon ?

N.B.: OpenCL 2.0 standard mandates of implementation to uniformly be withing several ULPs of the true value.

Hints from issue discussion:
The test should be changed. Note that we no longer recommend `assert_almost_equal` for these
```
    .. note:: 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.
```
Forgot to add that a PR would be welcome.
@charris Goes without saying. Will have it ready in a bit.

Commit message: MAINT: adjustments to test_ufunc_noncontigous

Patch to review:
diff --git a/numpy/core/tests/test_ufunc.py b/numpy/core/tests/test_ufunc.py
index 69fbc35e33d4..bb365aa76252 100644
--- a/numpy/core/tests/test_ufunc.py
+++ b/numpy/core/tests/test_ufunc.py
@@ -1933,4 +1933,13 @@ def test_ufunc_noncontiguous(ufunc):
             warnings.filterwarnings("always")
             res_c = ufunc(*args_c)
             res_n = ufunc(*args_n)
-        assert_equal(res_c, res_n)
+        dt = np.find_common_type([np.dtype(c) for c in out], [])
+        if np.issubdtype(dt, np.floating):
+            # for floating point results allow a small fuss in comparisons
+            # since different algorithms (libm vs. intrinsics) can be used
+            # for different input strides
+            res_eps = np.finfo(dt).eps
+            tol = 4*res_eps
+            assert_allclose(res_c, res_n, atol=tol, rtol=tol)
+        else:
+            assert_equal(res_c, res_n)


Code review comments:
diff_hunk: @@ -1933,4 +1933,13 @@ def test_ufunc_noncontiguous(ufunc):
             warnings.filterwarnings("always")
             res_c = ufunc(*args_c)
             res_n = ufunc(*args_n)
-        assert_equal(res_c, res_n)
+        dt = np.find_common_type([np.dtype(c) for c in out], [])
+        if np.issubdtype(dt, np.floating):
+            # for floating point results allow a small fuss in comparisons
+            # since different algorithms (libm vs. intrinsics) can be used
+            # for different input strides
+            res_eps = np.finfo(dt).eps
+            tol = 4*res_eps
+            assert_allclose(res_c, res_n, atol=tol, rtol=tol)
path: numpy/core/tests/test_ufunc.py
text: @user1:
It seems to me that `assert_array_almost_equal_nulp` is a better fit maybe? `nulp=4` seems pretty big, but I guess anything less can easily lead to false positives...

diff_hunk: @@ -1933,4 +1933,13 @@ def test_ufunc_noncontiguous(ufunc):
             warnings.filterwarnings("always")
             res_c = ufunc(*args_c)
             res_n = ufunc(*args_n)
-        assert_equal(res_c, res_n)
+        dt = np.find_common_type([np.dtype(c) for c in out], [])
path: numpy/core/tests/test_ufunc.py
text: @user1:
Hmmm, maybe 
```
if len(out) == 1:
    res_c = (res_c,)
    res_n = (res_n,)
```
and then, `for c, n in zip(res_c, res_n): ...` then you can avoid the whole common type magic?

@author:
I tried to implement this suggestion, as well as the one to use `nulp` variant only to discover that it does not like `NaN` in the array:

```
>>> import numpy as np
>>> x = np.full((5,), np.nan, 'd')
>>> np.testing.assert_array_almost_equal_nulp(x, x, nulp=2)
/tmp/numpy_upstream/numpy/testing/_private/utils.py:1566: RuntimeWarning: invalid value encountered in greater
  ref = nulp * np.spacing(np.where(ax > ay, ax, ay))
/tmp/numpy_upstream/numpy/testing/_private/utils.py:1567: RuntimeWarning: invalid value encountered in less_equal
  if not np.all(np.abs(x-y) <= ref):
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/tmp/git/numpy_upstream/numpy/testing/_private/utils.py", line 1573, in assert_array_almost_equal_nulp
    raise AssertionError(msg)
AssertionError: X and Y are not equal to 2 ULP (max is 0)
```

So I will perhaps stick with `assert_allclose`

@user2:
Sounds like that function needs fixing, I'll open an issue. I'm fine with assert_allclose for now.

@user1:
Should we use `atol=0` though, or does that fail for denormals, or maybe `atol` should be just much smaller? Oh well, maybe just leave it in any case, right...