Numpy 3D array transposed when indexed in single step vs two steps

Question

import numpy as np
x = np.random.randn(2, 3, 4)
mask = np.array([1, 0, 1, 0], dtype=np.bool)
y = x[0, :, mask]
z = x[0, :, :][:, mask]
print(y)
print(z)
print(y.T)

Why does doing the above operation in two steps result in the transpose of doing it in one step?

Answer 1

Here's the same behavior with a list index:

In [87]: x=np.arange(2*3*4).reshape(2,3,4)
In [88]: x[0,:,[0,2]]
Out[88]: 
array([[ 0,  4,  8],
       [ 2,  6, 10]])
In [89]: x[0,:,:][:,[0,2]]
Out[89]: 
array([[ 0,  2],
       [ 4,  6],
       [ 8, 10]])

In the 2nd case, x[0,:,:] returns a (3,4) array, and the next index picks 2 columns.

In the 1st case, it first selects on the first and last dimensions, and appends the slice (the middle dimension). The 0 and [0,2] produce a 2 dimension, and the 3 from the middle is appended, giving (2,3) shape.

This is a case of mixed basic and advanced indexing.

http://docs.scipy.org/doc/numpy/reference/arrays.indexing.html#combining-advanced-and-basic-indexing

In the first case, the dimensions resulting from the advanced indexing operation come first in the result array, and the subspace dimensions after that.

This is not an easy case to comprehend or explain. Basically there's some ambiguity as to what the final dimension should be. It tries to illustrate with an example x[:,ind_1,:,ind_2] where ind_1 and ind_2 are 3d (or together broadcast to that).

Earlier attempts to explain this are:

How does numpy order array slice indices?

Combining slicing and broadcasted indexing for multi-dimensional numpy arrays

===========================

A way around this problem is to replace the slice with an array - a column vector

In [221]: x[0,np.array([0,1,2])[:,None],[0,2]]
Out[221]: 
array([[ 0,  2],
       [ 4,  6],
       [ 8, 10]])
In [222]: np.ix_([0],[0,1,2],[0,2])
Out[222]: 
(array([[[0]]]), array([[[0],
         [1],
         [2]]]), array([[[0, 2]]]))
In [223]: x[np.ix_([0],[0,1,2],[0,2])]
Out[223]: 
array([[[ 0,  2],
        [ 4,  6],
        [ 8, 10]]])

Though this last case is 3d, (1,3,2). ix_ didn't like the scalar 0. An alternate way of using ix_:

In [224]: i,j=np.ix_([0,1,2],[0,2])
In [225]: x[0,i,j]
Out[225]: 
array([[ 0,  2],
       [ 4,  6],
       [ 8, 10]])

And here's a way of getting the same numbers, but in a (2,1,3) array:

In [232]: i,j=np.ix_([0,2],[0])
In [233]: x[j,:,i]
Out[233]: 
array([[[ 0,  4,  8]],

       [[ 2,  6, 10]]])