[Numpy-discussion] argsort and take along arbitrary axes

2007-05-14 Thread Zachary Pincus
Hello all,

I've got a few questions that came up as I tried to calculate various  
statistics about an image time-series. For example, I have an array  
of shape (t,x,y) representing t frames of a time-lapse of resolution  
(x,y).

Now, say I want to both argsort and sort this time-series, pixel- 
wise. (For example.)

In 1-d it's easy:
indices = a.argsort()
sorted = a[indices]

I would have thought that doing this on my 3-d array would work  
similarly:
indices = a.argsort(axis=0)
sorted = a.take(indices, axis=0)

Unfortunately, this gives a ValueError of dimensions too large.  
Now, I know that 'a.sort(axis=0)' works fine for the given example,  
but I'm curious about how to this sort of indexing operation in the  
general case.

Thanks for any insight,

Zach Pincus

Program in Biomedical Informatics and Department of Biochemistry
Stanford University School of Medicine

___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] argsort and take along arbitrary axes

2007-05-14 Thread Zachary Pincus
 I've got a few questions that came up as I tried to calculate various
 statistics about an image time-series. For example, I have an array
 of shape (t,x,y) representing t frames of a time-lapse of resolution
 (x,y).

 Now, say I want to both argsort and sort this time-series, pixel-
 wise. (For example.)

 In 1-d it's easy:
 indices = a.argsort()
 sorted = a[indices]

 I would have thought that doing this on my 3-d array would work
 similarly:
 indices = a.argsort(axis=0)
 sorted = a.take(indices, axis=0)

 Unfortunately, this gives a ValueError of dimensions too large.
 Now, I know that 'a.sort(axis=0)' works fine for the given example,
 but I'm curious about how to this sort of indexing operation in the
 general case.

 Unfortunately, argsort doesn't work transparently with take or  
 fancy indexing for multidimensional arrays. I am thinking of adding  
 a function argtake for this, and also for the results returned by  
 argmax and argmin, but at the moment you have to fill in the   
 values of the other indices and use fancy indexing. For now, it is  
 probably simpler, prettier, and faster to just sort the array.

Thanks Charles. Unfortunately, the argsort/sort buisness was, as I  
mentioned, just an example of the kind of 'take' operation that I am  
trying to figure out how to do. There are other operations that will  
have similarly-formatted 'indices' arrays (as above) that aren't  
generated from argsort...

As such, how do I fill in the values of the other indices and use  
fancy indexing? Even after reading the numpy book about that, and  
reading the docstring for numpy.take, I'm still vague on this. Would  
I use numpy.indices to get a list of index arrays, and then swap in  
(at the correct position in this list) the result of argsort (or the  
other operations), and use that for fancy indexing? Is there an  
easier/faster way?

Thanks again,

Zach





___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] argsort and take along arbitrary axes

2007-05-14 Thread Charles R Harris

On 5/14/07, Zachary Pincus [EMAIL PROTECTED] wrote:


 I've got a few questions that came up as I tried to calculate various
 statistics about an image time-series. For example, I have an array
 of shape (t,x,y) representing t frames of a time-lapse of resolution
 (x,y).

 Now, say I want to both argsort and sort this time-series, pixel-
 wise. (For example.)

 In 1-d it's easy:
 indices = a.argsort()
 sorted = a[indices]

 I would have thought that doing this on my 3-d array would work
 similarly:
 indices = a.argsort(axis=0)
 sorted = a.take(indices, axis=0)

 Unfortunately, this gives a ValueError of dimensions too large.
 Now, I know that 'a.sort(axis=0)' works fine for the given example,
 but I'm curious about how to this sort of indexing operation in the
 general case.

 Unfortunately, argsort doesn't work transparently with take or
 fancy indexing for multidimensional arrays. I am thinking of adding
 a function argtake for this, and also for the results returned by
 argmax and argmin, but at the moment you have to fill in the
 values of the other indices and use fancy indexing. For now, it is
 probably simpler, prettier, and faster to just sort the array.

Thanks Charles. Unfortunately, the argsort/sort buisness was, as I
mentioned, just an example of the kind of 'take' operation that I am
trying to figure out how to do. There are other operations that will
have similarly-formatted 'indices' arrays (as above) that aren't
generated from argsort...

As such, how do I fill in the values of the other indices and use
fancy indexing? Even after reading the numpy book about that, and
reading the docstring for numpy.take, I'm still vague on this. Would
I use numpy.indices to get a list of index arrays, and then swap in
(at the correct position in this list) the result of argsort (or the
other operations), and use that for fancy indexing? Is there an
easier/faster way?



I've attached  a quick, mostly untested, version of argtake. It's in Python,
probably not too fast, but see if it works for you.

Chuck
import numpy as N

def argtake(arr, ind, axis=-1) :
Take using output of argsort

*Description*

The usual output of argsort is not easily used with take when the relevent
array is multidimensional. This is a quick and dirty standin.

*Parameters*:

arr : ndarray
The array from which the elements are taken. Typically this will be
the same array to which argsort was applied.

ind : ndarray
Indices returned by argsort or lexsort.

axis : integer
Axis to which argsort or lexsort was applied. Must match the
original call. Defaults to -1.

*Returns*:

reordered_array : same type as arr
The input array reordered by ind along the axis.


if arr.shape != ind.shape :
raise shape mismatch
if arr.ndim == 1 :
return N.take(arr, ind)
else :
naxis = arr.shape[axis]
aswap = N.swapaxes(arr, axis, -1)
iswap = N.swapaxes(ind, axis, -1)
shape = aswap.shape
aswap = aswap.reshape(-1, naxis)
iswap = iswap.reshape(-1, naxis)
oswap = N.empty_like(aswap)
for i in xrange(len(aswap)) :
N.take(aswap[i], iswap[i], out=oswap[i])
oswap = oswap.reshape(shape)
oswap = N.swapaxes(oswap, axis, -1)
return oswap
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion