Repository: matplotlib/matplotlib
Pull request: MNT: pcolormesh robust underflow
Closes #27770

This looks to be a numpy masked array problem, but multiplication doesn't trigger the spurious underflow, and probably is a touch faster as well.  



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

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

### Code for reproduction

```python
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

Hints from issue discussion:
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.

Commit message: MNT: make pcolormesh more robust to masked array underflow

Patch to review:
diff --git a/lib/matplotlib/axes/_axes.py b/lib/matplotlib/axes/_axes.py
index b1aeb87e6b45..711d930a1253 100644
--- a/lib/matplotlib/axes/_axes.py
+++ b/lib/matplotlib/axes/_axes.py
@@ -5883,7 +5883,7 @@ def _pcolorargs(self, funcname, *args, shading='auto', **kwargs):
                 def _interp_grid(X):
                     # helper for below
                     if np.shape(X)[1] > 1:
-                        dX = np.diff(X, axis=1)/2.
+                        dX = np.diff(X, axis=1) * 0.5
                         if not (np.all(dX >= 0) or np.all(dX <= 0)):
                             _api.warn_external(
                                 f"The input coordinates to {funcname} are "


Code review comments:
diff_hunk: @@ -1517,6 +1517,22 @@ def test_pcolorargs():
         ax.pcolormesh(X, Y, Z, shading='auto')
 
 
+def test_pcolormesh_underflow_error():
+    """
+    Test that underflow errors don't crop up in pcolormesh.  Probably
+    a numpy bug, but...
+    """
+    oldsettings = np.geterr()
+    np.seterr(under="raise")
path: lib/matplotlib/tests/test_axes.py
text: @user1:
Worth using the [`np.errstate` context manager](https://numpy.org/doc/stable/reference/generated/numpy.errstate.html) instead to save a couple of lines?