docs: def choose_ancestral_states_mppa(tree, feature, states, force_joint=True):
      """
      Chooses node ancestral states based on their marginal probabilities using MPPA method.
  
      :param force_joint: make sure that Joint state is chosen even if it has a low probability.
      :type force_joint: bool
      :param tree: tree of interest
      :type tree: ete3.Tree
      :param feature: character for which the ancestral states are to be chosen
      :type feature: str
      :param states: possible character states in order corresponding to the probabilities array
      :type states: numpy.array
      :return: number of ancestral scenarios selected,
          calculated by multiplying the number of selected states for all nodes.
          Also modified the get_personalized_feature_name(feature, ALLOWED_STATES) feature of each node
          to only contain the selected states.
      :rtype: int
      """
      lh_feature = get_personalized_feature_name(feature, LH)
      allowed_state_feature = get_personalized_feature_name(feature, ALLOWED_STATES)
      joint_state_feature = get_personalized_feature_name(feature, JOINT_STATE)
  
      n = len(states)
      _, state2array = get_state2allowed_states(states, False)
  
      num_scenarios = 1
      unresolved_nodes = 0
      num_states = 0
  
      # If force_joint == True,
      # we make sure that the joint state is always chosen,
      # for this we sort the marginal probabilities array as [lowest_non_joint_mp, ..., highest_non_joint_mp, joint_mp]
      # select k in 1:n such as the correction between choosing 0, 0, ..., 1/k, ..., 1/k and our sorted array is min
      # and return the corresponding states
      for node in tree.traverse():
          marginal_likelihoods = getattr(node, lh_feature)
          marginal_probs = marginal_likelihoods / marginal_likelihoods.sum()
          if force_joint:
              joint_index = getattr(node, joint_state_feature)
              joint_prob = marginal_probs[joint_index]
              marginal_probs = np.hstack((np.sort(np.delete(marginal_probs, joint_index)), [joint_prob]))
          else:
              marginal_probs = np.sort(marginal_probs)
          best_k = n
          best_correstion = np.inf
          for k in range(1, n + 1):
              correction = np.hstack((np.zeros(n - k), np.ones(k) / k)) - marginal_probs
              correction = correction.dot(correction)
              if correction < best_correstion:
                  best_correstion = correction
                  best_k = k
  
          num_scenarios *= best_k
          num_states += best_k
          if force_joint:
              indices_selected = sorted(range(n),
                                        key=lambda _: (0 if n == joint_index else 1, -marginal_likelihoods[_]))[:best_k]
          else:
              indices_selected = sorted(range(n), key=lambda _: -marginal_likelihoods[_])[:best_k]
          if best_k == 1:
              allowed_states = state2array[indices_selected[0]]
          else:
              allowed_states = np.zeros(len(states), dtype=np.int)
              allowed_states[indices_selected] = 1
              unresolved_nodes += 1
          node.add_feature(allowed_state_feature, allowed_states)
  
      return num_scenarios, unresolved_nodes, num_states