NF - add ulp function (cf. matlab eps function)

Function to return gap between value and nearest representable value of
same type.

Use ``ulp`` instead of ``eps` because the wikipedia article pointed out
the different meanings possible for eps:

http://en.wikipedia.org/wiki/Machine_epsilon#Variant_definitions

Author: matthew-brett

nibabel/casting.py

@@ -177,9 +177,10 @@ def type_info(np_type):
     info : dict
         with fields ``min`` (minimum value), ``max`` (maximum value), ``nexp``
         (exponent width), ``nmant`` (significand precision not including
-        implicit first digit) ``width`` (width in bytes). ``nexp``, ``nmant``
-        are None for integer types. Both ``min`` and ``max`` are of type
-        `np_type`.
+        implicit first digit), ``minexp`` (minimum exponent), ``maxexp``
+        (maximum exponent), ``width`` (width in bytes). (``nexp``, ``nmant``,
+        ``minexp``, ``maxexp``) are None for integer types. Both ``min`` and
+        ``max`` are of type `np_type`.
 
     Raises
     ------
@@ -188,9 +189,10 @@ def type_info(np_type):
     Notes
     -----
     You might be thinking that ``np.finfo`` does this job, and it does, except
-    for PPC long doubles (http://projects.scipy.org/numpy/ticket/2077). This
-    routine protects against errors in ``np.finfo`` by only accepting values
-    that we know are likely to be correct.
+    for PPC long doubles (http://projects.scipy.org/numpy/ticket/2077) and
+    float96 on Windows compiled with Mingw. This routine protects against such
+    errors in ``np.finfo`` by only accepting values that we know are likely to
+    be correct.
     """
     dt = np.dtype(np_type)
     np_type = dt.type
@@ -200,13 +202,18 @@ def type_info(np_type):
     except ValueError:
         pass
     else:
-        return dict(min=np_type(info.min), max=np_type(info.max),
-                    nmant=None, nexp=None, width=width)
+        return dict(min=np_type(info.min), max=np_type(info.max), minexp=None,
+                    maxexp=None, nmant=None, nexp=None, width=width)
     info = np.finfo(dt)
     # Trust the standard IEEE types
     nmant, nexp = info.nmant, info.nexp
-    ret = dict(min=np_type(info.min), max=np_type(info.max), nmant=nmant,
-               nexp=nexp, width=width)
+    ret = dict(min=np_type(info.min),
+               max=np_type(info.max),
+               nmant=nmant,
+               nexp=nexp,
+               minexp=info.minexp,
+               maxexp=info.maxexp,
+               width=width)
     if np_type in (_float16, np.float32, np.float64,
                    np.complex64, np.complex128):
         return ret
@@ -223,14 +230,12 @@ def type_info(np_type):
         pass # these are OK
     elif vals in ((52, 15, 12), # windows float96
                   (52, 15, 16)): # windows float128?
-        # On windows 32 bit at least, float96 appears to be a float64 padded to
-        # 96 bits.  The nexp == 15 is the same as for intel 80 but nexp in fact
-        # appears to be 11 as for float64
-        return dict(min=np_type(info_64.min), max=np_type(info_64.max),
-                    nmant=info_64.nmant, nexp=info_64.nexp, width=width)
+        # On windows 32 bit at least, float96 is Intel 80 storage but operating
+        # at float64 precision. The finfo values give nexp == 15 (as for intel
+        # 80) but in calculations nexp in fact appears to be 11 as for float64
+        return type_info(np.float64).update(dict(width=width))
     elif vals == (1, 1, 16) and processor() == 'powerpc': # broken PPC
-        return dict(min=np_type(info_64.min), max=np_type(info_64.max),
-                    nmant=106, nexp=11, width=width)
+        ret = type_info(np.float64).update(dict(nmant=106, width=width))
     else: # don't recognize the type
         raise FloatingError('We had not expected type %s' % np_type)
     return ret
@@ -438,21 +443,32 @@ def floor_log2(x):
 
     Returns
     -------
-    L : int
-        floor of base 2 log of `x`
+    L : None or int
+        floor of base 2 log of `x`.  None if `x` == 0.
 
     Examples
     --------
     >>> floor_log2(2**9+1)
     9
     >>> floor_log2(-2**9+1)
     8
+    >>> floor_log2(0.5)
+    -1
+    >>> floor_log2(0) is None
+    True
     """
     ip = 0
     rem = abs(x)
-    while rem>=2:
-        ip += 1
-        rem //= 2
+    if rem > 1:
+        while rem>=2:
+            ip += 1
+            rem //= 2
+        return ip
+    elif rem == 0:
+        return None
+    while rem < 1:
+        ip -= 1
+        rem *= 2
     return ip
 
 
@@ -523,3 +539,41 @@ def able_int_type(values):
         if mn >= info.min and mx <= info.max:
             return ityp
     return None
+
+
+def ulp(val=np.float64(1.0)):
+    """ Return gap between `val` and nearest representable number of same type
+
+    This is the value of a unit in the last place (ULP), and is similar in
+    meaning to the MATLAB eps function.
+
+    Parameters
+    ----------
+    val : scalar, optional
+        scalar value of any numpy type.  Default is 1.0 (float64)
+
+    Returns
+    -------
+    ulp_val : scalar
+        gap between `val` and nearest representable number of same type
+
+    Notes
+    -----
+    The wikipedia article on machine epsilon points out that the term *epsilon*
+    can be used in the sense of a unit in the last place (ULP), or as the
+    maximum relative rounding error.  The MATLAB ``eps`` function uses the ULP
+    meaning, but this function is ``ulp`` rather than ``eps`` to avoid confusion
+    between different meanings of *eps*.
+    """
+    val = np.array(val)
+    if not np.isfinite(val):
+        return np.nan
+    if val.dtype.kind in 'iu':
+        return 1
+    aval = np.abs(val)
+    info = type_info(val.dtype)
+    fl2 = floor_log2(aval)
+    if fl2 is None or fl2 < info['minexp']: # subnormal
+        fl2 = info['minexp']
+    # 'nmant' value does not include implicit first bit
+    return 2**(fl2 - info['nmant'])

nibabel/tests/test_casting.py

@@ -4,7 +4,8 @@
 import numpy as np
 
 from ..casting import (float_to_int, shared_range, CastingError, int_to_float,
-                       as_int, int_abs, best_float)
+                       as_int, int_abs, floor_log2, able_int_type, best_float,
+                       ulp)
 
 from numpy.testing import (assert_array_almost_equal, assert_array_equal)
 
@@ -138,6 +139,18 @@ def test_int_abs():
         assert_array_equal(int_abs(in_arr), [e_mn, mx])
 
 
+def test_floor_log2():
+    assert_equal(floor_log2(2**9+1), 9)
+    assert_equal(floor_log2(-2**9+1), 8)
+    assert_equal(floor_log2(2), 1)
+    assert_equal(floor_log2(1), 0)
+    assert_equal(floor_log2(0.5), -1)
+    assert_equal(floor_log2(0.75), -1)
+    assert_equal(floor_log2(0.25), -2)
+    assert_equal(floor_log2(0.24), -3)
+    assert_equal(floor_log2(0), None)
+
+
 def test_able_int_type():
     # The integer type cabable of containing values
     for vals, exp_out in (
@@ -194,18 +207,33 @@ def test_best_float():
         assert_equal(best, np.longdouble)
 
 
-def test_eps():
-    assert_equal(eps(), np.finfo(np.float64).eps)
-    assert_equal(eps(1.0), np.finfo(np.float64).eps)
-    assert_equal(eps(np.float32(1.0)), np.finfo(np.float32).eps)
-    assert_equal(eps(np.float32(1.999)), np.finfo(np.float32).eps)
+def test_ulp():
+    assert_equal(ulp(), np.finfo(np.float64).eps)
+    assert_equal(ulp(1.0), np.finfo(np.float64).eps)
+    assert_equal(ulp(np.float32(1.0)), np.finfo(np.float32).eps)
+    assert_equal(ulp(np.float32(1.999)), np.finfo(np.float32).eps)
     # Integers always return 1
-    assert_equal(eps(1), 1)
-    assert_equal(eps(2**63-1), 1)
+    assert_equal(ulp(1), 1)
+    assert_equal(ulp(2**63-1), 1)
     # negative / positive same
-    assert_equal(eps(-1), 1)
-    assert_equal(eps(7.999), eps(4.0))
-    assert_equal(eps(-7.999), eps(4.0))
-    assert_equal(eps(np.float64(2**54-2)), 2)
-    assert_equal(eps(np.float64(2**54)), 4)
-    assert_equal(eps(np.float64(2**54)), 4)
+    assert_equal(ulp(-1), 1)
+    assert_equal(ulp(7.999), ulp(4.0))
+    assert_equal(ulp(-7.999), ulp(4.0))
+    assert_equal(ulp(np.float64(2**54-2)), 2)
+    assert_equal(ulp(np.float64(2**54)), 4)
+    assert_equal(ulp(np.float64(2**54)), 4)
+    # Infs, NaNs return nan
+    assert_true(np.isnan(ulp(np.inf)))
+    assert_true(np.isnan(ulp(-np.inf)))
+    assert_true(np.isnan(ulp(np.nan)))
+    # 0 gives subnormal smallest
+    subn64 = np.float64(2**(-1022-52))
+    subn32 = np.float32(2**(-126-23))
+    assert_equal(ulp(0.0), subn64)
+    assert_equal(ulp(np.float64(0)), subn64)
+    assert_equal(ulp(np.float32(0)), subn32)
+    # as do multiples of subnormal smallest
+    assert_equal(ulp(subn64 * np.float64(2**52)), subn64)
+    assert_equal(ulp(subn64 * np.float64(2**53)), subn64*2)
+    assert_equal(ulp(subn32 * np.float32(2**23)), subn32)
+    assert_equal(ulp(subn32 * np.float32(2**24)), subn32*2)

nibabel/tests/test_floating.py

@@ -28,7 +28,9 @@ def test_type_info():
     for dtt in np.sctypes['int'] + np.sctypes['uint']:
         info = np.iinfo(dtt)
         infod = type_info(dtt)
-        assert_equal(dict(min=info.min, max=info.max, nexp=None, nmant=None,
+        assert_equal(dict(min=info.min, max=info.max,
+                          nexp=None, nmant=None,
+                          minexp=None, maxexp=None,
                           width=np.dtype(dtt).itemsize), infod)
         assert_equal(infod['min'].dtype.type, dtt)
         assert_equal(infod['max'].dtype.type, dtt)
@@ -37,6 +39,7 @@ def test_type_info():
         infod = type_info(dtt)
         assert_equal(dict(min=info.min, max=info.max,
                           nexp=info.nexp, nmant=info.nmant,
+                          minexp=info.minexp, maxexp=info.maxexp,
                           width=np.dtype(dtt).itemsize),
                      infod)
         assert_equal(infod['min'].dtype.type, dtt)
@@ -54,12 +57,20 @@ def test_type_info():
                 (112, 15, 16), # real float128
                 (106, 11, 16)): # PPC head, tail doubles, expected values
         assert_equal(dict(min=info.min, max=info.max,
+                          minexp=info.minexp, maxexp=info.maxexp,
                           nexp=info.nexp, nmant=info.nmant, width=width),
                      infod)
     elif vals == (1, 1, 16): # bust info for PPC head / tail longdoubles
         assert_equal(dict(min=dbl_info.min, max=dbl_info.max,
+                          minexp=-1022, maxexp=1024,
                           nexp=11, nmant=106, width=16),
                      infod)
+    elif vals == (52, 15, 12):
+        exp_res = type_info(np.float64)
+        exp_res['width'] = width
+        assert_equal(exp_res, infod)
+    else:
+        raise ValueError("Unexpected float type to test")
 
 
 def test_nmant():