[Numpy-discussion] Indexing bug

[Nathaniel Smith]

On Sun, Mar 31, 2013 at 6:14 AM, Ivan Oseledets
<ivan.oseledets at gmail.com> wrote:
>> I am using numpy 1.6.1,
>> and encountered a wierd fancy indexing bug:
>>
>> import numpy as np
>> c = np.random.randn(10,200,10);
>>
>> In [29]: print c[[0,1],:200,:2].shape
>> (2, 200, 2)
>>
>> In [30]: print c[[0,1],:200,[0,1]].shape
>> (2, 200)
>>
>> It means, that here fancy indexing is not working right for a 3d array.
>>
>
> On Sat, Mar 30, 2013 at 11:01 AM, Ivan Oseledets
> <ivan.oseledets at gmail.com>wrote:
>
>> I am using numpy 1.6.1,
>> and encountered a wierd fancy indexing bug:
>>
>> import numpy as np
>> c = np.random.randn(10,200,10);
>>
>> In [29]: print c[[0,1],:200,:2].shape
>> (2, 200, 2)
>>
>> In [30]: print c[[0,1],:200,[0,1]].shape
>> (2, 200)
>>
>> It means, that here fancy indexing is not working right for a 3d array.
>>
> -->
> It is working fine, review the docs:
>
> http://docs.scipy.org/doc/numpy/reference/arrays.indexing.html#advanced-indexing
>
> In your return, item [0, :] is c[0, :, 0] and item[1, :]is c[1, :, 1].
>
> If you want a return of shape (2, 200, 2) where item [i, :, j] is c[i, :,
> j] you could use slicing:
>
>  c[:2, :200, :2]
>
> or something more elaborate like:
>
> c[np.arange(2)[:, None, None], np.arange(200)[:, None], np.arange(2)]
>
> Jaime
> --->
>
> Oh!  So it is not a bug, it is a feature, which is completely
> incompatible with other array based languages (MATLAB and Fortran). To
> me, I can not find a single explanation why it is so in numpy.
> Taking submatrices from a matrix is a common operation and the syntax
> above is very natural to take submatrices, not a weird diagonal stuff.
> i.e.,
>
> c = np.random.randn(100,100)
> d = c[[0,3],[2,3]]
>
> should NOT produce two numbers! (and you can not do it using slices!)
>
> In MATLAB and Fortran
> c(indi,indj)
> will produce a 2 x 2 matrix.
> How it can be done in numpy (and why the complications?)>
> So, please consider this message as a feature request.

Numpy's handling of such things is strictly more general than MATLAB
and Fortran's (AFAIK), and fits in with the rest of the system quite
nicely, I'd say.

The logic is: if you do
  a[row_coords, col_coords]
then it treats row_coords and col_coords as parallel arrays, where
each corresponding pair of entries in the two arrays gives the
coordinates of one entry in the result -- so you get [a[row_coords[0],
col_coords[0]], a[row_coords[1], col_coords[1]], ...].

This follows numpy's usual rules for arrays that are supposed to
align: just like if you did 'row_coords + col_coords', say, which
would give you [row_coords[0] + col_coords[0], row_coords[1] +
col_coords[1], ...].

AND, all of this works just the same if row_coords and col_coords are
arbitrary-dimensional arrays: the output of indexing (just like the
output of '+') will have the same dimensionality as its input. So if
they're both 2x2 arrays, then a[row_coords, col_coords] will be
  [[a[row_coords[0, 0], col_coords[0, 0]], a[row_coords[0, 1],
col_coords[0, 1]],
   [a[row_coords[1, 0], col_coords[1, 0]], a[row_coords[1, 1],
col_coords[1, 1]],
  ]

AND (here's the solution to your problem), this "aligning" uses the
same rules as "+" does, i.e., broadcasting applies:
 http://docs.scipy.org/doc/numpy/user/basics.broadcasting.html

So for your example, write
  c[[[0], [1]], :200, [[0, 1]]]
which by broadcasting is equivalent to
  c[[[0, 1], [0, 1]],
    :200,
    [[0, 0], [1, 1]]]
which is what you want.

The only problem then is that slicing indexing happens "inside" fancy
indexing ("fancy indexing" is the name for this
line-up-the-arrays-and-extract-coordinates thing). So what this means
is that when you use both inside a single indexing operation, the what
it does is for each set of corresponding fancy indexing coordinates,
it doesn't just extract a single element for the final array -- it
extracts an entire sub-array. The practical result is that the shape
of your output array will be the shape of your aligned fancy indexing
arrays (after broadcasting), and then with any sliced axes stuck on
the end. So you might expect that the above expression would give you
an array of shape
  (2, 200, 2)
but in fact it will be of shape
  (2, 2, 200)
because the (2, 2) part is the shape of the fancy index arrays, and
the (200,) is the shape of the slices. (Notice that there is no
relation at all between the shape of the input array and the shape of
the fancy indexes!) You can fix this up with a call to rollaxes(), or
by getting even fancier, using fancy indexes for all three axes, and
making sure they each broadcast to the shape (2, 200, 2):

coord1 = np.asarray([0, 1]).reshape((2, 1, 1))
coord2 = np.arange(200).reshape((1, 200, 1))
coord3 = np.asarray([0, 1]).reshape((1, 1, 2))
c[coord1, coord2, coord3]

-n