Merge pull request #134 from matthew-brett/better-longdouble-checks

Better longdouble checks

Make checks for longdouble implementations more specific.
Also consider POWER7 when checking for double pair implementation.
Work round some odd float128 / complex256 behavior on PPC.
Disable use of float128 nifti data type unless we have a real binary128 data type.

Author: matthew-brett

nibabel/casting.py

@@ -168,6 +168,14 @@ class FloatingError(Exception):
     pass
 
 
+def on_powerpc():
+    """ True if we are running on a Power PC platform
+
+    Has to deal with older Macs and IBM POWER7 series among others
+    """
+    return processor() == 'powerpc' or machine().startswith('ppc')
+
+
 def type_info(np_type):
     """ Return dict with min, max, nexp, nmant, width for numpy type `np_type`
 
@@ -244,13 +252,96 @@ def type_info(np_type):
         # at float64 precision. The finfo values give nexp == 15 (as for intel
         # 80) but in calculations nexp in fact appears to be 11 as for float64
         ret.update(dict(width=width))
-    elif vals == (1, 1, 16) and processor() == 'powerpc': # broken PPC
-        ret.update(dict(nmant=106, width=width))
-    else: # don't recognize the type
+    # Oh dear, we don't recognize the type information.  Try some known types
+    # and then give up. At this stage we're expecting exotic longdouble or their
+    # complex equivalent.
+    if not np_type in (np.longdouble, np.longcomplex) or width not in (16, 32):
         raise FloatingError('We had not expected type %s' % np_type)
+    if (vals == (1, 1, 16) and on_powerpc() and
+        _check_maxexp(np.longdouble, 1024)):
+        # double pair on PPC.  The _check_nmant routine does not work for this
+        # type, hence the powerpc platform check instead
+        ret.update(dict(nmant = 106, width=width))
+    elif (_check_nmant(np.longdouble, 52) and
+          _check_maxexp(np.longdouble, 11)):
+        # Got float64 despite everything
+        pass
+    elif (_check_nmant(np.longdouble, 112) and
+          _check_maxexp(np.longdouble, 16384)):
+        # binary 128, but with some busted type information. np.longcomplex
+        # seems to break here too, so we need to use np.longdouble and
+        # complexify
+        two = np.longdouble(2)
+        # See: http://matthew-brett.github.com/pydagogue/floating_point.html
+        max_val = (two ** 113 - 1) / (two ** 112) * two ** 16383
+        if np_type is np.longcomplex:
+            max_val += 0j
+        ret = dict(min = -max_val,
+                   max= max_val,
+                   nmant = 112,
+                   nexp = 15,
+                   minexp = -16382,
+                   maxexp = 16384,
+                   width = width)
+    else: # don't recognize the type
+        raise FloatingError('We had not expected long double type %s '
+                            'with info %s' % (np_type, info))
     return ret
 
 
+def _check_nmant(np_type, nmant):
+    """ True if fp type `np_type` seems to have `nmant` significand digits
+
+    Note 'digits' does not include implicit digits.  And in fact if there are no
+    implicit digits, the `nmant` number is one less than the actual digits.
+    Assumes base 2 representation.
+
+    Parameters
+    ----------
+    np_type : numpy type specifier
+        Any specifier for a numpy dtype
+    nmant : int
+        Number of digits to test against
+
+    Returns
+    -------
+    tf : bool
+        True if `nmant` is the correct number of significand digits, false
+        otherwise
+    """
+    np_type = np.dtype(np_type).type
+    max_contig = np_type(2 ** (nmant + 1)) # maximum of contiguous integers
+    tests = max_contig + np.array([-2, -1, 0, 1, 2], dtype=np_type)
+    return np.all(tests - max_contig == [-2, -1, 0, 0, 2])
+
+
+def _check_maxexp(np_type, maxexp):
+    """ True if fp type `np_type` seems to have `maxexp` maximum exponent
+
+    We're testing "maxexp" as returned by numpy. This value is set to one
+    greater than the maximum power of 2 that `np_type` can represent.
+
+    Assumes base 2 representation.  Very crude check
+
+    Parameters
+    ----------
+    np_type : numpy type specifier
+        Any specifier for a numpy dtype
+    maxexp : int
+        Maximum exponent to test against
+
+    Returns
+    -------
+    tf : bool
+        True if `maxexp` is the correct maximum exponent, False otherwise.
+    """
+    dt = np.dtype(np_type)
+    np_type = dt.type
+    two = np_type(2).reshape((1,)) # to avoid upcasting
+    return (np.isfinite(two ** (maxexp - 1)) and
+            not np.isfinite(two ** maxexp))
+
+
 def as_int(x, check=True):
     """ Return python integer representation of number
 
@@ -548,6 +639,13 @@ def best_float():
     return np.float64
 
 
