Repository: numba/numba
Pull request: Fix array iterator layout.
As title. Fixes #5692

<!--

Thanks for wanting to contribute to Numba :)

First, if you need some help or want to chat to the core developers, please
visit https://gitter.im/numba/numba for real time chat or post to the Numba
mailing list https://groups.google.com/a/continuum.io/forum/#!forum/numba-users.

Here's some guidelines to help the review process go smoothly.

0. Please write a description in this text box of the changes that are being
   made.

1. Please ensure that you have written units tests for the changes made/features
   added.

2. If you are closing an issue please use one of the automatic closing words as
   noted here: https://help.github.com/articles/closing-issues-using-keywords/

3. If your pull request is not ready for review but you want to make use of the
   continuous integration testing facilities please label it with `[WIP]` and
   then remove the label when you'd like it to be reviewed.

4. Once review has taken place please do not add features or make changes out of
   the scope of those requested by the reviewer (doing this just add delays as
   already reviewed code ends up having to be re-reviewed/it is hard to tell
   what is new etc!). Further, please do not rebase your branch on master/force
   push/rewrite history, doing any of these causes the context of any comments
   made by reviewers to be lost. If conflicts occur against master they should
   be resolved by merging master into the branch used for making the pull
   request.

Many thanks in advance for your cooperation!

-->


Problem statement:
issue with >1d array iteration
numba 0.49.0
numpy 1.18.2
Python 3.8.2 (default, May  7 2020, 20:00:49) [GCC 7.3.0] on linux

```python
import numba as nb
import numpy as np


@nb.njit()
def outer_product_sum(x: np.ndarray, y: np.ndarray) -> np.ndarray:
    assert x.ndim == 2, "dimensions of both arrays must be 2"
    assert y.ndim == 2, "dimensions of both arrays must be 2"
    assert x.shape[1] == y.shape[1], "The second dimension of both arrays must be equal"

    o = np.zeros((x.shape[0], y.shape[0]))
    for v1, v2 in zip(x.T, y.T):
        o += np.outer(v1, v2)
    return o

x = np.array([[1, 8, 9, 1], [4, 64, 27, 1]])
y = np.array([[1, 1/2, 1/3, 1],  [1/2, 1/4, 1/9, 1], [1/2, 1/4, 1/9, 1]])
expected = np.outer(x[:,0], y[:,0]) + np.outer(x[:,1], y[:,1]) + np.outer(x[:,2], y[:,2]) + np.outer(x[:,3], y[:,3])
actual = outer_product_sum(x, y)
assert np.isclose(expected, actual).all()
```

Without the decorator, this test passes, but with the decorator, it fails. On inspection, it is clear that the issue is not with the sum, but with the inner product. In this context, `np.outer` is not producing correct results. `v1`, `v2` are as expected, but `np.outer` no longer performs the outer product.

Hints from issue discussion:
Thanks for the report, I can reproduce. I don't think it's a problem with `outer`, but instead a problem with iteration over >1d arrays (which, IIRC, was added in 0.49). For example, this seems to work as expected:
```python
import numba as nb
import numpy as np


@nb.njit()
def outer_product_sum(x: np.ndarray, y: np.ndarray) -> np.ndarray:
    assert x.ndim == 2, "dimensions of both arrays must be 2"
    assert y.ndim == 2, "dimensions of both arrays must be 2"
    assert x.shape[1] == y.shape[1], "The second dimension of both arrays must be equal"

    o = np.zeros((x.shape[0], y.shape[0]))
    for i in range(x.shape[0]):
        v1 = x[:, i]
        v2 = y[:, i]
        o += np.outer(v1, v2)
    return o

x = np.array([[1, 8, 9, 1], [4, 64, 27, 1]])
y = np.array([[1, 1/2, 1/3, 1],  [1/2, 1/4, 1/9, 1], [1/2, 1/4, 1/9, 1]])

actual = outer_product_sum(x, y)
expected = outer_product_sum.py_func(x, y)
assert np.isclose(expected, actual).all()
```
Thanks, that fixed the issue for me too. Looks like it was an issue with iteration.
Narrowed it down to:
```python
import numba as nb
import numpy as np

x = np.array([[1, 2,], [3, 4]])

@nb.njit
def foo(x):
    z = zip(x.T, x.T)
    y1, _ = next(z)
    print(y1.ravel())

foo(x)
foo.py_func(x)
```
probably iteration, F order and `ravel()` related.

