Cyril Fontlup

El-Ghazali Talbi

Keywords: Reconstruction Phylogénétique, Maximum Parcimonie, Optimisation Combinatoire, Recuit Simulé Phylogenetics, Maximum Parsimony, Combinatorial Optimization, Simulated Annealing

This work was carried out with the economic support of CONACyT

List of Tables

Bioinformatics is an interdisciplinary research area where intervene: biology, statistics, applied mathematics, chemistry, biochemistry, artificial intelligence and computational science. Bioinformatics studies the development of computational methods and statistical techniques to solve practical and theoretical problems derived from the storage, extraction, handling and distribution of information related with biological macromolecules such as DNA (Deoxyribonucleic acid), RNA (Ribonucleic acid) and proteins. These macromolecules are represented by their primary structure 1 . Hereafter the term "biological macromolecules" will be used to refer to DNA, RNA and protein sequences.

In Bioinformatics there are two important fields. The first one develops databases and computational tools. The second one uses these databases and tools to extract biological knowledge in order to understand the living systems. It is important to note that these fields are complementary and one can benefit from the other. Some of the main research areas in bioinformatics are:

-Sequence alignment: this refers to the comparison procedure of two or more sequences of macromolecules by looking for a series of individual characters or character patterns that are in the same order in the sequences to identify regions of similarity that may be a consequence of functional, structural, or evolutionary relationships between the sequences. Identical or similar characters are placed in the same column, and non-identical characters can either be placed in the same column as a mismatch or opposite to a gap 2 in one of the other sequences. In an alignment, non-identical characters and gaps are placed so as to bring as many identical or similar characters as possible into columns. -Protein structure prediction: is the prediction of the three-dimensional structure of a protein from its amino acid sequence i.e., the prediction of its secondary 3 , tertiary 4 , and quaternary structure from its primary structure. It is very important in bioinformatics and theoretical chemistry and it is critical in medicine (for example in drug design) and biotechnology (for example in the design of new enzymes). -Gene prediction: is the problem of locating genes in a genomic sequence 5 , i.e. identifying the regions of genomic 1. The primary structure of a protein is the linear sequence of its amino acid structural units, and partly comprises its overall biomolecular structure.

2. Gap: a special character used to indent columns and that serves as a spacer.

3. Secondary structure: is the local spatial arrangement of a polypeptide's backbone atoms regardless of the conformations of its side chains. 4. Tertiary structure: refers to the three-dimensional structure of an entire polypeptide. 5. Genomic sequence, is the entire nucleic acid sequence necessary for the synthesis of a functional protein (or functional RNA). It includes regions that contain no information.

12 CHAPTER 1. INTRODUCTION DNA that encode genes. It includes protein-coding genes as well as RNA genes, but may also include prediction of other functional elements such as regulatory regions. Gene finding is one of the first and most important steps in understanding the genome of a species once it has been sequenced. -Prediction of genic expression: is the process by which information from a gene is used in the synthesis of a functional gene product. These products are often proteins, but in non-protein coding genes such as ribosomal RNA (rRNA), transfer RNA (tRNA) or small nuclear RNA (snRNA) genes, the product is a functional RNA. The process of gene expression is used by eukaryotes and prokaryotes to generate the macromolecular machinery for life. Several steps in the gene expression process may be modulated, including the transcription, RNA splicing, translation, and posttranslational modification of a protein. -Protein-protein interactions: occurs when two or more proteins bind together, often to carry out their biological function. Many of the most important molecular processes in the cell such as DNA replication are carried out by large molecular machines that are built from a large number of protein components organized by their protein-protein interactions. As protein-protein interactions are so important there are a multitude of methods to detect them [START_REF] Phizicky | Protein-protein interactions: Methods for detection and analysis[END_REF]. -Evolution modeling: evolution is the change in the inherited characteristics of biological populations over successive generations. Evolutionary processes give rise to diversity at every level of biological organization, including species, individual organisms and molecules such as DNA and proteins. Evolution modeling tries to model the changes that occur in biological populations. -Reconstruction of phylogenetic trees: a phylogenetic tree or evolutionary tree is a branching diagram or tree showing the inferred evolutionary relationships among various biological species or other entities based upon similarities and differences in their physical and/or genetic characteristics. The sequences (also called taxa, singular taxon) joined together in the tree are considered to descend from a common ancestor. The phylogenetic reconstruction tries to model the evolutionary relationships between a group of species in structures of trees (rooted or unrooted) using their genetic information represented by sequences of characters.

In the remaining of this work, we will focus on phylogenetic reconstruction also called Phylogenetics.

There are two main approaches used to the phylogenetic reconstruction. The first one known as comparative morphology, is based on the analysis of anatomic, physiologic and behavioral characteristics shared between the studied organisms. The second one, called comparative biochemistry, studies the common characters (amino acids and nucleotides) of macromolecules that belong to a group of analyzed individuals.

Different methods have been designed for the reconstruction of phylogenetic trees, Maximum Parsimony (MP) is one of the most interesting in the field of combinatorial optimization. MP is a reconstruction method based on the assumption that the most probable hypothesis is the one that requires the fewest explanations. MP is used to classify the living beings with the purpose of recovering their phylogenetic history. According to Cladism 6 , the species should be classified according to monophyletic groups 7 [START_REF] Hennig | Phylogenetic Systematic. Phylogeny[END_REF]. MP tries to maximize the evolutionary similitudes between species, this implies identifying a tree topology with the minimum tree length, i.e. the topology that requires the smallest number of evolutionary changes (or transformation in characters state) in order to explain the differences between the sequences [START_REF] Kluge | Quantitative phyletics and the evolution of anurans[END_REF].

Motivation

There are two main reasons for the study of phylogenetic tree reconstruction. The first one concerns biology where it is used for the generation of new vaccines, antibacterial, herbicide molecules and in the study of the dynamic of microbial communities [START_REF] Pace | A molecular view of microbial diversity and the biosphere[END_REF]. For example, in the pharmaceutical industry, it helps to the development of new drugs. It is worth mentioning, that in molecular biology, MP is used in the prediction of protein structures, determination of homology of sequences and the classification of proteins [START_REF] Murakami | Sharp2: Protein-protein interaction predictions using patch analysis[END_REF].

The second reason is related to Computer Science where MP is a NP-complete problem [START_REF] Garey | The rectilinear steiner tree problem is np-complete[END_REF], due to its equivalence to the problem of combinatorial optimization known as The problem of the Steiner tree in a hypercube which represents an important challenge due to its nature and the fact that the solutions are not vectors of values but trees. 6. Cladism: biology area that defines the evolutionary relationships between organisms based on derived similitudes. 7. A Monophyletic group is composed by the common ancestor and all its descendants. Indeed, the size of the search space, i.e. the number of possible rooted tree topologies, for this problem is given by the next expression [START_REF] Xiong | Essential Bioinformatics[END_REF]:

(2n -3)! 2 n-2 (n -2)! (1.1)
where n is the number of studied species. It is easy to observe that this expression grows in factorial way according to the value of n, and generates huge search spaces.

The main goal of our work is to develop algorithms that provide solutions of good quality in a reasonable amount of time.

It is also worth mentioning that the cladistic community has contributed to the resolution of MP during several years and has used techniques that are well-known from the combinatorial optimization community. Those techniques were tailored to the resolution of MP.

State of the art 2.1 History

The phylogenetics analysis finds its origin from the XVIII century, when evolutionary theories arise and used the laws of inheritance as a fundamental basis. For this reason, an historical review with important facts will be shown.

Stephen Jay Gould argues that the French biologist Jean-Baptiste Lamarck (1744 -1829) was the primary evolutionary theorist and his ideas and the way in which he structured his theory is the base of the most subsequent thinking in evolutionary biology, through to the present day [START_REF] Gould | The Structure of Evolutionary Theory[END_REF]. Lamarck constructed one of the first theoretical frameworks of organic evolution which was rejected during his lifetime [START_REF] Burkhardt | The Spirit of System: Lamarck and Evolutionary Biology: Now with" Lamarck in 1995[END_REF].

Lamarck believed in the inheritance of acquired characteristics and the use and disuse model by which organisms developed their characteristics. He incorporated this belief into his theory of evolution giving rise to what we know as Lamarckism (or Lamarckian inheritance).

Lamarckism is the idea that an organism can pass on characteristics that it acquired during its lifetime to its offspring (also known as heritability of acquired characteristics or soft inheritance). Lamarck incorporated two ideas into his theory of evolution which in his days were considered to be generally true [START_REF] Packard | Lamarck, the founder of evolution. his life and work, with translations of his writings on organic evolution[END_REF]:

use and disuse: individuals lose characteristics they do not require (or use) and develop characteristics that are useful.

inheritance of acquired traits: individuals inherit the traits of their ancestors. When Charles Darwin (1809 -1882) published his theory of evolution by natural selection in On the Origin of Species (Darwin, 1859), he gave credence to what Lamarck called use and disuse inheritance, but rejected other aspects of Lamarck's theories. He proposed the mechanisms that had given origin to the whole biodiversity. His work gave rise to the Theory of Evolution (Darwin, 1872). The ideas of Darwin are generally summed up in the sentence survival of the fittest, i.e., those that have adaptation abilities to the environment will have more chance to survive, and then reproduce. Their offspring will inherit their characters. In a nutshell we can say that, from Darwin's point of view, all living species have been created from previous generations, given that all the living beings share the same genetic code. We could conclude that all the inhabitants in the planet are descendants of a common ancestor (Darwin, 1859).

Gregor Mendel (1822Mendel ( -1884) ) made some experiments in a small garden, which consisted of crossing plants that should differ in at least two characters. After some time of observations, he proposed the laws of inheritance. In the year 1865 he published his results with the title Experiments on Plant Hybridization [START_REF] Mendel | Experiments in plant hybridization (1865)[END_REF]. Nevertheless, his work had not been recognized until the year 1900 when the scientific community began to take an interest in the mechanisms of genetic transmission.

Ernst Haeckel (1834Haeckel ( -1919) ) in the year 1866 supported by the Darwin idea, radically defended The Recapitulation theory, also called the biogenetic law. It is a biological hypothesis that in developing from embryo to adult, animals go through stages resembling or representing successive stages in the evolution of their remote ancestors. This theory is currently discarded (Haeckel, 1866). He anticipated the fact that the key of the hereditary factors reside in the cell core. He was the first one to establish the phylogenetic hypothesis of the biologic biodiversity adjusted to the evolution theory.

The German entomologist Willi Hennig (1913Hennig ( -1976) ) gave birth to the cladistic analysis referring it as phylogenetic systematics [START_REF] Hennig | Phylogenetic Systematic. Phylogeny[END_REF]. The use of the terms Cladistics and clade was popularized by other researchers.

Cladistics is a classification approach in which items are grouped together based on whether or not they have one or more shared characteristics that come from the group's common ancestor and are not present in more distant ancestors. Therefore, members of the same group are thought to share a common history and are considered to be more closely related [START_REF] Pearsall | The New Oxford Dictionary of English[END_REF][START_REF] Lagassé | The Columbia Encyclopedia[END_REF][START_REF] Shaban-Nejad | Ontology-inferred phylogeny reconstruction for analyzing the evolutionary relationships between species: Ontological inference versus cladistics[END_REF].

From the time of his original formulation until the end of the 1980s cladistics remained a minority approach to both phylogenetics and taxonomy. However, in the 1990s it rapidly became the dominant set of methods of phylogenetics in evolutionary biology, because computers made it possible to process large quantities of data about organisms and their characteristics.

Cladism is a method used to reconstruct the genealogy of organisms and is based on three principles:

1. a sequence is called a taxon (plural taxa) and represents taxonomic units in the biological system of classification of organisms, 2. every group is a descendant of a common ancestor (i.e. are monophyletic),

3. the tree with the highest probability to be the correct is the one that proves to be the most simple and has the minimum number of changes and is defined as the most parsimonious tree.

The Maximum Parsimony (MP) problem, i.e. the search for a tree that is most parsimonious, arises from the necessity to classify the different species, in order to exhibit their evolutionary history. Phylogenetic trees represent the possibility of representation of the evolutionary history of a group of species under the parsimony criterion.

Please note that there are different parsimony optimality criteria known as Fitch, Wagner, Camin-Sokal, Dollo or weighted parsimony. Those criteria determine the number of changes of a substitution from one site to another [START_REF] Felsenstein | Evolutionary trees from dna sequences: a maximum likelihood approach[END_REF]:

-Fitch also called unweighted or unordered parsimony [START_REF] Fitch | Toward defining the course of evolution : Minimum change for a specific tree topology[END_REF]) is a generalized form of Wagner parsimony where any mutation has the same cost of 1, this is the parsimony criterion that we are using in this thesis, -Wagner or ordered parsimony [START_REF] Wagner | Problems in the classification of ferns[END_REF] imposes some order on the mutations: characters are allowed to reverse so that change from A to C costs the same number of steps as C to A, but also characters are additive so that if A to C costs 1, and C to G costs 1, then A to G must cost 2, -Dollo parsimony [START_REF] Farris | Phylogenetic analysis under Dollo's Law[END_REF] is based on the probability that a reversal is higher than a forward change, for example A to C is allowed with a higher probability than C to A, -Camin-Sokal parsimony [START_REF] Camin | A method for deducing branching sequences in phylogeny[END_REF] requires that all character evolution is irreversible: A to C is possible but not C to A, -finally generalized parsimony [START_REF] Swofford | Paup*: Phylogenetic analysis using parsimony 4.0 beta[END_REF]) assigns a cost to each transformation between states. For this work, we are only interested in Fitch parsimony for which all substitutions are given the same weight of one unit. [START_REF] Edwards | The reconstruction of evolution[END_REF] were the first ones to mention the maximum parsimony criterion in Phylogeny. They introduced computational methods to compute evolutionary trees based on genetic information. Their most important contribution was the introduction of stochastic methods of estimation applied to stochastic models of evolution.

We have analyzed the historical events that were the basis to this research problem. In the next section, the methods used to solve the phylogenetic reconstruction problem (or phylogeny) will be explained.

Methods to solve phylogenetic reconstruction

Since the 60's they all are dependent of a preliminary operation called a multiple alignment of sequences. However, as pointed out before, some software like Clustal [START_REF] Thompson | Multiple Sequence Alignment Using ClustalW and ClustalX[END_REF] use distance methods with the objective to determine the order in which the sequences will be added to compose the multiple alignment.

Many methods have been developed for the phylogenetic tree reconstruction. The principal approaches are:

Distance based or Clustering methods

They rely on a measure of distance between species [START_REF] Felsenstein | Inferring phylogenies[END_REF], generally the Hamming distance is used. It corresponds to the number of differences between two sequences. Distance methods start by generating a first matrix that describes the distances between each pair of sequences. Then an iterative procedure groups the two closest sequences and the matrix distance is reduced and recomputed.

The complexity of distance methods is low in the order of O(n 2 ) which makes them suitable to analyze huge sets of data.

Grouping based

-WPGMA (Weighted Pair Group Method with Averaging): is a clustering method that uses a sequential clustering algorithm, in which a tree is built in a stepwise manner, by grouping sequences or groups of sequences (taxa) that are most similar to each other; i.e. for which the evolutionary distance is smallest. When two taxa are grouped, they are treated as a new single taxon (see Figure 2.1) and from among the new group of taxa, it is necessary again to identify the pair for which the similarity is the highest, and so on, until only two taxa remain (see Figure 2.1, below). In WPGMA, the averaging of the distances is not based on the total number of taxa in the respective clusters. For example, when taxa A, B (which have been grouped before) and C are grouped into a new node 'u', then the distance from node 'u' to any other node 'k' (e.g. grouping D and E) is simply computed as follows:

d(u, k) = d(A, B) + d(C, k) 2 (2.1)
-UPGMA (Unweighted Pair Group Method Using Arithmetic Average): [START_REF] Sokal | A statistical method for evaluating systematic relationships[END_REF] proposed this iterative process. It is similar to WPGMA where the equation to compute the new distance is:

d(u, k) = |AB|d(A, B) + |C|d(C, k) |AB| + |C| (2.2)
where |C| is the number of taxa in C.

-Neighbor-Joining: is a clustering method for the creation of phylogenetic trees, created by [START_REF] Saitou | The neighbor-joining method: a new method for reconstructing phylogenetic trees[END_REF].

Usually used for trees based on DNA or protein sequence data, the algorithm requires knowledge of the distance between each pair of taxa to build the tree. In Figure 2.2(A) we can see the starting with a star tree; then a matrix Q is calculated in order to choose a pair of nodes for joining, in this case f and g. They are joined as a newly created node u (see Figure 2.2(B)). The part of the tree shown as dotted lines is now fixed and will not be changed in subsequent joining steps. The distances from node u to the nodes ae are computed. This process is then repeated, using a matrix of just the distances between the nodes, a, b, c, d, e and u, and a Q matrix derived from it. In this case u and e are joined to the newly created v, as shown in Figure 2.2(C). Two more iterations lead first to (D), and then to (E), at which point the algorithm is done, as the tree is fully resolved. [START_REF] Fitch | Construction of phylogenetic trees[END_REF] proposed a method that selects the best tree taking into account the minimum deviation between the distances computed for all the branches in the tree and the distances of the set of original data. -Minimum evolution: it employs linear computation to build trees. This method searches a tree with the minimum number of mutations. [START_REF] Rzhetsky | Theoretical foundation of the minimum-evolution method of phylogenetic inference[END_REF] applied this method for the reconstruction of phylogenetic trees.

-

Characters based methods

-Maximum parsimony: also called minimum-steps or minimum evolutionary steps method [START_REF] Felsenstein | Maximum Likelihood and Minimum-Steps Methods for Estimating Evolutionary Trees from Data on Discrete Characters[END_REF]) is a character-based approach that searches the tree with the minimum number of mutations (i.e. it needs to minimize the total number of evolutionary steps), which explains in a better way the evolutionary relationships between the studied taxa represented on the leaves. MP uses a matrix of discrete phylogenetic characters to infer one or more optimal phylogenetic trees for a set of taxa and operates by evaluating candidate phylogenetic trees according to an optimality criterion. Occam's razor, a principle of theoretical parsimony suggested by William of Ockham in the 1320's, asserted that it is vain to give an explanation which involves more assumptions than necessary. Under this principle the best phylogenetic tree is the one that contains the least evolutionary changes to explain observed data [START_REF] Fitch | Toward defining the course of evolution : Minimum change for a specific tree topology[END_REF].

Trees are scored (evaluated) by using a simple algorithm to determine how many "steps" (number of mutations) are required to explain the distribution of each character. A step represents a change from one character state to another. The algorithm does not assign particular character states to internal nodes (hypothetical ancestors) on a tree: the least number of steps can involve multiple, equally costly assignments and distributions of evolutionary transitions. MP optimizes the total number of changes. -Maximum likelihood: evaluates a hypothesis about evolutionary history in terms of the probability that the proposed model and the hypothesized history would give rise to the observed data set. The maximum likelihood method uses standard statistical techniques for inferring probability distributions to assign probabilities to particular possible phylogenetic trees. The method requires a substitution model to assess the probability of particular mutations; roughly, a tree that requires more mutations at interior nodes to explain the observed phylogeny will be assessed as having a lower probability. This is broadly similar to the maximum-parsimony method, but maximum likelihood allows additional statistical flexibility by permitting varying rates of evolution across both lineages and sites. Maximum likelihood is thus well suited to the analysis of distantly related sequences, but because it formally requires search of all possible combinations of tree topology and branch length, it is computationally expensive to perform on more than a few sequences. The "pruning" algorithm, a variant of dynamic programming, is often used to reduce the search space by efficiently calculating the likelihood of subtrees [START_REF] Felsenstein | Inferring phylogenies[END_REF]. The method calculates the likelihood for each site in a "linear" manner, starting at a node whose only descendants are leaves (that is, the tips of the tree) and working backwards toward the "bottom" node in nested sets. Searching tree topologies defined by likelihood has not been shown to be NP-complete [START_REF] Felsenstein | Inferring phylogenies[END_REF]), but remains extremely challenging because branch-and-bound search is not yet effective for trees represented in this way. This exhaustive method tries to search all the possible tree topologies and considers every position in the alignment (not only the informative sites). Foster and Hickey (1999) applied this method.

The advantages of Maximum Likelihood are the following:

-has often lower variance than other methods (i.e. it is frequently the estimation method least affected by sampling error) -tends to be robust to many violations of the assumptions in the evolutionary model -even with very short sequences, it tends to outperform alternative methods such as parsimony or distance methods -it uses all the sequence information The disadvantages are:

-very CPU intensive and thus extremely slow -the result is dependent of the model of evolution used

In the remaining of this work we will focus on Maximum Parsimony as it is an interesting optimization problem and the subject of this thesis.

Problem statement

Let S = {S 1 , S 2 , . . . , S n } be a set of n taxa 1 of length k that represent existing species and that are defined over an alphabet A . In the case of DNA A is composed of four nucleotides, the gap symboland ? which represents a missing or unknown symbol: A = {-, A, C, G, T, ?}. Generally, S is the result of a multiple alignment that in its turn relies on distance methods (see section 2.2) to define the order in which taxa are chosen (for example with Clustal). A phylogenetic binary tree T = (V, E) (see Figure 2.3) which represents ancestral relationships between species is composed by a set of nodes V = {v 1 , . . . , v r } and a set of edges E ⊆ V × V = {{u, v}|u, v ∈ V }. The set of nodes V of size (2n -1) is divided in two subsets:

1. taxa is the plural of taxon which stands for taxonomic unit, we can also use the term OTU which stands for Operational Taxonomic Unit.

-L is a subset of n leaves, i.e., nodes without descendants or OTUs, -I represents the n -1 internal nodes (also called hypothetical ancestors or HTUs 2 ), each one has two descendants. Note that the parsimony score is also called tree length by the cladistic community. In order to compute the overall score of a tree, Fitch's algorithm [START_REF] Fitch | Toward defining the course of evolution : Minimum change for a specific tree topology[END_REF]) starts from the root and goes down to the leaves and then gradually moves back from the leaves to the root in order to compute hypothetical ancestral taxa. This is often referred to as the first-pass of the algorithm. At each position of an internal node a set of bases is assigned. If two descendants x 2. HTU: Hypothetical Taxonomic Unit and y of an internal node z have some bases in common they are assigned to the internal node z i = x i ∩ y i for the i-th residue of the taxon. Otherwise all bases of both descendants are assigned to the parent z i = x i ∪ y i and a cost of one unit is added to the overall score of the tree (see Figure 2.4) as it corresponds to a mutation. The second-pass of the algorithm, which starts from the root and reaches the leaves, enables to assign one nucleotide for a site if many possibilities exist, in order to obtain a hypothetical tree. Only the first-pass is necessary to obtain the parsimony score. More formally, the parsimony sequence z of every internal node z = (x, y) ∈ I whose descendants are x and y and are represented by the sequences

