Repository: numpy/numpy
Pull request: BUG: fix mgrid output for lower precision float inputs
Floats besides float64 were being coerced to integers in np.mgrid and complex step sizes for the index trick classes would fail for complex64 input. Fixes #16466

Before:
```python
>>> np.mgrid[np.float32(0.1):np.float32(0.33):np.float32(0.1),]
array([[0, 0, 0]])
>>> np.r_[0:10:np.complex64(3j)]
array([], dtype=complex128)
```
After:
```python
>>> np.mgrid[np.float32(0.1):np.float32(0.33):np.float32(0.1),]
array([[0.1, 0.2, 0.3]])
>>> np.r_[0:10:np.complex64(3j)]
array([ 0.,  5., 10.])
```

The input types do not affect the output dtype, which are always int64 or float64, but they can affect the output in other ways due to `arange` precision issues (#5808).
```python
>>> np.mgrid[np.float32(0.1):np.float32(0.3):np.float32(0.1)]
array([0.1, 0.2, 0.3])
>>> np.mgrid[np.float64(0.1):np.float64(0.3):np.float64(0.1)]
array([0.1, 0.2])
```

@seberg [mentioned](https://github.com/numpy/numpy/issues/16466#issuecomment-636928083) that the `mgrid` code could use a refactor. A lot of its functionality is duplicated in `AxisConcatenator` and `meshgrid`. I think the code could almost be replaced by `np.meshgrid(*[np.r_[k.start:k.stop:k.step] for k in key], sparse=self.sparse)`. I could verify this and clean it up if it sounds reasonable.

Problem statement:
BUG: mgrid result sensitive to data type
Using `np.mgrid' depends on slice dtype

### Reproducing code example:
```python
import numpy as np

a = np.mgrid[[slice(np.float32(0.1), np.float32(0.3), np.float32(0.1)) for _ in range(1)]]
# returns array([[0, 0, 0]])

b = np.mgrid[[slice(np.float64(0.1), np.float64(0.3), np.float64(0.1)) for _ in range(1)]]
# returns array([[0.1, 0.2]])

np.testing.assert_allclose(a, b)
```

### Numpy/Python version information:

<!-- Output from 'import sys, numpy; print(numpy.__version__, sys.version)' -->
1.18.4 3.6.7 (default, Mar 27 2020, 15:56:09) 
[GCC 6.5.0 20181026]

Hints from issue discussion:
Using tuple syntax is much simpler here:
```
np.mgrid[np.float32(0.1):np.float32(0.3):np.float32(0.1),]   # note the last comma
```
In any case, you are right. `mgrid` seems a bit funny, and could use some cleanup and refactor also. In particular, I think it should just call `np.arange` (or `np.linspace` for the complex step case) and then reshape/create a dense grid if necessary.

Maybe the float32 bug is so obvious, that we can simply make it work (by going through arange unfortunately, I guess)?
The problem seems to be that `isinstance(start, float)` returns False when `start` is a `np.float32`, so the output is initialized to use the default `int` instead. Oddly, `np.ogrid` avoids this problem because it outputs a list instead of an array, so the values just get overwritten, instead of coerced into `int` like `np.mgrid`.

I can fix this.

Commit message: BUG: fix mgrid output for lower precision float inputs

Floats besides float64 were being coerced to integers
and complex step sizes for the index trick classes
would fail for complex64 input. Fixes #16466

Patch to review:
diff --git a/numpy/lib/index_tricks.py b/numpy/lib/index_tricks.py
index e86c7a7bbc6a..8a73e35eddb2 100644
--- a/numpy/lib/index_tricks.py
+++ b/numpy/lib/index_tricks.py
@@ -154,15 +154,15 @@ def __getitem__(self, key):
                     start = 0
                 if step is None:
                     step = 1
-                if isinstance(step, complex):
+                if isinstance(step, (_nx.complexfloating, complex)):
                     size.append(int(abs(step)))
                     typ = float
                 else:
                     size.append(
                         int(math.ceil((key[k].stop - start)/(step*1.0))))
-                if (isinstance(step, float) or
-                        isinstance(start, float) or
-                        isinstance(key[k].stop, float)):
+                if (isinstance(step, (_nx.floating, float)) or
+                        isinstance(start, (_nx.floating, float)) or
+                        isinstance(key[k].stop, (_nx.floating, float))):
                     typ = float
             if self.sparse:
                 nn = [_nx.arange(_x, dtype=_t)
@@ -176,7 +176,7 @@ def __getitem__(self, key):
                     start = 0
                 if step is None:
                     step = 1
-                if isinstance(step, complex):
+                if isinstance(step, (_nx.complexfloating, complex)):
                     step = int(abs(step))
                     if step != 1:
                         step = (key[k].stop - start)/float(step-1)
@@ -194,7 +194,7 @@ def __getitem__(self, key):
             start = key.start
             if start is None:
                 start = 0
-            if isinstance(step, complex):
+            if isinstance(step, (_nx.complexfloating, complex)):
                 step = abs(step)
                 length = int(step)
                 if step != 1:
@@ -344,7 +344,7 @@ def __getitem__(self, key):
                     start = 0
                 if step is None:
                     step = 1
-                if isinstance(step, complex):
+                if isinstance(step, (_nx.complexfloating, complex)):
                     size = int(abs(step))
                     newobj = linspace(start, stop, num=size)
                 else:
diff --git a/numpy/lib/tests/test_index_tricks.py b/numpy/lib/tests/test_index_tricks.py
index 905165a99b89..b71356219ebc 100644
--- a/numpy/lib/tests/test_index_tricks.py
+++ b/numpy/lib/tests/test_index_tricks.py
@@ -249,6 +249,26 @@ def test_mgrid_size_none_handling(self, start, stop, step, expected):
         assert_equal(grid.size, expected[0])
         assert_equal(grid_small.size, expected[1])
 
+    def test_accepts_npfloating(self):
+        # regression test for #16466
+        grid64 = mgrid[0.1:0.33:0.1, ]
+        grid32 = mgrid[np.float32(0.1):np.float32(0.33):np.float32(0.1), ]
+        assert_(not isinstance(grid32.dtype, int))
+        assert_array_almost_equal(grid64, grid32)
+
+        # different code path for single slice
+        grid64 = mgrid[0.1:0.33:0.1]
+        grid32 = mgrid[np.float32(0.1):np.float32(0.33):np.float32(0.1)]
+        assert_(not isinstance(grid32.dtype, int))
+        assert_array_almost_equal(grid64, grid32)
+
+    def test_accepts_npcomplexfloating(self):
+        grid = mgrid[0.1:0.3:np.complex64(3j), ]
+        assert_array_almost_equal(grid, np.array([[0.1, 0.2, 0.3]]))
+
+        # different code path for single slice
+        grid = mgrid[0.1:0.3:np.complex64(3j)]
+        assert_array_almost_equal(grid, np.array([0.1, 0.2, 0.3]))
 
 class TestConcatenator:
     def test_1d(self):
@@ -270,6 +290,9 @@ def test_complex_step(self):
         g = r_[0:36:100j]
         assert_(g.shape == (100,))
 
+        g = r_[0:36:np.complex64(100j)]
+        assert_(g.shape == (100,))
+
     def test_2d(self):
         b = np.random.rand(5, 5)
         c = np.random.rand(5, 5)


Code review comments:
diff_hunk: @@ -249,6 +249,26 @@ def test_mgrid_size_none_handling(self, start, stop, step, expected):
         assert_equal(grid.size, expected[0])
         assert_equal(grid_small.size, expected[1])
 
+    def test_accepts_npfloating(self):
+        # regression test for #16466
+        grid64 = mgrid[0.1:0.33:0.1, ]
+        grid32 = mgrid[np.float32(0.1):np.float32(0.33):np.float32(0.1), ]
+        assert_(not isinstance(grid32.dtype, int))
path: numpy/lib/tests/test_index_tricks.py
text: @user1:
This test is impossible to make fail, `isinstance` isn't useful on dtypes

@user1:

```suggestion
        assert grid32.dtype == np.float32
```

@author:
whoops, it should be `np.float64` though as `mgrid` only outputs `float64` or `int64` currently.