[Bug]: pcolormesh issue with np.seterr(under='raise')

2sn commented 9 months ago

Bug summary

When np.seterr(under="raise") is set, the pcolormesh fails. Maybe this should be internally disabled?

Code for reproduction

import numpy as np
import matplotlib.pylab as plt

ca=array([0.5, 0.6, 0.7, 0.8, 0.9, 1. , 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7,
       1.8, 1.9, 2. ])
a3=array([0.75, 0.8 , 0.85, 0.9 , 0.95, 1.  , 1.05, 1.1 , 1.15, 1.2 , 1.25])
data=array([[1.2713495 , 1.27445031, 1.28460717, 1.29130818, 1.29799399,
        1.28663907, 1.31302497, 1.300941  , 1.30953053, 1.28866943,
        1.27942581],
       [1.26072153, 1.34242149, 1.39931996, 1.36029362, 1.27837626,
        1.26751629, 1.27348899, 1.29718153, 1.29149684, 1.29825976,
        1.29844045],
       [6.08656104, 1.30067448, 1.25214664, 1.26035875, 1.36715818,
        1.32552442, 1.39957926, 1.26443228, 1.27086277, 1.27747337,
        1.27804707],
       [3.00639509, 3.71437317, 6.09823091, 5.08496857, 1.33917387,
        1.275968  , 1.276208  , 5.85914639, 1.3672074 , 1.3054208 ,
        1.29196619],
       [1.40590936, 2.7645637 , 1.59388749, 3.00583827, 6.09541914,
        6.10693772, 1.43934882, 1.4218234 , 1.24456015, 1.2583217 ,
        4.73126518],
       [1.51797233, 4.07758869, 1.35024387, 6.16087726, 6.16087726,
        2.99797873, 1.4925019 , 6.15846305, 6.16087726, 6.17099357,
        6.16087726],
       [1.42225226, 6.20368688, 3.77458528, 1.63940466, 3.05123787,
        1.4034648 , 2.78257648, 6.1955168 , 1.63842382, 3.01200839,
        3.83237647],
       [6.12860987, 1.52526702, 6.12860987, 3.4192525 , 1.66872783,
        3.03455313, 1.32343997, 1.40334972, 2.73661128, 6.13119387,
        1.60614133],
       [2.91161999, 6.1955168 , 6.20368688, 1.56871188, 1.38913448,
        6.21100954, 3.38222841, 2.95083762, 3.06491493, 6.20368688,
        1.50694172],
       [1.74007769, 2.92087855, 6.22356158, 6.23099418, 6.22356158,
        1.53980699, 1.39476663, 6.22356158, 3.34994297, 1.57842062,
        2.96292181],
       [1.73497151, 1.7467381 , 2.92083986, 1.77700437, 6.20368688,
        6.20368688, 1.68301207, 1.48597638, 1.38461998, 3.04317972,
        3.42806827],
       [1.71963249, 1.72142379, 1.73780007, 1.72713968, 2.92199405,
        1.78053682, 6.14481403, 6.14481403, 1.55941905, 1.41043732,
        6.14481403],
       [2.91315155, 1.742551  , 1.72151811, 1.73169847, 1.7398712 ,
        2.92240585, 2.91794373, 6.15232143, 6.14481403, 6.14481403,
        1.54133219],
       [6.14481403, 6.15232143, 2.92461883, 1.73502069, 1.7412515 ,
        1.74323843, 2.92054641, 1.78489886, 1.75280085, 6.15232143,
        6.14481403],
       [6.27016475, 6.27016475, 6.27016475, 2.90887762, 2.9217294 ,
        2.92487936, 2.92513504, 2.92124614, 2.91207688, 2.90331518,
        2.86488642],
       [6.17730685, 6.17730685, 6.18773682, 6.17730685, 2.92253426,
        1.76412583, 1.76438894, 1.75630465, 1.76600374, 2.91196777,
        2.90819997]])

plt.pcolormesh(ca, a3, data.T)

Actual outcome

