RF: make floor_exact like np.floor, add ceil_exact

Floor exact previously floored towards 0. This was unlike np.floor,
which finds the nearest integer <= val.

Added ceil_exact function (nearest exact integer >= val).

Infs from floor / ceil_exact return infs (previously pulled infs down to
float max)

Consider -+ infs in shared_range calculation (now infs not removed by
floor / ceil_exact).

Author: matthew-brett

nibabel/casting.py

@@ -129,10 +129,10 @@ def shared_range(flt_type, int_type):
 
     Examples
     --------
-    >>> shared_range(np.float32, np.int32)
-    (-2147483648.0, 2147483520.0)
-    >>> shared_range('f4', 'i4')
-    (-2147483648.0, 2147483520.0)
+    >>> shared_range(np.float32, np.int32) == (-2147483648.0, 2147483520.0)
+    True
+    >>> shared_range('f4', 'i4') == (-2147483648.0, 2147483520.0)
+    True
     """
     flt_type = np.dtype(flt_type).type
     int_type = np.dtype(int_type).type
@@ -143,9 +143,15 @@ def shared_range(flt_type, int_type):
     except KeyError:
         pass
     ii = np.iinfo(int_type)
-    mn_mx = floor_exact(ii.min, flt_type), floor_exact(ii.max, flt_type)
-    _SHARED_RANGES[key] = mn_mx
-    return mn_mx
+    fi = np.finfo(flt_type)
+    mn = ceil_exact(ii.min, flt_type)
+    if mn == -np.inf:
+        mn = fi.min
+    mx = floor_exact(ii.max, flt_type)
+    if mx == np.inf:
+        mx = fi.max
+    _SHARED_RANGES[key] = (mn, mx)
+    return mn, mx
 
 # ----------------------------------------------------------------------------
 # Routines to work out the next lowest representable integer in floating point
@@ -335,7 +341,7 @@ def int_to_float(val, flt_type):
 
 
 def floor_exact(val, flt_type):
-    """ Get nearest exact integer to `val`, towards 0, in float type `flt_type`
+    """ Return nearest exact integer <= `val` in float type `flt_type`
 
     Parameters
     ----------
@@ -344,15 +350,14 @@ def floor_exact(val, flt_type):
         because large integers cast as floating point may be rounded by the
         casting process.
     flt_type : numpy type
-        numpy float type.  Only IEEE types supported (np.float16, np.float32,
-        np.float64)
+        numpy float type.
 
     Returns
     -------
     floor_val : object
-        value of same floating point type as `val`, that is the next excat
-        integer in this type, towards zero, or == `val` if val is exactly
-        representable.
+        value of same floating point type as `val`, that is the nearest exact
+        integer in this type such that `floor_val` <= `val`.  Thus if `val` is
+        exact in `flt_type`, `floor_val` == `val`.
 
     Examples
     --------
@@ -370,26 +375,74 @@ def floor_exact(val, flt_type):
 
     >>> floor_exact(2**24+1, np.float32) == 2**24
     True
+
+    As for the numpy floor function, negatives floor towards -inf
+
+    >>> floor_exact(-2**24-1, np.float32) == -2**24-2
+    True
     """
     val = int(val)
     flt_type = np.dtype(flt_type).type
-    sign = val > 0 and 1 or -1
-    aval = abs(val)
+    sign = 1 if val > 0 else -1
     try: # int_to_float deals with longdouble safely
-        faval = int_to_float(aval, flt_type)
+        fval = int_to_float(val, flt_type)
     except OverflowError:
-        faval = np.inf
+        return sign * np.inf
+    if not np.isfinite(fval):
+        return fval
     info = type_info(flt_type)
-    if faval == np.inf:
-        return sign * info['max']
-    if as_int(faval) <= aval: # as_int deals with longdouble safely
-        # Float casting has made the value go down or stay the same
-        return sign * faval
+    diff = val - as_int(fval)
+    if diff >= 0: # floating point value <= val
+        return fval
     # Float casting made the value go up
-    biggest_gap = 2**(floor_log2(aval) - info['nmant'])
+    biggest_gap = 2**(floor_log2(val) - info['nmant'])
     assert biggest_gap > 1
-    faval -= flt_type(biggest_gap)
-    return sign * faval
+    fval -= flt_type(biggest_gap)
+    return fval
+
+
+def ceil_exact(val, flt_type):
+    """ Return nearest exact integer >= `val` in float type `flt_type`
+
+    Parameters
+    ----------
+    val : int
+        We have to pass val as an int rather than the floating point type
+        because large integers cast as floating point may be rounded by the
+        casting process.
+    flt_type : numpy type
+        numpy float type.
+
+    Returns
+    -------
+    ceil_val : object
+        value of same floating point type as `val`, that is the nearest exact
+        integer in this type such that `floor_val` >= `val`.  Thus if `val` is
+        exact in `flt_type`, `ceil_val` == `val`.
+
+    Examples
+    --------
+    Obviously 2 is within the range of representable integers for float32
+
+    >>> ceil_exact(2, np.float32)
+    2.0
+
+    As is 2**24-1 (the number of significand digits is 23 + 1 implicit)
+
+    >>> ceil_exact(2**24-1, np.float32) == 2**24-1
+    True
+
+    But 2**24+1 gives a number that float32 can't represent exactly
+
+    >>> ceil_exact(2**24+1, np.float32) == 2**24+2
+    True
+
+    As for the numpy ceil function, negatives ceil towards inf
+
+    >>> ceil_exact(-2**24-1, np.float32) == -2**24
+    True
+    """
+    return -floor_exact(-val, flt_type)
 
 
 def int_abs(arr):

