[Numpy-discussion] Logical indexing and higher-dimensional arrays.
Wed Feb 8 12:25:09 CST 2012
On Wed, Feb 8, 2012 at 8:49 AM, Travis Oliphant <firstname.lastname@example.org> wrote:
> On Feb 8, 2012, at 8:29 AM, Sturla Molden wrote:
> > On 08.02.2012 06:01, Travis Oliphant wrote:
> >> Recall that the shape of the output with fancy indexing is determined
> by broadcasting together the indexing objects and using that as the shape
> of the output:
> >> x[ind1, ind2] will produce an output with the shape of "broadcast(ind1,
> ind2)" whose elements are selected by the broadcasted tuple.
> > Thank you for the explanation, Travis.
> > I think my main confusion is why we broadcast indices at all.
> > Broadcasting is something I would expect to happen for binary operations
> > between ndarrays, not for indexing.
> We broadcast because of the "zip-based" semantics of current
> fancy-indexing which is basically an element-by-element operation:
> x[ind1, ind2] will choose it's elements out of x by taken elements from
> ind1 as the first index and elements out of ind2 as the second.
> As an essentially element-by-element operation, if the shapes of the input
> arrays are not exactly the same, it is natural to broadcast them to the
> same shape. This is fairly intuitive if you are used to broadcasting in
> NumPy. The problem is that it does not coincide other intuitions in
> certain cases.
> This behavior actually allows easy-to-access "cross-product" behavior by
> broadcasting a 1-d ind1 shaped as (4,) with a ind2 shaped as (4,1). The
> numpy.ix_ function does the basic reshaping needed so that 0:2, 0:2
> matches [0,1],[0,1] indexing: x[ix_([0,1],[0,1])] is the same as
> There are also some very nice applications where you can select out of a
> 3-d volume a depth-surface defined by indexes like so:
> arr[ i[:,newaxis], j, depth]
> where arr is a 3-d array, i and j are 1-d index arrays: i =
> arange(arr.shape) and j = arange(arr.shape), and depth is a 2-d array
> of "depths". The selected result will be 2-d.
> When you understand what is happening, it is consistent and it does make
> sense and it has some very nice applications. I recognize that people
> get confused because their expectations are different, but the current
> behavior can be understood with a fairly simple mental model.
> I could support, however, a move to push this style of indexing to a
> method, and make fancy-indexing use cross-product semantics if:
> * we do it in a phased-manner with lot's of deprecation warnings
> and even a tool that helps you change code (the tool would have to "play"
> your code and make a note of the locations where this style of indexing was
> used --- because static code analysis wouldn't know which objects were
> arrays and which weren't).
> * we also make the ndarray object general enough so that
> fancy-indexing could be returned as a view on the array (the same as
> regular indexing)
> This sort of thing would take time, but is not out of the question in my
> mind because I suspect the number of users and use-cases of "broadcasted"
> fancy-indexing is small.
Speaking for myself, I've used both methods quite often. For the
broadcasting method, it is very useful for image processing scenarios where
you want to extract a small subregion with an arbitrary shape. It is also
extremely useful for this to work both ways, for example if I wanted to do
some processing on that subregion and then put it back in the larger image:
sub_reg_vals = img[ ind1, ind2 ]
img[ ind1, ind2 ] = sub_reg_vals
Of course this may require awareness of the coordinates versus flat index,
but the ravel/unravel_index functions do this for you. Note that IDL works
in this way, unlike MATLAB. The "quick and dirty" way I've been doing it in
sub_reg_vals = diag( img[ind1, ind2] )
Which works but is obviously a bit inefficient.
I have no convincing ideas about which method should be the "default", but
I think both methods get a fair bit of use, so as long as there are fast
and well documented methods for doing either, I would be happy. I actually
had no idea that np.ix_ did the other method in NumPy. (Thanks for
mentioning that! I had been using A[ind1,:][:,ind2]). There is no generic
terminology for these operations, so I think that generates a lot of
-------------- next part --------------
An HTML attachment was scrubbed...
More information about the NumPy-Discussion