Merge pull request #91 from matthew-brett/array-writers-ppc

Array writers ppc

After merging the array writers branch recently, I discovered that PPC min and max values for longdouble were wrong. This branch introduces a type_info function to return the min, max etc for a particular numpy type, and special-cases the broken PPC longdouble.

I've used the type_info function throughout for consistency (instead of np.finfo).

Author: matthew-brett

nibabel/arraywriters.py

@@ -25,7 +25,7 @@ def __init__(self, array, out_dtype=None)
 
 import numpy as np
 
-from .casting import int_to_float, as_int, int_abs
+from .casting import int_to_float, as_int, int_abs, type_info
 from .volumeutils import finite_range, array_to_file
 
 
@@ -329,13 +329,13 @@ def _range_scale(self):
         """ Calculate scaling based on data range and output type """
         mn, mx = self.finite_range() # These can be floats or integers
         out_dtype = self._out_dtype
+        info = type_info(out_dtype)
+        t_mn_mx = info['min'], info['max']
         if out_dtype.kind == 'f':
-            t_mn_mx = np.finfo(out_dtype).min, np.finfo(out_dtype).max
             # But we want maximum precision for the calculations. Casting will
             # not lose precision because min/max are of fp type.
             t_min, t_max = np.array(t_mn_mx, dtype = np.longdouble)
         else: # (u)int
-            t_mn_mx = np.iinfo(out_dtype).min, np.iinfo(out_dtype).max
             t_min, t_max = [int_to_float(v, np.longdouble) for v in t_mn_mx]
         if self._out_dtype.kind == 'u':
             if mn < 0 and mx > 0:
@@ -484,7 +484,8 @@ def _range_scale(self):
             mn2mx = int_to_float(as_int(mx) - as_int(mn), np.longdouble)
         if out_dtype.kind == 'f':
             # Type range, these are also floats
-            t_mn_mx = np.finfo(out_dtype).min, np.finfo(out_dtype).max
+            info = type_info(out_dtype)
+            t_mn_mx = info['min'], info['max']
         else:
             t_mn_mx = np.iinfo(out_dtype).min, np.iinfo(out_dtype).max
             t_mn_mx= [int_to_float(v, np.longdouble) for v in t_mn_mx]

nibabel/casting.py

@@ -1,6 +1,8 @@
 """ Utilties for casting floats to integers
 """
 
+from platform import processor
+
 import numpy as np
 
 
@@ -152,56 +154,75 @@ def shared_range(flt_type, int_type):
 except AttributeError: # float16 not present in np < 1.6
     _float16 = None
 
-# The number of significand digits in IEEE floating point formats, not including
-# the implicit leading 0.  See http://en.wikipedia.org/wiki/IEEE_754-2008
-_flt_nmant = {
-    _float16: 10,
-    np.float32: 23,
-    np.float64: 52,
-    }
-
 
 class FloatingError(Exception):
     pass
 
 
-def flt2nmant(flt_type):
-    """ Number of significand bits in float type `flt_type`
+def type_info(np_type):
+    """ Return dict with min, max, nexp, nmant, width for numpy type `np_type`
+
+    Type can be integer in which case nexp and nmant are None.
 
     Parameters
     ----------
-    flt_type : object
-        Numpy floating point type, such as np.float32
+    np_type : numpy type specifier
+        Any specifier for a numpy dtype
 
     Returns
     -------