nibabel/tests/test_floating.py

@@ -4,8 +4,8 @@
 
 import numpy as np
 
-from ..casting import (floor_exact, as_int, FloatingError, int_to_float,
-                       floor_log2, type_info)
+from ..casting import (floor_exact, ceil_exact, as_int, FloatingError,
+                       int_to_float, floor_log2, type_info)
 
 from nose import SkipTest
 from nose.tools import assert_equal, assert_raises
@@ -164,7 +164,8 @@ def test_floor_exact_16():
     # A normal integer can generate an inf in float16
     if not have_float16:
         raise SkipTest('No float16')
-    assert_equal(floor_exact(2**31, np.float16), type_info(np.float16)['max'])
+    assert_equal(floor_exact(2**31, np.float16), np.inf)
+    assert_equal(floor_exact(-2**31, np.float16), -np.inf)
 
 
 def test_floor_exact_64():
@@ -192,37 +193,50 @@ def test_floor_exact():
     # When numbers go above int64 - I believe, numpy comparisons break down,
     # so we have to cast to int before comparison
     int_flex = lambda x, t : as_int(floor_exact(x, t))
+    int_ceex = lambda x, t : as_int(ceil_exact(x, t))
     for t in to_test:
         # A number bigger than the range returns the max
         info = type_info(t)
-        assert_equal(floor_exact(2**5000, t), info['max'])
+        assert_equal(floor_exact(2**5000, t), np.inf)
+        assert_equal(ceil_exact(2**5000, t), np.inf)
+        # A number more negative returns -inf
+        assert_equal(floor_exact(-2**5000, t), -np.inf)
+        assert_equal(ceil_exact(-2**5000, t), -np.inf)
+        # Check around end of integer precision
         nmant =  info['nmant']
         for i in range(nmant+1):
             iv = 2**i
             # up to 2**nmant should be exactly representable
-            assert_equal(int_flex(iv, t), iv)
-            assert_equal(int_flex(-iv, t), -iv)
-            assert_equal(int_flex(iv-1, t), iv-1)
-            assert_equal(int_flex(-iv+1, t), -iv+1)
+            for func in (int_flex, int_ceex):
+                assert_equal(func(iv, t), iv)
+                assert_equal(func(-iv, t), -iv)
+                assert_equal(func(iv-1, t), iv-1)
+                assert_equal(func(-iv+1, t), -iv+1)
         # The nmant value for longdouble on PPC appears to be conservative, so
         # that the tests for behavior above the nmant range fail
         if t is np.longdouble and processor() == 'powerpc':
             continue
-        # 2**(nmant+1) can't be exactly represented
+        # Confirm to ourselves that 2**(nmant+1) can't be exactly represented
         iv = 2**(nmant+1)
         assert_equal(int_flex(iv+1, t), iv)
+        assert_equal(int_ceex(iv+1, t), iv+2)
         # negatives
-        assert_equal(int_flex(-iv-1, t), -iv)
+        assert_equal(int_flex(-iv-1, t), -iv-2)
+        assert_equal(int_ceex(-iv-1, t), -iv)
         # The gap in representable numbers is 2 above 2**(nmant+1), 4 above
         # 2**(nmant+2), and so on.
         for i in range(5):
             iv = 2**(nmant+1+i)
             gap = 2**(i+1)
             assert_equal(as_int(t(iv) + t(gap)), iv+gap)
-            for j in range(1,gap):
+            for j in range(1, gap):
                 assert_equal(int_flex(iv+j, t), iv)
                 assert_equal(int_flex(iv+gap+j, t), iv+gap)
+                assert_equal(int_ceex(iv+j, t), iv+gap)
+                assert_equal(int_ceex(iv+gap+j, t), iv+2*gap)
             # negatives
-            for j in range(1,gap):
-                assert_equal(int_flex(-iv-j, t), -iv)
-                assert_equal(int_flex(-iv-gap-j, t), -iv-gap)
+            for j in range(1, gap):
+                assert_equal(int_flex(-iv-j, t), -iv-gap)
+                assert_equal(int_flex(-iv-gap-j, t), -iv-2*gap)
+                assert_equal(int_ceex(-iv-j, t), -iv)
+                assert_equal(int_ceex(-iv-gap-j, t), -iv-gap)