Commit message: Add test for strided

Patch to review:
diff --git a/numba/core/types/iterators.py b/numba/core/types/iterators.py
index c66a706fdab..8c86c3a2a85 100644
--- a/numba/core/types/iterators.py
+++ b/numba/core/types/iterators.py
@@ -103,5 +103,6 @@ def __init__(self, array_type):
         elif nd == 1:
             yield_type = array_type.dtype
         else:
-            yield_type = array_type.copy(ndim=array_type.ndim - 1)
+            # iteration semantics leads to C order layout
+            yield_type = array_type.copy(ndim=array_type.ndim - 1, layout='C')
         super(ArrayIterator, self).__init__(name, yield_type)
diff --git a/numba/tests/test_array_iterators.py b/numba/tests/test_array_iterators.py
index 59303775b86..daf91ef1feb 100644
--- a/numba/tests/test_array_iterators.py
+++ b/numba/tests/test_array_iterators.py
@@ -174,7 +174,7 @@ def check_array_ndenumerate_sum(self, arr, arrty):
         self.check_array_unary(arr, arrty, array_ndenumerate_sum)
 
     def test_array_iter(self):
-        # Test iterating over a 1d array
+        # Test iterating over arrays
         arr = np.arange(6)
         self.check_array_iter_1d(arr)
         self.check_array_iter_items(arr)
@@ -190,6 +190,40 @@ def test_array_iter(self):
         self.check_array_iter_items(arr)
         self.check_array_iter_items(arr.T)
 
+    def test_array_iter_yielded_order(self):
+        # See issue #5692
+        @jit(nopython=True)
+        def foo(arr):
+            t = []
+            for y1 in arr:
+                for y2 in y1:
+                    t.append(y2.ravel())
+            return t
+
+        # 'F' ordered
+        arr = np.arange(24).reshape((2, 3, 4), order='F')
+        expected = foo.py_func(arr)
+        got = foo(arr)
+        np.testing.assert_allclose(expected, got)
+
+        # 'A' ordered, outer strided
+        arr = np.arange(64).reshape((4, 8, 2), order='F')[::2, :, :]
+        expected = foo.py_func(arr)
+        got = foo(arr)
+        np.testing.assert_allclose(expected, got)
+
+        # 'A' ordered, middle strided
+        arr = np.arange(64).reshape((4, 8, 2), order='F')[:, ::2, :]
+        expected = foo.py_func(arr)
+        got = foo(arr)
+        np.testing.assert_allclose(expected, got)
+
+        # 'A' ordered, inner strided
+        arr = np.arange(64).reshape((4, 8, 2), order='F')[:, :, ::2]
+        expected = foo.py_func(arr)
+        got = foo(arr)
+        np.testing.assert_allclose(expected, got)
+
     def test_array_view_iter(self):
         # Test iterating over a 1d view over a 2d array
         arr = np.arange(12).reshape((3, 4))


Code review comments:
diff_hunk: @@ -103,5 +103,6 @@ def __init__(self, array_type):
         elif nd == 1:
             yield_type = array_type.dtype
         else:
-            yield_type = array_type.copy(ndim=array_type.ndim - 1)
+            # iteration semantics leads to C order layout
+            yield_type = array_type.copy(ndim=array_type.ndim - 1, layout='C')
path: numba/core/types/iterators.py
text: @user1:
setting layout to `'A'` is probably the safest fix.

@author:
Done, also added the example as a test.