docs: def ml_acr(tree, character, prediction_method, model, states, avg_br_len, num_nodes, num_tips, freqs=None, sf=None,
             kappa=None, force_joint=True):
      """
      Calculates ML states on the tree and stores them in the corresponding feature.
  
      :param states: numpy array of possible states
      :param prediction_method: str, MPPA (marginal approximation), MAP (max a posteriori) or JOINT
      :param tree: ete3.Tree, the tree of interest
      :param character: str, character for which the ML states are reconstructed
      :param model: str, evolutionary model, F81 (Felsenstein 81-like), JC (Jukes-Cantor-like) or EFT (estimate from tips)
      :param avg_br_len: float, average non-zero branch length of the tree.
      :param freqs: numpy array of predefined frequencies (or None if they are to be estimated)
      :param sf: float, predefined scaling factor (or None if it is to be estimated)
      :return: dict, mapping between reconstruction parameters and values
      """
      n = len(states)
      state2index = dict(zip(states, range(n)))
      missing_data = 0.
      observed_frequencies = np.zeros(n, np.float64)
      for _ in tree:
          state = getattr(_, character, set())
          if state:
              num_node_states = len(state)
              for _ in state:
                  observed_frequencies[state2index[_]] += 1. / num_node_states
          else:
              missing_data += 1
      total_count = observed_frequencies.sum() + missing_data
      observed_frequencies /= observed_frequencies.sum()
      missing_data /= total_count
  
      logger = logging.getLogger('pastml')
      logger.debug('Observed frequencies for {}:{}{}.'
                   .format(character,
                           ''.join('\n\tfrequency of {}:\t{:.3f}'.format(state, observed_frequencies[
                               state2index[state]])
                                   for state in states),
                           '\n\tfraction of missing data:\t{:.3f}'.format(
                               missing_data) if missing_data else ''))
  
      if freqs is not None and model not in {F81, HKY}:
          logging.warning('Some frequencies were specified in the parameter file, '
                          'but the selected model ({}) ignores them. '
                          'Use F81 (or HKY for nucleotide characters only) '
                          'for taking user-specified frequencies into account.'.format(model))
      optimise_frequencies = model in {F81, HKY} and freqs is None
      if JTT == model:
          frequencies = JTT_FREQUENCIES
      elif EFT == model:
          frequencies = observed_frequencies
      elif model in {F81, HKY} and freqs is not None:
          frequencies = freqs
      else:
          frequencies = np.ones(n, dtype=np.float64) / n
  
      initialize_allowed_states(tree, character, states)
      alter_zero_tip_allowed_states(tree, character)
      if sf:
          optimise_sf = False
      else:
          sf = 1. / avg_br_len
          optimise_sf = True
      if HKY == model:
          if kappa:
              optimise_kappa = False
          else:
              optimise_kappa = True
              kappa = 4.
      else:
          optimise_kappa = False
  
      likelihood = get_bottom_up_likelihood(tree=tree, character=character, frequencies=frequencies, sf=sf, kappa=kappa,
                                            is_marginal=True, model=model)
      if not optimise_sf and not optimise_frequencies and not optimise_kappa:
          logger.debug('All the parameters are fixed for {}:{}{}{}{}.'
                       .format(character,
                               ''.join('\n\tfrequency of {}:\t{:.3f}'.format(state, frequencies[
                                   state2index[state]])
                                       for state in states),
                               '\n\tSF:\t{:.3f}, i.e. {:.3f} changes per avg branch'
                               .format(sf, sf * avg_br_len),
                               '\n\tkappa:\t{:.3f}'.format(kappa) if HKY == model else '',
                               '\n\tlog likelihood:\t{:.3f}'.format(likelihood)))
      else:
          logger.debug('Initial values for {} parameter optimisation:{}{}{}{}.'
                       .format(character,
                               ''.join('\n\tfrequency of {}:\t{:.3f}'.format(state, frequencies[
                                   state2index[state]])
                                       for state in states),
                               '\n\tSF:\t{:.3f}, i.e. {:.3f} changes per avg branch'
                               .format(sf, sf * avg_br_len),
                               '\n\tkappa:\t{:.3f}'.format(kappa) if HKY == model else '',
                               '\n\tlog likelihood:\t{:.3f}'.format(likelihood)))
          if optimise_sf:
              (_, sf, _), likelihood = optimize_likelihood_params(tree=tree, character=character, frequencies=frequencies,
                                                                  sf=sf, kappa=kappa,
                                                                  optimise_frequencies=False, optimise_sf=optimise_sf,
                                                                  optimise_kappa=False, avg_br_len=avg_br_len,
                                                                  model=model)
              if optimise_frequencies or optimise_kappa:
                  logger.debug('Pre-optimised SF for {}:{}{}.'
                               .format(character,
                                       '\n\tSF:\t{:.3f}, i.e. {:.3f} changes per avg branch'
                                       .format(sf, sf * avg_br_len),
                                       '\n\tlog likelihood:\t{:.3f}'.format(likelihood)))
          if optimise_frequencies or optimise_kappa:
              (frequencies, sf, kappa), likelihood = \
                  optimize_likelihood_params(tree=tree, character=character, frequencies=frequencies, sf=sf, kappa=kappa,
                                             optimise_frequencies=optimise_frequencies, optimise_sf=optimise_sf,
                                             optimise_kappa=optimise_kappa, avg_br_len=avg_br_len, model=model)
  
          logger.debug('Optimised {} values:{}{}{}{}'
                       .format(character,
                               ''.join('\n\tfrequency of {}:\t{:.3f}'.format(state, frequencies[
                                   state2index[state]])
                                       for state in states) if optimise_frequencies else '',
                               '\n\tSF:\t{:.3f}, i.e. {:.3f} changes per avg branch'
                               .format(sf, sf * avg_br_len),
                               '\n\tkappa:\t{:.3f}'.format(kappa) if HKY == model else '',
                               '\n\tlog likelihood:\t{:.3f}'.format(likelihood)))
  
      result = {LOG_LIKELIHOOD: likelihood, CHARACTER: character, METHOD: prediction_method, MODEL: model,
                FREQUENCIES: frequencies, SCALING_FACTOR: sf, CHANGES_PER_AVG_BRANCH: sf * avg_br_len, STATES: states,
                NUM_NODES: num_nodes, NUM_TIPS: num_tips}
      if HKY == model:
          result[KAPPA] = kappa
  
      results = []
  
      def process_reconstructed_states(method):
          if method == prediction_method or is_meta_ml(prediction_method):
              method_character = get_personalized_feature_name(character, method) \
                  if prediction_method != method else character
              convert_allowed_states2feature(tree, character, states, method_character)
              res = result.copy()
              res[CHARACTER] = method_character
              res[METHOD] = method
              results.append(res)
  
      def process_restricted_likelihood_and_states(method):
          alter_zero_tip_allowed_states(tree, character)
          restricted_likelihood = get_bottom_up_likelihood(tree=tree, character=character,
                                                           frequencies=frequencies, sf=sf, kappa=kappa,
                                                           is_marginal=True, model=model)
          unalter_zero_tip_allowed_states(tree, character, state2index)
          note_restricted_likelihood(method, restricted_likelihood)
          process_reconstructed_states(method)
  
      def note_restricted_likelihood(method, restricted_likelihood):
          logger.debug('Log likelihood for {} after {} state selection:\t{:.3f}'
                       .format(character, method, restricted_likelihood))
          result[RESTRICTED_LOG_LIKELIHOOD_FORMAT_STR.format(method)] = restricted_likelihood
  
      if prediction_method != MAP:
          # Calculate joint restricted likelihood
          restricted_likelihood = get_bottom_up_likelihood(tree=tree, character=character,
                                                           frequencies=frequencies, sf=sf, kappa=kappa,
                                                           is_marginal=False, model=model)
          note_restricted_likelihood(JOINT, restricted_likelihood)
          unalter_zero_tip_joint_states(tree, character, state2index)
          choose_ancestral_states_joint(tree, character, states, frequencies)
          process_reconstructed_states(JOINT)
  
      if is_marginal(prediction_method):
          initialize_allowed_states(tree, character, states)
          alter_zero_tip_allowed_states(tree, character)
          get_bottom_up_likelihood(tree=tree, character=character, frequencies=frequencies, sf=sf, kappa=kappa,
                                   is_marginal=True, model=model)
          calculate_top_down_likelihood(tree, character, frequencies, sf, kappa=kappa, model=model)
          unalter_zero_tip_allowed_states(tree, character, state2index)
          calculate_marginal_likelihoods(tree, character, frequencies)
          # check_marginal_likelihoods(tree, feature)
          result[MARGINAL_PROBABILITIES] = convert_likelihoods_to_probabilities(tree, character, states)
  
          choose_ancestral_states_map(tree, character, states)
          process_restricted_likelihood_and_states(MAP)
  
          if MPPA == prediction_method or is_meta_ml(prediction_method):
  
              if ALL == prediction_method:
                  pars_acr_results = parsimonious_acr(tree, character, MP, states, num_nodes, num_tips)
                  results.extend(pars_acr_results)
                  for pars_acr_res in pars_acr_results:
                      _parsimonious_states2allowed_states(tree, pars_acr_res[CHARACTER], character, state2index)
                      alter_zero_tip_allowed_states(tree, character)
                      restricted_likelihood = get_bottom_up_likelihood(tree=tree, character=character,
                                                                       frequencies=frequencies, sf=sf, kappa=kappa,
                                                                       is_marginal=True, model=model)
                      note_restricted_likelihood(pars_acr_res[METHOD], restricted_likelihood)
  
              result[NUM_SCENARIOS], result[NUM_UNRESOLVED_NODES], result[NUM_STATES_PER_NODE] = \
                  choose_ancestral_states_mppa(tree, character, states, force_joint=force_joint)
              result[NUM_STATES_PER_NODE] /= num_nodes
              result[PERC_UNRESOLVED] = result[NUM_UNRESOLVED_NODES] * 100 / num_nodes
              logger.debug('{} node{} unresolved ({:.2f}%) for {} by {}, '
                           'i.e. {:.4f} state{} per node in average.'
                           .format(result[NUM_UNRESOLVED_NODES], 's are' if result[NUM_UNRESOLVED_NODES] != 1 else ' is',
                                   result[PERC_UNRESOLVED], character, MPPA,
                                   result[NUM_STATES_PER_NODE], 's' if result[NUM_STATES_PER_NODE] > 1 else ''))
              process_restricted_likelihood_and_states(MPPA)
  
      return results