Repository: numpy/numpy
Pull request: BUG: Raise correct errors in boolean indexing fast path
Previously the logic assumed that an index was valid if the size was the same
as the indexed array, rather than the shape. This resulted in an index being
incorrectly allowed if the index was all False, like np.empty((3, 3))[np.full((1, 9),
False)], or incorrectly giving a ValueError about broadcasting if the index
was not all False (like np.empty((3, 3))[np.full((1, 9), True)]). Now these
examples both give IndexError with the usual error message.

Fixes #16997.

(not really sure about the NumPy commit message flags, or which one should apply to this. Is this considered an API break?)

<!-- Please be sure you are following the instructions in the dev guidelines
http://www.numpy.org/devdocs/dev/development_workflow.html
-->

<!-- We'd appreciate it if your commit message is properly formatted
http://www.numpy.org/devdocs/dev/development_workflow.html#writing-the-commit-message
-->


Problem statement:
"ValueError: operands could not be broadcast together with shapes" error on invalid boolean array index
<!-- Please describe the issue in detail here, and fill in the fields below -->

In some instances, when indexing an array with a boolean array index, instead of getting an `IndexError`, you instead get a `ValueError: operands could not be broadcast together with shapes ...` error.  It seems to happen when the product of the array and index shapes are the same (but the shapes to not match).

### Reproducing code example:

<!-- A short code example that reproduces the problem/missing feature. It should be
self-contained, i.e., possible to run as-is via 'python myproblem.py' -->

```python
>>> import numpy as np
>>> a = np.empty((3, 3))
>>> a[np.empty((1, 9), dtype=bool)]
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ValueError: operands could not be broadcast together with shapes (3,3) (1,9)
```

Compare the usual case:

```py
>>> a = np.empty((3, 4))
>>> a[np.empty((1, 9), dtype=bool)]
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
IndexError: boolean index did not match indexed array along dimension 0; dimension is 3 but corresponding boolean dimension is 1
```

Another example to back my theory that this only happens when the array sizes match

```py
>>> a = np.empty((2, 2, 4))
>>> a[np.empty((1, 4, 4), dtype=bool)]
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ValueError: operands could not be broadcast together with shapes (2,2,4) (1,4,4)
>>> a[np.empty((1, 2, 4), dtype=bool)]
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
IndexError: boolean index did not match indexed array along dimension 0; dimension is 2 but corresponding boolean dimension is 1
```

Obviously, `ValueError` is the incorrect error for an invalid index, but aside from this, an error mentioning broadcasting of the array shape and the index shape is very confusing since (unless I am mistaken), this never actually happens.

<!-- Remove these sections for a feature request -->

### Numpy/Python version information:

<!-- Output from 'import sys, numpy; print(numpy.__version__, sys.version)' -->

NumPy master:

```py
>>> import sys, numpy; print(numpy.__version__, sys.version)
1.20.0.dev0+8f60522 3.7.6 | packaged by conda-forge | (default, Mar 23 2020, 22:45:16)
[Clang 9.0.1 ]
```

Hints from issue discussion:
The reason is that the above case goes into a boolean array fast path and fails there, instead of the general advanced indexing path (which uses the `nonzero` logic first).  I agree the error could and probably should be changed to IndexError, it is currently being raised by the iterator, so would need to be checked ahead of time or replaced.  Or maybe just fall back to the slow path to give the same errors.
I guess I'm confused about boolean indexing semantics. If the index array is all False, the indexing works with no errors. In every other case, it fails:

```py
>>> a.shape
(3, 3)
>>> a[np.full((1, 9), False)]
array([], dtype=int64)
>>> a[np.array([[False, False, False, False, False, False, False, False,  True]])]
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ValueError: operands could not be broadcast together with shapes (3,3) (1,9)
>>> a[np.full((1, 100), False)]
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
IndexError: boolean index did not match indexed array along dimension 0; dimension is 3 but corresponding boolean dimension is 1
>>> a[np.nonzero(np.full((1, 9), False))]
array([], dtype=int64)
>>> a[np.nonzero(np.full((1, 100), False))]
array([], dtype=int64)
```

