Numpy: Multidimensional index. Row by row with no loop

Question

I have a Nx2x2x2 array called A. I also have a Nx2 array called B, which tells me the position of the last two dimensions of A in which I am interested. I am currently getting a Nx2 array, either by using a loop (as in C in the code below) or by using a list comprehension (as in D in the code below). I want to know whether there would be time gains from vectorization, and, if so, how to vectorize this task.

My current approach to vectorization (E in the code below), is using B to index each submatrix of A, so it does not return what I want. I want E to return the same as C or D.

Input:

A=np.reshape(np.arange(0,32),(4,2,2,2))
print("A")
print(A)
B=np.concatenate((np.array([0,1,0,1])[:,np.newaxis],np.array([1,1,0,0])[:,np.newaxis]),axis=1)
print("B")
print(B)
C=np.empty(shape=(4,2))
for n in range(0, 4):
    C[n,:]=A[n,:,B[n,0],B[n,1]]
print("C")
print(C)
D = np.array([A[n,:,B[n,0],B[n,1]] for n in range(0, 4)])
print("D")
print(D)
E=A[:,:,B[:,0],B[:,1]]
print("E")
print(E)

Output:

A
[[[[ 0  1]
   [ 2  3]]

  [[ 4  5]
   [ 6  7]]]


 [[[ 8  9]
   [10 11]]

  [[12 13]
   [14 15]]]


 [[[16 17]
   [18 19]]

  [[20 21]
   [22 23]]]


 [[[24 25]
   [26 27]]

  [[28 29]
   [30 31]]]]
B
[[0 1]
 [1 1]
 [0 0]
 [1 0]]
C
[[  1.   5.]
 [ 11.  15.]
 [ 16.  20.]
 [ 26.  30.]]
D
[[ 1  5]
 [11 15]
 [16 20]
 [26 30]]
E
[[[ 1  3  0  2]
  [ 5  7  4  6]]

 [[ 9 11  8 10]
  [13 15 12 14]]

 [[17 19 16 18]
  [21 23 20 22]]

 [[25 27 24 26]
  [29 31 28 30]]]

Divakar · Accepted Answer · 2015-06-02 23:34:30Z

The complicated slicing operation could be done in a vectorized manner like so -

shp = A.shape
out = A.reshape(shp[0],shp[1],-1)[np.arange(shp[0]),:,B[:,0]*shp[3] + B[:,1]]

You are using the first and second columns of B to index into the third and fourth dimensions of input 4D array, A. With it means is, basically you are slicing the 4D array, with the last two dimensions being fused together. So, you need to get the linear indices with that fused format using B. Of course, before doing all that, you need to reshape A to a 3D array with A.reshape(shp[0],shp[1],-1).

Verify results for a generic 4D array case -

In [104]: A = np.random.rand(6,3,4,5)
     ...: B = np.concatenate((np.random.randint(0,4,(6,1)),np.random.randint(0,5,(6,1))),1)
     ...: 

In [105]: C=np.empty(shape=(6,3))
     ...: for n in range(0, 6):
     ...:     C[n,:]=A[n,:,B[n,0],B[n,1]]
     ...:     

In [106]: shp = A.shape
     ...: out = A.reshape(shp[0],shp[1],-1)[np.arange(shp[0]),:,B[:,0]*shp[3] + B[:,1]]
     ...: 

In [107]: np.allclose(C,out)
Out[107]: True

Collectives™ on Stack Overflow

Numpy: Multidimensional index. Row by row with no loop

1 Answer 1

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

Comments

Your Answer

Sign up or log in

Post as a guest

Related