Kathrin Mme

Gieseler

M Yannick Vanpoulle

M Bruno Andrioletti

Keywords: Cophylogeny, Enumeration algorithms, Equivalence relation, Dynamic programming. Résumé Cophylogénie, Algorithmes d'énumération, Relation d'équivalence, Programmation dynamique

To understand the dynamics of interaction between groups of species, for example, in a hosts/parasites (or hosts/symbionts) system, it is essential to consider the process of coevolution. The cophylogeny reconciliation model formulates the coevolution question as an optimization problem. The set of optimal solutions represents different scenarios of coevolution which need to be analyzed separately. The main result of this thesis is a novel approach that addresses the issue of the often huge number of optimal solutions which makes it difficult to analyze the results. We introduce several biologically motivated equivalence relations, each one splitting the solution space into equivalence classes. We propose polynomial-delay algorithms that enumerate the equivalence classes, and the representative solutions, i.e., the solutions taken from each equivalence class. Experimental results are then presented to illustrate the practical benefits of considering equivalence classes as a means of efficiently exploring the solution space of cophylogeny reconciliation. Based on the algorithmic results, a software called Capybara has been developed as a practical tool for analyzing cophylogeny datasets. More generally, in the area of enumeration algorithms, the problem of efficiently enumerating representative solutions is of theoretical importance. In the second part of this thesis, we provide a general framework for the enumeration of the equivalence classes of solutions for a family of optimization problems. We show that this framework can be applied to various dynamic programming problems, of which the cophylogeny problem is a special case.

Résumé en français

Le lecteur avisé ne devrait pas constater avec surprise que, contrairement à ce qu'un titre de « dynamiques des interactions inter-espèces » aurait suggéré, le sujet principal de cette thèse concerne non pas l'écologie mais l'évolution. En effet, des approches mêlant écologie et évolution ont été développées notamment par des écologues qui s'intéressent à l'impact des mécanismes coévolutifs sur des communautés biologiques en étroite interaction.

Pour préserver la biodiversité ou pour maîtriser l'émergence de maladies infectieuses, les sciences de l'écologie doivent répondre au défi de la prédiction. Puisqu'il est impossible de concevoir un avenir sans une compréhension approfondie du passé, une science prédictive, quelle que soit l'échelle temporelle visée, s'appuie nécessairement sur les sciences du passé : géologie, paléontologie, archéologie, histoire, ou encore, la biologie évolutive. La cophylogénie, la juxtaposition des deux arbres phylogénétiques, consiste à reconstruire l'histoire évolutive conjointe de deux groupes d'organismes ayant à l'heure actuelle une certaine relation écologique entre eux, par exemple un groupe de mammifères et un groupe de leurs parasites. Il s'agit alors d'établir deux types de liens de façon simultanée : des liens écologiques entre ces deux groupes d'espèces, ce qui délimite le périmètre de l'étude, et des liens évolutifs entre le présent et le passé, ce qui accorde aux résultats un pouvoir informatif et prédictif.

Il est facile de représenter des liens entre des organismes vivants : connaissant quels parasites vivent chez quels hôtes, nous pouvons simplement écrire leurs noms sur un papier et les relier par des traits de crayon. Sur l'axe temporel, comment représenter des liens évolutifs entre le présent et le passé ?

La philosophie de la Grèce antique s'offre deux temporalités différentes : Chronos et i Kairos. Chronos fait référence au temps qui défile de manière linéaire et séquentielle, tandis que Kairos symbolise le temps du moment opportun pour un geste décisif. Un arbre phylogénétique, avec ses noeuds à chaque ramification, incarne Kairos, les points de basculement où de nouvelles espèces sont nées. De même, la cophylogénie cherche à établir des liens évolutifs qui ne s'inscrivent pas dans une temporalité chronologique, mais qui peuvent être identifiés par des événements coévolutifs, une sorte de transition critique durant laquelle se forment de nouvelles associations écologiques.

Le problème cophylogénétique est le suivant. Nous disposons de deux arbres phylogénétiques (qui représentent l'histoire évolutive de deux groupes d'espèces, comme Un algorithme, une suite finie et non ambiguë d'instructions qui peuvent s'exécuter sur un ordinateur, apporte une réponse à la question cophylogénétique au moyen d'un autre type de lien plus abstrait. Entre les liens écologiques inter-espèces et les liens évolutifs composés d'événements coévolutifs, un algorithme forme le maillon central qui relie deux espaces : l'un est l'espace d'hypothèses biologiquement plausibles de l'histoire coévolutive, l'autre est l'espace de solutions optimales du problème informatique.

La vraisemblance d'une hypothèse coévolutive se rapporte à une pluralité de réflexions. Or, un modèle mathématique ne peut prendre en compte qu'un petit nombre d'entre elles. En outre, l'inclusion de certaines contraintes biologiques pourrait rendre le problème informatique difficile à résoudre, ou, comme les informaticiens l'appellent, NP-difficile. La réponse apportée par un algorithme ne peut donc pas être une réponse directe à la question biologique. Si l'espace des hypothèses biologiques n'est jamais identique à celui des solutions du problème informatique, le maillon algorithmique entre ces deux espaces repose sur une intime conviction des bioinformaticiens : avec un modèle soigneusement élaboré, la vérité biologique est intégralement saisie, et il ne reste plus qu'à la localiser au sein de l'espace de solutions, c'est-à-dire, dans la sortie de l'algorithme. En d'autres termes, si un algorithme pouvait fournir une représentation exhaustive de l'espace de solutions, une analyse appropriée de sa sortie permettrait de dévoiler l'ensemble des scénarios coévolutifs plausibles.

Ainsi, des algorithmes d'énumération ont été développés pour le problème de cophylogénie. Un tel algorithme produit une liste contenant différentes manières d'associer les deux arbres phylogénétiques, toutes optimales par rapport aux contraintes mathématiques basées sur les événements coévolutifs. Un biologiste peut alors analyser chacune de ces solutions, appliquer des critères de sélection supplémentaires, et arriver à un espace d'hypothèses coévolutives appuyé par son expertise.

La notion d'espace de solutions joue un rôle clé dans notre investigation algorithmique de la cophylogénie et marque le point d'entrée de cette thèse. Bien qu'un algorithme d'énumération efficace semble donner une réponse satisfaisante à la question biologique, il y a toutefois un obstacle de taille : il s'agit de la taille de l'espace de solutions, c'est-à-dire du nombre de solutions optimales qui le composent. Afin de soumettre la liste de solutions à un biologiste qui procédera ensuite à l'analyse manuelle des résultats, ce nombre-là doit rester raisonnable. Au cas où il y aurait plusieurs milliers de solutions, la tâche d'analyse de résultats pourrait être confiée à un programme informatique. Cependant, aucun programme informatique ne pourra jamais traiter une liste contenant un si grand nombre d'objets qu'il dépasse le nombre d'atomes dans l'univers ! L'espace de solutions pour le problème de cophylogénie peut, malheureusement, Dans le Chapitre , le plus grand espace de solutions que nous allons rencontrer est de l'ordre de .

Le nombre d'atomes dans l'univers a un ordre de grandeur estimé entre et (une brève explication de cette estimation se trouve au lien suivant : https://physics.stackexchange.com/questions/47941/ dumbed-down-explanation-how-scientists-know-the-number-of-atoms-in-the-universe).

devenir incroyablement vaste, et surtout si on s'intéresse à un système d'interactions hôtes-parasites impliquant des milliers d'organismes. Que faire dans ce cas-là ?

Une stratégie possible consiste à affiner le modèle en y ajoutant des contraintes supplémentaires. Davantage de contraintes pourront effectivement réduire la taille de l'espace de solutions. Néanmoins, comme nous l'avons évoqué plus tôt, il se peut que ces contraintes rendent le problème bien plus difficile à résoudre, jusqu'au point où l'on se retrouve dans l'incapacité de concevoir un algorithme efficace. De plus, de nouvelles contraintes exigent de nouvelles informations biologiques (e.g., des connaissances biogéographiques qui excluent l'association entre des espèces spatialement séparées), qui peuvent s'avérer difficiles à obtenir ou insuffisamment fiables.

Le travail de cette thèse relève d'une autre stratégie : naviguer à travers un vaste espace de solutions de manière efficace. Lorsqu'il est irréaliste de faire un parcours exhaustif de la longue liste de solutions, un « parcours guidé » de quelques solutions représentatives laisserait peut-être déjà se dessiner le contour des scénarios coévolutifs les plus vraisemblables. v

Dans le

Acknowledgements

Throughout the preparation and the writing of this thesis, I have received a great deal of support and assistance from many people. As someone who is quite uncomfortable with names (even with my own!), I would avoid using names as much as possible (the same applies to the rest of the manuscript).

I would like to express my deep gratitude to my research supervisors, for their patient guidance, enthusiastic encouragement and invaluable critiques of this research work. I would also like to thank the members of the defense jury, the thesis committee, and the administrative units from the lab and the university.

I would like to acknowledge all the colleagues with whom I worked in the same office, the same team, the same lab, the same building, the same campus. I should also appreciate my confrères outside of academia, collaborators, family members, partners, companions, and those whom I met in one of the many serendipitous encounters. This thesis would not have been done in the way that it has been done without the wonderful human environment in which I continue to evolve.

vii ì e (φ)

Table of contents

Résumé en français i

The event vector of a reconciliation.

T(G)

The set of solution subtrees of an ad-AND/OR graph.

gch

The tuple of grandchildren of a node.

π(T(G))

The set of equivalence classes of a reconciliation graph.

Introduction

All my apologies to the reader who, intrigued by the title of dynamics of species interactions and eager to catch a glimpse of the complex mechanisms within ecological communities, will discover with surprise that most of this thesis is concerned with evolution rather than ecology. Yet, it is a hallmark of the modern-day ecology to consider the process of coevolution for the purpose of investigating biodiversity and certain aspects of disease ecology.

To understand the present and to predict the future, humankind possesses an essential source of knowledge: the past. Just as geology, paleontology, archaeology or history, evolutionary biology makes it possible to establish a link between the present time and a bygone time, either distant or nearby. In cophylogeny, evolution is studied with regard to two taxonomic groups assumed to have some kind of current relationship with each other, such as where one is a group of mammal species and the other is a group of parasites of those mammals. Cophylogeny thus tries to elucidate two links simultaneously: the ecological link between the two groups of species, which sets the perimeter of the investigation, and the evolutionary link between the present and the past, which endows the results with informative and predictive power.

To make the link between two groups of species is easy: after observing which parasites live in which mammals, we can simply write down their names on a paper and draw lines between them with a pencil. But how do we connect the present and the past?

The ancient Greeks had two words for time: chronos and kairos. Chronos refers to sequential time, while kairos signifies a proper time for a decisive act. A phylogenetic tree, with its nodes connecting separate branches, embodies kairos, the tipping points where new species are born. Likewise, the evolutionary link made by cophylogeny is not chronological in nature, but is identified by coevolutionary events, the critical transitions during which new ecological connections are formed.

The cophylogeny problem is the following. Our information thus far has three components: the two phylogenetic trees (representing respectively the evolutionary history of a group of host species and of a group of parasite species), and a set of associations that state which taxon in one tree is associated with which taxon in the other tree. Our objective is to find a way of associating the ancestral species in the parasite evolutionary tree with places in the host evolutionary tree, in a manner that makes the most sense, and we judge what makes sense in terms of the coevolutionary events that would account for those past associations.

Formulated as a computational problem, the phylogenetic tree reconciliation model translates the cophylogeny problem into a hunt for efficient algorithmic tools. Our investigation is therefore an algorithmic one.

Algorithms, finite sequences of well-defined and computer-implementable instructions, help answering the cophylogeny question by providing another, more abstract, kind of link. Built on top of the ecological link between species and the evolutionary link composed of coevolutionary events, the algorithmic link connects spaces: one is the space of biologically plausible hypotheses on associations and events, which forms the answer of the biological question, the other one is the space of optimal solutions of the computational problem, which is a mathematically well-defined set, and can constitute the output of a computer program.

For a hypothesis about the coevolutionary history to be biologically plausible, a broad range of factors should be taken into account, although only a select few can be encoded as mathematical constraints into the computational problem. Besides, there are biological constraints that, once added to the computational problem, render it intractable, or, as the computer scientists call it, NP-hard. Needless to say, an algorithm

does not directly answer the biological question: even with the most accurate mathematical model, the space of biological hypotheses is never equivalent to the space of solutions of the computational problem. By saying that there is an algorithmic link be-tween the two spaces, we express our belief that, given a sufficiently accurate model, biological insights should be hidden somewhere in the solution space, and hence could be found in the output of the algorithm. In other words, an algorithm that supplies the whole solution space should open the door to understanding different possibilities of coevolutionary scenarios. This is the case for enumeration algorithms that have been developed for the cophylogeny problem. Such an algorithm enumerates, that is, provides a list of different ways of associating the two phylogenetic trees, all of which are optimal with respect to mathematical constraints based on the occurrences of coevolutionary events. Then, a biologist can analyze each one of these solutions, apply additional selection criteria, and come up with her own space of biologically plausible coevolutionary hypotheses.

Indeed, the solution space plays a pivotal part in our algorithmic investigation of the cophylogeny problem, and is the starting point of this thesis. It may seem that using an efficient enumeration algorithm we can already give a satisfactory answer to the biological question. There is nonetheless a stumbling block: the size of the solution space. For a human expert to be able to analyze the solutions one by one, the number of solutions cannot be too large. When there are thousands of solutions, we can resort to computer programs that automatically perform the analysis. However, no computer program can ever process a list containing such a huge number of items that it exceeds the number of atoms in the observable universe!

For the cophylogeny problem, the size of the solution space does unfortunately get incredibly large, especially if we want to study a hosts/parasites system composed of thousands of organisms. What can we do in such cases?

One possible strategy is to refine the model by including more constraints. This new version of the computational problem will necessarily have a smaller solution space.

We already mentioned one major drawback of this approach: adding more constraints can make the problem harder to solve, and we may be unable to find efficient algo-

In Chapter , the largest solution space that we will encounter has a size in the order of . It is estimated that the there are between to atoms in the observable universe (a short explanation of this estimation can be found in https://physics.stackexchange.com/questions/47941/ dumbed-down-explanation-how-scientists-know-the-number-of-atoms-in-the-universe).

rithms. Those constraints also require further biological information (such as biogeographical information that makes associations between spatially separated species impossible), which can be difficult to obtain or insufficiently reliable.

Another strategy consists in exploring the solution space efficiently despite its large size. Many efforts have been undertaken in this direction. The work of the current thesis also lies within the scope of this approach.

In Chapter , along with a succinct presentation of the mathematical and biological background, we introduce our main algorithmic challenge: finding representative solutions of the cophylogeny problem. Chapter sets out the formal definition of the model and our theoretical results. First, we define equivalence relations on the set of solutions and explain the biological motivation of the definitions. Mathematically, an equivalence relation splits the solution space into parts of mutually equivalent subsets, called equivalence classes. In this setting, we find representative solutions through the search of solutions that are not equivalent to each other, i.e., solutions in different equivalence classes. We propose algorithms for enumerating the equivalence classes, for each notion of equivalence relation. In Chapter , we show the experimental results that illustrate the practical benefits of using the equivalence classes as a means of exploring the solution space. A software called Capybara is also presented in Chapter as a practical tool for analyzing cophylogeny datasets. Then, we conclude our work on the cophylogeny problem and present some perspectives.

Our algorithmic investigation of the cophylogeny problem has fostered two related investigations. In Chapter , we venture into the territory of the abstract: in a general framework where the solution space of some problem needs to explored, we propose an algorithm for enumerating the equivalence classes of solutions, with respect to an equivalence relation represented as a graph coloring. We show that the algorithm can be applied to various computational problems, of which the cophylogeny problem is a special case. In Chapter , we turn our attention to another problem in phylogenetic studies that has to do with the quantification of dissimilarity between phylogenetic trees. A fast and accurate computation of a distance measure between trees has many applications. For instance, in a hosts/parasites system where species have co-evolved together, the dissimilarity between the phylogenetic trees of the two groups of species should be low. We propose a new algorithm for approximating the subtree prune-and-regraft distance between two phylogenetic trees, and we perform a series of experiments to demonstrate its practical advantages in comparison with existing algorithms.

The main algorithmic results and some of the experimental results from Chap- As little is known about the relationships between extinct species, phylogenetic trees are reconstructed by using phylogenetic inference methods, where the past relationships are inferred through mathematical models from features of contemporary species. Early methods for reconstructing phylogenetic trees relied on morphological or structural characteristics. Today, advances in molecular biology, especially in sequencing technologies, have allowed the extensive use of DNA or protein sequence data to build phylogenetic trees.

Phylogenetic inference methods that use molecular data are based on the identification of homologous genes, that is, genes from different organisms that have a shared ancestry [Fit ]. The development of new sequencing technologies has given rise to massive datasets of ever increasing size, both in the number of genes and in the number of taxa (i.e., species) that can be considered [KYT ]. Consequently, more sophisticated methods have been developed to extract the most information out of those datasets, resulting in more accurate trees on larger and larger scales [Hug+ ; Bur+ ; CK ]. Figure . shows, for example, the tree of life constructed by Hug et al. [Hug+ ] using genomic data from organisms.

As the volume of the data is huge and complex methods are required to deal with the uncertainties of the data and systematic errors (these errors are generally due to heterogeneous rates of evolution across taxa or time, which violate the assumptions of the models), modern phylogenetic inference methods pose a heavy computational burden. The inference of the tree in Figure . required computational hours on a supercomputer.

Algorithmic efficiency is thus crucial. In Chapter , from an algorithmic standpoint, we will look at one problem that naturally arises in phylogenetics analysis, that is, the quantification of the difference between phylogenetic trees.

. . Cophylogeny

In cophylogeny, one is interested in reconstructing the coevolutionary history of groups of species that interact over a long period of time and in ecologically close environments based on their phylogenetic information. This prolonged and intimate interaction between different species is known as symbiosis. Symbioses are a common feature of life and can be found across all major phylogenetic lineages [Pou ]. Typically defined between a larger host organism and a smaller symbiont, symbiotic relationships may be mutualistic, where both partners benefit from the interaction, parasitic, where one partner benefits and the other suffers a cost, or commensalistic, where one partner benefits while the other receives neither benefit nor harm [LP ; EKM ]. Sto+ ]. In a nutshell, doing a phylogenetic tree reconciliation means to map the parasite tree to the host tree. Once an appropriate mapping between the two trees is chosen, one should be able to identify coevolutionary events. In general, four main evolutionary events are considered: cospeciation (when both parasite and host speciate), duplication (when the parasite speciates but not the host), loss (when the host speciates but not the parasite), and host-switch (when the parasite speciates and one of the new species goes to infect another unrelated host). A schematic representation of each of the four types of events is given in Figure . . Formulated as a computational optimization problem, the phylogenetic tree reconciliation model includes a cost assigned to each type of events. The goal is then to find a most parsimonious reconciliation, i.e., a mapping between the two trees that minimizes the total cost of the events. cospeciation duplication host-switch loss If timing information (i.e., the order in which the speciation events occurred in the host phylogeny) is not available (or sufficiently reliable to be used with enough confidence), as is usually the case, finding a most parsimonious reconciliation is a NP-hard problem [Ova+ ; THL ]. To deal with this, one can either rely on heuristics [Con+ ],

or accept solutions that may be biologically unfeasible, that is, to accept reconciliations in which some of the host-switches induce a contradictory time ordering on the host species [BAK ; Don+ ]. Following this second path, efficient algorithms have been proposed to enumerate optimal reconciliations. The output of those algorithm consists of not one but a list of reconciliations that all achieve the minimal cost within the model.

One issue that arises then is the size of the output. As the output of the enumeration algorithm can be composed of billions of solutions, there is no hope of analyzing each one of them separately. For instance, for the dataset of Wolbachia and their arthropod hosts, collected by the Biometry and Evolutionary Biology Laboratory [Sim+ ; Sim ], the number of optimal reconciliations is unrealistically large (e.g., in the order of , see Chapter ).

The good news is that, even when the number of solutions is huge, the set of so- In Chapter , we propose a new exact method that allows to efficiently enumerate representative solutions based on the biologically motivated notion of equivalence:

lutions
the output of the algorithm is then a list of reconciliations that are biologically nonequivalent to each other. In Chapter , we will see that in practice, by only considering the representative solutions, we effectively reduce the number of solutions that need to be analyzed while still maintaining important biological information about the solution space (for example, Is there a high incidence of cospeciations? In which branches of the parasite tree have there probably been host-switches?).

.

Mathematical background

. . Graphs, trees, forests

Graphs

Throughout the manuscript, we will use several basic notions in graph theory such as cyclicity and connectivity. The definitions can be found in any classical textbook such as [Die ]. For a directed graph G, we denote by V(G) and A(G) respectively the set of nodes and the set of arcs of G. An arc is denoted by an ordered pair of nodes (v, w).

The set of out-neighbors of a node v is called its set of children and is denoted ch(v). If w ∈ ch(v), the node w is a child of v. A directed acyclic graph (DAG) is a directed graph without any directed cycle. A directed tree is a DAG whose underlying undirected graph is acyclic and connected.

Trees

In this manuscript, a tree will refer to a rooted directed tree in which arcs are directed away from the root. Any such tree T has the following property: ( ) exactly one node, called the root of T and denoted by r(T ), does not have any in-neighbor; ( ) any other node v has exactly one in-neighbor in T ; this unique in-neighbor is called the parent of v and is denoted by p(v). A node without any child is called a leaf. We denote by L(T )

the set of leaf nodes of T . The non-leaf nodes, i.e., nodes which are not leaf nodes, are called the internal nodes of T .

In a tree T , if there exists a directed path from a node v to a node w, the node w is called a descendant of v, and v is called an ancestor of w; if moreover v ≠ w, we say that w is a proper descendant of v, and that v is a proper ancestor of w. If neither w is an ancestor of v nor w is an ancestor of v, we say that the two nodes are incomparable, and denote it v w w. When v and w are comparable (i.e., they are not incomparable), the number of arcs on the directed path between v and w is called the distance between v and w, and is denoted by d T (v, w) (it is a non-negative integer). The largest distance between a node v and a leaf of T is called the height of the node v. The leaf nodes have height zero. We denote by LCA (v, w) the lowest common ancestor of two nodes v and w, that is, the node having the smallest height among all nodes which are ancestors of both v and w. The height of a tree T is defined to be the height of its root node. A tree having one single node (its root is also a leaf) has height zero.

When T is a tree and v is a node of T , the subtree of T rooted at v is the subtree of T containing all descendants of v, and is denoted by T | v . Clearly, T | v is a rooted tree, and

its root is r(T | v ) = v.
This way of taking subtrees does not change the height of the nodes: for any node w of T | v , its height in T | v is the same as its height in the original tree T .

Forests

In this manuscript, a forest will always refer to a DAG whose connected components are rooted directed trees. A component having one single node is called an isolated node.

The order of a forest F, denoted by |F|, is its number of connected components. A tree is a forest of order one. All terms that are previously defined for trees can be naturally generalized to forests. In particular, a node v of a forest F is a leaf (resp. a root) if F has a component T such that v is a leaf (resp. the root) of T . If two nodes v and w are in two distinct components of F, the lowest common ancestor LCA (v, w) does not exist.

. . Phylogenetic trees and forests

A full rooted binary tree is a rooted tree in which every internal node has exactly two children. Let T be a full rooted binary tree having n := |V(T )| nodes. The following results are easy: ( ) T has n+ leaves and n-internal nodes; ( ) the maximum distance between any two comparable nodes in T (this is also called the diameter of T ) is bounded by n-= O(n).

An ordered full rooted binary tree is a full rooted tree that satisfies the following: for every internal node v, the set of children ch(v) has a fixed ordering; we will write ch(v) Let X be a fixed arbitrary finite set, called the set of taxa. An X-tree is a tree whose leaves are bijectively labeled by the set X. The element of X assigned to a node v by the bijection is called the label of v.

as
Definition . (Phylogenetic trees and forests). Let X be a fixed set of taxa. A phylogenetic tree on X is a tree that satisfies all of the following properties: ( ) it is a full rooted binary tree; ( ) it is ordered; ( ) it is an X-tree. A phylogenetic forest is a forest whose connected components are phylogenetic trees.

Biologically, two different phylogenetic trees according to Definition . can reflect the same evolutionary relationships between the taxa, because the ordering of the children of a node does not matter. We will use the notion of isomorphic trees to qualify phylogenetic trees that are biologically the same.

Definition . (Isomorphic phylogenetic forests). Two phylogenetic forests F and F on the same set X of taxa are isomorphic, denoted by F F ,if there exists a bijection between their node sets ϕ : V(F ) → V(F ) that satisfies all of the following conditions:

. For all pairs of nodes

(u, v) ∈ V (F ), (u, v) ∈ A(F ) ⇐⇒ (ϕ(u), ϕ(v)) ∈ A(F ). . For each node u ∈ V(F ), u is a root of F if and only if ϕ(u) is a root of F .
. For each leaf u ∈ L(F ), ϕ(u) is a leaf of F , and the labels of u and ϕ(u) are identical.

The first condition ensures that the two forests are isomorphic as directed graphs;

the second condition additionally ensures that they are isomorphic as rooted forests, in other words, there exists a bijection between the components of two forests that identifies isomorphic rooted trees. Together, the first two conditions guarantee that the forests are the same as unordered forests, in other words, they have the same branching pattern, the same topology. The third condition ensures that the isomorphism pre-serves the labels of the leaves. For two phylogenetic trees (i.e., forests of order one) T and T to be isomorphic, the second condition can be simplified to: ϕ(r(T )) = r(T ).

Remarks

For the sake of theoretical clarity, we have chosen a definition of the phylogenetic trees that is slightly different from the one that is commonly used in the literature; the latter does not consider the trees to be ordered. Our choice is motivated by the fact that, in many practical situations, it is more convenient or natural to consider ordered trees. In any such case, one must always bear in mind that two trees that "look different" may be biologically the same:

• Two phylogenetic trees that are biologically the same (isomorphic) can have different representations in Newick format, the standard format for representing trees in phylogenetic analysis (see Figure . for examples).

• For the purpose of visualizing phylogenetic trees on a two-dimensional graphic (we are doing a planar embedding), an ordering of the children of each node must be chosen. Visualization software such as [Yu+ ] allow the user to flip the left and right subtrees manually. Two isomormphic trees can have different graphical representations (see Figure . ).

• When performing algorithms on phylogenetic trees, one often wants to traverse all the nodes in the tree. Assuming a sequential model of computation (i.e., not parallel), a convenient way of traversing an unordered tree would be to first choose a fixed ordering for the children of each node (so it becomes ordered), and then apply the standard traversal techniques for ordered binary trees (e.g., do a preorder depth-first traversal: first the root, then the left subtree, finally the right subtree).

. . Concepts in algorithms and complexity

Dynamic programming

Dynamic programming (DP) is an easy and powerful method that is widely used for solving combinatorial optimization problems. First formally developed by Richard Bellman in the 's [Bel ; Bel ], it has found numerous applications in various fields.

In bioinformatics, DP is used for solving problems such as Sequence alignment [NW ],

Protein structure comparison [HS ], RNA structure prediction [RE ], etc.

In Chapter , we will present a formal model for characterizing problems that can be solved using DP or DP-like techniques. Here, we will only provide the necessary intuitions and present some terminologies that will be used in Chapter .

Simply put, dynamic programming is a general problem-solving method that tries to solve a problem by breaking it down into simpler subproblems in a recursive manner (DP-subproblem). To solve an optimization problem using a DP-algorithm, the value to be optimized should be computed via some recurrence relations using the optimal values of the subproblems. Usually, these optimal values of the subproblems are stored in a structure called the dynamic programming ?), (c) define the output (how to obtain the optimal value of the optimization problem from the DP-table?).

