Repository: numpy/numpy
Pull request: BUG: Raise when unravel_index, ravel_multi_index are given empty input
Fixes #7165. Replaces #7170.

Adds tests and fixes for `unravel_index`, `ravel_multi_index` when the input is an empty list. The strategy in #7170 relied on internal use of `PyArray_DTypeFromObject`. But we can proactively use it with a minimal `int` type (I used `NPY_BOOL` as the most basic type) to determine the dtype of the object while restricting the possibilities to `int` dtypes before calling `PyArray_FromAny`.

Also added a small refactor of privately used sequence_to_arrays to put output values at the end of the paramter list and use a goto to reduce duplicate code.

Problem statement:
np.unravel_index() does not accept empty inputs
I would expect a tuple of empty outputs back. I get the following:

```
>>> np.unravel_index([], (10, 12))
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
TypeError: Iterator operand 0 dtype could not be cast from dtype('float64') to dtype('int64') according to the rule 'same_kind'
```

This does not seem to be a regression, as I can see this with numpy v1.9.2, but it does seem silly.

Hints from issue discussion:
At the very least, the error message does not really make much sense in this context.

Digging into the code in multiarray/compiled_base.c, I am trying to figure out where this error is occurring. At line 1123, `NpyIter_New()` is called on the passed in indices, and one of the flags is `NPY_ITER_ZEROSIZE_OK`, so that can't be where it is failing (and it should be having a dtype of int64).

What's a good way for me to figure out where this is failing?

Indeed, looking further down in the function, there are a bunch of checks for empty iterators, so I suspect this function was supposed to work for empty inputs, but something bad is inadvertently happening.

Here are two ways I can think of:

1) IINM, you can use `gdb` to add breakpoints to the C code.

2) Create your own "breakpoints" by throwing exceptions at strategic points in the code.  Then you can zoom into where the code is breaking by whether or not your exception or that `TypeError` is thrown first.

The failure is in the iterator, but it has nothing do, at the time of failing, with it being an empty array, but with it being a float array.

The source of the issue is that [this call](https://github.com/numpy/numpy/blob/maintenance/1.11.x/numpy/core/src/multiarray/compiled_base.c#L1107), when `indices0` is an empty sequence, defaults to a `float64` dtype, and that's what throws the iterator off. Try `np.array([]).dtype` to see the same thing at work on your Python interpreter.

There may be more elegant ways, but a check for `PyArray_SIZE(indices) == 0` that triggers a forced cast to `NPY_INTP` is about as straightforward as you can go.

I think you can just pass in the dtype into FromAny here, since there is only one reasonable dtype for indices really (Yeah, I know if someone manages to put a huge 64bit number into it on a 32bit syste, that may get funny stuff, but....).

Ah ok, I guess if you want to avoid the copy, ...

I think I can put together a PR that does that check and forced cast. I
thought it was already getting forced cast, so I originally didn't think it
was the source of the problem, but I can see now what the problem is.

On Tue, Feb 2, 2016 at 10:46 AM, seberg notifications@github.com wrote:

> Ah ok, I guess if you want to avoid the copy, ...
> 
> —
> Reply to this email directly or view it on GitHub
> https://github.com/numpy/numpy/issues/7165#issuecomment-178647541.

May not be worth it, but that's one of the sweetnesses of the new iterator.

What seems ludicrous is that we have specialized branches and functions for C and Fortran order... There's a PR with a net negative line count of ~70 lines waiting there...

@seberg: the nice thing is that we even have a special type for indices
(intp), so if someone wants a 64 bit index on a 32 bit system we can just
point to the name of the type :-)
On Feb 2, 2016 8:00 AM, "Jaime" notifications@github.com wrote:

> May not be worth it, but that's one of the sweetnesses of the new iterator.
> 
> What seems ludicrous is that we have specialized branches and functions
> for C and Fortran order... There's a PR with a net negative line count of
> ~70 lines waiting there...
> 
> —
> Reply to this email directly or view it on GitHub
> https://github.com/numpy/numpy/issues/7165#issuecomment-178653513.

