docs: def solve_mbar_once(u_kn_nonzero, N_k_nonzero, f_k_nonzero, method="hybr", tol=1E-12, options=None):
      """Solve MBAR self-consistent equations using some form of equation solver.
  
      Parameters
      ----------
      u_kn_nonzero : np.ndarray, shape=(n_states, n_samples), dtype='float'
          The reduced potential energies, i.e. -log unnormalized probabilities
          for the nonempty states
      N_k_nonzero : np.ndarray, shape=(n_states), dtype='int'
          The number of samples in each state for the nonempty states
      f_k_nonzero : np.ndarray, shape=(n_states), dtype='float'
          The reduced free energies for the nonempty states
      method : str, optional, default="hybr"
          The optimization routine to use.  This can be any of the methods
          available via scipy.optimize.minimize() or scipy.optimize.root().
      tol : float, optional, default=1E-14
          The convergance tolerance for minimize() or root()
      verbose: bool
          Whether to print information about the solution method.
      options: dict, optional, default=None
          Optional dictionary of algorithm-specific parameters.  See
          scipy.optimize.root or scipy.optimize.minimize for details.
  
      Returns
      -------
      f_k : np.ndarray
          The converged reduced free energies.
      results : dict
          Dictionary containing entire results of optimization routine, may
          be useful when debugging convergence.
  
      Notes
      -----
      This function requires that N_k_nonzero > 0--that is, you should have
      already dropped all the states for which you have no samples.
      Internally, this function works in a reduced coordinate system defined
      by subtracting off the first component of f_k and fixing that component
      to be zero.
  
      For fast but precise convergence, we recommend calling this function
      multiple times to polish the result.  `solve_mbar()` facilitates this.
      """
      u_kn_nonzero, N_k_nonzero, f_k_nonzero = validate_inputs(u_kn_nonzero, N_k_nonzero, f_k_nonzero)
      f_k_nonzero = f_k_nonzero - f_k_nonzero[0]  # Work with reduced dimensions with f_k[0] := 0
      u_kn_nonzero = precondition_u_kn(u_kn_nonzero, N_k_nonzero, f_k_nonzero)
  
      pad = lambda x: np.pad(x, (1, 0), mode='constant')  # Helper function inserts zero before first element
      unpad_second_arg = lambda obj, grad: (obj, grad[1:])  # Helper function drops first element of gradient
  
      # Create objective functions / nonlinear equations to send to scipy.optimize, fixing f_0 = 0
      grad = lambda x: mbar_gradient(u_kn_nonzero, N_k_nonzero, pad(x))[1:]  # Objective function gradient
      grad_and_obj = lambda x: unpad_second_arg(*mbar_objective_and_gradient(u_kn_nonzero, N_k_nonzero, pad(x)))  # Objective function gradient and objective function
      hess = lambda x: mbar_hessian(u_kn_nonzero, N_k_nonzero, pad(x))[1:][:, 1:]  # Hessian of objective function
  
      with warnings.catch_warnings(record=True) as w:
          if method in ["L-BFGS-B", "dogleg", "CG", "BFGS", "Newton-CG", "TNC", "trust-ncg", "SLSQP"]:
              if method in ["L-BFGS-B", "CG"]:
                  hess = None  # To suppress warning from passing a hessian function.
              results = scipy.optimize.minimize(grad_and_obj, f_k_nonzero[1:], jac=True, hess=hess, method=method, tol=tol, options=options)
              f_k_nonzero = pad(results["x"])
          elif method == 'adaptive':
              results = adaptive(u_kn_nonzero, N_k_nonzero, f_k_nonzero, tol=tol, options=options)
              f_k_nonzero = results # they are the same for adaptive, until we decide to return more.
          else:
              results = scipy.optimize.root(grad, f_k_nonzero[1:], jac=hess, method=method, tol=tol, options=options)
              f_k_nonzero = pad(results["x"])
  
      # If there were runtime warnings, show the messages
      if len(w) > 0:
          can_ignore = True
          for warn_msg in w:
              if "Unknown solver options" in str(warn_msg.message):
                  continue
              warnings.showwarning(warn_msg.message, warn_msg.category,
                                   warn_msg.filename, warn_msg.lineno, warn_msg.file, "")
              can_ignore = False  # If any warning is not just unknown options, can ]not skip check
          if not can_ignore:
              # Ensure MBAR solved correctly
              w_nk_check = mbar_W_nk(u_kn_nonzero, N_k_nonzero, f_k_nonzero)
              check_w_normalized(w_nk_check, N_k_nonzero)
              print("MBAR weights converged within tolerance, despite the SciPy Warnings. Please validate your results.")
  
      return f_k_nonzero, results