docs: def optimize_likelihood_params(tree, character, frequencies, sf, kappa, avg_br_len,
                                 optimise_sf=True, optimise_frequencies=True, optimise_kappa=True,
                                 model=F81):
      """
      Optimizes the likelihood parameters (state frequencies and scaling factor) for the given tree.
  
      :param model: model of character evolution
      :type model: str
      :param avg_br_len: avg branch length
      :type avg_br_len: float
      :param tree: tree of interest
      :type tree: ete3.Tree
      :param character: character for which the likelihood is optimised
      :type character: str
      :param frequencies: array of initial state frequencies
      :type frequencies: numpy.array
      :param sf: initial scaling factor
      :type sf: float
      :param optimise_sf: whether the scaling factor needs to be optimised
      :type optimise_sf: bool
      :param optimise_frequencies: whether the state frequencies need to be optimised
      :type optimise_frequencies: bool
      :return: optimized parameters and log likelihood: ((frequencies, scaling_factor), optimum)
      :rtype: tuple
      """
      bounds = []
      if optimise_frequencies:
          bounds += [np.array([1e-6, 10e6], np.float64)] * (len(frequencies) - 1)
      if optimise_sf:
          bounds += [np.array([0.001 / avg_br_len, 10. / avg_br_len])]
      if optimise_kappa:
          bounds += [np.array([1e-6, 20.])]
      bounds = np.array(bounds, np.float64)
  
      def get_real_params_from_optimised(ps):
          freqs = frequencies
          if optimise_frequencies:
              freqs = np.hstack((ps[: (len(frequencies) - 1)], [1.]))
              freqs /= freqs.sum()
          sf_val = ps[(len(frequencies) - 1) if optimise_frequencies else 0] if optimise_sf else sf
          kappa_val = ps[((len(frequencies) - 1) if optimise_frequencies else 0) + (1 if optimise_sf else 0)] \
              if optimise_kappa else kappa
          return freqs, sf_val, kappa_val
  
      def get_v(ps):
          if np.any(pd.isnull(ps)):
              return np.nan
          freqs, sf_val, kappa_val = get_real_params_from_optimised(ps)
          res = get_bottom_up_likelihood(tree=tree, character=character, frequencies=freqs,
                                         sf=sf_val, kappa=kappa_val, is_marginal=True, model=model)
          return np.inf if pd.isnull(res) else -res
  
      for i in range(10):
          if i == 0:
              vs = np.hstack((frequencies[:-1] / frequencies[-1] if optimise_frequencies else [],
                              [sf] if optimise_sf else [],
                              [kappa] if optimise_kappa else []))
          else:
              vs = np.random.uniform(bounds[:, 0], bounds[:, 1])
          fres = minimize(get_v, x0=vs, method='L-BFGS-B', bounds=bounds)
          if fres.success and not np.any(np.isnan(fres.x)):
              return get_real_params_from_optimised(fres.x), -fres.fun