Commit message: BUG: fixes from review

Patch to review:
diff --git a/numpy/core/src/multiarray/compiled_base.c b/numpy/core/src/multiarray/compiled_base.c
index 1e0c1eb5ff8b..aebab4e43f62 100644
--- a/numpy/core/src/multiarray/compiled_base.c
+++ b/numpy/core/src/multiarray/compiled_base.c
@@ -832,17 +832,62 @@ arr_interp_complex(PyObject *NPY_UNUSED(self), PyObject *args, PyObject *kwdict)
     return NULL;
 }
 
+static const char *EMPTY_SEQUENCE_ERR_MSG = "indices must be integral: the provided " \
+    "empty sequence was inferred as float. Wrap it with " \
+    "'np.array(indices, dtype=np.intp)'";
+
+static const char *NON_INTEGRAL_ERROR_MSG = "only int indices permitted";
+
+/* Convert obj to an ndarray with integer dtype or fail */
+static PyArrayObject *
+astype_anyint(PyObject *obj) {
+    PyArrayObject *ret;
+
+    if (!PyArray_Check(obj)) {
+        /* prefer int dtype */
+        PyArray_Descr *dtype_guess = NULL;
+        if (PyArray_DTypeFromObject(obj, NPY_MAXDIMS, &dtype_guess) < 0) {
+            return NULL;
+        }
+        if (dtype_guess == NULL) {
+            if (PySequence_Check(obj) && PySequence_Size(obj) == 0) {
+                PyErr_SetString(PyExc_TypeError, EMPTY_SEQUENCE_ERR_MSG);
+            }
+            return NULL;
+        }
+        ret = (PyArrayObject*)PyArray_FromAny(obj, dtype_guess, 0, 0, 0, NULL);
+        if (ret == NULL) {
+            return NULL;
+        }
+    }
+    else {
+        ret = (PyArrayObject *)obj;
+        Py_INCREF(ret);
+    }
+
+    if (!(PyArray_ISINTEGER(ret) || PyArray_ISBOOL(ret))) {
+        /* ensure dtype is int-based */
+        PyErr_SetString(PyExc_TypeError, NON_INTEGRAL_ERROR_MSG);
+        return NULL;
+    }
+
+    return ret;
+}
+
 /*
  * Converts a Python sequence into 'count' PyArrayObjects
  *
- * seq       - Input Python object, usually a tuple but any sequence works.
- * op        - Where the arrays are placed.
- * count     - How many arrays there should be (errors if it doesn't match).
- * paramname - The name of the parameter that produced 'seq'.
+ * seq         - Input Python object, usually a tuple but any sequence works.
+ *               Must have integral content.
+ * paramname   - The name of the parameter that produced 'seq'.
+ * count       - How many arrays there should be (errors if it doesn't match).
+ * op          - Where the arrays are placed.
  */
-static int sequence_to_arrays(PyObject *seq,
-                                PyArrayObject **op, int count,
-                                char *paramname)
+static int int_sequence_to_arrays(PyObject *seq,
+                              char *paramname,
+                              int count,
+                              PyArrayObject **op
+                              )
 {
     int i;
 
@@ -854,29 +899,26 @@ static int sequence_to_arrays(PyObject *seq,
     }
 
     for (i = 0; i < count; ++i) {
+        PyArray_Descr *dtype = NULL;
         PyObject *item = PySequence_GetItem(seq, i);
         if (item == NULL) {
-            while (--i >= 0) {
-                Py_DECREF(op[i]);
-                op[i] = NULL;
-            }
-            return -1;
+            goto fail;
         }
-
-        op[i] = (PyArrayObject *)PyArray_FROM_O(item);
+        op[i] = astype_anyint(item);
+        Py_DECREF(item);
         if (op[i] == NULL) {
-            while (--i >= 0) {
-                Py_DECREF(op[i]);
-                op[i] = NULL;
-            }
-            Py_DECREF(item);
-            return -1;
+            goto fail;
         }
-
-        Py_DECREF(item);
     }
 
     return 0;
+
+fail:
+    while (--i >= 0) {
+        Py_XDECREF(op[i]);
+        op[i] = NULL;
+    }
+    return -1;
 }
 
 /* Inner loop for unravel_index */
