DOC+RF: more comments / rewrite of scaling checks

I found the scaling checks difficult to follow, so I added some comments
and refactored the scaling checks for the array_to_file routine.

Author: matthew-brett

nibabel/casting.py

@@ -14,7 +14,10 @@ class CastingError(Exception):
     pass
 
 
-#: Test for VC truncation when casting floats to uint64
+# Test for VC truncation when casting floats to uint64
+# Christoph Gohlke says this is so for MSVC <= 2010 because VC is using x87
+# instructions; see:
+# https://github.com/scipy/scipy/blob/99bb8411f6391d921cb3f4e56619291e91ddf43b/scipy/ndimage/tests/test_datatypes.py#L51
 _test_val = 2**63 + 2**11 # Should be exactly representable in float64
 TRUNC_UINT64 = np.float64(_test_val).astype(np.uint64) != _test_val
 

nibabel/tests/test_arraywriters.py

@@ -138,24 +138,35 @@ def test_no_scaling():
         with suppress_warnings():
             back_arr = round_trip(aw)
         exp_back = arr.copy()
+        # If converting to floating point type, casting is direct.
+        # Otherwise we will need to do float-(u)int casting at some point.
         if out_dtype in IUINT_TYPES:
             if in_dtype in CFLOAT_TYPES:
+                # Working precision is (at least) float
                 with suppress_warnings():
                     exp_back = exp_back.astype(float)
+                # Float to iu conversion will always round, clip
                 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
+                    # Clip to shared range of working precision
+                    exp_back = np.clip(exp_back,
+                                       *shared_range(float, out_dtype))
+            else:  # iu input and output type
+                # No scaling, never gets converted to float.
+                # Does get clipped to range of output type
                 mn_out, mx_out = _dt_min_max(out_dtype)
                 if (mn_in, mx_in) != (mn_out, mx_out):
+                    # Use smaller of input, output range to avoid np.clip
+                    # upcasting the array because of large clip limits.
                     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
+        elif in_dtype in COMPLEX_TYPES:
+            # always cast to real from complex
             with suppress_warnings():
                 exp_back = exp_back.astype(float)
         exp_back = exp_back.astype(out_dtype)

nibabel/tests/test_spm99analyze.py

@@ -343,25 +343,31 @@ def test_no_scaling(self):
             with suppress_warnings():  # invalid mult
                 back_arr = rt_img.get_data()
             exp_back = arr.copy()
+            # If converting to floating point type, casting is direct.
+            # Otherwise we will need to do float-(u)int casting at some point
             if out_dtype in IUINT_TYPES:
                 if in_dtype in FLOAT_TYPES:
                     # Working precision is (at least) float
                     exp_back = exp_back.astype(float)
+                    # Float to iu conversion will always round, clip
                     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
+                else:  # iu input and output type
+                    # No scaling, never gets converted to float.
+                    # Does get clipped to range of output type
                     mn_out, mx_out = _dt_min_max(out_dtype)
                     if (mn_in, mx_in) != (mn_out, mx_out):
+                        # Use smaller of input, output range to avoid np.clip
+                        # upcasting the array because of large clip limits.
                         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:
+                # always cast to real from complex
                 exp_back = exp_back.astype(out_dtype)
             else:
                 # Cast to working precision

nibabel/tests/test_utils.py

@@ -536,13 +536,12 @@ 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)
+        # 0 when scaled to output will also be the output value for NaN
         nan_fill = -intercept / divslope
         if out_dtype in IUINT_TYPES:
             nan_fill = np.round(nan_fill)
+        # nan2zero will check whether 0 in scaled to a valid value in output
         if (in_dtype in CFLOAT_TYPES and not mn_out <= nan_fill <= mx_out):
             assert_raises(ValueError,
                           array_to_file,
@@ -552,29 +551,36 @@ def test_a2f_scaled_unscaled():
                           divslope=divslope,
                           intercept=intercept)
             continue
-        with suppress_warnings():  # cast to real
+        with suppress_warnings():
             back_arr = write_return(arr, fobj,
                                     out_dtype=out_dtype,
                                     divslope=divslope,
                                     intercept=intercept)
             exp_back = arr.copy()
-            if iu_via_float:
-                if out_dtype in IUINT_TYPES:
+            if (in_dtype in IUINT_TYPES and
+                out_dtype in IUINT_TYPES and
+                (intercept, divslope) == (0, 1)):
+                # Direct iu to iu casting.
+                # Need to clip if ranges not the same.
+                # Use smaller of input, output range to avoid np.clip upcasting
+                # the array because of large clip limits.
+                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))
+            else: # Need to deal with nans, casting to float, clipping
+                if in_dtype in CFLOAT_TYPES and 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:
+                if (exp_back.dtype.type in CFLOAT_TYPES and
+                    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)