As an example, consider the following optimization problem, known as the -K : given a knapsack of capacity W max and n items, each of value v i and of weight w i , the goal is to select a subset of items to put in the knapsack such that the total weight does not exceed W max and the total value is maximized. In the case where the weights w , . . . , w n and W max are positive integers, there is a well-known DP-algorithm, which we will now describe by identifying the three components. (a) A subproblem, indexed by (i, j), asks to find the maximum value M(i, j) that can be attained using the first i items and with total weight at most j, where ≤ i ≤ n and

≤ j ≤ W max . (b)
The recurrence relations are:

M( , j) = , M(i, j) =            M(i -, j) if j < w i , max{M(i -, j), M(i -, j -w i ) + v i } otherwise. (c) The output is M(n, W max ).
A problem may be solved by several different DP-algorithms. Depending on the formulation, the time and space complexity of a DP-algorithm may or may not be a polynomial of the input size (the above algorithm for the -K has pseudopolynomial complexity).

Counting, enumeration

An optimization problem asks to find one solution, among all the feasible solutions, such that its value (or measure) maximizes or minimizes the objective function. For each optimization problem, there is a corresponding decision problem that asks whether there is a feasible solution for some particular value (a decision problem is simply a problem whose output set consists of only two values yes or no). This is called the decision variant, or the associated decision problem of the optimization problem. We will most often consider a NP optimization problem: every feasible solution has polynomial size, and its value can be computed in polynomial time. In this case, the associated decision problem will also be in the complexity class NP.

We can also define the counting and the enumeration problems that are associated with a given optimization problem. The first is concerned with establishing how many solutions are optimal. The second is concerned with enumerating, i.e., producing a list without duplicates of all optimal solutions. Counting and enumeration problems without an underlying optimization problem can be formulated in a similar manner.

If we have an algorithm for the enumeration problem, it can obviously be turned into an algorithm for the counting problem. However, in most cases, we do not have an upper bound for the number of optimal solutions that is polynomial in the size of the input. If we are interested in finding this number in polynomial time, the counting algorithm should not rely on enumeration. Just as for a decision problem in the class NP, we can either try to find a polynomial time counting algorithm, or try to show that it belongs to a subclass of "harder" problems (NP-complete for decision problems);

the complexity class that characterizes the "harder" counting problems is the class #P-complete.

For the enumeration problem, since the output size can be exponential in the input size, the complexity is expressed in an output-sensitive manner. The performance of an enumeration algorithm can be measured in different ways:

• We can measure the running time in terms of the total time that is required for producing all solutions.

• By looking into the intermediate steps of the algorithm, we can measure the running time more carefully in terms of ( ) the preprocessing time, i.e., the time required to produce the first solution, and ( ) the maximum time delay between two consecutive outputs.

An algorithm is said to be total polynomial time or output polynomial time if its total running time is a polynomial of the input size and the output size. A total polynomial time algorithm is said to be incremental polynomial time if the time delay between the k-th and the (k+ )-th outputs is a polynomial of k and the input size. A total polynomial time algorithm is said to be polynomial delay if the time delay between two consecutive outputs is a polynomial of the input size. In Chapters and , we will see several polynomial delay enumeration algorithms.

Approximation, heuristic

An approximation algorithm is an efficient algorithm (polynomial time) that finds an approximate solution to an optimization problem, with provable guarantees on the quality of the approximate solution: this is usually a multiplicative factor called the approximation ratio. For a minimization problem, a k-approximation algorithm, that is, an approximation algorithm with ratio k, guarantees to output a feasible solution whose value is within k times the optimal (minimum) value. Searching for approximation algorithms is a widely used approach for dealing with NP-hard optimization problems.

Some algorithms are known as heuristics as they provide an approximate solution efficiently but do not provide a guaranteed approximation ratio. Another type of algorithm, also called heuristic, provides exact solutions with no guaranteed running time bound. We will encounter both kinds of heuristics in this manuscript.

Chapter

Phylogenetic tree reconciliation: theoretical contributions

.

Introduction and definitions

. . The R

In this section, we will define formally the P (shortly, the R ). At the end of the section, we will discuss the biological relevance of this definition and provide some references.

Reconciliation, events, event vector

Let H and P be two phylogenetic trees (see Definition . ), respectively on the set of taxa of host and parasite species. Let σ be a function from the parasite leaves L(P) to the host leaves L(H), representing the parasite/host associations between the taxa (i.e., the extant species). The triple (H, P, σ) will be called a dataset.

A reconciliation is simply a function φ that assigns, for each parasite node p ∈ V(P), a host node φ(p) ∈ V(H). It should represent the parasite/host associations for the extant species as well as for all the ancestral parasite species in the tree. Naturally, the function φ must extend the leaf association function σ, that is, for each p ∈ L(P), we should have φ(p) = σ(p).

Now, for an internal node p, what are the constraints on the associated host φ(p)?

In the event-based model, a reconciliation must induce events. Concretely, for each internal node p under a reconciliation φ, the events at p will have two components:

( ) exactly one event type (denoted by a special symbol) among ℂ, , , meaning respectively cospeciation, duplication, and host-switch, and ( ) zero or more loss events (denoted by a nonnegative integer). The event type and the number of losses at a node p will depend on the relationship between three host nodes in the host tree: the host φ(p) assigned to p by φ, and the hosts φ(p ) and φ(p ) which are assigned to the two children (p , p ) of p. Those dependencies can be expressed as an event type function E and a loss number function ξ, both taking as input three host nodes. The definitions of E and ξ are given in Table . .

Case

Value of Using the functions E and ξ, we now give in Definition . the formal definition of a reconciliation. In Definition . , we also define the event type E φ (p) and the loss number ξ φ (p) of a parasite node p under a reconciliation φ, as these notions will be ubiquitous in the subsequent theoretical developments.

E(h, h , h ) Value of ξ(h, h , h ) Both h , h
Definition . (Reconciliation). Given two phylogenetic trees H and P, and a function σ : L(P) → L(H), a function φ : V(P) → V(H) is called a reconciliation of (H, P, σ) if it satisfies the following two conditions:

. For every leaf node p ∈ L(P), φ(p) = σ(p).

. For every internal node p ∈ V(P) \ L(P) with children (p , p ):

E(φ(p), φ(p ), φ(p )) ∈ {ℂ, , }.
Definition . (Event type and loss number of a parasite node under a reconciliation).

Let φ be a reconciliation of (H, P, σ). For a parasite node p ∈ V(P), the event type E φ (p)

and the loss contribution ξ φ (p) of p under φ are defined respectively as follows:

• If p is a leaf (here is a special symbol for denoting the terminal event, i.e., the event type for the leaves):

E φ (p) := , ξ φ (p) := .
• Otherwise, p is an internal node with children (p , p ):

E φ (p) := E(φ(p), φ(p ), φ(p )), ξ φ (p) := ξ(φ(p), φ(p ), φ(p )) .
The image of the node-wise event type function E φ is denoted by E := { , ℂ, , }.

The loss event is denoted by another special symbol . The values of E φ partitions the set of internal parasite nodes into three disjoint subsets according to their event type; these subsets are denoted by V ℂ (P), V (P), V (P). The number of occurrences of each of the three event types together with the number of losses make up the event vector of the reconciliation φ, as defined in Definition . . We also define in Definition . the starred event vector which is simply the event vector without the number of losses.

Definition . (Event vector). The event vector of a reconciliation φ is the quadruple of integers consisting of the total number of each type of events ℂ, , , , i.e. ì e (φ) := V ℂ (P) , V (P) , V (P) , p ∈V(P)

ξ φ (p) . ( . )
Definition . (Starred event vector). The starred event vector of a reconciliation φ is the triple of integers that correspond to the first three elements of the event vector, i.e., ì e * (φ) := V ℂ (P) , V (P) , V (P) .

( . )

The following result is trivial: the sum of the three elements of ì e * (φ) is equal to the number of internal parasite nodes, that is, V(P)|-.

In Figure . , for ≤ i ≤ , we show six different reconciliations φ i , on the same dataset (H, P, σ). On the right side of each reconciliation φ i , we also marked the event type E φ i (p) and the loss number ξ φ i (p) for each internal p, and the event vector ì e (φ i ).

E φ (p ) = ξ φ (p ) = E φ (p ) = ℂ ξ φ (p ) = ì e (φ ) = ( , , , ) p p p c p b h h h b h c p a h a φ E φ (p ) = ξ φ (p ) = E φ (p ) = ℂ ξ φ (p ) = ì e (φ ) = ( , , , ) p p p c p b h h h b h c p a h a φ E φ (p ) = ℂ ξ φ (p ) = E φ (p ) = ξ φ (p ) = ì e (φ ) = ( , , , ) p p p c p b h h h b h c p a h a φ E φ (p ) = ℂ ξ φ (p ) = E φ (p ) = ξ φ (p ) = ì e (φ ) = ( , , , ) p p p c p b h h h b h c p a h a φ E φ (p ) = ℂ ξ φ (p ) = E φ (p ) = ℂ ξ φ (p ) = ì e (φ ) = ( , , , ) p p p c p b h h h b h c p a h a φ E φ (p ) = ξ φ (p ) = E φ (p ) = ξ φ (p ) = ì e (φ ) = ( , , , ) p p p c p b h h h b h c p a h a φ

Cost of a reconciliation, R

A cost vector ì c := (c(ℂ), c( ), c( ), c( )) is a quadruple of real numbers representing the costs that we want to assign to each event type (including the loss event). Given a cost vector ì c, the cost of a reconciliation φ is defined as the dot product between its event vector and the cost vector: ì e (φ) • ì c.

We are now ready to define the combinatorial optimization problem which we call the R . Recall that a dataset is a triple (H, P, σ). An instance is a pair consisting of a dataset and a cost vector ((H, P, σ), ì c). Given an instance, the set of feasible solutions is the set of reconciliations of the dataset (H, P, σ). The goal of the optimization problem is to find an optimal feasible solution, that is, a reconciliation φ that minimizes the cost ì e (φ) • ì c. In the associated enumeration problem, the goal is to find all reconciliations of minimum cost.

For example, for the dataset depicted in Figure . , if the cost vector is ( , , , ), all reconciliations are optimal; if the cost vector is ( , , , ), only φ is optimal.

Remarks

The R is studied extensively in the literature [THL ; BAK ; Don+ ; Ma+ ]. Notice that, for the sake of clarity, we used a slightly different set of notations to define the event-based model (in particular, we rely on the functions E and ξ in Table . ). The computational problem is equivalent to the ones that can be found in those references.

There exists in the literature a non-equivalent formulation of the R which considers the host tree to be dated: the input also includes a time function that maps each host node to a positive number [Doy+ ; Sco+ ]. This different computational problem is not studied in the present manuscript.

For any given dataset, the set of feasible solutions is nonempty: a function φ that satisfies the leaf association constraint (condition . of Definition . ) and that maps every internal parasite node to the root of the host tree (for example, φ of Figure . ) is always a valid reconciliation.

. .

Time-feasibility

A phylogenetic tree T , which represents a possible scenario for the evolutionary history of species, implies the existence of a node ordering that is consistent with the ancestor/descendant relations in the tree. In practice, we can consider a time function τ : V(T ) → ℝ > that maps each node to a positive number and satisfies the following:

for each arc (v, w) ∈ A(T ), we have τ(v) < τ(w). By transitivity, if v is an ancestor of w, then τ(v) < τ(w).

For a dataset (H, P, σ) of the R , let τ H be the time function for the host tree. For a reconciliation φ to be biologically meaningful, the following property must be satisfied:

for each arc (p, p i ) ∈ A(P), we should have τ H (φ(p)) < τ H (φ(p i )).

( . )

