, M, Nk = a.shape[:axis], a.shape[axis], a.shape[axis+1:]
    J = indices.shape[axis]  # Need not equal M
    out = np.empty(Ni + (J,) + Nk)

    for ii in ndindex(Ni):
        for kk in ndindex(Nk):
            a_1d       = a      [ii + s_[:,] + kk]
            indices_1d = indices[ii + s_[:,] + kk]
            out_1d     = out    [ii + s_[:,] + kk]
            for j in range(J):
                out_1d[j] = a_1d[indices_1d[j]]

Equivalently, eliminating the inner loop, the last two lines would be::

            out_1d[:] = a_1d[indices_1d]

See Also
--------
take : Take along an axis, using the same indices for every 1d slice
put_along_axis :
    Put values into the destination array by matching 1d index and data slices

Examples
--------
>>> import numpy as np

For this sample array

>>> a = np.array([[10, 30, 20], [60, 40, 50]])

We can sort either by using sort directly, or argsort and this function

>>> np.sort(a, axis=1)
array([[10, 20, 30],
       [40, 50, 60]])
>>> ai = np.argsort(a, axis=1)
>>> ai
array([[0, 2, 1],
       [1, 2, 0]])
>>> np.take_along_axis(a, ai, axis=1)
array([[10, 20, 30],
       [40, 50, 60]])

The same works for max and min, if you maintain the trivial dimension
with ``keepdims``:

>>> np.max(a, axis=1, keepdims=True)
array([[30],
       [60]])
>>> ai = np.argmax(a, axis=1, keepdims=True)
>>> ai
array([[1],
       [0]])
>>> np.take_along_axis(a, ai, axis=1)
array([[30],
       [60]])

If we want to get the max and min at the same time, we can stack the
indices first

>>> ai_min = np.argmin(a, axis=1, keepdims=True)
>>> ai_max = np.argmax(a, axis=1, keepdims=True)
>>> ai = np.concatenate([ai_min, ai_max], axis=1)
>>> ai
array([[0, 1],
       [1, 0]])
>>> np.take_along_axis(a, ai, axis=1)
array([[10, 30],
       [40, 60]])
r(