From d3382a5e7e3c45a2fc9840509de0c8eadc25fa07 Mon Sep 17 00:00:00 2001 From: Stephannie Jimenez Gacha Date: Mon, 19 Jul 2021 14:01:40 -0500 Subject: [PATCH] ENH: Add smallest_normal and smallest_subnormal attributes to finfo (#18536) * ENH: Add smallest_normal and smallest_subnormal attributes to finfo class * ENH: Add smallest_normal and smallest_subnormal attributes to MachAr * BUG: Fix minor bug in the calc of smallest_subnormal attribute * ENH: Use nextafter to calculate smallest subnormal in finfo * ENH: Use nextafter for calculating all subnormal numbers * ENH: Add warning in case the subnormal number is zero and hardcode double double value * FIX: Fix docstring in new function * FIX: Add stacklevel parameter to warning * ENH: Move smallest_subnormal to a property and add a test * FIX: Only hardcode the subnormal value for double double * FIX: Fix subnormal warning test * FIX: Fix the hardcoded value for double double type * FIX: Fix failing test by supressing an overflow warning * ENH: Fix the value for smallest_normal and tiny in double double * BUG: Fix failing tests when tiny is NaN * BUG: Fix failing test when tiny is NaN * ENH: add review fixes * ENH: Raise TypeError for smallest_subnormal/tiny value in longdouble * FIX: Replace TypeError with warning * BUG: Fix tests * ENH: Add release note * ENH: Add review changes * BUG: Fix wrong value in smallest_subnormal Co-authored-by: Sebastian Berg * BUG: Fix wrong value in subnormal number * ENH: Bind the smallest_normal warning to the finfo class * BUG: Revert warning supression in tests * ENH: Add review changes * BUG: Fix failing test * BUG: Fix test identation * BUG: Fix typo * Update numpy/core/tests/test_getlimits.py Co-authored-by: Sebastian Berg --- .../upcoming_changes/18536.improvement.rst | 7 + numpy/core/getlimits.py | 214 ++++++++++++++---- numpy/core/machar.py | 14 +- numpy/core/tests/test_getlimits.py | 26 ++- numpy/core/tests/test_longdouble.py | 1 + numpy/core/tests/test_numeric.py | 17 +- numpy/core/tests/test_scalarmath.py | 6 +- numpy/core/tests/test_ufunc.py | 18 +- numpy/core/tests/test_umath.py | 10 +- 9 files changed, 249 insertions(+), 64 deletions(-) create mode 100644 doc/release/upcoming_changes/18536.improvement.rst diff --git a/doc/release/upcoming_changes/18536.improvement.rst b/doc/release/upcoming_changes/18536.improvement.rst new file mode 100644 index 000000000000..8693916dbdba --- /dev/null +++ b/doc/release/upcoming_changes/18536.improvement.rst @@ -0,0 +1,7 @@ +Add ``smallest_normal`` and ``smallest_subnormal`` attributes to `finfo` +------------------------------------------------------------------------- + +The attributes ``smallest_normal`` and ``smallest_subnormal`` are available as +an extension of `finfo` class for any floating-point data type. To use these +new attributes, write ``np.finfo(np.float64).smallest_normal`` or +``np.finfo(np.float64).smallest_subnormal``. diff --git a/numpy/core/getlimits.py b/numpy/core/getlimits.py index fcb73e8ba3a4..0f7031bac96a 100644 --- a/numpy/core/getlimits.py +++ b/numpy/core/getlimits.py @@ -9,8 +9,8 @@ from .overrides import set_module from . import numeric from . import numerictypes as ntypes -from .numeric import array, inf -from .umath import log10, exp2 +from .numeric import array, inf, NaN +from .umath import log10, exp2, nextafter, isnan from . import umath @@ -29,32 +29,96 @@ def _fr1(a): a.shape = () return a + class MachArLike: """ Object to simulate MachAr instance """ - - def __init__(self, - ftype, - *, eps, epsneg, huge, tiny, ibeta, **kwargs): - params = _MACHAR_PARAMS[ftype] - float_conv = lambda v: array([v], ftype) - float_to_float = lambda v : _fr1(float_conv(v)) - float_to_str = lambda v: (params['fmt'] % array(_fr0(v)[0], ftype)) - - self.title = params['title'] + def __init__(self, ftype, *, eps, epsneg, huge, tiny, + ibeta, smallest_subnormal=None, **kwargs): + self.params = _MACHAR_PARAMS[ftype] + self.ftype = ftype + self.title = self.params['title'] # Parameter types same as for discovered MachAr object. - self.epsilon = self.eps = float_to_float(eps) - self.epsneg = float_to_float(epsneg) - self.xmax = self.huge = float_to_float(huge) - self.xmin = self.tiny = float_to_float(tiny) - self.ibeta = params['itype'](ibeta) + if not smallest_subnormal: + self._smallest_subnormal = nextafter( + self.ftype(0), self.ftype(1), dtype=self.ftype) + else: + self._smallest_subnormal = smallest_subnormal + self.epsilon = self.eps = self._float_to_float(eps) + self.epsneg = self._float_to_float(epsneg) + self.xmax = self.huge = self._float_to_float(huge) + self.xmin = self._float_to_float(tiny) + self.smallest_normal = self.tiny = self._float_to_float(tiny) + self.ibeta = self.params['itype'](ibeta) self.__dict__.update(kwargs) self.precision = int(-log10(self.eps)) - self.resolution = float_to_float(float_conv(10) ** (-self.precision)) - self._str_eps = float_to_str(self.eps) - self._str_epsneg = float_to_str(self.epsneg) - self._str_xmin = float_to_str(self.xmin) - self._str_xmax = float_to_str(self.xmax) - self._str_resolution = float_to_str(self.resolution) + self.resolution = self._float_to_float( + self._float_conv(10) ** (-self.precision)) + self._str_eps = self._float_to_str(self.eps) + self._str_epsneg = self._float_to_str(self.epsneg) + self._str_xmin = self._float_to_str(self.xmin) + self._str_xmax = self._float_to_str(self.xmax) + self._str_resolution = self._float_to_str(self.resolution) + self._str_smallest_normal = self._float_to_str(self.xmin) + + @property + def smallest_subnormal(self): + """Return the value for the smallest subnormal. + + Returns + ------- + smallest_subnormal : float + value for the smallest subnormal. + + Warns + ----- + UserWarning + If the calculated value for the smallest subnormal is zero. + """ + # Check that the calculated value is not zero, in case it raises a + # warning. + value = self._smallest_subnormal + if self.ftype(0) == value: + warnings.warn( + 'The value of the smallest subnormal for {} type ' + 'is zero.'.format(self.ftype), UserWarning, stacklevel=2) + + return self._float_to_float(value) + + @property + def _str_smallest_subnormal(self): + """Return the string representation of the smallest subnormal.""" + return self._float_to_str(self.smallest_subnormal) + + def _float_to_float(self, value): + """Converts float to float. + + Parameters + ---------- + value : float + value to be converted. + """ + return _fr1(self._float_conv(value)) + + def _float_conv(self, value): + """Converts float to conv. + + Parameters + ---------- + value : float + value to be converted. + """ + return array([value], self.ftype) + + def _float_to_str(self, value): + """Converts float to str. + + Parameters + ---------- + value : float + value to be converted. + """ + return self.params['fmt'] % array(_fr0(value)[0], self.ftype) + _convert_to_float = { ntypes.csingle: ntypes.single, @@ -91,6 +155,7 @@ def _register_type(machar, bytepat): _KNOWN_TYPES[bytepat] = machar _float_ma = {} + def _register_known_types(): # Known parameters for float16 # See docstring of MachAr class for description of parameters. @@ -208,23 +273,27 @@ def _register_known_types(): # https://en.wikipedia.org/wiki/Quadruple-precision_floating-point_format#Double-double_arithmetic # These numbers have the same exponent range as float64, but extended number of # digits in the significand. - huge_dd = (umath.nextafter(ld(inf), ld(0)) - if hasattr(umath, 'nextafter') # Missing on some platforms? - else float64_ma.huge) + huge_dd = nextafter(ld(inf), ld(0), dtype=ld) + # As the smallest_normal in double double is so hard to calculate we set + # it to NaN. + smallest_normal_dd = NaN + # Leave the same value for the smallest subnormal as double + smallest_subnormal_dd = ld(nextafter(0., 1.)) float_dd_ma = MachArLike(ld, - machep=-105, - negep=-106, - minexp=-1022, - maxexp=1024, - it=105, - iexp=11, - ibeta=2, - irnd=5, - ngrd=0, - eps=exp2(ld(-105)), - epsneg= exp2(ld(-106)), - huge=huge_dd, - tiny=exp2(ld(-1022))) + machep=-105, + negep=-106, + minexp=-1022, + maxexp=1024, + it=105, + iexp=11, + ibeta=2, + irnd=5, + ngrd=0, + eps=exp2(ld(-105)), + epsneg=exp2(ld(-106)), + huge=huge_dd, + tiny=smallest_normal_dd, + smallest_subnormal=smallest_subnormal_dd) # double double; low, high order (e.g. PPC 64) _register_type(float_dd_ma, b'\x9a\x99\x99\x99\x99\x99Y<\x9a\x99\x99\x99\x99\x99\xb9\xbf') @@ -341,8 +410,13 @@ class finfo: The approximate decimal resolution of this type, i.e., ``10**-precision``. tiny : float - The smallest positive floating point number with full precision - (see Notes). + An alias for `smallest_normal`, kept for backwards compatibility. + smallest_normal : float + The smallest positive floating point number with 1 as leading bit in + the mantissa following IEEE-754 (see Notes). + smallest_subnormal : float + The smallest positive floating point number with 0 as leading bit in + the mantissa following IEEE-754. Parameters ---------- @@ -363,12 +437,12 @@ class finfo: impacts import times. These objects are cached, so calling ``finfo()`` repeatedly inside your functions is not a problem. - Note that ``tiny`` is not actually the smallest positive representable - value in a NumPy floating point type. As in the IEEE-754 standard [1]_, - NumPy floating point types make use of subnormal numbers to fill the - gap between 0 and ``tiny``. However, subnormal numbers may have - significantly reduced precision [2]_. - + Note that ``smallest_normal`` is not actually the smallest positive + representable value in a NumPy floating point type. As in the IEEE-754 + standard [1]_, NumPy floating point types make use of subnormal numbers to + fill the gap between 0 and ``smallest_normal``. However, subnormal numbers + may have significantly reduced precision [2]_. + References ---------- .. [1] IEEE Standard for Floating-Point Arithmetic, IEEE Std 754-2008, @@ -420,7 +494,7 @@ def _init(self, dtype): 'maxexp', 'minexp', 'negep', 'machep']: setattr(self, word, getattr(machar, word)) - for word in ['tiny', 'resolution', 'epsneg']: + for word in ['resolution', 'epsneg', 'smallest_subnormal']: setattr(self, word, getattr(machar, word).flat[0]) self.bits = self.dtype.itemsize * 8 self.max = machar.huge.flat[0] @@ -434,6 +508,8 @@ def _init(self, dtype): self._str_epsneg = machar._str_epsneg.strip() self._str_eps = machar._str_eps.strip() self._str_resolution = machar._str_resolution.strip() + self._str_smallest_normal = machar._str_smallest_normal.strip() + self._str_smallest_subnormal = machar._str_smallest_subnormal.strip() return self def __str__(self): @@ -446,6 +522,8 @@ def __str__(self): 'minexp = %(minexp)6s tiny = %(_str_tiny)s\n' 'maxexp = %(maxexp)6s max = %(_str_max)s\n' 'nexp = %(nexp)6s min = -max\n' + 'smallest_normal = %(_str_smallest_normal)s ' + 'smallest_subnormal = %(_str_smallest_subnormal)s\n' '---------------------------------------------------------------\n' ) return fmt % self.__dict__ @@ -457,6 +535,46 @@ def __repr__(self): return (("%(klass)s(resolution=%(resolution)s, min=-%(_str_max)s," " max=%(_str_max)s, dtype=%(dtype)s)") % d) + @property + def smallest_normal(self): + """Return the value for the smallest normal. + + Returns + ------- + smallest_normal : float + Value for the smallest normal. + + Warns + ----- + UserWarning + If the calculated value for the smallest normal is requested for + double-double. + """ + # This check is necessary because the value for smallest_normal is + # platform dependent for longdouble types. + if isnan(self.machar.smallest_normal.flat[0]): + warnings.warn( + 'The value of smallest normal is undefined for double double', + UserWarning, stacklevel=2) + return self.machar.smallest_normal.flat[0] + + @property + def tiny(self): + """Return the value for tiny, alias of smallest_normal. + + Returns + ------- + tiny : float + Value for the smallest normal, alias of smallest_normal. + + Warns + ----- + UserWarning + If the calculated value for the smallest normal is requested for + double-double. + """ + return self.smallest_normal + @set_module('numpy') class iinfo: diff --git a/numpy/core/machar.py b/numpy/core/machar.py index 55285fe5928f..04dad4d77d28 100644 --- a/numpy/core/machar.py +++ b/numpy/core/machar.py @@ -56,13 +56,19 @@ class MachAr: epsilon : float Same as `eps`. tiny : float - Same as `xmin`. + An alias for `smallest_normal`, kept for backwards compatibility. huge : float Same as `xmax`. precision : float ``- int(-log10(eps))`` resolution : float ``- 10**(-precision)`` + smallest_normal : float + The smallest positive floating point number with 1 as leading bit in + the mantissa following IEEE-754. Same as `xmin`. + smallest_subnormal : float + The smallest positive floating point number with 0 as leading bit in + the mantissa following IEEE-754. Parameters ---------- @@ -293,6 +299,8 @@ def _do_init(self, float_conv, int_conv, float_to_float, float_to_str, title): else: xmax = xmax * beta + smallest_subnormal = abs(xmin / beta ** (it)) + self.ibeta = ibeta self.it = it self.negep = negep @@ -316,6 +324,8 @@ def _do_init(self, float_conv, int_conv, float_to_float, float_to_str, title): self.epsilon = self.eps self.tiny = self.xmin self.huge = self.xmax + self.smallest_normal = self.xmin + self.smallest_subnormal = float_to_float(smallest_subnormal) import math self.precision = int(-math.log10(float_to_float(self.eps))) @@ -333,6 +343,8 @@ def __str__(self): 'negep =%(negep)s epsneg=%(_str_epsneg)s (beta**epsneg)\n' 'minexp=%(minexp)s xmin=%(_str_xmin)s (beta**minexp == tiny)\n' 'maxexp=%(maxexp)s xmax=%(_str_xmax)s ((1-epsneg)*beta**maxexp == huge)\n' + 'smallest_normal=%(smallest_normal)s ' + 'smallest_subnormal=%(smallest_subnormal)s\n' '---------------------------------------------------------------------\n' ) return fmt % self.__dict__ diff --git a/numpy/core/tests/test_getlimits.py b/numpy/core/tests/test_getlimits.py index bcf8cf659b79..de7b3e769394 100644 --- a/numpy/core/tests/test_getlimits.py +++ b/numpy/core/tests/test_getlimits.py @@ -1,6 +1,7 @@ """ Test functions for limits module. """ +import warnings import numpy as np from numpy.core import finfo, iinfo from numpy import half, single, double, longdouble @@ -47,7 +48,8 @@ def test_basic(self): for dt1, dt2 in dts: for attr in ('bits', 'eps', 'epsneg', 'iexp', 'machar', 'machep', 'max', 'maxexp', 'min', 'minexp', 'negep', 'nexp', - 'nmant', 'precision', 'resolution', 'tiny'): + 'nmant', 'precision', 'resolution', 'tiny', + 'smallest_normal', 'smallest_subnormal'): assert_equal(getattr(finfo(dt1), attr), getattr(finfo(dt2), attr), attr) assert_raises(ValueError, finfo, 'i4') @@ -112,6 +114,28 @@ def test_known_types(): assert_ma_equal(ld_ma, _float_ma[128]) +def test_subnormal_warning(): + """Test that the subnormal is zero warning is not being raised.""" + with np.errstate(all='ignore'): + ld_ma = _discovered_machar(np.longdouble) + bytes = np.dtype(np.longdouble).itemsize + with warnings.catch_warnings(record=True) as w: + warnings.simplefilter('always') + if (ld_ma.it, ld_ma.maxexp) == (63, 16384) and bytes in (12, 16): + # 80-bit extended precision + ld_ma.smallest_subnormal + assert len(w) == 0 + elif (ld_ma.it, ld_ma.maxexp) == (112, 16384) and bytes == 16: + # IEE 754 128-bit + ld_ma.smallest_subnormal + assert len(w) == 0 + else: + # Double double + ld_ma.smallest_subnormal + # This test may fail on some platforms + assert len(w) == 0 + + def test_plausible_finfo(): # Assert that finfo returns reasonable results for all types for ftype in np.sctypes['float'] + np.sctypes['complex']: diff --git a/numpy/core/tests/test_longdouble.py b/numpy/core/tests/test_longdouble.py index acef995f3f05..1a54e62d8b3d 100644 --- a/numpy/core/tests/test_longdouble.py +++ b/numpy/core/tests/test_longdouble.py @@ -8,6 +8,7 @@ ) from numpy.core.tests._locales import CommaDecimalPointLocale + LD_INFO = np.finfo(np.longdouble) longdouble_longer_than_double = (LD_INFO.eps < np.finfo(np.double).eps) diff --git a/numpy/core/tests/test_numeric.py b/numpy/core/tests/test_numeric.py index 1315364bcabd..9551855ad544 100644 --- a/numpy/core/tests/test_numeric.py +++ b/numpy/core/tests/test_numeric.py @@ -646,7 +646,7 @@ def test_floating_exceptions(self): if np.dtype(ftype).kind == 'f': # Get some extreme values for the type fi = np.finfo(ftype) - ft_tiny = fi.tiny + ft_tiny = fi.machar.tiny ft_max = fi.max ft_eps = fi.eps underflow = 'underflow' @@ -655,7 +655,7 @@ def test_floating_exceptions(self): # 'c', complex, corresponding real dtype rtype = type(ftype(0).real) fi = np.finfo(rtype) - ft_tiny = ftype(fi.tiny) + ft_tiny = ftype(fi.machar.tiny) ft_max = ftype(fi.max) ft_eps = ftype(fi.eps) # The complex types raise different exceptions @@ -664,10 +664,13 @@ def test_floating_exceptions(self): overflow = 'overflow' invalid = 'invalid' - self.assert_raises_fpe(underflow, - lambda a, b: a/b, ft_tiny, ft_max) - self.assert_raises_fpe(underflow, - lambda a, b: a*b, ft_tiny, ft_tiny) + # The value of tiny for double double is NaN, so we need to + # pass the assert + if not np.isnan(ft_tiny): + self.assert_raises_fpe(underflow, + lambda a, b: a/b, ft_tiny, ft_max) + self.assert_raises_fpe(underflow, + lambda a, b: a*b, ft_tiny, ft_tiny) self.assert_raises_fpe(overflow, lambda a, b: a*b, ft_max, ftype(2)) self.assert_raises_fpe(overflow, @@ -898,7 +901,7 @@ def test_promote_types_strings(self, swap, string_dtype): promote_types = np.promote_types S = string_dtype - + # Promote numeric with unsized string: assert_equal(promote_types('bool', S), np.dtype(S+'5')) assert_equal(promote_types('b', S), np.dtype(S+'4')) diff --git a/numpy/core/tests/test_scalarmath.py b/numpy/core/tests/test_scalarmath.py index 09a734284a76..5981225c1f05 100644 --- a/numpy/core/tests/test_scalarmath.py +++ b/numpy/core/tests/test_scalarmath.py @@ -670,8 +670,10 @@ def _test_abs_func(self, absfunc, test_dtype): x = test_dtype(np.finfo(test_dtype).max) assert_equal(absfunc(x), x.real) - x = test_dtype(np.finfo(test_dtype).tiny) - assert_equal(absfunc(x), x.real) + with suppress_warnings() as sup: + sup.filter(UserWarning) + x = test_dtype(np.finfo(test_dtype).tiny) + assert_equal(absfunc(x), x.real) x = test_dtype(np.finfo(test_dtype).min) assert_equal(absfunc(x), -x.real) diff --git a/numpy/core/tests/test_ufunc.py b/numpy/core/tests/test_ufunc.py index a47f1df49871..0251f21a9899 100644 --- a/numpy/core/tests/test_ufunc.py +++ b/numpy/core/tests/test_ufunc.py @@ -12,7 +12,7 @@ from numpy.testing import ( assert_, assert_equal, assert_raises, assert_array_equal, assert_almost_equal, assert_array_almost_equal, assert_no_warnings, - assert_allclose, HAS_REFCOUNT, + assert_allclose, HAS_REFCOUNT, suppress_warnings ) from numpy.compat import pickle @@ -583,7 +583,13 @@ def test_true_divide(self): else: tgt = float(x)/float(y) rtol = max(np.finfo(dtout).resolution, 1e-15) - atol = max(np.finfo(dtout).tiny, 3e-308) + # The value of tiny for double double is NaN + with suppress_warnings() as sup: + sup.filter(UserWarning) + if not np.isnan(np.finfo(dtout).tiny): + atol = max(np.finfo(dtout).tiny, 3e-308) + else: + atol = 3e-308 # Some test values result in invalid for float16. with np.errstate(invalid='ignore'): res = np.true_divide(x, y, dtype=dtout) @@ -596,7 +602,13 @@ def test_true_divide(self): dtout = np.dtype(tcout) tgt = complex(x)/complex(y) rtol = max(np.finfo(dtout).resolution, 1e-15) - atol = max(np.finfo(dtout).tiny, 3e-308) + # The value of tiny for double double is NaN + with suppress_warnings() as sup: + sup.filter(UserWarning) + if not np.isnan(np.finfo(dtout).tiny): + atol = max(np.finfo(dtout).tiny, 3e-308) + else: + atol = 3e-308 res = np.true_divide(x, y, dtype=dtout) if not np.isfinite(res): continue diff --git a/numpy/core/tests/test_umath.py b/numpy/core/tests/test_umath.py index ebd61d199612..8d7fa1e0a2a3 100644 --- a/numpy/core/tests/test_umath.py +++ b/numpy/core/tests/test_umath.py @@ -3463,8 +3463,14 @@ def test_nextafterl(): def test_nextafter_0(): for t, direction in itertools.product(np.sctypes['float'], (1, -1)): - tiny = np.finfo(t).tiny - assert_(0. < direction * np.nextafter(t(0), t(direction)) < tiny) + # The value of tiny for double double is NaN, so we need to pass the + # assert + with suppress_warnings() as sup: + sup.filter(UserWarning) + if not np.isnan(np.finfo(t).tiny): + tiny = np.finfo(t).tiny + assert_( + 0. < direction * np.nextafter(t(0), t(direction)) < tiny) assert_equal(np.nextafter(t(0), t(direction)) / t(2.1), direction * 0.0) def _test_spacing(t):