In [83]: plt.pcolormesh(ca, a3, data.T)
---------------------------------------------------------------------------
FloatingPointError                        Traceback (most recent call last)
Cell In[83], line 1
----> 1 plt.pcolormesh(ca, a3, data.T)

File ~/Python/lib/python3.11/site-packages/matplotlib/pyplot.py:3478, in pcolormesh(alpha, norm, cmap, vmin, vmax, shading, antialiased, data, *args, **kwargs)
   3465 @_copy_docstring_and_deprecators(Axes.pcolormesh)
   3466 def pcolormesh(
   3467     *args: ArrayLike,
   (...)
   3476     **kwargs,
   3477 ) -> QuadMesh:
-> 3478     __ret = gca().pcolormesh(
   3479         *args,
   3480         alpha=alpha,
   3481         norm=norm,
   3482         cmap=cmap,
   3483         vmin=vmin,
   3484         vmax=vmax,
   3485         shading=shading,
   3486         antialiased=antialiased,
   3487         **({"data": data} if data is not None else {}),
   3488         **kwargs,
   3489     )
   3490     sci(__ret)
   3491     return __ret

File ~/Python/lib/python3.11/site-packages/matplotlib/__init__.py:1465, in _preprocess_data.<locals>.inner(ax, data, *args, **kwargs)
   1462 @functools.wraps(func)
   1463 def inner(ax, *args, data=None, **kwargs):
   1464     if data is None:
-> 1465         return func(ax, *map(sanitize_sequence, args), **kwargs)
   1467     bound = new_sig.bind(ax, *args, **kwargs)
   1468     auto_label = (bound.arguments.get(label_namer)
   1469                   or bound.kwargs.get(label_namer))

File ~/Python/lib/python3.11/site-packages/matplotlib/axes/_axes.py:6289, in Axes.pcolormesh(self, alpha, norm, cmap, vmin, vmax, shading, antialiased, *args, **kwargs)
   6286 shading = shading.lower()
   6287 kwargs.setdefault('edgecolors', 'none')
-> 6289 X, Y, C, shading = self._pcolorargs('pcolormesh', *args,
   6290                                     shading=shading, kwargs=kwargs)
   6291 coords = np.stack([X, Y], axis=-1)
   6293 kwargs.setdefault('snap', mpl.rcParams['pcolormesh.snap'])

File ~/Python/lib/python3.11/site-packages/matplotlib/axes/_axes.py:5873, in Axes._pcolorargs(self, funcname, shading, *args, **kwargs)
   5870     return X
   5872 if ncols == Nx:
-> 5873     X = _interp_grid(X)
   5874     Y = _interp_grid(Y)
   5875 if nrows == Ny:

File ~/Python/lib/python3.11/site-packages/matplotlib/axes/_axes.py:5852, in Axes._pcolorargs.<locals>._interp_grid(X)
   5849 def _interp_grid(X):
   5850     # helper for below
   5851     if np.shape(X)[1] > 1:
-> 5852         dX = np.diff(X, axis=1)/2.
   5853         if not (np.all(dX >= 0) or np.all(dX <= 0)):
   5854             _api.warn_external(
   5855                 f"The input coordinates to {funcname} are "
   5856                 "interpreted as cell centers, but are not "
   (...)
   5859                 "edges, in which case, please supply "
   5860                 f"explicit cell edges to {funcname}.")

File ~/Python/lib/python3.11/site-packages/numpy/ma/core.py:4275, in MaskedArray.__truediv__(self, other)
   4273 if self._delegate_binop(other):
   4274     return NotImplemented
-> 4275 return true_divide(self, other)

File ~/Python/lib/python3.11/site-packages/numpy/ma/core.py:1171, in _DomainedBinaryOperation.__call__(self, a, b, *args, **kwargs)
   1169 domain = ufunc_domain.get(self.f, None)
   1170 if domain is not None:
-> 1171     m |= domain(da, db)
   1172 # Take care of the scalar case first
   1173 if not m.ndim:

