On Wed, Jul 16, 2008 at 17:12, Jack.Cook@shell.com wrote:

Robert,

I can understand how this works if K is a constant time value but in my case K varies at each location in the two-dimensional slice. In other words, if I was doing this in a for loop I would do something like this

for i in range(numI): for j in range(numJ): k = slice(i,j) trace = cube(i,j,k-half_width:k+half_width) # shove trace in sub volume

What am I missing?

Ah, okay. It's a bit tricky, though. Yes, you need to use fancy indexing. Since axis you want to be index fancifully is not the first one, you have to be more explicit than you might otherwise want. For example, it would be great if you could just use slices for the first two axes:

cube[:,:,slice + numpy.arange(-half_width,half_width)]

but the semantics of that are a bit different for reasons I can explain later, if you want. Instead, you have to have explicit fancy-index arrays for the first two axes. Further, the arrays for each axis need to be broadcastable to each other. Fancy indexing will iterate over these broadcasted arrays in parallel with each other to form the new array. The liberal application of numpy.newaxis will help us achieve that. So this is the complete recipe:

In [29]: import numpy In [30]: ni, nj, nk = (10, 15, 20)

# Make a fake data cube such that cube[i,j,k] == k for all i,j,k. In [31]: cube = numpy.empty((ni,nj,nk), dtype=int) In [32]: cube[:,:,:] = numpy.arange(nk)[numpy.newaxis,numpy.newaxis,:]

# Pick out a random fake horizon in k. In [34]: kslice = numpy.random.randint(5, 15, size=(ni, nj)) In [35]: kslice Out[35]: array([[13, 14, 9, 12, 12, 11, 8, 14, 11, 13, 13, 13, 8, 11, 8], [ 7, 12, 12, 6, 10, 12, 9, 11, 13, 9, 14, 11, 5, 12, 12], [ 7, 5, 10, 9, 6, 5, 5, 14, 5, 6, 7, 10, 6, 10, 11], [ 6, 9, 11, 14, 7, 11, 10, 6, 6, 9, 9, 11, 5, 5, 14], [12, 8, 11, 6, 10, 8, 5, 9, 8, 10, 7, 5, 9, 9, 14], [ 9, 8, 10, 9, 10, 12, 10, 10, 6, 10, 11, 6, 8, 7, 7], [11, 12, 7, 13, 5, 5, 8, 14, 5, 14, 9, 10, 12, 7, 14], [ 7, 7, 7, 12, 10, 6, 13, 13, 11, 13, 8, 11, 13, 14, 14], [ 6, 13, 13, 10, 10, 14, 10, 8, 9, 14, 13, 12, 9, 9, 5], [13, 14, 10, 8, 11, 11, 10, 6, 12, 11, 12, 12, 13, 11, 7]])

In [36]: half_width = 3

# These two replace the empty slices for the first two axes. In [37]: idx_i = numpy.arange(ni)[:,numpy.newaxis,numpy.newaxis] In [38]: idx_j = numpy.arange(nj)[numpy.newaxis,:,numpy.newaxis]

# This is the substantive part that actually picks out our window. In [41]: idx_k = kslice[:,:,numpy.newaxis] + numpy.arange(-half_width,half_width+1) In [42]: smallcube = cube[idx_i,idx_j,idx_k] In [43]: smallcube.shape Out[43]: (10, 15, 7)

# Now verify that our window is centered on kslice everywhere: In [47]: smallcube[:,:,3] Out[47]: array([[13, 14, 9, 12, 12, 11, 8, 14, 11, 13, 13, 13, 8, 11, 8], [ 7, 12, 12, 6, 10, 12, 9, 11, 13, 9, 14, 11, 5, 12, 12], [ 7, 5, 10, 9, 6, 5, 5, 14, 5, 6, 7, 10, 6, 10, 11], [ 6, 9, 11, 14, 7, 11, 10, 6, 6, 9, 9, 11, 5, 5, 14], [12, 8, 11, 6, 10, 8, 5, 9, 8, 10, 7, 5, 9, 9, 14], [ 9, 8, 10, 9, 10, 12, 10, 10, 6, 10, 11, 6, 8, 7, 7], [11, 12, 7, 13, 5, 5, 8, 14, 5, 14, 9, 10, 12, 7, 14], [ 7, 7, 7, 12, 10, 6, 13, 13, 11, 13, 8, 11, 13, 14, 14], [ 6, 13, 13, 10, 10, 14, 10, 8, 9, 14, 13, 12, 9, 9, 5], [13, 14, 10, 8, 11, 11, 10, 6, 12, 11, 12, 12, 13, 11, 7]])