From the definition of a reconciliation ( ter, we will concentrate on the R , while keeping in mind that the solutions that we will explore may or may not be time-feasible. In Chapter , Section . , we will see that, for many instances, by simply enumerating the optimal reconciliations and filtering out those that are not time-feasible, we can already obtain time-feasible reconciliations. We will also present in Chapter some more methods for obtaining time-feasible reconciliations.

. . The reconciliation graph

Putting aside the time-feasibility and returning to our R , (or the reason for its popularity) comes primarily from the fact that, even though the number of optimal solutions can grow exponentially with the size of the input trees, many features of the solution space can be effectively understood by exploring this graph structure of polynomial size [Ma+ ]. In this section, we will describe the DPalgorithm for constructing the graph that we call the reconciliation graph. This graph will serve as the basis for a number of problems that we will tackle later in the chapter.

Dynamic programming

We will first describe a DP-algorithm that finds the optimal value, i.e., the minimum cost of any reconciliation. Using the terminologies from Section . . , given an instance ((H, P, σ), ì c), we need to: (a) define the DP-subproblems and the DP-table, (b) write down the recurrence relations for filling up the DP-table, (c) define the output.

Recall that E := { , ℂ, , } is the set of possible event types for a node. Given a parasite node p ∈ V(P), we denote by P| p the subtree of P rooted at node p, and by σ| L(P| p ) the restriction of σ to the leaves of P| p . For the purpose of describing the DP-subproblems, for a fixed dataset (H, P, σ), we will take the subtree of P at a particular node p ∈ V(P).

In The output, i.e., the minimum cost of any reconciliation of (H, P, σ) is given by min h∈V(H),e∈E M(r(P), h, e).

• If p is a leaf, M(p, h, e) =            if h = σ(p)
( . )

ad-AND/OR graphs

Before describing how to construct the reconciliation graph, we will first define a more general graph structure: the acyclic decomposable AND-OR graph (ad-AND/OR graph), which is known for having an intimate relationship with dynamic programming on a tree. The reconciliation graph will be a particular ad-AND/OR graph with some additional properties. In Chapter , we will study the relationship between general DPalgorithms and ad-AND/OR graphs and provide some references.

Definition . (ad-AND/OR graph).

A directed graph G is an acyclic decomposable AND/OR graph (shortly, ad-AND/OR graph) if it satisfies the following:

• G is a acyclic (it is a DAG).

• G is bipartite: its node set V(G) can be partitioned into (A, O) so that all arcs of G are between these two sets. Nodes in A are called AND nodes; nodes in O are called OR + nodes.

• Every AND node has in-degree at least one and out-degree at least one. The set of nodes of out-degree zero is then a subset of O and is called the set of goal nodes;

the remaining OR + nodes are simply the OR nodes. The subset of OR nodes of indegree zero is the set of start nodes.

• G is decomposable: for any AND node, the sets of nodes that are reachable from each one of its child nodes are pairwise disjoint.

An ad-AND/OR graph represents the space of solutions for some problem instance.

A solution will be a particular kind of subgraph which is a rooted tree.

Definition . (Solution subtree).

A solution subtree T of an ad-AND/OR graph G is a subgraph of G which: ( ) contains exactly one start node; ( ) for any OR node in T it contains one of its child nodes in G, and for any AND node in T it contains all its children in G.

The set of solution subtrees of G is denoted by T(G). It is immediate to see that a solution subtree is indeed a subtree of G: it is a rooted tree, the root of which is a start node. If we would drop the requirement of G being decomposable, the object defined in Definition . would not be guaranteed to be a tree.

It is easy to see that the following four-step procedure allows to obtain one solution subtree: ( ) start at any start node, ( ) for any visited OR node, visit one child, ( ) for any visited AND node, visit all children, ( ) stop when the goal nodes are reached. In fact, a well-known folklore result states that the set T(G) of solution subtrees of any ad-AND/OR graph G can be enumerated efficiently: after a pre-processing step in time linear in the size of G, the delay between outputting two consecutive solutions is linear in the size of the solution. In Chapter we will discuss this result again.

We will define now a notion of subgraph for ad-AND/OR graphs that we will use extensively later in this chapter as well as in Chapter . Intuitively, the relationship between such a subgraph and the entire graph correspond to the relationship between a DP-subproblem and the "full" optimization problem. ii. Create a single OR node, connect it to every AND node created in the previous step. Its label is (p, h, e), and its value is the minimum of the values of its children.

. For each (r(P), h, e) ∈ U, remove the OR node labeled by that cell unless its value is equal to the optimal cost (see Equation ( . )). For each OR node w, remove the arc to its child AND node w i if the value of w i is not equal to the value of w. Finally, remove recursively all AND nodes without incoming arcs.

It can be checked that the reconciliation graph is indeed an an-AND/OR graph as defined in Definition . . An OR + node labeled by (p, h, e) is a start node if and only if p = r(P), and is a goal node if and only if p ∈ L(P). It is also immediate to see that each AND node in the reconciliation graph has exactly one in-neighbor and exactly two children.

We will consider the two children of each AND node to be ordered. More precisely, for an AND node v, if its in-neighbor is labeled by (p, h, e) where p has two ordered children (p , p ), we know that the two children of v are respectively labeled by cells of the form (p , h , e ) and (p , h , e ): we will call the first one the left child w and the second one the right child w of v, and we will write ch(v) = (w , w ). However, the set of children of an OR node w is unorderd, and we will write ch(w) = {v i }.

For an OR node, we will typically be interested not in its set of children but in its set of "grandchildren", hence, we introduce here a new notation. If w is an OR node, we call the grandchild couples of w, denoted by gch(w), the union of the children of its child AND nodes (it is a set of couples of OR + nodes): gch(w) The space required for storing the reconciliation graph is thus O(mn ). Practical algorithms do not construct the reconciliation graph in the same way as we just described, as the functions E and ξ are not trivial to compute. Nevertheless, the reconciliation graph can be constructed using O(mn ) time and space complexity [Don+ ].

:= v i ∈ch(w) ch(v i ).

Solution subtrees of the reconciliation graph are optimal reconciliations

As announced in the beginning of this section, the reconciliation graph is a compact representation of the set of optimal reconciliations. Formally, if G is the reconcilia-tion graph for a given instance of the R , then the set T(G)

of solution subtrees of G corresponds bijectively to the set of optimal solutions of that instance, i.e., the set of reconciliations of minimum cost. Moreover, we have the following result for the subgraph G/{w} starting from an OR + node w labeled by (p, h, e):

the set T(G/{w}) of solution subtrees of the subgraph corresponds bijectively to the set of optimal solution of the DP subproblem indexed by (p, h, e).

In practice, to convert a solution subtree T ∈ T(G) into a reconciliation φ, we only need to look at the labels (p, h, e) of the OR + nodes in T (notice that a reconciliation can simply be viewed as a collection of triples of the form (p, h, e)). We will henceforth use interchangeably the terms solution subtrees of the reconciliation graph and optimal reconciliations of the problem instance.

.

Defining equivalence relations

In this section, we will start by giving the formal definitions for a number of different equivalence relations on the set of reconciliations of a given dataset, then we will explain the motivation for such definitions, show some preliminary results, and state the computational problems that are related to those equivalence relations. Each one of these equivalence relation splits the set of optimal reconciliations of a given instance into equivalence classes, i.e., subsets of pairwise equivalent reconciliations. One representative of an equivalence class is simply a reconciliation in the corresponding subset. We will abuse the terminology and call equivalence classes the objects that best represent the common property of the reconciliations in that subset (e.g., a V-equivalence class is simply a vector). A reconciliation in an equivalence class will then be a reconciliation satisfying that property.

p , h b , p , h , ℂ p , h a , p , h c , p , h , p , h b , p c , h c , p a , h a , p , h , ℂ p b , h b , p , h c ,
We will now introduce the notations for denoting the equivalence classes and making statements such as "a reconciliation is in an equivalence class". Recall that, if G is the reconciliation for a given instance, then T(G) denotes the set of optimal reconciliations. For a fixed equivalence relation, we consider a function π that maps each recon-

ciliation φ ∈ T(G) to the equivalence class π(φ) of φ (in other words, φ is in the equiva- lence class π(φ)).
For the E-, EL-, and CD-equivalence relations, the "common property" of a class will be expressed via some local properties on each parasite node; we define in Definition . this localized property that we call the color of a node. Then, the definitions of π for each equivalence relation are given in Definition . .

Definition .

(Color of a parasite node under a reconciliation). Let φ be a reconciliation of (H, P, σ). For each parasite node p ∈ V(P), the color κ φ (p) of p under φ is defined in Table . .

Equivalence relation

The color κ φ (p) of p under φ (it is either a couple or a triple)

E-equivalence relation (p, E φ (p)) EL-equivalence relation                    (p, E φ (p)) if E φ (p) ≠ (p, L ) if (p, p ) is the host-switch arc at p (p, R ) otherwise, (p, p ) is the host-switch arc
where L and R are two special symbols CD-equivalence relation

           (p, φ(p), E φ (p)) if E φ (p) ≠ p, ?, otherwise, E φ (p) =
where ? is a special symbol (denoting unknown host)

Table . : For each of E-, EL-, and CD-equivalence relations, definition of the color of a parasite node under a reconciliation.

Definition . (Equivalence class)

. Consider a fixed instance on a dataset (H, P, σ) and denote by T(G) the set of optimal reconciliations. The equivalence class π(φ) of a reconciliation φ is the object defined in Table . .

Equivalence relation

The equivalence class π(φ) of φ

V-equivalence relation ì e (φ) (it is a vector of four integers) V * -equivalence relation ì e * (φ) (it is a vector of three integers) E-, EL-, CD-equivalence relations π(φ) is a function that maps each p ∈ V(P) to its color κ φ (p)
Table . : For each equivalence relation, definition of the equivalence class of a reconciliation.

For a fixed equivalence relation, the set of equivalence classes of the instance, denote by π(T(G)), is defined as follows:

π(T(G)) := φ∈T(G) π(φ). ( . )
. .

Motivations

The first and foremost motivation of defining equivalence relations is the need of capturing useful biological information from the set of optimal reconciliations, when this set is too large for manual analyses or for exhaustive enumeration. While the number of optimal reconciliations can grow exponentially fast with the size of the input trees, we will see in Section . . that the number of V-or V * -equivalence classes has a polynomial upper bound. The set of V-equivalence classes (i.e., the event vectors) provides already a first information about the co-evolutionary history of the hosts and their parasites. Indeed, a high number of cospeciations may indicate that hosts and parasites evolved together, while a high number of host-switches may indicate that the parasites are able to infect different host species. Under the scope of the E-equivalence relation, we are also interested in which parasites are associated to each event type (disregarding losses): by looking at the E-equivalence classes, we can get knowledge of the number as well as the locality of the events inside the parasite tree. The information from the E-equivalence classes can be further refined using the EL-equivalence classes, as they tell us not only where the host-switches happen but also which arcs of the parasite tree are involved (this information is relevant in the analysis of timefeasibility; a practical application can be found in Chapter , Section . . ).

The intuition behind the CD-equivalence relation is that, when a host-switch happens, there may be various hosts that can be selected as the parasite's "landing site".

In this case, we choose to consider as equivalent those reconciliations for which, while the hosts that receive the switching parasites may differ, all the other parasite/host associations (not corresponding to a host-switch) are the same. These reconciliations are similar and often indistinguishable without additional biological information. Indeed, take the two reconciliations φ and φ in Figure . : they are identical except for one switching parasite p , which is mapped to h by φ and to h a by φ . Since h and h a are two sibling nodes sharing the same parent in the host tree, without further information, there is no good way to tell apart the two reconciliations φ and φ , hence we consider them as equivalent.

Among the five equivalence relations, the CD-equivalence relation is the only one that captures the precise parasite/host associations for some internal parasite nodes.

At first glance, it might seem too restrictive, and one might fear that the number of CD-equivalence classes would be too large, which limits the practical usability. In Chapter , Section . , we will see that, for real-world datasets, the number of CDequivalence classes is almost always much smaller compared to the number of optimal reconciliations. In fact, there is also a theoretical foundation for allowing the host-switch associations to be "movable": as we will see in Lemma . , finding one optimal reconciliation means that (a potentially very large number of) certain other reconciliations must also be optimal, all of them differ only by the associations of the host-switch nodes.

Lemma . . (The proof is straightforward and can be found in my Master thesis). Let φ be

an optimal reconciliation of (H, P, σ) with a cost vector ì c satisfying c(ℂ) < c( ) + c( ) (this is verified by any reasonable cost vector used in practice). Let p ∈ V(P) be an internal node with children (p , p ). Suppose E φ (p ) = and suppose wlog that (p , p ) is the host-switch

arc at p . Suppose that φ(p ) ≠ φ(p ) (so necessarily φ(p ) is a proper ancestor of φ(p ), see Figure .
). Take any node h on the path from φ(p ) to φ(p ) in H, excluding φ(p ) (there is at least one choice for h which is φ(p )). Consider φ ′ :

V(P) → V(H) such that φ ′ (p ) = h, and ∀p ∈ V(P) \ {p }, φ ′ (p) = φ(p)
(φ ′ is identical to φ except for the association at p ). Then φ ′ is an optimal reconciliation and is CD-equivalent to φ. . .

p 0 p 2 p 1 φ(p 0 ) φ(p 2 ) d H ⩾ 1

Relationship between the equivalence relations

From the Definitions . -. , we can easily deduce the following:

. If φ and φ are V-equivalent, then they are V * -equivalent.

. If φ and φ are E-equivalent, then they are V * -equivalent.

. If φ and φ are EL-equivalent, then they are E-equivalent.

. If φ and φ are CD-equivalent, then they are E-equivalent.

The next result is only guaranteed when φ and φ are optimal reconciliations for a cost vector ì c satisfying c( ) ≠ (nonzero loss cost):

. If φ and φ are V * -equivalent, then they are V-equivalent.

The converse of the above five statements is generally not true; counterexamples can be found easily.

In practice, we will often look at the number of equivalence classes of a fixed instance (i.e., the size |π(T(G))| of the set of equivalence classes, see Definition . ) for several different equivalence relations. If we know that being A-equivalent implies being B-equivalent, where A and B are two equivalence relations, then we know that the number of A-equivalence classes is no bigger than the number of B-equivalence classes. If the experimental observation contradicts this, there is a strong indication that something went wrong (for example, we did not compute the things that we thought we did). Therefore, knowing the relationship between the equivalence relations, apart from the obvious usefulness of helping us to better understand the meaning of equivalence classes, is also useful during experiments as a safeguard measure.

. .

Computational problems

Equipped with our definitions of equivalence classes (see Definition . ), we aim at studying the features of the set of optimal reconciliations by enumerating the equivalence classes (see Chapter for some other approaches, and a comparison with the equivalence-based approach). Naively, one would iterate through every reconciliation and record their properties, then report the equivalence classes, and, only at the end, report the statistics of the reconciliations in each equivalence class. However, when the number of reconciliations is too large, for example, in the order of (see Chapter , Section . . and [Wan+ ]), the naive method is not feasible.

The challenge is then to enumerate directly the equivalence classes without enumerating the optimal reconciliations explicitly. Concretely, the set of optimal reconciliations will be represented implicitly as T(G), the set of solution subtrees of a reconciliation graph G. Given a reconciliation graph as input, we will tackle the following problems:

• Count the number of equivalence classes.

• Enumerate the equivalence classes.

• Study a particular equivalence class. That is, given an equivalence class, -Count the number of reconciliations in that class, -Find one representative (i.e., one optimal reconciliation) of that class, -Enumerate all reconciliations in that class.

. Dealing with the V-and V * -equivalence relations

. . V-equivalence relation

The enumeration of the V-equivalence classes (i.e., event vectors) can be achieved by modifying the dynamic programming algorithm in such a way that the set of V-equivalence classes for each DP subproblem is kept along the way. This basic idea can, in theory, be applied to any equivalence relation. However, such a modified algorithm will require explicitly storing all the equivalence classes. Fortunately, for the V-equivalence relation, the space complexity will remain a polynomial of the size of the input trees, as we will see next. possible combinations. The number of losses ξ φ (p) for each parasite node p under any reconciliation φ is at most twice the diameter of the host tree (i.e., twice the maximum distance between two nodes), so the fourth element of any event vector is bounded by O(nm). Therefore, the number of event vectors is bounded by O(nm ).

Given an instance and its reconciliation graph G, we are interested in the following two problems: listing all event vectors, and, given a particular event vector, listing one (or all) optimal reconciliations in that V-equivalence class. We will now describe how to modify the DP algorithm for constructing the reconciliation graph G to solve both problems.

Recall that if w is an OR + node of the reconciliation graph G, the set of solution subtrees T(G/{w}) of the subgraph starting from w corresponds to the optimal reconciliations of the DP subproblem identified by the cell (p, h, e) with which w is labeled. We define the set EV(w) of an OR + node w to be the set of V-equivalence classes of T(G/{w}) (i.e., the set π(T(G/{w})), where π maps a reconciliation to its event vector, see Definition . ). Then, the sets EV can be computed as follows (for simplicity, we will identify an OR + node with the cell (p, h, e) with which it is labeled):

• For each goal node (p, h, ), EV(p, h, ) := {( , , , }.

• For each OR node (p, h, e), let { (p i , h i , e i ), (p i , h i , e i ) } ≤i≤k be its set of grandchild couples, then EV(p, h, e) can be computed as

≤i≤k ì u ∈EV(p i ,h i ,e i ) ì u ∈EV(p i ,h i ,e i )                    ì u + ì u + ( , , , ξ(h, h , h )) +                    ( , , , ) if e = ℂ ( , , , ) if e = ( , , , ) otherwise, e =                    . ( . )
The set of event vectors of T(G) that we seek is the union w EV(w) taken over the set of start nodes of G, i.e., the OR + nodes labeled with a cell of the form (r(P), h, e).

For each of the O(nm) OR + nodes of the reconciliation graph, we need to keep an extra set of size O(nm ). The space complexity is therefore O(n m ). For each OR node and for each of its O(n ) grandchild couples, we need to compute the Cartesian sum of two sets of EV's of size O(nm ) each; this can be done naively in time O(n m ) (to improve this, see, e.g., [FJ ]). The overall time complexity is O(n m ).

Given an event vector, finding one optimal reconciliation in that V-equivalence class can be done using a standard technique called backtracking. Let us define a function

Backtrack that takes two parameters: an OR + node w of the reconciliation graph and a vector ì v satisfying ì v ∈ EV(w). The function outputs an optimal subproblem reconcil-

iation φ w ∈ T(G/{w}) such that ì e (φ w ) = ì v.
We choose to represent a reconciliation as a sequence of triples of the form (p, h, e), and the function will output one triple in each recursive step. The function Backtrack(w, ì v) can be implemented as follows:

. Let (p, h, e) be the cell with which w is labeled. Output the triple (p, h, e). If w is a goal node, stop. Otherwise, go to step .

. Let { (p i , h i , e i ), (p i , h i , e i ) } ≤i≤k be the grandchild couples of w. Find any index i such that there exist ì u ∈ EV(p i , h i , e i ) and ì u ∈ EV(p i , h i , e i ) such that the sum inside the big braces of Equation ( .) is equal to ì v (such an i necessarily exists). Choose any such ì u and ì u . Then, call Backtrack((p i , h i , e i ), ì u ) and

Backtrack((p i , h i , e i ), ì u ).
Given a start node w and an event vector ì v ∈ EV(w), it suffices to call Backtrack(w, ì v)

to get one representative of the V-equivalence class ì v. Finally, if we replace "any" by "all" in step . of Backtrack, we can easily adapt the algorithm so that it enumerates all reconciliations, or counts the number of reconciliations in a V-equivalence class.

. . V * -equivalence relation

It is easy to see that the method described in the previous section for dealing with the V-equivalence relation can be applied directly to the V * -equivalence relation. For finding all starred event vectors, the extra spaced needed for storing those is only O(nm ),

as the total number of starred event vectors is bounded by O(m ); the time complexity is O(n m ).

.

Dealing with the E-, EL-, and CD-equivalence relations

.

. E-equivalence relation

By Definition . , an E-equivalence class is a function that maps each parasite node p to its color of the form (p, e), where e is an event type e ∈ E := { , ℂ, , }. In this section, we will represent an E-equivalence class as a set T of colors of the form (p, e).

In the same manner, a reconciliation φ ∈ T(G) will be written as a set of triples of the form (p, h, e). Using this notation, a reconciliation φ is in an E-equivalence class T (i.e.,

π(φ) = T ) if
∈ T(G) such that π(φ) = T .
The problem of studying a particular E-equivalence class is easy. Given an E-equivalence class T , the reconciliation graph G can be pruned in such a way that its set of solution subtrees corresponds precisely to the reconciliations in the class T (we simply need to remove all OR nodes unless its label (p, h, e) corroborates the given class, i.e., (p, e) ∈ T ).

Counting and enumerating the reconciliations in that particular E-equivalence class can then be done using exactly the same method as for the "full" reconciliation graph.

Counting and enumerating the E-equivalence classes are, however, more challenging problems. At present, we will first concentrate on the enumeration of the set π(T(G))

of E-equivalence classes.

Enumerating the E-equivalence classes

The algorithm is based on the simple idea of traversing the reconciliation graph in a top-down fashion (a similar approach can be used in an algorithm that enumerates all the solution subtrees). In order to obtain a polynomial time delay algorithm, during the traversal, we can no longer consider the nodes one by one; the sets of nodes that are contained in E-equivalent reconciliations must be traversed together. To make this clear, we will define conveniently the color of an OR + node in the reconciliation graph; the algorithm will then traverse nodes of the same color together.

Definition . (Color of an OR + node. Color couple).

• Let w be an OR + node labeled by a cell (p, h, e) ∈ U, we say that w is colored by the ordered pair κ(w) := (p, e) ∈ V(P) × E.

• Let w and w be two OR + nodes colored respectively by (p , e ) and by (p , e ). The color couple of the couple of nodes (w , w ) is the couple of colors ((p , e ), (p , e )).

To The pitfall of the naive approach is that not every combination between the Eequivalence classes of the reconciliations of the two child subtrees is valid. Our algorithm, shown in Algorithm , can be viewed as an improved version of the naive algorithm in which particular care has been taken to ensure that only valid combinations are outputted. Along with each E-equivalence class T , it also outputs a set O which is a subset of the input set O: it is equal the union of the root OR + nodes of all solution subtrees φ ∈ T(G/O) such that π(φ) = T . Notice that in Algorithm we employ both the return and the yield statements for the output, the difference being that the latter does not halt the algorithm.

Before the proof of correctness, let us recall that we use interchangeably the terms optimal reconciliations and solution subtrees. We denote by r(φ) the root node of the solution subtree φ. Observe that, if φ is a reconciliation of (H, P, σ), then r(φ) is a start node; if φ is an optimal solution of the DP-subproblem indexed by a cell (p, h, e) (that is, a reconciliation of P| p ), then r(φ) is an OR + node colored by (p, e).

Lemma . . The function

O = φ {r(φ) | π(φ) = T , φ ∈ T(G/O)}. ( . )
Proof. The proof is by induction on the height h p of the subtree P| p . We use the fact that the pre-condition in the Require statement in Algorithm is true for all recursive calls of Enumerate (easy induction). When h p = , p is a leaf and {(p, σ(p), )} is the

Counting the E-equivalence classes

The next algorithm is a heuristic: it computes the exact value of |π(T(G))| without a guaranteed upper bound on the running time or on the space. The idea is to build an ad-AND/OR graph with a structure similar to that of the reconciliation graph and whose solution subtrees correspond bijectively to the E-equivalence classes; after this step, counting and enumerating the E-equivalence classes can be done using the standard techniques for ad-AND/OR graphs. Just as the reconciliation graph, such a graph can have a relatively small size even when the number of solution subtrees is huge (say, arcs in the graph versus solution subtress). In Chapter , Section . . , we will see that this heuristic method does allow us to compute some huge numbers of E-equivalence classes.

Exactly like the reconciliation graph, the graph is constructed in a bottom-up fashion using the dynamic programming technique. The DP-table is denoted by D. Each step is indexed by a cell (p, h, e) ∈ U. The space of cells is traversed following the postorder of the parasite tree. Notice that the input of the algorithm will be the reconciliation graph. For simplicity, we use the notation M(p, h, e) to denote the unique OR + node in the reconciliation graph labeled by (p, h, e) (see Section . . , and also Equations ( . ) and ( . )). At the end of the step (p, h, e), the value D(p, h, e) of the DP-table

will also be an OR + node: instead of being labeled by a triple (p, h, e), this node will be labeled by a couple (p, e) (i.e., a color, see Definition . ).

We will work on OR + nodes while keeping in mind the following: an OR + node w induces a complete ad-AND/OR graph structure containing all nodes that are reachable from w and in which w is the unique start node. We denote by T(w) the set of solution subtrees of the ad-AND/OR graph induced by the OR + node w. This "new" notation should be intuitive: if w is an OR + node in a given graph G, the set T(w) was denoted by T(G/{w}). When we consider an OR + nodes for which the set of solution subtrees corresponds to a set of reconciliations (that is, an OR + node of the reconciliation graph), we denote by π(T(w)) the set of E-equivalence classes (see also Equation ( . ), Definition . ):

π(T(w)) := φ∈T(w) π(φ). ( . )
Either for a reconciliation or for an E-equivalence classes, we will use its correspondence with the label set (or a set of colors) of a solution subtree: by looking at the labels of the OR + nodes in each solution subtree in T(w), we obtain a set of labels of the form (p, h, e) or (p, e). As we mentioned in the beginning of this section, a set of colors T is a valid E-equivalence class for the subproblem (p, h, e), i.e., T ∈ π(T(M(p, h, e))), if and only if there exists φ ∈ T(M(p, h, e)) such that π(φ) = T .

Using this notation, we know that the OR + node M(p, h, e) of the reconciliation graph satisfies the following:

• T(M(p, h, e)) corresponds bijectively to the set of optimal reconciliations of the DP-subproblem indexed by (p, h, e) (those are reconciliations φ of the subtree P| p satisfying φ(p) = h and E φ (p) = e)).

To obtain an ad-AND/OR graph as desired, the OR • (Pre-condition ) Every w i has the same label,

• (Pre-condition ) For every grandchild couple (w i , w i ) ∈ gch(w i ), both T(w i ) and T(w i ), as multisets of label sets, have no redundancies.

It returns a single OR node w m such that

• (Post-condition )
We have the equality between the two sets of label sets

T(w m ) = k i= T(w i ) , ( . )
• (Post-conditions ) T(w m ), as a multiset of label sets, has no redundancies. Now, we are ready to describe how to fill up the DP-table D. Since the space of cells is the same, and the cells are filled in the same order, the algorithm is very similar to the construction of the reconciliation graph (see Section . . ). For each cell (p, h, e) ∈ U such that M(p, h, e) exists, in the post-order of V(P):

• If p is a leaf, D(p, h,
π(φ i ) = T i , for i = , . Let T := T ∪ T ∪ {(p, e)}.
It is clear that π(φ) = T . For T to be a solution subtree, i.e., for T to be in T(D(p, h, e)), it suffices to find an OR node w i created during the construction step (before applying Merge, see also Lemma . ) such that T ∈ T(w i ). Such a w i can be found, since there exists a child AND node of M(p, h, e)

having (M(p , h , e ), M(p , h , e )) as children (it is the child AND node of M(p, h, e) in the solution subtree φ).

It remains to describe the Merge function. The presentation is complicated, and is deferred to the end of this chapter in Section . as Supplementary material.

. .

CD-equivalence relation

. . EL-equivalence relation

To apply the same methods as for the E-equivalence relation, we can build a slightly modified version of the reconciliation graph which keeps tracks of the host-switch arcs. Consider a new set of event types

E ′ := { , ℂ, , L , R } (see Definition . ).
Given a reconciliation graph G, we will modify G so that each OR + node, instead of being labeled by (p, h, e) ∈ V(P) × V(H) × E, will now be labeled by

(p, h, e ′ ) ∈ V(P) × V(H) × E ′ .
To do this, for every OR node w in G labeled by a cell (p, h, e) such that e = , in the post-order of V(P), we perform the following operations:

. (to duplicate a node means to copy all the incoming and outgoing arcs). Remove the arcs between w L and some of its children until gch(w L ) = gch L . Remove the arcs between w R and some of its children until gch(w R ) = gch R .

After this step (which takes linear time in the size of the graph), any method for the E-equivalence relation can be used directly without modification.

.

Supplementary material: the Merge function

The goal of this section is to present the Merge function mentioned in Lemma . .

In this section, we do not work on the reconciliation graph but on some general ad-AND/OR graph represented by its start node w, and such that each of its OR + node is labeled by a color of the form (p, e). We will be interested in questions such as "does T(w), as a multiset of label sets, have redundancies?" In order to clarify and simplify this statement, we start by introduce some new notations.

. . Notations and definitions

If T ∈ T(w) is a solution subtree, we denote by C(T ) the set of colors of the OR + nodes in T (i.e., the label set). What we have been calling the multiset of label sets is denoted

as the multiset C(w) := {C(T ) | T ∈ T(w)}.
Instead of talking about label sets or sets of colors, we will consider colorings, which we define next. Indeed, for any w that we will encounter in the algorithm, a label set in C(w) always coincides with a coloring of a parasite subtree.

We say that an OR + node w colors a node p ∈ V(P) if its color κ(w) = (p, e) for some e ∈ E. We will also say that an OR + node colors a parasite subtree, using the following recursive definition:

Definition . . Let p ∈ V(P). Let w be an OR + node that colors p. We say that w colors the subtree P| p if it satisfies the following:

• If p is a leaf, then w is a goal node.

• If p is an internal node with children (p , p ), then w is an OR node, and for each grandchild couples (w , w ) ∈ gch(w), w colors P| p and w colors P| p .

Given a tree and a fixed set of colors, a coloring (or node-coloring) of the tree is a function that maps each node of the tree to a color. If w colors P| p , then each solution subtree T ∈ T(w) corresponds to a coloring of the subtree P| p ; this coloring, which can be written as a label set, is also denoted by C(T ). Notice that each coloring in the multiset C(w) has the same root color, i.e., the color assigned to the root p of the subtree P| p (the root color is simply κ(w)). The next important definition replaces the statement "the multiset C(w) has no redundancies" by a single adjective:

Definition . . Let p ∈ V(P) and let w be an OR + node that colors P| p . We say that the node w is neat if no two colorings in the multiset C(w) are the same, in other words, for any two solution subtrees T and T ∈ T(w),

C(T ) = C(T ) ⇐⇒ T = T .
The following result is easy: if an OR node w is neat, and if (w , w ) ∈ gch(w), then both OR + nodes w and w are neat.

Using the new notations, Lemma . can be re-written as follows:

Lemma . (Reformulation of Lemma . ). The Merge function takes as input a list w , . . . , w k of OR + nodes satisfying:

• (Pre-condition ) Every w i colors the same subtree P| p with the same root color, for some node p ∈ V(P),

• (Pre-condition ) Every w i is neat.

It returns a single OR + node w m such that

• (Post-condition )
We have the equality between the two sets of colorings

C(w m ) = k i= C(w i ). ( . ) • (Post-condition ) w m is neat.
Notice that the union operation in Equation ( .) is the union for sets; by the precondition , every C(w i ) is in fact a set (multiset whose elements are of multiplicity one); the left hand side is also a set thanks to the post-condition .

. . Description of Merge

Merge

Informally speaking, the Merge function has two objectives: completeness (we must capture all colorings of the input) and neatness (we must not produce duplicates). While the second objective seems nontrivial, there is an easy method for achieving the first objective, or Equation ( . ) (it is clear that any candidate output w m should have the same color as any input w i ; we will omit the color of w m and focus on its children or its grandchild couples; we also omit the trivial case where the function applies to goal nodes instead of OR nodes):

• (Step of Merge) Create an OR node w m and connect it to every one of the child AND nodes of the input w i 's, so that the set of children of w m is the union of the children of the input w i :

ch(w m ) = k i= ch(w i ).
This is exactly what we decided to do in Merge, followed by a second step:

• (Step of Merge) Call MakeNeat(w m ) so that w m becomes neat.

We still need to describe MakeNeat, but we are one step closer, because, instead of trying to capture all the colorings of the (multiple) input graphs, we now only need to modify one single graph while maintaining the same set of distinct colorings.

Lemma . . The MakeNeat function takes as input a single OR node w satisfying:

• (Pre-condition ) w colors a subtree P| p , for some node p ∈ V(P).

• (Pre-condition ) For every grandchild couple (w , w ) ∈ gch(w ), both w and w are neat.

It returns a single OR node w m such that

• (Post-condition ) w m is neat. • (Post-condition ) C(w m ) is equal to the set of distinct colorings in C(w ).

Notations

To give a characterization of the neatness of a node, we need some more notations.

Notice that a grandchild couple (w , w ) ∈ gch(w ) corresponds to an AND node in ch(w ). Since an AND node is uniquely identified by its two child OR + node, in this section, we simply denote an AND node by a couple of OR + nodes, and we also write (w , w ) ∈ ch(w ). The algorithm will create AND nodes that are not children of some already fixed OR nodes; their common properties are related to colorings of a parasite subtree.

Definition . . Let p ∈ V(P) \ L(P) with children (p , p ). Let (w , w ) be an AND node.

We say that (w , w ) neatly colors the subtree P| p if w colors P| p and w colors P| p , and both w and w are neat.

If an AND node (w , w ) neatly colors P| p , we denote by w × w the set

w × w := C(w ) × C(w ).
An element C of w × w is called a coloring of P| p with undetermined root color. In a similar spirit as for the Cartesian products of colorings, we next introduce shorthands for the intersection and the difference of sets of colorings. If w and w are two neat OR + nodes that color P| p , we denote, respectively, by w ∩ w and w \ w the sets

w ∩ w := C(w ) ∩ C(w ), w \ w := C(w ) \ C(w ).

MakeNeat

We can now give a characterization of the neatness of a node and present a strategy for the MakeNeat function. Let w be the input of MakeNeat. It is an OR + that colors P| p for some p ∈ V(P). If p is a leaf, then w is necessarily neat. Otherwise, that is, if p is an internal node:

• If, for every two child AND nodes of w , indexed by i and j, and denoted respectively by (w i , w i ) and (w j , w j ) ∈ ch(w ) (by the pre-conditions, all child AND nodes neatly color P| p ), we have

i ≠ j =⇒ w i × w i ∩ w j × w j = ∅, ( .
)
then w is neat.

To put it short, an OR node is neat if the sets of colorings (with undetermined root color) of its child AND nodes have pairwise empty intersections.

A straightforward strategy for MakeNeat is to repeatedly find an "intersecting" pair of child AND nodes and to make them "intersection-free":

• Function MakeNeat(w ): while there is a pair of AND nodes (w i , w i ) and (w j , w j ) in ch(w ) such that w i × w i ∩ w j × w j ≠ ∅, replace those two nodes by MakeInter-sectionFree((w i , w i ), (w j , w j )).

Lemma . . The MakeIntersectionFree function takes as input two AND nodes (w i , w i ) and (w j , w j ) such that both input nodes neatly color P| p for some p ∈ V(P), and

w i × w i ∩ w j × w j ≠ ∅. ( . )
It returns a set A of AND nodes (w k , w k ) k , all of which neatly color P| p , and such that

k w k × w k = w i × w i ∪ w j × w j , ( .
)
and the sets of colorings w k × w k of nodes in A have pairwise empty intersections.

Table . shows what we have chosen as the output of MakeIntersectionFree depending on the cases (some easy set-theoretical computation will suffice for checking that the output verifies Equation ( . ) and the condition of pairwise empty intersection). The cases represent the different causes for two Cartesian products of sets to be have nonempty intersection (see Equation ( .)). The first two cases (a) and (b) utilizes the Merge function (here, Merge is used recursively, as it is called on OR + nodes that color a subtree P| p ′ of smaller height than P| p , the latter is colored by the input nodes of the Merge call from which the current call of MakeIntersectionFree originates); the output in these two cases is only one AND node. In cases (c) -(f), the output consists of two AND nodes (however, one is potentially empty, see the paragraph after Lemma . below). Notice that the four cases (c) -(f) are not not mutually exclusive (all the cases are checked sequentially until the first applicable case is found). In the most general case, that is, case (g), all four set differences w i \ w ), which we will present in the next section.

Case

Output AND nodes Notes

(a) C(w i ) = C(w j ) w i , Merge(w i , w j ) (b) C(w i ) = C(w i ) Merge(w i , w j ), w i (c) C(w i ) ⊆ C(w j ) (w i , w diff ) and (w j , w j ) C(w diff ) = w i \ w j (d) C(w j ) ⊆ C(w i ) (w i , w i ) and (w j , w diff ) C(w diff ) = w j \ w i (e) C(w i ) ⊆ C(w j ) (w diff , w i ) and (w j , w j ) C(w diff ) = w i \ w j (f) C(w j ) ⊆ C(w i ) (w i , w i ) and (w diff , w j ) C(w diff ) = w j \ w i (g) none of the above applies (w diff , w i ) and (w inter , w diff ) and (w j , w j ) C(w diff ) = w i \ w j C(w inter ) = w i ∩ w j C(w diff ) = w i \ w j
C(w inter ) = w ∩ w , C(w diff ) = w \ w . ( . )
In the case of empty intersection or difference, the output is a special symbol Empty.

In cases (c) -(f), it can happen that the "diff" node in the output is Empty. This corresponds to the situation where both cases (c) and (e) or both (d) and (f) apply. Intuitively, we can say that one input AND node is entirely included in the other one as sets of colorings. In this case, the output is just one single AND node, and it is equal to one of the input nodes (the "bigger" one, or the one that includes the other).

Notice that the MakeIntersectionFree function needs to test for the inclusion and the equality between sets of colorings. The equality is checked by double inclusion (an improved method is discussed at the end of Section . . ). The inclusion can be checked by computing the intersection and the difference:

C(w ) ⊆ C(w ) ⇐⇒ (w ∩ w ≠ ∅) ∧ (w \ w = ∅).
. .

Intersection and difference

In this section, an AND node can either be denoted by a single letter v or by a couple (w , w ) of OR + nodes that represents its children. If v := (w , w ) neatly colors P| p , we denote its set of colorings with undetermined root color by

C(v) = w × w := C(w ) × C(w ).

Preliminaries

˙ v∈A inter C(v) should be equal to ˙ v i ∈A C(v i ) ∩ C(v j ) = ˙ v i ∈A C(v i ) ∩ C(v j ) , ˙ v∈A diff C(v) should be equal to ˙ v i ∈A C(v i ) \ C(v j ) = ˙ v i ∈A C(v i ) \ C(v j ) .
( . )

All the union operations in the above equations are decorated with a small dot, which

indicates that those are unions on pairwise disjoint sets. To go from the left hand side to the right hand side of the equal signs in Equation ( . ), we have used the right distributive property of the intersection and the difference over the union.

AND-InterDiff

The function AND-InterDiff(A, v j ) can be implemented as follows:

• If A = ∅, return two empty sets.

• If |A| > , for each v i ∈ A, call AND-InterDiff on the singleton set {v i } and v j .

Return the union of all the outputted A inter , and the union of all the outputted

A diff (we are using the right distributive property, see Equation ( .)). Since all the union operations are on pairwise disjoint sets that contain pairwise disjoint AND nodes, a union can simply be implemented as putting AND nodes to a list.

• In the remaining case, let v i be the only element of A. We have two AND nodes v i := (w i , w i ) and v j := (w j , w j ) such that both neatly color P| p for some p ∈ V(P),

and we want to compute their intersection and difference. To do this, we will use the InterDiff function for two OR + nodes. Let w inter and w diff be the output of InterDiff(w i , w j

). Let w inter and w diff be the output of InterDiff(w i , w j

).

-If neither w inter nor w inter is Empty, return A inter := {(w inter , w inter )}. Otherwise, return A inter := ∅.

-( ) If w diff is not Empty and w diff is Empty, return A diff := {(w diff , w i )}.

( ) If w diff is Empty and w diff is not Empty, return A diff := {(w i , w diff )}.

( ) Otherwise, that is, neither w diff nor w diff is Empty, ( a) If either w inter or w inter is Empty (or both), return

A diff := {(w i , w i )}. ( b) Otherwise, return A diff := (w diff , w i ), (w inter , w diff ) .
Similar to what we did in MakeIntersectionFree, the different cases of AND-InterDiff can be understood by asking in which situations two Cartesian products have nonempty set difference; the correctness of the output can be verified by some easy set-theoretical computation. Notice that only the case ( b) requires us to check carefully that the two AND nodes in A diff have empty intersection.

InterDiff

Algorithm shows the function InterDiff. The idea that we mentioned in the beginning of this section corresponds to the for loop at Lines -. Every occurrence of the ∪ operator corresponds to the union between disjoint sets of AND nodes, and the AND nodes in the same set always have pairwise empty intersection. The variable A i records the difference between v i , the i-th child AND node of w and the union of the first j child AND nodes of w (here, the sets of children have an arbitrary ordering), and is equal to A diff i,j during the inner loop (Line ), for ≤ j ≤ |ch(w )|. We have, at Line , if j > :

C(A diff i,j ) = C(A diff i,j-) \ C(v j ) = C(v i ) \ ≤k≤j- C(v k ) \ C(v j ) = C(v i ) \ ≤k≤j C(v k ) . If j = , C(A diff i,j ) = C(v i ) \ C(v j ).
The intersection between the i-th child v i of w and the j-th child of w can also be expressed using the intersection between v i and the first j children of w (before computing the intersection with v j , we first subtract some previous child AND nodes that are disjoint with v j ; this does not change the intersection).

We have, at Line , if j > :

C(A inter i,j ) = C(A diff i,j-) ∩ C(v j ) = C(v i ) \ ≤k≤j- C(v k ) ∩ C(v j ) = C(v i ) ∩ C(v j ). If j = , C(A inter i,j ) = C(v i ) ∩ C(v j ).
At the end of the inner loop, the following sets are added (using unions of disjoint sets) to A inter and A diff that collects all the intersection and difference between v i , the i-child of w , and all the children of w :

|ch(w )| j= C(A inter i,j ) = |ch(w )| j= C(v i ) ∩ C(v j ) = C(v i ) ∩ |ch(w )| j= C(v j ) C A diff i,|ch(w )| = C(v i ) \ |ch(w )| j= C(v j ) .
The above observations, together with the right distributive property of intersection and union, show that we indeed get the correct sets of colorings for A inter and A diff , at the end of the outer loop (Line ):

C(A inter ) = |ch(w )| i= |ch(w )| j= C(A inter i,j ) = |ch(w )| i= C(v i ) ∩ |ch(w )| j= C(v j ) C(A diff ) = |ch(w )| i= C A diff i,|ch(w )| = |ch(w )| i= C(v i ) \ |ch(w )| j= C(v j ) .
These two sets consist of colorings with undetermined root color. At the end of Algorithm , we add the correct root color and create the two OR nodes having respectively A inter and A diff as sets of children. We then obtain the desired result as stated in Lemma . , Equation ( . ):

C(w inter ) = w ∩ w , C(w diff ) = w \ w .

. . Another operation on a set of AND nodes

We have already finished the description of all the components of the Merge function.

The idea was quite simple: after a "naive merge" of the input AND nodes, we remove all the duplicates in the set of colorings by making every pair of AND nodes "intersectionfree". This basic operation, called MakeIntersectionFree, can be easily implemented

A i ← {v i }
for each child AND node v j ∈ ch(w ) do

Let A inter i,j and A diff i,j be the output of AND-InterDiff(A i , v j ) after repeatedly applying MakeIntersectionFree so that there are no intersecting pair of AND nodes, RemovePartners can be used to reduce the number of children of the output OR node.

A inter ← A inter ∪ A inter i,j A i ← A diff i,j end A diff ← A diff ∪ A i end
Given a set A of AND nodes, the RemovePartners operation repeatedly finds a pair of AND nodes satisfying certain conditions, called a partner pair, and replaces them by a single AND node, and it stops when no partner pair can be found. By assumption, any two distinct AND nodes (w i , w i ) and (w j , w j ) in A satisfy the condition of being disjoint as sets of colorings, that is,

w i × w i ∩ w j × w j = ∅. ( . )
Table . shows the two cases in which two AND nodes in A form a partner pair, and

gives the new AND node that replaces them. The ∪ symbol means to perform a "naive merge" between two OR nodes of the same color by taking the union of their child AND nodes (this is the same as the first step of Merge; the second step of Merge is unnecessary here, because the sets of colorings of the two OR nodes are necessarily disjoint).

Notice that the first two cases of MakeIntersectionFree also deal with similar situations (there, however, the second step of Merge cannot be omitted).

Case

To be replaced by

C(w i ) = C(w j ) and κ(w i ) = κ(w j ) (w i , w i ∪ w j ) C(w i ) = C(w j ) and κ(w i ) = κ(w j ) (w i ∪ w j , w i )
Table . : The two cases in which two AND nodes (w i , w i ) and (w j , w j ) satisfying Equation ( . ) form partner pair, and the corresponding replacement node in the RemovePartners function.

In RemovePartners, we need to check for the equality between the sets of colorings of two OR nodes. Unlike MakeIntersectionFree which checks the equality using the double inclusion criterion (by computing the intersection and the difference between several pairs of OR + nodes), RemovePartners, being an optional operation by nature, does not require an exact equality check with double inclusion. A sufficient condition for C(w ) = C(w ) is that the two graphs are isomorphic with an isomorphism that preserves the colors. This condition can be checked in linear time in the size of the graph (or faster in practice using hashing techniques), as we can see in the following recursive definition.

Definition . . Given two neat OR + nodes w and w that color P| p for some p ∈ V(P),

we say that they are sufficiently equal if κ(w ) = κ(w ) and, either

. p is a leaf (so w and w are goal nodes), or

. w and w have the same number of children, and ch(w ) and ch(w ) can be ordered in such a way that for each i = , . . . , |ch(w )|, the i-th child (w , w ) of w and the i-th child (w , w ) of w satisfy

• w and w are sufficiently equal, and

• w and w are sufficiently equal.

In the implementation, MakeIntersectionFree can also gain an improved performance by checking the equality (for the first two cases) in two phases: first we check the sufficient equality, and only if the first test fails, we check the double inclusion.

Chapter

Phylogenetic tree reconciliation: practical contributions

.

Experimental results

We implemented in Python the algorithms for enumerating and counting the equivalence classes from Chapter . In this section, we apply them on biological datasets and give some experimental results that demonstrate their usefulness in practice. The experiments are run on a laptop PC with Intel i -M CPU ( . GHz, cores) and GB RAM.

. . Instances

An instance of the R consists of a dataset, that is a triple (H, P, σ), where H and P are phylogenetic trees and σ : L(P) → L(H) is a function, and a cost vector ì c := (c(ℂ), c( ), c( ), c( )) that gives the cost value that we associate with each event type among cospeciation, duplication, host-switch, and loss.

Each dataset that we used is identified by an alphanumeric code. The code, the size (measured as the number of nodes in the trees H and P), and the source of each dataset are given in Table . . All the datasets have already been used in the literature (see [BAK ; Don+ ]) and together cover different situations (sizes and topologies of the trees) and different contexts: the COG datasets are from a genes-species context, while all the other datasets are from a hosts-parasites or hosts-symbionts context.

For each dataset, we formed five instances using five different cost vectors: (-, , , ), for maximizing the cospeciation; ( , , , ), for minimizing the events that lead to incongruencies between the tree topologies; ( , , , ) and ( , , , ), where host-switches are more penalized; ( , , , ), which is a vector chosen only for theoretical purposes and does not penalize cospeciations and losses.

With datasets and cost vectors, we have a total of instances. Table . : All the datasets that were used in the experiments and their references.

Dataset

Introduction

The first set of experiments is concerned with the number of equivalence classes of the set of optimal reconciliations, for each equivalence relation (this number is denoted by |π(T(G))| in Chapter ). As we mentioned in Chapter , Section . . , for an equivalence relation to be of practical interest, we typically expect the number of equivalence classes to be significantly smaller than the number of optimal reconciliations, that is,

|π(T(G))| ≪ |T(G)|,
where ≪ denotes "much smaller than". Indeed, if the number of equivalence classes were close to the number of optimal optimal reconciliations, analyzing one representative per equivalence class would not be advantageous, as analyzing all optimal reconciliations would require a comparable level of human effort but provides more information.

Notice that, even without performing the enumeration of equivalence classes or the enumeration of the representatives, the number of equivalence classes itself can still be used as a diversity measure. Consider for example the following situation: if we have a dataset that yields two drastically different numbers of equivalence classes (one is very small, while the other one is very large), say the numbers of E-equivalence classes, for two "similar" cost vectors, say ( , , , ) and ( , , , ), this may indicate that the solution space (i.e., the set of optimal reconciliations) for this dataset is sensitive to the choice of the cost vector. The issue of sensitivity to event costs has been studied literature, in particular in [BAK ], where the authors have also been interested in the consistency of event type assignments when different cost vectors are used (those assignments correspond to the E-equivalence classes). A small number of E-equivalence classes would indicate that the consistency is high. In [BAK ], the authors tried to discover different event type assignments (i.e., different E-equivalence classes) by sampling the solution space uniformly at random. In comparison, our algorithm for enumerating the E-equivalence classes can answer the consistency question exactly.

Discussion

First, let us look at the third column of Table . , corresponding to the numbers of Vand V * -equivalence classes. We have shown in Section . . that these numbers are bounded by a polynomial of the input size. They are indeed always very small; the only exception is that, when the cost vector is ( , , , ), the number of V-equivalence classes can get large: this can be understood easily since, by setting the loss cost c( ) to zero, we allow more variability in the number of losses of the solutions. One can argue that, if the loss cost is set to zero, it is more reasonable to consider the V * -equivalence relation rather than the V-equivalence relation.

The numbers of E-and CD-equivalence classes are also very small except when the cost vector is ( , , , ), in this case those numbers can be relatively larger. For example, for the dataset COG and cost vector , , , ), there are optimal reconciliations but only CD-equivalence classes. This result suggests that the associations that are involved in host-switches (in this case, horizontal gene transfer events) can account for almost all the diversity of the solution space.

