BF: work round horrible MSVC uint64 casting

MSVC casting floats to uint64 casts floats > 2**63 to 2*63, at least on
32-bit Windows.  Platform standard compiler.

Author: matthew-brett

nibabel/arraywriters.py

@@ -411,15 +411,16 @@ def _do_scaling(self):
     def _iu2iu(self):
         # (u)int to (u)int scaling
         mn, mx = self.finite_range()
-        if self._out_dtype.kind == 'u':
+        out_dt = self._out_dtype
+        if out_dt.kind == 'u':
             # We're checking for a sign flip.  This can only work for uint
             # output, because, for int output, the abs min of the type is
             # greater than the abs max, so the data either fits into the range
             # (tested for in _do_scaling), or this test can't pass. Need abs
             # that deals with max neg ints. abs problem only arises when all
             # the data is set to max neg integer value
-            imax = np.iinfo(self._out_dtype).max
-            if mx <= 0 and int_abs(mn) <= imax: # sign flip enough?
+            o_min, o_max = shared_range(self.scaler_dtype, out_dt)
+            if mx <= 0 and int_abs(mn) <= as_int(o_max): # sign flip enough?
                 # -1.0 * arr will be in scaler_dtype precision
                 self.slope = -1.0
                 return
@@ -563,15 +564,19 @@ def _iu2iu(self):
         # range may be greater than the largest integer for this type.
         # as_int needed to work round numpy 1.4.1 int casting bug
         out_dtype = self._out_dtype
-        t_min, t_max = np.iinfo(out_dtype).min, np.iinfo(out_dtype).max
-        type_range = as_int(t_max) - as_int(t_min)
+        # Options in this method are scaling using intercept only.  These will
+        # have to pass through ``self.scaler_dtype`` (because the intercept is
+        # in this type).
+        o_min, o_max = [as_int(v)
+                        for v in shared_range(self.scaler_dtype, out_dtype)]
+        type_range = o_max - o_min
         mn2mx = mx - mn
         if mn2mx <= type_range: # might offset be enough?
-            if t_min == 0: # uint output - take min to 0
+            if o_min == 0:  # uint output - take min to 0
                 # decrease offset with floor_exact, meaning mn >= t_min after
                 # subtraction.  But we may have pushed the data over t_max,
                 # which we check below
-                inter = floor_exact(mn - t_min, self.scaler_dtype)
+                inter = floor_exact(mn - o_min, self.scaler_dtype)
             else: # int output - take midpoint to 0
                 # ceil below increases inter, pushing scale up to 0.5 towards
                 # -inf, because ints have abs min == abs max + 1
@@ -581,8 +586,8 @@ def _iu2iu(self):
                 inter = floor_exact(midpoint, self.scaler_dtype)
             # Need to check still in range after floor_exact-ing
             int_inter = as_int(inter)
-            assert mn - int_inter >= t_min
-            if mx - int_inter <= t_max:
+            assert mn - int_inter >= o_min
+            if mx - int_inter <= o_max:
                 self.inter = inter
                 return
         # Try slope options (sign flip) and then range scaling

nibabel/casting.py

@@ -14,6 +14,11 @@ class CastingError(Exception):
     pass
 
 
