docs: 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)