Is this intended behavior? It doesn't seem to be documented. I guess it's because `np.nonzero` gives empty arrays when there are no Trues, so that the shape doesn't matter? But why does it allow shapes that don't actually match (only the product matches)?

Also it doesn't work in every case:

```py
>>> np.empty((2, 4, 6))[np.full((2*2*2, 1, 6), False)]
array([], dtype=float64)
>>> np.empty((2, 4, 6))[np.full((2*2, 2, 6), False)]
array([], dtype=float64)
>>> np.empty((2, 4, 6))[np.full((2, 4*6), False)]
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
IndexError: boolean index did not match indexed array along dimension 1; dimension is 4 but corresponding boolean dimension is 24
```

It seems to work as long as the array and index have the same number of dimensions and the same size, and the index is all False.
Sounds like a buggy fast path for the all-false case. I am a bit surprised as to why we would have such a fast path, but...

Note that all the weird cases completely bypass any "nonzero" logic, since all of them occur in the fast-path for a single boolean array of identical size as the input array (which should be identical shape).

Fixing that fast path to only accept identical *shape* should automatically also switch that behaviour you just note.  I think I would be fine with just doing it, it seems strange enough and more likely to flush out serious bugs then anything else.
If you are up for it, you could try modifying this line: https://github.com/numpy/numpy/blob/e20dded7c3bf28c5767c645f074920d8017eeba9/numpy/core/src/multiarray/mapping.c#L557

So that it compares the shape exactly, which will have the effect that I described, of just using the "slow" path (which gives all the errors you would like to see).
By the way, I found both this issue and https://github.com/numpy/numpy/issues/15753 using hypothesis to test my own library. I noticed NumPy already has a couple of hypothesis tests, but not very many. I would suggest adding a lot more, especially for things that don't rely on testing finicky floating point arithmetic (like indexing semantics).
Can you make any sense of the comment above that line? It seems to say to me that this is intentional, but I'm probably misreading it. 
Sure, new tests are welcome.  The comment makes sense, but is outdated, since the paths which it intentionally lets through are now all errors, actually reversing the situation: Making it more strict ensures the same errors are raised.  But of course, as the comment states, you can also delete the shape line entirely, and instead raise the IndexError(s) you have been talking about explicitly in the boolean array case.
OK, I guess I'm missing the history of how indexing semantics used to work differently?
Changing `PyArray_SIZE` to `PyArray_SHAPE` on that line seems to fix it. I'm running the tests to see if anything else breaks. 
@asmeurer the semantics had a huge amount of weird corner cases at some point, most are fixed, but as you see some corner cases of corner cases are still around.  You need to use `PyArray_CompareLists` (or delete the line and fix up/clarify the errors in the boolean fast path explicitly, which may even be a bit nicer).

Commit message: Update the comments around boolean arrays in mapping.c