The numbers of EL-equivalence classes are generally larger than the other equivalence classes. For some of the instances, we can observe that while the number of E-equivalence classes is small, the number of EL-equivalence classes can be large and close to (or even equal to) the number of optimal reconciliations. For example, for the dataset PP and cost vector (-, , , ), there are only E-equivalence classes whereas all the optimal reconciliations are in different EL-equivalence classes. This indicates that, compared to the E-equivalence classes that tell us which nodes are involved in host-switches, the EL-equivalence classes, which distinguish which arcs of the parasite tree are involved in host-switches, can convey much more information.

To better measure and visualize how small the number of equivalence classes are relative to the number of optimal reconciliations, we computed, for instances having between and optimal reconciliations, a value that we call Reduction: it is equal to the number of equivalence classes divided by the number of optimal reconciliations. In Figure . , each x coordinate corresponds to an instance, where the instances are sorted in the increasing order of the number of optimal reconciliations (the denominator of Reduction); for each instance, we plotted three points that correspond to the Reduction values for the E-, EL-, and CD-equivalence relations. . .

Lists of equivalence classes or representatives

As we have seen in the previous section, for many of the biological instances, the number of equivalence classes is small enough so that we can enumerate them all. We will now see that the list of equivalence classes can provide useful information about the solution space. For this purpose, we use the dataset WOLB that has the following particularity: all the cost vectors that we used lead to a number of optimal reconciliations that is at least , which is too large for any exhaustive enumeration method. To gain insight into the solution space of this dataset, apart from the methods based on random sampling [BAK ; Don+ ], only a limited number of the preexisting exact methods can be applied in practice (e.g., computing the diameter of the solution space

[Haa+ ], computing the distribution of the pairwise distances [SML ], computing a median reconciliation [Ngu+ ], etc.). The enumeration of the equivalence classes or of the representatives of the equivalence classes is an addition to the list of exact methods, and we will see that it can be quite informative. Consider a metric m on the solution space that measures the differences in the event type assignments between reconciliations:

m(φ , φ ) := {p ∈ V(P) | E φ (p) ≠ E φ (p)} . ( . )
Then, knowing the list of E-equivalence classes and the size of each class, we can compute the diameter, the distribution of pairwise distances, as well as a median reconciliation of the solution space with respect to the metric m.

For the experiments, we choose the instance with the WOLB dataset and the cost vector ( , , , ). First, we look at the list of V-equivalence classes (i.e., the event vectors), shown in Table . . The table also indicates the size of each class as proportions of the solution space (i.e., the proportion of optimal reconciliations in each V-equivalence class among all optimal reconciliations).

From the list of event vectors, one can see that the dataset can be explained by a large number of host-switches and cospeciations, and that there have probably been no duplication. The diversity of the solution space in terms of the event vectors is low, not only because the number of different event vectors is small, but also because these event vectors are very similar to each other. Furthermore, if we view the proportions of the event vectors as probability densities, the distributions of the solutions in the spaces of all four dimensions of the event vectors are unimodal. For example, along the axis of the number of cospeciations, solutions that have (the most common value)

cospeciations or a number close to are the most probable; in contrast, we are less likely to find solutions that have a number of cospeciations which is far away from (extreme values). Now, we turn to the list of the E-equivalence classes. Recall that an E-equivalence class is a function π that associates each parasite node p ∈ V(P) with a color of the form (p, e), where e ∈ E is an event type. One way to measure how dissimilar two Eequivalence classes are from each other is to take the sum of differences:

m(π , π ) := |{p ∈ V(P) | π (p) ≠ π (p)}| . ( . )
Let s i denote the size of an E-equivalence class π i (i.e., the number of optimal reconciliations in that class). Equation ( . ) defines a measure on the set of E-equivalence classes which can be used to compute, for instance, the average pairwise distance of the solution space with respect to the metric in Equation ( . ) by noticing the following relationship:

i≠j m(φ i , φ j ) = i≠j m(π i , π j ) • s i • s j , ( . )
where the first sum is over the set of optimal reconciliations, the second sum is over the set of E-equivalence classes. We computed the distances between the E-equivalence classes. The maximum distance is . The average pairwise distance of the solution space (see Equation ( . )) is around . Informally, it means that if we randomly select two optimal reconciliations, we expect to find two differences in their node-wise event type assignments. Those distance values are extremely small compared to the theoretical maximum of (the number of internal nodes of the parasite tree). Another observation that can be made from the list of E-equivalence classes is about the parasite nodes that can receive different event types across the solution space. Instead of looking at the total number of differences, we can ask which nodes contribute to the differences. There are only of them, in other words, all other internal nodes receive a consistent event type across the entire solution space. We have further confirmed that the diversity of the solution space is low: not only the event vectors are similar, the distributions of the events on the nodes of the parasite tree are also similar.

For the CD-equivalence classes, we can ask the similar question of which nodes contribute to the differences. Here, a parasite node p contributes to the differences be-

. . Finding time-feasible reconciliations using EL-equivalence classes

In Chapter , Section . . , we mentioned that a biologically meaningful reconciliation should satisfy the requirement of time-feasibility. Unfortunately, the T -, i.e., finding a time-feasible reconciliation that minimizes the cost, is NP-hard [Ova+ ; THL ]. A practical way of obtaining time-feasible reconciliations is to first enumerate the optimal solution of R (where the time-feasible requirement is entirely dropped), and then filter out those that are not time-feasible. If we only want to obtain one time-feasible optimal reconciliation or check if there is any, we can stop the enumeration immediately as soon as one such solution is found. However, if the number of optimal reconciliations is large and the time-feasible reconciliations are sparse (or nonexistent), this approach would require a long computation time before producing any output.

Based on the work of [Nøj+ ], we imagined a method for searching for time-feasible optimal reconciliations that utilizes the enumeration of EL-equivalence classes. Precisely, given an instance, it outputs a time-feasible optimal reconciliation if there is any, or "no" if there is none. In • Given an EL-equivalence class π (see Definition . ), is there a reconciliation φ in the class π that is time-feasible?

If the answer is yes, it also constructs such a reconciliation φ. Our new method is thus straightforward: first we enumerate the EL-equivalent classes, then, for each class, we apply the algorithm of [Nøj+ ] to output a time-feasible reconciliation if there is any. There are several pros and cons to this method:

• Pros:

-If there are no time-feasible reconciliations among the optimal ones, with this method we can expect to obtain the "no" answer faster than with the naive method (that is, by enumerating all optimal reconciliations and checking their time-feasibility).

-If the time-feasible reconciliations are sparse, we can expect to discover one such solution faster.

-Very often, we want to save the reconciliations to a file for subsequent analysis. In other words, we separate the enumeration and the time-feasibility check into two steps. The space needed for storing the EL-equivalence classes can be much smaller than the space required for storing all optimal reconciliations. The EL-equivalence classes file can be used for the task of finding time-feasible reconciliations as well as other tasks (such as visualization, computation of the average distance or a median solution).

• Cons:

-In an EL-equivalence class, there can be multiple time-feasible reconciliations. This method can output only one of them. It should not be used if we want to find all time-feasible optimal reconciliations.

-If the number of EL-equivalence classes is close to the total number of optimal reconciliations (we have seen in Table . that these two numbers can even be equal), there is little or no advantage in computation time or in storage space compared to the naive method.

As we have seen in Section . . , for the majority of biological datasets and cost vectors that we have chosen, the number of EL-equivalence classes is quite close to the number of optimal reconciliations, so there is little advantage in terms of the speed or space. Moreover, after enumerating the optimal reconciliations and computing the proportion of those that are time-feasible, we observed that this proportion is almost always high, therefore, we are rarely in a case where the time-feasible reconciliations are sparse or nonexistent. In fact, as we hinted in Chapter , Section . . , in practical situations, by applying the naive method, we are usually able to obtain time-feasible optimal reconciliations.

Nevertheless, there are still a few instances (for example, the dataset COG with the cost vector ( , , , )) in which the number of EL-equivalence classes is significantly smaller than the number of optimal reconciliations, and the latter is very large;

there are also a few instances (for example, the dataset COG with the cost vector ( , , , )) where none of the optimal reconciliations is time-feasible. Take the example of the dataset COG with the cost vector (-, , , ), there are optimal reconciliations and EL-equivalence classes (corresponding to a Reduction of . , see Section . . ). The space needed for storing the EL-equivalence classes is thus approximately . times the space needed for storing all optimal reconciliations, or times smaller. It turned out that none of the optimal reconciliations is time-feasible.

We compared the time needed for obtaining the "no" answer using the two methods: ( ) apply the O(n + m) time time-feasibility check of [THL ] to the list of optimal reconciliations, ( ) apply the (m log(n)) time algorithm of [Nøj+ ] to the list of ELequivalence classes. Roughly, the second method produces the answer times faster than the first one ( seconds versus seconds). In the future, it is possible that we come across more cases where finding time-feasible reconciliations by means of the EL-equivalence classes would be clearly superior to the naive method.

.

Capybara

. . Overview

