[Matrix-SIG] Creative slicing

[Gary Strangman]

> Try this code. It is not general for the number of dimensions, but
> general for the length of the dimensions. Perhaps it is a start...

<snip>

A good start indeed.  Thanks!  I just thought I'd pass along the following
code to any others who might be interested.  (In the FWIW category, I
agree that this kind of slicing would be an excellent addition to the
NumPy core.)

> x=arange(10*4*5)
> x.shape = (10,4,5) 
> y = ones((4,5))*2
> y[0,0]=3
> y[0,1]=4
> y[1,0]=5
> y.shape = (4,5)   # filled with ints in the 0-9 range 
>  
> # These lines need to be adapted for more dimensions
> yi,xi=indices(y.shape)
> linindex = y*product(x.shape[1:])+yi*y.shape[-2]+xi

This wasn't quite right ... but very close (an off-by-one error, probably
a mail typo).  Instead use ... 

linindex = y*(product(x.shape[1:]))+yi*y.shape[1]+xi

> z=take(ravel(x),ravel(linindex))
> z.shape=(4,5)

Starting from your idea, I modified the code to work in the N-D case. 
It's not exactly pretty (and I didn't bother to find a way around the
for-loop) but it seems to work.  This general type of workaround is not
going to be fast for large arrays ... which is what would make it nice to
have such a function ... a[indices] = vals ... in the numpy core. ;-)

Gary

<------CUT BELOW------>

def elementslice(a, indexarray, dimension=0):
    """
Usage:   elementslice(a, indexarray, dimension=0)
Returns: a[indices] along 'dimension', resulting shape = indices.shape
"""
    if dimension <> 0:  # for convenience, put operating dimension first
        newdims = range(len(a.shape))
        del newdims[dimension]
        newdims = [dimension] + newdims
        a = N.transpose(a,newdims)
    else:
        newdims = a.shape
    ind = N.indices(indexarray.shape)
    linindex = indexarray*N.product(a.shape[1:])  # get offsets for 1st dim
    for i in range(len(ind)-1):  # for remaining dims but one, add in offsets
        multindices = N.product(indexarray.shape[i+1:])
        linindex = linindex + ind[i]*multindices
    linindex = linindex + ind[-1]                 # add in last dim's offsets
    slice = N.take(N.ravel(a),N.ravel(linindex))  # get vals from ravel'd array
    return N.reshape(slice,indexarray.shape)      # reshape properly


--------------------------------------------------------------
Gary Strangman, PhD           |  Neural Systems Group
Office: 617-724-0662          |  Massachusetts General Hospital
Home:   xxx-xxx-xxxx          |  13th Street, Bldg 149, Room 9103
strang at nmr.mgh.harvard.edu |  Charlestown, MA  02129
http://www.nmr.mgh.harvard.edu/Neural_Systems_Group/gary/