NF - add able_int_type function

Function selects integer type capable of representing all passed values

Numpy doesn't itself try and minimize the integer type so able_int_type
has to do this the slow way.

Confirm able_int_type corresponds to numpy casting.

Author: matthew-brett

nibabel/casting.py

@@ -1,4 +1,7 @@
-""" Utilties for casting floats to integers
+""" Utilties for casting numpy values in various ways
+
+Most routines work round some numpy oddities in floating point precision and
+casting.  Others work round numpy casting to and from python ints
 """
 
 from platform import processor
@@ -483,3 +486,40 @@ def ok_floats():
 
 
 OK_FLOATS = ok_floats()
+
+
+def able_int_type(values):
+    """ Find the smallest integer numpy type to contain sequence `values`
+
+    Prefers uint to int if minimum is >= 0
+
+    Parameters
+    ----------
+    values : sequence
+        sequence of integer values
+
+    Returns
+    -------
+    itype : None or numpy type
+        numpy integer type or None if no integer type holds all `values`
+
+    Examples
+    --------
+    >>> able_int_type([0, 1]) == np.uint8
+    True
+    >>> able_int_type([-1, 1]) == np.int8
+    True
+    """
+    if any([v % 1 for v in values]):
+        return None
+    mn = min(values)
+    mx = max(values)
+    if mn >= 0:
+        for ityp in np.sctypes['uint']:
+            if mx <= np.iinfo(ityp).max:
+                return ityp
+    for ityp in np.sctypes['int']:
+        info = np.iinfo(ityp)
+        if mn >= info.min and mx <= info.max:
+            return ityp
+    return None

nibabel/tests/test_casting.py

@@ -136,3 +136,76 @@ def test_int_abs():
         assert_equal(int_abs(mx), mx)
         assert_equal(int_abs(mn), e_mn)
         assert_array_equal(int_abs(in_arr), [e_mn, mx])
+
+
+def test_able_int_type():
+    # The integer type cabable of containing values
+    for vals, exp_out in (
+        ([0, 1], np.uint8),
+        ([0, 255], np.uint8),
+        ([-1, 1], np.int8),
+        ([0, 256], np.uint16),
+        ([-1, 128], np.int16),
+        ([0.1, 1], None),
+        ([0, 2**16], np.uint32),
+        ([-1, 2**15], np.int32),
+        ([0, 2**32], np.uint64),
+        ([-1, 2**31], np.int64),
+        ([-1, 2**64-1], None),
+        ([0, 2**64-1], np.uint64),
+        ([0, 2**64], None)):
+        assert_equal(able_int_type(vals), exp_out)
+
+
+def test_able_casting():
+    # Check the able_int_type function guesses numpy out type
+    types = np.sctypes['int'] + np.sctypes['uint']
+    for in_type in types:
+        in_info = np.iinfo(in_type)
+        in_mn, in_mx = in_info.min, in_info.max
+        A = np.zeros((1,), dtype=in_type)
+        for out_type in types:
+            out_info = np.iinfo(out_type)
+            out_mn, out_mx = out_info.min, out_info.max
+            B = np.zeros((1,), dtype=out_type)
+            ApBt = (A + B).dtype.type
+            able_type = able_int_type([in_mn, in_mx, out_mn, out_mx])
+            if able_type is None:
+                assert_equal(ApBt, np.float64)
+                continue
+            # Use str for comparison to avoid int32/64 vs intp comparison
+            # failures
+            assert_equal(np.dtype(ApBt).str, np.dtype(able_type).str)
+
+
+def test_best_float():
+    # Finds the most capable floating point type
+    # The only time this isn't np.longdouble is when np.longdouble has float64
+    # precision.
+    best = best_float()
+    end_of_ints = np.float64(2**53)
+    # float64 has continuous integers up to 2**53
+    assert_equal(end_of_ints, end_of_ints + 1)
+    # longdouble may have more, but not on 32 bit windows, at least
+    end_of_ints = np.longdouble(2**53)
+    if end_of_ints == (end_of_ints + 1): # off continuous integers
+        assert_equal(best, np.float64)
+    else:
+        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)
+    # Integers always return 1
+    assert_equal(eps(1), 1)
+    assert_equal(eps(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)