A A G A G G C C G A C G A G +1 +1 +1 a) first-pass from leaves to root A A G G C A C A A b) assignment after second-pass
S x = {x 1 , • • • , x k } and S y = {y 1 , • • • , y k } is computed with the next expression: ∀i, 1 ≤ i ≤ k, z i = x i ∪ y i , if x i ∩ y i = ∅ x i ∩ y i , otherwise (2.3)
Consequently, the parsimony score of the sequence z (or number of mutations) is defined as:

φ(z) = k i=1 C i where C i = 1, if x i ∩ y i = ∅ 0, otherwise (2.4)
and the parsimony score of the tree T is computed with the next equation:

φ(T ) = ∀z∈I φ(z) (2.5)
i.e. the sum of all mutations.

In contrast, the large Maximum Parsimony problem, which is far more complex, consists in finding a tree topology T * where φ(T * ) is minimum, indeed it is a minimization problem: Definition 2 (Large parsimony problem or Maximum Parsimony problem) Given a set S of n sequences of length k find a most parsimonious tree T , i.e. a tree with minimum score.

φ(T * ) = min{φ(T ) : T ∈ T } (2.6)
where T is the set of all possible tree topologies or search space.

In the following we will use the Newick notation to display trees which is a standard for representing trees in computerreadable format using parenthesis and commas. For example, the trees on Figure 2.5 will be represented respectively by : (A, B), ((A, B), C) and(A, (B, C)).

A B A B C A B C Figure 2.

-Newick tree representation

As an example, we provide the graphics of figures 2.6, 2.7, 2.8. The initial tree as a parsimony score of 4 (green rectangle) and we put into red circles the local parsimony score of each internal node. The first move (degraph, regraph) (see Figure 2.7) enables to obtain a tree of score 3. The last move (see Figure 2.8) will generate a tree with an optimal parsimony score of 2.

Computational complexity

The cost of a tree can be computed in polynomial time (see [START_REF] Fitch | Toward defining the course of evolution : Minimum change for a specific tree topology[END_REF]). A rooted binary tree of n leaves has n -1 internal nodes, thus the complexity of the small parsimony problem is O(n × k). Indeed, we need to compute the hypothetical sequences of n -1 internal nodes. However, the search for an optimal tree is computationally intractable: the large parsimony problem is extremely difficult to solve since it is equivalent to the NP-complete Steiner problem in a hypercube [START_REF] Foulds | The Steiner tree problem in phylogeny is NP-complete[END_REF]. This is why, as we shall see later on, heuristics methods constitute the main alternative in order to obtain near-optimal trees with reasonable computation time [START_REF] Goloboff | Character optimization and calculation of tree lengths[END_REF][START_REF] Nixon | The parsimony ratchet, a new method for rapid parsimony analysis[END_REF].

Complexity analysis

The MP problem is NP-complete [START_REF] Gusfield | Algorithms on Strings, Trees, and Sequences: Computer Science and Computational Biology[END_REF]. This is easily proved because it is equivalent to the problem of Steiner tree in a hypercube, which is known in combinatorial optimization as a NP-complete problem [START_REF] Garey | The rectilinear steiner tree problem is np-complete[END_REF]. Definition 3 (Hypercube) A hypercube d is a undirected graph with nodes 2 d bijectively tagged, i.e., each tag belongs to only one element in the graph and all elements in the graph belong at least to one element of the group of tags, where the tags are represented by the integer numbers between 0 and 2 d -1. Two nodes in the hypercube are adjacent if and only if, the binary equivalent of their tags differs by only one bit.

In order to easily understand the reduction of the MP problem to the problem of Steiner tree in a hypercube, we will take as a basis the work of [START_REF] Goëffon | Nouvelles Heuristiques de Voisinage et Mémétiques pour le Probléme Maximum de Parcimonie[END_REF], where a phylogenetic tree could be seen as a set of points connected in a sequence inside a hypercube. We will consider two binary sequences S 1 and S 2 , with A = {0, 1} and 3 points in simultaneous way, where every sequence is composed by one or two characters (see the Now, imagine an hypercube of m dimensions, whose nodes are represented by each one of 2 m different binary sequences, which in accordance with the definition, are considered adjacent only if their Hamming distance is one.

It should be noticed that with two taxa, only one phylogenetic tree could be constructed. On Figure 2.9 we represent a tree in a hypercube, where the biggest circles correspond to the sequences of the problem and are tagged with each one of the 2 mn possible sequences. In this context the MP problem consists in finding the shortest tree in the hypercube that at the same time contains all nodes corresponding to the specific data. Every edge in the tree represents a state change, and the parsimony score of the tree is equal to its number of edges (length). The resulting tree is not always a binary tree, but it is easy to transform it. First the nodes of grade 3 are inserted where some nodes in the hypercube have a superior grade. Then, particular nodes are extracted from the hypercube to be transformed in leaves if their grade is 2. Finally, in order to complete the binary tree, it is sufficient to replace all the pathways (x i , . . . , x i+l ) where only x i and x i+l have a different grade of 2, by a simple edge (x i , x i+l ) tagged with the length l of the pathway. In this way, a solution without root is obtained and it is possible to insert a root in any part in the tree.

In Figure 2.9, the sample 1 represents a hypercube of dimension one, which only have a solution of the problem. In the sample 2, almost of the trivial topology of a tree with two leaves, it is represented the shortest tree that optimizes the result of parsimony. In the same way in the sample 3, there are two possible explanations: an hypothetical ancestor tagged with 01 and another one with 10. Definition 4 (Steiner tree.) Let G = (V, E) be an undirected graph. V is the set of nodes and E ⊆ V 2 are the edges. In each edge (v i , v j ) ∈ E a non negative weight is associated

ω ij ∈ N * . We consider X ⊆ V is a sub-set of nodes of V . A Steiner tree ST G (X) = (V ′ , E ′ ) is a tree such that X ⊆ V ′ ⊆ V and E ′ ⊆ E. The nodes V ′ \ X are Steiner points. The weight ST G (X),
that is denoted as: Ω(ST G (X)), is defined by:

Ω(ST G (X)) = i,j (vi,vj )∈E ′ ω ij (2.7)
The Steiner problem called in this way by its creator the Swiss mathematician Jakob Steiner (1796-1863), consists in finding the tree with the minimum length given G and X. The Steiner problem is generated when the edges of G do not have an associated value. In this case, ω ij = 1 for all i, j such that

(v i , v j ) ∈ E ′ .
When the MP problem is mapped to binary sequences it is equivalent to the problem of the Tree Steiner in a hypercube. A hypercube of dimension m has a register for each one of the 2 m string of m bits. Two sequences separated by a state change, are represented by two adjacent vertexes in the hypercube. It is easy to see how the Steiner tree is the shortest, and is equivalent to the most parsimonious tree, tolerating the possibility of a non binary tree. The number of edges in the Steiner tree is equal to the minimum parsimony score. [START_REF] Sankoff | Locating the vertices of a Steiner tree in an arbitrary metric space[END_REF] In order to solve the Steiner problem, we will consider an hypercube H = (V, E) of dimension 4 and a set of five nodes X = {0000, 0011, 0101, 0110, 1001}. In Figure 2.10.a H is represented in the five squares denoted by the dark circles of higher size, as well as a Steiner tree of length 6. This tree is showed in In practice there is no need to go beyond, because the tree in Figure 2.10.b, by the elimination of the unnecessary Steiner points (level 2), shows all the necessary information. However, to achieve greater rigor and satisfy the Definitions 3 and 4, the tree is represented most clearly in Figure 2.10.c. We could find a score of 6.

In accordance with [START_REF] Foulds | The Steiner tree problem in phylogeny is NP-complete[END_REF], the Steiner problem in a hypercube is NP-complete, therefore the MP problem is NP-complete too.

Analysis of search space

The size of the search space for the MP problem is given by the number of the possible tree topologies for the set of n analyzed species. The next expression gives an estimation of the size [START_REF] Xiong | Essential Bioinformatics[END_REF]:

N R = (2n -3)! 2 n-2 (n -2)! (2.8)
Where N R is the total number of possible tree topologies with root. For example, if we apply the last equations for a set of n = 6 species the result would be 945 different topologies which is relatively small. But, for an instance with n = 125 species, handled commonly in the literature, the total number of possible rooted tree topologies is N R = 1.62129384E+243.

In Figure 2.11 where we use a logarithmic scale we could see that the number of possible topologies of phylogenetic trees with root grows in function of the number n of species and the search space for this problem increases in factorial way in function of the value of n.

Resolution of MP

A thorough and historical introduction of methods to solve the Maximum Parsimony problem is available in Giribet (2007).

Exact methods

Exact methods explore the entire search space and provide the best solutions although they represent unfeasible methods to solve the problems due to the size of the search space which grows in an exponential way in function of the size of the instances. For the MP problem the size of an instance is given by the number of taxa studied n.

Exhaustive search

This method constructs all possible tree topologies then computes their parsimony score. Due to the prohibitive size of the search space (see Section 2.4.2), this method is impractical and may only be used to solve instances with less than ten taxa (Swofford, Olsen, Waddell, and Hillis, 1996a), from which all possible tree topologies are built.

Branch and bound algorithm

Branch and bound is a strategy to avoid exhaustive search. Theoretically, a branch and bound algorithm can not ensure polynomial time complexity in the worst case. [START_REF] Hendy | Branch and bound algorithms to determine minimal evolutionary trees[END_REF] were the first to implement a branch and bound algorithm which was designed to construct a minimum evolutionary tree of 11 species.

A first tree T 0 is generated (e.g. by a random process) and its score serves as an upper bound. Then a new tree T is built with 2 sequences. All remaining sequences are added iteratively on all possible branches creating new trees that are kept only if their score is lower than T 0 . If by adding a taxon on a branch a tree T i has a score greater than T 0 it is discarded. [START_REF] Wu | Approximation and exact algorithms for constructing minimum ultrametric trees from distance matrices[END_REF] used branch and bound to build trees with instances between 12 and 25 taxa. At first they build a n × n distance matrix M . Next the taxa are relabeled (1, 2, . . . , n) using a method to build an efficient permutation called MaxMin permutation. This method orders the taxa such that the taxon that has more characters in common with the others is in the first place and the last one is the taxon that has less characters in common.

Then they start to create the root v of the BBT (Branch and Bound Tree) such that v represents the only topology with leaves 1 and 2, and then the first tree is created using a modification of the algorithm UPGMA called UPGMM (Unweighted Pair Group Method with Maximum) and its score is stored in the variable U B. The building of the BBT then consists in adding one by one every taxon and if its score exceeds the value of U B all the nodes are deleted and they then build a new tree with the taxa that have not been used. They show that the efficiency of their algorithm with respect to the traditional method implemented by [START_REF] Hendy | Branch and bound algorithms to determine minimal evolutionary trees[END_REF] depends on the use of the MaxMin permutation.

Nevertheless, the fact that MP is a highly combinatorial problem makes exact methods impractical when the number of taxa exceeds 25 taxa for efficiency and storage reasons.

Heuristic methods

Heuristic methods try to generate optimal or near-optimal solutions to the problem in an iterative way. These methods are used when exact methods are unable to provide a solution in a reasonable amount of time. Although, the heuristic methods are generally more efficient than exact methods, they do not guarantee that generated solutions are optimal.

The first heuristic method of tree construction is the algorithm proposed by [START_REF] Wagner | Problems in the classification of ferns[END_REF] and implemented by Farris in 1970. Such trees were originally calculated by hand and used as final results to interpret a phylogeny. But it became evident that in the presence of homoplasy 3 Wagner trees were suboptimal solutions.

Greedy algorithms

Greedy algorithms are used in phylogenetic reconstruction under the name stepwise addition (Swofford, Olsen, Waddell, and Hillis, 1996b), Random Addition Sequence (RAS) [START_REF] Goloboff | [END_REF] or Wagner trees. They use an optimization criterion: a taxon is inserted on a branch that minimizes the parsimony score of the resulting tree. It is possible that there are many positions with a minimum score; generally one is chosen randomly. [START_REF] Andreatta | Heuristics for the phylogeny problem[END_REF] proposed a comparison between three greedy algorithms of different complexity and efficiency: -1stRotuGbr: selects a taxon randomly at every iteration and inserts it in the position that minimizes the parsimony score of the tree. -Gstep_wR: selects a taxon and a branch such that when the taxon is added in that position, the increase in the parsimony score is minimum. -GRstep: is a variant of Gstep_wR where the selected branch and taxon do not increase the parsimony score more than ten From their experiments they concluded that Gstep_wR is slightly more efficient than 1stRotuGbr but its complexity (O(n 3 ) vs O(n 2 )) increases dramatically the runtime. GRstep is useful only when it is combined with local search. [START_REF] Weber | Representation in stochastic search for phylogenetic tree reconstruction[END_REF] explored the representation based on a greedy decoder, a parameterization of the greedy variant of stepwise addition. They consider that any permutation of taxa determines a tree, by execution of a greedy procedure, the permutation itself can be taken to be a representation of that tree. The score of a permutation is calculated by greedily decoding it into the corresponding tree and calculating the score of that tree. The random variant of the stepwise addition constructs a tree under the direct representation. Random trees under the greedy decoder representation can be generated by decoding randomly selected 3. homoplasy: correspondence between parts or organs acquired as the result of parallel evolution or convergence permutation of the taxa. They compared both greedy and randomly built permutation to construct a greedy tree and they get a considerable difference between the two methods getting the best tree scores with the greedy decoder.

Stochastic local search (LS)

In Combinatorial Optimization, an instance of a problem consists of a set of feasible solutions and a cost function used to evaluate the quality of solutions. The cost function is also called objective function, score function or fitness measure. The goal is to find a solution with the optimal cost among all feasible solutions. Generally the problems addressed are computationally intractable, thus approximation algorithms have to be used. One class of approximation algorithms that have been surprisingly successful in spite of their simplicity are local search methods [START_REF] Bibliography | Local search in Combinatorial Optimization[END_REF][START_REF] Ch | Combinatorial Optimization: Algorithms and Complexity[END_REF].

The MP problem has a huge search space; for that reason it was empirically proved that the methods of stochastic local search (or neighborhood search) are suitable to solve this problem [START_REF] Ganapathy | On contract-and-refine transformations between phylogenetic trees[END_REF]. This kind of methods starts with an initial solution s 0 , that is improved through an iterative search process. This process continues until a stop condition is satisfied, which could be the number of solutions examined.

Local search is based on the concept of neighborhood. A neighborhood N of a solution s is a set of solutions that are in some sense close to s, for example because they can be easily computed from s or because they share a significant amount of structure with s. We will say that s ′ is a neighbor of s by using the following notation: s ′ ∈ N (s).

Definition 5 (Local optimum)

A local optimum of a combinatorial optimization problem is a solution that is optimal (either maximal or minimal) within a neighboring set of solutions. This is in contrast to a global optimum.

Definition 6 (Global optimum) A global optimum is the optimal solution among all possible solutions.

In this sense the neighborhood has much influence on a search as it may, or may not, be able to generate the optimal solution.

Descent or branch-swapping.

The descent algorithm (also called hill-climbing in the case of maximization (see Algorithm 1) generates an initial solution s 0 , then searches a close solution s ′ in the neighborhood of the current solution s, which should be better than s for the cost function. This process is repeated until there is no more improvement of s. In the case of MP, the descent is also called a round of branch-swapping by biologists from the Clasdistics community. A neighborhood function is then a branch-swapping for them.

Algorithm 1: Descent algorithm 2005) made a comparison between different methods to construct phylogenetic trees using a permutation either random or greedy for the initial solution, then they compared their algorithm called Greedy High Climbing (GHC) using a new perturbation method that exchanges the positions of two pairs of leaves. They compared their results against the efficiency of the TBR (see next Section) neighborhood function and showed that this new method of perturbation is not always more efficient than TBR.

1 Descent(s 0 , N ) input: s 0 : initial solution, N : a neighborhood output: best solution found 2 s ← s 0 3 iteration ← 1 4 while ∃s ′ ∈ N (s), s.t. f itness(s ′ ) < f itness(s) do 5 s ← s ′ 6 iteration ← iteration + 1 7 return s Weber et al. (
The descent algorithms are the base of many methods that efficiently solve the MP problem.

Neighborhood functions.

We now come to the description of the main neighborhood functions found in the literature relative to trees:

Nearest neighbor interchange (NNI) NNI is a technique proposed by [START_REF] Waterman | On the similarity of dendrograms[END_REF] which consists in exchanging two branches separated by an internal node. The size of the neighborhood is O(n). Each tree has 2n -6 NNI neighbors, n -3 internal branches and two possible movements by branch [START_REF] Robinson | Comparison of labeled trees with valency three[END_REF]. Figure 2.12 shows how NNI works. From an initial tree ((A, B), (C, D)) we can generate two neighbors ((A, C), (B, D)) and ((A, D), (B, C)). NNI can be considered as a small size neighborhood. TBR (Swofford and Olsen, 1990) consists in breaking the tree in two subtrees which will be reconnected from one of their branches. It can be considered as a double SPR where all possibilities of reconnection are tried: t is re-inserted on Tt and Tt is reinserted on t. From a given tree, the TBR neighborhood induces at most (2n -3)(n -3) 2 neighbor trees or O(n 3 ) (Allen and Steel, 2001) and is a large size neighborhood. See Figure 2.14.

The following property holds for the neighborhoods: N N I ⊆ SP R ⊆ T BR. A more thorough study shows that NNI is too simple to generate some topologies and will let the descent algorithm get stuck quickly. TBR as a large size neighborhood can have more influence and will sometimes enable to find high quality solutions. SPR represents a trade-off between efficiency and quality.

Single step (STEP)

STEP (Swofford and Olsen, 1990;[START_REF] Andreatta | Heuristics for the phylogeny problem[END_REF]) takes out a taxon (i.e. a leaf) from the current solution and puts it back into another branch of the tree (see Figure 2.15). Since each taxon may be reinserted into 2(n -1) -3 -1 different branches; each solution has 2n(n -3) neighbors. According to [START_REF] Mladenović | Variable neighborhood search[END_REF], VNS is a metaheuristic which systematically performs the procedure of neighborhood change, both in descent to local minimum and in escape from the valleys which contain them.

The variable neighborhood descent algorithm proposed by Ribeiro and Vianna (2005) (see Algorithm 2) follows this model to solve the MP problem, it has a good behavior but its runtime is important. [START_REF] Ganapathy | On contract-and-refine transformations between phylogenetic trees[END_REF] implemented this algorithm to solve the MP problem. They describe a new tree-rearrangement operation called p-ECR move, for p-Edge-Contract-and-Refine. Their algorithm computes the best 2-ECR neighbors of a given tree, based upon a simple data structure which also allows to calculate the best neighbor under NNI, SPR and TBR. They shows CHAPTER 2. STATE OF THE ART

Algorithm 2: Variable Neighborhood Search Algorithm 1 VNSDescent(s 0 , {N 1 , N 2 , . . . , N l }) input: s 0 : initial solution, {N 1 , N 2 , . . . , N l } : set of neighborhoods output: best solution found s ⋆ 2 s ← s 0 3 s ⋆ ← s 4 i ← 1 5 while i l do 6 s ′ ← Descent(s, N i ) 7 if f itness(s ′ ) < f itness(s ⋆ ) then 8 s ⋆ ← s ′ 9 s ← s ′ 10 i ← i + 1 11 return s ⋆
that the use of 2-ECR in conjunction with TBR and/or NNI moves, may be a more effective technique for exploring tree space than TBR alone. The algorithm uses a preprocessing step, in which they assign three labels to each node in the tree in order to assign three optimal state-assignment labels to each internal node.

Progressive neighborhood

Based on the observation that important topological modifications of the tree are only performed at the beginning of the descent, [START_REF] Goëffon | Progressive tree neighborhood applied to the maximum parsimony problem[END_REF] have proposed a Progressive Neighborhood (PN), which contrary to VNS starts with a medium size neighborhood (SPR) and is iteratively reduced to NNI. Results show that PN needs a smaller number of iterations than traditional SPR searches to obtain solutions of the same quality by limiting the evaluation of relevant configurations. Recently PN has been used by [START_REF] Pirkwieser | Finding consensus trees by evolutionary, variable neighborhood search, and hybrid algorithms[END_REF] to find consensus trees of high quality.

In order to evolve neighborhoods, a topological distance on trees is defined in [START_REF] Goëffon | Progressive tree neighborhood applied to the maximum parsimony problem[END_REF]) that enables to build a distance matrix for a set of taxa given a tree topology. This distance is also used to control the size of the neighborhood (i.e. the distance between a pruned edge and its inserted edge is at most equal to a given limit).

Definition 7 (Topological distance) Let i and j be two taxa of a tree T . The topological distance δ T (i, j) between i and j is defined as the number of edges of the path between parents of i and j, minus 1 if the path contains the root of the tree.

B C D E F A B C D E 0 1 1 3 3 2 3 3 2 0 2 2 1 1 1 δ T A B f k Tree T D E F j h g C Distance Figure 2.16 -Example of topological distance δ T
For example, on Figure 2.16, A and B have the same parent f , so δ T (A, B) = 0, and δ T (A, D) = 3, because the number of edges between f and g is 4 (f ↔ k ↔ j ↔ h ↔ g), and as we pass through the root node k, we decrease the value by one unit. Note that for the topological distance, we consider trees as unrooted, this is why we remove one unit when passing through the root node. Progressive Neighborhood based on the topological distance was implemented in the software Hydra [START_REF] Goëffon | Progressive tree neighborhood applied to the maximum parsimony problem[END_REF]. The process used in Hydra to reduce the size of the neighborhood takes into account a parameter M which corresponds to a maximum number of Local Search iterations. A parameter d is introduced to control the size of the neighborhood and is defined as the maximal distance between a pruned edge and the edge where it is reinserted (i.e. distance δ between their two descendant nodes). As such, changing d leads to neighborhoods of different sizes which are explored with a descent algorithm.

Iterated local search (ILS)

There are different ways to obtain a new solution of a neighborhood when the descent algorithm is stuck in a local optimum. Some algorithms generate a new initial solution and explore new neighborhoods while other such as Iterated local search (ILS) [START_REF] Charon | The noising methods: a survey[END_REF] disturb the optimum before to restart the descent [START_REF] Lourenco | Iterated local search[END_REF].

The ILS metaheuristic iteratively constructs a sequence of solutions by embedding a local search. ILS has four essential elements: initial solution, local search, perturbation and acceptance criterion of the local search (see Algorithm 3).

Algorithm 3:

Iterated Local Search algorithm for MP 1 IteratedLocalSearch(s 0 , N , stop_condition) input: s 0 : initial solution, N : a neighborhood output: best solution found s ⋆ 2 s ← s 0 3 while not(stop_condition) do 4 s ′ ← Descent(s, N ) 5 if f itness(s ′ ) < f itness(s ⋆ ) then 6 s ⋆ ← s ′ 7 s ← perturbation(s ′ ) 8 return s ⋆
ILS starts from a solution s 0 and performs a first descent until a local optimum s * is found. Then it applies a perturbation method in order to generate an intermediate state. The next step consists in applying local search to the new solution. The algorithm iterates until a stop condition is satisfied.

Vázquez-Ortiz and Rodríguez-Tello (2011) applied the ILS algorithm to the MP problem and compared it against a simulated annealing (SA) algorithm. Their experiments showed that the SA provides best results than ILS.

Simulated annealing (SA)

Simulated Annealing is a metaheuristic based on the work of [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF] in the field of statistical thermodynamic, where the process of annealing is modeled by the simulation of the energy changes in a system of particles. The temperature decreases until it converges to a freezing state. [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF] and [START_REF] Cerny | A thermodynamic approach to the traveling salesman problem: An efficient simulation algorithm[END_REF] worked independently and showed how this process could be applied to optimization problems, associating key concepts of the simulation process with elements of combinatorial optimization.

Algorithm 4 shows the basic method of simulated annealing for minimization problems with a Markov Chain to sample solutions iteratively and stochastically.

Parameter t represents the temperature. A solution with an increase ∆f in the cost function will be accepted with probability e (-∆f /t) .

In phylogenetic reconstruction this metaheuristic was applied by [START_REF] Barker | LVB: Parsimony and simulated annealing in the search for phylogenetic trees[END_REF] in his software named LVB 4 . It starts with a tree of random topology improved by a descent using alternately NNI and SPR (Swofford and Olsen, 1990). The algorithm uses the temperature as an arbitrary control parameter varying between 0 and 1.

During the search, changes that do not improve the tree score are accepted with probability p = e -(1/T )∆f , where T is the current temperature and ∆f is the change (increase of the score).

Vázquez-Ortiz and Rodríguez-Tello (2011) have implemented a simulated annealing algorithm analyzing carefully the different possibilities for the key components in order to find the combination that could offer the best quality solution to the problem at a reasonable computational effort. They analyzed two initialization methods (greedy and random), different combinations of neighborhoods, cooling schedules and stop conditions. They compared the best combinations of parameters with the best results 4. LVB: Ludwig Van Beethoven, because the author was listening to Beethoven when he was implementing the software.

Algorithm 4: Simulated Annealing algorithm 1 SimulatedAnnealing(s 0 , N , T 0 , T f , α) input: N : neighborhood, f : fitness function, M CL : Markov Chain length, α : cooling scheme, T i : initial temperature, T f : final temperature output: best solution found s ⋆ 2 s 0 ← GenerateInitialSolution() Another so-called implementation of SA is called Tree-Drifting. It is described in [START_REF] Goloboff | Techniques for analyzing large data sets[END_REF] as a TBR descent able to accept suboptimal trees. The key component of the method is the function designed to determine the probability of acceptance, which is based on both the absolute step difference and a measure of character conflict: the relative fit difference, which is the ratio of steps gained and saved in all characters, between the two trees. However there is no algorithm given for this method in the literature. Tree-Drifting is embedded in the software TNT (Goloboff, Farris, and Nixon, 2008a), but used alone, it does not provide results of good quality.

3 s ← s 0 4 s * ← s 0 5 t ← T i 6 while t > T f do 7 i ← 0 8 while i < M CL do 9 s ′ ← GenerateN eighbor(s, N ) 10 ∆f ← f (s ′ ) -f (s) 11 generate a random u ∈ [0, 1] 12 if (∆f < 0) or (e -∆f /t > u) then 13 s ← s ′ 14 if f (s ′ ) < f (s * ) then 15 s ⋆ ← s ′ 16 i ← i + 1 17 t ← αt

Tabu search (TS)

Tabu search is a local search metaheuristic algorithm created by [START_REF] Glover | The General Employee Scheduling Problem. An Integration of MS and AI[END_REF] and formalized in [START_REF] Glover | Tabu searchâpart i[END_REF] and [START_REF] Glover | Tabu Search, volume 3 of Handbook of Combinatorial Optimization[END_REF].

Tabu search uses local search or neighborhood functions to move from one potential solution s 0 to an improved solution s. To avoid being stuck in poor-scoring areas and explore regions of the search space that would be left unexplored by other local search, TS uses memory structures known as tabu list, a set of rules and banned solutions used to filter which solutions will be admitted to be explored by the search. A tabu list is a short-term set of the solutions that have been visited in the recent past (see Algorithm 5). the number of previous solutions to be stored). [START_REF] Lin | A tabu search algorithm for maximum parsimony phylogeny inference[END_REF] designed TABUPARS, that seems to be the only work in the literature to solve the MP problem using TS.

