analysis

phyllotaxis_analysis.analysis

This module contains the functions used to detect permutations in sequences of divergence angles.

code_sequence

code_sequence(seq, permutation_block_max_size, canonical_angle, borders_dict, kappa, threshold)

Make a list of n-admissible trees coding sequence of measured angles.

Given a sequence of measured angles, this function generates a list of n-admissible trees that encode the sequence. The trees are constructed by analyzing segments of the sequence, with each tree corresponding to a specific segment defined by its starting and ending indices.

Parameters:

Name	Type	Description	Default
seq	list[float]	Sequence of measured angles to be analyzed.	required
permutation_block_max_size	int	Maximum number of organs involved in permutations.	required
canonical_angle	float	Canonical divergence angle used in the analysis.	required
borders_dict	dict	Dictionary mapping theoretical angles to their corresponding intervals.	required
kappa	float	Concentration parameter for the analysis.	required
threshold	float	Statistical threshold used to prune the n-admissible trees.	required

Returns:

Type	Description
list[tuple[object, int]]	A list of tuples where each tuple contains an n-admissible tree and its corresponding index in the sequence. Each tree encodes a segment of the sequence from seq[(index - 1): index].

Notes

• The function processes the sequence in segments, starting from the beginning and moving forward.
• For each segment, it constructs an n-admissible tree using the make_n_admissible_tree function.
• The trees and additional data are saved to a file named 'nadm.pkl' using pickle.
• The function assumes the existence of helper functions theoretical_divergence_angles and make_n_admissible_tree.

Source code in src/phyllotaxis_analysis/analysis.py

def code_sequence(seq, permutation_block_max_size, canonical_angle, borders_dict, kappa, threshold):
    """
    Make a list of `n`-admissible trees coding sequence of measured angles.

    Given a sequence of measured angles, this function generates a list of `n`-admissible trees that encode the sequence. The trees are constructed by analyzing segments of the sequence, with each tree corresponding to a specific segment defined by its starting and ending indices.

    Parameters
    ----------
    seq : list[float]
        Sequence of measured angles to be analyzed.
    permutation_block_max_size : int
        Maximum number of organs involved in permutations.
    canonical_angle : float
        Canonical divergence angle used in the analysis.
    borders_dict : dict
        Dictionary mapping theoretical angles to their corresponding intervals.
    kappa : float
        Concentration parameter for the analysis.
    threshold : float
        Statistical threshold used to prune the `n`-admissible trees.

    Returns
    -------
    list[tuple[object, int]]
        A list of tuples where each tuple contains an `n`-admissible tree and its corresponding index in the sequence. Each tree encodes a segment of the sequence from `seq[(index - 1): index]`.

    Notes
    -----
    - The function processes the sequence in segments, starting from the beginning and moving forward.
    - For each segment, it constructs an `n`-admissible tree using the `make_n_admissible_tree` function.
    - The trees and additional data are saved to a file named 'nadm.pkl' using `pickle`.
    - The function assumes the existence of helper functions `theoretical_divergence_angles` and `make_n_admissible_tree`.
    """
    d_n, d_n_coeff = theoretical_divergence_angles(permutation_block_max_size, canonical_angle)
    d_n_coeff_dict = dict((d_n[i], d_n_coeff[i]) for i in range(3 * permutation_block_max_size - 2))
    sequence_length = len(seq)
    trees_indexes = list()
    index = 0
    while index <= sequence_length - 1:
        if index == 0:
            print("Beginning of sequence")
            tree, all_trees_to_pickle = make_n_admissible_tree(seq[index:], permutation_block_max_size, canonical_angle,
                                                               borders_dict, threshold, d_n, d_n_coeff, d_n_coeff_dict,
                                                               kappa,
                                                               seq_begin=True)
        else:
            tree, all_trees_to_pickle = make_n_admissible_tree(seq[index:], permutation_block_max_size, canonical_angle,
                                                               borders_dict, threshold, d_n, d_n_coeff, d_n_coeff_dict,
                                                               kappa,
                                                               d_n_coeff_dict)
        index += tree.get_max_level() + 1
        trees_indexes.append([tree, index - 1])

    return trees_indexes

extract_sequences

extract_sequences(trees_indexes, seq, permutation_block_max_size, canonical_angle, prob_check=True)

Extract n-admissible sequences from a list of n-admissible trees.

Given a sequence of measured angles and corresponding n-admissible trees, this function extracts the theoretical angles that best explain the sequence, identifies unexplained angles, and classifies permutations as valid or invalid.

