RF: use best_float instead of np.longdouble

Refactor to use float64 instead of np.longdouble on systems (?only
windows) in which np.longdouble is no more capable than float64.

Author: matthew-brett

nibabel/arraywriters.py

@@ -25,7 +25,8 @@ def __init__(self, array, out_dtype=None)
 
 import numpy as np
 
-from .casting import int_to_float, as_int, int_abs, type_info
+from .casting import (int_to_float, as_int, int_abs, type_info, floor_exact,
+                      best_float)
 from .volumeutils import finite_range, array_to_file
 
 
@@ -331,12 +332,13 @@ def _range_scale(self):
         out_dtype = self._out_dtype
         info = type_info(out_dtype)
         t_mn_mx = info['min'], info['max']
+        big_float = best_float()
         if out_dtype.kind == 'f':
             # 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)
+            t_min, t_max = np.array(t_mn_mx, dtype = big_float)
         else: # (u)int
-            t_min, t_max = [int_to_float(v, np.longdouble) for v in t_mn_mx]
+            t_min, t_max = [int_to_float(v, big_float) for v in t_mn_mx]
         if self._out_dtype.kind == 'u':
             if mn < 0 and mx > 0:
                 raise WriterError('Cannot scale negative and positive '
@@ -472,27 +474,28 @@ def _range_scale(self):
             self.inter = mn
             return
         # Straight mx-mn can overflow.
+        big_float = best_float() # usually longdouble except in win 32
         if mn.dtype.kind == 'f': # Already floats
             # float64 and below cast correctly to longdouble.  Longdouble needs
             # no casting
-            mn2mx = np.diff(np.array([mn, mx], dtype=np.longdouble))
+            mn2mx = np.diff(np.array([mn, mx], dtype=big_float))
         else: # max possible (u)int range is 2**64-1 (int64, uint64)
             # int_to_float covers this range.  On windows longdouble is the same
             # as double so mn2mx will be 2**64 - thus overestimating slope
             # slightly.  Casting to int needed to allow mx-mn to be larger than
             # the largest (u)int value
-            mn2mx = int_to_float(as_int(mx) - as_int(mn), np.longdouble)
+            mn2mx = int_to_float(as_int(mx) - as_int(mn), big_float)
         if out_dtype.kind == 'f':
             # Type range, these are also floats
             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]
+            t_mn_mx= [int_to_float(v, big_float) for v in t_mn_mx]
         # We want maximum precision for the calculations. Casting will
         # not lose precision because min/max are of fp type.
         assert [v.dtype.kind for v in t_mn_mx] == ['f', 'f']
-        scaled_mn2mx = np.diff(np.array(t_mn_mx, dtype = np.longdouble))
+        scaled_mn2mx = np.diff(np.array(t_mn_mx, dtype = big_float))
         slope = mn2mx / scaled_mn2mx
         self.inter = mn - t_mn_mx[0] * slope
         self.slope = slope

nibabel/casting.py

@@ -469,3 +469,17 @@ def best_float():
     if type_info(np.longdouble)['nmant'] > type_info(np.float64)['nmant']:
         return np.longdouble
     return np.float64
+
+
+def ok_floats():
+    """ Return floating point types sorted by precision
+
+    Remove longdouble if it has no higher precision than float64
+    """
+    floats = sorted(np.sctypes['float'], key=lambda f : type_info(f)['nmant'])
+    if best_float() != np.longdouble:
+        floats.remove(np.longdouble)
+    return floats
+
+
+OK_FLOATS = ok_floats()

nibabel/volumeutils.py

@@ -16,7 +16,7 @@
 import numpy as np
 
 from .py3k import isfileobj, ZEROB
-from .casting import shared_range, type_info
+from .casting import (shared_range, type_info, as_int, best_float, OK_FLOATS)
 
 sys_is_le = sys.byteorder == 'little'
 native_code = sys_is_le and '<' or '>'
@@ -1057,7 +1057,7 @@ def best_write_scale_ftype(arr, slope = 1.0, inter = 0.0, default=np.float32):
     try:
         return _ftype4scaled_finite(arr, slope, inter, 'write', default)
     except ValueError:
-        return FLOAT_TYPES[-1]
+        return OK_FLOATS[-1]
 
 
 def better_float_of(first, second, default=np.float32):
@@ -1111,14 +1111,17 @@ def _ftype4scaled_finite(tst_arr, slope, inter, direction='read',
     """ Smallest float type for scaling of `tst_arr` that does not overflow
     """
     assert direction in ('read', 'write')
-    def_ind = FLOAT_TYPES.index(default)
+    if not default in OK_FLOATS and default is np.longdouble:
+        # Omitted longdouble
+        return default
+    def_ind = OK_FLOATS.index(default)
     # promote to arrays to avoid numpy scalar casting rules
     tst_arr = np.atleast_1d(tst_arr)
     slope = np.atleast_1d(slope)
     inter = np.atleast_1d(inter)
     warnings.filterwarnings('ignore', '.*overflow.*', RuntimeWarning)
     try:
-        for ftype in FLOAT_TYPES[def_ind:]:
+        for ftype in OK_FLOATS[def_ind:]:
             tst_trans = tst_arr.copy()
             slope = slope.astype(ftype)
             inter = inter.astype(ftype)