+#: Test for VC truncation when casting floats to uint64
+_test_val = 2**63 + 2**11 # Should be exactly representable in float64
+TRUNC_UINT64 = np.float64(_test_val).astype(np.uint64) != _test_val
+
+
 def float_to_int(arr, int_type, nan2zero=True, infmax=False):
     """ Convert floating point array `arr` to type `int_type`
 
@@ -151,6 +156,8 @@ def shared_range(flt_type, int_type):
     mx = floor_exact(ii.max, flt_type)
     if mx == np.inf:
         mx = fi.max
+    elif TRUNC_UINT64 and int_type == np.uint64:
+        mx = min(mx, flt_type(2**63))
     _SHARED_RANGES[key] = (mn, mx)
     return mn, mx
 

nibabel/tests/test_arraywriters.py

@@ -135,17 +135,29 @@ def test_no_scaling():
         kwargs = (dict(check_scaling=False) if awt == ArrayWriter
                   else dict(calc_scale=False))
         aw = awt(arr, out_dtype, **kwargs)
-        with suppress_warnings():  # cast to real from cplx
+        with suppress_warnings():
             back_arr = round_trip(aw)
-            exp_back = arr.astype(float)
+        exp_back = arr.copy()
         if out_dtype in IUINT_TYPES:
-            with np.errstate(invalid='ignore'):
-                exp_back = np.round(exp_back)
-            if hasattr(aw, 'slope') and in_dtype in FLOAT_TYPES:
-                # Finite scaling sets infs to min / max
-                exp_back = np.clip(exp_back, 0, 1)
-            else:
-                exp_back = np.clip(exp_back, *shared_range(float, out_dtype))
+            if in_dtype in CFLOAT_TYPES:
+                with suppress_warnings():
+                    exp_back = exp_back.astype(float)
+                with np.errstate(invalid='ignore'):
+                    exp_back = np.round(exp_back)
+                if hasattr(aw, 'slope') and in_dtype in FLOAT_TYPES:
+                    # Finite scaling sets infs to min / max
+                    exp_back = np.clip(exp_back, 0, 1)
+                else:
+                    exp_back = np.clip(exp_back, *shared_range(float, out_dtype))
+            else: # iu input type
+                mn_out, mx_out = _dt_min_max(out_dtype)
+                if (mn_in, mx_in) != (mn_out, mx_out):
+                    exp_back = np.clip(exp_back,
+                                       max(mn_in, mn_out),
+                                       min(mx_in, mx_out))
+        elif in_dtype in COMPLEX_TYPES:  # always cast to real from cplx
+            with suppress_warnings():
+                exp_back = exp_back.astype(float)
         exp_back = exp_back.astype(out_dtype)
         # Sometimes working precision is float32 - allow for small differences
         assert_allclose_safely(back_arr, exp_back)

nibabel/tests/test_spm99analyze.py

@@ -343,15 +343,29 @@ def test_no_scaling(self):
             with suppress_warnings():  # invalid mult
                 back_arr = rt_img.get_data()
             exp_back = arr.copy()
-            if in_dtype not in COMPLEX_TYPES:
-                exp_back = arr.astype(float)
             if out_dtype in IUINT_TYPES:
-                with np.errstate(invalid='ignore'):
-                    exp_back = np.round(exp_back)
-                exp_back = np.clip(exp_back, *shared_range(float, out_dtype))
-                exp_back = exp_back.astype(out_dtype).astype(float)
-            else:
+                if in_dtype in FLOAT_TYPES:
+                    # Working precision is (at least) float
+                    exp_back = exp_back.astype(float)
+                    with np.errstate(invalid='ignore'):
+                        exp_back = np.round(exp_back)
+                    if in_dtype in FLOAT_TYPES:
+                        # Clip to shared range of working precision
+                        exp_back = np.clip(exp_back,
+                                           *shared_range(float, out_dtype))
+                else:  # iu input type
+                    # Clipping only to integer range
+                    mn_out, mx_out = _dt_min_max(out_dtype)
+                    if (mn_in, mx_in) != (mn_out, mx_out):
+                        exp_back = np.clip(exp_back,
+                                           max(mn_in, mn_out),
+                                           min(mx_in, mx_out))
+                # exp_back = exp_back.astype(out_dtype)
+            if out_dtype in COMPLEX_TYPES:
                 exp_back = exp_back.astype(out_dtype)
+            else:
+                # Cast to working precision
+                exp_back = exp_back.astype(float)
             # Allow for small differences in large numbers
             with suppress_warnings():  # invalid value
                 assert_allclose_safely(back_arr,

nibabel/tests/test_utils.py

@@ -536,6 +536,9 @@ def test_a2f_scaled_unscaled():
         mn_in, mx_in = _dt_min_max(in_dtype)
         nan_val = np.nan if in_dtype in CFLOAT_TYPES else 10
         arr = np.array([mn_in, -1, 0, 1, mx_in, nan_val], dtype=in_dtype)
+        iu_via_float = (in_dtype in CFLOAT_TYPES or
+                        (intercept, divslope) != (0, 1) or
+                        out_dtype in CFLOAT_TYPES)
         mn_out, mx_out = _dt_min_max(out_dtype)
         nan_fill = -intercept / divslope
         if out_dtype in IUINT_TYPES:
@@ -555,20 +558,24 @@ def test_a2f_scaled_unscaled():
                                     divslope=divslope,
                                     intercept=intercept)
             exp_back = arr.copy()
-            if out_dtype in IUINT_TYPES:
-                exp_back[np.isnan(exp_back)] = 0
-            if in_dtype not in COMPLEX_TYPES:
-                exp_back = exp_back.astype(float)
-            if intercept != 0:
-                exp_back -= intercept
-            if divslope != 1:
-                exp_back /= divslope
-            if out_dtype in IUINT_TYPES:
-                exp_back = np.round(exp_back).astype(float)
-                exp_back = np.clip(exp_back, *shared_range(float, out_dtype))
-                exp_back = exp_back.astype(out_dtype)
-            else:
-                exp_back = exp_back.astype(out_dtype)
+            if iu_via_float:
+                if out_dtype in IUINT_TYPES:
+                    exp_back[np.isnan(exp_back)] = 0
+                if in_dtype not in COMPLEX_TYPES:
+                    exp_back = exp_back.astype(float)
+                if intercept != 0:
+                    exp_back -= intercept
+                if divslope != 1:
+                    exp_back /= divslope
+                if out_dtype in IUINT_TYPES:
+                    exp_back = np.round(exp_back).astype(float)
+                    exp_back = np.clip(exp_back, *shared_range(float, out_dtype))
+            else: # Not via float, direct iu to iu casting
+                if (mn_in, mx_in) != (mn_out, mx_out):
+                    exp_back = np.clip(exp_back,
+                                    max(mn_in, mn_out),
+                                    min(mx_in, mx_out))
+            exp_back = exp_back.astype(out_dtype)
         # Allow for small differences in large numbers
         assert_allclose_safely(back_arr, exp_back)