-    nmant : int
-        Number of digits in the signficand
+    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`.
+
+    Raises
+    ------
+    FloatingError : for floating point types we don't recognize
+
+    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.
     """
-    try:
-        return _flt_nmant[flt_type]
-    except KeyError:
+    dt = np.dtype(np_type)
+    np_type = dt.type
+    width = dt.itemsize
+    try: # integer type
+        info = np.iinfo(dt)
+    except ValueError:
         pass
-    fi = np.finfo(flt_type)
-    nmant, nexp = fi.nmant, fi.nexp
-    # Assuming the np.float type is always IEEE 64 bit
-    if flt_type is np.float and (nmant, nexp) == (52, 11):
-        return 52
-    # Now we should be testing long doubles
-    assert flt_type is np.longdouble
-    if (nmant, nexp) == (63, 15): # 80-bit intel type
-        return 63 # Not including explicit first digit
-    # We test the declared nmant by stepping up and down.  These tests assume a
-    # binary format
-    i_end_contig = 2**(nmant+1) # int
-    f_end_contig = flt_type(i_end_contig)
-    # We need as_int here because long doubles do not necessarily convert
-    # correctly to ints with int() - see
-    # http://projects.scipy.org/numpy/ticket/1395
-    if as_int(f_end_contig-1) == (i_end_contig-1): # still representable
-        if as_int(f_end_contig+1) == i_end_contig: # Rounding down
-            return nmant
-    raise FloatingError('Cannot be confident of nmant value for %s' % flt_type)
+    else:
+        return dict(min=np_type(info.min), max=np_type(info.max),
+                    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)
+    if np_type in (_float16, np.float32, np.float64,
+                   np.complex64, np.complex128):
+        return ret
+    if dt.kind == 'c':
+        assert np_type is np.longcomplex
+        vals = (nmant, nexp, width / 2)
+    else:
+        assert np_type is np.longdouble
+        vals = (nmant, nexp, width)
+    if vals in ((112, 15, 16), # binary128
+                (63, 15, 12), (63, 15, 16)): # Intel extended 80
+        pass # these are OK
+    elif vals == (1, 1, 16) and processor() == 'powerpc': # broken PPC
+        dbl_info = np.finfo(np.float64)
+        return dict(min=np_type(dbl_info.min), max=np_type(dbl_info.max),
+                    nmant=106, nexp=11, width=width)
+    else: # don't recognize the type
+        raise FloatingError('We had not expected this type')
+    return ret
 
 
 def as_int(x, check=True):
@@ -342,14 +363,14 @@ def floor_exact(val, flt_type):
         faval = int_to_float(aval, flt_type)
     except OverflowError:
         faval = np.inf
+    info = type_info(flt_type)
     if faval == np.inf:
-        return sign * np.finfo(flt_type).max
+        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
     # Float casting made the value go up
-    nmant = flt2nmant(flt_type)
-    biggest_gap = 2**(floor_log2(aval) - nmant)
+    biggest_gap = 2**(floor_log2(aval) - info['nmant'])
     assert biggest_gap > 1
     faval -= flt_type(biggest_gap)
     return sign * faval

nibabel/tests/test_arraywriters.py

@@ -37,7 +37,7 @@ def __init__(self, array, out_dtype=None, order='F')
                             WriterError, ScalingError, ArrayWriter,
                             make_array_writer, get_slope_inter)
 
-from ..casting import int_abs
+from ..casting import int_abs, type_info
 
 from ..volumeutils import array_from_file, apply_read_scaling
 
@@ -354,8 +354,8 @@ def test_io_scaling():
                                    (np.float32, np.int16, None)):
         out_dtype = np.dtype(out_type)
         arr = np.zeros((3,), dtype=in_type)
-        info = np.finfo(in_type) if arr.dtype.kind == 'f' else np.iinfo(in_type)
-        arr[0], arr[1] = info.min, info.max
+        info = type_info(in_type)
+        arr[0], arr[1] = info['min'], info['max']
         aw = SlopeInterArrayWriter(arr, out_dtype, calc_scale=False)
         if not err is None:
             assert_raises(err, aw.calc_scale)
@@ -429,17 +429,16 @@ def test_writers_roundtrip():
 
 
 def test_to_float():
+    start, stop = 0, 100
     for in_type in NUMERIC_TYPES:
-        if in_type in IUINT_TYPES:
-            info = np.iinfo(in_type)
-            mn, mx, start, stop, step = info.min, info.max, 0, 100, 1
-        else:
-            info = np.finfo(in_type)
-            mn, mx, start, stop, step = info.min, info.max, 0, 100, 0.5
+        step = 1 if in_type in IUINT_TYPES else 0.5
+        info = type_info(in_type)
+        mn, mx = info['min'], info['max']
         arr = np.arange(start, stop, step, dtype=in_type)
         arr[0] = mn
         arr[-1] = mx
         for out_type in CFLOAT_TYPES:
+            out_info = type_info(out_type)
             for klass in (SlopeInterArrayWriter, SlopeArrayWriter,
                           ArrayWriter):
                 if in_type in COMPLEX_TYPES and out_type in FLOAT_TYPES:
@@ -451,8 +450,7 @@ def test_to_float():
                 arr_back = round_trip(aw)
                 assert_array_equal(arr.astype(out_type), arr_back)
                 # Check too-big values overflowed correctly
-                out_min = np.finfo(out_type).min
-                out_max = np.finfo(out_type).max
+                out_min, out_max = out_info['min'], out_info['max']
                 assert_true(np.all(arr_back[arr > out_max] == np.inf))
                 assert_true(np.all(arr_back[arr < out_min] == -np.inf))
 
@@ -498,8 +496,8 @@ def test_writer_maker():
 def test_float_int_min_max():
     # Conversion between float and int
     for in_dt in FLOAT_TYPES:
-        finf = np.finfo(in_dt)
-        arr = np.array([finf.min, finf.max], dtype=in_dt)
+        finf = type_info(in_dt)
+        arr = np.array([finf['min'], finf['max']], dtype=in_dt)
         for out_dt in IUINT_TYPES:
             try:
                 aw = SlopeInterArrayWriter(arr, out_dt)