Parameters:

Name	Type	Description	Default
trees_indexes	list[tuple[Tree, int]]	A list of tuples where each tuple contains an n-admissible tree and its starting index. Each tree encodes the subsequence seq[(index - 1): index].	required
seq	list[float]	Sequence of measured angles to analyze.	required
permutation_block_max_size	int	Maximum number of organs involved in permutations.	required
canonical_angle	float	Canonical divergence angle used in the analysis.	required

Other Parameters:

Name	Type	Description
prob_check	bool	If True, includes probability information in the output. Default is True.

Returns:

Name	Type	Description
prediction	list[float]	List of theoretical angles coding the input sequence.
not_explained_list	list[int]	Indices of angles in the original sequence that were not explained by the model.
valid_permutations	list[list[int]]	List of valid permutations found in the sequence.
invalids	list[list[int]]	List of invalidated permutations.

Source code in src/phyllotaxis_analysis/analysis.py

def extract_sequences(trees_indexes, seq, permutation_block_max_size, canonical_angle, prob_check=True):
    """
    Extract n-admissible sequences from a list of n-admissible trees.

    Given a sequence of measured angles and corresponding n-admissible trees, this function
    extracts the theoretical angles that best explain the sequence, identifies unexplained angles,
    and classifies permutations as valid or invalid.

    Parameters
    ----------
    trees_indexes : list[tuple[Tree, int]]
        A list of tuples where each tuple contains an n-admissible tree and its starting index.
        Each tree encodes the subsequence `seq[(index - 1): index]`.
    seq : list[float]
        Sequence of measured angles to analyze.
    permutation_block_max_size : int
        Maximum number of organs involved in permutations.
    canonical_angle : float
        Canonical divergence angle used in the analysis.

    Other Parameters
    ----------------
    prob_check : bool, optional
        If True, includes probability information in the output. Default is True.

    Returns
    -------
    prediction : list[float]
        List of theoretical angles coding the input sequence.
    not_explained_list : list[int]
        Indices of angles in the original sequence that were not explained by the model.
    valid_permutations : list[list[int]]
        List of valid permutations found in the sequence.
    invalids : list[list[int]]
        List of invalidated permutations.
    """
    theoretical_angles, theoretical_coeff_angles = theoretical_divergence_angles(permutation_block_max_size,
                                                                                 canonical_angle)
    d_n_coeff_dict = dict(
        (theoretical_angles[i], theoretical_coeff_angles[i]) for i in range(3 * permutation_block_max_size - 2))
    index = 0
    best_seq = []
    best_prob = 0
    for tree_index in trees_indexes:
        best_path = []
        seq_list = []
        max_prob = None
        for path_all in tree_index[0].all_paths_all():
            result, inversion_list = is_n_admissible(path_all[0], permutation_block_max_size, canonical_angle,
                                                     d_n_coeff_dict)
            if not result:
                result, inversion_list = is_n_admissible(path_all[0][:-1], permutation_block_max_size, canonical_angle,
                                                         d_n_coeff_dict)
            new_inversion_list = list()
            for inversion in inversion_list:
                new_inversion = [j + index for j in inversion]
                new_inversion_list.append(new_inversion)
            seq_list.append(path_all)
            if max_prob == None:
                max_prob = path_all[3]
            elif max_prob < path_all[3]:
                max_prob = path_all[3]
        for path_all in seq_list:
            if path_all[3] == max_prob:
                best_path.append(path_all)
                break  # ???
        best_prob += max_prob
        best_seq.append(best_path)
        index = tree_index[1] + 1
    ind = 0
    not_explained_list = []
    prediction = []
    inversions = []
    all_prob = []
    for s_list in best_seq:
        for path_all in s_list:
            prediction.extend(path_all[0])
            if prob_check:
                all_prob.extend(path_all[1])
            not_explained_list.extend([j[1] + ind for j in path_all[2]])
            result, inversion_list = is_n_admissible(path_all[0], permutation_block_max_size, canonical_angle,
                                                     d_n_coeff_dict)
            if not result:
                result, inversion_list = is_n_admissible(path_all[0][:-1], permutation_block_max_size, canonical_angle,
                                                         d_n_coeff_dict)
            new_inversion_list = list()
            for inversion in inversion_list:
                new_inversion = [j + ind for j in inversion]
                new_inversion_list.append(new_inversion)
            inversions.extend(new_inversion_list)
            ind += len(path_all[0])
    not_explained_angles = [seq[j] for j in not_explained_list]
    d_n, d_n_com = theoretical_divergence_angles(permutation_block_max_size, canonical_angle)
    code = dict((d_n[j], d_n_com[j]) for j in range(len(d_n)))
    coded_seq_now = [code[j] for j in prediction]

    ss_list = []
    ends = []
    for l in trees_indexes:
        ends.append(l[1])
    if prediction[0: ends[0] + 1] != []:
        sub_seqs = [prediction[0: ends[0] + 1]]
    else:
        sub_seqs = []
    for i in range(len(ends) - 1):
        sub_seqs.append(prediction[ends[i] + 1: ends[i + 1] + 1])
    if prediction[ends[-1] + 1: len(prediction)] != []:
        sub_seqs.append(prediction[ends[-1] + 1: len(seq)])

    invalid_counter = 0
    ind = 0
    for ss in sub_seqs:
        res = is_n_admissible_first_angles(ss, permutation_block_max_size, canonical_angle, d_n_coeff_dict)
        if res[0]:
            ss_list.append([ss, res[1], []])
        elif len(ss) > 1:
            res = is_n_admissible_first_angles(ss[:-1], permutation_block_max_size, canonical_angle, d_n_coeff_dict)
            invalides = sub_invalidate_last(res[1], len(ss) - 1)
            for inv in invalides:
                invalid_counter += len(inv)
            ss_list.append([ss, res[1], invalides])
        else:
            ss_list.append([ss, res[1], []])

        ind += len(ss)
    invalid_counter += len(not_explained_list)
    invalids = invalidate_seq(coded_seq_now, not_explained_list, permutation_block_max_size)
    valid_permutations = lists_complement_as_set(new_inversion_list, invalids)
    return prediction, not_explained_list, valid_permutations, invalids