+def have_binary128():
+    """ True if we have a binary128 IEEE longdouble
+    """
+    ti = type_info(np.longdouble)
+    return (ti['nmant'], ti['maxexp']) == (112, 16384)
+
+
 def ok_floats():
     """ Return floating point types sorted by precision
 

nibabel/nifti1.py

@@ -20,6 +20,7 @@
 from .quaternions import fillpositive, quat2mat, mat2quat
 from . import analyze # module import
 from .spm99analyze import SpmAnalyzeHeader
+from .casting import have_binary128
 
 # Needed for quaternion calculation
 FLOAT32_EPS_3 = -np.finfo(np.float32).eps * 3
@@ -76,13 +77,12 @@
 header_dtype = np.dtype(header_dtd)
 
 # datatypes not in analyze format, with codes
-try:
-    _float128t = np.float128
-except AttributeError:
+if have_binary128():
+    # Only enable 128 bit floats if we really have IEEE binary 128 longdoubles
+    _float128t = np.longdouble
+    _complex256t = np.longcomplex
+else:
     _float128t = np.void
-try:
-    _complex256t = np.complex256
-except AttributeError:
     _complex256t = np.void
 
 _dtdefs = ( # code, label, dtype definition, niistring

nibabel/tests/test_arraywriters.py

@@ -87,9 +87,11 @@ def test_arraywriters():
             # Byteswapped is OK
             bs_arr = arr.byteswap().newbyteorder('S')
             bs_aw = klass(bs_arr)
-            assert_array_equal(bs_arr, round_trip(bs_aw))
+            # assert against original array because POWER7 was running into
+            # trouble using the byteswapped array (bs_arr)
+            assert_array_equal(arr, round_trip(bs_aw))
             bs_aw2 = klass(bs_arr, arr.dtype)
-            assert_array_equal(bs_arr, round_trip(bs_aw2))
+            assert_array_equal(arr, round_trip(bs_aw2))
             # 2D array
             arr2 = np.reshape(arr, (2, 5))
             a2w = klass(arr2)
@@ -508,6 +510,10 @@ def test_float_int_min_max():
     for in_dt in FLOAT_TYPES:
         finf = type_info(in_dt)
         arr = np.array([finf['min'], finf['max']], dtype=in_dt)
+        # Bug in numpy 1.6.2 on PPC leading to infs - abort
+        if not np.all(np.isfinite(arr)):
+            print 'Hit PPC max -> inf bug; skip in_type %s' % in_dt
+            continue
         for out_dt in IUINT_TYPES:
             try:
                 aw = SlopeInterArrayWriter(arr, out_dt)

nibabel/tests/test_floating.py

@@ -1,14 +1,13 @@
 """ Test floating point deconstructions and floor methods
 """
-from platform import processor
-
 import numpy as np
 
 from ..casting import (floor_exact, ceil_exact, as_int, FloatingError,
-                       int_to_float, floor_log2, type_info)
+                       int_to_float, floor_log2, type_info, _check_nmant,
+                       _check_maxexp, ok_floats, on_powerpc, have_binary128)
 
 from nose import SkipTest
-from nose.tools import assert_equal, assert_raises
+from nose.tools import assert_equal, assert_raises, assert_true, assert_false
 
 IEEE_floats = [np.float32, np.float64]
 try:
@@ -80,6 +79,25 @@ def test_nmant():
         assert_equal(type_info(np.longdouble)['nmant'], 63)
 
 
+def test_check_nmant_nexp():
+    # Routine for checking number of sigificand digits and exponent
+    for t in IEEE_floats:
+        nmant = np.finfo(t).nmant
+        maxexp = np.finfo(t).maxexp
+        assert_true(_check_nmant(t, nmant))
+        assert_false(_check_nmant(t, nmant - 1))
+        assert_false(_check_nmant(t, nmant + 1))
+        assert_true(_check_maxexp(t, maxexp))
+        assert_false(_check_maxexp(t, maxexp - 1))
+        assert_false(_check_maxexp(t, maxexp + 1))
+    # Check against type_info
+    for t in ok_floats():
+        ti = type_info(t)
+        if ti['nmant'] != 106: # This check does not work for PPC double pair
+            assert_true(_check_nmant(t, ti['nmant']))
+        assert_true(_check_maxexp(t, ti['maxexp']))
+
+
 def test_as_int():
     # Integer representation of number
     assert_equal(as_int(2.0), 2)
@@ -214,7 +232,7 @@ def test_floor_exact():
                 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':
+        if t is np.longdouble and on_powerpc():
             continue
         # Confirm to ourselves that 2**(nmant+1) can't be exactly represented
         iv = 2**(nmant+1)
@@ -240,3 +258,12 @@ def test_floor_exact():
                 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)
+
+
+def test_usable_binary128():
+    # Check for usable binary128
+    yes = have_binary128()
+    abf = np.longdouble(2) ** 16383
+    assert_equal(yes,
+                 abf.dtype.itemsize == 16 and
+                 np.isfinite(abf))

nibabel/tests/test_nifti1.py

@@ -14,7 +14,7 @@
 
 import numpy as np
 
-from ..casting import type_info
+from ..casting import type_info, have_binary128
 from ..tmpdirs import InTemporaryDirectory
 from ..spatialimages import HeaderDataError
 from ..affines import from_matvec
@@ -309,6 +309,14 @@ def test_binblock_is_file(self):
         hdr.write_to(str_io)
         assert_equal(str_io.getvalue(), hdr.binaryblock + ZEROB * 4)
 
+    def test_float128(self):
+        hdr = self.header_class()
+        if have_binary128():
+            hdr.set_data_dtype(np.longdouble)
+            assert_equal(hdr.get_data_dtype().type, np.longdouble)
+        else:
+            assert_raises(HeaderDataError, hdr.set_data_dtype, np.longdouble)
+
 
 class TestNifti1Image(tana.TestAnalyzeImage):
     # Run analyze-flavor spatialimage tests

nibabel/tests/test_scaling.py

@@ -224,6 +224,11 @@ def check_int_a2f(in_type, out_type):
     this_min, this_max = info['min'], info['max']
     if not in_type in np.sctypes['complex']:
         data = np.array([this_min, this_max], in_type)
+        # Bug in numpy 1.6.2 on PPC leading to infs - abort
+        if not np.all(np.isfinite(data)):
+            if DEBUG:
+                print 'Hit PPC max -> inf bug; skip in_type %s' % in_type
+            return
     else: # Funny behavior with complex256
         data = np.zeros((2,), in_type)
         data[0] = this_min + 0j