File ~/Python/lib/python3.11/site-packages/numpy/ma/core.py:858, in _DomainSafeDivide.__call__(self, a, b)
    856 a, b = np.asarray(a), np.asarray(b)
    857 with np.errstate(invalid='ignore'):
--> 858     return umath.absolute(a) * self.tolerance >= umath.absolute(b)

FloatingPointError: underflow encountered in multiply

Expected outcome

a plot

Additional information

seems to be generic, but specific to this dataset

Operating system

Fedora 39

Matplotlib Version

3.8.2

Matplotlib Backend

GTK4Agg

Python version

3.11.8

Jupyter version

IPython 8.21.0

Installation

pip

jklymak commented 9 months ago

This seems a bug upstream in numpy:

np.seterr(under="raise")
x=np.arange(0, 3, 0.1)
dX = np.diff(x)
X = np.ma.array(x)
dX = np.diff(X)
dX = dX / 2.0

returns

Traceback (most recent call last):
  File "/Users/jklymak/matplotlib/testit.py", line 9, in <module>
    dX = dX / 2.0
         ~~~^~~~~
  File "/Users/jklymak/mambaforge/envs/mpl-dev/lib/python3.11/site-packages/numpy/ma/core.py", line 4275, in __truediv__
    return true_divide(self, other)
           ^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/jklymak/mambaforge/envs/mpl-dev/lib/python3.11/site-packages/numpy/ma/core.py", line 1171, in __call__
    m |= domain(da, db)
         ^^^^^^^^^^^^^^
  File "/Users/jklymak/mambaforge/envs/mpl-dev/lib/python3.11/site-packages/numpy/ma/core.py", line 858, in __call__
    return umath.absolute(a) * self.tolerance >= umath.absolute(b)
           ~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~
FloatingPointError: underflow encountered in multiply

Note that you need to divide by 2 to get the error.

A work around for us is to multiply by 0.5 instead.

tacaswell commented 9 months ago

Digging a bit more, this definitly looks like a numpy bug as the multiple that is failing is part of its internal business.

In [15]: dX / 2
---------------------------------------------------------------------------
FloatingPointError                        Traceback (most recent call last)
Cell In[15], line 1
----> 1 dX / 2

File /usr/lib/python3.11/site-packages/numpy/ma/core.py:4275, in MaskedArray.__truediv__(self, other)
   4273 if self._delegate_binop(other):
   4274     return NotImplemented
-> 4275 return true_divide(self, other)

File /usr/lib/python3.11/site-packages/numpy/ma/core.py:1171, in _DomainedBinaryOperation.__call__(self, a, b, *args, **kwargs)
   1169 domain = ufunc_domain.get(self.f, None)
   1170 if domain is not None:
-> 1171     m |= domain(da, db)
   1172 # Take care of the scalar case first
   1173 if not m.ndim:

File /usr/lib/python3.11/site-packages/numpy/ma/core.py:858, in _DomainSafeDivide.__call__(self, a, b)
    856 a, b = np.asarray(a), np.asarray(b)
    857 with np.errstate(invalid='ignore'):
--> 858     return umath.absolute(a) * self.tolerance >= umath.absolute(b)

FloatingPointError: underflow encountered in multiply

In [16]: %debug
> /usr/lib/python3.11/site-packages/numpy/ma/core.py(858)__call__()
    856         a, b = np.asarray(a), np.asarray(b)
    857         with np.errstate(invalid='ignore'):
--> 858             return umath.absolute(a) * self.tolerance >= umath.absolute(b)
    859
    860

ipdb> a
self = <numpy.ma.core._DomainSafeDivide object at 0x77c30803cc90>
a = array([0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
       0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
       0.1, 0.1, 0.1])
b = array(2)
ipdb> self.tolerance
2.2250738585072014e-308
ipdb> umath.absolute(b)
2
ipdb> umath.absolute(a)
array([0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
       0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
       0.1, 0.1, 0.1])
ipdb> umath.absolute(a) * self.tolerance
*** FloatingPointError: underflow encountered in multiply
ipdb> umath.absolute(a) * self.tolerance

matplotlib / matplotlib