##############################################
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):

.. nbplot::

    >>> # 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])
    [[  0.   4.   8.]
     [ 12.  16.  20.]]
    >>> print(py_arr[:, :, 1])
    [[  1.   5.   9.]
     [ 13.  17.  21.]]

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.

.. code-block:: matlab

    >> % 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:

.. nbplot::

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

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

.. was:

    .. runblock:: octave

.. code-block:: matlab

    >> 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``:

.. nbplot::

    >>> # First axis first fetching, like MATLAB
    >>> first_1_reshaped = numbers.reshape((2, 3, 4), order='F')
    >>> print(first_1_reshaped[:, :, 0])
    [[0 2 4]
     [1 3 5]]
    >>> print(first_1_reshaped[:, :, 1])
    [[ 6  8 10]
     [ 7  9 11]]