Capybara [Wan+ ] (see https://capybara-doc.readthedocs.io/) is a cross-platform desktop application that help the user to perform various tasks around the R -. After loading an input file encoding a dataset (H, P, σ), the user can specify a cost vector that makes up the instance of interest, then select one of the following tasks:

• Counting or enumerating all optimal reconciliations.

• Applying the time-feasibility filter during the enumeration of reconciliations.

• Counting or enumerating the V-, E-, and CD-equivalence classes, or enumerating one representative in each equivalence class.

• Enumerating the Best-K reconciliations (see Section . . ).

• Converting a file containing a list of E-or CD-equivalence classes to a format readable by Capybara Viewer, a Web visualization tool that specifically handles these two types of equivalence classes (see Section . . ).

A subset of the features of Capybara is distributed as a Python package, which allows for customized manipulations of the equivalence classes during the enumeration (see Section . . ).

. . Best-K enumeration

Although at least one time-feasible reconciliations exists for any instance, it does not have necessarily have the minimum cost among all possible reconciliations. When we say that there is no time-feasible optimal reconciliation, it implies that all the timefeasible reconciliations of this instance are suboptimal, that is, they have a cost that is strictly larger than the minimum. A natural strategy of searching for time-feasible reconciliations consists, in this case, of enumerating suboptimal reconciliations and then applying the time-feasibility filter.

As early as , finding suboptimal solutions to an optimization problem has been studied for its connection with dynamic programming [BK ; Glu ; WB ]. The Best-K enumeration is an enumeration method that aims at finding a set of k solutions that are better than any solution not in this set (they can thus be optimal or suboptimal), and has been used extensively in the problem of finding k shortest paths in a graph [Yen ; KIM ; Epp ; HMS ; AL ]. Intuitively, for an instance of a minimization problem, we can consider a partial order on the set of feasible solutions such that one solution φ precedes another solution φ if and only if the cost of φ is strictly smaller than the cost of . The set of feasible solutions can then be viewed as an ordered list, where the order between solutions of the same cost are chosen arbitrarily.

A Best-K enumeration algorithm should output the first K elements of this list, or, if K is larger than the number of feasible solutions, all elements of this list.

Capybara offers the option of Best-K enumeration and its combination with the time-feasibility filter. The algorithm is based on a simple modification of the dynamic programming algorithm for enumerating all optimal reconciliations and was discussed in my Master thesis. It should be noticed that, whether the user chooses the optimal or the Best-K enumeration options, Capybara will check the time-feasibility criteria of [Sto+ ], which are more restrictive than the criteria of [THL ].

. . Visualization of E-and CD-equivalence classes

Given an E-equivalence class, we have the knowledge of the color for each parasite node bara Viewer; the user can either visualize one particular frame by selecting its index, or visualize all frames sequentially as an animation. Common interaction options such as zooming, panning, and scaling are also available.

For the CD-equivalence classes, we have the additional knowledge of the parasite/host associations, for each parasite node with an event type other than host-switch.

Using the same frames as in the visualization of E-equivalence classes, Capybara Viewer displays the label of the host node when the mouse hovers over a parasite node (if the parasite node is of event type host-switch, the host label is the symbol ?, in conformity with our Definition . ).

. . Python package

As we have seen in Section . , counting or enumerating equivalence classes can be a useful method for obtaining biological information from an instance. If we want to include the analysis of equivalence classes as one step of an experimental pipeline, Capybara's Python package can then be a preferred choice over the graphical version of Capybara. The Python package is called capybara-cophylogeny and can be found in the repository Python Package Index (PyPI). It is easy to install and to update the package to the newest version using pip, the standard Python package manager. The documentation page provides several examples illustrating how to effectively use the package for counting and enumerating equivalence classes in practical situations, such as performing the same task on a large number of instances using multiprocessing, or speeding up the computation using an alternative Python interpreter.

Generally, while the graphical version offers more features and is more suitable for users who are not technically oriented, the Python package offers more flexibility for running multiple tasks automatically and simultaneously.

.

Conclusion and perspectives

In Chapter , we introduced the R , defined equivalence relations on reconciliations, and proposed algorithms for counting and enumerating equivalence classes. In particular, we are able to enumerate, in polynomial delay, one representative reconciliation each equivalence class. In Chapter , we applied our algorithms to real datasets from the literature and analyzed their space of optimal reconciliations using equivalence classes. Even when the number of optimal reconciliations is huge (e.g., ), we managed to obtain useful information about the solution space through the list of equivalence classes. Finally, we presented Capybara, a user-friendly software which performs various computational tasks based on our algorithms and provides a tool for visualizing equivalence classes.

As a future direction, we would like to extend our algorithms to other biologically motivated definitions of equivalence relations. A future version of Capybara will include more visualization options. How to efficiently visualize several equivalence classes together in a compact and static fashion (i.e., without animation) is an especially interesting and challenging question.

The problem of efficiently enumerating equivalence classes is interesting in its own right. In Chapter , we will see a generalization of Algorithm to solve the E for general ad-AND/OR graphs; more perspectives will be presented there in the scope of enumeration algorithms.

Chapter

Equivalence classes of solutions and ad-AND/OR graphs

.

Introduction

An AND/OR graph is a well-known structure in the field of Logic and Artificial Intelligence (AI) that represents problem solving and problem decomposition (see, for example [MM ; Nil ]). In this chapter, we consider a particular flavor of AND/OR graphs: the acyclic decomposable AND/OR graphs or ad-AND/OR graphs, also known as the explicit AND/OR graphs for trees [DM ]. In Chapter , the same structure served as the basis for describing the reconciliation graph, a graph that represents the solution space of an instance of the R . More generally, ad-AND/OR graphs have an intimate relationship with problems that can be solved by dynamic programming (DP-problems).

As we suggested at the end of Chapter , enumerating equivalence classes is an .

Background and motivation

Whereas classical optimization problems require to find a single best solution, enumeration problems require to build an entire set of solutions. The study of enumeration problems offers many theoretical challenges as well as perspectives in practical applications. On one hand, enumeration problems can give rise to theoretical questions that are not applicable in the more traditional setting of decision or optimization problems, such as the output-sensitive complexity analysis (see Section . . ). For instance, the enumeration of minimal traversals of a hypergraph is a major problem that is not known to be solved in polynomial total time [FGS ]. On the other hand, enumeration problems often provide more satisfactory answers to real-world questions compared to optimization problems, as mathematical formulations (such as the sum of weights) may not adequately capture all the desired properties of the best solution.

This can be due to various reasons: the quality criteria may be complicated (for example, including the time-feasibility criterion in the R would yield a NP-hard problem); the choice may be determined by data that is not yet available (for example, finding a shortest path in a transportation network subject to real time disturbances); the decision may involve subjective factors that cannot be formally modeled. In these situations, an enumeration algorithm is able to provide a list of candidate solutions (optimal, or of high quality), so that the user can apply more complicated criteria, wait for data to become available or make full use of human expertise before making a choice.

While real-world applications motivate the formulation of enumeration problems, the practical usability of enumeration algorithms can be limited. The main reason for this is the large number of solutions. The set of candidate solutions can be too slow to be analyzed by the user or requires complex data analysis methods. One way of dealing with this problem is the generation of representative solutions; similar formulations include the generation of solutions that are far apart form each other, or the generation of a subset of solutions of high diversity. In theses kinds of approaches, we need to carefully define a new problem to better answer the original enumeration problem.

Regardless of the meaning of representative solutions, notice immediately that, for the method to be applicable, one either requires the set of candidate solutions to be small enough to be explored to a sufficient extend, or, when that is not the case (for example, when there are candidate solutions), one needs to find the representative solutions directly from the original enumeration problem instance.

In our work, we chose to find representative solutions through the enumeration of equivalence classes; this is by no means the only method of defining and obtaining representative solutions. Depending on the data types, it can be more straightforward to define an equivalence relation rather than a distance measure or a diversity measure on the set of solutions (or the other way around, for example, when the data consists of vectors of numerical values). Mathematically, a distance or diversity measure can be more costly or complicated to compute, as its computation typically involves multiple objects that need to be treated together, whereas the computation of equivalence classes is usually done by mapping each single object to a canonical object of the class to which it belongs.

Although there are similarities in the rationale between methods based on equivalence relations and methods based on distance measures, there are important differences that should not be ignored, especially when interpreting of the results. Along with a chosen distance or diversity measure, the user usually also needs to specify a number k as the size of the output (or an upper bound of the size), in such a way that the level of detail or the granularity of the representative solutions in the output can be explicitly controlled. The number of equivalence classes is directly determined by the definition of the equivalence relation, and consequently, the choices of equivalence relations are limited in practice to those that do not yield too much or too little equivalence classes. Another potential issue is that, in the enumeration output, equivalence classes may be presented without any good way of assessing their relative disparities (two solutions are simply non equivalent to each other, but cannot be "more or less" non equivalent). When the size of each equivalence class (or the fraction of the solution space that it represents) is of importance to the user, this information should also be included in the output to prevent any cognitive bias. Indeed, according to the applications, the user may be more interested in solutions that are the most representative, in the sense that they are similar to a large portion of the whole set of solutions, or she may be more interested in solutions that are outliers. Notice that, in data analysis methods based on distance measures, the preference between community detection and anomaly detection is usually directed encoded in the design of the algorithm. An algorithm that generates a list of equivalence classes does not necessarily need to consider the size of each equivalence class. Later in the chapter, we solve the problem of enumerating equivalence classes and the problem of finding the size of an equivalence class in our framework using two entirely different algorithms.

Computational problems that explore the idea of equivalence relation or equiva- .

Colored classes of solution subtrees of an ad-AND/OR graph . . Basic definitions

First, recall the following notations. If w is a node in a directed graph, its out-neighbors are called its children and the set of children is denoted by ch(w). The root node of a rooted tree T is denoted by r(T ). Next, let us restate the definitions of ad-AND/OR graphs given in Chapter .

Definition . (ad-AND/OR graph).

A directed graph G is an acyclic decomposable AND/OR graph (shortly, ad-AND/OR graph) if it satisfies the following:

• G is a acyclic (it is a DAG).

• G is bipartite: its node set V(G) can be partitioned into (A, O) so that all arcs of G are between these two sets. Nodes in A are called AND nodes; nodes in O are called OR + nodes.

• Every AND node has in-degree at least one and out-degree at least one. The set of nodes of out-degree zero is then a subset of O and is called the set of goal nodes;

the remaining OR + nodes are simply the OR nodes. The subset of OR nodes of indegree zero is the set of start nodes.

• G is decomposable: for any AND node, the sets of nodes that are reachable from each one of its child nodes are pairwise disjoint.

An example of an ad-AND/OR graph is given in Figure . .

Definition . (Solution subtree).

A solution subtree T of an ad-AND/OR graph G is a subgraph of G which: ( ) contains exactly one start node; ( ) for any OR node in T it contains one of its child nodes in G, and for any AND node in T it contains all its children in G.

The set of solution subtrees of G is denoted by T(G). It is immediate to see that a solution subtree is indeed a subtree of G: it is a rooted tree, the root of which is a start node. If we would drop the requirement of G being decomposable, the object defined in Definition . would not be guaranteed to be a tree (in general AND/OR graphs, the solutions are not called subtrees but instead subgraphs). In Figure . , one solution subtree is shown in bold. The following four-step procedure allows to obtain one solution subtree: ( ) start at any start node, ( ) for any visited OR node, visit one child, ( ) for any visited AND node, visit all children, ( ) stop when the goal nodes are reached. The set T(G) of solution subtrees of any ad-AND/OR graph G can be enumerated efficiently: after a preprocessing step in time linear in the size of G, the delay between outputting two consecutive solutions is linear in the size of the solution. This result is folklore and the algorithm is sometimes redesigned to fit the needs of a particular problem. Indeed, for a specific dynamic programming problem, the enumeration of the solutions can be done in an ad hoc manner, and the graph structure used by the algorithm is not explicitly called an ad-AND/OR graph. For instance, in [Don+ ], the authors devised an enumeration algorithm for the R that works not only on reconciliation graphs but on any ad-AND/OR graph. See [Rot ] for another example.

In the literature that explicitly studies general AND/OR graphs, the effort is naturally focused on the question of enumerating solution subgraphs in AND/OR graphs that are neither acyclic nor decomposable, which is much more challenging (see, e.g., [MM ; MB ; JT ; Gho+ ]).

We will define now a notion of subgraph for ad-AND/OR graphs that we will be used extensively. Intuitively, the relationship between such a subgraph and the entire graph correspond to the relationship between a DP-subproblem and the "full" optimization problem. To see this, let N(w) denote the number of solution subtrees of G/w for each OR + node w. If w is a goal node, then N(w) = . Otherwise, N(w)

Definition

= v i ∈ ch(w) w k ∈ ch(v i ) N(w k ).
In a nutshell, the number of solution subtrees of a subgraph starting from an OR and (respectively, an AND node) node is the sum (respectively, the product) of that of its children.

. . Color classes: definitions

Let G be an ad-AND/OR graph. Let ℂ be an ordered set of colors. The colors will be used to express the equivalence relation on the set of solution subtrees of G. Intuitively, two OR + nodes having the same color represent two alternative ways of solving the problem that can be considered equivalent.

Definition . (e-coloring

). An ad-AND/OR graph G is e-colored if its OR + nodes are colored in such a way that for any AND node all its children have distinct colors.

Notice here that only OR + nodes are colored and AND nodes are not colored. If w is a OR + node of G, we denote by c(w) its color. If w is an AND node, we denote by c(w)

the tuple of colors of the children of w sorted in increasing order of the colors. From the definition of an e-coloring, the color tuple c(w) of an AND node w is necessarily of the same size as the set of children ch(w). The intuition is that the color tuple c(w)

distinguishes each child OR + node of an AND node w, and gives a specific ordering between them.

If T ∈ T(G) is a solution subtree of G, we use the notation π(T ) for the result of contracting the AND nodes in T : for each OR node w of T , contract the only child node of w in T (i.e., remove the child and connect w to each one of its "grandchildren").

Definition . (Color class

. . Color classes: enumeration

Given an ad-AND/OR graph G that is e-colored with a fixed ordered color set ℂ, we propose a polynomial delay algorithm to enumerate all equivalence classes of G.

Total ordering over the set of color classes

We define a total ordering ≺ (called smaller than) over C(G), the set of color classes of (C ′ , . . . , C ′ ℓ )) be the child subtrees of r(C) (resp. of r(C ′ )) sorted recursively with respect to ≺. We then say that C is smaller than

C ′ if the tuple (C , . . . , C k ) is lexicographically smaller than (C ′ , . . . , C ′ ℓ ), i.e., if C i ≺ C ′ i with i being the smallest index such that C i ≠ C ′ i .
We also assume that ∅ is smaller than any tree, and therefore a single-node tree colored with a color c ∈ ℂ comes before any other tree whose root is colored with c in ≺.

Definitions and notations

Given an OR node w, we denote by ct (w) the set of color tuples of its children, i. Given a set O of OR + nodes, we consider the set ct (O in lexicographical order. For an index r ≤ |ct (O)|, let t r be the r-th color tuple of ct (O). We denote by P r (O) the subset of child AND nodes of O whose color tuple is t r , i.e.,

P r (O) := {v ∈ ch(O) : c(v) = t r } .
When the set O is clear from the context, it is omitted from the notation and we simply write P r Algorithm , we will only work with a fixed set O). The subsets P , . . . , P |ct (O| form a partition of ch(O).

Finally, given a color tuple t r := (c , . . . , c j ) ∈ ct (O), for each i ≤ j, by Definition . , each AND node in P r has exactly one child OR + node colored with c i . We denote by P r i the set of child OR + nodes of P r colored with c i , i.e.,

P r i = {w ∈ ch(v) : v ∈ P r , c(w) = c i } .
In summary, the key notation P r i that we will use in Algorithm can be described as follows: it is the set of grandchild OR + nodes of O colored with the i-th color c r of the r-th color tuple t r ∈ ct (O).

In the left panel of Figure . , an example graph is shown where each node is labeled by an integer. The colors are, in increasing order, w, x, y, and z. For all i ≤ j, let O i ⊆ P r i be the set of nodes in P r i compatible with (C , . . . , C i-), the first (i -) subtrees of r(C)

Let ℓ be the largest index i ≤ j such that Next(c i , O i , C i ) ≠ ⊥ if such index exists. Otherwise, C ← ∅ and go to line end C ℓ ← Next(c ℓ , O ℓ , C ℓ ) for ℓ < i ≤ j do Let O i ⊆ P r
i be the set of nodes in P r i compatible with (C , . . . , C i-)

C i ← Next(c i , O i , ∅)
end Return A tree with root color c and root child subtrees (C , . . . , C j )

In this second case, we turn our attention to the child subtrees (C , . . . , C j ) of r(C), ordered according to their root colors. Observe that, since C is a color class of G/O, for each i ≤ j, the i-th child subtree C i is a color class of G/P r i . It then suffices to replace the last child subtree C j by its successor C ′ j if there exists one, and we obtain a solution C ′ whose subtrees at the root are (C , . . . , C j-, C ′ j ), which is by definition of the successor of C. If C j has no successor, we replace if possible C j-by its successor C ′ jand we replace Cj by the smallest admissible color class (i.e., the successor of ∅). In general, we select at each step the greatest index ℓ such that C ℓ has a successor with respect to ≺ (Line ), and we replace it by its successor C ′ ℓ (Line ), and we take the smallest admissible color class for every ℓ < i ≤ j (Lines -).

Without further care, the above described procedure would output color classes C ′ whose child subtrees (C ′ , . . . , C ′ j ) at the root correspond to the elements of the Cartesian product of C(G/P r i ), for i ≤ j. However, while it is true that if C is a color class of G/O, its child subtree C i is a color class of G/P r i for all i ≤ j, the converse is not true. Indeed, not all elements of C(G/P r ) × • • • × C(G/P r j ) lead to an admissible color class (an example is given in the right panel of Figure . ). In order to find an admissible color class, we should guarantee that the choice of a given C i is compatible with the previous choices (C , . . . , C i-). This is done by selecting the subset of OR + nodes O i ⊆ P r i that are compatible with (C , . . . , C i-) (Lines and ). An admissible choice of C i will then be any color class of G/O i . The definition of compatible nodes will be given later. There are two important remarks with respect to the addition of this compatibility check in the algorithm. First, for each i ≤ j, the subset O i can be easily computed (this can be seen from Definition . ). Secondly, O i is never empty, i.e., there is always a choice for C i that is compatible with the previous choices of (C , . . . , C i-) (there is at least one choice that corresponds to the current color class). Notice that if this property were not satisfied, the algorithm would not have a polynomial delay complexity since we may spend exponential time without reaching any color class. With this property, we are guaranteed to be always able to extend a partial tuple of subtrees (C , . . . , C i ) until we reach a complete tuple (C , . . . , C j ) that will form an admissible color class.

Compatible nodes

Given Definition . . Let O be a set of OR + nodes of color c. Let t r =: (c , . . . , c j ) ∈ ct (O) be its r-th color tuple. Let k be an integer such that ≤ k < j. For i ≤ k, let C i be a color class in C(G/P r i ). We say that a OR + node w ∈ P r k+ is compatible with the k-tuple (C , . . . , C k ) if there exists an AND-node v ∈ P r that satisfies the following: ( ) w is a child of v, ( ) for

all i ≤ k, the set r(π -(C i )) contains a child of v.
In Algorithm , we extended Definition . to the case where k = , and any node in P r is compatible with the empty tuple.

The next lemma, essential to the proof of correctness, states the key property of compatible nodes: choosing a compatible node guarantees that a partial tuple of subtrees can be "safely" extended, without the risk of creating any non-admissible color class, that is, a color class that does not correspond to any solution subtree of G.

Lemma . . Let O be a set of OR + nodes G all colored with c. Let t r =: (c , . . . , c j ) be the r-th color tuple in ct (O). Let C be a color class of G/O, and let C , . . . , C k , be its first k child subtrees, for some fixed k < j, where C i ∈ C(G/P r i ) for all i ≤ k. Let O k+ ⊆ P r k+ be the set of OR + nodes compatible with (C , . . . , C k ). Given C k+ ∈ C G/P r k+ , there exists a color class

C ∈ C(G/O) whose first (k + ) child subtrees at the root are precisely (C , . . . , C k , C k+ ) if and only if C k+ ∈ C(G/O k+ ).
Proof. (First direction) Assume that there exists a color class C of G/O whose first (k+ ) child subtrees at the root are (C , . . . , C k , C k+ ), and let T be a solution subtree of G/O such that π(T ) = C. Let v be the unique child AND node of the r(T ). Notice that since we have C i ∈ P r i for all i ≤ k, the color tuple c(v) of v is necessarily the r-th color tuple t r ∈ ct (O), and so v ∈ P r . Let (w , . . . , w j ) be the child OR + of v in T , ordered in increasing order of their color. For all i ≤ k we have

w i ∈ r(π -(C i ), the node w k+ is thus compatible with (C , . . . , C k ). Therefore, C k+ ∈ C (G/O k+ ) since C k+ ∈ C(G/w k+ ), and
w k+ ∈ O k+ . (Second direction) Assume now that C k+ ∈ C (G/O k+ ).
Then there exists a compatible OR + node w k+ ∈ O k+ and a solution subtree T k+ of G/w k+ rooted at w k+ that satisfies π(T k+ ) = C k+ . Since w k+ ∈ O k+ , by the Definition . of compatible nodes, there exists an AND node v ∈ P r such that w k+ is a child of v, and, for all i ≤ k, there exists an OR + node w i ∈ r(π -(C i )) such that c(w i ) = c i (i.e., w i is the unique child of v colored with the color c i ). Therefore, for all i ≤ k, there exists a solution subtree T i of for all i ≤ j, C ′ i is the smallest color class of G/O i , where O i ⊆ P i is the set of nodes of P i compatible with C ′ , . . . , C ′ i-. Thus, Next(c, O, ∅) will return C ′ at Line , since r will receive at Line , C will receive Next(c , P , ∅) at Line which is equal to C ′ by the induction hypothesis, and for all i ≤ j, C i will receive Next(c i , O i , ∅) at Line which is equal to C ′ i by the induction hypothesis.

G/w i such that π(T i ) = C i .
Proof of Lemma . . Let (C ′ , . . . , C ′ j ) be the child subtrees of r(C) and let t r := (c , . . . , c j ) be the root color tuple of C, i.e., the color tuple with which the roots of the subtrees (C ′ , . . . , C ′ j ) are colored. Let O i ⊆ P r i be the set of nodes of P r i compatible with (C ′ , . . . , C ′ i-) for all i ≤ j. Let C ′ be the color class of G/O that follows C with respect to ≺. Let (C ′′ , . . . , C ′′ j ) be the child subtrees of r(C ′ ). Notice that the children of r(C ′ ) are either colored with the same color tuple t r , or colored with t r+ , the next color tuple of

ct (O), if C is the largest color class of G/O with the root color tuple t r .
Assume first that the children of r(C ′ ) are colored with the same color tuple t r . Let ℓ be the smallest index such that C ′ ℓ ≠ C ′′ ℓ . We claim that ℓ is also the largest index such that C ′ ℓ has a successor in C(G/O ℓ ) with respect to ≺, and thus it corresponds to the ℓ chosen by the algorithm at Line . Indeed, assume that there exist i < k ≤ j and a larger color class

C k ∈ C(G/O k ) such that C ′ k ≺ C k . By Lemma .
, there exists a color class of G/O whose first k child subtrees at the root are (C ′ , . . . , C ′ k-, C k ). However, in this case such a color class would be greater than C and smaller than C ′ with respect to ≺, and it would be in contradiction with the fact that C ′ immediately follows C in ≺.

Now, C ℓ will receive Next(c ℓ , O ℓ , C ′ ℓ )
at Line , which is by the induction hypothesis the color class that succeeds C ′ ℓ in C(G/O ℓ ). Since we assumed that ℓ is the smallest index such that C ′ ℓ ≠ C ′′ ℓ , by definition of ≺ and by Lemma . , C ′′ ℓ is the color class that succeeds C ′ ℓ in ℂ(G/O)ℓ), and so C ℓ will receive C ′′ ℓ at Line . Since for all i ≤ ℓ we have

C ′ i = C ′′ i
, and since (C , . . . , C ℓ-) are not modified by the algorithm, at the end of it, (C , . . . , C ℓ-, C ℓ ) will be equal to

(C ′ , . . . , C ′ ℓ-, C ′′ ℓ ) = (C ′′ , . . . , C ′′ ℓ-, C ′′ ℓ )
. It now remains to show that C i will be equal to C ′′ i for all ℓ < i ≤ j. Again, by Lemma . , for all ℓ < i ≤ j,

C ′′ i is the smallest color class of G/O ′ i ),
where O ′ i is the set of nodes of P r i compatible with (C ′′ , . . . , C ′′ i-), since otherwise, another color class of G/O greater than C and smaller than C ′ could be built. Thus, at Line , C i will receive Next(c i , O ′ i , ∅), which is equal to C ′′ i by the induction hypothesis, and C ′ will be returned at Line . Assume now that the children of r(C ′ ) are colored with the next color tuple t r+ . In this case, C is the greatest color class of C(G/O) with the root color tuple t r . We claim that each child subtree C ′ i of r(C) is also the greatest color class of G/O i for all i ≤ j. Indeed, assume otherwise that, for some i, there exists a color class

C i ∈ O(G/O i ) with C ′ i ≺ C i .
Then, by Lemma . , there exists a color class C ∈ C(G/O) with the same root color tuple t r that is greater than C, contradicting the fact that C is the greatest color class with the root color tuple t r . By the induction hypotheses, Next(c i , O i , C ′ i ) will therefore return ⊥ for all i ≤ j. So C will receive ∅ at Line and the algorithm will go to Line , and the next color tuple t r+ will be selected at Line . Using now similar to the ones used in the proof of Lemma . , the smallest color class with the root color tuple t r+ will be returned. For the complexity, notice that at most one recursive call is performed by each node of the next color class. More precisely, if Next(c, O, C) returns C ′ , there will be exactly one recursive call per node in C ′ that is not in C, and thus at most s recursive calls will be performed.

Proof of

In each recursive call, the set of color tuples ct (O) and the sets P r i , where r ≤ |ct (O| and i ≤ j := |t r |, can all be computed in O(n) time in total. It remains to show that the sets of compatible nodes O i , i ≤ j, can be computed in O(n) time in total which will conclude the proof. To do this, we should be able to compute the sets r(π

-(C i )) for all i ≤ j. If Next(c, O, C) = C ′ ,

. . Restricting the graph to a color class

After the enumeration of the color classes, it might be interesting to go back to the solution subtrees that belong to each color class. In particular, in practical applications, one might want to use the number of solution subtrees as a measure for the "importance" or "significance" of a color class. In Algorithm , we present an algorithm that,

given an e-colored ad-AND/OR graph G and a color class C ∈ C(G), constructs the subgraph G C of G restricted to C, that is, a subgraph of G of which the solution subtrees are exactly the ones of belonging to the color class C:

T(G C ) = π - G (C) := {T ∈ T(G) : π(T ) = C} .
Once the graph G C is obtained, the following questions can be answered (by applying the same method as for any ad-AND/OR graph): counting the number of, and enumerating the solution subtrees belonging to the color class C.

Algorithm relies on two recursive functions VisitOR + and VisitAND, both taking as input a node in G and a node in C. The Require statements are used to specify the preconditions that the two parameters of the two Visit functions must verify; it can be checked by inspection that these conditions are always satisfied whenever the functions are called. The algorithm performs an operation called Mark on the nodes in G.

All nodes are initially unmarked; the Mark operation changes the state of a node into marked.

Recall that π(T ) transforms a solution subtree

T ∈ T(G) into a color class C ∈ C(G)
by contracting the AND nodes in T . For a fixed T , we extend this notation and write π(w) = π(v) := x for every OR node w in T with its unique AND-child node v in T that are identified with the node x in C under the transformation.

Lemma . . Let T be a solution subtree of G belonging to the color class C. For every node w of T , there is a call to either the function VisitOR + or to VisitAND of Algorithm with parameters w and π(w), depending on whether w is an OR + node or an AND node.

Proof. By top-down induction. The start node of T is visited at Line since it has the correct color. In the induction step we separate two cases. For an OR + node w of T that is not a start node, suppose that the parent v of w in T is visited in a call VisitAND(v , π(v )). Then w is visited (Line ), and the second parameter is π(w). In the other case, for an AND node v of T , suppose that the parent w of v in T is visited in a call

VisitOR + (w , π(w )). Since we have π(v) = π(v ) and c(π(w )) = c(v), the condition at Line is satisfied and v is visited in a call VisitAND(v, π(v)).
The correctness of Algorithm is shown in Theorem . . We omit the analysis of complexity as the algorithm clearly requires a running time that is linear in the size of the input graph G.

Require: c(w) = c(x)
if w is a goal node and x is a leaf then

Mark(w) return end for each child AND node v i of w in G such that c(v i ) = c(x) do VisitAND (v i , x) end if at least one child of w is marked then Mark(w) end Function VisitAND(v, x): Require: c(v) = c(x)
for each child OR + node w i of v do

x i ← the unique child of x such that c(w i ) = c(x i ) VisitOR + (w i , x i ) end if all children of v are marked then Mark(v) end Proof. (First direction)
We show that any solution subtree of G C is also a solution subtree of C, and that it belongs to the color class C. For every marked OR node, at least one child is marked (Line ); for every marked AND node, all its children are marked (Line ). A solution subtree of G C is thus a solution subtree of G. Let T be a solution subtree of G C , consider the recursion tree of the Visit function calls during which the nodes in T are marked. By the preconditions of the Visit functions, that is, c(w) = c(x)

and c(v) = C(x), the tree π(C) is equal to the tree formed by the colored nodes that are used as the second parameter x in this recursion tree. The latter is simply equal to C (we start with the root of C, then visit each child of the current node), so we have π(T ) = C.

Therefore, every solution subtree of G C belongs to the color class C.

(Second direction) Let T be a solution subtree of G such that π(T ) = C, we show that every node in T is marked by bottom-up induction. By Lemma . , any goal node w in T is visited in a call VisitOR + (w, π(w)) so w is marked (Line ) because π(w) is necessarily a leaf. For the induction step we separate two cases. Let w be an OR node in T , and suppose that all nodes in the subtree T | w are already marked. Again, by Lemma . , w is visited. Then w is marked at Line , since exactly one child of w is in T and is thus marked. In the other case, let v be an AND node in T and suppose that all nodes in T | v are marked. By the lemma, v is visited. Then v is marked at Line , because all children of v are in T and are thus marked. This completes the proof.

. In this model, a tree-sequential problem can be specified by a finite (bottom-up) tree automaton A = (Q, Σ, δ, q , Q F ), where Q is a finite set of states, Σ is a ranked alphabet, δ is a set of transition rules of the form (q , . . . , q n , a, q) where q , . . . , q n , q ∈ Q and a ∈ Σ, q ∈ Q is the initial state, Q F ⊆ Q is a set of final states. The problem specification also includes a cost function. The set L(A) of trees accepted by the tree automaton A defines the set of feasible solutions. The minimization problem seeks to minimize the cost function over the set L(A) of feasible solutions.

Application to dynamic programming

We will consider the simple case of a positive additive cost function that always equals zero in the initial state. An additive cost function can be defined via an incremental cost function I : Q * × Σ → ℝ, where Q * consists of tuples of states in Q of the form (q , . . . , q n ). I(q , . . . , q n , a) can be viewed as the cost of attaching n child subtrees In this case, the answer of the minimization problem can be shown to be equal to min q∈Q F D(q), where D : Q → ℝ ≥ is defined by the following recurrence equations:

D(q ) = ,
for q ≠ q , D(q) = min (q ,...,q n ,a,q)∈δ ≤i≤n D(q i ) + I(q , . . . , q n , a) .

( . )

A dynamic programming algorithm for the minimization problem corresponds to an algorithm that computes D; the function D is commonly called a dynamic programming table (a DP-tabled, also called a DP-array, or a DP-matrix). In general, to be able to write down the recurrence relations does not imply that there exists an efficient algo-rithm to compute the function. Indeed, such an algorithm does not exist in general for given arbitrary tree automata and cost functions [Iba ].

Using an algebraic approach, Gnesi and Montanari [GMM ] have shown that solving the functional Equation ( . ) corresponds to finding the solution subtrees of a general AND/OR graph. An important special case in which DP-algorithms exist is when the underlying AND/OR graph is acyclic.

When a fixed tree is given as an input to the problem, the underlying AND/OR graph is acyclic and decomposable (that is, it is an ad-AND/OR graph

Explicit construction of the ad-AND/OR graph for DP on a fixed tree

Due to its usefulness for the examples that we will develop next, in the case of DP on a fixed tree, an explicit construction of the ad-AND/OR graph from Equation ( . ) is described below. The construction is done in two steps. In the first step, we build a graph in which every node retains an additional attribute, its value, and every OR + node is labeled by a state q ∈ Q. In the second step, we prune the graph by removing nodes that do not yield optimal values.

. For each (q , a, q) ∈ δ, create a goal node of value labeled by q. Then, for each q ≠ q in post-order, i. For each (q , . . . , q n , a, q) ∈ δ, create an AND node, connect it to the n OR + nodes labeled by q , . . . , q n . Its value is equal to the sum of the values of its children, plus I(q , . . . , q n , a).

ii. Create a single OR node, connect it to every AND node created in the previous step. Its label is q, and its value is the minimum of the values of its children.

. For each q ∈ Q F , remove the OR node labeled by q unless its value is equal to min q∈Q F G(q). For each OR node s, remove the arc to its AND-child node s i if the value of s i is not equal to the value of s. Finally, remove recursively all AND nodes without incoming arcs.

A similar procedure was described in Chapter , Section . . for constructing the reconciliation graph, which is an ad-AND/OR graph with some additional properties (one such property is that every AND node has exactly two ordered children).

. . Examples

We will give several examples of DP-problems for which the enumeration of equivalence classes of solutions can be done by enumerating the color classes of the underlying ad-AND/OR graph.

Optimal tree coloring problem

Description A prototypical problem that fits into the framework of DP on a fixed tree is O , that is, finding an optimal node-coloring of the input tree. Many problems of practical interest can reduce to O ;

three concrete examples are given later in this section.

If T is the input (rooted, ordered) tree and ℂ is the set of colors, such a problem seeks a coloring φ : V(T ) → ℂ that minimizes the cost function. There can be many constraints on the coloring function: some nodes of T may be forced to have a certain color, the possible colors of a node may depend on the colors of its descendants. In our tree-sequential dynamic programming formalism, a tree automaton and a cost function are given as part of the input. The tree automaton defines the set L(A) of feasible coloring functions satisfying all those constraints. A state q can be interpreted as a colored subtree of T with a particular root color; the unique initial state is an empty coloring and transitions into a colored leaf of T ; a final state corresponds to a fully colored T with a particular root color. A commonly used form of cost functions considers the (possibly weighted) sum over the edges of the tree of the cost of putting two colors on each end of an edge, that is, an incremental cost function I of the form I(q , . . . , q n , a, p) = ≤i≤n p(a i , a) where a i is the color of the root of the subtree in state q i and p : C → ℝ ≥ is a function that gives the cost of putting two colors at each end of an edge.

Equivalence relations on the set of solutions A possible strategy to define equivalence classes on the solution space of the O problem is to consider some colors to be locally equivalent on a node. In practical applications, the space of colors can be quite large (this is the case for the R from Chapter ). Even though the precise colors of each node are necessary for correctly computing the cost function, when the solutions are analyzed by a human expert, it can be desirable to omit the colors and just look at whether the color of a node belongs to some group of colors. Therefore, this kind of equivalence relations is natural in many situations. Our Definition . of color classes of an e-colored ad-AND/OR graph deals exactly with equivalence relations of this type.

Let e be a function that maps a node u and a color c to the "color group" of that node, denoted by e(u, c). Two solutions of the O problem φ , φ : V(T ) →

ℂ are said to be equivalent if ∀u ∈ V(T ), e(u, φ (u)) = e(u, φ (u)). Let G be the ad-AND/OR graph associated with this instance. Then G can be e-colored as follows. For each OR + node w of G labeled with the state q, where q is interpreted as a colored subtree of T | u of T rooted at a node u colored with c, we color the node w with e(u, c). After this, G is e-colored, and its set of color classes C(G) corresponds to the set of equivalence classes of the solutions of the instance. Notice that the constraint we had on the e-coloring of an ad-AND/OR graph is naturally satisfied by any meaningful function e because in a DP setting we only consider ordered trees: the i-th and j-th children of a node of T cannot be in the same color group unless i = j.

Concrete examples of tree coloring problems

In each of the following problems, the connection with the O problem is straightforward. The aim here is to show that our framework allows to consider equivalence relations that are natural or appropriate for those problems.

Example

In the R that we studied in Chapter , the input tree is the parasite tree, the set of colors is the set of nodes of the host tree, and the notion of reconciliation defines the set of feasible colorings. The E-, EL-, and CDequivalence relations, all of which can be shown to yield an e-coloring of the reconciliation graph, are motivated by biological considerations and by practical needs. In Chapter , we have seen that useful information can be extracted from the equivalence classes to help us better understand the solution space, which can be huge for some datasets.

Example

This example is related to the alignment of gene sequences on a phylogenetic tree, known as the T [San ]. In this problem, one seeks to infer the ancestral sequences, knowing the sequence on each taxa, in a way that is most parsimonious (i.e., minimizing the edit distance between adjacent sequences).

The general version of the problem considers sequences of arbitrary length and is NP- The goal is to find a coloring φ : V(G) → C such that a ∀v ∈ V(G), φ(v) ∈ S(v), and

∀(u, v) ∈ E(G), (φ(u), φ(v)) ∉ T .
While this problem is hard in general, it can be solved by DP when the underlying graph is a tree. In this case, we can enumerate colorings of the input trees without any forbidden pair of colors on the edges. Defining equivalence relations by grouping some of the colors together (for example, frequencies that could be considered similar), and enumerating only one solution per equivalence class can be a practical way of reducing the size of the output.

Graph optimization problems and DP on tree decomposition

Graph optimization problems can model a large variety of real-world problems. Many hard graph problems can be solved in polynomial time with a dynamic programming algorithm when the input graph has bounded treewidth (see, for example, [Bod ]).

The underlying idea is that, given a tree decomposition of a graph, the DP-algorithm traverses the nodes (bags) of the decomposition and consecutively solves the respective sub-problems.

A well-studied type of graph problems is the vertex subset optimization problems, which ask for a subset of vertices of the input graph that is optimal with respect to some graph-theoretical properties. Given a bag X, a DP-algorithm generally computes for each subset Z ⊆ X the optimal solution of the DP-subproblem whose intersection with X is Z. In this context, we could define two solutions to be equivalent if they intersect each bag of the decomposition in an "equivalent" way. The equivalence relation on the solutions is then defined by an equivalence relation over the subsets of each bag, and two solutions S and S are equivalent if for all bags X of the decomposition, S ∩ X is equivalent to S ∩ X.

One of the simplest examples is to consider that all the nonempty subsets of ver-tices of a bag are equivalent. Thus, what we are interested in is whether a solution "hits" a bag (i.e., whether it has a nonempty intersection with the vertices in the bag).

Consequently, two solutions would be considered equivalent if they hit the same bags.

We can also consider two subsets of a bag to be equivalent if they have the same size.

In this case, two solutions would be equivalent if each bag contains the same number of vertices in the two solutions.

The above two examples give an idea on how equivalence relations which fit into our framework can be defined. We believe that considering solutions in the way they are distributed along the tree decomposition of a graph could give a good overview of the diversity of the solution space.

Conclusion and perspectives

In this chapter, we provided a framework for the enumeration of solutions in polynomial delay that can be applied to a whole family of problems, and we showed several examples of applications in dynamic programming problems. A few key ingredients are required: one is that the solutions must be represented in an AND/OR graph that is In these applications, one typically wishes to understand the data, which is the graph itself, by exhibiting its patterns, represented as subgraphs satisfying some desired properties (e.g., connectivity, density). For the purpose of discovering distinctive patterns, an equivalence relation can be introduced, and the new question is then how to efficiently enumerate subgraphs that satisfy the desired properties and that are not equivalent to each other.

In the last part of Section . . , we hinted that our approach can be used in the context of graph optimization problems, together with the tree decomposition technique and DP-algorithms on tree decomposition. For a specific graph problem such as D -

, one might want to define an equivalence relation that is more tailored to the problem and not based on tree decomposition. More generally, the most suitable equivalence relation on the set of solutions of a problem could be based on some non-local properties, and our framework would not apply. A recursive algorithm such as Algorithm works well when the property of interest is local in some sense (here, in the form of an e-coloring). In order to efficiently enumerate equivalence classes that are based on non-local properties, we have to develop new techniques. Indeed, in Chapter , the enumeration of V-equivalence classes, which is based on the non-local property of event vector, was done using an entirely different method.

Chapter

Maximum agreement forest

.

Introduction

A natural question that arises in phylogenetic analysis is the quantification of dissimilarity between different (i.e., non isomorphic, see Chapter , Definition . ) phylogenetic trees on the same set of taxa. The difference in the tree topology can be a consequence of different inference methods, or different data sources being used (e.g., morphological data, behavioral data, genetic data, etc). Furthermore, regardless of the inference method or the data source, the constructed tree may not correctly represent evolutionary relationships, because not all groups of species follow a simple tree-like evolutionary pattern. The non tree-like evolutionary processes, such as hybridiza- In this chapter, we define formally the MAF and present a heuristic for it, that is, an algorithm that provides an approximate solution without a guaranteed approximation ratio, then we compare this algorithm experimentally with the existing approximation algorithms and show that it outperforms all but one of them in terms of the quality of the approximation. We conjectured that this algorithm also achieves the approximation ratio of . If this were true, we would have yet another -approximation, again based on entirely different ideas from the two existing ones, that is easier to understand and to implement and performs better in practice.

.

Definitions

Basic definitions about trees and forests can be found in Chapter , Section . . . In this chapter, the central object that we study is a X-forest, or a leaf-labeled rooted binary forest. It is a graph whose connected components are rooted binary trees, and its leaves are bijectively labeled by a set X. To be more consistent with the terminology commonly found in the literature (where the vocabulary associated with undirected trees and forests is most often used), in this chapter, the set of arcs of a forest F will be denoted by E(F); nodes and arcs will be called respectively vertices and edges. Recall from Chapter that if v is a vertex in a forest, p(v) denotes its parent.

The edge incident on a leaf is called the pendant edge of that leaf. An isolated vertex in a rooted binary forest is a leaf that is also the root of a component (i.e., an out-degree zero and in-degree zero vertex). A pendant subtree of an X-forest is a subtree that can be Definition . (Leaf pair). Let (F , F ) be a pair of X-forests, i.e., two rooted binary forests whose leaves are labeled by the same set X. A leaf pair (v , v ) of (F , F ) is an ordered pair consisting of a leaf v of F and a leaf v of F having the same label.

Definition . (Cutting an edge from an X-forest). Let F be an X-forest. Let E ⊆ E(F) be a set of edges of F. Then C (F, E) is an X-forest obtained from F by applying the following operations:

. Remove from E(F) the edges in E.

. Repeatedly contract each out-degree one and in-degree one vertex, i.e., remove it and then add an edge between its parent and its unique child.

. Remove all out-degree one and in-degree zero vertices (their unique child then becomes the new root of that component).

When E contains only one edge e, C (F, E) can be simply written as C (F, e). The operation of cutting an edge is also called making a cut.

Definition . (AF, MAF

). An agreement forest (AF) of two X-forests F and F is an Xforest F such that there exist E ⊆ E(F ) and E ⊆ E(F ) satisfying

F C (F , E ) C (F , E ) .
By choosing E = E(F ) and E = E(F ), we can always obtain an AF of F and F whose components are all isolated vertices.

A maximum agreement forest (MAF) is an AF whose order is minimum.

For a pair of X-trees (T , T ) (i.e., two X-forests of order ), m(T , T ) denotes the order of a MAF minus one, that is, m(T , T ) := min{|F| -: F is an AF of T and T } . ( . )

Since each cut increases the number of components by one, the number m(T , T ) can also be interpreted as the number of edges that one needs to cut in each tree to produce an agreement forest.

The MAF can be stated as follows: given two X-trees T and T , find the value of m(T , T ).

It is known that m(T , T ) is equal to the rSPR distance if the following preprocessing operations have been applied to each of the two input X-trees: ( ) add a dummy leaf labeled by ρ (a symbol that is not already in the label set X), ( ) add ρ to the label set X, and ( ) add a dummy root whose two children are the dummy leaf ρ and the original root of the tree [BS ]. As our interest in solving the MAF originates from the computation of the rSPR distance, we will assume that this preprocessing step has already been applied. For convenience, we will refer to the optimal value of an instance of the MAF simply as the distance. Therefore, the value of m(T , T ) defined in Equation ( .) is called the distance between and T .

Even though an instance of the MAF is a pair of trees, for the purpose of analyzing our algorithm, it is worthwhile to extend the notion of distance to a pair of forests.

Definition . (Distance between a pair of X-forests). Let (F , F ) be a pair of X-forests. The following easy lemmas are useful for expressing the quality of a cut in terms of the distance. As the distance can be intuitively understood as the number of cuts that still need to be made before reaching a MAF, a cut in one forest decreases (the corresponding element of) the distance by either zero or one; a cut is effective if it indeed decreases the distance by one.

Let

Lemma . (Making cuts does not increase the distance). Let (F , F ) be a pair of Xforests. Let E ⊆ E(F ) and E ⊆ E(F ) be some subsets of edges of F and F . Then .

The algorithm

. . Reduction rules

Given a pair of X-trees T and T , we seek an approximate solution of their distance, that is, an upper bound of m(T , T ). The general strategy of the algorithm is straightforward. At any step, we have a pair of X-forests and we cut edges in both of them until the pair becomes isomorphic; the result is an AF and the final number of cuts is an upper bound of the distance. Notice that not all algorithms in the literature follow this approach, for example, the -approximation of [WZ ] only cuts edges in one of the two input trees.

The first important component of the algorithm is the reduction rules. These rules tell us how to transform a pair of forests into another one with strictly fewer leaves while preserving the distance. Intuitively, these rules correspond to the observations of the type "it is forbidden (or mandatory) to cut such an edge". those structures can be replaced by simpler ones. The fourth rule is of a different nature and identifies an edge that must be cut in exactly one of the two forests.

• Reduction rule (Subtree reduction): Replace any subtree of more than one leaf that occurs in both forests by a single leaf with a new label.

• Reduction rule (Chain reduction): Replace any chain of subtrees that occurs identically and with the same relative orientation in both forests by three new leaves with new labels correctly oriented to preserve the direction of the chain.

• Reduction rule (Root-chain reduction): Replace any chain of subtrees that occurs identically and with the same relative orientation in both forests and starts from the root of its component by a single leaf with a new label.

•

(C (F , e ), F )) = (d -, d ) (resp. m(F , C (F , e ))) = (d , d -)) where (d , d ) = m(F , F ).
Proof. The proof is trivial: if a leaf appears as an isolated vertex in either F or F , it must appear as an isolated vertex in every MAF of F and F .

During the algorithm, the reduction rules will be applied at the very beginning as well as at the end of each iteration. The remaining steps of an iteration thus only consider a pair of forests that is irreducible.

Definition . (Irreducible pair of forests)

. A pair of X-forests is called irreducible if none of the reduction rules -can be applied.

The next lemma says that a pair of forests can be made irreducible efficiently, and the size of an irreducible pair of forests is bounded linearly by their distance. The proof is omitted and follows from the similar statements in [AS ; BS ] (the addition of the two new reduction rules does not affect the analysis).

Lemma . . A pair of X-forests (F , F ) of size n (i.e., the number of vertices in each forest) can be transformed into an irreducible pair (F ′ , F ′ ) of size n ′ in O(n) time, and the size n ′ is bounded linearly by m(F , F ).

. .

Common continuous subtree (CCS) Definition

The central idea of the algorithm is to identify subtrees that are "almost" identical in the pair of forests. This is in a similar spirit as the first reduction rule, which identifies subtrees that are exactly identical. In our definition of "almost identical subtrees", we allow at most one cut in either forest, that is to say, after cutting at most one edge, that subtree will occur in both forests and the reduction rule applies. The next few definitions express this formally.

Definition . (CS).

Let F be a X-forest. A continuous subtree (CS) of F is a subtree H of F having at least two leaves, and in which at most one edge (u, v) is marked as the cut edge and satisfies u ≠ r(H) (i.e., the cut edge is not incident to the root of the CS).

When no edge is marked as the cut edge, the cut edge of the CS is defined to be the empty set ∅. Accordingly, we define C (F, ∅) := F.

Definition . (CCS)

. Let (F , F ) be a pair of X-forests. A common continuous subtree

(CCS) of (F , F ) is an ordered pair (H , H ) such that • H is a CS of F , H is a CS of F .
• Let e and e be the cut edges of (H , H ). Then C (H , e ) C (H , e ).

The number of leaves in C (H , e ) or equivalently in C (H , e ) is called the size of the CCS.

From the definition of a CS, after performing the cut, it must still remain at least two leaves in the CS, hence, the size of a CCS is at least two. Consequently, C (H , e ) and C (H , e ) share at least one common cherry.

Classification

In an irreducible pair of forests, the possible sizes or topologies of a CCS are very limited. As we just pointed out, after performing the cuts in a CCS, a common cherry appears in both forests. This common cherry must not exist before the cuts, otherwise the first reduction rule would apply. It implies that, in at least one of the two forests, the cut edge of the CCS is the edge that "separates" the two leaves in that cherry. The next lemma provides a list of properties of CCSes in an irreducible pair of forests, and classifies all possible CCSes into three types in relation to the cut edge that "separates the new common cherry".

Lemma . . Let (H , H ) be a CCS of a pair of irreducible X-forests (F , F ). Then . Th cut edges of H and H cannot both be empty.

. For either H or H or for both, the cut edge (u, v) is incident to the pendant edge of a leaf a, that is, u = p(a).

A CS satisfying this property is called a leading CS of that CCS.

. Assume wlog that H is a leading CS of the CCS. Let a be leaf adjacent to the cut edge . The maximum size of a CCS of type one is . The maximum size of a CCS of type two is . The maximum size of a CCS of type three is .

Proof.

. Suppose that both cut edges are empty, then H and H are common subtrees of F and F . Since the reduction rule does not apply, any common subtree can only have one leaf. By definition, a CS has at least two leaves.

. Suppose that neither of the two CSes is a leading CS, then there exists a common cherry in C (H , e ) and C (H , e ) that is unchanged by performing the cuts.

This cherry is a common subtree of F and F and contradicts the irreducibility.

. Trivial from the definition of a CCS.

. The limits on the size of a CCS come from the fact that neither reduction rule nor reduction rule applies. The assertion can be checked by listing exhaustively all possible configurations of a CCS. In Figure . , we provided such a list. The

Best CCS and best pre-CCS leaf pair

At each iteration, the algorithm identifies all CCSes of the irreducible pair of forests.

Two situations can then occur. If there are more than one CCS, it selects one of them.

If, on the other hand, there exists no CCS, it in turn tries to find a pre-CCS leaf pair, which we will define later; as the name suggests, after cutting out such a leaf pair from the forests, some CCSes will emerge as a result. In both situations, the choice of the CCS or the pre-CCS leaf pair is guided by a partial order between CCSes, taking into the type and the size.

Definition . (Best CCS). For a given set of CCSes, consider the following partial order (reflexive, antisymmetric, and transitive) "better than":

• A CCS of type one is better than a CCS of type two.

• A CCS of type two is better than a CCS of type three.

• For two CCSes of the same type, a CCS of larger size is better than the smaller one.

A CCS is called a best CCS if it is a minimal element, i.e., there exists no other better CCS in the set.

When the pair of forests does not have any CCS, the algorithm chooses to cut a leaf pair, or more precisely, the pendant edges of a leaf pair. Indeed, in the special case of an irreducible pair of forests without any CCS, this kind of cut preserves the irreducibility. This allows the algorithm to "look one step ahead" and to identify new CCSes immediately after the cut, i.e., without applying the reduction rules and going to the next iteration.

Lemma . . Let (F , F ) be a pair of irreducible X-forests that does not have any CCS. Let (v , v ) be a leaf pair of (F , F ). Then, after cutting the pendant edges e and e of v and v , the pair of forests (C (F , e ), C (F , e )) are still irreducible.

Proof. It is easy to check that the reduction rules and cannot apply to the pair (C (F , e ), C (F , e )). To see that the Reduction rule does not apply, suppose that (C (F , e ), C (F , e )) has one common subtree of at least two leaves. Let a and b be two leaves that form a common cherry of C (F , e ) and C (F , e ). The leaves a and b cannot appear as a cherry in both forests F and F the pair (F , F ) is irreducible.

Thus at least one edge e or e must separate the two leaves. Assume wlog that e is incident to the pendant edge of a in F . Let H be the CS of F rooted at the parent of b with e marked as the cut edge. Depending on the position of the edge e relative to the cherry (a ′ , b ′ ) in F (three cases: outside of the cherry, incident to the pendant edge of b ′ , or incident to the pendant edge of a ′ ), we can choose the CS H accordingly (each case corresponds to a type of CCS) such that (H , H ) is a CCS of (F , F ) (a and b form the target cherry of that CCS), contradicting the assumption that the pair (F , F ) does have any CCS.

Definition . (Pre-CCS leaf pair). Let (F , F ) be a pair of irreducible X-forests that does not have any CCS. A leaf pair (v , v ) of (F , F ) is called a pre-CCS leaf pair if, after cutting the pendant edges e and e of v and v , the new pair of forests (C (F , e ), C (F , e )) has at least one CCS.

Definition . (Best pre-CCS leaf pair). Let (F , F ) be a pair of irreducible X-forests that does not have any CCS. Consider the set of CCSes

H := (e ,e ) H , H ) : (H , H ) is a best CCS of (C (F , e ), (F , e ))
where e and e are the pendant edges of v and v , and the union is taken over all pre-CCS leaf pairs (v , v ) of (F , F ). A leaf pair (v , v ) for which a best CCS (H , H ) of (C (F , e ), C (F , e )) is a minimal element of H for the relation "better than" (see Definition . ) is called a best pre-CCS leaf pair associated with a best CCS (H , H ).