Patch to review:
diff --git a/numpy/core/src/multiarray/mapping.c b/numpy/core/src/multiarray/mapping.c
index c27e0c391ea9..91630b522ce8 100644
--- a/numpy/core/src/multiarray/mapping.c
+++ b/numpy/core/src/multiarray/mapping.c
@@ -539,22 +539,23 @@ prepare_index(PyArrayObject *self, PyObject *index,
             /*
              * There are two types of boolean indices (which are equivalent,
              * for the most part though). A single boolean index of matching
-             * dimensionality and size is a boolean index.
-             * If this is not the case, it is instead expanded into (multiple)
-             * integer array indices.
+             * shape is a boolean index. If this is not the case, it is
+             * instead expanded into (multiple) integer array indices.
              */
             PyArrayObject *nonzero_result[NPY_MAXDIMS];
 
             if ((index_ndim == 1) && allow_boolean) {
                 /*
-                 * If ndim and size match, this can be optimized as a single
-                 * boolean index. The size check is necessary only to support
-                 * old non-matching sizes by using fancy indexing instead.
-                 * The reason for that is that fancy indexing uses nonzero,
-                 * and only the result of nonzero is checked for legality.
+                 * If shapes match exactly, this can be optimized as a single
+                 * boolean index. The reason for that is that fancy indexing
+                 * uses nonzero, and only the result of nonzero is checked for
+                 * legality. Otherwise, if the shapes don't match and it isn't
+                 * the scalar case checked below, it should be an error.
                  */
                 if ((PyArray_NDIM(arr) == PyArray_NDIM(self))
-                        && PyArray_SIZE(arr) == PyArray_SIZE(self)) {
+                        && PyArray_CompareLists(PyArray_DIMS(arr),
+                                                PyArray_DIMS(self),
+                                                PyArray_NDIM(arr))) {
 
                     index_type = HAS_BOOL;
                     indices[curr_idx].type = HAS_BOOL;
diff --git a/numpy/core/tests/test_indexing.py b/numpy/core/tests/test_indexing.py
index 1069cbe8dd3a..c243be936672 100644
--- a/numpy/core/tests/test_indexing.py
+++ b/numpy/core/tests/test_indexing.py
@@ -7,8 +7,8 @@
 from numpy.core._multiarray_tests import array_indexing
 from itertools import product
 from numpy.testing import (
-    assert_, assert_equal, assert_raises, assert_array_equal, assert_warns,
-    HAS_REFCOUNT,
+    assert_, assert_equal, assert_raises, assert_raises_regex,
+    assert_array_equal, assert_warns, HAS_REFCOUNT,
     )
 
 
@@ -1239,6 +1239,33 @@ def test_boolean_indexing_weirdness(self):
         assert_raises(IndexError, lambda: a[False, [0, 1], ...])
 
 
+    def test_boolean_indexing_fast_path(self):
+        # These used to either give the wrong error, or incorrectly give no
+        # error.
+        a = np.ones((3, 3))
+        idx1 = np.array([[False]*9])
+        idx2 = np.array([[False]*10])
+        idx3 = np.array([[False]*8 + [True]])
+        assert_raises_regex(IndexError,
+            "boolean index did not match indexed array along dimension 0; "
+            "dimension is 3 but corresponding boolean dimension is 1",
+            lambda: a[idx1])
+        assert_raises_regex(IndexError,
+            "boolean index did not match indexed array along dimension 0; "
+            "dimension is 3 but corresponding boolean dimension is 1",
+            lambda: a[idx2])
+        assert_raises_regex(IndexError,
+            "boolean index did not match indexed array along dimension 0; "
+            "dimension is 3 but corresponding boolean dimension is 1",
+            lambda: a[idx3])
+        a = np.ones((1, 1, 2))
+        idx = np.array([[[True], [False]]])
+        assert_raises_regex(IndexError,
+            "boolean index did not match indexed array along dimension 1; "
+            "dimension is 1 but corresponding boolean dimension is 2",
+            lambda: a[idx])
+
+
 class TestArrayToIndexDeprecation:
     """Creating an an index from array not 0-D is an error.
 


Code review comments:
diff_hunk: @@ -539,22 +539,23 @@ prepare_index(PyArrayObject *self, PyObject *index,
             /*
              * There are two types of boolean indices (which are equivalent,
              * for the most part though). A single boolean index of matching
-             * dimensionality and size is a boolean index.
-             * If this is not the case, it is instead expanded into (multiple)
-             * integer array indices.
+             * shape is a boolean index. If this is not the case, it is
+             * instead expanded into (multiple) integer array indices.
              */
             PyArrayObject *nonzero_result[NPY_MAXDIMS];
 
             if ((index_ndim == 1) && allow_boolean) {
                 /*
-                 * If ndim and size match, this can be optimized as a single
-                 * boolean index. The size check is necessary only to support
-                 * old non-matching sizes by using fancy indexing instead.
-                 * The reason for that is that fancy indexing uses nonzero,
-                 * and only the result of nonzero is checked for legality.
+                 * If shapes match exactly, this can be optimized as a single
+                 * boolean index. The reason for that is that fancy indexing
+                 * uses nonzero, and only the result of nonzero is checked for
+                 * legality. Otherwise, if the shapes don't match and it isn't
+                 * the scalar case checked below, it should be an error.
path: numpy/core/src/multiarray/mapping.c
text: @user1:
The nonzero reference is not useful here, and I doubt its correct. Historically, it *was* correct, but that is a different issue. Maybe:

```suggestion
                 * boolean index. When the dimensions are identical but the shapes are not,
                 * this is always an error. The check ensures that these errors are raised
                 * and match those of the generic path.
```