def make_monotonic(values, amplitudes, result):
    """
    Make an ndarray monotonically increasing over its last axis

    This function is a numpy universal function which defines the operation
    applied for every dimension.

    The result argument can be ignored when calling this function.
    It is filled inplace and returned by this function.

    Current method: Assume centre is well behaved. Iterate outward in both
    directions from centre and linearly interpolate values that don't follow
    the monotonically increasing trend. If there is no second value for
    interpolation (that is following the montonically increasing trend),
    then set to last behaved value.

    Parameters
    ----------
    values : ndarray
        The values to ensure are monotonic
        Shape: (*tf_dimensions, n_amplitudes)
    amplitudes : list or ndarray
       Amplitudes for the X axis of the TF
       Shape: (n_amplitudes)
    result : ndarray
        Return argument for ufunc (ignore)
        Returns the resulting monotonic array

    Returns
    -------
    result : ndarray
        Input values transformed to be monotonic
        Shape: (*tf_dimensions, n_amplitudes)
    """
    def interp(x, x0, x1, y0, y1):
        return y0 + (x - x0) * (y1 - y0) / (x1 - x0)

    result[:] = values
    n_entries = values.size
    mid = n_entries // 2
    # First half
    for i in range(mid, 0, -1):
        ii = i - 1
        if result[ii] < result[i]:
            continue
        else:
            for j in range(i, 0, -1):
                jj = j - 1
                if result[jj] < result[i]:
                    result[ii] = interp(
                        amplitudes[ii], amplitudes[i], amplitudes[jj],
                        result[i], result[jj]
                    )
                    break
            else:
                result[ii] = result[i]
    # Second half
    for i in range(mid, n_entries-1, 1):
        ii = i + 1
        if result[ii] > result[i]:
            continue
        else:
            for j in range(i, n_entries-1, 1):
                jj = j + 1
                if result[jj] > result[i]:
                    result[ii] = interp(
                        amplitudes[ii], amplitudes[i], amplitudes[jj],
                        result[i], result[jj]
                    )
                    break
            else:
                result[ii] = result[i]