7
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?

user2133814
  • 2,431
  • 1
  • 24
  • 34

1 Answers1

6

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]]])
Community
  • 1
  • 1
hpaulj
  • 221,503
  • 14
  • 230
  • 353
  • I still don't understand it, but I appreciate the answer and that its hard to explain, so I'll accept it. Even harder to understand is why this is the desired behavior, independent of how the code works internally. The docs you linked to aren't very complete IMHO. My lesson is to just not combine basic and advanced indexing. – user2133814 Jan 26 '16 at 22:04
  • I've added examples of replacing the slice with indexing. – hpaulj Jan 27 '16 at 02:14
  • 1
    To be clear, this is not desired behavior, although there are good (historical) reasons for it. We (numpy devs) have preliminary plans to deprecate it in a future release, and provide easier to understand functionality that works as you might expect. See https://github.com/numpy/numpy/pull/6256 for more details. – shoyer Jan 27 '16 at 05:11
  • I was struggling to justify the behavior in this case with a scalar 1st index. It's not as ambiguous as more complex mixed indexing examples. But changing something as fundamental as indexing is fraught with backward compatibility issues. – hpaulj Jan 27 '16 at 07:26