For their experiments they used a data set obtained from nuclear ribosomal DNA sequences (instances of size between 10 and 20, which is relatively small). In order to prevent the search process from revisiting a tree that has been visited recently, they use a tabu list which is composed of two arrays L and N . L contains the leaves of the tree ordered according to their position in the tree (from left to right) and N contains the internal nodes ordered according to their positions. As neighborhood function they use the method of leaf swapping (exchange two leaves) and profile change (exchange two nodes).

5. SAMPARS: Simulated Annealing for Maximum PARSimony.

⋆ 2 s ← s 0 3 s ⋆ ← s 4 i ← 1 5 T abuList ← ∅ 6 while not(stop_condition) do 7 s ′ ← Descent(s, N ) 8 if s ′ ∈ T abuList then 9 if f itness(s ′ ) < f itness(s ⋆ ) then s ⋆ ← s s ← s ′ T abuList ← add(T abuList, s) return s ⋆ 2.5.2.

Greedy randomized adaptive search procedure (GRASP)

This procedure is divided in two phases. During the first one; a greedy algorithm is applied where a taxon is added step by step in the position that reduces the increase of the parsimony score. The hybrid method of addition suggests an exchange of branches every time that a taxon is added to the tree, in order to do a correction that improves the parsimony score.

During the second phase, local search is applied (see Algorithm 6). [START_REF] Ribeiro | A GRASP/VND heuristic for the phylogeny problem using a new neighborhood structure[END_REF] applied the GRASP+VND algorithm to solve the MP problem, using some instances taken from the literature and other randomly generated for their tests, proving that the new heuristic showed best performance in the solution quality with respect to the algorithms reported in the literature at that time.

Algorithm 6: GRASP Algorithm 1 GRASP(max_iterations, N ) input: max_iterations : integer, N : neighborhood output: best solution found s ⋆ 2 s ⋆ ← ∅ // f itness(∅) = +∞ 3 i ← 1 4 while i max_iterations do 5 s 0 ← greedy_random_generation()

6 s ′ ← Descent(s 0 , N ) 7 if f itness(s ′ ) < f itness(s ⋆ ) then 8 s ⋆ ← s ′ 9 i ← i + 1 return s ⋆ 2.5.2.

Genetic and memetic algorithms

The concept of genetic algorithms (GA) was introduced by Holland (1975) and popularized by [START_REF] Goldberg | Genetic Algorithms in Search, Optimization and Machine Learning[END_REF]. Genetic algorithms belong to the class of evolutionary algorithms (EA), which generate solutions to optimization problems using techniques inspired by natural evolution, such as inheritance, mutation, selection, and crossover. They are modeled loosely on the principles of the evolution via natural selection. They rely on a population of individuals (candidate solutions) that undergo selection in the presence of variation operators such as mutation and recombination (crossover). A fitness function is used to evaluate individuals, and reproductive success varies with fitness (see Algorithm 7).

The evolution usually starts from a population of randomly generated individuals. In each iteration (generation), the fitness of every individual in the population is evaluated. The fittest individuals are stochastically selected from the current population, and each individual is modified (recombined and possibly randomly mutated) to form a new generation. The new generation of candidate solutions is then used in the next iteration of the algorithm. Commonly, the algorithm terminates when either a maximum number of generations has been produced, or a satisfactory fitness level has been reached for the population.

According to [START_REF] Hoos | Stochastic Local Search, Foundations and Applications[END_REF] in most cases a genetic algorithm is not effective enough, because the crossover and mutation operators do not increase efficiency of the search.

For this reason GA have been hybridized with local search, in which a local search method replaces the mutation operator or is applied after it. This drives the search in different areas while guiding it with the crossover and selection mechanisms. This kind of hybridized methods are known as memetic algorithms [START_REF] Moscato | Memetic algorithms: A short introduction[END_REF] Memetic algorithms have had a few applications in phylogenetic reconstruction. To date only the works of [START_REF] Matsuda | Protein phylogenetic inference using maximum likelihood with a genetic algorithm[END_REF] and [START_REF] Lewis | A genetic algorithm for maximum likelihood phylogeny inference using nucleotide sequence data[END_REF] have been registered.

The most recent application of this algorithm to solve the MP problem was performed by Richer et al. ( 2009), with the software Hydra which is based on an integration of an effective local search operator with a specific topological tree crossover operator.

The crossover operator is based on a topological distance defined as follows: given i and j two taxa of a Tree T , the topological distance δT (i, j) is defined as the number of edges of the path between parents of i and j, minus 1 if the path contains the root of the tree. This topological distance tries to keep the properties of the parent trees. The local search improves the quality of each created offspring by a descent algorithm using a neighborhood called Progressive Neighborhood (PN). The initial population is generated with the greedy algorithm 1stRotuGbr. For the crossover they use the Distance-Based Information Preservation tree crossover (DiBIPX) which takes into account the topological distance of parents trees, the key idea is to preserve the topological distance of parent trees. Experimentations on a number of real benchmark instances from TreeBase show that Hydra competes very well when compared to TNT. Hydra is able to find phylogenetic trees of better parsimony score with much fewer evaluations of candidate trees. Results show that Hydra improved the results of TNT on 6 instances and matched TNT on 5 instances and got worse in only one.

Parallel algorithms

A parallel algorithm or concurrent algorithm is an algorithm which can be executed a piece at a time on many different processing devices, and combined together again at the end to get the correct result [START_REF] Blelloch | Parallel algorithms[END_REF].

The cost or complexity of serial algorithms is estimated in terms of the space (memory) and time (processor cycles) that they take. Parallel algorithms need to optimize an additional resource: the communication between different processors. There are two ways parallel processors communicate, shared memory or message passing.

In order to solve the MP problem in efficient way, some authors [START_REF] Snell | Parallel phylogenetic inference[END_REF] and [START_REF] Du | Reconstruction of large phylogenetic trees: a parallel approach[END_REF] considered the parallel algorithms.

- [START_REF] Blazewicz | Adaptive memory programming: Local search parallel algorithms for phylogenetic tree construction[END_REF] proposed a parallel adaptive neighborhood search method. The algorithm achieves super-linear speedup and find solutions of good quality. Their main goal is to develop a strategy of cooperation of distributed computational nodes using an adaptive memory heuristic. In this paper, parallel adaptive memory programming (AMP) algorithms for parsimonious phylogenetic tree construction are proposed. The term adaptive memory programming was introduced by Glover and, it was originally related to tabu search [START_REF] Glover | Tabu searchâpart i[END_REF][START_REF] Glover | Tabu Search, volume 3 of Handbook of Combinatorial Optimization[END_REF]. Three versions of the parallel local search algorithms were created. Each one of the versions was combined with three neighborhood (i.e. NNI, SPR and TBR).The solution space is explored by the local search procedure with an adaptive memory mechanism. The algorithms try to find a local optimum using NNI. When found, such an optimum becomes a starting point for the local search procedure combined with the SPR. Next, the solution becomes a starting point for the local search procedure combined with TBR. The parallel algorithms have the master-slave structure. At the beginning of computation, the master process generates 30 starting solutions, which are then sent on the request to slave processes. When a local optimum has been found by a slave, it is sent to the master process, which collects the local optima that are used to create new starting solutions for slave process.

The motivation for such strategy is the hypothesis that local optima can contain fragments of the global optimum. They considered two stop criteria. The first was a value of the solution obtained while the second was based on a runtime limit. They compare their results with [START_REF] Ribeiro | A GRASP/VND heuristic for the phylogeny problem using a new neighborhood structure[END_REF]

Quartets concept

The meaning of this concept applied to phylogenetic trees is to break a big problem in different small problems and to solve every one independently. At the end we join all the partial solutions in order to get only one solution of the whole problem. In this method a phylogenetic tree is divided in quartets (topology on any four labels from the tree S) (see Figure 2.17). Then for each quartet we could notice 3 different topologies, where the topology that reduces the parsimony score is selected. The key idea is to consider small subsets of taxa, one at a time, and infer the phylogenies for these subsets. Given a complete set of correct quartets, a tree can be constructed in polynomial time.

Ben-Dor et al. (1998) use as input a list of weighted quartets over n taxa. Each quartet is a subtree on four taxa, and its weight represents a confidence level for the specific topology. The goal is to construct a binary tree with n leaves such that the total weight of the satisfied quartets is maximized. The first approach is based on geometric ideas, they embed the n points on the ndimensional unit sphere, while maximizing an objective function, which depends of the Euclidean distance between the four points and reflects the quartet topology. Given the embedding, a binary tree is constructed by performing geometric clustering, similar to neighbor joining, with the difference that the update phase retains geometric meaning: When two neighbors are joined together, their common ancestor is taken to be the center of mass of the original points. The second is based on dynamic programming, and it is guaranteed to find the optimal tree (with respect to the given quartets). They implemented both algorithms and ran them on real data for n = 15 taxa.

Reconfigurable computing

Reconfigurable computing is a computer architecture combining some of the flexibilities of software with the high performance of hardware by processing with very flexible high speed computing fabrics like field-programmable gate arrays (FPGAs). The principal difference when compared to ordinary microprocessors is the ability to make substantial changes to the data-path itself in addition to the control flow [START_REF] Estrin | Reconfigurable computer origins: the ucla fixed-plus-variable (f+ v) structure computer[END_REF].

Reconfigurable computing has been used to solve the MP problem by [START_REF] Kasap | High performance phylogenetic analysis with maximum parsimony on reconfigurable hardware. Very Large Scale Integration (VLSI) Systems[END_REF] who used Maxwell, a FPGA based supercomputer developed by FPGA High Performance Computing Alliance (FHPCA) in Scotland. They use the Sankoff's algorithm to calculate the tree length. They take into account the levels of the tree; the first level is the closest to the leaves. The sub-trees in each level are computed in parallel until the root is reached. This is the first FPGA implementation but it supports a maximum of 12 taxa.

Solving MP as a multi-objective problem

New ideas have emerged in order to solve the MP problem under a multi-objective concept where one tries to optimize both the Maximum Parsimony and Maximum Likelihood versions of the problem.

As explained in the Section 2.2 different criteria have been employed to evaluate possible solutions to get the best phylogenetic tree in order to guide a search algorithm towards the best tree. However, these criteria may lead to distinct phylogenies, which are often conflicting among them. In this context [START_REF] Cancino | A multi-objective evolutionary approach for phylogenetic inference[END_REF] proposed a multi-objective evolutionary algorithm, named phyloMOEA, which employs the maximum parsimony and likelihood criteria to evaluate solutions. They consider that this approach can be useful since it could produce a spectrum of equally optimal trees (Pareto front) according to all criteria. They tested phyloMOEA using 4 datasets. The algorithm finds for all datasets a Pareto front representing a trade-off between the criteria. The main objective of phyloMOEA is the computation of a set of non-dominated solutions (trees), which represents a trade-off between parsimony and likelihood scores.

phyloMOEA is a multi-objective evolutionary algorithm that uses the following parameters:

-Initial population: creates random trees -Fitness Evaluation: the parsimony and likelihood scores are calculated using Fitch and Felsenstein algorithms, respectively. The rank is calculated using a non-dominated sorting algorithm applied to a population for all generations -Selection: tournament selection, picks two individuals at random and choose the best one -Crossover operator: prune a subtree s T from parent T 1 and remove all leaves of s T from T 2. The offspring T ′ 1 results from a regraph of s T in T 2 -Mutation Operator: NNI At each iteration i, the algorithm divides R = P i ∪ Q i (where P is an initial population and Q the current population) into several frontiers, denoted by F 1 , F 2 , ..., F j . The first frontier (F 1 ) is formed by non-dominated solutions from R. Thus the second frontier F 2 is by non-dominated solutions from R -F 1 . This process is repeated to R -F 1 -F 2 , and so on, until R is empty. At the end of a PhyloMOEA execution, duplicate trees are removed from the final population. Finally, the Pareto optimal solutions are calculated.

Unfortunately phyloMOEA uses only a set of 4 instances of different size including the instance called zilla (500 taxa of 759 residues).

Other techniques

Some other techniques have been developed to tackle the problem of Maximum Parsimony. They were integrated in the sofware TNT (Goloboff, Farris, and Nixon, 2008a).

-Ratchet [START_REF] Nixon | The parsimony ratchet, a new method for rapid parsimony analysis[END_REF]: it consists in temporarily modifying the problem by selecting some sites or by changing the weight of some sites. Consequently the objective function is modified and it helps escaping from a local optimum, -Tree fusing (TF) [START_REF] Goloboff | Techniques for analyzing large data sets[END_REF]: given many trees of same score, this method uses sub-groups with an identical taxa composition but different topologies and switch them, the resulting trees are then evaluated; it can be considered as a kind of crossover operator, -Sectorial search (SS) [START_REF] Goloboff | Techniques for analyzing large data sets[END_REF]: consists in optimizing a subtree t of the tree T , in other words to find a topology of the subtree that can give the whole tree a lower score.

Software

Table 2.2 presents a list of software related to the resolution of Maximum Parsimony. We did not put in this list old software like NONA [START_REF] Goloboff | [END_REF] or platforms like MEGA. A complete list can be found on the web page of Joe Felsenstein 6 . POY (Wheeler, Lucaroni, Hong, Crowley, and Varón, 2015) for example is close to TNT.

Conclusion

Maximum Parsimony is a challenging problem and much work has been devoted to its resolution that range from a simple descent to memetic algorithms. However, it is sometimes difficult to compare the different implementations as there is no standard benchmark like BaliBase (Nuin, Wang, and Tillier, 2006) for multiple sequence alignments methods. In Chapter 5 we describe a database of more than 12,000 problems that we have gathered in order to predict the MP score of a problem and that could be used to assess the performance of the search and the quality of the solutions obtained by parsimony software.

Simulated annealing for the Maximum Parsimony problem 3.1 Introduction

Simulated Annealing (SA) is a general-purpose stochastic optimization technique that has proved to be an effective tool for the approximation of global optimal solutions to many NP-hard optimization problems. However, it is well known that the design of an effective SA algorithm requires a careful implementation of some essential components and an appropriate tuning of the parameters used (Johnson, Aragon, McGeoch, andSchevon, 1989, 1991).

In this chapter we present an improved implementation of a SA algorithm (see Algorithm 4), that we called SAMPARS to find tree topologies (phylogenies) with near-optimal parsimony costs under Fitch optimality criterion. The algorithm is a simple cooling schedule with a Markov Chain process.

Simulated annealing

The simulated annealing algorithm is based on the analogy between the simulation of the annealing of solids and the problem of solving large combinatorial optimization problems [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF][START_REF] Cerny | A thermodynamic approach to the traveling salesman problem: An efficient simulation algorithm[END_REF].

In the field of statistical thermodynamic, annealing denotes a physical process in which the temperature of a solid is increasing to a maximum value at which all particles of the solid randomly arrange themselves in the liquid phase, followed by a cooling process where the temperature is decreasing slowly. During the process from hot to cold the particles are modeled to reach an optimal form [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF]. [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF] and [START_REF] Cerny | A thermodynamic approach to the traveling salesman problem: An efficient simulation algorithm[END_REF] worked independently and they showed how this process could be applied to optimization problems, associating key concepts of the simulation process with elements of combinatorial optimization.

In order to implement an optimized version of simulated annealing algorithm to solve the phylogenetic reconstruction (see Algorithm 9), we have identified the principal components of this algorithm and we have selected different values to find the combinations that could offer the best quality solutions with reasonable computational effort. We call this algorithm and the resulting software SAMPARS for Simulated Algorithm for the Maximum PARSimony problem.

In the next section we will present all the implementation details of the SAMPARS algorithm, each of its components and the possible combinations of parameters to obtain the best quality solution.

Algorithm 9: SAMPARS algorithm 1 SimulatedAnnealing(s 0 , N , T 0 , T f , α) input: N : neighborhood, f : fitness function, M CL : Markov Chain length, α : cooling scheme, T i : initial temperature, T f : final temperature output: best solution found s ⋆ 2 s 0 ← GenerateInitialSolution()

3 s ← s 0 4 s * ← s 0 5 t ← T i 6 while t > T f do 7 i ← 0 8 while i < M CL do 9 s ′ ← GenerateN eighbor(s, i, N ) 10 ∆f ← f (s ′ ) -f (s) 11 generate a random u ∈ [0, 1] 12 if (∆f < 0) or (e -∆f /t > u) then 13 s ← s ′ 14 if f (s ′ ) < f (s * ) then 15 s ⋆ ← s ′ 16 i ← i + 1 17 t ← αt 18 return s ⋆

Initial solution

The initial solution s 0 is the starting phylogenetic tree used for the algorithm to begin the search for better configurations in the search space T . SAMPARS can create the starting solution using a random or greedy procedure. A greedy initial solution is supposed to guarantee a better quality for the final solution. Other initialization methods could be used:

upgma, nj: distance methods like UPGMA [START_REF] Sneath | Numerical Taxonomy. The Principles and Practice of Numerical Classification[END_REF] or Neighbor-Joining (NJ) [START_REF] Saitou | The neighbor-joining method: a new method for reconstructing phylogenetic trees[END_REF][START_REF] Gascuel | On the optimization principle in phylogenetic analysis and the minimum-evolution criterion[END_REF] as they have a low complexity nj + fitch: use of the NJ distance method where the distance is Fitch evaluation φ; at each step of the NJ algorithm we recompute the distance matrix using Fitch scoring function nj + greedy: we have also designed a NJ and greedy procedure which uses a random selection of half of the sequences to generate a first NJ tree, and adds the second half to this tree in a greedy manner