From this definition, if there exist two distinct best pre-CCS leaf pairs associated respectively with the best CCSes (H , H ) and (H ′ , H ′ ), then both best CCSes are of the same type and of the same size, although they are CCSes in two different pairs of forests.

. . Summary of the algorithm

The algorithm, shown in Algorithm , is called NewCCS. The name acknowledges the fact that it originated from ideas of algorithms that were previously explored in the team and in which the term CCS, though defined differently, was first introduced.

At each iteration, the algorithm updates a pair of forests (F , F ) and an integer k that keeps track of the number of cuts made so far in F . The edges that are cut from F by the algorithm come from three sources: either it is the result of applying the reduction rule (Line ), or it is the pendant edge of a leaf (Lines and ), or it is the cut edge of a CCS (Lines and ). As the cut edge e of a CCS may not exist (defined to be ∅ in this case), we use the notation |e| that equals zero if the edge is empty and equals one otherwise. At least one edge is cut in one of the forests in each iteration (by Lemma . , the two cut edges of a CCS cannot both be empty), so the while loop necessarily stops and ends with a pair of isomorphic forests. As the reduction rules We end this section with the following conjecture:

Conjecture . . Algorithm is a -approximation for computing m(T , T ).

.

Experimental results

Algorithm has been implemented in Python. Because there is some randomness, the implementation performs random runs and outputs the best In this section, the name NewCCS only refers to the best-out-of-ten implementation.

To compare with our NewCCS algorithm, we use five existing algorithms that have a publicly available implementation.