In [50]: (smallcube[:,:,3] == kslice).all() Out[50]: True

Clear as mud? I can go into more detail if you like particularly about how newaxis works.

Hi Robert

2008/7/17 Robert Kern robert.kern@gmail.com:

In [42]: smallcube = cube[idx_i,idx_j,idx_k]

Fantastic -- a good way to warm up the brain-circuit in the morning! Is there an easy-to-remember rule that predicts the output shape of the operation above? I'm trying to imaging how the output would change if I altered the dimensions of idx_i or idx_j, but it's hard.

It looks like you can do all sorts of interesting things by manipulation the indices. For example, if I take

In [137]: x = np.arange(12).reshape((3,4))

I can produce either

In [138]: x[np.array([[0,1]]), np.array([[1, 2]])] Out[138]: array([[1, 6]])

or

In [140]: x[np.array([[0],[1]]), np.array([[1], [2]])] Out[140]: array([[1], [6]])

and even

In [141]: x[np.array([[0],[1]]), np.array([[1, 2]])] Out[141]: array([[1, 2], [5, 6]])

or its transpose

In [143]: x[np.array([[0,1]]), np.array([[1], [2]])] Out[143]: array([[1, 5], [2, 6]])

Is it possible to separate the indexing in order to understand it better? My thinking was

cube_i = cube[idx_i,:,:].squeeze() cube_j = cube_i[:,idx_j,:].squeeze() cube_k = cube_j[:,:,idx_k].squeeze()

Not sure what would happen if the original array had single dimensions, though.

Back to the original problem:

In [127]: idx_i.shape Out[127]: (10, 1, 1)

In [128]: idx_j.shape Out[128]: (1, 15, 1)

In [129]: idx_k.shape Out[129]: (10, 15, 7)

For the constant slice case, I guess idx_k also have been (1, 1, 7)?

The construction of the cube could probably be done using only

cube.flat = np.arange(nk)

Fernando is right: this is good food for thought and excellent cookbook material!

Regards Stéfan

2008/7/17 Robert Kern robert.kern@gmail.com:

So the way fancy indexing interacts with slices is a bit tricky, and this is why we couldn't use the nicer syntax of cube[:,:,idx_k]. All axes with fancy indices are collected together. Their index arrays are broadcasted and iterated over. *For each iterate*, all of the slices are collected, and those sliced axes are *added* to the output array. If you had used fancy indexing on all of the axes, then the iterate would be a scalar value pulled from the original array. If you mix fancy indexing and slices, the iterate is the *array* formed by the remaining slices.

So if idx_k is shaped (ni,nj,3), for example, cube[:,:,idx_k] will have the shape (ni,nj,ni,nj,3). So smallcube[:,:,i,j,k]==cube[:,:,idx_k[i,j,k]].

Is that clear, or am I obfuscating the subject more?

Crystal, thank you for taking the time to explain! This is such valuable information; we should consider adding a section to numpy.doc.indexing or wherever is more appropriate.

...

For the constant slice case, I guess idx_k also have been (1, 1, 7)?

The construction of the cube could probably be done using only

cube.flat = np.arange(nk)

[...]

If the RHS cannot be broadcasted to the right shape (in this case because it is not the same length as the final axis of the LHS), an error is raised. I find the reuse of the broadcasting concept to be more memorable, and robust over the (mostly) ad hoc use of plain repetition with .flat.

I've become used to exploiting the repeating property of .flat, and forgot its dangers. Thanks for the reminder!

Cheers Stéfan