.. highlight:: none

    Implements the Nelder-Mead algorithm described in Lagarias et al. (1998)
    modified to maximize instead of minimizing. JIT-compiled in `nopython`
    mode using Numba.

    Parameters
    ----------
    fun : callable
        The objective function to be maximized.
            `fun(x, *args) -> float`
        where x is an 1-D array with shape (n,) and args is a tuple of the
        fixed parameters needed to completely specify the function. This
        function must be JIT-compiled in `nopython` mode using Numba.

    vertices : ndarray(float, ndim=2)
        Initial simplex with shape (n+1, n) to be modified in-place.

    args : tuple, optional
        Extra arguments passed to the objective function.

    ρ : scalar(float), optional(default=1.)
        Reflection parameter. Must be strictly greater than 0.

    χ : scalar(float), optional(default=2.)
        Expansion parameter. Must be strictly greater than max(1, ρ).

    γ : scalar(float), optional(default=0.5)
        Contraction parameter. Must be stricly between 0 and 1.

    σ : scalar(float), optional(default=0.5)
        Shrinkage parameter. Must be strictly between 0 and 1.

    tol_f : scalar(float), optional(default=1e-10)
        Tolerance to be used for the function value convergence test.

    tol_x : scalar(float), optional(default=1e-10)
        Tolerance to be used for the function domain convergence test.

    max_iter : scalar(float), optional(default=1000)
        The maximum number of allowed iterations.

    Returns
    ----------
    results : namedtuple
        A namedtuple containing the following items:
        ::

            "x" : Approximate solution
            "fun" : Approximate local maximum
            "success" : 1 if successfully terminated, 0 otherwise
            "nit" : Number of iterations
            "final_simplex" : The vertices of the final simplex

def _nelder_mead_algorithm(fun, vertices, bounds=np.array([[], []]).T,
                           args=(), ρ=1., χ=2., γ=0.5, σ=0.5, tol_f=1e-8,
                           tol_x=1e-8, max_iter=1000):
    """
    .. highlight:: none

    Implements the Nelder-Mead algorithm described in Lagarias et al. (1998)
    modified to maximize instead of minimizing. JIT-compiled in `nopython`
    mode using Numba.

    Parameters
    ----------
    fun : callable
        The objective function to be maximized.
            `fun(x, *args) -> float`
        where x is an 1-D array with shape (n,) and args is a tuple of the
        fixed parameters needed to completely specify the function. This
        function must be JIT-compiled in `nopython` mode using Numba.

    vertices : ndarray(float, ndim=2)
        Initial simplex with shape (n+1, n) to be modified in-place.

    args : tuple, optional
        Extra arguments passed to the objective function.

    ρ : scalar(float), optional(default=1.)
        Reflection parameter. Must be strictly greater than 0.

    χ : scalar(float), optional(default=2.)
        Expansion parameter. Must be strictly greater than max(1, ρ).

    γ : scalar(float), optional(default=0.5)
        Contraction parameter. Must be stricly between 0 and 1.

    σ : scalar(float), optional(default=0.5)
        Shrinkage parameter. Must be strictly between 0 and 1.

    tol_f : scalar(float), optional(default=1e-10)
        Tolerance to be used for the function value convergence test.

    tol_x : scalar(float), optional(default=1e-10)
        Tolerance to be used for the function domain convergence test.

    max_iter : scalar(float), optional(default=1000)
        The maximum number of allowed iterations.

    Returns
    ----------
    results : namedtuple
        A namedtuple containing the following items:
        ::

            "x" : Approximate solution
            "fun" : Approximate local maximum
            "success" : 1 if successfully terminated, 0 otherwise
            "nit" : Number of iterations
            "final_simplex" : The vertices of the final simplex

    """
    n = vertices.shape[1]
    _check_params(ρ, χ, γ, σ, bounds, n)

    nit = 0

    ργ = ρ * γ
    ρχ = ρ * χ
    σ_n = σ ** n

    f_val = np.empty(n+1, dtype=np.float64)
    for i in range(n+1):
        f_val[i] = _neg_bounded_fun(fun, bounds, vertices[i], args=args)

    # Step 1: Sort
    sort_ind = f_val.argsort()
    LV_ratio = 1

    # Compute centroid
    x_bar = vertices[sort_ind[:n]].sum(axis=0) / n

    while True:
        shrink = False

        # Check termination
        fail = nit >= max_iter

        best_val_idx = sort_ind[0]
        worst_val_idx = sort_ind[n]

        term_f = f_val[worst_val_idx] - f_val[best_val_idx] < tol_f

        # Linearized volume ratio test (see [2])
        term_x = LV_ratio < tol_x

        if term_x or term_f or fail:
            break

        # Step 2: Reflection
        x_r = x_bar + ρ * (x_bar - vertices[worst_val_idx])
        f_r = _neg_bounded_fun(fun, bounds, x_r, args=args)

        if f_r >= f_val[best_val_idx] and f_r < f_val[sort_ind[n-1]]:
            # Accept reflection
            vertices[worst_val_idx] = x_r
            LV_ratio *= ρ

        # Step 3: Expansion
        elif f_r < f_val[best_val_idx]:
            x_e = x_bar + χ * (x_r - x_bar)
            f_e = _neg_bounded_fun(fun, bounds, x_e, args=args)
            if f_e < f_r:  # Greedy minimization
                vertices[worst_val_idx] = x_e
                LV_ratio *= ρχ
            else:
                vertices[worst_val_idx] = x_r
                LV_ratio *= ρ

        # Step 4 & 5: Contraction and Shrink
        else:
            # Step 4: Contraction
            if f_r < f_val[worst_val_idx]:  # Step 4.a: Outside Contraction
                x_c = x_bar + γ * (x_r - x_bar)
                LV_ratio_update = ργ
            else:  # Step 4.b: Inside Contraction
                x_c = x_bar - γ * (x_r - x_bar)
                LV_ratio_update = γ

            f_c = _neg_bounded_fun(fun, bounds, x_c, args=args)
            if f_c < min(f_r, f_val[worst_val_idx]):  # Accept contraction
                vertices[worst_val_idx] = x_c
                LV_ratio *= LV_ratio_update

            # Step 5: Shrink
            else:
                shrink = True
                for i in sort_ind[1:]:
                    vertices[i] = vertices[best_val_idx] + σ * \
                        (vertices[i] - vertices[best_val_idx])
                    f_val[i] = _neg_bounded_fun(fun, bounds, vertices[i],
                                                args=args)

                sort_ind[1:] = f_val[sort_ind[1:]].argsort() + 1

                x_bar = vertices[best_val_idx] + σ * \
                    (x_bar - vertices[best_val_idx]) + \
                    (vertices[worst_val_idx] - vertices[sort_ind[n]]) / n

                LV_ratio *= σ_n

        if not shrink:  # Nonshrink ordering rule
            f_val[worst_val_idx] = _neg_bounded_fun(fun, bounds,
                                                    vertices[worst_val_idx],
                                                    args=args)

            for i, j in enumerate(sort_ind):
                if f_val[worst_val_idx] < f_val[j]:
                    sort_ind[i+1:] = sort_ind[i:-1]
                    sort_ind[i] = worst_val_idx
                    break

            x_bar += (vertices[worst_val_idx] - vertices[sort_ind[n]]) / n

        nit += 1

    return results(vertices[sort_ind[0]], -f_val[sort_ind[0]], not fail, nit,
                   vertices)