Influence of initial solution

We have tried to determine which initialization procedure was the most useful. The initial solution can sometimes influence the search. For example, with a descent algorithm (local search) that employs a NNI (Nearest Neighbor Interchange) neighborhood, if one starts from a solution with a high parsimony score then the final solution will generally be of poor quality as the algorithm will converge prematurely to a local optimum [START_REF] Goëffon | Progressive tree neighborhood applied to the maximum parsimony problem[END_REF].

Tests performed on set1 (see Table 5.3 in chapter 5 and Table 3.1) show that on average random > upgma > nj > nj + greedy > nj + f itch > greedy. In other words, the random initialization procedure produces initial solutions of worst quality while a greedy algorithm gives the best configurations. For set2 (results not reported here), we get slightly different results: random > upgma > nj > greedy > nj > greedy > nj + f itch.

Neighborhood functions

The most common practice in the reported metaheuristics for the MP problem (Andreatta and Ribeiro, 2002; Ribeiro andVianna, 2005, 2009) is to employ one of the three classical neighborhood functions: NNI, SPR, TBR (see Figure 3.1). [START_REF] Barker | LVB: Parsimony and simulated annealing in the search for phylogenetic trees[END_REF][START_REF] Barker | LVB homepage[END_REF] alternates the use of the NNI and SPR neighborhood functions at each iteration of the search process. In SAMPARS both the SPR and the TBR neighborhood relations are implemented but can only be used one at a time. However, from our preliminary experiments it has been observed that the separate use of these neighborhood functions is not sufficient to reach the best-known solutions, because both of them are highly disruptive. In order to achieve a better performance for SAMPARS, we have decided to use a third complementary neighborhood structure. It is based on a stochastic descent algorithm with a best-improvement scheme (see Algorithm 10) which is occasionally applied to the neighboring solution s ′ prior to returning it. Our neighborhood function is inspired by the ideas reported in [START_REF] Lü | Neighborhood analysis: A case study on curriculum-based course timetabling[END_REF], where the advantage of using this approach is well documented. The BestImproveDescent consists in a descent for which all neighbors of the current solution are processed and only the one that gives the best result is kept.

In Algorithm 10, every 25 iterations of the Markov Chain we improve the new configuration with a descent algorithm. This mechanism seems to greatly improve the parsimony scores. However, the frequency used (here 25) ideally should not be a constant but rather a proportion of the iterations per Markov chain (e.g. 1 out of 4 iterations) since the parameterization of SA is usually dependent on each problem instance it seems more sensible to use a ratio for frequency rather than an absolute value. But in order to to study the influence of the frequency of improvements we have decided to keep it a fixed value that we make

Experiments and results obtained for the neighborhood function

Results obtained for different instances (see Tables 3.2 and 3.3 for instance tst12 in particular) show that:

-the computation time increases with a small frequency (f req) as there are more descents to perform: for example for instance tst12, with a frequency of 5, the computation time is around 72 minutes, while with a frequency of 10 it is only 31 minutes, -if the frequency is too low (< 20) then it is more difficult to reach the best known optimal solution as we get stuck in local optima. For f req = 10, the average score of the solutions is 1226.50 and the minimum score found is 1222 while for f req = 20 the average score is 1219 and the minimum score found reaches 1213, -if the frequency is too high (> 200) then it is more difficult to reach a good solution, -depending on the problem, the interval for which we can obtain good results is 20 ≤ f req ≤ 100.

Cooling schedule

A cooling schedule is defined by the following parameters: an initial temperature T i , a final temperature T f or a stopping criterion, the number of solutions generated at each temperature (Markov Chain Length = M CL), and a rule for decrementing the temperature. The cooling schedule governs the convergence of the SA algorithm. At the beginning of the search, when the temperature has a high value, the probability of accepting solutions of worse quality than the current solution (uphill moves) is high. It allows the algorithm to escape from local minimum. The probability to accept such moves is gradually decreased as the temperature decreases to zero.

