With readability, np.choose definitely looks great.

If performance is of essence, you can calculate the linear indices and then use np.take or use a flattened version with .ravel() and extract those specific elements from val_arr. The implementation would look something like this -

def linidx_take(val_arr,z_indices):

    # Get number of columns and rows in values array
     _,nC,nR = val_arr.shape

     # Get linear indices and thus extract elements with np.take
    idx = nC*nR*z_indices + nR*np.arange(nR)[:,None] + np.arange(nC)
    return np.take(val_arr,idx) # Or val_arr.ravel()[idx]

Runtime tests and verify results -

Ogrid based solution from here is made into a generic version for these tests, like so :

In [182]: def ogrid_based(val_arr,z_indices):
     ...:   v_shp = val_arr.shape
     ...:   y,x = np.ogrid[0:v_shp[1], 0:v_shp[2]]
     ...:   return val_arr[z_indices, y, x]
     ...: 

Case #1: Smaller datasize

In [183]: val_arr = np.random.rand(30,30,30)
     ...: z_indices = np.random.randint(0,30,(30,30))
     ...: 

In [184]: np.allclose(z_indices.choose(val_arr),ogrid_based(val_arr,z_indices))
Out[184]: True

In [185]: np.allclose(z_indices.choose(val_arr),linidx_take(val_arr,z_indices))
Out[185]: True

In [187]: %timeit z_indices.choose(val_arr)
1000 loops, best of 3: 230 µs per loop

In [188]: %timeit ogrid_based(val_arr,z_indices)
10000 loops, best of 3: 54.1 µs per loop

In [189]: %timeit linidx_take(val_arr,z_indices)
10000 loops, best of 3: 30.3 µs per loop

Case #2: Bigger datasize

In [191]: val_arr = np.random.rand(300,300,300)
     ...: z_indices = np.random.randint(0,300,(300,300))
     ...: 

In [192]: z_indices.choose(val_arr) # Seems like there is some limitation here with bigger arrays.
Traceback (most recent call last):

  File "<ipython-input-192-10c3bb600361>", line 1, in <module>
    z_indices.choose(val_arr)

ValueError: Need between 2 and (32) array objects (inclusive).


In [194]: np.allclose(linidx_take(val_arr,z_indices),ogrid_based(val_arr,z_indices))
Out[194]: True

In [195]: %timeit ogrid_based(val_arr,z_indices)
100 loops, best of 3: 3.67 ms per loop

In [196]: %timeit linidx_take(val_arr,z_indices)
100 loops, best of 3: 2.04 ms per loop

If I understand correctly, you are looking for advances indexing of first dimension of B. You can use np.indices to create the indices required for the other two dimensions of B and use advanced indexing:

idx = np.indices(A.shape)
C = B[A,idx[0],idx[1]]

Example:

B = np.random.rand(10,20,30)
A = np.array([[   1,    0,    0,    1,    0],
       [   0,    2,    0,    2,    4],
       [   1,    7,    0,    5,    5],
       [   2,    1,    0,    1,    7],
       [   0,    7,    0,    1,    4],
       [   0,    0,    3,    6,    0]])

print(C[4,1]==B[7,4,1])
#True

In [107]: arr_3d = np.arange(2*3*4).reshape(2,3,4)                                                           
In [108]: arr_2d = np.array(([3,2,0], [2,3,2]))                                                              
In [109]: arr_2d.shape                                                                                       
Out[109]: (2, 3)
In [110]: arr_3d[[[0],[1]],[0,1,2],arr_2d]                                                                   
Out[110]: 
array([[ 3,  6,  8],
       [14, 19, 22]])

[[0],[1]],[0,1,2] broadcast with each other to index a (2,3) block, the same size as `arr_2d.

ix_ can be used to construct those 2 indices:

In [114]: I,J = np.ix_(range(2), range(3))                                                                   
In [115]: I,J                                                                                                
Out[115]: 
(array([[0],
        [1]]), array([[0, 1, 2]]))
In [116]: arr_3d[I, J, arr_2d]                                                                               
Out[116]: 
array([[ 3,  6,  8],
       [14, 19, 22]])

The question is somewhat ambiguous, but if the intent is to increment some unknown array B at indices (0,0,0), (1,0,1), (3,1,0), and (3,1,1), then the following should be fine:

B[(A.ravel(), ) + np.unravel_index(np.arange(np.prod(A.shape)), A.shape)] += increment

For example:

A = np.array([[0,1],[3,3]])
B = np.zeros((4,2,2), dtype=int)
increment = 1

B[(A.ravel(), ) + np.unravel_index(np.arange(np.prod(A.shape)), A.shape)] += increment

>>> B
array([[[1, 0],
        [0, 0]],

       [[0, 1],
        [0, 0]],

       [[0, 0],
        [0, 0]],

       [[0, 0],
        [1, 1]]])

Another way of doing the same thing is:

w, h = A.shape
indices = (A.ravel(),) + tuple(np.mgrid[:w, :h].reshape(2, -1))

# then
B[indices] += increment

If your matrices get much bigger than 3x3, to the point that memory involved in np.ogrid is an issue, and if your indexes remain binary, you could also do:

np.where(a, b[1], b[0])

But other than that corner case (or if you like code golfing one-liners) the other answer is probably better.

You can use take array method:

import numpy

a = numpy.array([[5,6,7,8,9],[10,11,12,14,15]])

idx = numpy.random.randint(5, size=(3, a.shape[0], a.shape[1]))

print a.take(idx)

To get component #0, use

data[:, 0]

i.e. we get every entry on axis 0 (samples), and only entry #0 on axis 1 (components), and implicitly everything on the remaining axes.

This can be easily generalized to

data[:, indices]

to select all relevant components.

But what OP really wants is just the diagonal of this array, i.e. (data[0, indices[0]], (data[1, indices[1]]), ...) The diagonal of a high-dimensional array can be extracted using the diagonal function:

>>> np.diagonal(data[:, indices])
array([[7, 7, 4, 8, 5],
       [4, 3, 5, 2, 8]])

(You may need to transpose the result.)

This works:

sub = arr[[0,1], :, [1,3]]

Output:

>>> sub
array([[ 1,  5,  9],
       [15, 19, 23]])

A more dynamic version by @Psidom:

>>> sub = arr[np.arange(len(arr)), :, [1,3]]
array([[ 1,  5,  9],
       [15, 19, 23]])