Multidimensional matrix permutation Julia vs. Python disagreement

I have noticed a difference in behavior between python's numpy.permute_dims and julia's Base.permutedims.

On an input 3x3x3 matrix containing elements 0:26, inclusive in both languages, they agree for the axes argument (1,2,0) but disagree for (0,2,1).

As far as I can tell from the docs, these functions should be equivalent. There's a note about permutedims being non-recursive, but I don't see why that should have behavioral consequences.

Julia is also column-major order, but again I don't see why that should effect the overall behavior.

Python code:

arr = np.array([[[0, 1, 2],
                 [3, 4, 5],
                 [6, 7, 8]],
                [[9, 10, 11],
                 [12, 13, 14],
                 [15, 16, 17]],
                [[18, 19, 20],
                 [21, 22, 23],
                 [24, 25, 26]]])
    
arr_perm = np.permute_dims(arr, axes=[0,2,1])

Output:

array([[[ 0,  3,  6],
        [ 1,  4,  7],
        [ 2,  5,  8]],

       [[ 9, 12, 15],
        [10, 13, 16],
        [11, 14, 17]],

       [[18, 21, 24],
        [19, 22, 25],
        [20, 23, 26]]])

Julia code:

arr = [
    0 1 2 
    3 4 5 
    6 7 8;;; 
    9 10 11 
    12 13 14 
    15 16 17;;;
    18 19 20 
    21 22 23 
    24 25 26
]

arr_perm = permutedims(arr, [1,3,2])

Output:

3×3×3 Array{Int64, 3}:
[:, :, 1] =
 0   9  18
 3  12  21
 6  15  24

[:, :, 2] =
 1  10  19
 4  13  22
 7  16  25

[:, :, 3] =
 2  11  20
 5  14  23
 8  17  26

It looks like your problem comes from the fact that numpy does not adhere to the "natural" mental model for indexing. In the natural mental model, if you want to address the number 20 in a 3d matrix like this

    a = np.arange(2*3*4).reshape(2, 3, 4)
    array([[[ 0,  1,  2,  3],
            [ 4,  5,  6,  7],
            [ 8,  9, 10, 11]],

           [[12, 13, 14, 15],
            [16, 17, 18, 19],
            [20, 21, 22, 23]]])

you instinctively think (counting starts from zero): «row=2, column=0, plane=1» but, in numpy, if you issue the command a[2, 0, 1], you get an error.
The correct mental model to apply with numpy is the Boxes model. The boxes model works like boxes containing boxes which in turn contain boxes and so on.

If you want address the number 20, ask yourself: «Starting from the outer box, which element in the box level 0 contains the number 20? It's the number 1», then: «Which element in box level 1 contains the number 20? it's the number 2», finally: «Which element in box level 2 contains the number 20? it's the number 0»

    a[1, 2, 0]
    20

How to address the same element when you transpose the matrix? Just reverse the sequence of indexes:

    a[1, 2, 0] == a.T[0, 2, 1]
    True

Note that when working in 2d (and only in 2d), the natural mental model and the boxes models comes to coincidence, so a[row, column] translates to

box level 0 = row

box level 1 = column

the numpy expression is a[row, column]

If you work in 3d in natural model, you have a[row, column, plane] and the numpy translation is:

box level 0 = plane

box level 1 = row

box level 2 = column

the numpy expression is a[plane, row, column]

Returning to your example, the Julia address (row, column, plane) containing the number 20 in arr is: row = 1, column = 3, plane = 3.

To transform to numpy coordinates, apply:

box level 0 = 2 # plane=3 - 1 because Python starts counting from 0

box level 1 = 0 # row=1 - 1 because Python starts counting from 0

box level 2 = 2 # column=3 - 1 because Python starts counting from 0

so the Julia indexing expression arr[1, 3, 3] becomes the numpy indexing expression arr[2, 0, 2).

Any operation you apply to Julia arrays, can be translated into numpy applying the above transformation.

Hope this helps.