[Bug]: fill_between gives IndexError with numpy 1.24.0.dev

[pllim]

Bug summary

There is a test in astropy that started breaking when using matplotlib + numpy 1.24.0.dev0+896.g5ecaf36cd . Same error using stable matplotlib and matplotlib from the nightly wheel.

The actual test is at https://github.com/astropy/astropy/blob/75f9b60b5521ed8dae08611ddf766c95ce421801/astropy/visualization/tests/test_units.py#L24 but I will give reproducible example without astropy below.

Example log with the astropy failure: https://github.com/astropy/astropy/actions/runs/3198205490/jobs/5222440278

Code for reproduction

import matplotlib
import matplotlib.pyplot as plt

matplotlib.use('agg')
plt.figure()
plt.plot([1, 2, 3], [3, 4, 5], label='label')
plt.plot([1.05, 2.10, 3.15], [3.050, 3.025, 3.010])
plt.fill_between([1, 3], [3, 5], [3.050, 3.010])  # IndexError

Actual outcome

File .../matplotlib/pyplot.py:2514, in fill_between(x, y1, y2, where, interpolate, step, data, **kwargs)
   2510 @_copy_docstring_and_deprecators(Axes.fill_between)
   2511 def fill_between(
   2512         x, y1, y2=0, where=None, interpolate=False, step=None, *,
   2513         data=None, **kwargs):
-> 2514     return gca().fill_between(
   2515         x, y1, y2=y2, where=where, interpolate=interpolate, step=step,
   2516         **({"data": data} if data is not None else {}), **kwargs)

File .../matplotlib/__init__.py:1423, in _preprocess_data.<locals>.inner(ax, data, *args, **kwargs)
   1420 @functools.wraps(func)
   1421 def inner(ax, *args, data=None, **kwargs):
   1422     if data is None:
-> 1423         return func(ax, *map(sanitize_sequence, args), **kwargs)
   1425     bound = new_sig.bind(ax, *args, **kwargs)
   1426     auto_label = (bound.arguments.get(label_namer)
   1427                   or bound.kwargs.get(label_namer))

File .../matplotlib/axes/_axes.py:5335, in Axes.fill_between(self, x, y1, y2, where, interpolate, step, **kwargs)
   5333 def fill_between(self, x, y1, y2=0, where=None, interpolate=False,
   5334                  step=None, **kwargs):
-> 5335     return self._fill_between_x_or_y(
   5336         "x", x, y1, y2,
   5337         where=where, interpolate=interpolate, step=step, **kwargs)

File .../matplotlib/axes/_axes.py:5263, in Axes._fill_between_x_or_y(self, ind_dir, ind, dep1, dep2, where, interpolate, step, **kwargs)
   5259 ind, dep1, dep2 = np.broadcast_arrays(
   5260     np.atleast_1d(ind), dep1, dep2, subok=True)
   5262 polys = []
-> 5263 for idx0, idx1 in cbook.contiguous_regions(where):
   5264     indslice = ind[idx0:idx1]
   5265     dep1slice = dep1[idx0:idx1]

File .../matplotlib/cbook/__init__.py:1333, in contiguous_regions(mask)
   1330     return []
   1332 # Find the indices of region changes, and correct offset
-> 1333 idx, = np.nonzero(mask[:-1] != mask[1:])
   1334 idx += 1
   1336 # List operations are faster for moderately sized arrays

IndexError: too many indices for array: array is 0-dimensional, but 1 were indexed

Expected outcome

No error.

Additional information

Only with numpy 1.24.0.dev0+896.g5ecaf36cd but not with stable numpy.

Operating system

Debian or Ubuntu

Matplotlib Version

stable or nightly

Matplotlib Backend

agg

Python version

3.10 and 3.11

Jupyter version

N/A

Installation

pip

[tacaswell]

This seems to be due to a change in how masked arrays work in numpy. Under some condition we are getting masks which are a single bool rather than an array of bools.

Looking at %debug you can see

ipdb> p functools.reduce(np.logical_or, map(np.ma.getmask, [ind, dep1, dep2]))
False
ipdb> ind
masked_array(data=[1, 3],
             mask=False,
       fill_value=999999)
ipdb> dep1
masked_array(data=[3, 5],
             mask=False,
       fill_value=999999)
ipdb> dep2
masked_array(data=[3.05, 3.01],
             mask=False,
       fill_value=1e+20)

which means that we pass a bool into

-> 5281         for idx0, idx1 in cbook.contiguous_regions(where):

which also fails on older versions of numpy (which seems both reasonable and correct).

That is far as made it debugging this. Not sure if this is a regression upstream in numpy we should report or if this is intentional change we need to adapt to.

@melissawm as our local numpy dev.

[QuLogic]

This bisects to numpy/numpy@44c8da9, which coincidentally fixes numpy/numpy#19332 which is a bug I reported. np.ma.masked_invalid now mirrors np.ma.masked_where as documented.

The np.ma.masked_invalid call is a couple lines earlier, but strangely, we don't use copy=False. I guess this is just another inconsistency between the two.

[QuLogic]

Something like this might work?

diff --git a/lib/matplotlib/axes/_axes.py b/lib/matplotlib/axes/_axes.py
index 64ee4a5122..cd18072fe4 100644
--- a/lib/matplotlib/axes/_axes.py
+++ b/lib/matplotlib/axes/_axes.py
@@ -5274,8 +5274,8 @@ default: :rc:`scatter.edgecolors`
         where = where & ~functools.reduce(
             np.logical_or, map(np.ma.getmask, [ind, dep1, dep2]))
 
-        ind, dep1, dep2 = np.broadcast_arrays(
-            np.atleast_1d(ind), dep1, dep2, subok=True)
+        ind, dep1, dep2, where = np.broadcast_arrays(
+            np.atleast_1d(ind), dep1, dep2, where, subok=True)
 
         polys = []
         for idx0, idx1 in cbook.contiguous_regions(where):

[QuLogic]

Or actually, since we definitely want an array here, this is probably better:

diff --git a/lib/matplotlib/axes/_axes.py b/lib/matplotlib/axes/_axes.py
index 64ee4a5122..e6d8a73689 100644
--- a/lib/matplotlib/axes/_axes.py
+++ b/lib/matplotlib/axes/_axes.py
@@ -5272,7 +5272,7 @@ default: :rc:`scatter.edgecolors`
                 raise ValueError(f"where size ({where.size}) does not match "
                                  f"{ind_dir} size ({ind.size})")
         where = where & ~functools.reduce(
-            np.logical_or, map(np.ma.getmask, [ind, dep1, dep2]))
+            np.logical_or, map(np.ma.getmaskarray, [ind, dep1, dep2]))
 
         ind, dep1, dep2 = np.broadcast_arrays(
             np.atleast_1d(ind), dep1, dep2, subok=True)

[pllim]

Thanks for the quick fix! 🙇‍♀️