index_reshape.Rmd

---
jupyter:
  jupytext:
    text_representation:
      extension: .Rmd
      format_name: rmarkdown
      format_version: '1.2'
      jupytext_version: 1.11.5
---

# Index ordering and reshape in NumPy and MATLAB

Let’s make a 3D array by taking a 1D array from 0 through 23, and filling the
3D array by depth (the last axis), and then by column (the second axis) and
then by row (the first axis):

```{python}
# Fill in a NumPy array
import numpy as np
numbers = np.arange(24)
py_arr = np.zeros((2, 3, 4))
n_index = 0  # Python has 0-based indices
for i in range(2):  # row index changes slowest
    for j in range(3):  # then column index
        for k in range(4):  # depth index changes fastest
            py_arr[i, j, k] = numbers[n_index];
            n_index += 1

print(py_arr[:, :, 0])
print(py_arr[:, :, 1])
```

We can do the same thing in MATLAB, or its open-source version, Octave:

<!-- We previously did this live, with
.. runblock:: octave

    >> % Fill in a MATLAB / Octave array

etc.  But then we need octave for every build, which seemed a bit much. -->
```
>> % Fill in a MATLAB / Octave array
>> numbers = 0:23;
>> m_arr = zeros(2, 3, 4);
>> n_index = 1;
>> for i = 1:2  % row index changes slowest
       for j = 1:3  % then column index
           for k = 1:4  % depth index changes fastest
               m_arr(i, j, k) = numbers(n_index);
               n_index = n_index + 1;
           end
       end
   end
>> m_arr(:, :, 1)
ans =

    0    4    8
   12   16   20

>> m_arr(:, :, 2)
ans =

    1    5    9
   13   17   21
```

Remember that MATLAB and Octave have 1-based indices.  That is, the index to
the first element on an axis is 1.  Python has 0-based indices.  Given that,
these two arrays are the same.  By “the same” I mean that that `m_arr[i, j,
k] == py_arr[i-1, j-1, k-1]` for any `i, j, k`.  You can see this from the
printout above of the first two planes in each array in Python and MATLAB.

So far, we see that NumPy / Python and MATLAB indexing are the same, apart
from the 0-based / 1-based difference.  They differ in their default way of
getting elements when doing a reshape.

When NumPy or MATLAB reshapes one array into another array, it takes the
elements from the first array in some order, and puts them into the new array
using the same order.  The default order in NumPy is to take the elements off
the last axis first, then the second to last, back to the first axis.  This is
the same order as the loop above:

```{python}
numbers = np.arange(24)
reshaped = np.reshape(numbers, (2, 3, 4))
# This is exactly the same as the original ``py_arr``
np.all(reshaped == py_arr)
```

MATLAB does a reshape using the opposite order, taking the elements off the
first axis first:

<!-- was:

.. runblock:: octave -->
```
>> m_reshaped = reshape(0:23, [2 3 4]);
>> m_reshaped(:, :, 1)
ans =

   0   2   4
   1   3   5

>> m_reshaped(:, :, 2)
ans =

    6    8   10
    7    9   11
```

If you prefer this ordering, you can ask NumPy to do the same, by using the
`order` parameter to `reshape`:

```{python}
# First axis first fetching, like MATLAB
first_1_reshaped = np.reshape(numbers, (2, 3, 4), order='F')
print(first_1_reshaped[:, :, 0])
print(first_1_reshaped[:, :, 1])
```