Why does numpy have this quirk when indexing with slices vs. lists?-CodePudding

One would think that indexing with slices and equivalent lists is equivalent in the result, and it mostly is:

>>> b = np.array([[0,1,2],[3,4,5]])
>>> b[0:2,0:2] # slice & slice
array([[0, 1],
    [3, 4]])
>>> b[0:2,[0,1]] # slice & list
array([[0, 1],
    [3, 4]])
>>> b[[0,1],0:2] # list & slice
array([[0, 1],
    [3, 4]])

However:

>>> b[[0,1],[0,1]] # list & list
array([0, 4])

Okay, there is a work-around:

>>> b[[0,1],:][:,[0,1]]
array([[0, 1],
    [3, 4]])

But why is this necessary?

(Note: I originally refered to slices as "ranges".)

CodePudding user response：

One would think that indexing with ranges and equivalent lists is equivalent in the result

One would be wrong. When you say "ranges" you actually mean "slices". In any case, indexing with a slice is considered "basic indexing", and indexing with a list (or array) is "advanced indexing".

The rules to how these two mix together and the shape of the resulting array are pretty arcane. However, your examples are the simplest case, a single advanced index and a basic index (the slice), from the docs:

In the simplest case, there is only a single advanced index. A single advanced index can for example replace a slice and the result array will be the same, however, it is a copy and may have a different memory layout. A slice is preferable when it is possible.

So, you'll notice, advanced indexing will always create a copy, if you only use basic indexing, you get a view. This is an important distinction (1).

But that explains why the slicing seem equivalent in that specific case. However, in the case of:

>>> b = np.array([[0,1,2],[3,4,5]])
>>> b[[0,1],[0,1]]
array([0, 4])

This is just purely integer array advanced indexing

When the index consists of as many integer arrays as the array being indexed has dimensions, the indexing is straight forward, but different from slicing.

Advanced indexes always are broadcast and iterated as one result[i_1, ..., i_M] == x[ind_1[i_1, ..., i_M], ind_2[i_1, ..., i_M], ..., ind_N[i_1, ..., i_M]]

(1) The difference between views and copies, basic vs advanced indxeing

>>> b = np.array([[0,1,2],[3,4,5]])
>>> view = b[0:2,0:2]
>>> copy = b[0:2,[0,1]]
>>> view
array([[0, 1],
       [3, 4]])
>>> copy
array([[0, 1],
       [3, 4]])
>>> view[0,0] = 1337
>>> b
array([[1337,    1,    2],
       [   3,    4,    5]])
>>> view
array([[1337,    1],
       [   3,    4]])
>>> copy
array([[0, 1],
       [3, 4]])

CodePudding user response：

As an addenda to Juanpa's answer, https://stackoverflow.com/a/70212077/901925 you can use two arrays to get the block - but they have to broadcast against each other.

In [348]: b = np.array([[0,1,2],[3,4,5]])
In [349]: b[[[0],[1]], [0,1]]
Out[349]: 
array([[0, 1],
       [3, 4]])

np.ix_ is a convenience function for creating this kind of indexing set:

In [350]: np.ix_([0,1],[0,1])
Out[350]: 
(array([[0],
        [1]]),
 array([[0, 1]]))
In [351]: b[np.ix_([0,1],[0,1])]
Out[351]: 
array([[0, 1],
       [3, 4]])

The first index array is (n,1) shape, the second (1,m); together they define a (n,m) space.

You specified a (n,) and (n,) pair, which broadcast to (n,). Think of this as the 'diagonal' of [351].

In some languages like MATLAB pairs of indexing arrays like this do specify the block. But to get the diagonal, you have use an extra step, something like sub2ind to convert the pairs to a 1d indexing array. In numpy the block index requires the extra ix_ like step, but logically it's all explainable by broadcasting.