¿Cuál es la forma más elegante de acceder a una matriz de n dimensiones con una matriz de dimensiones (n-1) a lo largo de una dimensión dada como en el ejemplo ficticio

a = np.random.random_sample((3,4,4))
b = np.random.random_sample((3,4,4))
idx = np.argmax(a, axis=0)

¿Cómo puedo acceder ahora conidx a para obtener los máximos ena como si hubiera usadoa.max(axis=0)? o cómo recuperar los valores especificados poridx enb?

Pensé en usarnp.meshgrid Pero creo que es una exageración. Tenga en cuenta que la dimensiónaxis puede ser cualquier eje útil (0,1,2) y no se conoce de antemano. ¿Hay alguna forma elegante de hacer esto?