make_n_admissible_tree

make_n_admissible_tree(seq, permutation_block_max_size, canonical_angle, borders_dict, threshold, d_n, theoretical_coeff_angles, d_n_coeff_dict, kappa, seq_begin=False)

make n-admissible tree coding sequence of measured angles. It stops when it can not code an angle, i.e. at a not explained angle.

Parameters:

Name	Type	Description	Default
seq	list	Sequence of measured angles.	required
permutation_block_max_size	int	Maximum number of organs involved in permutations.	required
canonical_angle	float	Canonical divergence angle.	required
borders_dict	dict	Dictionary of theoretical angles and the corresponding intervals.	required
kappa	float	Concentration parameter.	required
threshold	float	A statistical threshold to prune the n-admissible trees.	required
d_n	list	List of theoretical angles.	required
theoretical_coeff_angles	list	List of coefficients of theoretical angles.	required
d_n_coeff_dict	dict	Dictionary of theoretical angles and their coefficients.	required
seq_begin	bool	Indicates whether permutations at the beginning of the sequence should be taken into account. Default is False.	False

Returns:

Type	Description
Nadmissibletree	The constructed n-admissible tree.
list	The list of all intermediate trees generated during the process.

Source code in src/phyllotaxis_analysis/analysis.py

def make_n_admissible_tree(seq, permutation_block_max_size, canonical_angle, borders_dict, threshold, d_n,
                           theoretical_coeff_angles, d_n_coeff_dict, kappa, seq_begin=False):
    """
    make `n`-admissible tree coding sequence of measured angles. It stops when it can not code an angle, i.e. at a not explained angle.

    Parameters
    ----------
    seq : list
        Sequence of measured angles.
    permutation_block_max_size : int
        Maximum number of organs involved in permutations.
    canonical_angle : float
        Canonical divergence angle.
    borders_dict : dict
        Dictionary of theoretical angles and the corresponding intervals.
    kappa : float
        Concentration parameter.
    threshold : float
        A statistical threshold to prune the `n`-admissible trees.
    d_n : list
        List of theoretical angles.
    theoretical_coeff_angles : list
        List of coefficients of theoretical angles.
    d_n_coeff_dict : dict
        Dictionary of theoretical angles and their coefficients.
    seq_begin : bool, optional
        Indicates whether permutations at the beginning of the sequence should be taken into account. Default is ``False``.

    Returns
    -------
    phyllotaxis_analysis.n_admissible_tree.Nadmissibletree
        The constructed `n`-admissible tree.
    list
        The list of all intermediate trees generated during the process.
    """
    print("inside make_n_admissible_tree")

    all_trees_to_pickle = []

    sequence_length = len(seq)
    n_adm_tree = Nadmissibletree()
    index = 0  # the index of seq
    # Initializing tree

    # possible = list()
    candidates = list()
    horizon = min(sequence_length - index, permutation_block_max_size)
    for i in range(0, horizon):

        # all of the theoretical angles that can correspond to the measured subsequence
        possible = candidate_angles_list(seq[index: index + i + 1], borders_dict)  # remove 1 ?
        print("subsequence: ", seq[index: index + i + 1])

        if seq_begin:
            for pos in possible:
                nadm_list = is_n_admissible_first_angles(pos, permutation_block_max_size, canonical_angle,
                                                         d_n_coeff_dict)
                if nadm_list[0]:
                    candidates.append([pos, nadm_list[2]])
        else:
            for pos in possible:
                nadm_list = is_n_admissible_not_first_angles(pos, permutation_block_max_size, canonical_angle,
                                                             d_n_coeff_dict)
                if nadm_list[0]:
                    candidates.append([pos, nadm_list[2]])
    print("candidates : ", candidates)
    cand_len = len(candidates)
    not_subseq = np.ones(cand_len, dtype=bool)
    for j in range(cand_len):
        for i in range(j + 1, cand_len):
            if is_sub_sequence(candidates[j][0], candidates[i][0]):
                not_subseq[i] = False
                print("to remove : ", candidates[j][0], "is subsequence of ", candidates[i][0])

    print(not_subseq)
    new_candidates = [candidates[i] for i in range(cand_len) if not_subseq[i]]

    print("new_candidates : ", new_candidates)

    print(80 * "=")
    candidates_prob = list()
    for pos in new_candidates:
        for i in range(len(pos[0])):
            if prob_of_angle(seq[index + i], pos[0][i], kappa, d_n) < threshold:
                break
        else:
            candidates_prob.append(pos)
    print("candidates_prob :", candidates_prob)

    if candidates_prob == []:
        for affected_angle in max_prob_angle(seq[index], kappa, d_n):
            pb_angle = prob_of_angle(seq[index], affected_angle, kappa, d_n)
            log_pb_angle = log_of_prob_of_angle(seq[index], affected_angle, kappa, d_n)
            #            tree.addMarkedNodeProb(tree.treeRoot().id ,[ affected_angle, index, pb_angle, log_pb_angle, "null"])
            n_adm_tree.add_child(parent=n_adm_tree.root, value=affected_angle, level=index, pb_angle=pb_angle,
                                 log_pb_angle=log_pb_angle, question_mark=True, order_index="null")
            # question_mark == True marks unexplained angles
            all_trees_to_pickle.append(copy.deepcopy(n_adm_tree))
        return n_adm_tree, all_trees_to_pickle
    else:
        for pos in candidates_prob:
            tree_pos = list()  # tree_pos is a list of affected values and their level in the tree.
            current_id = n_adm_tree.root
            for i in range(len(pos[0])):
                pb_angle = prob_of_angle(seq[index + i], pos[0][i], kappa, d_n)
                log_pb_angle = log_of_prob_of_angle(seq[index + i], pos[0][i], kappa, d_n)
                order_index = pos[1][i + 1]
                tree_pos.append([pos[0][i], index + i, pb_angle, log_pb_angle, order_index])
            n_adm_tree.add_list_value_level_prob(n_adm_tree.root, tree_pos)

    print("vertices : ", n_adm_tree.vertices())
    print("properties : ", n_adm_tree.properties())

    # return n_adm_tree
    # Initializing tree finished
    #     print("Initializing tree finished")
    counter = 0
    while True:

        all_trees_to_pickle.append(copy.deepcopy(n_adm_tree))

        leaves_not_last = n_adm_tree.leaves_not_last(
            sequence_length - 1)  # The leaves whose level are less than sequence_length - 1

        if len(leaves_not_last) == 0:  # all leaves correspond to the last angle. The search is over.

            return n_adm_tree, all_trees_to_pickle

        leaves_not_last_not_question = [vertex for vertex in leaves_not_last if
                                        (not n_adm_tree.get_vertex_property(vertex)["question_mark"])]
        print("leaves_not_last_not_question", leaves_not_last_not_question)

        if leaves_not_last_not_question == []:  # all leaves whose levels are less than sequence_length - 1 have question marks
            leaves = n_adm_tree.leaves()
            if n_adm_tree.leaves_last(sequence_length - 1) == []:

                max_level = n_adm_tree.get_max_level()
                for feuille in leaves:

                    if n_adm_tree.get_vertex_property(feuille)["level"] != max_level:
                        ver = feuille
                        while n_adm_tree.is_leaf(ver):
                            pid = n_adm_tree.parent(ver)
                            n_adm_tree.remove_leaf(ver)
                            ver = pid
                    all_trees_to_pickle.append(copy.deepcopy(n_adm_tree))
                return n_adm_tree, all_trees_to_pickle
            else:
                for vertex in leaves_not_last:
                    ver = copy.deepcopy(vertex)
                    while n_adm_tree.is_leaf(ver):
                        pid = n_adm_tree.parent(ver)
                        n_adm_tree.remove_leaf(ver)
                        ver = pid
                    all_trees_to_pickle.append(copy.deepcopy(n_adm_tree))
                return n_adm_tree, all_trees_to_pickle

        for leaf in leaves_not_last_not_question:
            all_trees_to_pickle.append(copy.deepcopy(n_adm_tree))

            index = n_adm_tree.get_vertex_property(leaf)["level"] + 1
            pid = leaf
            possible = list()
            candidates = list()
            horizon = min(sequence_length - index, permutation_block_max_size)
            for i in range(0, horizon):
                possible = candidate_angles_list(seq[index: index + i + 1], borders_dict)
                if index <= permutation_block_max_size:
                    # if True:
                    # print("VIP 0", index, permutation_block_max_size)
                    for pos in possible:
                        path_root = n_adm_tree.path_to_root(leaf)
                        new_pos = list(path_root)
                        new_pos.extend(pos)
                        if seq_begin:  # here a voir
                            nadm_list = is_n_admissible_first_angles(new_pos, permutation_block_max_size,
                                                                     canonical_angle,
                                                                     d_n_coeff_dict)
                        else:
                            nadm_list = is_n_admissible_not_first_angles(new_pos, permutation_block_max_size,
                                                                         canonical_angle,
                                                                         d_n_coeff_dict)
                        if nadm_list[0]:
                            candidates.append([pos, nadm_list[2][len(path_root):]])
                else:
                    # elif False:
                    print("VIP 1", index, permutation_block_max_size)
                    for pos in possible:
                        new_pos = list(pos)
                        print("new_pos[0], 1", new_pos[0])
                        new_pos[0] = (new_pos[0] + (
                                n_adm_tree.get_vertex_property(leaf)["order_index"] - index) * canonical_angle) % 360
                        print("new_pos[0], 2", new_pos[0])
                        nadm_list = is_n_admissible_not_first_angles(new_pos, permutation_block_max_size,
                                                                     canonical_angle,
                                                                     d_n_coeff_dict)
                        if nadm_list[0]:
                            candidates.append([pos, [index + i for i in nadm_list[2]]])
            cand_len = len(candidates)
            not_subseq = np.ones(cand_len, dtype=bool)
            for j in range(cand_len):
                for i in range(j + 1, cand_len):
                    if is_sub_sequence(candidates[j][0], candidates[i][0]):
                        not_subseq[i] = False
            new_candidates = [candidates[i] for i in range(cand_len) if not_subseq[i]]
            candidates_prob = list()
            for pos in new_candidates:
                for i in range(len(pos[0])):
                    if prob_of_angle(seq[index + i], pos[0][i], kappa, d_n) < threshold:
                        break
                else:
                    candidates_prob.append(pos)
            if candidates_prob == []:
                for affected_angle in max_prob_angle(seq[index], kappa, d_n):
                    pb_angle = prob_of_angle(seq[index], affected_angle, kappa, d_n)
                    log_pb_angle = log_of_prob_of_angle(seq[index], affected_angle, kappa, d_n)
                    n_adm_tree.add_child(parent=pid, value=affected_angle, level=index, pb_angle=pb_angle,
                                         log_pb_angle=log_pb_angle, question_mark=True, order_index="null")
            else:
                for pos in candidates_prob:
                    tree_pos = list()  # tree_pos is a list of affected values and their level in the tree.
                    for i in range(len(pos[0])):
                        pb_angle = prob_of_angle(seq[index + i], pos[0][i], kappa, d_n)
                        log_pb_angle = log_of_prob_of_angle(seq[index + i], pos[0][i], kappa, d_n)
                        order_index = pos[1][i + 1]
                        tree_pos.append([pos[0][i], index + i, pb_angle, log_pb_angle, order_index])
                    n_adm_tree.add_list_value_level_prob(pid, tree_pos)