The literature offers a number of different cooling schedules, see for instance ( dynamic. In a static cooling schedule, the parameters are fixed and cannot be changed during the execution of the algorithm. With a dynamic cooling schedule the parameters are adaptively changed during the execution.

In the SAMPARS implementation we preferred a geometrical cooling scheme mainly for its simplicity, with a reheating system that makes it dynamic. It starts at an initial temperature T i that can either be defined by the user or automatically computed using the following formula: 3 (k + n) which generates values under 6.0 for most of the tested benchmark instances, which was found to give good results. Then, this temperature is decremented at each round by a factor α = 0.99 using the relation t = αt. A reheat mechanism has also been implemented. If the best-so-far solution is not improved during 40 consecutive temperature decrements, the current temperature t is increased by a factor β = 1.4 using the function t = βt. In our implementation this reheat mechanism can be applied at most max_reheat = 4 times, since it represents a good trade-off between efficiency and quality of solution found.

For each temperature t, the maximum number of visited neighboring solutions is M CL. It depends directly on the parameters n and k of the studied instance, since we have observed that more moves are required for bigger trees (Vázquez-Ortiz and Rodríguez-Tello, 2011). The three different values that M CL can take and that we call effort were empirically decided, are:

-small:

M CL = 15 × (n + k) -medium: M CL = 23 × (n + k) -large: M CL = 40 × (n + k)
The parameter values presented in this section were chosen based on the parsimony score in a preliminary experimentation (Vázquez-Ortiz and Rodríguez-Tello, 2011).

Influence of the reheat mechanism and effort level

In order to determine the influence of the reheat mechanism we ran a total of 26,000 tests for medium and large effort levels, different initialization methods (random, nj, greedy, nj + greedy, nj + f itch), with different improvement frequencies (20, 30, 60, 90) and reheat values that range from 0 (no reheat) to 4 (reheat is performed at most 4 times) and a SPR neighborhood.

Results are shown in Table 3.4. For each problem we provide the average parsimony score obtained, the minimum score that could be reached and the number of times it was reached out of 20 executions.

We can see that problem tst10 does not seem to be too complex because we could find the minimum score of 720 for every configuration. However, we can see that permitting up to 4 reheats enables the minimum score to be found 38 times over all tests performed, i.e. a total of 400 (5 methods, 4 frequencies and 20 executions). Problem tst20 is a more complex problem as the minimum score of 659 could only be found with a large effort level with a reheat of 1 or 4.

The results confirm our former observations: a large effort level will lead to results of better quality compared to a medium effort level and a reheat of 1 or 4 can provide better results.

In a second experimentation (Table 3.5) we have tested some methods to get the initial solution combined with effort level and a given number of reheats in order to assess if the initial solution has some influence on the final result. The answer, obtained from the results, is that the nj and nj + greedy methods will on average provide better results and help to reach a maximum number of good solutions (see results in boldface).

Computational experiments

This section reports on three main experiments were conducted to evaluate the performance of the SAMPARS algorithm. The objective of the first experiment is to determine both a component combination, and a set of parameter values which enables SAMPARS to attain the best trade-off between solution quality and computational effort. The purpose of the second experiment is to carry out a performance comparison of SAMPARS with respect to an existing SA algorithm called LVB [START_REF] Barker | LVB: Parsimony and simulated annealing in the search for phylogenetic trees[END_REF][START_REF] Barker | LVB homepage[END_REF]. The third experiment is devoted to assess the performance of SAMPARS with respect to three representative state-of-the-art procedures: GA+PR+LS (Ribeiro and Vianna, 2009), TNT (Goloboff et al., 2008a) and Hydra [START_REF] Goëffon | Nouvelles Heuristiques de Voisinage et Mémétiques pour le Probléme Maximum de Parcimonie[END_REF].

For these experiments SAMPARS was coded in C++ and compiled with g++ using the optimization flag -O3. It was run on a cluster composed of Xeon X5650 CPUs under Linux 64 bits operating system. Due to the non-deterministic nature of the studied algorithms, 30 independent runs were executed for each of the selected benchmark instances in each experiment presented in this section.

Benchmark instances and performance assessment

The test-suites that we have used in our experiments are the same proposed by Ribeiro and Vianna (Ribeiro andVianna, 2003, 2005) and later employed in other works [START_REF] Goëffon | Nouvelles Heuristiques de Voisinage et Mémétiques pour le Probléme Maximum de Parcimonie[END_REF][START_REF] Ribeiro | A hybrid genetic algorithm for the phylogeny problem using path-relinking as a progressive crossover strategy[END_REF]. It consists of two sets (see section 5.2). A first set (set1 ) of 8 benchmarks obtained from real data. A second set (set2 ) of 20 randomly generated instances (tst01 to tst20) with a number of sequences (n) ranging from 45 to 75 whose length (k) varies from 61 to 159. The set2 is composed of some problems for which it is hard to reach the global optimum while for set1 the instances are quite easy to solve. The criteria used for evaluating the performance of the algorithms are the same as those used in the literature: the best parsimony cost found for each instance (smaller values are better) and the CPU time in seconds.

Components and parameters tunning

Optimizing parameter settings is an important task in the context of algorithm design. Different procedures have been proposed in the literature to find the most suitable combination of parameter values (de Landgraaf et al., 2007; [START_REF] Gunawan | Fine-tuning algorithm parameters using the design of experiments[END_REF]. In this work we employ a tuning methodology based on Combinatorial Interaction Testing (CIT) [START_REF] Cohen | The combinatorial design approach to automatic test generation[END_REF]. We have decided to use CIT, because it allows to significantly reduce the number of tests (experiments) needed to determine the best parameter settings of an algorithm. Instead of exhaustively testing all the parameter value combinations of the algorithm, it only analyzes the interactions of t (or fewer) input parameters by creating interaction test-suites that include at least once all the t-way combinations between these parameters and their values.

Covering arrays (CAs) are combinatorial designs which are extensively used to represent those interaction test-suites. A covering array CA(N ; t, k, v), of size N , strength t, degree k, and order v, is an N × k array on v symbols, such that, every N × t sub-array includes at least once all the ordered subsets from v symbols of size t (t-tuples) [START_REF] Colbourn | Combinatorial aspects of covering arrays[END_REF]). The minimum N for which a CA(N ; t, k, v) exists is the covering array number and it is defined according to the following expression:

CAN(t, k, v) = min{N : ∃ CA(N ; t, k, v)}.
CAs are used to represent an interaction test-suite as follows. In an algorithm we have k input parameters. Each of these has v values or levels. An interaction test-suite is an N × k array where each row is a test case, and each column represents an input parameter. The value at row i, column j is the particular setting of the j t h parameter in the i t h test case. This test-suite covers all the t-way combinations of input parameter values at least once. Thus, the costs of tuning the algorithm can be substantially reduced by minimizing the number of test cases N in the covering array.

In practice, algorithms' input parameters do not have exactly the same number of values (levels). To overcome this limitation of CAs, mixed level covering arrays (MCAs) are used.

A MCA(N ; t, k, (v 1 , v 2 , • • • , v k )
) is an N × k array on v symbols, where each column i (1 ≤ i ≤ k) of this array contains only elements from a set S i , with |S i | = v i . This array has the property that the rows of each N × t sub-array cover all t-tuples of values from the t columns at least once. Next, we present the details of the tuning process, based on CIT, for the particular case of our SAMPARS algorithm. First, we have identified k = 8 input parameters used for SAMPARS: neighborhood function N , maximum number of visited neighboring solutions M CL, stop condition SC, reheating factor β, cooling factor α, initial temperature T i , initial solution procedure IS and number of reheats N R. Based on some preliminary experiments, certain reasonable values were selected for each one of those input parameters (shown in Table 3.6).

√ n 0.95 - - - - 3 NJ + Fitch 3 √ n - - - - - 4 - - - - - - -
The smallest possible mixed level covering array MCA(240; 4, 8, (2, 2, 2, 2, 3, 4, 4, 5)), shown (transposed) in Table 3.7, was obtained by using the Memetic Algorithm reported in (Rodriguez-Tello and Torres-Jimenez, 2010). As reference, we can map 0 in the last column (the first line in Table 3.7) to SPR and 1 to TBR. The resulting interaction test-suite contains, thus, 240 test cases (parameter settings) which include at least once all the 8-way combinations between SAMPARS's input parameters and their values 1 Each one of those 240 test cases was used to execute the SAMPARS algorithm 30 times over the 20 instances of the test-suite described in Section 3.3.1.

From the results obtained we have selected the 5 test cases which yield the best results. Their average parsimony cost and the average CPU time in seconds are presented in Table 3.8. This table allowed us to observe that the parameter setting giving the best trade-off between solution quality and computational effort corresponds to the test case number 59 (shown in bold). The best average parsimony cost with an acceptable speed is thus reached with the following input parameter values: number of reheat N R = 4, initial solution procedure IS = N J, initial temperature T i = 6.0, cooling factor α = 0.99, cooling factor β = 1.6, stop condition SC = T f , maximum number of visited neighboring solutions M CL = Big, neighborhood function N = SP R. These values are thus used in the experimentation reported next.

SAMPARS compared to an existing SA implementation

For this experiment a subset of six representative benchmark instances, taken from the test-suite described Section 3.3.1, was selected (comparable results were obtained with all the other tested instances). Then, the latest version of LVB was obtained, compiled and executed on our computational platform using the input parameters suggested by the author [START_REF] Barker | LVB homepage[END_REF]. Table 3.9 displays the detailed computational results produced by this experiment. The first three columns in the table indicate the name of the instance as well as its number of taxa (n) and length (k). For each compared algorithm the best (B), average (Avg.), and standard deviation (Dev.) of the parsimony cost attained in 30 independent executions and its average CPU time in seconds are listed in columns 4 to 11. A statistical significance analysis was performed for this experiment. First, D'Agostino-Pearson's omnibus K 2 test was used to evaluate the normality of data distributions. For normally distributed data, either ANOVA or the Welch's t parametric tests were used depending on whether the variances across the samples were homogeneous (homoskedasticity) or not. This was investigated using the Bartlett's test. For non-normal data, the nonparametric Kruskal-Wallis test was adopted. A significance level of 0.05 has been considered. The resulting P -value is presented in Column 12.

From Table 3.9 we can observe that SAMPARS is the most time-consuming algorithm, since it uses an average of 539.57 seconds for solving these six instances. On the contrary, LVB employs only 288.16 seconds. However, we can also remark that SAMPARS can take advantage of its longer executions. Indeed it is able to consistently improve the best results found by LVB, 1. In contrast, with an exhaustive testing which contains 5( 42 ) × 3(2 4 ) = 3840 test cases. 189 190 191 192 193 194 195 196 197 198 199 200

3 0 0 4 1 0 3 4 1 0 4 0 0 2 3 2 3 2 3 2 1 2 2 2 0 3 3 1 1 0 0 0 2 0 3 0 0 3 2 3 0 0 2 2 2 1 1 0 3 2 3 3 0 1 1 0 0 1 2 0 1 2 0 1 0 2 2 0 1 0 1 0 2 0 0 1 1 1 0 0 1 0 1 0 1 0 1 1 0 1 1 0 1 1 1 0 0 0 0 0 0 1 0 1 1 0 0 0 0 1 1 1 0 0 1 0 0 0 1 0 0 0 0 0 1 1 1 1 0 1 0 1 1 0 1 0 0 0 0 1 1 1 0 0 1 1 0 1 1 1 0 0 1 0 1 0 1 1 1 1 1 0 0 0 1 1 0 0 0 1 1 1 0 1 1 1 1 1 0 1
3 3 0 0 3 2 1 3 1 0 4 3 2 4 3 2 1 0 3 0 1 1 0 1 2 1 0 1 3 3 3 0 1 2 3 3 0 2 0 0 0 0 1 2 1 0 1 3 2 1 3 2 2 3 3 0 3 1 1 0 0 0 1 3 1 2 1 1 2 2 0 1 1 2 1 2 2 1 2 0 1 0 1 0 0 0 1 2 1 0 1 1 0 0 1 0 1 1 0 1 0 1 1 1 1 0 0 1 1 1 0 1 1 1 0 0 1 0 0 0 1 1 1 0 1 1 0 0 1 1 0 1 0 0 0 1 1 0 0 0 1 1 0 1 1 0 0 1 0 0 1 0 0 1 1 1 1 1 0 0 0 1 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 1 1 0 0 1 0
3 0 0 3 1 2 0 0 2 0 0 2 1 1 1 1 3 1 0 3 3 1 2 3 2 1 2 2 1 3 1 1 3 2 3 0 3 3 1 3 3 2 0 0 0 0 1 2 0 0 3 2 1 1 1 0 3 3 1 1 1 1 2 0 2 1 2 2 0 1 1 1 1 2 1 2 0 2 2 0 2 0 2 1 1 0 0 1 0 0 1 1 1 1 0 0 0 0 1 1 0 1 0 0 0 0 1 0 0 1 1 0 0 1 0 0 0 0 0 1 0 1 0 0 1 1 1 1 1 0 0 0 1 1 0 0 1 0 0 0 1 1 0 1 1 0 1 1 1 1 1 0 1 0 0 0 1 0 1 1 0 1 0 0 0 0 0 1 0 0 1 0 1 1 0 1 1 1 1 0 1
3 3 0 3 0 1 3 0 0 2 1 3 0 3 2 2 0 2 0 3 3 0 2 0 0 2 2 0 3 2 0 1 1 1 1 2 2 2 2 2 1 1 2 0 2 3 2 0 2 0 2 0 0 2 1 0 2 1 0 0 0 2 2 2 0 0 1 1 1 1 0 2 0 1 1 0 1 0 0 0 1 1 0 0 0 0 1 1 1 1 0 1 1 1 0 1 0 0 1 1 1 1 0 1 0 0 0 0 0 1 1 0 0 1 1 0 1 0 0 0 0 0 0 1 1 0 1 0 1 1 1 0 1 1 1 0 1 0 1 0 0 1 0 1 0 0 1 1 1 0 0 0 1 1 1 1 1 0 1 0 0 1 0 0 0
2 3 0 1 0 3 1 4 3 4 2 3 4 1 3 2 1 2 2 3 1 1 2 0 1 1 3 0 2 3 0 0 1 0 0 2 0 1 2 1 1 3 3 2 3 3 3 1 1 3 0 3 1 3 2 1 1 0 2 2 1 2 2 1 2 1 0 1 2 1 2 0 1 1 1 0 1 0 0 0 1 0 1 0 1 0 1 1 0 1 1 0 1 1 1 1 1 0 1 0 1 1 1 0 0 1 0 0 1 1 0 1 1 0 0 1 0 1 1 0 1 1 1 0 0 0 0 0 1 0 1 0 1 0 0 0 0 0 1 1 1 0 0 0 1 0 0 1 1 0 0 1
0 4 3 1 3 0 4 4 0 1 2 0 3 1 2 1 4 2 3 3 0 2 0 0 2 2 0 1 1 0 3 2 1 1 3 1 3 1 3 3 2 1 1 2 3 1 3 3 2 0 1 0 2 1 3 1 2 1 2 2 2 1 0 2 0 2 0 2 0 0 0 1 1 2 0 1 2 2 0 0 0 0 1 1 0 1 1 0 1 1 0 0 0 1 0 0 0 0 1 1 1 0 0 1 0 0 1 0 1 1 1 1 1 0 0 1 0 0 0 1 0 0 1 1 1 0 0 0 0 0 0 0 0 0 1 0 0 0 1 1 0 0 0 1 0 1 1 0 0 1 0 1 0 1 0 0 1
3 3 3 1 1 2 2 3 0 2 0 1 3 0 2 3 3 1 1 2 0 2 3 0 3 0 1 1 0 3 1 0 0 3 0 0 3 2 0 0 2 2 0 0 0 1 2 0 1 2 0 2 1 1 0 1 1 2 1 2 0 1 0 1 0 1 0 1 0 0 1 1 1 1 0 0 1 1 1 0 0 1 1 0 1 1 0 0 1 0 1 0 1 1 1 0 0 0 0 1 0 1 1 0 1 0 1 0 0 0 0 0 1 0 0 0 0 1 1 0 1 1 1 0 1 1 0 0 1 0 1 0 0 1 1 1 0
2 2 0 1 2 1 0 1 0 1 2 4 3 0 1 2 3 1 1 2 2 2 1 3 1 0 3 2 1 0 0 3 1 3 3 2 3 3 0 0 1 0 2 2 0 1 1 2 3 3 2 0 2 0 3 2 1 1 3 2 2 2 1 1 0 2 1 0 0 0 0 1 0 0 1 1 1 2 2 0 1 0 0 1 0 0 0 0 1 0 1 1 1 1 0 0 1 0 1 0 1 0 1 1 1 1 0 0 1 1 1 1 1 1 1 0 0 1 1 1 0 1 1 1 0 0 1 0 1 1 0 1 0 1 1 0 1 0 0 1 0 0 1 0 0 0 1 1 0 0 0 1 0 1 1 1 1 0 0 185 186 187
0 0 2 4 1 2 3 0 4 2 0 1 2 2 0 2 0 3 0 2 0 3 1 2 1 1 1 2 1 0 1 1 0 3 3 3 2 1 0 0 2 0 3 3 0 3 1 2 2 1 2 0 0 0 2 1 2 0 1 0 2 2 2 1 2 0 0 1 0 1 0 1 1 1 1 1 1 1 0 0 1 1 1 0 0 1 1 0 0 0 1 0 1 0 0 1 1 1 0 1 0 0 0 1 0 1 1 0 0 0 1 0 1 1 1 1 0 0 1 0 1 0 0 0 1 1 1 1 1 0 1 0 1 0 1 1 0 0 1 0 1 1 1 1 1 1 1 0 0
2 3 3 1 1 2 0 1 2 2 1 3 2 2 2 3 0 3 2 1 1 0 1 3 2 0 2 0 2 3 2 3 1 2 3 0 3 1 3 3 2 0 0 0 1 2 2 0 0 2 0 2 2 2 0 0 1 1 1 0 1 1 0 1 0 1 0 0 1 1 0 1 1 0 1 0 0 0 1 1 1 0 0 0 1 0 0 0 0 1 0 1 1 1 1 1 1 0 0 1 0 1 1 1 0 1 0 0 0 1 1 0 1 1 0 0 1 0 1 0 1 1 1 1 0 0 1 1 1 0 1 1 0 1 0 0 0 1 0 1 0 0 1
0 2 3 0 3 3 1 1 3 1 0 2 2 0 3 3 0 3 1 0 0 0 2 0 1 2 2 2 1 1 0 2 2 1 2 2 0 0 1 0 1 1 0 0 1 1 1 1 1 1 1 1 0 0 0 0 0 1 1 0 1 0 1 1 0 0 1 0 1 1 0 1 0 1 0 1 0 1 0 0 0 0 0 1 1 0 1 1 0
obtaining in certain instances, like tst08, an important decrease in the parsimony cost (up to -15 = 852 -867). Furthermore, the solutions found by SAMPARS present a relatively small standard deviation (see column Dev.). It is an indicator of the algorithm's precision and robustness since it shows that in average the performance of SAMPARS does not present important fluctuations. We have noticed that LVB has a erratic behavior (see Table 3.10). For example, the resolution of problem tst04 can take from 84 seconds to 5 hours and the time spent for the resolution is not correlated to the finding of a better solution. The statistical analysis presented in the last two columns of Table 3.9 confirms that there exist a statistically significant increase in performance achieved by SAMPARS with regard to LVB on the six studied instances. Thus, we can conclude that SAMPARS is more effective than the existing SA algorithm reported in [START_REF] Barker | LVB: Parsimony and simulated annealing in the search for phylogenetic trees[END_REF][START_REF] Barker | LVB homepage[END_REF].

Another implementation called Tree-Drifting is described in [START_REF] Goloboff | Techniques for analyzing large data sets[END_REF] as a TBR descent able to accept suboptimal trees. The key component of the method is the function designed to determine the probability of acceptance, which is based on both the absolute step difference and a measure of character conflict: the relative fit difference, which is the ratio of steps gained and saved in all characters, between the two trees. A pseudocode of the algorithm is given in (De [START_REF] Laet | Pseudocode for some tree search algorithms in phylogenetics[END_REF]. Tests performed using only tree-drifting in TNT (Goloboff, Farris, and Nixon, 2008a) gave results far from the best known values (results not presented here).

Comparison of SAMPARS with the State-of-the-art procedures

In this experiment a performance comparison of the best solutions achieved by SAMPARS with respect to those produced by GA+PR+LS [START_REF] Ribeiro | A hybrid genetic algorithm for the phylogeny problem using path-relinking as a progressive crossover strategy[END_REF], TNT (Goloboff, Farris, and Nixon, 2008a) and Hydra [START_REF] Goëffon | Nouvelles Heuristiques de Voisinage et Mémétiques pour le Probléme Maximum de Parcimonie[END_REF] was carried out over the test-suite described in section 3.3.1.

TNT (Tree analysis using New Technology) is probably the fastest and one of the most effective parsimony analysis program for the MP problem. TNT is known for finding better trees several thousands times faster than other software. TNT uses many search strategies [START_REF] Goloboff | Techniques for analyzing large data sets[END_REF] coming from genetic algorithms, local search and supertrees. During the local search phases through the use of a tuning methodology based on Combinatorial Interaction Testing [START_REF] Cohen | The combinatorial design approach to automatic test generation[END_REF], to yield the best solution quality in a reasonable computational time.

An extensive experimentation was conducted to investigate the performance of SAMPARS over a set of 20 well-known benchmark instances. In these experiments our algorithm was carefully compared with an existing Simulated Annealing implementation (LVB) [START_REF] Barker | LVB: Parsimony and simulated annealing in the search for phylogenetic trees[END_REF][START_REF] Barker | LVB homepage[END_REF], and other three state-of-the-art algorithms GA+PR+LS, TNT and Hydra. The results show that there exists a statistically significant increase in performance achieved by SAMPARS with respect to LVB. SAMPARS is in fact able to consistently improve the best results produced by LVB, obtaining in certain instances important reductions in the parsimony cost. Regarding GA+PR+LS [START_REF] Ribeiro | A hybrid genetic algorithm for the phylogeny problem using path-relinking as a progressive crossover strategy[END_REF], we have observed that on average this algorithm returns worse quality solutions than SAMPARS. Compared with the state-of-the-art algorithm called Hydra [START_REF] Goëffon | Nouvelles Heuristiques de Voisinage et Mémétiques pour le Probléme Maximum de Parcimonie[END_REF] our SAM-PARS algorithm was able to improve 9 previous best-known solutions and to equal the results on the other 11 selected benchmark instances. Furthermore, it was observed that the solution cost found by SAMPARS presents a relatively small standard deviation, which indicates the precision and robustness of the proposed approach. These experimental results confirm the practical advantages of using our algorithm for solving the MP problem.

A Bottom-up implementation of Path-relinking

Introduction

Path Relinking (PR), has proved unusually effective for solving a wide variety of optimization problems from both classical and real world settings. PR operates with a population of solutions, rather than with a single solution at a time, and employs procedures for combining these solutions to create new ones [START_REF] Glover | Fundamentals of scatter search and path relinking[END_REF].

In this chapter we describe a bottom-up implementation of Path-Relinking for phylogenetic trees. This bottom-up implementation is compared to two versions of an existing top-down implementation. We show that our implementation is more efficient, more interesting to compare trees and to give an estimation of the distance between two trees in terms of the number of transformations.

Path-relinking

Path-relinking (PR) is a metaheuristic with an intensification strategy used to explore elite solutions. It was defined by [START_REF] Glover | Fundamentals of scatter search and path relinking[END_REF] and it is closely related to Tabu Search. The main objective of PR is to generate paths between and beyond a group of solutions by combining them, where solutions on such paths also serve as sources for generating additional paths. In order to get good paths, the solutions to combine must be selected in a neighborhood space, rather than in Euclidean space [START_REF] Glover | A template for scatter search and path relinking[END_REF].

PR has a group of guiding solutions with weighted attributes to determine which moves are given higher priority. The generation of such paths in neighborhood space characteristically relinks previous points in ways not achieved in the previous search history, hence giving the approach its name.

Given two solutions called source and guiding, PR consists in transforming the source solution into the guiding solution by applying a series of modifications. After each modification an exploration phase is performed on a copy of the current solution. The aim of PR is to generate a path from the source to the guiding solution in order to possibly find a better solution.

PR has been used in different domains for the resolution of combinatorial optimization problems (Wang et Ribeiro and Vianna (2009) describe an implementation of PR tailored to trees for the resolution of the Maximum Parsimony problem (see algorithm 11). The symmetric difference ∆(x, y) computes the set of moves to go from solution x to solution y. At each step the algorithm selects m ⋆ the move that minimizes the fitness function from the current solution s i to t. If the new solution s i+1 = s i ⊕ m ⋆ is better than s ⋆ then s ⋆ is updated. On randomly generated instances they observed that using GRASP + PRm they could always find better solutions than a simple GRASP algorithm.

Complexity of PR-Bottom-up

The complexity of the algorithm can be computed as follows: given n, the number of taxa of the problem, the reordering of the guiding and source trees needs n -1 comparisons and the computation of Ω g can be done in 2n -1 operations. The main loop will be executed a certain number of times, let's say p times, and we will need to compute Ω s , find a missing subtree (the maximum will be n comparisons) and perform a degraph and regraph (1 transformation). This adds up to:

2 × (n -1) + 2n -1 + p × (2n -1 + n + 1) = 4n -3 + p × (3n) ≃ 3n × p (4.1)
In the worst case p = n, so the worst case complexity is O(n 2 ) for the bottom-up implementation.

For the top-down implementation, the computation of the complexity is more difficult to establish as the process is recursive and the size of the subtrees changes. At each step we must compute the sets of leaves on the left and right subtrees of the source and guiding solutions. Then move the leaves in the source tree from left to right or from right to left. The number of transformations (see the tables in the results section) can give us some insight into the complexity of this implementation.

Experimentation and results

Benchmark

We used the instance called zilla (500 sequences of 759 DNA residues) that was originally obtained from the chloroplast gene rbcL (Chase, Soltis, Olmstead, Morgan, Les, Mishler, Duvall, Price, Hills, Qiu, Kron, Rettig, Conti, Palmer, Manhart, Sytsma, and al, 1993) and the best parsimony score of 16,218 was first found by TNT. In Table 4.2 we report the percentage of common subtrees that we could find between two trees that are being compared. The number of common subtrees can be thought as a measure of similarity between trees. Those different trees were obtained using SAMPARS.

For example two trees of score 16,218 (namely 16,218 and 16,218 b ) have 74.75% subtrees in common which means that they are topologically close. On the contrary the tree of score 21,727 has 0% subtrees in common with the tree of score 16,218. This means they are very far from each other and have only the leaves in common. bottom-up: iterative implementation explained in this chapter, based on subtree difference, top-down without minimization: recursive implementation without optimization on regraph, i.e. when a leaf is regraphed it is placed at the root of the subtree where it should appear (Ribeiro and Vianna, 2009), top-down with minimization: recursive implementation with optimization on regraph, i.e. when a leaf is regraphed all possible branches are tested and we keep the first one that minimizes the score of the tree (Ribeiro and Vianna, 2009).

Experiments

The experimentations were performed on an Intel Core i5 4570 and the program was coded in Java 1.7, it is part of a software called Arbalet1 . In Table 4.3 we report for each implementation the number of transformations (degraph + regraph), the execution time in seconds, the number of times the source tree had a score inferior or equal to the guiding tree (#Equal) and the number of times the source tree had a score strictly inferior to the guiding tree (#Less) during the generation of the path. Note that the source tree has a higher score than the guiding tree. This is not necessary for the bottom-up method for which we can invert the trees. However this is required by the top-down with minimization implementation. In Table 4.3, with the bottom-up implementation the path between the two trees of score 16,218 (called 16,218 and 16,218 b ) was built with 24 transformations in 0.1 seconds. As mentioned before those trees are topologically very close. During the transformation process the source tree, when modified, had a best score (of 16,218) 13 times among the 24 transformations. With the top-down with minimization algorithm (see results of Table 4.3) the construction of the path of the tree of score 21,727 into the tree of score 16,611 has needed 1528 transformations and took 22.4 seconds. It also lead to the generation of 586 trees of score under 16,611. This means that the top-down with minimization algorithm can sometimes help find a tree of lower score than the guiding tree.

Execution time

Although the bottom-up implementation requires to recompute the set of subtrees of the source tree for each iteration it is the fastest method. The top-down with minimization implementation takes much more time because of the optimization phase on regraph. However for the top-down implementations the number of transformations in terms of nodes to degraph and regraph is very important compared to the bottom-up version. For example for the search of a path from tree 21,727 to 16,218, the number of recursive calls is equal to 499. The number of leaves moved for the first 10 steps of the recursion are: 194, 143, 21, 3, 15, 34, 40, 4, 43, 3, ... for a total of 1881 transformations.

Comparison of trees

In the bottom-up implementation the number of transformations is proportional to the topological modifications of the tree. For example, on Table 4.3 the number of transformations from the tree of score 16,219 to 16,218 is 32, while the number of transformations from the tree of score 21,727 to 16,218 is 455. For the top-down implementation without minimization the transformation from the tree of score 16,401 to 16,218 is nearly equal to the number of transformations of the tree of score 21,727 to 16,218. Based on this results we could say that our implementation is more suitable to compare trees than the top-down of [START_REF] Ribeiro | A hybrid genetic algorithm for the phylogeny problem using path-relinking as a progressive crossover strategy[END_REF].

There exists different metrics and methods to compare trees. One of the most famous is the Robinson-Foulds metric (RF) [START_REF] Robinson | Comparison of phylogenetic trees[END_REF] initially designed for unrooted tree but which has a variant for rooted trees based on clusters (Lin, Rajan, and Moret, 2012). RF is equal to the number of different splits in compared trees. A split A|B of a set L is an unordered pair (ie, A|B = B|A) of its subsets, such that L = A ∪ B and A ∩ B = ∅. Our method can to some extent be brought closer to the Matching Cluster metric for rooted trees (Lin, Rajan, and Moret, 2012) but provides different results.

Ability to find better solutions

When PR does not use the exploration phase it is unable to find a tree with a score lower than the guiding tree (see columns #Less and #Equal), except for the top-down with minimization version but only for trees that are far from the best known solution.

Nevertheless it is possible to add a minimization phase to the bottom-up algorithm: the new subtree t = (X, Y ) obtained from line 11 of Algorithm 12 can be degraphed and regraphed somewhere on the tree in order to minimize the parsimony score. However t must not be regraphed on some previously modified subtree in order to avoid any cycle of transformation that would cause a modification to be undone and that would be performed again at the next iteration. We have tested this method but it does not lead to any improvement and we could not find a score inferior to the score of the guiding tree.

Conclusion

Path-Relinking in its bottom-up implementation represents an interesting tool to compare the topologies of the source and guiding trees. The bottom-up iterative implementation what we have described is faster than the top-down recursive implementation and can serve as a measure of distance between trees and could be applied to any other context when one has to compare trees: for example when dealing with a genetic algorithm, the replacement of some individual in the population by a new individual could only be possible if the new individual is not too close to the individuals that constitute the population.

Prediction of the parsimony score

Introduction

Given a solver it is always difficult to know the combination of parameters that will provide a result close to the optimum solution that we want to find, or in other words, how could we know how far or close is the final solution we obtain from the global optimum ? In this chapter we present a predictor designed to estimate the optimum score of an instance based on a database of solved instances. As we shall see later on, knowing that the estimate of the best score is far from the random or greedy score, we have the information that the problem is probably difficult to solve and that it would be preferable to solve it using powerful techniques. On the contrary if we know that the problem is easy to solve, it is not necessary to use resolution methods and techniques that will cost a lot of computation and consequently a lot of useless time and efforts. This work is still under way.

The main idea

In the first instance, the idea of designing a predictor is not new. It comes from the fact that it seems easier to estimate a solution than to solve the entire problem. However the design of a predictor for MP is a new concept, to which little (or no) attention has yet been paid and came from an idea we had in mind and that was partially covered by a small study we made using the R mathematical tool.

We took a set of twenty problems that we call set1 (see Table 5.1 for the description of the sets of problems) for which we had computed different initialization trees and their scores. We had generated a table with the following numerical data for each problem:

-K the number of sequences 1 -L the length of a sequence or number of residues (all sequences have the same length) -R the average score (over 100 runs) of initial solutions generated randomly -N the score of the tree obtained from the neighbor-joining algorithm (only one tree) -G the average score (over 100 runs) of initial solutions generated with a greedy algorithm -B the best known score, considered as the optimum score obtained by using TNT or SAMPARS We have decided to use a multiple linear regression (MLR) model [START_REF] Xin | Linear Regression Analysis: Theory and Computing[END_REF], mainly for its simplicity, in order to establish a relationship between (K, L, R, N, G) and B. The principle of MLR is to have as input a matrix M and a vector Y of observed data for which we want to determine if M and Y are correlated. We want to find a vector X and constants α and ǫ such that :

α + M.X + ǫ = Y ′ ≃ Y (5.1)
Generally linear regression models are often fitted using the least squares approach where it is necessary to minimize

|Y -Y ′ | 2 = i (y i -y ′ i ) 2
in order to have Y ′ as close as possible of Y . To be able to compute the MLR, the matrix M (u × p) will contain u instances for which we have p properties and u ≥ p.

From the data we had gathered and computed (see Table 5.1) we could obtain a model with a R-square (R 2 ) of 0.9999 and a p-value less than 2.2e-16, which means the data are strongly correlated.

Following is the output of the program R, where V2 is K, ..., V6 is G and V7 is B: Residual standard error: 4.497 on 14 degrees of freedom Multiple R-squared: 0.9999,Adjusted R-squared: 0.9999 F-statistic: 4.205e+04 on 5 and 14 DF, p-value: < 2.2e-16

> data = read.
The first column (Estimate) gives us the coefficients of the formula. The last column (Pr(>|t|)) is the p-value and helps test the null hypothesis, i.e. check whether the corresponding variable set to zero would affect or not the result. A variable that has a low p-value is likely to be a meaningful addition to the model. A large p-value suggests that changes in the variable are not associated with changes in the response.

Although the analysis of the p-values, R 2 and F-Test are very important to find out if the data are correlated, we are only interested by the forecasting, i.e. prediction from the MLR and we will not delve into statistical details.

The average error of the prediction was of 2.78 points with a minimum of 0.14 for instance tst05 and a maximum of 10.13 for instance tst02 i.e. 0,74 % of its score. For the sake of simplicity, in the remaining of the chapter, we will report real numbers with two decimal places.

The formula obtained is Y ′ = E(K, L, R, N, G) = 65.65 -1.13 × K -0.52 × L -0.17 The fact that we could obtain a formula which gives a good precision probably comes from the relative similarity of those problems: they were all generated using the same software, are composed of a small number of sequences and a small number of residues (see K and L columns).

We then tried to add more problems to the initial set of twenty problems and recomputed the MLR, but it seemed obvious that it would not work for two reasons: we needed more parameters and we needed to take into consideration problems that were close to each other.

For example if we take into account the problems of set1 and set2 for a total of 28 problems, then the average error rises to 7.61 with a minimum of 0.64 for tst20 and a maximum of 22.01 for instance tst15, i.e. 2.92 % of its optimum score.

Database of problems

We have created a database of more than 12,000 problems for which we have calculated all properties cited above. The main goal is to supply enough problems close to each other in order to be able to perform the MLR. The determination of the optimum score B was made using TNT 2 (Goloboff, Farris, and Nixon, 2008a) and the results obtained with SAMPARS. We consider that the best score obtained by the software represent the global optimum. This might not be true for all instances but we take it for granted.

The database or knowledge base is composed of the following sets:

-set1: eight real-life instances by [START_REF] Luckow | An empirical comparison of numerical wagner computer programs[END_REF] -set2: twenty randomly generated instances presented in [START_REF] Ribeiro | A GRASP/VND heuristic for the phylogeny problem using a new neighborhood structure[END_REF] 2. the parameters used to solve the instances are rseed 0, bground, mxram 5096, nstates num, nstates nogaps, rep+1, mult:tbr, mult=ho3, mult=ratchet, sectsch:rss -set3: problems from Treebase 3 -set4: problems from [START_REF] Morrison | Increasing the efficiency of searches for the maximum likelihood tree in a phylogenetic analysis of up to 150 nucleotide sequences[END_REF] -set5: problems taken from (White and Holland, 2011) -set6: problems from (Goloboff, Carpenter, Arias, and Esquivel, 2008b) especially zilla -set7: 12126 problems from TreeFam 4 -set8: 2 subproblems of ribosomal RNA sequences 23S 5 from the Empirical Datasets with Reference Topologies 6 -set9: 40 subproblems of 16S 7 ribosomal RNA trees of FastTree 8 Table 5.2 reports the different global information about each set of problems:

-the number of instances (#instances) -the minimum number of taxa (K min ) -the maximum number of taxa (K max ) -the minimum length for a taxon (L min ) -the maximum length for a taxon (L max ) -the minimum similarity (S min ) -the maximum similarity (S max )

Prediction with the best score

In this section we will explain the method that we have designed in order to predict with a good accuracy the best score of a problem but taking into account the instance in the prediction formula. It will prepare a simple foundation on which to extend to a more complex method when the best score of the instance to predict is not known and will serve:

1. to prove if the MLR is a valid model for the prediction 2. to detect which part of the MLR has more influence for the prediction.

3. treebase.org 4. www.treefam.org/ 5. Ribosomal RNAs (rRNA) perform critical functions in the ribosome that allow protein synthesis to occur. The genes that encode rRNAs evolve in a very specific manner that makes them excellent markers to trace evolutionary history and powerful tools to identifying species from sequence data.

6. http://www.cs.utexas.edu/users/phylo/datasets/phylogeny-topology.html 7. 16S rDNA are used in reconstructing phylogenies because of the slow rates of evolution of this region of the gene. 8. http://www.microbesonline.org/fasttree/#16S

Problem properties

Let's consider a problem P for which we have the following information:

-K the number of sequences -L the length of a sequence (all taxa have the same length)

we can compute a first set of properties:

-R the average score (over 100 runs) of initial solutions generated randomly -N the score of the tree obtained from the neighbor-joining algorithm -G the average score (over 100 runs) of initial solutions generated with a greedy algorithm -B the best score found by TNT, SAMPARS or any over parsimony-based software -S the similarity of sequences computed as the number of equal pairs over the total number of pairs of residues and a second set of properties (see Figure 5.1):

- On Table 5.3 we provide some information about set1 and set2 especially for the properties RN , RG. The set1 of problems is quite easy to solve, we can reach a value close to the best score very quickly. On the contrary the problems of set2 are very difficult to solve. We can note the difference by comparing columns RN and RB:

RN
-for set1, RN ranges from 39 to 68 and RB from 43 to 70 -for set2, RN ranges from 10 to 21 and RB from 17 to 28

The column RN B which stands for the formula 100 -(RN/RB × 100) indicates if we are close to the optimum B when we start from N knowing R 9 . When this quantity is small it means that the tree generated by neighbor-joining is close to B and it could be considered as one of the measures of the difficulty of a problem. For set1 it ranges from 3 to 10 %, and for set2 it ranges from 22 to 41 %. Among all problems of the database it concerns 68 instances if we chose RN B > 10. Note that we could also use |RN -RB| = N B, which is the last column of the

The method

The problem to predict will be represented by a tuple P 0 = (K 0 , L 0 , ..., GB 0 ). The algorithm is composed of three steps, namely:

1. problems selection The results obtained prove that the multiple linear regression model is a good candidate for the prediction of the optimum score of the problem but only if we have problems close to the instance to predict and that, in some way, we already know the solution we want to predict when we integrate P 0 into the MLR.

Prediction without the best score

We follow the same scheme as before but in this case we do not know B 0 so we could not use the selection criterion related to RB. The MLR will concern the properties K, L, N , R, G, S, RN , RG.

Problems selection

default selection

We use the same selection criteria as in Section 5.4.2.1 and we select from the database all problems P = (K, L, ..., GB) that are close to P 0 . The results are of poorer quality (see

Modification of problem selection

In order to get a prediction of better quality we decided to modify the selection criteria, from the database we select instances closer to the problem P 0 .

In this case the three following criteria are met:

-RN = RN 0 ± δ -RG = RG 0 ± δ -S = S 0 ± γ
Initially δ = 0.5 and γ = 5. If the selection does not provide 10 problems we increase δ by a factor of 0.5 and γ by a factor of 5. The results show an increase of the prediction (see Table 5 By studying some predictions we could observe that some of the problems selected to compute the MLR model are sometimes not retained in the final list of 10 problems that constitutes the model. For example we could have a selection of 50 instances but we will retain only 10 out of those 50 instances to compute the MLR. In order to have close instances to P 0 we have decided to define a distance from the selected instances to the problem to predict that will be used to order the problems.

If P and P 0 are close in terms of properties then the RN GS distance will be close to 0. This distance is simply the product of the differences of the properties of P and P 0 and was determined empirically. We needed to add some constant so that if two properties are equal then the overall expression would not be 0 even if the other properties differ.

rngs(P, P

0 ) = (|RN -RN 0 | + cst1) × (|RG -RG 0 | + cst1) × (|N G -N G 0 | + cst1) × (|S -S 0 | + cst2) × (|K -K 0 | + cst2) × (|L -L 0 | + cst2)
where cst1 = 0.001 and cst2 = 0.01.

We then sort the problems according to their increasing RN GS distance. The RN GS distance needs probably to be refined. For example the range of the length of the taxa L is important (from 36 to 152160), it is obvious that the influence of |L -L 0 | can sometimes be too important. Maybe it will be more interesting to use a logarithm to decrease the influence of L. We need to study the influence of each parameter in order to determine a formula with a better precision.

Estimation by weighted average (WA)

In order to improve the quality of the prediction we have decided to use a method different from the MLR. As we have the RB, N B and GB values of the selected problems we have decided to use them in order to predict the optimal score of P 0

We simply compute an average of the scores predicted using RB, N B, GB. Given 10 instances close to P 0 we compute:

predicted RB = 1 10 10 i=1 R 0 - (1 -RB i ) 100 predicted N B = 1 10 10 i=1 N 0 - (1 -N B i ) 100 predicted GB = 1 10 10 i=1 G 0 - (1 -GB i ) 100
For example, predicted RB is the average of the predictions using the RB property. The final score predicted is given by the weighted average:

predicted = 2 × predicted GB + 2 × predicted N B + predicted RB 5
Here it is obvious that GB and N B being closer to B will have less inaccuracy than RB. So we increase the weight of GB and N B.

Results

Table 5.8 shows the predictions obtained by this new method using:

1. the enhanced selection using the RN GS distance to order the selected problems 2. the weighted average to compute the prediction of the score With this new method we get results of very good quality with 96.23% of the instances predicted at ± 1% from the best score.

However, with an accuracy of 5% the problems not predicted are the following: CARP GOLO GRIS tst01 tst02 tst04 tst05 tst06 tst07 tst08 tst10 tst11 tst12 tst14 tst15 tst16 tst19 tst20 m0972 astr bracon dinos kearney norell Those problems are predicted with a precision between 5% and 9%. This is quite strange and needs further investigation because the problems tst01 to tst20 problems formed the basis of the predictor at first. We could observe that there are many TF problems (from set7) which appear as selected and retained to establish the prediction. ------------------------------------------------------------------------------------------- We have included both methods:

% of

-the multiple linear regression with a R 2 = 1.000 (see variable score_mlr) -the weighted average (see variable score_wa) It appears that in this case both methods have a high accuracy but the MLR gives the best prediction. However we use the method based on weighted average to provide the result. The WA works fine because for HIV, we have:

-RB 0 = 42.50 and the average of RB for the selected instances is 42.51 -N B 0 = 1.14 and the average of N B for the selected instances is 1.26 -GB 0 = 1.14 and the average of GB for the selected instances is 1.12

Example that does not work: tst07

In this case both prediction methods (MLR and WA) provide results far from the optimum of 1269. Here we can observe that the values of RB, N B and GB differ for tst07 and the selected instances:

-RB 0 = 20.19 and the average of RB for the selected instances is 14.73 -N B 0 = 9.42 and the average of N B for the selected instances is 2.58 -GB 0 = 6.96 and the average of GB for the selected instances is 1.41

Conclusion

In this chapter we have tried to design the key components of a predictor for the Maximum Parsimony problem based on a knowledge base of solved instances and on a multiple linear regression model. We observed that the MLR was interesting only in the case when we knew what we had to predict. If we removed the RB property from the MLR the percentage of the prediction falls below 86% with the modified problem selection.

In order to get an improved accuracy we had to substitute the MLR by a weighted average which proved to give very interesting results with a percentage of prediction close to 96%. Given that we have all data of the instances, we can deduce that the predictor does not work if the average of RB, N B and GB of the selected instances used for the prediction are far from the properties RB 0 , N B 0 and GB 0 of the instance to predict.

In the general case, for an instance that does not belong to the database we do not know the values of the RB 0 , N B 0 and GB 0 properties as B 0 is unknown. However, the method can provide the user with some useful information that can be harnessed to determine if the instance will be difficult to solve or not and consequently use the appropriate set of parameters for a software. We need to find more accurate measures to estimate the difficulty of the problem as for example the zilla instance (500 taxa, 759 residues) has a RN B of 1.8 which makes it a problem easy to solve while it is not the case because the best optimal score known of 16218 is difficult to reach. However it is not too difficult to reach a score of 16221 or 16220 for example by using SAMPARS. The difficulty of zilla probably comes from the important number of taxa. We also need to refine the method to analyze the cases that lead to a bad prediction. This study nonetheless is very encouraging and is spurring us on to greater efforts.

Improvement of the evaluation of the objective function 6.1 Introduction

The objective or fitness function is one of the key elements for the successful implementation of metaheuristic algorithms because it is in charge of guiding the search process toward good solutions in a combinatorial search space. Previous works for solving the MP problem have commonly evaluated the quality of a potential solution using the parsimony cost depicted in equation (2.5) [START_REF] Andreatta | Heuristics for the phylogeny problem[END_REF][START_REF] Barker | LVB: Parsimony and simulated annealing in the search for phylogenetic trees[END_REF]Ribeiro andVianna, 2005, 2009;[START_REF] Richer | A memetic algorithm for phylogenetic reconstruction with maximum parsimony[END_REF]. In terms of performance the objective function is called every time a new solution must be evaluated and should have the shortest execution time.

In this chapter we introduce two implementations to compute the objective function of MP, the first one is the conventional way to compute it using the CPU (processor) and the second one is a GPU version. This objective function is the most time consuming [START_REF] Bader | ExactMP: An efficient parallel exact solver for phylogenetic tree reconstruction using maximum parsimony[END_REF] but exhibits much parallelism. Generally GPUs are used to perform floating point calculations but in this case we will only use integers and take advantage of the many cores / threads of the GPU.

Background

The evaluation function is used during different phases:

-when building an initial solution in a greedy manner, the objective function is used to evaluate the addition of a taxon

x to any node z of a subtree in order to minimize the increase in parsimony cost φ(z) (see 2.3); during this step, the algorithm counts the number of unions (or mutations) of the new hypothetical taxon, -when we generate a new neighbor using NNI, SPR or TBR, a subtree is degraphed and regraphed on another node of the initial tree, it is then necessary to recompute the score of the new tree obtained.

The techniques that we will talk about in this chapter concern what is called multi-character optimization techniques.

A set of methods [START_REF] Goloboff | Character optimization and calculation of tree lengths[END_REF][START_REF] Gladstein | Efficient incremental character optimization[END_REF][START_REF] Ronquist | Fast fitch-parsimony algorithms for large data sets[END_REF][START_REF] Yan | Fast character optimization in parsimony phylogeny reconstruction[END_REF] fall into the category of fast character optimization techniques, i.e. a set of shortcuts that helps decrease the computation time by not recalculating the whole tree each time a SPR or TBR modification is applied. Those techniques are particularly effective when an important number of SPR or TBR neighbors has to be evaluated. In the case of Fitch's parsimony, characters are considered as unordered and in C or C++ for example.

On the point of view of NVidia the CPU is considered as a multi-core device (with 2 to 20 cores 3 ) designed for general purpose treatments while the GPU is a many-core device (with hundreds or thousands of cores 4 ) designed for massively parallel treatments. However the extra cost of memory transfer from central memory (DDR) to the memory of the GPU (GDR) can sometimes lower down or annihilate the speed-up gained by the many cores.

A good introduction to learn GPU computing with CUDA is the book of Sanders and Kandrot (2010).

Programming model

We use the NVidia CUDA programming model (see Figure 6.2) where:

1. the data stored in the global memory (called DDR) which are processed by the CPU are copied to the RAM of the GPU (called GDR); the kernel is also copied to the GPU, 2. the data in GDR are treated by a kernel executed as many threads dispatched on the GPU cores and the data in the GDR memory are generally modified by the threads, 3. the data in the GDR is then copied back to the DDR to be eventually processed by the CPU.

GPU implementation

We now come to the presentation of the GPU implementation of φ(t). We use a static or flat representation of the tree to speed up the treatments and to obtain coalesced 5 memory transaction when all of the threads in a warp 6 or half-warp access global memory at the same time.

A node is represented by a structure of four integers, namely: L, R, P, N which are the index of the Left and Right offspring and the index of the Parent. The N field is used only for data alignment in memory but also contains the index of the node (see Figure 6.3).

The first n (here n = 3) nodes contain the leaves (for example: S1, S2, S3) that have no offspring so the fields L and R are set to -1. The rest of the n -1 nodes are internal nodes. The root node (here I2) has no parent so the field P is set to -1. The 3. for example the Intel Xeon E5 2670 of taurus, the cluster of the LERIA, has 10 cores with Hyper-Threading. 4. for example the NVidia GTX Titan X with 3072 CUDA cores. 5. coalesced: a NVidia term for contiguous. 6. warp: a group of 32 threads which is the minimum size of the data processed in SIMD fashion by a CUDA multiprocessor.

S1 S2 I1 S3 I2 L R P N S1 S2 S3 I1 I2 -1 -1 0 1 2 3 4 3 -1 -1 -1 -1 0 1 3 2 3 4 4 -1 Figure 6.
3 -Matrix representation of a binary tree (flat tree) with three leaves and two internal nodes two basic operations on the tree which are the degraph (cut of a subtree) and the regraph (paste of a subtree on a branch) must be implemented carefully in order to avoid a false interpretation of the score. For example if we put I2 before I1 on the table of Figure 6.3, the result will be wrong because I1 appears before I2 in the tree when we go from the leaves to the root.

We store the initial and the hypothetical sequences in a matrix M of (2n -1) × k characters or integers called data on the GPU. The first n rows will not be modified as they contain the initial sequences (OTUs), the other rows will be modified by the kernel to compute the hypothetical sequences or HTUs (see Figure 6.4). The code of the kernel consists in computing each internal sequence by assigning to each thread a column of the matrix M . We start from the first internal sequence at index n and then we move back towards the root (see Figure 6.5). The input parameters are N (= n) the number of leaves of the tree, K (= k) the length of the sequences and data (previously described). The output parameter is an array of k integers called mutations which will record the number of mutations for each column of the sequences. This array is stored on the GPU. The flat tree is obtained from constant memory (variable gpu_tree) for efficiency reason as the tree is a read-only structure.

The call to the kernel is presented on Figure 6.6. First, the tree to evaluate called cpu_tree is copied to the constant memory area on the GPU gpu_tree. After the execution of the kernel, the mutations are copied from the GPU to the CPU and are then summed on the CPU to obtain the parsimony score φ(t).

for large instances. This is also the case with the Tesla K40 (results not reported here). Indeed the Tesla K40 offers 40 percent higher performance than its predecessor, the Tesla K20 and has a boost technology analogous to the Turbo Boost seen in Intel processors. We can also note that it takes more time for Tesla K20 to compute the parsimony score than on the GTX 770 for small lengths (n = 64, k ≃ 1024). This is probably due to the occupancy factor which is small for small lengths. The occupancy is a metric for quantifying the ability of a given kernel to supply enough parallel work to achieve reasonable performance.

Conclusion

In this chapter we have presented the implementation details of the evaluation of the parsimony score of a tree on a CPU and a GPU architectures. The results obtained show that for a small number of sequences of short length the CPU is faster than the GPU. But for an important number of sequences with a long length the GPU becomes much faster than the CPU. This will be very interesting for phylogenies based on multi-genes or whole genomes where sequences can have from thousands to millions of residues.

Another area of research would be to evaluate many regraphs at the same time on the GPU on the same principle of (Robilliard, Marion-Poty, and F., 2008): instead of having a parallel version of the objective function we could have a parallel version of the SPR or TBR neighborhoods on the GPU. This work needs some amount of time to study the feasibility of the solution, design the algorithm and implement the code.

Parallelism can also be used with multi-core CPUs for example when we use the NNI, SPR or TBR neighborhoods using OpenMP 7 . It enables to evaluate different regraphs of subtrees in parallel.

7. This was a project of students of the Master 2 Professional of the Faculty of Science of Angers in 2015.

Conclusion and perspectives 7.1 Analysis of the work

In this thesis we have tried to address the Maximum Parsimony problem by developing new methods and ideas to solve it:

• We have designed a new Simulated Annealing algorithm much more interesting and accurate than the LVB implementation: the resolution time is not erratic like in the case of LVB and we provide solutions of better quality. Although relatively slow (compared to TNT) the software SAMPARS could benefit from an implementation close to the one of TNT, but the source code of TNT is not available. Recently, Pablo Goloboff the author of TNT has released Oblong [START_REF] Goloboff | Oblong, a program to analyse phylogenomic data sets with millions of characters, requiring negligible amounts of ram[END_REF] which source code is available and is close to the one of TNT. We are currently designing a new implementation to handle trees in a static way which should be more efficient than the present implementation and closer to the one of Oblong.

• From the experiments we have carried out the Path-Relinking approach does not seem to generate improvements most of the time. However the bottom-up implementation we have designed could be used as a measure of distance between trees. It is used in the Arbalet 1 Java software than can compare trees given in newick format. We now offer the possibility to the user to print the tree obtained by SAMPARS in newick format or to visualize it calling Arbalet.

• The prediction of the score of some instance is a very interesting tool that enables the experimenter or the software designer to decide what kind of parameters he has to use in order to reach a good solution in a reasonable time based on the measure of difficulty of the problem. Much work remains to be done to sharpen the accuracy of the prediction in particular for the instances that are predicted incorrectly at this time.

• The GPU implementation could be used to help decrease the computation time for large instances and could be selected instead of an assembly function. The GPU integration needs to be finalized. An interesting alternative to tackle the problem of long sequences is Oblong [START_REF] Goloboff | Oblong, a program to analyse phylogenomic data sets with millions of characters, requiring negligible amounts of ram[END_REF] which was designed for parsimony analysis of long sequences storing the data into the disk as blocks. However as pointed out by the author, using disk files to store the data slows down searches by a factor 4 to 5 times.

Our ideas have been integrated into the biosbl library 2 , a set of classes and programs used to handle input/outputs of bioinformatics data files and basic operations as multiple alignment and phylogenetic reconstruction with distance methods and Maximum Parsimony (see Appendix A).

The Maximum Parsimony problem is a very interesting optimization problem. What makes it all the more interesting, but also difficult to solve is the fact that the solutions of the problem are not vectors of values but trees. The static representation 1. www.info.univ-angers.fr/pub/richer/ur.php?arbalet 2. https://sourceforge.net/projects/biosbl/ CHAPTER 7. CONCLUSION AND PERSPECTIVES of trees is supposed to decrease the computation time as the tree can reside in the L2 cache of the processor and that the data are contiguous in memory. However the dynamic representation with pointers to left and right sons for binary trees is easier to handle. A trade-off between efficiency and ease of use must be found. All those implementation details, although important, were not taking into account as we have decided to write a generic code to test our ideas, especially to represent trees and to write the metaheuristic algorithms.

Further reflection 7.2.1 Reduction of the problem

Another approach, not developed here, and that could be very useful is a kind of multi-level or reduction resolution. This idea is not implemented in TNT but TNT uses the notion of cluster, where a cluster is a subtree of close taxa. TNT tries to optimize clusters, generally at the end of the search, but it would also be interesting to start the search by replacing a cluster by a unique sequence which is the taxa of the root of the cluster and solve the modified version of the problem during some iterations. As the number of taxa has an important influence on the resolution time it could enable us to decrease the overall computation time. When we reach some local optimum we can decide to come back to the original version of the problem by replacing the root sequence by the subtree that represents the cluster.

Operator to escape from local optimum

The methods used to escape a local optimum are the following:

-use a different neighborhood, -change the objective function (like the Ratchet technique), -perturb the current solution, -start the search from another initial solution.

However the perturbation of a solution is sometimes too simple and after a few moves we come back to the previous local optimum. It would be interesting to detect nodes that are potentially cause of a high contribution to the score of the tree and that could be degraphed and processed to build a new tree far from the last local optimum.

Predictor

The predictor could be integrated to many software or as a web service in order to help determine which set of parameters is suitable for the resolution of a particular instance. For the moment the database of instances is integrated in the code of biosbl library as a C array. It would be interesting to make it a database available by Internet and constitute a set of hard instances that could help test the efficiency and accuracy of Maximum Parsimony solvers.

Another interesting study would be to determine if the RN , RG and N G values computed for all instances are correlated. In other words, knowing RN , RG and N G can we infer the values of RB, N B and GB. If it is the case we could increase the accuracy of the predictor.

Multi-agent resolution

As the biosbl package now integrates the work of [START_REF] Goëffon | Nouvelles Heuristiques de Voisinage et Mémétiques pour le Probléme Maximum de Parcimonie[END_REF]; [START_REF] Goëffon | Progressive tree neighborhood applied to the maximum parsimony problem[END_REF] and my own work it could be possible to start in parallel the resolution of a problem using different metaheuristics (Simulated Annealing, Memetic Algorithm, Iterated Local Search, ...) and let the programs cooperate and exchange information during the search. For example when a new improving solution is found it could be broadcasted to the other software and they could take into account this solution in the next iteration of the search.

A Software A.1 Introduction

The ideas and concepts we have developed have been integrated in the biosbl library1 , a set of classes used to handle input/outputs of bioinformatics data files and basic operations like multiple alignment and phylogenetic reconstruction with distance methods and maximum parsimony.

We have designed two tree implementations: a dynamic one that is based on pointers and a static one based on an array representation. Those implementations give rise to different versions of a software. It should probably be more clever and less troublesome for the user to design a unique version which hides the implementation details and to allow the user to chose the implementation using command line parameters. However this is the choice of the last version of the implementation (version 3.1).

The different problems used to create the database are under the benchmarks/phylo directory. The files have been stored under different formats: fasta, phylip, nexus or sequence group.

A.2 Programs

A.2.1 Sampars SAMPARS, which stands for Simulated Annealing for Maximum PARSimony is the main software that helps solve the Maximum Parsimony problem. There are two versions following the tree implementation: sampar_dynamic.exe and sampar_static.exe.

The set of command line parameters can be obtained using the -h or --help option: where value=random, nj, greedy method used to generate initial solution --initial-temperature=FLOAT within [1..7] initial temperature of simulated annealing, default 0 means it will be evaluated --input-file=STRING or -i STRING input file (same as --taxa-file), use @file_name to look for the file in the current directory and subdirectories --markov-chain-length=NATURAL maximum number of Markov chain iterations (see also --effort) --max-reheat=NATURAL maximum number of times a reheat is performed --max-threads=NATURAL maximum number of threads used by OpenMP, use this parameter if you are using the 'spr_omp' neighborhood --neighborhood=VALUE where value=spr, spr_1, spr_omp, spr_pr, tbr neighborhood used by the local search algorithm --output=STRING or -o STRING output file, used to redirect output --output-newick=STRING or -n STRING print final tree into newick to file --random-generator-seed=NATURAL set random number generator seed --silent silent mode --stuck-limit=NATURAL number of iterations that causes a reheat if no improvement is found --taxa-file=STRING input file (same as --input-file), use @file_name to look TaxLabels chlorant amborell myristic hamameli canellac nymphaea tetracen caryophy sargento berberid winterac dioscore rafflesi himantan ceratoph gomorteg trichopo aristolo chenopod >convalla monimiac trochode piperace austroba hydnorac papavera cabombac dillenia saururac paeoniac sterculi magnolia degeneri lactorid illiciac eupomati phytolac ranuncul schizand lardizab nelumbon euptelea fumariac calycant annonace lauracea polygona menisper trimenia ; END; MATRIX a00000 AAAA?CAAAC?CAAAA??ACAA?AAACAAAAAAAACCCAAAA?CAAAAAAAACCAAAAA a00001 AAAA??AACAAC?AAAAAAAAAAAAAC?C?AA?AAAAAAAAAAA?AAAA????????AA ... a00048 AAAA?CAACAACAAAAAAAAAAAAAACC??AA?AACCAAAAAAA?AAAA????????AA

2. 1

 1 Three samples with two binary sequences S 1 and S 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Non exhaustive list of parsimony software . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Score of initial solution for different initialization methods for set 1 of real instances . . . . . . . . . . . . . . . 3.2 Parsimony score for instance tst12 with nj+greedy for f req = 5, . . . , 70 . . . . . . . . . . . . . . . . . . . . . . 3.3 Parsimony score for instance tst12 with nj+greedy for f req = 80, . . . , 600 . . . . . . . . . . . . . . . . . . . . . 3.4 Performance comparison for the reheat mechanism for problems tst10 to tst20 . . . . . . . . . . . . . . . . . . . 3.5 Performance comparison for the initialization methods for tst10 to tst20 . . . . . . . . . . . . . . . . . . . . . . 3.6 Input parameters of the SAMPARS algorithm and their selected values. . . . . . . . . . . . . . . . . . . . . . . 3.7 Mixed level covering array MCA(240; 4, 8, (2, 2, 2, 2, 3, 4, 4, 5)) representing an interaction test-suite for tuning SAMPARS (transposed). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8 Results from the 5 best parameter test cases in the tuning experiment. . . . . . . . . . . . . . . . . . . . . . . . 3.9 Comparison between SAMPARS and LVB for a subset of six selected instances. All tests were significant at the 0.05 level . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.10 Erratic results of LVB on instance tst04 for 5 executions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.11 Performance comparison among SAMPARS, GA+PR+LS, TNT and Hydra over 20 standard benchmark instances. 4.1 Sets of subtrees in the source and guiding trees ordered by number of leaves . . . . . . . . . . . . . . . . . . . . 4.2 Percentage of common subtrees for Parsimony trees of zilla used for Path-Relinking . . . . . . . . . . . . . . . . 4.3 Results of Path-Relinking for the different implementations (Times in seconds) . . . . . . . . . . . . . . . . . . 5.1 Multiple linear regression model for set2, E(K, L, R, N, G) = 65.65 -1.13 × K -0.52 × L -0.17 × R + 0.09 × N + 1.08 × G . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Properties of the problems of the knowledge base . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Data for set1 and set2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.6 Prediction without best score but with same criteria as in section 5.4.2.1 . . . . . . . . . . . . . . . . . . . . . . 5.7 Prediction without best score but with an increased selectivity . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.8 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1 Results in seconds for different implementations and number of threads for 64 bits architecture on an Intel i5-4570 CPU @ 3.20GHz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 Results in seconds for different number of taxa (n) and lengths (k) for 64 bits architecture . . . . . . . . . . . . .

Figure 2

 2 Figure 2.1 -Example of WPGMA Method

Figure 2 . 2 -

 22 Figure 2.2 -Example of Neighbor-Joinig Method

Figure 2 . 3 -

 23 Figure 2.3 -Characteristics of a binary tree with 5 leaves and 4 internal nodes

Figure 2 . 4 -

 24 Figure 2.4 -First pass and second pass for a tree of score of 3 under Fitch's optimality criterion

Figure 2 . 6 - 3 Figure 2 . 7 -Figure 2 . 8 -

 2632728 Figure 2.6 -Initial tree from which the right subtree (A,C) will be degraphed and regraphed on the left subtree (A,C)

Figure 2 . 9 -

 29 Figure 2.9 -Example of the MP problem representation of in a hypercube

Figure 2 .Figure 2 .

 22 Figure 2.10 -Steiner problem in a hypercube dimension 4

Figure 2 .

 2 Figure 2.11 -Total number of tree topologies with root in function of n species.

Figure 2 . 12 -

 212 Figure 2.12 -Example of neighborhood NNI

Figure 2 .Figure 2 .

 22 Figure 2.13 -Example of neighborhood SPR

Figure 2 . 15 -

 215 Figure 2.14 -Example of neighborhood TBR

  18 return s ⋆ known in the literature for a set of instances. They could improve 2 best-known solutions and match the results of 13 instances over 20.The most recent work to solve MP with the simulated annealing algorithm was SAMPARS5 implemented by[START_REF] Richer | Improvevement of fitch function for maximum parsimony in phylogenetic reconstruction with intel avx2 assembler instructions[END_REF] as the combination of the work of[START_REF] Richer | A memetic algorithm for phylogenetic reconstruction with maximum parsimony[END_REF] and Vázquez-Ortiz and Rodríguez-Tello (2011). SAMPARS provides the best results reported for the set of instances used in this work.

Algorithm 5 :

 5 Tabu Search Algorithm 1 TabuSearch(s 0 , N ) input: s 0 : initial solution, N : neighborhood output: best solution found s

Algorithm 7 :

 7 Genetic Algorithm 1 GeneticAlgorithm(P , N , stop_condition) input: P : initial population, N : neighborhood output: best solution found s ⋆ 2 generation ← 1 3 while not(stop_condition) do 4 (f ather, mother) ← ParentSelection(P ) 5 child ← Crossover(f ather, mother) 6 child ← Mutation(child) 7 Replace(child, P ) 8 generation ← generation + 1 9 s ⋆ ← best of P 10 return s ⋆

Figure 2 . 17 -

 217 Figure 2.17 -The three possible quartets for the species set {a, b, c, d}.

Figure 3

 3 Figure 3.1 -NNI, SPR and TBR neighborhoods

Algorithm 10 : 4 5

 104 GenerateNeighbor function input: s : current solution, iteration : iteration of Markov Chain, N : neighborhood output: best solution in the neighborhood s ′ 1 randomly select s ′ ∈ N (s) 2 if iteration is a multiple of 25 then 3 s ′ ← BestImproveDescent(s) return s ′

Table 3 . 4 -

 34 Performance comparison for the reheat mechanism for problems tst10 to tst20

Figure 4 . 2 -

 42 Figure 4.2 -Example of Bottom-Up Path-Relinking with source and guiding trees

Figure 5 . 1 -

 51 Figure 5.1 -Percentage of decrease between random (R), neighbor-joining (N), greedy(G) and best (B).

Figure 6 . 2 -

 62 Figure 6.2 -NVidia programming model

Figure 6 . 4 -

 64 Figure 6.4 -Data representation with sequences in one array and flat tree

5 -

 5 NAME sampars_dynamic.exe -Simulated Annealing for Maximum Parsimony SYNOPSIS sampars_dynamic.exe [OPTIONS] ... DESCRIPTION --alphabet=VALUE or -a VALUE where value=na, aa alphabet of input taxa file, use 'na' for nucleic acid and 'aa' for amino acid --cooling-factor=FLOAT within [0.96..0.9999] cooling factor (alpha) of simulated annealing, default value is 0.99 --descent-type=VALUE where value=strict, loose type of descent --effort=VALUE where value=small, medium, large will determine the number of iterations of the markov chain (see also --max-mc-iterations) --final-temperature=FLOAT within [1e-06..3.8] final temperature of simulated annealing, default value is 1e--heating-factor=FLOAT within [1.00001..1.8] reheat factor (beta) of simulated annealing, default value is 1.4 --help or -h help flag, print synopsis --improvement-frequency=NATURAL number of step in the Markov chain after which a descent is performed to improve neighbor --initial-solution=VALUE

75 D 26 H 31 L 33 N 28 O 37 T 17 V 29 W

 7526313328371729 Unified Device Architecture . . . . . . . . . . . . Descent, Branch-swapping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 E effort . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 F FPGA, Field-Programmable Gate Arrays . . . . . . . . . . . . . . . . . 36 G GA, Genetic Algortihm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 gap . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 GPU, Graphics Processing Unit . . . . . . . . . . . . . . . . . . . . . . . . . 75 GRASP, Greedy Randomized Adaptive Search Procedure . . 33 Greedy Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . HTU, Hypothetical Taxonomic Unit . . . . . . . . . . . . . . . . . . . . . 20 I ILS, Iterated Local Search . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Local Search . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 M MA, Memetic Algortihm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . NJ, Neighbor-Joining . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 NNI, Nearest Neighbor Interchange . . . . . . . . . . . . . . . . . . . . . . OTU, Operational Taxonomic Unit . . . . . . . . . . . . . . . . . . . . . . 19 P parsimony, optimality criteria . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 PN, Progressive Neighborhood . . . . . . . . . . . . . . . . . . . . . . . . . . 30 primary structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 R Ratchet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 S SA, effort . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 SA, Simulated Annealing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 secondary structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 SPR, Subtree Pruning Regrafting . . . . . . . . . . . . . . . . . . . . . . . . 28 SS, Sectorial Search . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . TBR, Tree Bisection Reconnection . . . . . . . . . . . . . . . . . . . . . . 28 Tertiary structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 TF, Tree Fusing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 TNT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 Topological distance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 TS, Tabu Search . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 U UPGMA, Unweighted Pair Group Method Using Arithmetic Average . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . VNS, Variable neighborhood search . . . . . . . . . . . . . . . . . . . . . Wagner tree . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 WPGMA, Weighted Pair Group Method with Averaging . . . 17

table 2

 2

		.1).
	Sample 1	Sample 2	Sample 3
	S 1 0	S 1 00	S 1 00
	S 2 1	S 2 01	S 2 11

Table 2 .

 2 1 -Three samples with two binary sequences S 1 and S 2

  were the first to establish a link between the MP problem and the

	Steiner problem
	Back to the example proposed by Goëffon (2006), we consider five sequences of length 4 :
	Sequence 1	0000
	Sequence 2	0011
	Sequence 3	0101
	Sequence 4	0110
	Sequence 5	1001

  or genetic local search algorithms (Ulder, Aarts, Bandelt, van Laarhoven, and Pesch, 1991) (see Algorithm 8).

	Algorithm 8: Memetic Algorithm

1 MemeticAlgorithm(P , N , stop_condition) input: P : initial population, N : neighborhood output: best solution found s ⋆ 2 generation ← 1 3 while not(stop_condition) do 4 (f ather, mother) ← ParentSelection(P ) 5 child ← Crossover(f ather, mother) 6 child ← Mutation(child) 7 child ← LocalSearch(child) 8 Replace(child, P ) 9 generation ← generation + 1 10 s ⋆ ← best of P 11 return s ⋆

6 .

 6 Parsimony software list: http://evolution.genetics.washington.edu/phylip/software.html# Parsimony

	Name of Software		Author	Year Metaheuristic
	PHYLIP			J. Felsenstein	1980 greedy + descent with NNI
	POY			W. Wheeler	1997 RAS + TBR
	PAUP Phylogenetic Analysis	D. L. Swofford	2003 many new algorithms
	Using Parsimony
	LVB	(Ludwig	Van	D. Barker	2003 Simulated Annealing
	Beethoven)
	GRASP+VND		Ribeiro and Vianna	2005 GRASP
	TABUPARS		Lin et al.	2007 Tabu Search
	TNT			P. Goloboff	2008 RAS + SPR/TBR, SS,
					Ratchet, TF, TD
	Hydra			Goëffon, Richer, Hao	2009 Memetic Algorithm
	XMP			White and Holland	2011 Branch and Bound
	SAMPARS		Vazquez-Ortiz, Richer, Le-	2013 Simulated Annealing
				saint, Rodriguez-Tello

Table 2 . 2

 22

-Non exhaustive list of parsimony software

Table 3 .

 3 1 -Score of initial solution for different initialization methods for set 1 of real instances

	problem	random	upgma	nj	greedy	nj + greedy	nj + f itch	best
	ANGI	425.80	240.00	234.00	229.42	228.82	226.00	226.00
	CARP	1546.92	628.00	592.00	591.18	591.68	597.00	591.18
	ETHE	898.58	399.00	385.00	388.84	386.32	382.00	382.00
	GOLO	903.82	567.00	535.00	536.32	534.16	533.00	533.00
	GRIS	476.50	200.00	195.00	182.22	184.48	186.00	182.22
	ROPA	731.28	363.00	353.00	345.72	346.00	347.00	345.72
	SCHU	2627.02	880.00	842.00	803.22	812.64	805.00	803.22
	TENU	1604.18	748.00	728.00	710.40	712.00	719.00	710.40

Table 3 .

 3 2 -Parsimony score for instance tst12 with nj+greedy for f req = 5, . . . , 70

	frequency	5	10	20	30	40	50	60	70
	time (min)	72	31	17	13	10	9	8	8
	average	1227.70	1226.50	1219	1219.90	1218.70	1218.50	1216.10	1214
	minimum	1222	1222	1213	1215	1213	1210	1211	1208

range from 5 iterations or within 10 to 100 iterations (by step of 10) or 200 to 600 iterations (by step of 100).

Table 3 .

 3 [START_REF] Aarts | Statistical cooling: A general approach to combinatorial optimization problems[END_REF][START_REF] Van Laarhoven | Simulated annealing: Theory and applications[END_REF][START_REF] Abramson | Simulated annealing cooling schedules for the school timetabling problem[END_REF][START_REF] Rodriguez-Tello | An effective two-stage simulated annealing algorithm for the minimum linear arrangement problem[END_REF]. They can be divided into two main categories: static and 3 -Parsimony score for instance tst12 with nj+greedy for f req = 80, . . . , 600

	frequency	80	90	100	200	300	400	500	600
	time (min)	8	8	7	6	6	5	5	4
	average 1217.40 1218,30 1218.40 1219,30 1219,80 1217,40 1220,30 1220,10
	minimum	1215	1216	1215	1212	1216	1213	1213	1215

Table 3 .

 3

		5 -Performance comparison for the initialization methods for tst10 to tst20
	effort level			medium					large
	max reheats	0	1	2	3	4	0	1	2	3	4	total
	random	0	1	0	0	0	3	4	2	2	4	43
	greedy	0	1	0	0	1	1	4	1	4	3	47
	nj + f itch	0	1	0	1	1	0	3	2	1	4	48
	nj	1	1	2	0	0	2	3	3	2	6	63
	nj + greedy	0	2	0	2	0	2	3	5	5	3	66
	allmethods	0	0	0	0	0	0	5	1	2	5	156

Table 3 .

 3 6 -Input parameters of the SAMPARS algorithm and their selected values.

	N R	IS	T i	α	β	SC	M CL	N
	0 1	Random Greedy	6 √ n	0.99 1.4 0.97 1.6 # max. of solutions that don't improve Final temperature	Medium SPR Big TBR
	2	NJ	2.5

Table 3 .

 3 7 -Mixed level covering array MCA(240; 4, 8, (2, 2, 2, 2, 3, 4, 4, 5)) representing an interaction testsuite for tuning SAMPARS (transposed).

	1 0	2 2	3 4	4 3	6 2	7 3	8 0	9 2	10 1	11 3	12 3	13 0	14	15	16	17	18	19	20	21	22	23	24	25

Table 3 .

 3 8 -Results from the 5 best parameter test cases in the tuning experiment.

	Num.		Test case	Avg. parsimony	Avg. CPU time
	59		43001010		1004.06	3512.51
	233		40010010		1004.14	3511.35
	217		43000100		1005.00	2058.93
	37		30001000		1005.02	2047.52
	118		30300101		1005.29	3136.40
	Table 3.9 -Comparison between SAMPARS and LVB for a subset of six selected instances. All tests were
	significant at the 0.05 level
				LVB		SAMPARS
	Instance n k	B	Avg. Dev. T ime		B	Avg. Dev. T ime P -value
	tst01	45 61	549 553.87 2.47 85.57	545 545.83 0.87 295.80 1.80E-11
	tst02	47 151	1367 1375.33 4.77 23.63	1354 1356.13 1.33 479.50 5.08E-21
	tst03	49 111	840 850.83 4.86 68.02	833 834.00 1.05 577.12 1.38E-18
	tst08	57 119	867 879.80 5.01 922.61	852 854.53 2.37 665.50 2.12E-11
	tst09	59 93	1153 1160.77 4.60 58.53	1143 1145.50 1.11 719.24 3.48E-18
	tst10	60 71	727 738.00 5.59 570.62	720 721.27 0.78 500.28 1.89E-16
	Avg.		917.17 926.43 4.55 288.16	907.83 909.54 1.25 539.57

Table 3 .

 3 10 -Erratic results of LVB on instance tst04 for 5 executions

	time (s)	score	number of trees
	84	609	16
	13194 (>3h)	595	4747
	19532 (>5h)	600	30
	1960 (32m)	605	92
	307 (5m)	609	510

  al., 2012; Seridi et al., 2013; Zhou et al., 2013; Resende and Ribeiro, 2014).

Table 4 .

 4 2 -Percentage of common subtrees for Parsimony trees of zilla used for Path-Relinking

	source / guiding % common subtrees
	16218 b / 16218	74.75%
	16219 / 16218	73.15%
	16250 / 16218	61.12%
	16401 / 16218	45.69%
	16611 / 16218	36.27%
	21727 / 16218	0.00%
	16250 / 16219	61.72%
	16401 / 16219	48.70%
	16401 / 16250	45.49%
	16611 / 16250	44.89%
	21727 / 16250	0.20%
	16611 / 16401	36.67%
	21727 / 16611	0.20%
	Three different implementations were compared:

Table 4 .

 4 3 -Results of Path-Relinking for the different implementations (Times in seconds)

			Bottom-Up		Top-Down No Minimization	Top-Down With Minimization
	source / guiding	Trans. Time	#Equal	#Less	Trans. Time	#Equal	#Less	Trans.	Time	#Equal	#Less
	16218 b / 16218	24	0.10	13	0	118	0.33	16	0	74	0.85	21	0
	16219 / 16218	32	0.14	2	0	462	1.36	1	0	343	4.08	3	0
	16250 / 16218	97	0.40	3	0	1770	4.79	1	0	1519 25.32	2	0
	16401 / 16218	151	0.60	2	0	1891	5.16	1	0	1474 25.36	2	0
	16611 / 16218	186	0.75	2	0	1903	5.19	1	0	1225 19.57	2	0
	21727 / 16218	446	1.68	2	0	1881	5.17	1	0	1241 18.70	2	0
	16250 / 16219	92	0.37	2	0	1722	4.76	1	0	1418 24.70	1	0
	16401 / 16219	152	0.61	2	0	1867	6.17	1	0	1390 23.74	1	0
	16401 / 16250	162	0.63	1	0	1707	4.71	1	0	1203 19.74	1	0
	16611 / 16250	144	0.56	1	0	1746	4.79	1	0	1233 35.78	3	0
	21727 / 16250	449	1.71	1	0	1712	4.77	1	0	1196 19.65	1	0
	16611 / 16401	202	0.79	3	0	1602	4.42	1	0	1216 47.75	2	1
	21727 / 16611	455	1.84	2	0	1842	5.16	1	0	1528 22.40	593	586

table (

 (

		"data.rdata")
	> model <-lm(V7 ~V2 + V3 + V4 + V5 + V6, data)
	> summary(model)
	Call:
	lm(formula = V7 ~V2 + V3 + V4 + V5 + V6, data = data)
	Residuals:
	Min	1Q Median	3Q	Max
	-9.1111 -1.7844 -0.6903 2.1100 10.1303
	Coefficients:
	Estimate Std. Error t value Pr(>|t|)
	(Intercept) 65.65223	10.49157	6.258 2.10e-05 ***
	V2	-1.13401	0.31890 -3.556 0.003163 **
	V3	-0.52265	0.11382 -4.592 0.000419 ***
	V4	-0.16641	0.11465 -1.451 0.168685
	V5	0.09717	0.13636	0.713 0.487797
	V6	1.08315	0.19320	5.606 6.47e-05 ***
	---
	Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 ' ' 1

Table 5 .

 5 × R + 0.09 × N + 1.08 × G, where E is the estimation of the best score Y = B. The constant (65.65), K, L and G are reported to have a high significance. 1 -Multiple linear regression model for set2, E(K, L, R, N, G) = 65.65 -1.13 × K -0.52 × L -0.17 × R + 0.09 × N + 1.08 × G

	problem K tst01 45 61 724 588 580 545 L R N G B	E B -E 547.62 -2.62
	tst02	47 151 1629 1454 1422 1354 1343.87	10.13
	tst03	49 111 1064	908 894 833	831.57	1.43
	tst04	50 97 789 674 645 587	591.08	-4.08
	tst05	52 75 1012	845 838 789	788.86	0.14
	tst06	54 65 814 651 642 596	593.62	2.38
	tst07	56 143 1586 1401 1365 1269 1278.11	-9.11
	tst08	57 119 1108	951 929 852	853.08	-1.08
	tst09	59 93 1450 1216 1215 1144 1143.02	0.98
	tst10	60 71 994 799 785 720	722.99	-2.99
	tst11	62 63 772 611 598 540	541.04	-1.04
	tst12	64 147 1547 1357 1309 1208 1208.50	-0.50
	tst13	65 113 1917 1629 1606 1515 1511.69	3.31
	tst14	67 99 1519 1257 1247 1160 1157.98	2.02
	tst15	69 77 1066	832 833 752	752.87	-0.87
	tst16	70 69 778 614 594 527	523.79	3.21
	tst17	71 159 2975 2597 2578 2450 2451.67	-1.67
	tst18	73 117 1956 1656 1630 1521 1522.66	-1.66
	tst19	74 95 1380 1126 1110 1012 1014.14	-2.14
	tst20	75 79 956 774 738 659	654.79	4.21

Table 5 .

 5 set #instances K min K max L min L max S min S max 2 -Properties of the problems of the knowledge base

	set1	8	47	117	59	179	71	88
	set2	20	45	75	61	159	68	88
	set3	11	116	299	355	8245	64	100
	set4	20	13	158	411 120762	67	97
	set5	14	10	36	64 10539	56	96
	set6	29	49	500	36	759	62	91
	set7	12126	10	400	147 152160	37	97
	set8	2	144	263 8619	10305	74	87
	set9	40	50	200 1287	1287	74	94

  table but it seems less obvious to infer the difficulty of an instance with this quantity.

	problem	K	L	R	N	G	B S RN RG RB RNB NB
	ANGI	49 59 421 234 229 216 77 44 46 49	9	5
	CARP	117 110 1520	592 590 548 82 61 61 64	5	3
	ETHE	58 86 878 385 388 372 71 56 56 58	3	2
	GOLO	77 97 883 535 535 497 84 39 39 44	10	5
	GRIS	47 93 470 195 183 172 88 59 61 63	8	4
	ROPA	75 82 719 353 347 326 82 51 52 55	7	4
	SCHU	113 146 2559	810 803 759 74 68 69 70	3	2
	TENU	56 179 1568	728 710 682 76 54 55 57	5	3
	tst01	45 61 723 594 586 545 68 18 19 25	28	7
	tst02	47 151 1631 1465 1426 1354 75 10 13 17	40	7
	tst03	49 111 1067	899 895 833 82 16 16 22	28	6
	tst04	50 97 787 663 645 587 88 16 18 25	38	9
	tst05	52 75 1015	855 844 789 68 16 17 22	29	6
	tst06	54 65 816 666 649 596 75 18 20 27	32	9
	tst07	56 143 1590 1401 1364 1269 82 12 14 20	41	8
	tst08	57 119 1114	945 929 852 87 15 17 24	35	9
	tst09	59 93 1447 1235 1216 1158 68 15 16 20	27	5
	tst10	60 71 991 795 787 734 75 20 21 26	24	6
	tst11	62 63 771 620 601 552 82 20 22 28	31	8
	tst12	64 147 1546 1357 1310 1245 87 12 15 19	37	7
	tst13	65 113 1912 1617 1608 1535 68 15 16 20	22	7
	tst14	67 99 1524 1266 1252 1174 76 17 18 23	26	6
	tst15	69 77 1062	839 834 768 82 21 21 28	24	7
	tst16	70 69 776 617 600 555 88 20 23 28	28	8
	tst17	71 159 2970 2588 2581 2475 68 13 13 17	23 10
	tst18	73 117 1962 1666 1637 1551 76 15 17 21	28	6
	tst19	74 95 1382 1124 1112 1041 82 19 20 25	24	6
	tst20	75 79 955 769 741 688 88 19 22 28	30	9
				Table 5.3 -Data for set1 and set2

Table 5 . 5

 55 2. multiple linear regression model computation9. RNB: note that RN B is different from N B

	% of error	within	outside	not		error
		range	range predicted min max	avg
	1%	8,307 (67.70%)	3960	3	0	4260 115.83
	2%	9,754 (79.49%)	2513	3	0	8755 140.76
	3%	10,427 (84.98%) 1840	3	0	8755 155.25
	4%	10,804 (88.05%) 1463	3	0	13278 164.84
	5%	11,056 (90.11%) 1211	3	0	13278 173.22
	10%	11,604 (94.57%)	663	3	0	41082 199.01

Table 5 .

 5

	6)

Table 5 .

 5 6 -Prediction without best score but with same criteria as in section 5.4.2.1

  .7).

	% of error	within	outside	not		error
		range	range predicted min max	avg
	1%	10627 (86.61%) 1635	8	0	4279	99.81
	2%	11543 (94.07%)	719	8	0	7289 134.07
	3%	11782 (96.02%)	480	8	0	18043 147.17
	4%	11909 (97.06%)	353	8	0	18043 156.69
	5%	11985 (97.68%)	277	8	0	18043 164.95
	10%	12109 (98.69%)	153	3	0	27572 175.87
	Table 5.7 -Prediction without best score but with an increased selectivity
	5.5.1.3 Enhanced problem selection

  Here is a detailed output of the prediction for instance HIV of best score 73, 537. We could predict a score of 73, 500.17 which corresponds to an error of 36.83, i.e. 0.05% of the score. The closest instances selected are at least 207. They were sorted according to the RN GS distance and only the 10 closest to HIV were kept.

	for HIV, K=94.00, L=10922.00 RN=41.83, RG=41.83, S=79.00
	selection criterion RN : [41.33..42.33]
	selection criterion RG : [41.33..42.33]
	selection criterion S : [69.00..89.00]
	instances=207
	....
	selected
	problem;	K ;	L ;	R ;	N ;	G ;	B ; S ;	RB ;	NB ;	GB ;	RNGS
			error		within		outside	not		error
					range		range predicted min max	avg
			1%	11808 (96.23%)	454	8	0.00 90.81 4579.12
			2%	12194 (99.38%)	68	8	0.00 96.74 4679.34
			3%	12221 (99.60%)	41	8	0.00 96.94 4679.34
			4%	12229 (99.67%)	33	8	0.00 96.91 4679.34
			5%	12238 (99.74%)	24	8	0.00 96.85 4679.34
			10%	12262 (99.93%)	0	8	0.00 96.76 4679.34
							Table 5.8
	5.5.3.1 Example that works: HIV

Table 6 .

 6

		2 -Results in seconds for different number of taxa (n) and lengths (k) for 64 bits architecture
	taxa (n)	length (k)	i5 4570	GTX 770	Tesla K20	speed up
			8 bits data -DNA, RNA
	64	1024	0.03	1.43	1.89	-
		32768	1.59	5.41	5.07	-
		100000	12.50	15.01	11.81	1.05
		200000	36.80	28.90	22.45	1.63
	512	1024	0.42	6.14	9.56	-
		32768	38.70	26.36	26.76	1.44
		100000	125.26	82.50	68.49	1.82
		200000	252.39	161.75	131.77	1.91
			32 bits data, PROTEINS
	64	1024	0.17	1.43	1.86	-
		32768	17.48	6.78	6.71	-
		100000	58.83	20.33	18.10	3.25
		200000	120.76	39.89	35.35	3.41
	512	1024	1.69	6.55	11.05	-
		32768	162.60	38.23	43.09	3.77
		100000	508.90	127.55	122.06	4.15
		200000	1022.30	261.94	243.17	4.2

http://www.info.univ-angers.fr/pub/richer/ur.php?arbalet

Note that throughout the thesis we used N or n for the number of sequences, but in this case N will be used to identify the neighbor-joining score.

int fitch(char * x, char * y, char * z, int k) {

https://sourceforge.net/projects/biosbl/

Acknowledgments

and technical and administrative support of the LERIA laboratory of the University of Angers.

List of Figures

based on SPR and TBR, TNT can visit millions of trees in a very short time as it is based on Gladstein's incremental down-pass optimization [START_REF] Gladstein | Efficient incremental character optimization[END_REF].

The results from this experiment are depicted in Table 3.11. Columns 1 to 3 indicate the instance and its size in terms of taxa (n) and length (k). The best solutions found by GA+PR+LS , TNT and Hydra, in terms of parsimony cost Φ are listed in the next three columns. Note that for TNT this corresponds to the best result over 30 runs. Columns 7 to 9 present the best (B), average (Avg.) and standard deviation (Dev.) of the parsimony cost attained by SAMPARS in 30 independent executions, as well as its average CPU time in seconds (Column 10). Finally, the difference (δ) between the best result produced by the SAMPARS algorithm and the best-known solution produced by either GA+PR+LS or Hydra is shown in the last column. The analysis of the data presented in Table 3.11 lead us to the following observations. First, we clearly see that the procedure GA+PR+LS [START_REF] Ribeiro | A hybrid genetic algorithm for the phylogeny problem using path-relinking as a progressive crossover strategy[END_REF] consistently returns poorer quality solutions than Hydra and SAMPARS. Second, the best solutions attained by the proposed SAMPARS algorithm are very competitive with respect to that produced by the existing state-of-the-art procedure Hydra [START_REF] Goëffon | Nouvelles Heuristiques de Voisinage et Mémétiques pour le Probléme Maximum de Parcimonie[END_REF], since on average SAMPARS provides solutions whose parsimony costs are smaller (compare columns 6 and 7). In fact, it is able to improve 9 previous best-known solutions produced by Hydra and to reach the results of the other 11 benchmark instances.

Thus, as this experiment confirms, our SAMPARS algorithm is an effective alternative for solving the MP problem, compared with the three representative state-of-art algorithms: GA+PR+LS, TNT and Hydra in terms of the quality of the solution.

Conclusions

In this chapter we have presented an improved Simulated Annealing algorithm called SAMPARS to find solutions for the MP problem under the optimality criterion of Fitch. SAMPARS's components and parameter values were carefully determined, Algorithm 11: Path Relinking from [START_REF] Ribeiro | A hybrid genetic algorithm for the phylogeny problem using path-relinking as a progressive crossover strategy[END_REF] 1 PathRelinking(s, t) input: s, t : trees such that f itness(s) >= f itness(t) output: best solution found s ⋆ 2 s ⋆ ← s (by definition)

The implementation of [START_REF] Ribeiro | A hybrid genetic algorithm for the phylogeny problem using path-relinking as a progressive crossover strategy[END_REF] can be called top-down: it starts from the root of the tree and recursively explores each left and right subtree. For each iteration they compare the taxa in the left and right subtrees of the source and guiding solutions. All taxa of the source left (resp. right) subtree that are in the right (resp. left) subtree of the guiding solution are moved to the right (resp. left) subtree of the source solution. The major drawback of this implementation is that it requires a lot of modifications and moves of the taxa from left to right (resp. right to left) subtrees. When a taxon is degraphed it is then regraphed on a branch of the sibling subtree that minimizes the parsimony score of the tree. We can also put the degraphed taxon at the top of the sibling subtree to avoid the search of a minimum tree. If we remove the exploration phase from the PR algorithm we then have an algorithm that transforms one tree into another which can be used to compare trees and define a measure of distance between them. This measure of distance seems interesting because it is closer to the topology of the trees to compare.

In the next section we describe the method that we have designed to perform Path-Relinking between two trees based on the difference of subtrees. Then we show the results and analysis of our experimentation.

The Bottom-up implementation

Ribeiro and Vianna (2009) implemented PR to solve the MP problem by using a recursive algorithm that starts from the root of the tree and compares the left and right subtrees of the source and guiding solutions. All taxa of the left subtree of the source that are present in the right subtree of the guiding solution are moved to the right subtree of the source solution and conversely (see Figure 4.1). The major drawback of this implementation is that it requires a lot of modifications and moves of the taxa from left (resp. right) to right (resp. left) subtrees.

In contrast, we have implemented a bottom-up iterative solution which compares the subtrees present in the source and guiding solutions (see Figure 4.2). For this the subtrees of each solution are ordered by their number of leaves and we start to compare subtrees of size 2, then subtrees of size 3, and so on (see Algorithm 12).

When a subtree of the guiding solution t = (X, Y ) is not found in the source solution, we have to transform the subtree in the source solution into the one in the guiding solution. Here X is the left subtree of t and Y its right subtree.

Algorithm 12: Path-Relinking with bottom-up iterative implementation input: s: source tree, g : guiding tree output: number of transformations 1 reorder(s) ; 2 reorder(g) ; 3 transf ormations ← 0; 4 Ω g ← ordered set of subtrees of guiding tree g; 5 change ← true; 6 while change do 7 Ω s ← ordered set of subtrees of source tree s;

degraph Y and regraph on X in s; transf ormations ← transf ormations + 1; return transf ormations Table 4.1 -Sets of subtrees in the source and guiding trees ordered by number of leaves

#

Source Guiding

Consider the example of Figure 4.2 (a) where we can see the source and guiding trees. A preliminary modification of the trees (line 1 and 2 of the algorithm) consists in reordering the trees by the lexicographic order of the leaves in order to be able to efficiently compare the subtrees. In fact the subtrees (A, B) and (B, A) are equal and we do not want to have to check both cases, so we will only allow subtrees of the form (A, B).

On each node, one of the sequences on the left subtree must be inferior to all the sequences on the right subtree. For example, the node (B, A) will be changed by swapping the leaves in order to obtain (A, B). The subtree (C, (A, B))) will be reordered as ((A, B), C) because on the right node, A and B are inferior to C. The subtree ((C, E), (A, F )) will be reordered as ((A, F ), (C, E)) because A that initially appears on the right subtree is inferior to C and E.

The subtrees present in the source and guiding solutions are represented in Table 4.1 in Newick notation. After reordering the source and guiding trees we can enter the main loop of the algorithm (lines 6 to 12).

The subtree (A, B) present in the guiding tree is not present in the source tree so we degraph B from the source tree, and regraph it on A (see Figure 4.2(b)). On Figure 4.2(b) the subtree (E, F ) present in the guiding tree is not present in the source tree so we degraph F from the source tree and regraph it on E.

On Figure 4.2(c) the subtree ((A, B), C) present in the guiding tree is not present in the source tree, we have to degraph C from the source tree and regraph it on (A, B).

The subtree (D, (E, F )) present in the guiding tree is not present in the source tree: degraph (E, F ) from the source tree and regraph it on D (see Figure 4.2(e)).

Finally the source and guiding trees are equal so we can stop the algorithm.

3. estimation using the model

Problems selection

We select from the database all problems P = (K, L, ..., GB) that are close to P 0 and for which the four following criteria are met:

Initially δ = 0.5 and we retain the first 10 problems that satisfy the constraints and we add P 0 to this set. If the selection does not provide 10 problems we increase δ by a factor of 0.5. In some cases we have a lack of problems close to P 0 .

Linear regression

We compute lm a multiple linear regression model taking into account the properties K, L, R, N , G, S, RN , RG, RB in order to determine B. For this we use the GSL (GNU Scientific Library) and especially the function gsl_multifit_linear. We can then determine the average error over all the problems.

Estimation

We compute E 0 which is the score estimated from the lm model using the data of P 0 .

Results and analysis

If we consider a prediction with an accuracy of 1%, i.e. the score we predict E is such that E = B ± 1%, then we get the following results (see Table 5.4):

-12,132 instances were predicted successfully, so that means 98.88% of the database -135 instances have an accuracy greater than 1% -3 problems could not be predicted because we could not get enough close instances to use the MLR. If we remove the instance to be predicted from the MLR then we get results of bad quality (see Table 5.5) with only 67.7% of the problems that could be predicted with an error of 1%. CHAPTER 6. IMPROVEMENT OF THE EVALUATION OF THE OBJECTIVE FUNCTION multistate: they can transform from one state to another independently. As a consequence, an unrooted tree may be rooted on any branch with no modification of the parsimony score, which means there is a potential root node for any branch.

% of error

In [START_REF] Goloboff | Character optimization and calculation of tree lengths[END_REF], the author proposed a method for indirect calculation of the parsimony score which uses two passes. This method needs only to compare the root of the clipped tree with the potential root of the target tree to obtain the score of a potential new tree for a SPR search. Gladstein (1997) also proposed an algorithm which is exact and correct. In [START_REF] Yan | Fast character optimization in parsimony phylogeny reconstruction[END_REF] a two passes algorithm is described which has the same complexity of Goloboff's and is faster than the incremental method of Gladstein.

However those fast character optimization techniques are difficult to implement and are effective with basic C implementation.

CPU implementation

In this section we briefly describe some implementation details to efficiently compute the parsimony score φ(t) of a tree using the SIMD (Single Instruction on Multiple Data) units of the CPU (Central Processing Unit). This version will serve as a reference for our benchmarks. The most time consuming function in a search algorithm for the resolution of MP in our case is the function of Fitch (see Algorithm 13).

Algorithm 13: Fitch algorithm to compute hypothetical sequence z from two sequences x and y and evaluate the number of mutations input: x, y: array[k] of sets of characters output: z: array[k] of sets of characters, mutations: number of mutations 1 mutations ← 0;

This function will be called to compute all the hypothetical sequences z for each node of the tree. The sum of all mutations will constitute the parsimony score of a tree.

In order to efficiently perform the union and intersection of Algorithm 13, each character is represented by a power of 2. For nucleotides the different symbols are (-, A, C, G, T, ?). We represent the gap symbolby 2 0 = 1, until T 2 4 = 16. The undefined character ? which can represent any other character is then coded by the value 31 = 1 + 2 + • • • + 16. With this representation the union can be performed by the binary-OR (|) and the intersection by the binary-AND (&). This gives rise to the C version of the previous algorithm (see Figure 6.1).

Although quite simple this function can not be vectorized by C++ compilers like GNU g++ or Intel icpc. Nevertheless, we can implement an efficient version of the Fitch function in assembly by taking full advantage of some relevant features offered by modern x86 processors. More precisely, the core of modern x86 processors has a SSE (SIMD Streaming Extension) or AVX (Advanced Vector Extensions) unit which enables to treat data as vectors. The use of vectors let us perform the same operation on different data at the same time. Recent Intel processors offer on a 32-bits architecture a set of 8 SSE registers of 128 bits or 8 AVX registers of 256 bits long. If we represent a nucleotide with one byte, in the case of DNA, then a SSE register can store and handle 16 bytes (nucleotides) at a time (resp. 32 with AVX). In the case of proteins it would be necessary to use 32 bits integers to represent the 20 different amino acids, so a SSE register would handle 4 integers (resp. 8 with AVX).

The experiments that were carried out by [START_REF] Richer | Three new techniques to improve phylogenetic reconstruction with maximum parsimony[END_REF] show that the vectorization of Fitch's function gives a 90% speedup on Intel Core 2 Duo processors compared to the basic version of Figure 6.1, while other architectures (Pentium II/III/4, Pentium-M, Athlon 64, Sempron) provide 70 to 80% improvement. This improvement enables to divide the overall computation time of 2 int mutations = 0; 3 for (int i=0; i<k; ++i

return mutations; 10 } a program by a factor of 3 to 4. A first pseudo-code was given in [START_REF] Ronquist | Fast fitch-parsimony algorithms for large data sets[END_REF] for PowerPC processors.

Finally, note that the most recent processors (Intel Core i5 or i7, AMD Phenom) introduce the SSE4.2 instructions set that contains the popcnt instruction which counts the number of bits set to one in a general purpose register. This instruction is used essentially to determine the number of mutations that occur when we perform the union between two SSE registers. By replacing a C implementation of popcnt by the native SSE4.2 assembly instruction, the experiments carried out show an overall improvement of 95% (on an Intel Core i7 860 processor) compared to the basic implementation. The introduction of the AVX and then AVX2 instructions sets in 2013 by Intel enables to use vectors of 256 bits. Rewriting the SSE4.2 version into AVX2 could only lead of a gain of 1% to 3% compared to the SSE4.2 version on Intel Haswell processors. For more details please see [START_REF] Richer | Improvevement of fitch function for maximum parsimony in phylogenetic reconstruction with intel avx2 assembler instructions[END_REF].

Some recent tests we made on a descent algorithm using the SPR neighborhood on a new implementation of trees based on an array representation gives the following results (see Table 6.1): by using SSE4.2 we can decrease the execution time by 96% compared to the C language implementation.

GPU programming

GPU programming also called GP (for General Purpose) GPU programming consists in using a graphics card to perform very intensive and parallel computations. The GPU (Graphics Processing Unit) can be viewed as a co-processor although this analogy is far from reality because the GPU can also be considered like a computer inside the computer as it is formed of cores and memory.

NVidia which sells GPUs since 1993 has developed CUDA 1 a parallel computing platform and application programming interface (API). The CUDA platform is designed to work with programming languages such as C, C++ and Fortran 2 and enables the programmer to write code that will be executed on the GPU, but that will be called from a program executed on the CPU. The interesting point is that we do not have to learn a new language because the code of the GPU, called a kernel, can be written 1.

Results

We have implemented a simple benchmark in order to compare the CPU and GPU implementations. The program takes as input the number of sequences and the number of residues of the sequences. The sequences are randomly generated. We then randomly generate 50 different trees and perform 1000 evaluations of the parsimony score for each tree. Results were obtained on the latest Intel CPU architecture (Haswell -Core i5-4570 running at 3.20 GHz with AVX2), a common GPU (GTX 770) and a high end GPU (Tesla K20).

From Table 6.2 we can see that depending of the number of sequences the GPU will become more efficient as long as the number of sequences increases. For 8 bits data, when n = 64, the GTX 770 is more efficient than the Core i5 when k is greater than millions residues. When n = 512 then it is also the case but with k ≃ 32768. The Tesla K20 seems particularly efficient for the file in the current directory and subdirectories --verbose-level=NATURAL or -v NATURAL verbose level, chose between 0 (silent) and 5, default is 1 --view=VALUE where value=text, arbalet type of visualization For example to solve problem tst01 of set2, we can use: build/obj/v3.1/release/gnu/programs/sampars_static.exe --input-file=benchmarks/phylo/set2/tst01.fasta

We can use the script bin/links.sh that will create symbolic links to the binaries. So we can after calling the script, write:

./sampars_static.exe --input-file=benchmarks/phylo/set2/tst01.fasta If we need to solve a hard instance we could use the TBR neighborhood with an important Markov Chain Length and using an initial solution generated by a greedy process:

./sampars_static.exe --input-file=benchmarks/phylo/set2/tst01.fasta --neighborhood=tbr --effort=large --initial-solution=greedy

A.2.2 Predictor

In order to predict the best score of some instance we can use the parsimony_predictor program with the following options:

--deduce=percentage which tries to deduce relationships with multiple linear regression from the problems of the database and takes into account predictions less than the given percentage knowing the best score. Use for example 1 for 1% accuracy, --predict=percentage which tries to predict the parsimony score of the instances of the database without knowing the optimal score of the instance to predict. It takes into account predictions less than the given percentage without knowing the best score. Use for example 1 for 1% accuracy, --export to export the database as a text file --input-file=file to try to predict the best score of some problem not in the database Finally the following programs are also available:

• fitch_score which given an instance and a newick tree computes the parsimony score of the tree • tree_comparator enables to compare two trees using Path-Relinking

A.3 File formats

We show here the different phylogenetic file formats for some instances.

A.3.1 Sequence group .sg

This file format is very simple and contains one taxon per line. The line starts with the name of the taxon followed by one or more spaces and the residues either in IUPAC 3 nomenclature or using numbers where 0=A, 1=C, ... austroba 0000?100100100000000000000100000000000000000?0000????????00 annonace 00010000000110001000000000100000000000000000000000010000000 myristic 0000?10000010000100000?001000000?00111000000100000001000010 ... convalla 11101000??1000001001000001000000000100100000100010000010000

A.3.2 Fasta .fasta

The fasta format identifies each taxon by a line that starts with the symbol > followed by the name of the taxon, a new line and the residues. >austroba aaaa?caacaacaaaaaaaaaaaaaacaaaaaaaaaaaaaaaaa?aaaa????????aa >annonace aaacaaaaaaaccaaacaaaaaaaaacaaaaaaaaaaaaaaaaaaaaaaaacaaaaaaa ... >convalla cccacaaa??caaaaacaacaaaaacaaaaaaaaacaacaaaaacaaacaaaaacaaaa

A.3.3 Phylip .phy

The Phylip format of Joe Felsenstein starts with a line made of two numbers: the number of taxa and the length of a taxon. It is followed by the taxa defined by their names and the residues. This file format is not suitable for taxa that have a long name because the residues must be found on column 11 of the line. The first 10 characters of each line contain the name of the taxon.

49 59 austroba AAAA?CAAAC?CAAAA??ACAA?AAACAAAAAAAACCCAAAA?CAAAAAAAACCAAAAA annonace AAAA??AACAAC?AAAAAAAAAAAAAC?C?AA?AAAAAAAAAAA?AAAA????????AA ... convalla AAAA?CAACAACAAAAAAAAAAAAAACC??AA?AACCAAAAAAA?AAAA????????AA A.3.4 Nexus .nexus, .nx

The Nexus format [START_REF] Maddison | Nexus: An extensible file format for systematic information[END_REF]) is very complex and can contain taxa and eventually trees. In this example we use an identifier for each taxon (a00000 to a00048) and we can establish a correspondence between the identifiers and the names of the taxa using the TaxLabels keyword. Dans cette thèse nous proposons différentes méthodes pour s'attaquer à la nature combinatoire de ce problème NP-complet. Premièrement, nous présentons un algorithme de Recuit Simulé compétitif qui nous a permis de trouver des solutions de meilleure qualité pour un ensemble de problèmes. Deuxièmement, nous proposons une nouvelle technique de Path-Relinking qui semble intéressante pour comparer des arbres mais pas pour trouver des solutions de meilleure qualité. Troisièmement, nous donnons le code d'une implantation sur GPU de la fonction objectif dont l'intérêt est de réduire le temps d'exécution de la recherche pour des instances dont la longueur des séquences est importante. Finalement, nous introduisons un prédicteur capable d'estimer le score optimum pour un vaste ensemble d'instances avec une très grande précision.

Abstract

Phylogenetic reconstruction is considered a central underpinning of diverse fields like ecology, molecular biology and physiology where genealogical relationships of species or gene sequences represented as trees can provide the most meaningful insights into biology. Maximum Parsimony (MP) is an important approach to solve the phylogenetic reconstruction based on an optimality criterion under which the tree that minimizes the total number of genetic transformations is preferred.

In this thesis we propose different methods to cope with the combinatorial nature of this NP-complete problem. First we present a competitive Simulated Annealing algorithm which helped us find trees of better parsimony score than the ones that were known for a set of instances. Second, we propose a Path-Relinking technique that appears to be suitable for tree comparison but not for finding trees of better quality. Third, we give a GPU implementation of the objective function of the problem that can reduce the runtime for instances that have an important number of residues per taxon. Finally, we introduce a predictor that is able to estimate the best parsimony score of a huge set of instances with a high accuracy.