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

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