ENH: add assert_array_max_ulp comparison function.

This new comparison raises an error if the number of representable
numbers between two arrays exceeds a tolerance.

2 files changed, 65 insertions, 1 deletions

diff --git a/numpy/testing/tests/test_utils.py b/numpy/testing/tests/test_utils.py
index 332a1df78..8d154be82 100644
--- a/numpy/testing/tests/test_utils.py
+++ b/numpy/testing/tests/test_utils.py
@@ -367,5 +367,58 @@ class TestSpacing(unittest.TestCase):
             x = np.array([1e-5, 1, 1000, 10500], dtype=dt)
             assert_array_almost_equal(spacing(x), ref[dt], decimal=dec)
 
+class TestULP(unittest.TestCase):
+    def test_equal(self):
+        x = np.random.randn(10)
+        assert_array_max_ulp(x, x, maxulp=0)
+
+    def test_single(self):
+        # Generate 1 + small deviation, check that adding eps gives a few UNL
+        x = np.ones(10).astype(np.float32)
+        x += 0.01 * np.random.randn(10).astype(np.float32)
+        eps = np.finfo(np.float32).eps
+        assert_array_max_ulp(x, x+eps, maxulp=20)
+
+    def test_double(self):
+        # Generate 1 + small deviation, check that adding eps gives a few UNL
+        x = np.ones(10).astype(np.float32)
+        x += 0.01 * np.random.randn(10).astype(np.float64)
+        eps = np.finfo(np.float64).eps
+        assert_array_max_ulp(x, x+eps, maxulp=200)
+
+    def test_inf(self):
+        for dt in [np.float32, np.float64]:
+            inf = np.array([np.inf]).astype(dt)
+            big = np.array([np.finfo(dt).max])
+            assert_array_max_ulp(inf, big, maxulp=200)
+
+    def test_nan(self):
+        # Test that nan is 'far' from small, tiny, inf, max and min
+        for dt in [np.float32, np.float64]:
+            if dt == np.float32:
+                maxulp = 1e6
+            else:
+                maxulp = 1e12
+            inf = np.array([np.inf]).astype(dt)
+            nan = np.array([np.nan]).astype(dt)
+            big = np.array([np.finfo(dt).max])
+            tiny = np.array([np.finfo(dt).tiny])
+            zero = np.array([np.PZERO]).astype(dt)
+            nzero = np.array([np.NZERO]).astype(dt)
+            self.failUnlessRaises(AssertionError,
+                                  lambda: assert_array_max_ulp(nan, inf,
+                                                               maxulp=maxulp))
+            self.failUnlessRaises(AssertionError,
+                                  lambda: assert_array_max_ulp(nan, big,
+                                                               maxulp=maxulp))
+            self.failUnlessRaises(AssertionError,
+                                  lambda: assert_array_max_ulp(nan, tiny,
+                                                               maxulp=maxulp))
+            self.failUnlessRaises(AssertionError,
+                                  lambda: assert_array_max_ulp(nan, zero,
+                                                               maxulp=maxulp))
+            self.failUnlessRaises(AssertionError,
+                                  lambda: assert_array_max_ulp(nan, nzero,
+                                                               maxulp=maxulp))
 if __name__ == '__main__':
     run_module_suite()
diff --git a/numpy/testing/utils.py b/numpy/testing/utils.py
index 3e789272e..c5be4c49d 100644
--- a/numpy/testing/utils.py
+++ b/numpy/testing/utils.py
@@ -14,7 +14,8 @@ __all__ = ['assert_equal', 'assert_almost_equal','assert_approx_equal',
            'assert_array_almost_equal', 'assert_raises', 'build_err_msg',
            'decorate_methods', 'jiffies', 'memusage', 'print_assert_equal',
            'raises', 'rand', 'rundocs', 'runstring', 'verbose', 'measure',
-           'assert_', 'spacing', 'assert_array_almost_equal_nulp']
+           'assert_', 'spacing', 'assert_array_almost_equal_nulp',
+           'assert_array_max_ulp']
 
 verbose = 0
 
@@ -1110,6 +1111,16 @@ def assert_array_almost_equal_nulp(x, y, nulp=1):
         raise AssertionError("X and Y are not equal to %d ULP "\
                              "(max is %d)" % (nulp, max_nulp))
 
+def assert_array_max_ulp(a, b, maxulp=1, dtype=None):
+    """Given two arrays a and b, check that every item differs in at most N
+    Unit in the Last Place."""
+    import numpy as np
+    ret = nulp_diff(a, b, dtype)
+    if not np.all(ret <= maxulp):
+        raise AssertionError("Arrays are not almost equal up to %d ULP" % \
+                             maxulp)
+    return ret
+
 def nulp_diff(x, y, dtype=None):
     """For each item in x and y, eeturn the number of representable floating
     points between them.