@@ -1023,11 +1065,10 @@ arr_ravel_multi_index(PyObject *self, PyObject *args, PyObject *kwds)
     }
 
     /* Get the multi_index into op */
-    if (sequence_to_arrays(coords0, op, dimensions.len, "multi_index") < 0) {
+    if (int_sequence_to_arrays(coords0, "multi_index", dimensions.len, op) < 0) {
         goto fail;
     }
 
-
     for (i = 0; i < dimensions.len; ++i) {
         op_flags[i] = NPY_ITER_READONLY|
                       NPY_ITER_ALIGNED;
@@ -1222,15 +1263,9 @@ arr_unravel_index(PyObject *self, PyObject *args, PyObject *kwds)
 
     unravel_size = PyArray_MultiplyList(dimensions.ptr, dimensions.len);
 
-    if (!PyArray_Check(indices0)) {
-        indices = (PyArrayObject*)PyArray_FROM_O(indices0);
-        if (indices == NULL) {
-            goto fail;
-        }
-    }
-    else {
-        indices = (PyArrayObject *)indices0;
-        Py_INCREF(indices);
+    indices = astype_anyint(indices0);
+    if (indices == NULL) {
+        goto fail;
     }
 
     dtype = PyArray_DescrFromType(NPY_INTP);
diff --git a/numpy/lib/tests/test_index_tricks.py b/numpy/lib/tests/test_index_tricks.py
index 3246f68ff41d..028bba37d2e4 100644
--- a/numpy/lib/tests/test_index_tricks.py
+++ b/numpy/lib/tests/test_index_tricks.py
@@ -77,6 +77,26 @@ def test_basic(self):
             [[3, 6, 6], [4, 5, 1]])
         assert_equal(np.unravel_index(1621, (6, 7, 8, 9)), [3, 1, 4, 1])
 
+    def test_empty_indices(self):
+        msg1 = 'indices must be integral: the provided empty sequence was'
+        msg2 = 'only int indices permitted'
+        assert_raises_regex(TypeError, msg1, np.unravel_index, [], (10, 3, 5))
+        assert_raises_regex(TypeError, msg1, np.unravel_index, (), (10, 3, 5))
+        assert_raises_regex(TypeError, msg2, np.unravel_index, np.array([]),
+                            (10, 3, 5))
+        assert_equal(np.unravel_index(np.array([],dtype=int), (10, 3, 5)),
+                     [[], [], []])
+        assert_raises_regex(TypeError, msg1, np.ravel_multi_index, ([], []),
+                            (10, 3))
+        assert_raises_regex(TypeError, msg1, np.ravel_multi_index, ([], ['abc']),
+                            (10, 3))
+        assert_raises_regex(TypeError, msg2, np.ravel_multi_index,
+                    (np.array([]), np.array([])), (5, 3))
+        assert_equal(np.ravel_multi_index(
+                (np.array([], dtype=int), np.array([], dtype=int)), (5, 3)), [])
+        assert_equal(np.ravel_multi_index(np.array([[], []], dtype=int),
+                     (5, 3)), [])
+
     def test_big_indices(self):
         # ravel_multi_index for big indices (issue #7546)
         if np.intp == np.int64:


Code review comments:
diff_hunk: @@ -854,29 +899,26 @@ static int sequence_to_arrays(PyObject *seq,
     }
 
     for (i = 0; i < count; ++i) {
+        PyArray_Descr *dtype = NULL;
path: numpy/core/src/multiarray/compiled_base.c
text: @user1:
    numpy/core/src/multiarray/compiled_base.c:902:24: warning: unused variable ‘dtype’

@author:
fixed