• WZ (https://github.com/cwhidden/rspr), the -approximation by [WZ ] which runs in linear time.

• SVN (https://nolver.net/maf/), the -approximation by [SZS ]. It runs in quadratic time.

• CMW, CHN, CombMCTS (http://rnc.r.dendai.ac.jp/rsprHN.html), respectively, the ( / )-approximation by [CMW ] (quadratic time), the -approximation by

[Che+ ] (cubic time), the Monte Carlo tree search algorithm by [YCW ], which is a practical improvement of the previous two (no complexity analysis).

. . Practical issues

Ideally, we would like to run all the above algorithms on all the available datasets. This was however not possible due to various runtime problems that we encountered when running the software.

• CMW can output a number that is smaller than the real distance. For example, for the following instance with distance , the output is .

( ((((((((L1,L2),L3),L4),L5),L6),L7),L8),L9),L10); (((((((((L1,L2),L4),L8),L9),L10),L5),L6),L3),L7);

• CMW can output a number that is more than / times the real distance, while it is a ( / )-approximation. For example, for the following instance with distance , the output is .

(((((((((L1,L2),L3),L4),L5),L6),L7),L8),L9),L10);

(((((((((L2,L4),L5),L1),L6),L7),L3),L8),L9),L10);

• CHN can output a number that is smaller than the real distance. For example, for the following instance with distance , the output is (however, the output is if the order of two input trees is switched).

(((L5,L4),L3),(((L9,L7),((L10,L8),(L6,L2))),L1));

((((L8,((((L9,L7),(L10,L2)),L1),L5)),L4),L6),L3); Any of the above problems will be referred to as an error. It should be clear that the possibility of running into errors limits the experiments and the analyses that we are able to perform.

We will now proceed to the comparison of NewCCS with the five existing algorithms. The datasets that we use are either from the literature or from our own construction. For any instance in the datasets, the exact distance is known (either by construction, or pre-computed using one of the available software, e.g., [WBZ ]). The experiments are run on a laptop PC with Intel i -M CPU ( . GHz, cores) and GB RAM.

. . Comparison using randomly generated datasets

For the first experiment, we use the three datasets R , R , and R from [YCW ], each containing instances. In R b a , each instance (T , T ) is obtained by first generating a random phylogenetic tree T with a leaves and then obtaining T by applying b random rSPR operations. For example, an instance in R has size , and the distance is at most .

On these three datasets, all six algorithms run fast enough, and without encountering any error, so we can compare their approximation ratio as the authors did in [YCW ]. For each dataset, we computed the average and the maximum approximation ratios among the instances. The results are shown in Table . . . : For each of the six algorithms and for each dataset, the average and maximum approximation ratio on the instances of the randomly generated datasets.

Average Ratio Max Ratio R R R R R
On the randomly generated datasets, the -approximation WZ produces slightly better ratios than the -approximation SVN. Those two algorithms are significantly outperformed by the ( / )-approximation CMW and the -approximation CHN. The lowest average approximation ratio is achieved by CombMCTS, which was designed to improve the approximation ratio in practice by means of a Monte Carlo tree search technique. The NewCCS algorithm yields the second best average approximation ratio on all three datasets.

Instead of looking at the quality of the output of each algorithm on average on the whole dataset, we can ask, for each individual instance, how well is the output of one algorithm compared to the output of another. For this, we can compute the average difference of the approximation ratio, that is, the difference between the two output values divided by the real distance, averaged over the instances. Take the dataset R . On average, NewCCS achieves an approximation ratio that is . better than CMW, .

better than CHN, and . worse than CombMCTS.

. . Comparison using a dataset of caterpillars of fixed size

We have constructed a dataset consisting of pairs of trees with distances between and , which correspond to all possible pairs of caterpillars with leaves (a rooted caterpillar is a rooted binary tree with only one cherry). Since the number of instances in this dataset is quite large, we randomly selected a subset of instances (with distances also spanning from to ). On this subset, we ran all six algorithms and measured the running time. We also computed the average and the maximum approximation ratios. If an error occurs (see Section . . ), the instance is excluded from the computation of the average and the max approximation ratios (for that algorithm only). The results are shown in Table . .

. . Comparison using a special family of instances

For the last experiment, we used a dataset that we call HUGE. It consists of pairs of trees of sizes between and . These instances belong to the I ′ family introduced by [RSW ] and for which the distances are known exactly. Precisely, for an instance of size q + in this dataset (where q is an integer), the distance is q + . For every instance (T , T ) in this dataset, another instance (T , T ) is also included, that is, the same pair of trees with the order switched.

We ran all six algorithms on HUGE and recorded the approximation ratio and the running time. The results are plotted in Figures . -. Table . : Behaviors of the approximation ratio of each algorithm with increasing instance sizes in the dataset HUGE, for each of the two possible orderings of the input.

Interestingly, apart from SVN and NewCCS, the other four algorithms all exhibit significantly different behaviors in terms of the approximation ratio when the order of the input trees is switched. Such observations are summarized in Table . . For two of those algorithms, namely CMW and CombMCTS, the running time also follows two distinctive growth patterns, i.e., the running time grows more rapidly with the instance size for one of two orderings of the input, and grows more slowly for the other ordering.

Among the six algorithms, WZ and CHN are the fastest but achieve the worst approximation ratios; SVN always produces an approximation ratio near the worst factor ; the performance of CMW is most sensitive to the ordering of the input. Only

CombMCTS and NewCCS consistently produce near optimal solutions (approximation ratio < . ) for both orderings of the input. However, CombMCTS can be more than a thousand times slower than NewCCS (e.g., min versus s), making it impractical on this family of instances as soon as the instance size becomes larger than a few hundreds. We have confirmed again the observation that we made in the previous experiment: NewCCS produces results of good quality, is fast enough, and is robust against special cases (no errors, no sensitivity to the ordering of the input).

.

Conclusion and perspectives

We proposed a new algorithm NewCCS for finding an approximate solution of the M -. Although its approximation ratio is currently unknown, it performs well in practice, beating existing approximation algorithms either in terms of the quality of the solution or in terms of speed. The immediate perspective is to prove or disprove Conjecture . . There are a number of small known instances on which NewCCS achieves the factor . Based only on few simple ideas,

NewCCS is easier to understand and to implement than the two algorithms that achieve the current best approximation ratio of , and it has the potential to be improved in the future. Finding a better approximation algorithm for computing the rSPR distance is an interesting question both theoretically and in practical applications. Any improvement in the approximation factor or in the running time can possibly lead to better exact algorithms for computing the rSPR distance (e.g., by using a bounded search technique), as well as better exact or approximation algorithms for related problems such as the computation of the hybridization number.

Conclusion

The interaction between organisms living in closely related environments, such as humans and their gut microbiota, can present complex dynamics. To understand it through the lens of coevolution, the field of computational biology provides powerful methods for finding the most parsimonious reconciliations between phylogenetic trees. In the case where efficient enumeration algorithms exist, that is, when we have put behind us the computational complexity of the problem itself, what we are facing now is the biological complexity of the solutions of the problem: in a hosts/parasites system consisting of more than a thousand species, any of the current event-based models would generate a solution space that is so large that it precludes the analysis or the interpretation of the results [Bau+ ]. There are two natural ways of addressing this issue: either we refine the model by including more constraints, hoping that the output size will be significantly reduced (and hoping that we can still solve the problem efficiently), or, we cling to this huge solution space but try to summarize it somehow, so that we get an output of reasonable size without losing important biological information. This thesis is concerned with the second approach.

In Chapter , we defined different notions of equivalence between reconciliations.

To summarize the solution space of the phylogenetic tree reconciliation problem then means to provide a list of equivalence classes. In each equivalence class, any solution can be chosen as the representative, as the biological characteristics that we want to capture with the equivalence relation are necessarily shared by any solution in the class. We designed practical algorithms for enumerating the equivalence classes and the representative solutions, and in Chapter , applied those algorithms to biological datasets and made them available to end users through the Capybara software.

Our approach is remarkable from several points of view. First of all, for the phylogenetic tree reconciliation problem, the practical benefit of considering the equivalence classes when analyzing the solution space is striking. For example, in Chapter , Section . . , we managed to reduce the number of solutions that need to be considered from the order of to only . Beyond the reconciliation problem, many more problems in bioinformatics or computer science also suffer from a large size of the solution space. Compared to other possible approaches that can drastically reduce the size of the output, such as agglomerative clustering methods that produce consensus solutions or centroid solutions [Ozd+ ; ML ; San+ ; AQE ], the equivalence-based approach has quite a few advantages.

The first one is that, since we only need to define an adequate notion of equivalence, the equivalence-based approach can be used in a wider context than clustering methods. The latter require some kind of distance or diversity measure between solutions which, depending on their mathematical nature (numerical values, strings of symbols, graphs, functions on graphs, etc.), can be difficult to define or costly to compute (see also Chapter , Section . ).

The second advantage is the interpretability of the results. One reason for this is that the enumeration (and the counting) of equivalence classes, as we formulated it, is not a data analysis method but instead a computational problem that could be solved exactly. Clustering methods will perform poorly when the solution space is not wellseparated into clusters with respect to the chosen distance measure, and, when there are indeed some distinctive clusters in the distribution of the solutions, such methods would fail to discover the pattern if the number of "real" clusters is too large (this is because the running time usually depends heavily on the size of the output, i.e., the final number of clusters chosen by the user). Quite unlike the clustering methods, an algorithm that counts or enumerates the equivalence classes will invariably discover any "distinctive pattern" with respect to the chosen equivalence relation, as a pattern here corresponds simply to an equivalence class. The equivalence relation also provides a precise characterization of the solutions inside an equivalence class. In contrast, solutions inside a cluster in the clustering output are described by their similarities or dissimilarities with the representative solution (the centroid of the cluster): these are numerical values which can be less straightforward for the user to interpret.

For all these reasons, we believe that, outside of the context of the reconciliation problem, there are more problems that can benefit from considering the equivalence classes of solutions. In Chapter , we proposed a general framework for enumerating equivalence classes of solutions that can be useful for a variety of problems.

This general framework has a strong connection with the key Algorithm of Chapter that enumerates E-equivalence classes in polynomial delay for the reconciliation problem. First, while Algorithm is described in terms of the reconciliation graph which might sound limited to that particular problem, we pointed out early on that this graph has the same structural properties as ad-AND/OR graphs, a type of graph that naturally emerges in problems that can be solved by dynamic programming algorithms. Moreover, since Algorithm only employed very general techniques such as recursive depth-first search and backtracking, it is not at all surprising that the idea can be generalized to tackle a wider range of problems.

Making sense of the output by efficiently summarizing the solution space is only one way of increasing the practical usability of an algorithm. In Chapter , we saw that even for the well-studied problem of computing the rooted subtree prune-and-regraft distance between two phylogenetic trees, current algorithms stumble upon various practical issues: quality of the solutions, speed, or the robustness of the implementation when tested on special datasets (see Section . ). We thus proposed a new algorithm that is easy to understand and to implement, is fast enough, and consistently produces solutions of high quality on randomly generated instances as well as on special families of instances.

Other than the open mathematical and computational challenges mentioned at the end of each chapter, there remains several avenues of expansion for the work of the present thesis. It will be interesting to see how well the results of equivalence classes enumeration match up with the biological understanding of the coevolutionary history of the organisms in other datasets, and maybe under new notions of equivalence that yet need to be defined. Then, based on the feedback from end users, new features can be added to Capybara to include more computational and visualization options. It will also be interesting to apply the general framework of Chapter to study some different problems in order to see how well it works in practice. We provided a few examples that can be a good starting point, and we would like to find more applications. For problems whose solutions do not have a tree structure, the enumeration of equivalence classes as a means of understanding the solution space can still be worth investigating, although new techniques will be necessary.

  celui d'hôtes et celui de parasites), et d'un ensemble d'associations qui indiquent quel organisme d'un arbre interagit au temps présent avec quel organisme de l'autre arbre. L'objectif est d'associer les espèces ancestrales entre les deux arbres, d'une façon qui soit la plus judicieuse possible, et ce, en termes d'événements coévolutifs représentés par ces associations. Le modèle de réconciliation d'arbres, proposé par des bioinformaticiens, permet de transformer le problème cophylogénétique en un problème d'optimisation combinatoire. Pour trouver une solution optimale d'un problème d'optimisation, le développement d'algorithmes est un outil essentiel. S'attaquant à ce problème informatique, notre investigation est donc de nature algorithmique.

  Chapitre , je présente brièvement quelques concepts biologiques et mathématiques et introduis le défi de trouver des solutions représentatives du problème de cophylogénie. Le Chapitre est consacré aux principaux résultats algorithmiques. D'abord, je définis plusieurs relations d'équivalence sur l'ensemble de solutions optimales et explique la motivation biologique de ces définitions. Mathématiquement, une relation d'équivalence partitionne l'espace de solutions en des sous-parties appelées « classes d'équivalence ». Dans ce contexte, on considère comme étant représentatives les solutions, non équivalentes entre elles, qui appartiennent à chacune des classes d'équivalence. Je présente ensuite des algorithmes permettant d'énumérer les classes d'équivalence pour chaque relation d'équivalence. Afin d'illustrer les avantages de considérer les classes d'équivalence en pratique, des résultats expérimentaux sur des jeux de données biologiques sont donnés dans le Chapitre . Ce chapitre présente aussi le logiciel Capybara, un outil d'aide à l'analyse des données de cophylogénie qui dérive de notre méthode d'exploration de l'espace de solutions par le biais des classes d'équivalence. Notre investigation algorithmique du problème cophylogénétique nous mène à réaliser deux autres investigations liées à la première. Dans le Chapitre , on se propose de s'attaquer à la question d'énumérer les classes d'équivalence de solutions dans un contexte plus général, où une solution est identifiée par un sous-arbre dans un graphe, et la relation d'équivalence est donnée par un coloriage de ce graphe. Je présente un algorithme d'énumération de classes d'équivalence et montre que celui-ci pourrait s'appliquer à une famille de problèmes différents, dont fait partie le problème de cophylogénie. Le Chapitre aborde la question de la quantification des différences entre des arbres phylogénétiques. De nombreuses applications nécessitent un calcul rapide et précis d'une distance entre des arbres. Par exemple, dans un système d'interactions hôtes-parasites soumis à une forte pression favorisant la coévolution, on s'attend à constater une faible dissimilarité entre les deux arbres phylogénétiques. Je propose un nouvel algorithme pour approximer la distance dite Subtree prune and regraft entre deux arbres, et je réalise une série d'expériences visant à illustrer ses avantages pratiques par rapport aux algorithmes existants. Les principaux résultats algorithmiques et une partie des résultats expérimentaux des Chapitres et se trouvent dans l'article « Making sense of a cophylogeny output : Efficient listing of representative reconciliations », publié dans le Workshop on Algorithms in Bioinformatics (WABI) [Wan+ b]. Le logiciel Capybara a été présenté dans l'article « Capybara : equivalence ClAss enumeration of coPhylogenY event-BAsed ReconciliAtions », publié sous forme d'une Application Note dans la revue Bioinformatics [Wan+ ]. Le contenu du Chapitre constitue l'article « A general framework for enumerating equivalence classes of solutions », publié dans l'European Symposium on Algorithms (ESA) [Wan+ a].

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . List of Tables . Definition of the event type function and the loss number function. . . . . . Definition of the color of a parasite node under a reconciliation. . . . . . . . Definition of the equivalence class of a reconciliation. . . . . . . . . . . . . . . A function MakeIntersectionFree that satisfies Lemma . . . . . . . . . . . . Two cases for a partner pair, and the corresponding replacement. . . . . . . All the datasets that were used in the experiments and their references. . . The number of optimal reconciliations and the number of equivalence classes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The list of V-equivalence classes of the instance WOLB-( , , , ). . . . . . . The average and maximum approximation ratio for each of the randomly generated dataset. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The running time, average and maximum approximation ratio for the datasets of pairs of caterpillars with leaves. . . . . . . . . . . . . . . . Behaviors of the approximation ratio with increasing instance sizes in the dataset HUGE. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv List of symbols and abbreviations ch The set or tuple of children of a node. DAG Directed acyclic graph. v w w Two nodes are incomparable. d T (v, w) Distance between two nodes in a tree.LCA (v, w) Lowest common ancestor of two nodes in a tree.

  ters and are included in the paper "Making sense of a cophylogeny output: Efficient listing of representative reconciliations", published in the Workshop on Algorithms in Bioinformatics (WABI) [Wan+ b]. The cophylogeny software Capybara has been presented in "Capybara: equivalence ClAss enumeration of coPhylogenY event-BAsed ReconciliAtions", published as an application note in the journal Bioinformatics [Wan+ ]. The material in Chapter forms the paper "A general framework for enumerating equivalence classes of solutions", published in European Symposium on Algorithms in the study of evolution is the understanding of evolutionary relationships between organisms, both living and extinct. To describe such relationships, a conceptual tool, dating back to the mid-nineteenth century, is the so-called Tree of Life. Charles Darwin in his famous book The Origin of Species presented a first graphical representation of evolutionary relationship among species, in the form of a tree [Dar ]. In evolutionary biology, phylogenies, or phylogenetic trees, have been a model of choice, on the tenet that contemporary species all share a common history through their ancestry [DBP ].

Figure . :

 . Figure . : The tree of life constructed from genomic data by Hug et al. [Hug+ ].

Figure . :

 . Figure . : Illustration of the four types of coevolutionary events in a host/parasite reconciliation. The tube represents the host tree; the dotted line represents the parasite tree.

  can be compactly represented in a graph structure. By carefully examining this graph, we can gain insight into the solution space and obtain information about the coevolutionary history of the species. Various approaches have been developed in the literature that aim at efficiently exploring the solution space or computing some characteristics on the set of solutions, including data analysis methods (e.g., sampling solutions uniformly at random [BAK ; Don+ ], clustering [Ozd+ ; ML ; San+ ]) and exact methods (e.g., computing the diameter of the solution space [Haa+ ], computing the distribution of the pairwise distances [SML ], computing a median reconciliation [Ngu+ ]).

  an ordered pair instead of a set. The two children of v are respectively called the left child and the right child of v. When (v , v ) are the children of v, the two subtrees T | v and T | v are called respectively the left subtree and the right subtree of T at node v.

Figure . :

 . Figure . : Four phylogenetic trees and their representations in Newick format. T and T are on the same set of taxa {A, B, C} but are not isomorphic. T and T are on the same set of taxa {A, B, C, D} and are isomorphic.

Figure . :

 . Figure . : Example of six different reconciliations on the same dataset. For each reconciliation, the parasite and the host trees are drawn side by side. The solid curves at the bottom indicate the parasite/host associations for the parasite leaves; the dashed lines indicate the associations for the internal parasite nodes.

  we are interested in finding reconciliations of minimum cost. It turns out that this problem can be solved efficiently using the dynamic programming technique [BAK ; Don+ ]. The DP-algorithm also solves the enumeration version of the problem: in fact, it produces a graph structure which is a compact representation of all reconciliations of minimum cost [Don+ ; Ma+ ]. One can argue that the power of this model

  Figure . ). Let n = |V(H)| and m = |V(P)|. It is clear that the reconciliation graph has O(mn) OR + nodes. Each OR node has O(n ) children (or, equivalently, O(n ) grandchild couples).

  Let n = |V(H)| and m = |V(P)|. As we mentioned, the reconciliation graph can be constructed in O(mn ) time. After the construction, using the standard techniques for ad-AND/OR graphs, the total number of optimal reconciliations can be easily computed, and the solutions subtrees can be enumerated with a linear time delay, that is O(m). Therefore, there is an algorithm with a O(mn ) time pre-processing step and O(m) time delay for enumerating all optimal reconciliations.

Figure . shows a

  Figure . shows a reconciliation graph based on the same dataset (H, P, σ) as in Figure . with nine solution subtrees. Among these nine reconciliations, four have event vector ( , , , ), two have ( , , , ) (φ and φ of Figure . ), two have ( , , , ) (φ and φ of Figure . ), and one has ( , , , ) (φ of Figure . ). The reconciliation shown in bold is φ of Figure . .

Figure . :

 . Figure . : Example of a reconciliation graph for the dataset (H, P, σ) of Figure . . Crossed circles are AND nodes. Rectangles are OR + nodes. The cells with which the OR + nodes are labeled are written inside. One solution subtree is shown in bold (it is φ of Figure . ).

For

  example, φ and φ (or φ and φ ) of Figure.are V-, V * -, E-, and CD-equivalent but not EL-equivalent; φ and φ are V * -equivalent, but they are not equivalent under any other equivalence relation.

Figure . :

 . Figure . : Illustration of Lemma . . Circles are parasite nodes; rectangles with rounded corners are host nodes; dashed lines represent the parasite/host associations; the zigzag line indicates that (p , p ) is the host-switch arc at p under φ.

  Let n = |V(H)| and m = |V(P)|. The first three elements of any event vector necessarily sum up to m-, the number of internal parasite nodes, hence there are only O(m )

  enumerate the E-equivalence classes by a top-down recursive traversal of the reconciliation graph, the algorithm should achieve the following goal: given a set O of OR + nodes of the same color (p, e), enumerate π(T(G/O)), i.e., all E-equivalence classes of the subgraph G/O. Any such class will include the color (p, e). If p is not a leaf, the possible events of the two children of the node p are dictated by the color couples of the grandchild couples gch(O) (by extension, gch of a set of nodes is the union of gch of every node in the set). A naive algorithm can be described as follows: for each color couple ((p , e ), (p , e )) of gch(O), first take the union O of the first grandchildren of color (p , e ) and the union O of the second grandchildren of color (p , e ), then call the algorithm on O and independently on O , and finally combine the results together, that is, perform a Cartesian product between π(T(G/O )) and π(T(G/O )).

  only reconciliation in T(G/O), therefore, {(p, e)} is the only E-equivalence class. The outputted set O contains in this case the unique goal node of G labeled by (p, σ(p), ). Now we assume h p > . (First direction) Consider a fixed pair of T := T ∪ T ∪ {(p, e)} and O outputted at Line , and take a node w in O. We show that there exists a reconciliation φ ∈ T(G/O) such that w = r(φ) and π(φ) = T (in other words, T is a valid E-equivalence class in π(T(G))). By the induction hypotheses, T is an E-equivalence class so there exists a reconciliation φ of P| p such that π(φ ) = T . Let w := r(φ ). Take a node w ∈ O such that (w , w ) ∈ gch(w). By the induction hypotheses, there exists a reconciliation φ of P| p such that r(φ ) = w and π(φ ) = T . Define φ := φ ∪ φ ∪ {(p, h, e)}, where (p, h, e) is the label of w. Then φ is a valid reconciliation in T(G/O) (notice that φ is a solution subtree of G/O if and only if (w , w ) ∈ gch(w)), and satisfies π(φ) = T . (Second direction) Consider an E-equivalence class T ∈ π(T(G/O)), and take a reconciliation φ ∈ T(G/O) such that π(φ) = T . We show that T is outputted exactly once at Line together with a set O containing the root node of φ. Assume that the root node w := r(φ) is labeled with the triple (p, h, e), then φ can be uniquely written as the union φ ∪ φ ∪ {(p, h, e)} where φ and φ are respectively reconciliations of P| p and P| p . Furthermore, T can be uniquely written as the union T ∪ T ∪ {(p, e)} where T = π(φ ) and T = π(φ ). Notice that T and T do not depend on the choice of φ; for T to be outputted exactly once, it suffices to show that each of T and T is outputted exactly once. For i = , , let w i := r(φ i ) and let (p i , e i ) := κ(w i ) be the color of w i . At Line , we only need to consider the iteration corresponding to the color couple ((p , e ), (p , e )), as no other iteration can output T or T from a recursive call. Since w ∈ O and φ ∈ T(G/O ), by the induction hypotheses, T is outputted exactly once in Line together with a set O containing w . For this pair of T and O , the set O computed at Line contains the node w . Hence, by applying again the induction hypotheses to φ ∈ T(G/O ), T is outputted exactly once in Line together with O containing w . It remains to check that the set O outputted together with T does contain the node w. As we have (w , w ) ∈ O × O , this is straightforward from the computation of O (Line ). Theorem . . Using Algorithm , the set of E-equivalence classes of a reconciliation graph can be enumerated in O(mn ) time delay, where m = |V(P)| and n = |V(H)|. Proof. To obtain all E-equivalence classes in π(T(G)), it suffices to first partition the set of start nodes of the reconciliation graph according to their colors, then, for each subset O i of start nodes of color (p, e), make one call of Enumerate(p, e, O). By Lemma . , we output every E-equivalence class of T(G/O) exactly once. Since any E-equivalence class of T(G) is an E-equivalence class of T(G/O k ) for a unique k, we output every Eequivalence class of T(G) exactly once. For the complexity, consider the recursion tree formed by the recursive calls of Enumerate. Notice that each node p of the parasite tree corresponds to exactly one recursive call, the size of the recursion tree is thus O(m). In each recursive call, the partitioning of gch(O) and the computation of the sets O , O , and O can all be done in time linear in the size of gch(O), which is O(n ). Therefore, O(mn ) time is needed in the worst case between outputting two E-equivalence classes.

  the output consists of three AND nodes. In cases (c) -(g), one of the two input AND nodes does not undergo any modification and is outputted directly; the modification made on the other input node (in case (g), it is also split in two) is based on the computation of intersection and difference of OR + nodes (last column of Table .

  Let w and w be the two input nodes of InterDiff. If their colors are different, it is clear that C(w ) ∩ C(w ) = ∅ and C(w ) \ C(w ) = C(w ). Suppose now that they have the same color. If they are both goal nodes, the sets of colorings C(w ) and C(w ) have only one element, that is, their common color. In this case, C(w ) ∩ C(w ) = C(w ) and C(w ) \ C(w ) = ∅. Suppose now that w and w are both OR nodes. The idea is to express the intersection and the difference using recursive calls of InterDiff on OR + nodes that can be found in the grandchild couples of w and w . This idea turns out to work well: we are able to find a practical algorithm. Before presenting the algorithm, it is convenient to introduce a subroutine, called AND-Inter-Diff. It is a function that takes two parameters, a set A of AND nodes and a single AND node v j , all of them neatly color P| p for some p ∈ V(P). It returns two sets of AND nodes A inter and A diff . All the sets of AND nodes (one in the input and two in the output) satisfy additionally the property of having pairwise empty intersections as sets of colorings (see Equation (.)). The output should satisfy the following:

  given that we know how to compute the intersection and the difference between two OR + nodes. This is done by the InterDiff function which builds the set of child AND nodes of the two output nodes (intersection and difference nodes) in a nested loop over the children of the input nodes.The InterDiff function is therefore the most basic operation and is needed at least O(n ) times, where n is the number of child AND nodes in all the input nodes of Merge.The only nontrivial operation in the nested loop of InterDiff is the subroutine called AND-InterDiff which applies InterDiff recursively while iterating through a set of AND nodes. The performance of InterDiff is sensitive to the size of those sets of AND nodes. Any set of AND nodes that are considered in the algorithm contains only AND nodes that are pairwise disjoint as sets of colorings, and it is constructed in two ways, either ( ) by initializing with a small number of AND nodes, or ( ) by performing the union operation between pairwise disjoint sets. For the purpose of improving the performance and reducing the memory required for storing the result AND/OR graph, we Algorithm : Intersection and difference between two OR + nodes Input: Two neat OR + nodes that color P| p for some p ∈ V(P).Output: Two neat OR + nodes w inter and w diff . Function InterDiff(w , w ): if w and w do not have the same color, that is, κ(w ) ≠ κ(w ) then return w inter := Empty and w diff := w end Let p ∈ V(P) such that w and w color P| p if p is a leaf then return w inter := w and w diff := Empty end Let (p , p ) be the children of p Let A inter and A diff be empty sets for each child AND node v i ∈ ch(w ) do

  Let w inter be an OR node with κ(w inter ) := κ(w ) and ch(w inter ) := A inter Let w diff be an OR node with κ(w diff ) := κ(w ) and ch(w diff ) := A diff /* If A inter or A diff is empty, the corresponding output is set to Empty */ return w inter and w diff now introduce another operation, called RemovePartners, that decreases the size of a set of AND nodes (without changing, of course, the sets of colorings). It can be used anywhere in AND-InterDiff and InterDiff, immediately after a union operation. It can also be used in the function MakeNeat which also operates on a set of AND nodes:

FD

  |V(H)| |V(P)| Notes and references EC Encyrtidae (parasitic wasps) and Coccidae [Den+ ] Lice, the trees are generated using either a maximum likelihood (L) or a maximum parsimony (P) approach [Hug+ ] trees from the Clusters of Orthologous Groups of proteins (COG) database. The unrooted gene trees are rooted in such a way that the minimum cost of the R is minimized among all choices of root. [Tat+ ; DA ] WOLB Wolbachia (bacteria) and their arthropod hosts [Sim+ ; Sim ]

  Figure . : X-axis:The instances in increasing order of the number of optimal reconciliations. Y-axis: In logarithmic scale, the Reduction value that is equal to the number of equivalence classes over the number of optimal reconciliations. For each instance, three points are plotted: the blue circle, the red triangle, and the black X, corresponding respectively to the number of E-, EL-, and CD-equivalence classes.

  [Nøj+ ], the authors presented an algorithm that answers the following question in O(m log(n)) time, where m = |V(P)| and n = |V(H)|:

p

  Figure . : Visualization of the four E-equivalence classes the six reconciliations from Figure . . The colors corresponding to , ℂ, , event types are respectively white, orange, green, and blue.

  interesting and challenging topic in the field of enumeration algorithms. Unlike the other chapters of this thesis which all have to do with phylogenetic trees, the current chapter is concerned with the problem of enumerating equivalence classes of solutions in the context of DP-problems. The algorithm that we will propose can be seen as a generalization of Algorithm of Chapter ; this generalization is made from a partic-ular DP-problem (the R ) to a whole family of DP-problems.Generalizing an algorithm to a wider range of problems is not a trivial task. It is crucial to have a perspicuous formulation of the generalized version of the equivalence classes enumeration problem. To this end, we conveniently use the framework of ad-AND/OR graphs to represent the common structure of those DP-problems, and we define a specific kind of node-coloring on ad-AND/OR graphs, with the aim of characterizing the equivalence relations that we are allowed to consider on the solutions of the DP-problems.The remainder of the chapter is organized as follows. After presenting the background and the motivation, we will first define an intermediate problem, namely the C on ad-AND/OR graphs, and give a polynomial delay algorithm. Then, as a second step, we show the connection between the colored classes of an ad-AND/OR graph and the equivalence classes of solutions of a DP-problem, and provide several concrete examples to illustrate how our algorithm can be applied to a wide variety of DP-problems.

  lence classes have been identified in various areas, such as Genome rearrangements [Bra+ ; BS ], Artificial intelligence [AMP ], Pattern matching [Blu+ ; Nar+ ], or the study of RNA shapes [GVR ]. The efficient enumeration of representative solutions to a problem was listed as an important open problem in a Dagstuhl workshop in on Algorithmic Enumeration: Output-sensitive, Input-Sensitive, Parameterized, Approximative [FGS ]. This challenge has been addressed in the literature in a few specific problems and using different approaches, such as diversity measures [AK ], super-solutions [MSS ], and equivalence classes [Bra+ ; Mor ]. To the best of our knowledge, no generalized method for finding representative solutions has been designed to deal with a family of problems.

Figure . :

 . Figure . : An acyclic decomposable AND/OR graph with four start nodes. Squares are OR + nodes (OR nodes or goal nodes); crossed circles are AND nodes. One solution subtree of size is shown in bold.

  . (Subgraph starting from a set of nodes). Let G be an ad-AND/OR graph. Let O be a set of OR + nodes of G. The subgraph of G starting from O, denoted by G/O, is the subgraph obtained from G by setting O as the new set of start nodes (i.e., by removing all nodes are unreachable from O through directed paths). When the set O consists of a single node w, the notation G/O simplifies to G/w. It may seem trivial but let us mention that the number of solution subtrees |T(G)| of an ad-AND/OR graph G can be computed easily, without enumerating any of them.

  Figure . .

GFigure . :

 . Figure . : An e-colored ad-AND/OR graph and its five color classes. Crossed circles are AND nodes; squares are OR + nodes. The colors of the OR nodes are written inside the squares.

  e., ct (w) := c(v) : v ∈ ch(w) . From the definition, a color tuple (c , . . . , c j ) belongs to ct (w) if w has a child AND node whose children are colored with (c , . . . , c j ). If we consider a color class C of C(G/w), the tuples in ct (w) are precisely the possible colorings of the children of r(C) (in other words, the root colors of the child subtrees of r(C)). Indeed, if the child AND node v ∈ ch(w) is chosen in a solution subtree, then c(v) will be the colors of the children of w in the color class in which that solution subtree belongs. Notice that several child AND nodes of w may have the same color tuple. If w is a goal node, we define ct (w) := ∅. We extend the notation ct to a set O of OR + nodes with ct (O) := w∈O ct (w). In the same way, (c , . . . , c j ) is a color tuple of ct (O) if and only if there exists a color class C of C(G/O) such that the children of r(C) are colored with (c , . . . , c j ).

Figure . :

 . Figure . : Left panel: An e-colored ad-AND/OR graph, where each node is labeled by an integer. Right panel: For O = { , } and t = (x, y), we have P = { , } and P = { , }. There are four combinations between C G/P and C G/P ; only two of them are admissible.

  a set of OR + nodes O all colored with the same color c and a color class C of C(G/O), we denote by r(π -(T )) := {r(T ) : T ∈ π -(C)} the subset of OR + nodes in O, each one of which is the root of a solution subtree belonging to the color class C. The following definition formalizes the notion of compatible nodes mentioned previously.

  Now consider any solution subtree T of G/O containing the nodes v, w , . . . , w k , and w k+ . Then its first (k+ ) child subtrees at v are (T , . . . , T k , T k+ ), and therefore, the first (k+ ) child subtrees at the root of its color class π(T ) are precisely (C , . . . , C k , C k+ ). Analysis Lemma . . Let C be a color class of G/O for a set O of OR + nodes of G, all colored with the same color c. Then, the function Next of Algorithm satisfies the following: . Next(c, O, ∅) returns the smallest equivalence class of G/O with respect to ≺. . If C is not the last color class of G/O, then Next(c, O, C) returns the successor of C. . If C is the last color class of G/O, then Next(c, O, C) returns ⊥. Proof. Let us define the height of the subgraph G/O, denoted by h(G/O), to be the maximum height of a color class of G/O, i.e., the number of OR nodes in a longest path from O to a goal node minus . The proof will be done by induction on h(G/O). Assume first that h(G/O) = , i.e., O contains only goal nodes. Then G/O has only one color class C, which is the single-node tree of color c. The call of Next(c, O, ∅) will output C at Line of the algorithm, and the call of Next(c, O, C) will return ⊥ at Line since |ct (O)| = . Assume now that h(G/O) > . Proof of Lemma . . If O contains goal nodes, then the smallest color class of G/O with respect to ≺ is the single-node tree colored with c and Next(c, O, ∅) outputs it at Line . Otherwise, let C ′ be the smallest color class of G/O and let (C ′ , . . . , C ′ j ) be the its child subtrees at the root. By definition of ≺, the roots of (C ′ , . . . , C ′ j ) are colored with the lexicographically smallest color tuple t := (c , . . . , c j ) of ct (O), and by Lemma . ,

  Lemma . . Assume now that C is the last equivalence class of G/O. Notice that in this case, the children of r(C) are colored with t r where r = |ct (O)|. As previously, let C ′ i be the i-th child subtree at r(C). Then, for all i ≤ j := |t r |, C ′ i is the largest color class of G/O i . By the induction hypothesis, Next(c i , O i , C ′ i ) would return ⊥ for all i ≤ j. Therefore, C will receive ∅ at Line and the algorithm will return to Line . Since r = |ct (O)|, r will receive r + at Line and ⊥ will be returned at Line . Theorem . . Given an e-colored ad-AND/OR graph G, the set C(G) can be enumerated with delay O(n • s) where n is the number of nodes of G and s is the maximum size of a color class. Proof. To enumerate C(G), we first partition the set of start nodes of G according to their colors. For each subset O i of color c i , starting with C := ∅, we repeatedly assign the output of Next(c i , O i , C) to C and output it until C = ⊥. By Lemma . , this guarantees that we output every color of O(G/O i ) exactly once. Since any color class of G belongs to C(G/O i ) for a given i ≤ k, every color class will be outputted exactly once.

  the easiest way is to return the set r(π -(C )) together with C ′ when the call Next(c, O, C) returns. To compute this, observe that if C ∈ C(G/O) where O is a set of goal nodes all having the same color c, then r(π -(C)) = O, and if r(C) has child subtrees (C , . . . , C j ), then r(π -(C)) is the set of nodes of O that has at least one child AND node v such that ch(v) contain exactly one node in r(π -(C i )) for each i ≤ j, which can be found in O(n) time. Thus only O(n) time is necessary at each recursive call to return r(π -(C ′ )) in addition to C ′ .

Theorem . .

 . The set of solution subtrees of G C , the graph returned by Algorithm , is equal to π - G (C), i.e., the set of solution subtrees of G belonging to the color class T . Algorithm : Restricting the graph to a color class Input: An e-colored ad-AND/OR graph G and a color class C Output: The subgraph G C of G restricted to C Function Main(G, C): for each start node w of G such that c(w ) = c(r(C)) do VisitOR + (w , r(C)) end return G C obtained from G by removing all unmarked nodes Function VisitOR + (w, x):

  to a new root of symbol a. While it might seem restrictive to require an additive structure on the cost function, this simple case does cover many important problems admitting a DP-algorithm, for instance, T[Bel ; HK ], K[KPP ], or L [WF ].

  hard[WJ ; War ]. Here we look at the special case of sequences of length (which also solves the problem of aligning, on a tree, sequences of a given fixed length). The input is a tree T , a set Σ of letters (DNA alphabet or protein alphabet), a function that labels each leaf node of T with a letter, and a distance function d : Σ → ℝ ≥ between two letters. The goal is to extend the leaf labeling to a full labeling φ : V(T ) → Σ such that the sum of the distances over the edges of is minimized. Defining equivalence relations of the solutions based on a grouping of the letters is uncontrived in this problem: for the DNA or protein alphabet, the letters can be subdivided into structurally similar groups. In practice, we have a fixed length n for the sequences on each node, and we solve the single-letter problem independently for each of the positions and then combine the solutions together. It is clear that, by combining n sets of solutions, the number of solutions can quickly become huge when the length grows. Summarizing the solution space by means of enumerating the equivalence classes is then a useful technique for analyzing the output.ExampleThe F problems are a family of problems that naturally arise in telecommunication networks, and that have been extensively stud-ied in graph theory as a generalization of graph coloring known as the T-coloring problem [Hal ; Rob ; Tes ]. In the variant called the list T-coloring, the input is a graph G representing the interference between radio stations, a set C of colors, a function S that gives for each vertex v ∈ V(G) a set S(v) ⊆ C of colors (possible frequencies for a station), and a set T ⊆ C of forbidden pairs of colors (interfering frequencies).

  acyclic and decomposable, the other is that the equivalence relation must be encoded into an e-coloring of the graph. It would be interesting to ask the same question for a general colored AND/OR graph, that is, whether we can efficiently enumerate color classes of solution subgraphs, if one or more of those ingredients is missing. In particular, it remains open whether the problem is hard without assuming the decomposability of the structure of the solutions.Putting together equivalence relations and graph enumeration problems can give rise to new research directions. We worked on a particular flavor of AND/OR graphs due to their close relationship with dynamic programming problems. Outside of a DP setting, in the area of graph algorithms, the enumeration of subgraphs of a graph satisfying some particular properties, of which the enumeration of color classes of an ad-AND/OR graph is a special case, is also a challenging problem, and has stimulated interesting development -both theoretical and experimental -in Network biology [ASA ; Geo+ ], Data mining [Mou+ ], Data bases theory [Zen+ etc.

  obtained from the forest by deleting at most one edge. It is necessarily a rooted binary tree. Unless explicitly stated otherwise, we simply use the term subtree to designate a pendant subtree. A cherry is a subtree consisting of only three vertices: two leaves and their common parent. A cherry containing two leaves a and b is denoted by (a, b), or equivalently, (b, a).

  (d , d ) be an ordered pair of integers satisfying d = d + |F | -|F |. • We say that (d , d ) is an upper bound of the distance, denoted by m(F , F ) ≤ (d , d ), if there exists an AF F of F and F such that |F| ≤ d + |F | = d + |F |. Equivalently, m(F , F ) ≤ (d , d ) if there exists E ⊆ E(F ) of size at most d and E ⊆ E(F ) of size at most d such that C (F , E ) C (F , E ). • We say that (d , d ) is a lower bound of the distance, denoted by m(F , F ) ≥ (d , d ), if for all AF F of F and F , the order of F satisfies |F| ≥ d + |F | = d + |F |. • If (d , d ) is an upper bound and a lower bound of the distance, we simply call it the distance between F and F , and write m(F , F ) = (d , d ).It is clear that if|F | = |F | = , the distance (d , d ) satisfies d = d = |F| -for any MAF F of F of F .The definition of m(F , F ) coincides with the one given in Equation ( . ) for a pair of trees, but with an important subtlety: the distance between a pair of forests is a pair of integers, while the distance between a pair of trees is a single integer. As we will see in the next section, the algorithm needs to work on pairs of forests with "unbalanced" numbers of components, i.e., (F , F ) with |F | ≠ |F |, hence, different numbers of cuts are needed in each forest before reaching a MAF.

  m(C (F , E ), C (F , E )) ≥ (d -|E |, d -|E |) where (d , d ) = m(F , F ). Proof. Any AF F of C (F , E ) and C (F , E ) is also an AF of F and F . So |F| ≥ d + |F | = d -|E | + |C (F , E )| and |F| ≥ d + |F | = d -|E | + |C (F , E )|. Lemma . (A cut is effective if it is used by a MAF). Let (F , F ) be a pair of X-forests with m(F , F ) = (d , d ). • (Cutting in only one forest) Let e be an edge of F . If there exists a MAF F of F and F such that e ∈ E , where E is the subset of edges of F such that F C (F , E ), then m(C (F , e ), F ) = (d -, d ). • (Cutting in both forests) Let e be an edge of F and e be an edge of F . If there exists a MAF F of F and F such that e ∈ E and e ∈ E , where E and E are the subsets of edges of F and F such that F C (F , E ) C (F , E ), then m(C (F , e ), C (F , e )) = (d -, d -). Proof. (For the first statement only; the proof of the second statement is similar) By Lemma . , m(C (F , e ), F ) ≥ (d -, d ). For the other direction, since e ∈ E , the forest F is an AF of C (F , e ) and F . We have |F| = d + |F | = (d -) + |C (F , e )|, so m(C (F , e ), F ) ≤ (d -, d ).

Figure . :

 . Figure . : Illustration of the reduction rule . Only the relevant and common structure of the two forests is drawn.

Figure . :

 . Figure . : Illustration of the reduction rule . Only the relevant and common structure of the two forests is drawn. Common subtrees have already been replaced by leaves after applying the reduction rule .

Proof.

  It suffices to consider two leaf pairs (a , a ) and (b , b ) (labeled respectively by a and b) such that p(b ) = p(p(a )) and p(b ) is the root of a component of F , and p(b ) = p(p(a )) and p(b ) is the root of a component of F (see Figure . ). Let (ab) be the label of the new leaf in (F ′ , F ′ ).

First

  , let F ′ be a MAF of F ′ and F ′ . If the leaf (ab) is isolated in F ′ , replace it by a cherry of two leaves a and b; otherwise, the parent of (ab) is the root of a component, and we can replace (ab) by a leaf b, and add a new root and connect it to the old root p(a) and to a leaf a. The result is an AF of F of F , hence, m(F , F ) ≤ m(F ′ , F ′ ).

For

  the other direction, let be a MAF of F and F We first show that the leaves labeled by a and b are in the same component of F. Assume that a and b are in different components. Then a must be an isolated vertex; b is either an isolated vertex or p(b) is the root of a component. In the first case, connect the two isolated vertices a and b to form a cherry. In the second case, add a new root to the component containing b, connecting it to the old root p(b) and to the isolated vertex a. In both cases, we get an AF of F of F of order strictly smaller than F, contradicting the minimality of the order of F. Therefore, a and b are in the same component in F. Replace the chain containing a and b by a new leaf labeled by (ab). The result is an AF of F ′ and F ′ , hence, m(F ′ , F ′ ) ≤ m(F , F ). Lemma . (The cut made by the reduction rule is effective). Let (F , F ) be a pair of X-forests. If the reduction rule cuts an edge e in F (resp. e in F ), then m

  of H , and let b be the leaf such that p(b) = p(p(a)) (such a leaf b exists because there must be a cherry (a, b) after cutting the cut edge of H ). Let a ′ and b ′ be the leaves of F having the same labels as a and b. Then exactly one of following is true: (a) (a ′ , b ′ ) is a cherry of H . In this case, (H , H ) is called a CCS of type one. (b) H has a cut edge (u, v) such that u = p(b ′ ) and p(u) = p(a ′ ). In this case, (H , H ) is called a CCS of type two.(c) H has a cut edge (u, v) such that u = p(a ′ ) and p(u) = p(b ′ ). In this case, (H , H ) is called a CCS of type three.In any of the three cases, the leaves a and b in F and the leaves a ′ and b ′ in F are said to be forming the target cherry of the CCS.

  figure should be read as follows. The two leaves forming the target cherry are labeled by a and b. The row indexed by shows the possible configurations of the leading CS of a CCS. Once a leading CS has been chosen in row , in the same column, the rows -provide possible configurations of the other CS (leading or not), such that these two CSes together form a correct CCS. Each column, indexed from to , shows the possible configurations for a given size of a CCS (starting from column , the row splits into two sub-rows). A possible configuration of a CCS is a combination of one configuration from row with one configuration from rows -, both taken from the same column. The choice between rowsdetermines whether the CCS is of type one, two, or three.

Algorithm:

  (NewCCS) Approximate m(T , T ) Input: A pair of X-trees T and T Output: An integer k such that m(T , T ) ≤ k F ← T ; F ← T ; k ← Make (F , F ) irreducible by applying the reductions rules while F and F are not isomorphic do if there exists a CCS of (F , F ) then Let (H , H ) be a best CCS of (F , F ) Let e and e be the cut edges of (H , H ) F ← C (F , e ); F ← C (F , e ); k ← k + |e | else if there exists a pre-CCS leaf pair of (F , F ) then Let (v , v ) be a best pre-CCS leaf pair of (F , F ) associated with a best CCS (H , H ) Let e and e be the pendant edges of v and v F ← C (F , e ); F ← C (F , e ); k ← k + Let e and e be the cut edges of (H , H ) F ← C (F , e ); F ← C (F , e ); k ← k + |e | else // i.e., there is no CCS and no pre-CCS leaf pair Let (v , v ) be a random leaf pair of (F , F ) Let e and e be the pendant edges of and v F ← C (F , e ); F ← C (F , e ); k ← k + end Make (F , F ) irreducible by applying the reductions rules Increment k with the number of edges the reduction rule has cut in F end return k preserve the distance, exactly k cuts are needed in either forests to arrive at an AF, therefore, k is an upper bound of the distance. The number of iterations is O(n) and each iteration requires O(n) time (to see this, notice that the number of leaf pairs is O(n) and cutting a leaf pair only creates O( ) new CCSes), hence the next lemma. Lemma . . Algorithm returns an upper bound of m(T , T ) in O(n ) time.

  . The X-coordinates correspond to the sizes of the instances. For each size, there are two points (blue circle and the red triangle) which correspond to the two instances that only differ by the ordering of the two input trees.

Figure . :Figure . :Figure . :

 ... Figure . : The approximation ratio and the running time (s) of WZ on HUGE.

FigureFigure . :Figure . :

 .. Figure . : The approximation ratio and the running time (s) of CHN on HUGE.

  are descendants of h,

	and h	w h and h = LCA (h , h ).		ℂ	d H (h, h ) + d H (h, h ) -
	Both h , h are descendants of h, and the previous case does not apply.	d H (h, h ) + d H (h, h )
	Either ( ) h is a descendant of h and h	w h,	Either ( ) d H (h, h )
	or ( ) h is a descendant of h and h	w h.	or ( ) d H (h, h )
	Any other case, that is,
	either ( ) both h	w h and h	w h,		nonvalid	nonvalid
	or ( ) h is a proper descendant of h or h .
	Table . : Definition of the event type function and the loss number function, both
	taking as input three host nodes.

  Table . and Definition . ), we can see that this desired property is already ensured (by noticing that φ(p) is an ancestor of φ(p i )), except for the arcs (p, p i ) involved in a host-switch event, that we define next:

	the cost among all time-feasible reconciliations. Unfortunately, T	-	-
		is NP-hard [Ova+ ; THL ]. For the remainder of this chap-
	Definition . (Host-switch arc). Let φ be a reconciliation of (H, P, σ). Let p ∈ V(P) be an
	internal parasite node of event type under φ, i.e., E φ (p) = , then the arc (p, p i ) ∈ A(P)
	is called the host-switch arc of p under φ, where p i is the unique child of p in P such that
	φ(p i )	w φ(p).
	The presence of the host-switch arcs can introduce incompatibilities in time, and
	as a result, a function τ H satisfying Equation ( . ) can no longer exist. Such a recon-
	ciliation, while being a feasible reconciliation in our model, is called time-unfeasible
	or cyclic (examples can be found, for instance, in [THL ]). There are two important
	observations with respect to this concept: ( ) given a reconciliation φ, whether φ is
	time-feasible can be decided in polynomial time [THL ]; ( ) for any dataset (H, P, σ),
	the set of time-feasible reconciliations is nonempty (the reconciliation that maps ev-
	ery internal parasite node to the root of the host tree does not have any host-switches).
	In order to only obtain biologically meaningful reconciliations, we can modify our
	R	into T	-	: given
	an instance, an algorithm for the latter should output a reconciliation that minimizes

  this context, a reconciliation of the dataset H, P| p , σ| L(P| p ) is simply called a reconciliation of P| p .

	The set U := V(P) × V(H) × E is called the space of cells of the DP-table. A triple
	(p, h, e) ∈ U is called a cell of the DP-table. A DP-subproblem, indexed by a cell (p, h, e),
	asks to find a reconciliation φ of P| p that minimizes the cost, and satisfies φ(p) = h
	and E φ (p) = e. The DP-table stores at each cell the optimal value of the subproblem,
	i.e., the minimum cost, denoted by M(p, h, e). The recurrence relations for computing
	M(p, h, e) are as follows:

  and e = ,

			( . )
			+∞ otherwise.
	• Otherwise, p is an internal node with children (p , p ). In this case,
	M(p, h, e) =	min e ,e ∈E h ,h ∈V(H) E(h,h ,h )=e	M(p , h , e ) + M(p , h , e ) + c(e) + c( ) ξ(h, h , h ) . ( . )

Constructing the reconciliation graph, first properties

  Definition . (Subgraph starting from a set of nodes). Let G be an ad-AND/OR graph.

	Let O be a set of OR + nodes of G. The subgraph of G starting from O, denoted by G/O, is the
	subgraph obtained from G by setting O as the new set of start nodes (i.e., by removing
	all nodes are unreachable from O through directed paths).
	We can finally describe now how to construct the reconciliation graph of a given in-
	stance of the R	from Equations ( . ) and ( . ), the recur-
	rence relations for filling up the DP-table.
	The construction is done in two steps. In the first step, we build a directed graph
	in which every node retains an additional attribute, its value, and every OR + node is
	uniquely labeled by a cell (p, h, e) ∈ U of the DP table. In the second step, we prune the
	graph by removing nodes that do not yield optimal values.
	. For each (p, h, e) ∈ U such that p is a leaf, create a goal node labeled by (p, h, e); its
	value is equal to if h = σ(p) and +∞ otherwise (see Equation ( . )). Then, for
	each (p, h, e) ∈ U in the post-order of V(P) \ L(P), let (p , p ) be the children of p,
	i. For each (p , h , e ) and each (p , h , e ) such that E(h, h , h ) = e, create an
	AND node, connect it to the two OR + nodes respectively labeled by (p , h , e )
	and (p , h , e ). Its value is equal to the sum of the values of its two children,
	plus c(e) + c( ) ξ(h, h , h ) (see Equation ( . )).

  and only if, for each triple (p, h, e) ∈ φ, there exists a unique couple (p, e) ∈ T . A set T of couples of the form (p, e) is an E-equivalence class of the instance (i.e.,

T ∈ π(T(G))) if and only if there exists φ

  is an OR node labeled by (p, e) and is equal to Merge(w , . . . , w k ), where the w i 's are the OR nodes created in the previous step. Merge are satisfied). We will prove the first part of Lemma . , i.e., Equation ( . ), by induction on the height h p of the subtree P| p . When h p = , p is a leaf and {(p, h, e)} is the only reconciliation in T(M(p, h, e)), thus, {(p, e)} is the only E-equivalence class. Now we assume h p > .

	Proof. (of Lemma . ) The second part, i.e., T(D(p, h, e)) has no redundancies, is a direct
	result of Lemma . , post-condition (it can be checked by easy induction that the
	pre-conditions of

e) is a goal node labeled by (p, e). • Otherwise, p is an internal node with children (p , p ). -For each child AND node of M(p, h, e), let (p , h , e ) and (p , h , e ) be the labels of its two child OR + nodes. Create an OR node w i labeled by (p, e) with one single child AND node. Connect this AND node to the two OR + nodes D(p , h , e ) and D(p , h , e ). -D(p, h, e) (First direction) Take T ∈ T(D(p, h, e)). It can be written as T ∪ T ∪ {(p, e)} where T i ∈ T(D(p i , h i , e i )), for i = , , for previously filled cells (p i , h i , e i ). By the induction hypotheses, T and T are valid E-equivalence classes of reconciliations of the subtrees P| p and P| p . Now, the set of colors T is a valid E-equivalence class if and only if there exists an OR + node w colored by (p, e) in the reconciliation graph that has a grandchild couple (w , w ) ∈ gch(w) of color couple ((p , e ), (p , e )) (see Definition . ). This is true; it suffices to choose w := M(p, h, e). (Second direction) Take φ ∈ T(M(p, h, e)). It can be written as φ ∪ φ ∪ {(p, h, e)} where φ i ∈ T(M(p i , h i , e i )), for i = , , for previously filled cells (p i , h i , e i ). By the induction hypotheses, there exist T ∈ T(D(p , h , e )) and T ∈ T(D(p , h , e )) such that

  To apply the methods for enumerating and counting the E-equivalence classes to the CD-equivalence relation, we only need to adapt the Definition . of the color of an OR + node so that it correctly reflects the color of a parasite node under a reconciliation for the CD-equivalence relation (seeDefinitions . and . ). Concretely, if w is an OR

+ node labeled by (p, h, e), its color is either (p, h, e), if e ≠ , or (p, ?, ) otherwise.

Table . :

 . A function MakeIntersectionFree that satisfies Lemma . . It uses the Merge function (see Lemma . ) and the InterDiff function (see Lemma . ).

Lemma . . The InterDiff function takes as input two OR + nodes w and w , such that both nodes are neat, and color P| p for some p ∈ V(P). It outputs two nodes w inter and w diff such that both nodes are neat, and

A formalism for tree-sequential dynamic programming

	able by DP or DP-like techniques [Hel ; BDI ].
	The framework of equivalence classes enumeration that we consider is situated
	within the immediate generalization of the mSDP model, i.e., generalizing finite au-
	tomata (regular expressions, paths in DAGs) to finite tree automata (regular tree gram-
	mars, solution trees of general AND/OR graphs). Further generalizations exist (from
	trees to graphs of treewidth > ); the collection of these methods is known as Non-
	serial dynamic programming [BB ].
	Since its introduction by Karp and Held [KH ], monotone sequential decision processes
	(mSDP) have been the classical model for problems solvable by dynamic programming
	(DP-problems). This formalism is based on finite-state automata. The solutions of
	DP-problems are thus equivalent to languages of regular expressions, or to paths in
	directed graphs. It is known that Bellman's principle of optimality [Bel ] also applies
	to problems for which the solutions are not sequential but tree-like [Bon ]. Various
	generalizations have been proposed to characterize broader classes of problems solv-

. .

  ). Such problems are hence naturally solvable by DP-algorithms. These algorithms are known in folklore under the name Dynamic programming on a tree. Many graph-theoretical problems (e.g., maximum matching, longest path) can be solved optimally on trees by DP-algorithms.

	Numerous real-world applications also rely on DP-algorithms on trees; examples can
	be found, for instance, in Data Science [RM ], Computer Vision [FH ; Vek ], and
	Computational Biology [BAK ; Don+ ].

  tion, recombination, and horizontal gene-transfer, are collectively known as reticulation events[HRS ]. Due to reticulation events, phylogenetic trees representing the evolutionary history of different genes found in the same set of species may differ.

	ulation events required to transform one tree to the other, which is a simple explana-
	tion for the difference between trees [BH ].
	The M	(MAF		) is closely related
	to the problem of computing the SPR distance [Hei+ ]. In the case of two rooted di-
	rected binary trees, the optimal value of the MAF	(using a subtle redefinition
	by [BS ]) coincides with the value of the rSPR distance, that is, the rooted variant of
	the SPR distance. In this chapter, we will focus on the MAF	.
	For the MAF	, several standard approaches for dealing with NP-hard prob-
	lems have been employed, including approximation algorithms, fixed-parameter al-
	gorithms [WZ ; WBZ ], and integer linear programming [Wu ]. Much of the lit-
	erature has been devoted to the search of better approximation algorithms, leading to
	a succession of algorithms with ever-improving approximation ratios. The first cor-
	rect approximation algorithm is a -approximation [Bon+ ]. Subsequently, several
	-approximations have been proposed [BMS ; RSW ; WZ ], each one improving
	on the running time. The approximation ratio was later improved to . by [Shi+ ]
	and then to / by [CMW ]. To date, the best approximation ratio of is achieved by
	[SZS ] and [Che+ ] which independently gave two algorithms using entirely differ-
	ent methods. Recently, [YCW ] made use of Monte Carlo tree search, improving on the
	-approximation of [Che+ ] to achieve a better practical result.
	Several distance metrics are commonly used to quantify the difference between
	phylogenetic trees. The Robinson-Foulds distance is a popular metric that can be cal-
	culated in linear time [RF ; Day ]. The subtree prune-and-regraft (SPR) distance
	and the hybridization number [Bar+ ] are more biologically meaningful, but are NP-
	hard to compute [AS ; BS ; BS ] (and even APX-hard [Hei+ ; Rod ]). The SPR
	distance is particularly interesting as it provides a lower bound on the number of retic-

Table . :

 . For each of the six algorithms, the running time, the average and maximum approximation ratio on the pairs of caterpillars with leaves.On this caterpillar dataset, both CMW and CHN run into some type of error. Out of the four remaining algorithms, SNV and WZ produce poorer average approximation ratios; CombMCTS gives the best approximation but is significantly slower; NewCCS is the second fastest and yields the second best average ratio. If we run the algorithms on the full dataset of instances, the estimated running time of CombMCTS is .days, while NewCCS only needs . hour. Hence, for this dataset consisting of a large number of instances of a very small size, NewCCS is of strong practical interest, as it offers a good compromise between the speed and the quality of the approximation.

		Time		Average Ratio Max Ratio	Remarks
	WZ		s	.	.
	SVN		s	.	.
	CMW	min	s	.	.	errors ( ratio > / , ratio < )
	CHN	min	s	.	.	errors (ratio < )
	CombMCTS	min	s	.	.
	NewCCS		s	.	.

Acknowledgements

Results

Dataset

Table . : The number of optimal reconciliations and the number of equivalence classes for each equivalence relation and for each instance. For numbers larger than one billion, only the order of magnitude is shown (as powers of ). For the instances with the cost vector ( , , , ), the numbers of V-and V * -equivalence classes do not coincide (see Section . . ); the latter is shown in parentheses. The * symbol indicates that the counting of the equivalence classes exceeded the memory limit of GB. . %

Event vector

Table . : The V-equivalence classes of the instance WOLB-( , , , ) and their size, as proportions of the solution space, sorted in the decreasing order of the size.

[AQE ] Nuraini Aguse, Yuanyuan Qi, and Mohammed El-Kebir. "Summarizing the solution space in tumor phylogeny inference by multiple consensus trees". In: Bioinformatics . (July ), pp. i -i . :

-.

: 10.1093/bioinformatics/btz312.

[ ----.

: 10.1007/978-3-319-42634-1_38.

[CK ] JaeJin Choi and Sung-Hou Kim. "Whole-proteome tree of life suggests a deep burst of organism diversity". In: . ( ), pp.

-. :

-. : 10.1073/pnas.1915766117.

[Con+ ] Chris Conow, Daniel Fielder, Yaniv Ovadia, and Ran Libeskind-Hadas.

"Jane: a new tool for the cophylogeny reconstruction problem". In: Algorithms for molecular biology (Feb.

), pp. -. :

-.

: 10.1186/1748-7188-5-16.

[Dar ] Charles Darwin. -.

: 10.1007/s10115-013-0625-y.

[Nar+ ] Kazuyuki Narisawa, Shunsuke Inenaga, Hideo Bannai, and Masayuki Takeda.

" -. :

-. : 10.1006/clad.1994.1010.

[PPG ] Adrian Paterson, Ricardo Palma, and Russell Gray. "Drowning on arrival, missing the boat, and x-events: how likely are sorting events". In: Tan- ), pp.

-. : -. : 10.1093/molbev/msn234.

[Ref+ ] Guislaine Refrégier, Mickaël Le Gac, Florian Jabbour, Alex Widmer, Jacqui A. Shykoff, Roxana Yockteng, Michael E. Hood, and Tatiana Giraud. "Cophylogeny of the anther smut fungi and their caryophyllaceous hosts:

Prevalence of host shifts and importance of delimiting parasite species for inferring cospeciation". In: BMC Evolutionary Biology . (Mar. ), p. . : -. : 10.1186/1471-2148-8-100.

[RE ] Elena Rivas and Sean R. Eddy. "A dynamic programming algorithm for RNA structure prediction including pseudoknots Edited by I. Tinoco".

In: Journal of Molecular Biology . ( : 10.1016/0012-365X(91)90258-4.

[RF ] D.F. ----.

: 10.1007/978-3-642-04241-6_32.

[Wu ] Yufeng Wu. "A practical method for exact computation of subtree prune and regraft distance". In: Bioinformatics . (Nov. ), pp.

-.

:

-. : 10.1093/bioinformatics/btn606.

[YCW ] Kohei Yamada, Zhi-Zhong Chen, and Lusheng Wang. : 10.1287/mnsc.17.11.712.

[Yu+ ]

Guangchuang Yu, David K. Smith, Huachen Zhu, Yi Guan, and Tommy Tsan-Yuk Lam. "ggtree: an r package for visualization and annotation of phylogenetic trees with their covariates and other associated data". In: