Guillaume Blin
email: gblin@univ-mlv.fr

Jens Stoye
email: stoye@techfak.uni-bielefeld.de

Finding Nested Common Intervals Efficiently

In this paper, we study the problem of efficiently finding gene clusters formalized by nested common intervals between two genomes represented either as permutations or as sequences. Considering permutations, we give several algorithms whose running time depends on the size of the actual output rather than the output in the worst case. Indeed, we first provide a straightforward O(n 3 ) time algorithm for finding all nested common intervals. We reduce this complexity by providing an O(n 2 ) time algorithm computing an irredundant output. Finally, we show, by providing a third algorithm, that finding only the maximal nested common intervals can be done in linear time. Considering sequences, we provide solutions (modifications of previously defined algorithms and a new algorithm) for different variants of the problem, depending on the treatment one wants to apply to duplicated genes.

Introduction and related work

Computational comparative genomics is a recent and active field of bioinformatics. One of the problems arising in this domain consists in comparing two or more species by seeking for gene clusters between their genomes. A gene cluster refers to a set of genes appearing, in spatial proximity along the chromosome, in at least two genomes. Genomes evolved from a common ancestor tend to share the same varieties of gene clusters. Therefore, they may be used for reconstructing recent evolutionary history and inferring putative functional assignments for genes.

The genome evolution process, including -among others -fundamental evolutionary events such as gene duplication and loss [START_REF] Ohno | Evolution by gene duplication[END_REF], has given rise to various genome models and cluster definitions. Indeed, genomes may be either represented as permutations (allowing one-to-one correspondence between genes of different genomes) or sequences -where the same letter (i.e. gene) may occur more than once (a more realistic model but with higher complexity). In both those models, there may exist, or not, genes not shared between two genomes (often called gaps).

Moreover, when modeling genomes for gene order analysis, one may consider either two or multiple genomes, seeking for exact or approximate occurrences, finding all or just non-extensible (i.e. maximal ) occurrences.

There are numerous ways of mathematical formalizations of gene clusters. Among others, one can mention common substrings (which require a full conservation), common intervals [START_REF] Bergeron | Computing common intervals of k permutations, with applications to modular decomposition of graphs[END_REF][START_REF] Didier | Character sets of strings[END_REF][START_REF] Schmidt | Quadratic time algorithms for finding common intervals in two and more sequences[END_REF][START_REF] Uno | Fast algorithms to enumerate all common intervals of two permutations[END_REF] (genes must occur consecutively, regardless of their order), conserved intervals [START_REF] Bergeron | On the similarity of sets of permutations and its applications to genome comparison[END_REF] (common intervals, framed by the same two genes), gene teams [START_REF] Zhang | Gene team tree: A compact representation of all gene teams[END_REF][START_REF] Bergeron | The algorithmic of gene teams[END_REF][START_REF] He | Identifying conserved gene clusters in the presence of homology families[END_REF] (genes in a cluster must not be interrupted by long stretches of genes not belonging to the cluster), and approximate common intervals [START_REF] Böcker | Computation of median gene clusters[END_REF][START_REF] Rahmann | Integer linear programs for discovering approximate gene clusters[END_REF] (common intervals that may contain few genes from outside the cluster). For more details, please refer to [START_REF] Bergeron | Formal models of gene clusters[END_REF].

In this article, we focus on another model -namely the nested common intervals -which was mentioned in [START_REF] Hoberman | The incompatible desiderata of gene cluster properties[END_REF]. In this model, an additional constraint -namely the nestedness -(observed in real data [START_REF] Kurzik-Dumke | Identification of a novel Drosophila melanogaster gene, angel, a member of a nested gene cluster at locus 59F4,5[END_REF]) is added to the cluster definition. Hoberman and Durand [START_REF] Hoberman | The incompatible desiderata of gene cluster properties[END_REF] argued that, depending on the dataset, if the nestedness assumption is not excluding clusters from the data, then it can strengthen the significance of detected clusters since it reduces the probability of observing them by chance.

As far as we know, [START_REF] Hoberman | The incompatible desiderata of gene cluster properties[END_REF] was the only attempt to take into account the nestedness assumption in a gene cluster model (namely gene teams) and yields to a quadratic-time greedy bottom-up algorithm. In fact, no explicit algorithmic analysis is given in [START_REF] Hoberman | The incompatible desiderata of gene cluster properties[END_REF], which might be a bit dangerous in view of the fact that genomes may consist of up to 25, 000 and more genes. In the following, we will give some efficient algorithms to find all gene clusters -considering nested common intervals -between two genomes.

Let π 1 and π 2 be our genomes, represented as permutations over N := {1, . . . , n}. For any i ≤ j, π[i, j] will refer to the sequence of elements

(π[i], π[i + 1], . . . , π[j]). Let CS(π[i, j]) := {π[k] | k ∈ [i, j]} denote the character set of the interval [i, j] of π. A subset C ⊆ N is called a common interval of π 1 and π 2 if and only if there exist 1 ≤ i 1 < j 1 ≤ n and 1 ≤ i 2 < j 2 ≤ n such that C = CS(π 1 [i 1 , j 1 ]) = CS(π 2 [i 2 , j 2 ]
). Note that this definition purposely excludes common intervals of size one since they would not be considered in the more general nested common interval definition. The intervals [i 1 , j 1 ] and [i 2 , j 2 ] are called the locations of C in π 1 and π 2 , respectively.

Given two common intervals C and

C ′ of π 1 and π 2 , C contains C ′ if and only if C ′ ⊆ C. This implies that the location of C ′ in π 1 (resp. π 2 ) is included in the location of C in π 1 (resp. π 2 ). A common interval C is called a nested common interval of π 1 and π 2 if either |C| = 2, or if |C| > 2 and it contains a nested common interval of size |C| -1.
Note that this recursive definition implies that for any nested common interval C there exists a series of nested common intervals such that

C 2 ⊆ C 3 ⊆ • • • ⊆ C with |C i | = i.
A nested common interval of size ℓ is maximal if it is not contained in a nested common interval of size ℓ + 1. A maximal nested common interval can however still be contained in a larger nested common interval. For example, considering π 1 := (3, 1, 2, 4, 5, 6) and π 2 := (1, 2, 3, 4, 5, 6), the maximal nested common interval [START_REF] Bergeron | On the similarity of sets of permutations and its applications to genome comparison[END_REF][START_REF] Böcker | Computation of median gene clusters[END_REF] in π 1 is contained in [START_REF] Bergeron | Computing common intervals of k permutations, with applications to modular decomposition of graphs[END_REF][START_REF] Böcker | Computation of median gene clusters[END_REF].

The general Nested Common Intervals problem may be defined as follows: Given two genomes, find all their nested common intervals. One can then consider genomes either as permutations or sequences and might also be interested in finding only the maximal nested common intervals and/or allowing gaps. In the two following sections, we will give efficient algorithms for both permutations and sequences but will leave the case considering gaps as an open problem.

In [START_REF] Bergeron | Computing common intervals of k permutations, with applications to modular decomposition of graphs[END_REF], Bergeron et al. proposed a theoretical framework for computing common intervals based on a linear space basis. Of importance here is the technique proposed in [START_REF] Bergeron | Computing common intervals of k permutations, with applications to modular decomposition of graphs[END_REF] in order to generate the PQ-tree [START_REF] Landau | Using pq trees for comparative genomics[END_REF] corresponding to a linear space basis for computing all the common intervals of K permutations. Generating this basis can be done in O(n) time for two permutations of size n. Then one can, by a browsing of the tree, generate all the common intervals in O(n + z) time where z is the size of the output. One can adapt this algorithm in order to find nested common intervals in O(n + z) time.

In this work, we did not follow that approach, since (1) PQ-trees are a heavy machinery and quite space consuming in practice, and (2) our aim was to provide easy-to-implement algorithms with small constants in the O-notation. Moreover, algorithms based on PQ-trees do not easily generalize to sequences.

Nested common intervals on permutations

As described in [START_REF] Didier | Character sets of strings[END_REF][START_REF] Schmidt | Quadratic time algorithms for finding common intervals in two and more sequences[END_REF][START_REF] Uno | Fast algorithms to enumerate all common intervals of two permutations[END_REF], when considering permutations, both common substrings and common intervals can be found in optimal, essentially linear time. As we will show, not surprisingly, finding nested common intervals on permutations can also be done efficiently.

Finding all nested common intervals

First, one has to notice that the number of nested common intervals can be quadratic in n (e.g. when π 1 = π 2 ). However, in many practical cases the number of nested common intervals may be much smaller, such that one can still achieve lower time complexity by developing methods whose running time depends on the size of the actual output and not of the output in the worst case.

In the following, we will w.l.o.g. assume that π 1 is the identity permutation and rename π 2 by π for ease of notation. A naive bottom-up algorithm, inspired from the one given in [START_REF] Hoberman | The incompatible desiderata of gene cluster properties[END_REF] and straightforwardly following the definition of nested common intervals, can be defined as in Algorithm 1.

Algorithm 1 Find all nested common intervals

1: for i ← 1, . . . , n do 2: l ← i, r ← i 3: repeat 4: l ′ ← l, r ′ ← r 5: if π[l ′ -1] = min(CS(π[l ′ , r ′ ])) -1 or π[r ′ + 1] = max(CS(π[l ′ , r ′ ])) + 1 then 6: while π[l -1] = min(CS(π[l, r ′ ])) -1 do l--done 7: while π[r + 1] = max(CS(π[l ′ , r])) + 1 do r++ done 8: else 9: while π[l -1] = max(CS(π[l, r ′ ])) + 1 do l--done 10: while π[r + 1] = min(CS(π[l ′ , r])) -1 do r++ done 11: end if 12: report all intervals [l ′′ , r ′′ ] with l ≤ l ′′ ≤ l ′ and r ′ ≤ r ′′ ≤ r except [l ′ , r ′ ] 13:
until l = l ′ and r = r ′ 14: end for Clearly such an algorithm requires O(n + z) time to report all nested common intervals where z is the size of the output. However, the output may be highly redundant as several intervals will be identified more than once. The worst case is when one considers π = (1, 2, . . . , n). More precisely, in this case some of the O(n 2 ) nested common intervals will be reported up to n times, giving a total worst-case runtime of O(n 3 ).

Therefore, one may be interested in computing an irredundant output. The main improvement we propose consists in a simple preprocessing step that will speed up our algorithm for nested common intervals. Let us define a run of two permutations π 1 and π 2 as a pair of intervals

([i 1 , j 1 ], [i 2 , j 2 ]) such that π 1 [i 1 , j 1 ] = π 2 [i 2 , j 2 ] or π 1 [i 1 , j 1 ] = ←----- π 2 [i 2 , j 2 ] where ← - x := (x k , x k-1 , . . . , x 1
) denotes the reverse of sequence x = (x 1 , x 2 , . . . , x k ). A run is maximal if it cannot be extended to the left or right. Since a run can also be of size one, two permutations can always be decomposed into their maximal runs with respect to each other. For example, in the following the maximal runs are underlined: π 1 = (1, 2, 3, 4, 5, 6, 7, 8, 9) and π 2 = (4, 3, 2, 1, 5, 9, 6, 7, 8).

Given a decomposition of two permutations into their maximal runs with respect to each other, a breakpoint will refer to any pair of neighboring elements that belong to different runs. In the above example, π 1 (resp. π 2 ) contains three breakpoints (4, 5), [START_REF] Blin | Comparing genomes with duplications: a computational complexity point of view[END_REF][START_REF] Böcker | Computation of median gene clusters[END_REF] and (8, 9) (resp. [START_REF] Bergeron | Computing common intervals of k permutations, with applications to modular decomposition of graphs[END_REF][START_REF] Blin | Comparing genomes with duplications: a computational complexity point of view[END_REF], [START_REF] Blin | Comparing genomes with duplications: a computational complexity point of view[END_REF][START_REF] Hoberman | The incompatible desiderata of gene cluster properties[END_REF] and [START_REF] Hoberman | The incompatible desiderata of gene cluster properties[END_REF][START_REF] Böcker | Computation of median gene clusters[END_REF]). By definition, the number of breakpoints is one less than the number of runs. When considering one of π 1 and π 2 as being the identity permutation then a run may be defined as a single interval.

Algorithm 2, hereafter defined, computes irredundant output by making use of the two following simple observations. Lemma 1. All subintervals of length at least 2 in a run are nested common intervals.

Lemma 2. In the procedure of constructing incrementally a nested common interval by extension, when one reaches up to one end of a run such that the begin/end of this run can be included, then the whole run can be included and all subintervals ending/beginning in this run can be reported as nested common intervals.

Proof. By definition, the elements of a run [i, j] in π with respect to the identity permutation are consecutive integers strictly increasing or decreasing. Therefore, in an incremental construction by extension of a nested common interval nc that has a run [i, j] as its right (resp. left) neighbor, if one may extend nc by i (resp. j) then, by definition, π[i] (resp. π[j]) is the minimal or maximal element among the elements of the extended interval. Thus, all the elements of the run [i, j] may be added one by one, each leading to a new nested common interval.

⊓ ⊔

Consequently, by identifying the runs in a preprocessing step (which can be easily done in linear time), whenever during an extension a border of a run is included, the whole run is added at once and all sub-intervals are reported. The details are given in Algorithm 2 which employs a data structure end defined as follows: if i is the index of the first or last element of a run in π then end[i] is the index of the other end of that run; end[i] = 0 otherwise.

It is easily seen that the time complexity remains in O(n + z) with z being the output size, except that this time the output is irredundant (i.e. each nested common interval is reported exactly once). When applied to π := (1, 2, 3, 6, 4, 5) for example, one may check that Algorithm 2 will report locations [START_REF] Bergeron | Computing common intervals of k permutations, with applications to modular decomposition of graphs[END_REF][START_REF] Bergeron | The algorithmic of gene teams[END_REF], [START_REF] Bergeron | Computing common intervals of k permutations, with applications to modular decomposition of graphs[END_REF][START_REF] Bergeron | Formal models of gene clusters[END_REF], [START_REF] Bergeron | The algorithmic of gene teams[END_REF][START_REF] Bergeron | Formal models of gene clusters[END_REF] when i = 1, locations [START_REF] Blin | Comparing genomes with duplications: a computational complexity point of view[END_REF][START_REF] Böcker | Computation of median gene clusters[END_REF], [START_REF] Bergeron | On the similarity of sets of permutations and its applications to genome comparison[END_REF][START_REF] Böcker | Computation of median gene clusters[END_REF], [START_REF] Bergeron | Computing common intervals of k permutations, with applications to modular decomposition of graphs[END_REF][START_REF] Böcker | Computation of median gene clusters[END_REF],

Algorithm 2 Find all nested common intervals, irredundant version 1: decompose π into maximal runs w.r.t. id and store them in a data structure end 2: for i ← 1, . . . , n -1 do 3:

if

end[i] > i then 4: l ← i, r ← end[i] 5:
report all intervals [l ′′ , r ′′ ] with l ≤ l ′′ < r ′′ ≤ r 6: repeat 7:

l ′ ← l, r ′ ← r 8: if π[l ′ -1] = min(CS(π[l ′ , r ′ ]))-1 or π[r ′ +1] = max(CS(π[l ′ , r ′ ]))+1 then 9: if π[l -1] = min(CS(π[l, r ′ ])) -1 then l ← end[l -1] end if 10: if π[r + 1] = max(CS(π[l ′ , r])) + 1 then r ← end[r + 1] end if 11: else 12: if π[l -1] = max(CS(π[l, r ′ ])) + 1 then l ← end[l -1] end if 13: if π[r + 1] = min(CS(π[l ′ , r])) -1 then r ← end[r + 1] end if 14: end if 15: report all intervals [l ′′ , r ′′ ] with l ≤ l ′′ ≤ l ′ and r ′ ≤ r ′′ ≤ r except [l ′ , r ′ ] 16:
until l = l ′ and r = r ′ 17:

end if 18: end for [START_REF] Bergeron | The algorithmic of gene teams[END_REF][START_REF] Böcker | Computation of median gene clusters[END_REF], [START_REF] Bergeron | Formal models of gene clusters[END_REF][START_REF] Böcker | Computation of median gene clusters[END_REF] when i = 5, and nothing for the other values of i. Let us prove this interesting property of Algorithm 2. Proposition 1. In Algorithm 2, any breakpoint is considered (i.e. passed through) at most once during the extension procedure described from lines 8 to 14.

Proof. Assume bp = (π[x], π[y]

) is a breakpoint in π with respect to the identity permutation. Then y = x + 1 and there are only two possibilities for bp to be passed through twice: either (1) once from the left and once from the right, or (2) twice from the same side.

Let us first consider the case where bp is passed through from both left and right. Therefore assume that

π = (. . . , π[X], . . . , π[x], π[y], . . . , π[Y ], . . .)
where C 1 := CS(π[X, x]) and C 2 := CS(π[y, Y ]) are nested common intervals (otherwise bp would not be passed through more than once). Further assume w.l.o.g. that bp is passed in the extension of interval [X, x] (a similar proof can be easily provided for the extension of [y, Y ]) whose maximum (or minimum) element is M := max(C 1 ) (resp. m := min(C 1 )).

Then, in order for [X, x] to be extensible, either π[y] = M + 1 or π[y] = m-1. Let us assume that π[y] = M +1 (the case where π[y] = m-1 can be shown similarly). We will show that bp cannot be passed through in the other direction, i.e. in an extension of [y, Y ]. Since π[y] = M + 1 and each of the two intervals [X, x] and [y, Y ] consists of consecutive integers, we have that all elements in C 1 are smaller than any element in C 2 . Thus, for an extension in the left direction across bp, the largest element of C 1 must be at its right end, i.e. π[x] = M . However, if this was the case, then bp = (π[x], π[y]) = (M, M + 1) would not be a breakpoint, a contradiction. Now let us consider the case where bp is passed through twice from the same side, starting with different runs. Therefore assume that

π = (. . . , π[X ′ ], . . . , π[X], . . . , π[x], π[y], . . .)
where C 1 := CS(π[X, x]) and C 2 := CS(π[X ′ , x]) are nested common intervals derived from different runs (i.e. reported from two different values of i), one in the interval [X, x] and the other in the interval [X ′ , X -1]. Further assume w.l.o.g. that bp is passed through in the extension of [X, x] (a similar proof can be easily provided for the extension of [X ′ , x]) whose maximum (or minimum) element is M := max(C 1 ) (resp. m := min(C 1 )).

Then, in order for [X, x] to be extensible, either π[y] = M + 1 or π[y] = m-1. Let us assume that π[y] = M +1 (the case where π[y] = m-1 can be shown similarly). We will show that bp cannot be passed again in this direction, i.e. in an extension of [X ′ , x]. Since π[y] = M + 1 and, by construction, C 1 ⊂ C 2 , we have that all elements in CS(π[X ′ , X -1]) are smaller than any element in C 1 . Moreover, since bp is a breakpoint, we have that π[x] = M and, more precisely, π[x] < M . Then, any extension of [X ′ , X -1] would not be able to include M since at least one necessary intermediate element (namely π[x]) would not have been previously included. Thus, all cases are covered and the proposition is proved.

⊓ ⊔

Irredundancy of the locations of nested common intervals returned by Algorithm 2 follows immediately. Proposition 2. In Algorithm 2, two different runs cannot yield the same nested common interval.

Proof. In order to be possibly reported twice, an interval would have to be a superinterval of two different runs. However, in order to yield the interval, the breakpoint(s) between these two runs would have to be passed through twice, which is not possible by Proposition 1.

⊓ ⊔

Finding all maximal nested common intervals

As previously mentioned, one might also be interested in finding only the maximal nested common intervals in optimal time O(n+z) where z is the number of maximal nested common intervals of π 1 and π 2 , since there will be fewer. In fact, we will first prove that the number of maximal nested common intervals is in O(n) leading to an overall linear time algorithm.

Proposition 3. Every element of π is a member of at most three different maximal nested common intervals.

Proof. This follows immediately from the correctness of Proposition 1. Indeed, according to Proposition 1 each position can be reached from at most two directions. Thus, the only case where an element of π may be member of exactly three different maximal nested common intervals

nc 1 = [i 1 , j 1 ], nc 2 = [i 2 , j 2 ] and nc 3 = [i 3 , j 3 ] is when i 1 ≤ i 2 ≤ i 3 ≤ j 1 ≤ j 2 ≤ j 3 .
For example, considering π 1 = (2, 1, 3, 4, 6, 5) and π 2 = (1, 2, 3, 4, 5, 6), the element 3 in π 1 is member of (2, 1, 3, 4), [START_REF] Bergeron | Formal models of gene clusters[END_REF][START_REF] Bergeron | On the similarity of sets of permutations and its applications to genome comparison[END_REF] and [START_REF] Bergeron | Formal models of gene clusters[END_REF][START_REF] Bergeron | On the similarity of sets of permutations and its applications to genome comparison[END_REF][START_REF] Böcker | Computation of median gene clusters[END_REF][START_REF] Blin | Comparing genomes with duplications: a computational complexity point of view[END_REF].

⊓ ⊔

In order to get only the locations of maximal nested common intervals, one has to modify Algorithm 2 such that only the locations at the end of an extension are reported. To do so, one has simply to (1) remove from Algorithm 2 lines 5 and 15 and (2) report the unique interval [l, r] -which is by definition maximal -just after the end of the repeat ... until loop (currently line 16).

Clearly, the time complexity of this slightly modified version of Algorithm 2 is unchanged; that is O(n + z) where z is the size of the output. Proposition 3 implies that the number of maximal nested common intervals is in O(n), leading to an overall linear time.

Nested common intervals on sequences

In this section, we will give algorithms (mainly ideas due to space constraints) to handle genomes represented as sequences (i.e. genes may be duplicated). In the following, we will assume that our genomes, denoted by S 1 and S 2 , are defined over a bounded integer alphabet Σ = {1, . . . , σ} and have maximal length n. The precise definition of nestedness in sequences is subtle. Therefore, we propose three different variants of the problem, depending on the treatment one wants to apply when, during the extension of an interval, an element that is already inside the interval is met once again.

First, one may just extend the interval "for free", only caring about the "innermost occurrence"; all other occurrences are considered as not contributing to the cluster content. This definition follows the same logic as earlier ones used for common intervals [START_REF] Didier | Character sets of strings[END_REF][START_REF] Schmidt | Quadratic time algorithms for finding common intervals in two and more sequences[END_REF] and for approximate common intervals [START_REF] Böcker | Computation of median gene clusters[END_REF].

A slight modification of our naive Algorithm 1 leads to an O(n 3 ) algorithm. Indeed, since the sequences may contain duplicates, one has to start the extension procedure with all possible pairs (S 1 [i], S 2 [j]) where 1 ≤ i ≤ |S 1 | and 1 ≤ j ≤ |S 2 |. Moreover, after each extension step all genes that are already members of the cluster have to be "freely" included. This can be tested efficiently by storing the elements belonging to the current cluster in a bit vector c[1, . . . , σ].

The resulting Algorithm 3 -which clearly runs in O(n 3 ) time as each pair of index positions (i, j) is considered at most once and for each of them the extension cannot include more than n steps -only reports maximal gene clusters as previously done for permutations.

Second, one may forbid the inclusion of a second copy of a gene in a nested common interval. This problem can also be solved easily, by a quite similar algorithm which stops any extension when a gene already contained in the common interval is encountered.

Finally, one may be interested in finding a bijection (sometimes called matching in the computational comparative genomics literature) where, inside a nested common interval, each gene occurrence in S 1 must match a unique gene occurrence in S 2 from the same gene family. Surprisingly, the nestedness constraint leads to a polynomial time algorithm whereas for many other paradigms, considering matching and duplicates leads to hardness [START_REF] Blin | Comparing genomes with duplications: a computational complexity point of view[END_REF].

For this last variant, we unfortunately only know a very inefficient algorithm described hereafter. The main idea is to, first, construct a directed acyclic graph G whose vertices correspond to pairs of intervals ([i 1 , j 1 ], [i 2 , j 2 ]), one from S 1 and one from S 2 , that contain the same multiset of characters. In G, an edge is drawn from a

vertex v = ([i 1 , j 1 ], [i 2 , j 2 ]) to a vertex v ′ = ([i ′ 1 , j ′ 1 ], [i ′ 2 , j ′ 2 ]
) iff the corresponding interval pairs differ by one in length -i.e. |(j 1i 1 ) -(j ′ 1i ′ 1 )| = 1 -and the shorter one is contained in the longer one, i.e. ((i 1 = i ′ 1 ) or (j 1 = j ′ 1 )) and ((i

2 = i ′ 2 ) or (j 2 = j ′ 2 )
). An illustration is given in Figure 1. Since, for a given multiset of cardinality ℓ there are at most (n-ℓ+1) 2 vertices in the graph, the total number of vertices in G is bounded by O(n 3 ). Moreover, by definition, each vertex has an output degree of at most four, hence the number of edges is also bounded by O(n 3 ). Finally, G can clearly be constructed in polynomial time. One can easily see that there is a correspondence between nested gene clusters and directed paths in G starting from vertices corresponding to multisets of size 2. Indeed, a path (ci 1 , ci 2 , . . . ci k ) in this DAG, where ci 1 , ci 2 , . . . ci k are common intervals, induces that ci 1 ⊆ ci 2 ⊆ . . . ⊆ ci k and ∀1 ≤ j < k we have Algorithm 3 Find all maximal nested common intervals in two sequences

1: for i ← 1, . . . , |S1| do 2:
for each occurrence j of S1

[i] in S2 do 3: for each k ← 1, . . . , σ do c[k] ← (k = S1[i]) done 4: l1 ← i, r1 ← i 5: l2 ← j, r2 ← j 6: repeat 7: while c[S1[l1 -1]] = true do l1 ← l1 -1 done 8: while c[S1[r1 + 1]] = true do r1 ← r1 + 1 done 9: while c[S2[l2 -1]] = true do l2 ← l2 -1 done 10: while c[S2[r2 + 1]] = true do r2 ← r2 + 1 done 11: l ′ 1 ← l1, r ′ 1 ← r1 12: if S1[l1 -1] = S2[l2 -1] or S1[r1 + 1] = S2[r2 + 1] then 13: while S1[l1 -1] = S2[l2 -1] do 14: l1--, l2--, c[S1[l1]] ← true 15: while c[S1[l1 -1]] = true do l1--done 16: while c[S2[l2 -1]] = true do l2--done 17:
end while 18:

while |ci j | + 1 = |ci j+1 |. Therefore, since any common interval of size 2 is a nested common interval, in any such path, if ci 1 is of size 2 then, by definition, any common interval of this path is a nested common interval. Thus, the nested gene clusters can be reported in polynomial time. Indeed, building the DAG can be done in O(n 3 ) time. Then, one has to browse any path starting from a vertex corresponding to a common interval of size 2. Since, there are at most n -1 such common intervals, and each vertex has an output degree of at most four, on the whole the number of such paths is bounded by (n -1) • 4(n -2); i.e. O(n 2 ).

S1[r1 + 1] = S2[r2 + 1] do 19: r1++, r2++, c[S1[r1]] ←

Conclusion

In this article, we proposed a set of efficient algorithms considering the nestedness assumption in the common intervals model of gene clusters for genomes represented both as permutations and as sequences. Two main questions remain open: (1) finding a more efficient algorithm for the last variant of nested common intervals on sequences and (2) allowing clusters to be interrupted by up to g consecutive genes from outside the cluster, as in [START_REF] Hoberman | The incompatible desiderata of gene cluster properties[END_REF].

(Fig. 1 .

 1 Fig. 1. Graph G for sequences S1 = (1, 3, 2, 3, 4, 5, 6) and S2 = (2, 3, 1, 3, 2, 4, 6, 3, 5).

Acknowledgments

The authors wish to thank Ferdinando Cicalese, Martin Milanič and Mathieu Raffinot for helpful suggestions on the complexity of finding nested common intervals on sequences and on PQ-trees aspects.


Michel Chilowicz

Etienne Duris

Gilles Roussel

Towards a multi-scale approach for source code approximate match report

Keywords: D.2.8 [Software Engineering]: Metrics; I.5.3 [Computing Methodologies]: Pattern Recognition Source code similarity, clones, software plagiarism

Finding exact clones in source code can be efficiently handled using classical exact substring or subtree pattern matching techniques inspired from genomics applications. These methods may be wisely employed as a foundation to sketch new techniques highlighting duplicated code chunks presenting minor edits or more extensive modifications at a higher structural scale. The main goal is to improve recall of small near matches and to aggregate them into larger ones to provide a more global view of similarities with a reasonable complexity. These concerns are essential to be able to address a large database of source code projects.

Most popular clone detection tools, either in a code reengineering or in a plagiarism highlighting context, consider sequences of tokens [START_REF] Prechelt | Finding plagiarisms among a set of programs with jplag[END_REF][START_REF] Kamiya | Ccfinder: A multilinguistic token-based code clone detection system for large scale source code[END_REF] or syntax trees [START_REF] Baxter | Clone detection using abstract syntax trees[END_REF] as an abstracted view of the source code. For token sequence representations, groups of occurrences of exact repeated token factors can be found with suffix indexing structures like the suffix arrays [START_REF] Manber | Suffix arrays: a new method for on-line string searches[END_REF] or the suffix trees [START_REF] Weiner | Linear pattern matching algorithm[END_REF] in linear time. Zones of high density of similar n-grams [START_REF] Schleimer | Winnowing: Local algorithms for document fingerprinting[END_REF] can also be used to spot similarities. Nevertheless, when local edits are introduced between code clones, looking for exact repeated factors is inappropriate: in this case, local alignment algorithms [START_REF] Irving | Plagiarism and collusion detection using the Smith-Waterman algorithm[END_REF] could be used to retrieve approximate substring matches containing small unmatched gaps. However, these techniques require a detrimental quadratic running cost in number of tokens: it is preferable to use them on carefully selected zones containing suspected clones.

FUNCTIONAL CALL GRAPHS OF TOKEN SEQUENCES

To cope with multiple occurrences of similarities and to have a global view of them, we proposed in [START_REF] Chilowicz | Finding similarities in source code through factorization[END_REF] a method allowing similarity report at a kind of function level. Using suffix structures, we split each original function of the source code into different components that are either token substrings shared by other functions or unmatched token sequences. Shared token substrings may themselves be decomposed with the help of smaller nested similarities found among them. Following this approach, we transform individual function call graphs of the compared projects into a global call graph introducing new synthetic functions that represent the granularity of shared chunks of code across the projects. We do not use these call graphs for refactoring purpose, but rather to define several metrics based on the amount of code represented by shared nodes between projects. This method allows similarity retrieval at a synthetic-function level despite extensive edits involving removal, transposition, addition of source code or even inlining or outlining of functions (clones of type 3 according to [START_REF] Bellon | Comparison and evaluation of clone detection tools[END_REF] or types 3-4 for [START_REF] Roy | Scenario-based comparison of clone detection techniques[END_REF]).

Thanks to this factorization approach, comparing elementary unmatched token sequences using an alignment algorithm allows retrieval of similarities that would be missed otherwise thanks to resiliency towards local edits (like removing or adding neutral token sequences, for instance -0 or * 1 in arithmetic expressions like in figure 1). Identifying tiny local similarities may also help to counter to the development of identifiers into expressions (like the development of identifier diff seen in lines 8-9 of stdDev2) but also to select pairs of token sequences containing high densities of tiny clones that will be tested through alignment methods.

CLASSIFYING ABSTRACT SYNTAX SUBTREES

Concerning methods based on syntax tree representations, most accurate approaches [START_REF] Zhang | Simple fast algorithms for the editing distance between trees and related problems[END_REF] involve extensive and very costly dynamic programming comparison of all subtrees comparing to linear token approaches. Preclustering potential similar subtrees, using degraded hashing strategies or metrics [START_REF] Jiang | Deckard: Scalable and accurate tree-based detection of code clones[END_REF], appears essential to reduce the number of subtree We explore in [START_REF] Chilowicz | Syntax tree fingerprinting: a foundation for source code similarity detection[END_REF] some of these hashing strategies based on an abstracted view of the syntax tree associated with different hash functions. We could then introduce hash values of subtrees according to several abstraction profiles (confounding all types, discarding small subtrees above a given size, considering commutativity of operators...) and infer similarity metrics on subtree pairs according to these abstraction profiles. Thus, duplicated subtrees involving most trivial local edits are detected and their distance quantified.

STRUCTURAL CONSOLIDATION OF LOCAL CLONES

The second challenge related to clone detection based on syntax trees focuses on merging simple subtree clones to form higher level clones. This problem does not involve only gathering sibling subtree clones that can be handled with suffix indexation algorithms [START_REF] Chilowicz | Syntax tree fingerprinting for source code similarity detection[END_REF] but also more distant clone subtrees like cousins. It allows detection of large clones composed of smaller ones that are transposed or flooded with other clones that are too heavily transformed to be recognized.

Gathering close clones in the syntax tree using various fast heuristics is an interesting issue to be explored. The call graph can also be used to assemble clones that are potentially reachable from a common function. It leads then to a more concise view of the similarities in source code at a higher level that may be later refined according to the user will. For example, in figure 1, reporting a single match between the complete functions stdDev1 and its copied counterpart stdDev2 appears more appropriate that two matches involving mean and stdDev1 respectively linked to the start and the end of stdDev2.

Figure 1 :

 1 Figure 1: An original standard deviation computation function stdDev1 in Java and a transformed copy stdDev2 involving inlining, identifier renaming and local edits


Frédéric Fauberteau
email: fauberte@univ-mlv.fr

Serge Midonnet
email: midonnet@univ-mlv.fr

Laurent George
email: lgeorge@ieee.org

Robust Partitioned Scheduling for Real-Time Multiprocessor Systems

Keywords: Fit}-Decreasing partitioning Real-time Scheduling, Partitioned Scheduling, Robustness

In this paper, we consider the problem of fixed-priority partitioned scheduling of sporadic real-time tasks for homogeneous processors. We propose a partitioning heuristic that takes into account possible Worst Case Execution Time (WCET) overruns. Our goal is to maximize the duration a task can be allowed to exceed its WCET without compromising the timeliness constraints of all the tasks. This duration is denoted in the paper the allowance of the task and is computed with a sensitivity analysis. The partitioning heuristic we propose, assigns the tasks to the processors in order (i) to maximize the allowance of the tasks and (ii) to tolerate bounded execution duration overruns. Property (ii) is important as real-time applications are often prone to be subject to OS approximations or software faults that might result in execution duration overruns. We show with performance evaluations that Allowance-Fit-Decreasing partitioning improves the temporal robustness of real-time systems w.r.t. classical {First-Fit/Best-Fit/Next

Introduction

Fixed-priority scheduling of recurring real-time tasks has been largely studied for uniprocessors. In such a scheduling, a Priority Assignment (PA) assigns a fixed priority to each job of the task. For instance, Rate-Monotonic (RM) is an optimal PA for periodic tasks with implicit-deadlines (deadlines equal to periods) [START_REF] Liu | Scheduling algorithms for multiprogramming in a hard real-time environment[END_REF]. Optimality implies that if a feasible PA over a taskset exists, then the optimal PA is also feasible. A feasible taskset is a taskset such that a scheduling algorithm exists which can schedule this taskset. We focus on the more general model of tasks with constrained-deadlines (deadlines less than or equal to periods) for which Deadline-Monotonic (DM) is an optimal PA [START_REF] Leung | On the complexity of fixed-priority scheduling of periodic, real-time tasks[END_REF]. Recently, the optimal Robust Priority Assignment (RPA) [START_REF] Davis | Robust priority assignment for fixed priority real-time systems[END_REF] has been proposed to find the PA which maximizes the interference that a tasks system can support. These interferences can be handled by the tasks by allowing WCETs overruns while the timeliness constraints of all the tasks are respected. These tolerated WCETs overruns are denoted allowance of the tasks. In the same way, our motivation is to propose a robust multiprocessor scheduling which maximizes the allowance.

The two most studied approaches to schedule real-time tasks on a multiprocessor are partitioned and global scheduling. The first one does not allow tasks to migrate whereas the second one allows unrestricted migrations. Recent architecture have reduced the cost of migration. Nevertheless, taking into account the cost of migration in the feasibility conditions of global scheduling is still an open issue. A recent performance evaluation of partitioned and global schedulings show that partitioned scheduling outperform global scheduling, in the current state-of-the-art of feasibility conditions [START_REF] Bertogna | Evaluation of existing schedulability tests for global EDF[END_REF]. We therefore focus on the partitioned approach. Several algorithms for fixed-priority partitioned scheduling have been proposed [START_REF] Burchard | New strategies for assigning realtime tasks to multiprocessor systems[END_REF][START_REF] Oh | Allocating fixed-priority periodic tasks on multiprocessor systems[END_REF][START_REF] Andersson | Preemptive multiprocessor scheduling anomalies[END_REF][START_REF] Fisher | The partitioned scheduling of sporadic tasks according to static-priorities[END_REF][START_REF] Lakshmanan | Partitioned fixed-priority preemptive scheduling for multi-core processors[END_REF]. The aim of the authors is to propose algorithms which improve the worst-case utilization bound. The worst-case utilization bound for a scheduling algorithm A is defined as the minimum utilization for which any implicit-deadline taskset is schedulable according to algorithm A. The utilization of an implicit-deadline taskset is the sum of the processor utilization (formally defined in Sect.2) of each task composing this taskset.

In this paper, our motivation is slightly different since we want to design a partitioned scheduling which improves the temporal robustness of a system i.e. to improve its capability to support variations on the system parameters at run time (WCET overruns for e.g.). Such events should not lead to a deadline violation in a hard real-time application. We focus on the WCET parameter and we propose an algorithm which allocates the tasks on the processor having the greatest capability to support WCETs overrun by maximizing the minimum allowance of all the tasks.

The rest of this paper is organized as follows. In Section 2, we introduce the terminology used in the rest of this paper. In Section 3, we give a definition of robustness in context of this paper. In Section 4, we discuss two manners to compute the allowance of the execution duration which is the criterion of our partitioning algorithm for the assignment of real-time tasks on the processors. In Section 5, we present our heuristic and describe how it works. In Section 6, we compare on some simulations the performance of the partitioning schedulings and we explain the benefits of our approach. We summarize the contributions of this paper in Section 7 and we give direction for our future work.

Terminology

In this paper, we consider an application built from a set τ = {τ 1 , τ 2 , . . . , τ n } of n sporadic real-time tasks. A sporadic task is a recurring task for which only a upper bound on the separation between release times of the jobs is known. Each task τ i is characterized by a minimum interarrival time T i (also denoted period), a worst-case execution time C i and a relative deadline D i . This application runs on a platform Π = {π 1 , π 2 , . . . , π m } of m identical processors (homogeneous case). We consider a fixed-priority scheduling on each processor. A fixed-priority scheduler assigns a priority to each task and all jobs of a task is released with the fixed priority of this task. We assume that tasks are indexed by decreasing priority: ∀i = 1, . . . , n -1, task τ i has a higher priority than task τ i+1 . A partitioning algorithm produces a partition P art(τ ) = {τ 1 , τ 2 , . . . , τ m } of m disjointed subsets where each subset τ j of real-time tasks is executed on processor π j . The subset τ j composed by n j tasks is also denoted by τ j = {τ j 1 , τ j 2 , . . . , τ j nj }. In the rest of this paper, we refer to τ i when the considered task is taken independently and to τ j i when it is considered assigned on processor π j . We define u i as the utilization of task τ i : u i ≡ Ci Ti and U j as the utilization of the taskset τ j : U j ≡ τ j k ∈τ j u k . On the processor π j , we denote lp j (i) (respectively hp j (i)) the subset of real-time tasks assigned to π j which have a priority lower than (respectively higher than or equal to) τ i . The response time of the task τ i is denoted R i . We denote R k i the k th iteration in the response time computation of the task τ i .

Robustness

We consider the robustness in the real-time systems as the capacity of the system to handle WCET overruns faults when the WCET are estimated. If the WCET of all the tasks of the system has been well defined, a feasibility analysis shows wheter the system is feasible. But in pratice, it may possible that a task makes a fault or that the time constraints has been miscalculated. Some real-time specifications -such as Real-Time Specification for Java [START_REF] Dibble | Jsr 1: Real-time specification for java[END_REF] -provide mechanism to handle cost overruns and deadline misses in the case of estimated WCET.

In this work, we consider the robustness as the capacity of a system to meet all the deadlines. We can guarantee that the system stay feasible if and only if we know the execution duration during which a task can exceed its WCET without any deadline is missed. This duration is denoted allowance and the more each task allowance has, the more robust the system regarding to our definition is.

Allowance concept

The allowance of a task is used as a criterion for allocating a task on a processor by our heuristic. We define the allowance A j i of a task τ j i on the processor π j as follows :

Definition 1. Let τ j be a given set of tasks assigned on processor π j . The allowance A j i of a task τ j i of τ j is the maximum duration which can be added to the WCET C i of τ j i such as all tasks of τ j meet their deadlines. We identified in the literature two approaches to compute the allowance: one based on a (Worst Case Response Time) WCRT computation and one based on a sensitivity analysis on the WCETs. The response time of a task is the duration between the time this task has been released and the time it has been finished. The WCRT of a task is the response time of this task in the worst activation scheme. We use the taskset given in Tab.1 to describe the two allowance computation methods in the following subsections.

Allowance Computed from WCRT

One of the available approach to compute the allowance of the execution duration has been proposed by Bougueroua et al. [START_REF] Bougueroua | Dealing with execution-overruns to improve the temporal robustness of real-time systems scheduled FP and EDF[END_REF]. For a given value of allowance A j i , this method consists in checking that the system remains schedulable when the execution duration of a task τ i is equal to

C i = C i + A j i .
In other words, this method consists in checking that the WCRT of all the tasks remains less than or equal to their deadline when their WCET is extended to C i +A j i . The 3 following equations perform this check for the task τ j i on a processor π j if τ j i was assigned on τ j .

U j + A j i T i ≤ 1 (1)
R k+1 i = C i + A j i + τ h ∈hp j (i) R k i T h C h ≤ D i (2)
∀τ l ∈ lp j (i), R k+1 l = C l + τ h ∈hp j (l) R k l T h C h + R k l T i A j i ≤ D l (3)
The value of allowance of a real-time task τ i is found by a binary search. Equation ( 1) tests if the utilization U j of the system when the WCET of τ j i is extended to C i +A j i does not exceed processor utilization. An upper bound A j i,up on the allowance of the task τ i can be found from (1):

A j i,up = (1 -U j ) • T i (4)
Equation ( 4) allows to bound the binary search in [0, A j i,up ]. For the task τ j 1 of our example, A j 1,up = (1 -0.58) • 70 = 29. We can carry out a binary search with 0 ≤ A 1 ≤ 29. Equation (2) tests if the response time R 1 of τ 1 , when its WCET has been extended to C 1 + A 1,j , does not exceed its deadline D 1 . Equation ( 3) tests if the response times R l of all tasks τ j l of lower priority than τ j 1 don't exceed their deadlines D l when the WCET of τ j 1 is extended by A 1,j . For the value 29/2 = 14, we must check that this value satisfies (2) and (3). For A 1,j = 14, the response time of

τ j 1 is R 1 = 24 ≤ D 1 .
We also obtain R 2 = 39 ≤ D 2 , R 3 = 69 ≤ D 3 and R 4 = 177 ≤ D 4 . Then A 1,j = 14 is a valid value of allowance for τ j 1 on the processor π j . We continue the binary search until A 1,j = 21, then A 1 = 21 is the maximum value of allowance for τ 1,j .

The complexity of this approach is pseudo-polynomial due to the Response Time Analysis (RTA) in [START_REF] Leung | On the complexity of fixed-priority scheduling of periodic, real-time tasks[END_REF] and in (3). This complexity is in O(n 2 ) where n is the number of tasks. Indeed, for a task τ i , a RTA is performed in O(n) and for each task of lower priority than τ i , a RTA is performed in O(n). In the worst case, there is n -1 tasks of lower priority than τ 1 , thus the complexity is O(n 2 ).

Allowance Computed from Sensitivity Analysis

Another approach to compute the allowance of the execution duration is the sensitivity analysis. This approach has been introduced by Bini et al. [START_REF] Bini | Sensitivity analysis for fixed-priority real-time systems[END_REF]. This approach is attractive compared to the previous one because no iterative computation (such as WCRT computation) is needed. The authors propose to consider the system only at time corresponding to the activation time of the highest priority tasks in [0, D i ] union time {D i }. The maximum allowance A j i of a task τ j i on the processor π j is computed by the following equations:

Sens(k) = max t∈schedP k t -   C k + h∈hp(k) t T h C h   t/T i (5)
A j i = min

k∈lp(i) Sens(k) (6)
where schedP k is the set of scheduling points defined by schedP k = P i-1 (D k ) and P k (t) is defined by :

P 0 (t) = {t} P k (t) = P k-1 ( t T k T k ) ∪ P k-1 (t) (7)
For the task τ 1 of our example, schedP 1 = D 1 = {60}. In the same way, 5. The allowance M 1 = min k∈lp j (i) Sens(k) = 21 is obtained from ( 5) and [START_REF] Oh | Allocating fixed-priority periodic tasks on multiprocessor systems[END_REF]. The complexity of this approach is exponential because |schedP n |, denoted as the size of the scheduling points set computed by ( 7) is 2 n-1 in the worst-case. But in practice, the size of schedP n << 2 n-1 for a great value of n. We notice that the size of schedP n is highly sensitive to the range of task periods as seen in appendix of [START_REF] Davis | Efficient exact schedulability tests for fixed priority real-time systems[END_REF].

Task Partitioning Problem

The task partitioning problem consists of finding a partition of a taskset τ in m subsets τ j , 1 ≤ j ≤ m such that each subset is feasible on processor π j . Since it has been proved that bin-packing problem (NP-hard in the strong sense) can be reduced in polynomial-time to a task partitioning problem [START_REF] Garey | Computers and Intractability: A Guide to the Theory of NP-Completeness[END_REF], no optimal algorithm exists to decide in polynomial-time if a given taskset is feasible. Fortunately, approximation algorithms and heuristics exist to find solutions for the task partitioning problem in polynomial-time. Heuristics for the tasks partitioning problem exist and are versions of the heuristics proposed for the bin-packing problem. The more cited in the literature are First-Fit [START_REF] Fisher | The partitioned scheduling of sporadic tasks according to static-priorities[END_REF] (FF), Best-Fit [START_REF] Oh | Allocating fixed-priority periodic tasks on multiprocessor systems[END_REF] (BF) and Next-Fit (NF) [START_REF] Andersson | Preemptive multiprocessor scheduling anomalies[END_REF]. These heuristics have been initially designed to minimize the number of bins (respectively the number of processors) for the bin-packing problem (respectively the tasks partitioning problem). Another heuristic Worst-Fit [START_REF] Coffman | Approximation algorithms for bin packing: A survey[END_REF], is rarely used because it provides poor performance to solve the bin-packing problem. On the problem of task partitioning, this heuristic allocates tasks to processors where utilization is the lowest. This approach is relevant because we want the best allocation of tasks to maximize the allowance of tasks. We propose in the next subsection a heuristic which allocates tasks to processors that has the greatest allowance rather than processors that has the lowest utilization.

Allowance-Fit-Decreasing

We propose a heuristic, denoted Allowance-Fit-Decreasing, to solve the task partitioning problem. We want to tolerate bounded WCET overruns, a property not considered in classical heuristics. WCETs overruns can be due to OS approximations, faults of the task or WCET under-estimation. By definition, the allowance of a processor is the minimum allowance for all task allocated to the processor. Our goal is to propose a partitioning scheme that assigns a task to the processor whose allowance is maximum.

We describe the Allowance-Fit-Decreasing heuristic with the pseudo-code given in Alg.1. The tasks are first sorted according to their utilization by function sort task by decreasing utilization() (line 1). For each task τ i of the taskset τ (iteration loop at lines 2-17), the proc parameter, denoting the processor on which τ i is allocated (at line 3), is initialized with a null value. The minimum value of allowance for the entire system (variable A min at line 4) is first initialized to minus infinity. We then consider all processors in Π and find the processor that maximizes the processor allowance. For each processor π j (iteration loop at lines [START_REF] Burchard | New strategies for assigning realtime tasks to multiprocessor systems[END_REF][START_REF] Oh | Allocating fixed-priority periodic tasks on multiprocessor systems[END_REF][START_REF] Andersson | Preemptive multiprocessor scheduling anomalies[END_REF][START_REF] Fisher | The partitioned scheduling of sporadic tasks according to static-priorities[END_REF][START_REF] Lakshmanan | Partitioned fixed-priority preemptive scheduling for multi-core processors[END_REF][START_REF] Dibble | Jsr 1: Real-time specification for java[END_REF][START_REF] Bougueroua | Dealing with execution-overruns to improve the temporal robustness of real-time systems scheduled FP and EDF[END_REF], our heuristic finds the minimum value of processor allowance A j min computed by function compute allowance(π j , τ j i ) when τ i is allocated to π j with the method described in Sect.4.2. If after the iteration loop, A min is greater than or equal to 0 then, τ i can be assigned to a processor, the one that maximizes the processor allowance, by construction. We then proceed with the other tasks until all the tasks until either all tasks have been assigned to a processor (we then return that the task set is schedulable in line 18) or one task is declared not schedulable (line 15).

Partitioned Scheduling Algorithm

A partitioned scheduling algorithm is the combination of a task partitioning algorithm with a schedulability condition. We build two partitioned scheduling algorithms, one from Worst-Fit and one from Allowance-Fit-Decreasing. For the first one, the schedulability condition is a necessary and sufficient condition implicitly given by Allowance-Fit-Decreasing. Indeed, our heuristic computes with the function compute allowance(π j ,τ i ) the value of A j min , the minimum allowance for all the tasks assigned to processor π j , including τ i . If this function returns a negative value, then τ i cannot be assigned on the processor π j . A j min is computed from the sensitivity analysis given in Sect.4.2. For the second one, we combine Worst-Fit with the necessary and sufficient schedulability condition RTA [START_REF] Audsley | Applying new scheduling theory to static priority pre-emptive scheduling[END_REF].

During the allocation, the tasks are taken in order of their decreasing utilization. In other words, the tasks with the greater utilization are allocated first. We consider a fixed-priority assignment and we use DM priority assignment since this PA is an optimal one when the considered tasks have constrained-deadlines (∀i, D i ≤ T i ) [START_REF] Leung | On the complexity of fixed-priority scheduling of periodic, real-time tasks[END_REF].

Methodology

Our simulations is based on randomly generated tasksets. Because we focus on tasks with constrained-deadlines, we consider tasksets such that for any task τ i , α = Di Ti ). We randomly generate 10 sets for α ∈ [0.1, 0.2, . . . , 1.0]. Each set is built from 100,000 randomly generated tasksets. Each taskset is composed by 24 tasks. A taskset is generated using the UUniFast algorithm [START_REF] Bini | Biasing effects in schedulability measures[END_REF] which produces a uniformly distributed set of task utilizations and which avoids bias in the generated tasksets. For a given task τ i , the period T i is generated with a uniform distribution between 100 and 100,000 ms. The deadline D i is given by D i = α • T i and the WCET C i is given by C i = u i • T i . We consider an homogeneous processor composed by 8 identical processors.

Simulation Results

We show in Fig. 1 the average number of iterations during the computation of allowance.

R k+1 i = C i + τ h ∈hp j (i) R k i T h C h (8)
We implement a function iteration() which computes the value given by [START_REF] Fisher | The partitioned scheduling of sporadic tasks according to static-priorities[END_REF].

The computation of the allowance based on the WCRT computation calls this function in [START_REF] Leung | On the complexity of fixed-priority scheduling of periodic, real-time tasks[END_REF] and in [START_REF] Davis | Robust priority assignment for fixed priority real-time systems[END_REF]. The computation of the allowance based on the sensitivity analysis calls this function in [START_REF] Burchard | New strategies for assigning realtime tasks to multiprocessor systems[END_REF]. For the two implementations of the allowance computation, we count the number of calls to the function iteration() for each randomly generated taskset and we keep the average of the allowance over all the taskset. We notice that despite the fact of the complexity of the sensitivity analysis seems greater than the complexity of the allowance computation based on the WCRT, the number of iterations for sensitivity analysis is well below the number of iterations for computation based on the WCRT. We therefore choose to compute allowance by sensitivity analysis [START_REF] Bini | Sensitivity analysis for fixed-priority real-time systems[END_REF]. We show in Fig. 2, 3 and 5 the comparison between First-Fit-Decreasing (FFD), Worst-Fit-Decreasing (WFD) and Allowance-Fit-Decreasing (AFD). Decreasing means that these heuristics assign the tasks in order of their decreasing utilization. We voluntarily omit to show results concerning BF and NF because their behavior are very similar to FF. AFD and WFD use a necessary and sufficient condition of schedulability. Therefore we used the necessary and sufficient condition RTA [START_REF] Audsley | Applying new scheduling theory to static priority pre-emptive scheduling[END_REF] for the heuristics FFD.

We show in Fig. 2 the number of partition found by the heuristics FFD, WFD and AFD. We notice that FFD provides a slightly better schedulability for α ≥ 0.4. This result is explained by the fact that FFD is one of the best heuristic for the task partitioning problem in terms schedulability. But the gap between FFD and the other two heuristics is not very large and may be acceptable if we want more robustness to WCETs overruns. We compare in Fig. 3 the minimum allowance obtained by the three heuristics on 4 processors and in Fig. 4 on 8 processors. Minimum allowance A min guarantees that any task of the system can bear an interference during A min without any deadline is missed. We show that AFD and WFD outperforms largely FFD. Indeed, AFD and WFD distributes the tasks among the processors instead of fills up all the first processors. We note that AFD is slightly better than WFD. We show in Fig. 5 the comparison between the computational time of the three heuristics. AFD offers better results than WFD in terms of minimum allowance. But the computation time of AFD is 6 to 10 greater than the computation time of WFD. For a robust allocation to the WCETs overruns, it is interesting to use AFD. But when tasks must be accepted online, WFD is a preferable choice.

Conclusion

We have proposed a fixed-priority partitioned scheduling for homogeneous processors which maximizes the allowance of the execution duration. This scheduling is more robust than the others based on FF, BF or NF because during allocation of the tasks, the processor offering the greatest value of allowance is chosen. In terms of maximization of allowance, Allowance-Fit is slighty efficient than Worst-Fit. But in terms of computation time, Worst-Fit is largelly better than Allowance-Fit. Thus Worst-Fit is a good heuristic to maximize the robustness of a partitioned system of real-time tasks. In a future work, we will extend this approach to the class of restricted migration scheduling to improve the schedulability of our solution. In such a scheduling, the different jobs of a recurring task can migrate from a processor to another, but no migration is allowed during the execution of the job.

schedP 2 =

 2 {70, 85}, schedP 3 = {70, 100, 140, 190} and schedP 4 = {140, 200, 210, 260}. The values Sens(1) = 50, Sens(2) = 45, Sens(3) = 33.33 and Sens(4) = 21.66 are computed by using Equation

Algorithm 1 : 1 foreach τi in τ do 2 proc = None; 3 Amin = -Inf; 4 foreach πj in Π do 5 A 6 if A j min > Amin then 7 proc = πj; 8 Amin = A j min ; 9 end 10 end 11 if Amin ≥ 0 then 12 assign τi to proc 13

 112345678910111213 Allowance-Fit-Decreasing sort task by decreasing utilization() j min = compute allowance(πj, τ j i )

Fig. 1 .Fig. 2 .

 12 Fig. 1. Comparison of computation time for the allowance computation approaches.

Fig. 3 .

 3 Fig. 3. Comparison of minimum allowance for the tasks partitioning heuristics on 4 processors.

Fig. 4 .Fig. 5 .

 45 Fig. 4. Comparison of minimum allowance for the tasks partitioning heuristics on 8 processors.

Table 1 .

 1 System of 4 sporadic real-time tasks.

	Ci Di Ti
	τ1 10 60 70
	τ2 15 85 100
	τ3 30 190 210
	τ4 45 260 320


Guillaume Blin
email: gblin@univ-mlv.fr

Romeo Rizzi
email: romeo.rizzi@uniud.it

Florian Sikora
email: sikora@univ-mlv.fr

Stéphane Vialette
email: vialette@univ-mlv.fr

Minimum Mosaic Inference of a Set of Recombinants

Keywords: Minimum Mosaic problem, SNP, Haplotype

investigate the central problem of finding recombination events

Introduction

Given any two unrelated people, their DNA sequences will only differ on about 0.1%. This small genetic variability is of particular importance since it influences how people differ in their risk of disease or their response to drugs. A main challenge is to discover the DNA variants that contribute to common disease risk. These variations mostly arise on specific single positions called Single Nucleotide Polymorphisms (SNPs) where the corresponding nucleotides, called alleles, differ. For example, given three DNA fragments . . . GGACCTG . . ., . . . GGACATG . . . and . . . GGACTTG . . . from different individuals of a population, the SNP is composed of the three alleles C, A and T. Furthermore, each person has two almost identical copies of all chromosomes except the sex chromosomes. Given a population, an haplotype refers to a combination of alleles at different positions on a chromosome. Therefore, each individual has two haplotypes (one maternal and one paternal) for each chromosome. The genotype corresponds to the set of alleles that a person has thereby defining at each SNP the alleles of the two haplotypes. Unfortunately, the current methods of data acquirement do not provide inheritance information.

As a simple illustration, consider the following two partial haplotypes (a chromosome region where only Alleles of an SNP are called heterozygous if they differ (e.g. A/C), and homozygous otherwise (e.g. A/A). Since most SNPs are composed of only two alleles (among the 16 possibilities) -that occur in a large percentage of the population -haplotypes are usually represented by binary sequences (one character for each of the two possible alleles per SNP) whereas genotypes are usually represented by ternary sequences (0 and 1 denote the two homozygote alleles and 2 denotes the heterozygote one). Notice that it is common to denote missing values by the extra symbol "-". Therefore, the sequences are built on the alphabet {0, 1, 2, -}.

Genetic variation within species is mostly induced by a process called recombination. Given two equal length sequences, a recombination generates a third sequence of the same length by concatening the prefix of one sequence with the suffix of the other sequence [START_REF] Koivisto | Recombination systems[END_REF]). In the resulting sequence, the assembly point is referred as a breakpoint. An illustration is given in Figure 1. Finding recombination events has become a central problem in computational biology [START_REF] Kececioglu | Reconstructing a history of recombinations from a set of sequences[END_REF][START_REF] Ukkonen | Finding founder sequences from a set of recombinants[END_REF][START_REF] Schwartz | Inferring piecewise ancestral history from haploid sequences[END_REF][START_REF] Koivisto | Recombination systems[END_REF][START_REF] Rastas | Haplotype inference via hierarchical genotype parsing[END_REF][START_REF] Wu | Improved algorithms for inferring the minimum mosaic of a set of recombinants[END_REF]. In most combinatorial models, a present population is assumed to be a descendent of a small number of specific sequences called founders. Due to recombination, a present sequence, called a recombinant, is thus composed of blocks from the founders. The term mosaic is often used to denote the mosaic-like structure of DNA induced by the recombinations (Ukkonen 2002) (c.f. Figure 2).

A major question related to founder sequences is the so-called Minimum Mosaic problem: using the natural parsimony criterion for the number of recombinations, find the "best" founders. More formally, the Minimum Mosaic problem, in its natural decisional form, is defined as follows: Given a set D of m This paper is organized as follows. In Section 2 we present the related state-of-the-art, and we introduce in Section 3 the needed material. Section 4 is devoted to computational complexity. In Section 5, we propose some exact exponential-time algorithms for the Minimum Mosaic problem.

Related works and known results

Ukkonen first formulated an optimization form of the problem based on the mosaic model and parsimony [START_REF] Ukkonen | Finding founder sequences from a set of recombinants[END_REF]. He considered two criteria: the number of founders and the number of recombinations in a solution. Ukkonen gave two polynomial-time algorithms: an O(n(m + K 3 )) time algorithm for the Minimum Segmentation problem and an O(mn) time algorithm for the Minimum Mosaic problem for K = 2 (Wu et al. gave a similar result [START_REF] Wu | Improved algorithms for inferring the minimum mosaic of a set of recombinants[END_REF]). Also, Ukkonen designed an O(nK O(m) ) time dynamic programming algorithm for the general (i.e. K ≥ 2) Minimum Mosaic problem (the time complexity was guessed from the description given in the paper and will be justified in Section 5). According to the authors, this latter algorithm does not perform well for a moderate number of founders and/or recombinants.

A different (no restriction on the number of founders) but related problem was introduced by Gusfield [START_REF] Gusfield | Haplotyping as perfect phylogeny: conceptual framework and efficient solutions[END_REF]. In the Haplotype Inference problem, we are given a set of n genotype sequences and the goal is to find a set of n pairs of haplotype sequences (one pair per genotype) that is a good "explanation" of the given genotype sequences. For biological pertinence, a solution must be compatible with (or guided by) a given perfect phylogeny. Rastas et al. [START_REF] Rastas | Haplotype inference via hierarchical genotype parsing[END_REF] observed that this latter problem is equivalent for K > √ 2m to the Minimum Mosaic for genotype recombinants without missing values (i.e., D ⊂ {0, 1, 2} n ), that was proved to be NP-complete in [START_REF] Lancia | Haplotyping populations by pure parsimony: Complexity of exact and approximation algorithms[END_REF].

In [START_REF] Rastas | Haplotype inference via hierarchical genotype parsing[END_REF]), Rastas et al. considered the Minimum Mosaic problem with missing values and mutations (i.e. mismatches between founders and recombinants) leading to a different parsimony criterion. More precisely, for each recombinant r ∈ D, a score k + k c is computed where (1) k is the number of breakpoints of a recombinant r such that the Hamming distance between r and r is k (i.e. the number of mutations) and (2) c is the relative weight for a mutation compared to recombinations. Rastas et al. improved the complexity of Minimum Segmentation problem by providing an O(nmK) algorithm.

They also proved that the Minimum Mosaic problem for haplotype recombinants with possible missing values (i.e. D ⊂ {0, 1, -} n ) is NP-complete. Finally, they proved that even when missing values are forbidden, the Minimum Mosaic problem for haplotype recombinants (i.e. D ⊂ {0, 1} n ) is NP-complete even for K = 2 since it becomes equivalent to the NP-complete Hypercube Segmentation problem [START_REF] Kleinberg | Segmentation problems[END_REF]. One has to notice that this latter result is a bit "artificial", in the sense that the extra condition that the mutation cost c = 1 nm (necessary in the proof) roughly corresponds to forbid (and thus ignore) recombinations. Indeed, nm mutations becomes less expensive than a single recombination event. This is precisely the reason why Ukkonen [START_REF] Ukkonen | Finding founder sequences from a set of recombinants[END_REF] and Wu et al. [START_REF] Wu | Improved algorithms for inferring the minimum mosaic of a set of recombinants[END_REF]) have been able to find a polynomial-time algorithm for the problem without mutation. One of the purpose of our contribution is to provide a better complexity insight of the problem by proving its complexity without relying on a somewhat unrealistic prohibitive mutation cost function.

Several heuristics and lower bounds have been proposed [START_REF] Roli | Tabu search for the founder sequence reconstruction problem: A preliminary study[END_REF][START_REF] Wu | Bounds on the minimum mosaic of population sequences under recombination[END_REF]). Wu et al. [START_REF] Wu | Improved algorithms for inferring the minimum mosaic of a set of recombinants[END_REF]) designed an O(nm) time algorithm for the Minimum Mosaic problem for genotype recombinants without missing values (i.e. D ⊂ {0, 1, 2} n ) and K = 2. Furthermore, they gave an exact algorithm for the general case (haplotype recombinants without missing values). Notice that the time complexity of this latter algorithm is not given in the paper and the authors claim it to be practical for moderate n and m.

Finally, in [START_REF] Zhang | Genotype sequence segmentation: Handling constraints and noise[END_REF]), Zhang et al. investigated the Minimum Segmentation problem for genotype recombinants and provided two dynamic programming algorithms with time complexity O(nK 4 ) and O(nK + P K 4 ), where P is the number of rows of the dynamic table.

Notations

In the rest of the paper, we do not consider missing values (i.e., haplotype are defined on {0, 1} n ). Given any string s = s 1 s 2 . . . s n , and two integers i and j, 1 ≤ i ≤ j ≤ n, we denote by s[i, j] the substring s i s i+1 . . . s j of s. Given a set of haplotype founders f 1 , f 2 , . . . , f K , each in {0, 1} n , and an haplotype recombinant r of size n, a segmentation of r out from

f 1 , f 2 , . . . , f K is a partition of the in- terval [1 . . . n] = {1, 2, . . . , n} into consecutive inter- vals I 1 , I 2 , . . . , I k such that, for each 1 ≤ i ≤ k, we have r[I i ] = f j [I i ] for some j ∈ {1, 2, . . . , K}. The cost of the segmentation is k -1. For exam- ple, for r = 001011, f 1 = 000000 and f 2 = 111111, {I 1 = f 1 [1 . . . 2], I 2 = f 2 [3 . . . 3], I 3 = f 1 [4 . . . 4], I 4 = f 2 [5 . . . 6]} is a segmentation of r out from f 1 and f 2 of cost 3. The cost of segmenting r out from f 1 , f 2 , . . . , f K is denoted cost(r; f 1 , f 2 , . . . , f K ).
(the cost of a best segmentation can be found in O(nK) time [START_REF] Rastas | Haplotype inference via hierarchical genotype parsing[END_REF]).

The Minimum Mosaic problem can be defined as follows: Given a set of m recombinants D = {r 1 , r 2 , . . . , r m }, r i ∈ {0, 1} n for 1 ≤ i ≤ m, and an integer K, find a set of

K founders F = {f 1 , f 2 , . . . , f K }, f i ∈ {0, 1} n for 1 ≤ i ≤ K, such that m i=1 cost(r i ; f 1 , f 2 , . . . , f K ) is minimized.

Hardness result for the Minimum Mosaic problem

In this section, we prove the Minimum Mosaic problem for an unbounded K to be NP-complete without relying on unrealistic prohibitive mutation costs (as done in [START_REF] Rastas | Haplotype inference via hierarchical genotype parsing[END_REF]). For the sake of presentation, we first generalize the problem to arbitrary strings. The following lemma proves that this can be done safely.

Lemma 1 If the Minimum Mosaic problem on arbitrary strings is NP-hard, then so is the Minimum Mosaic problem on binary strings.

Proof: Assume to be given a natural K and a set of m recombinants D = {r 1 , r 2 . . . r m } ⊂ Σ n where Σ = {σ 1 , σ 2 . . . σ k } is any alphabet on k symbols. Then, take any encoding of the symbols in Σ by binary strings of length log 2 k . In other words, let δ : Σ → {0, 1} k be any injection from Σ to {0, 1} k with k = log 2 k . Once such an encoding has been fixed, we can extend δ to get an injective function which maps every string

r i over Σ into a bi- nary string δ(r i ) of length |δ(r i )| = k |r i |.
Notice that any feasible solution f 1 , f 2 . . . f K for the instance K, r 1 , r 2 . . . r m naturally maps into the feasible solution δ(f 1 ), δ(f 2 ), . . . , δ(f K ) for the instance K, δ(r 1 ), δ(r 2 ) . . . δ(r m ) , and the cost remains unaffected.

Let us show that the converse is also true. To do so, we claim that, given any feasible solution f 1 , f 2 . . . f K for the instance K, δ(r 1 ), δ(r 2 ) . . . δ(r n ) , we can always modify it, without increasing its cost, in such a way that each of the f i is actually the binary encoding of some string over Σ.

In other words, we can assume that

f 1 = δ(f 1 ), f 2 = δ(f 2 ) . . . f K = δ(f K ).
In this way, the converse would directly follow. In order to prove our claim, we can actually act on the f i 's and "clean" them out one by one. Assume f i is not cleaned, that is, f i does not belong to the image of map δ. Then there exists some 1 ≤ j ≤ n such that f i [k (j -1) + 1, k j] does not belong to the codewords set ∆(Σ) := {δ(σ) : σ ∈ Σ }. In this case, where σ is any symbol in Σ such that δ(σ) has a longest suffix in common with f i [k (j -1) + 1, k j], we modify f i precisely on the interval [k (j -1) + 1, k j], and, more precisely, by replacing it by the string δ(σ). Considering the way the procedure for finding a minimum cost segmentation operates (e.g. Algorithm 1), it is easy to see that this modification does not increase the cost.

Algorithm 1 Algorithm for finding a minimum cost segmentation : F S(r; f 1 , f 2 , ..., f K )

1: Let i = |r| 2: if i = 0 then 3: return 4: end if 5: Let j be the smallest positive natural such that r[j..i] = f k [j..i] for some k 6: if j > i then 7: return "No production exists" 8: end if 9: return F S(r[1, j -1]; f 1 [1, j -1], f 2 [1, j - 1], . . . , f K [1, j -1]) + [j..i]
In order to facilitate the understanding of the proof, let us further strengthen the formulation of our problem: consider the more general formulation where K of the K founders comprising the solution (and to be given in output) are actually given as specified in advance in the input. We will refer to those specific founders as forced founders and will first show how one can force a part of the founders.

Let us first remark that it is easy to provide a polynomial-time reduction from any instance containning forced founders to an instance without ones. Indeed, for each forced founder f f , add n m copies of f f in the recombinants set D. It is clear that in a solution, one has to include f f in the founder set, otherwise, it will induce at least n m breakpoints which is the maximal number of breakpoints one can get in the original instance. Note that given an objective of at most B breakpoints, it suffices to add B copies of f f to get this property.

Thus, for the purpose of the reduction, we can without loss of generality assume to have instances of the form K, D, F f such that

D = {r 1 , r 2 . . . r m } ⊂ Σ n , F f = {f f 1 , f f 2 . . . f f m } ⊂ Σ n and m = m + m . Let us now provide a reduction from the NP-complete Vertex-Cover problem: Given a graph G = (V, E) such that |V | = n G , |E| = m G
and an integer k G , decide whether there exists a subset V of V such that each edge of G is incident to at least one vertex of V and |V | ≤ k G . Our reduction begins by giving an arbitrary orientation to each edge of G, that is, for each e j ∈ E, let t j and h j be indices such that v tj is the tail and v hj is the head of the arbitrary oriented edge e j .

We consider the alphabet Σ := {W, Z} ∪ {X i : i = 1, 2, . . . n G }. Informally, there is one letter for each vertex of V , while W and Z act like separators letters. We define m = C n G + 3 recombinants, each of length n = 6 m G built as follows (C is a constant defined later):

• r 1 = (W W ZZW W ) m G , • r 2 = m G j=1 (ZZX tj X tj ZZ) and r 3 = m G j=1 (ZZX hj X hj ZZ), • r j i = (X i X i X i X i X i X i ) m G , for each 1 ≤ i ≤ n G , 1 ≤ j ≤ C.
We then define a set F f of K = 1+2m G n G forced founders as follows.

• F f 1 = (ZZZZZZ) m G , • F f i,t = Z 3 t X i X i X i Z 3 (2m G -t-1) , for each 0 ≤ t ≤ 2m G -1, and 1 ≤ i ≤ n G . Finally, we set K = 2m G n G +k G +2; that is asking for K = k G + 1 founders.
Lemma 2 If the graph G = (V, E) admits a vertex cover of size at most k G then the corresponding built instance of our problem admits a solution of cost at most

C(n G -k G )(2m G -1) + 6m G . Proof: Let V ⊂ V be a vertex cover of size k G of G.
For each 1 ≤ j ≤ m G , we let c j be an index such that v cj ∈ V ∩ e j (such an index exists since V is a vertex cover). By definition, c j ∈ {t j , h j }. Let c j = {t j , h j }\c j . We construct our solution with the following set F of K = k G + 1 founders:

• F 0 = m G j=1 (W W X cj X cj W W ), • F i = (X i X i X i X i X i X i ) m G , for each v i ∈ V .
Let us now compute the corresponding cost:

• cost(r 1 ; F ∪ F f ) = 2 m G since we have to switch from F 0 to F f 1 and from F f 1 to F 0 for each m G blocks of length 6, • cost(r 2 ; F ∪ F f ) = 2 m G and cost(r 3 ; F ∪ F f ) = 2 m G .
For each 1 ≤ j ≤ m G , since V is a vertex cover, F cj and F 0 will prevent a breakpoint between the X tj 's (and between the X hj 's) since c j ∈ {t j , h j } and c j = {t j , h j }\c j . Therefore, by switching from F f 1 to F cj or F 0 and back to F f 1 for each m G blocks of length 6, it only induces 2 breakpoints for each block,

• for 1 ≤ j ≤ C, cost(r j i ; F ∪ F f ) = 0 for each v i ∈ V whereas cost(r j i ; F ∪ F f ) = C(2 m G -1) for each v i / ∈ V . Considering each recombinants r j i such that v i / ∈ V , by switching from F f i ,2t-2
to F f i ,2t-1 for the t th block of length 6, it will only induce 2 breakpoints for each block (except the last one). Since

|V | = k G , there are n G -k G recombinants which costs (2m G -1) each.
On the whole,

2 i=1 cost(r i ; F ∪ F f ) + n G i=1 C j=1 cost(r j i ; F ∪ F f ) = 2m G + (2 × 2m G ) + C(n G -k G )(2m G -1) = C(n G -k G )(2 m G -1) +6 m G .
We now turn to considering the reverse direction. Lemma 3 Given a graph G = (V, E) and the corresponding built instance of our problem, if the latter one admits a solution of cost at most C(n Gk G )(2m G -1) + 6m G , for any C > 6m, then G admits a vertex cover of size at most k G . Proof: Let us prove that the cost of recombinants

r j i , 1 ≤ i ≤ n G , 1 ≤ j ≤ C, for any solution is greater than C((n G -k G )(2m G -1)). Note first that, only considering the set of forced founders F f , each r j i , 1 ≤ i ≤ n G , 1 ≤ j ≤ C, has a cost(r j i , F f ) = C(n G (2m G -1)). Indeed, given a recombinant r j i , one has to switch from F f i,t to F f i,t+1 for 0 ≤ t ≤ 2m G -2. Considering now both F f and F , if the set of K = k G +1 founders is composed of exact copies of K dif- ferent recombinants {r j1 i1 , r j2 i2 . . . r j k G +1 i k G +1 } (for some 1 ≤ i 1 , i 2 . . . i k G +1 ≤ n G and 1 ≤ j 1 , j 2 . . . j k G +1 ≤ C), then n G i=1 C j=1 cost(r j i ; F ∪ F f ) = C(n G -(k G + 1))(2m G -1).
However, r 1 is built over some W letters; which does not belong to any forced founders. Therefore, there should be founders among the K including a W in positions 6t + 1, 6t + 2, 6t + 5, 6t + 6, for each 0 ≤ t ≤ m G -1. Then, for any 0 ≤ t ≤ m G -1, there will exist among the recombinants {r j1 i1 , r j2 i2 . . . r

j k+1
i k+1 } one recombinant -say r j1 i1 -that will induce a switch (for each of its C copies) from a founder of F to F f i1,t+1 (resp. F f i1,t+2 ) due to the W W at positions 6t+1, 6t+2 (resp. 6t+5, 6t+6) in the corresponding founder. Therefore, on the whole, one will end-up with an extra cost of 2Cm G . In the best case, one should only "sacrifice" one of the K founders (which was only allowing a gain of 2C(m G -1) breakpoints). On the whole,

n G i=1 C j=1 cost(r j i ; F ∪ F f ) ≥ C(n G -k G )(2m G -1). Moreover, considering the objective cost (i.e. C(n G -k G )(2m G -1) + 6m G )
and that C > 6m, the following is enforced:

1. Considering any position, each of the K founders has a different letter, 2. A founder with a letter X u in a given position contains only letters X u in the next positions until a letter W or the end of the founder is encountered, 3. All the W letters occur in the same founder.

Let us now prove each of these points.

Point 1. By contradiction, suppose there are two founders with the same letter in a given position. Then, there is at least one more breakpoint for each copy of a given recombinant of {r j1 i1 , r j2 i2 . . . r

j k G +1
i k G +1 } which leads to a cost above the objective.

Point 2. By contradiction, in an optimal solution, suppose there is a founder with a letter X u at position k followed by a letter X v . Then, due to Point 1, there should not exist another founder such that X u is at position k nor X v at position k + 1. Therefore, there is a breakpoint between positions k and k + 1 in all r j u (resp. r j v ); that is 2C breakpoints. Consider now changing X v into X u in the corresponding founder. Then it will induce at least C less breakpoints, leading to a better solution; a contradiction.

Point 3. Let us first prove that W letters at position 6t + 1 and 6t + 2 belong to the same founder. By contradiction, suppose it is not the case; e.g. there is a W letter at position 6t + 1 in founder F i and a W letter at position 6t + 2 in founder F j , i = j. By Point 1, there is no W letter at position 6t + 2 (resp. 6t + 1) in founder F i (resp. F j ). Thus, there is a breakpoint in the recombinant r 1 between positions 6t + 1 and 6t + 2. Consider now swapping F i [1..6t + 1] and F j [1..6t + 1]. This will decrease the number of breakpoints at least by one since it only changes the cost due to positions 6t+1 and 6t+2. Same argument holds for the other W positions.

We now know that the W letters appears by block of size 4 (except for the first and last blocks; which are of size 2). It also appears that if the blocks of letters W are not on the same founder, it creates at least C more breakpoints. Figure 3 illustrates this configuration. Indeed, in an optimal solution, consider two founders F i and F j such that there are W letters at positions 6t+1 and 6t+2 (resp. 6t+5 and 6t+6) in F i (resp. F j ). Assume, moreover, that a X u (resp. X v ) appears in position 6t + 3 (resp. 6t + 4) in F i (resp. F j ). Then considering the forced founders, F i and F j , there will be one breakpoint between positions 6t + 3 and 6t + 4 in all r x u and r x v , for 1 ≤ x ≤ C and one breakpoint between positions 6t+6 and 6t+7 in all r x v , for 1 ≤ x ≤ C. Consider now swapping F i [1..6t + 4] and F j [1..6t + 4] and replacing the X i 's of F j by X j 's. It will induce at least C less breakpoints, leading to a better solution; a contradiction.

In the following, let us assume that F 0 denote the founder containning the W letters; the other founders are referred as F 1 , F 2 , ..., F K . Now that we know the precise topology of the founders of any optimal solution, we can now prove that {F 1 , F 2 . . . F K } corresponds to a vertex cover in G. Let V ⊆ V be defined as follows. For each founder F i , 1 ≤ i ≤ K , if F i is the exact copy of the recombinant r y

x , then add v x in V . Let us prove that V is a vertex cover of G. By contradiction, suppose it is not. Then, it exists at least one y such that (v hy , v ty ) ∈ E and both {v hy , v ty } ∩ V = ∅. Therefore, there is at least a breakpoint between X ty and X ty in r 2 or between X hy and X hy in r 3 since neither h y nor t y appears in {F 1 , F 2 . . . F K }, and F 0 can contain at most one of the two in the corresponding positions. It is worth noting that none of the forced founders cannot prevent theses breakpoints since there never have the same letters at position 6t + 3 and 6t + 4.

Combining Lemma 2 and Lemma 3 we obtain the following result. Proposition 4 Minimum Mosaic given haplotype recombinants with no missing values is NP-Complete.

ments among the K founders. This ensures that when considering I x y , the choice of the founder for I x y will not interfere with placement of I x y . Indeed, either in the solution (1) they were on the same founder and thus will have a common prefix or (2) were not in the same founder. This leads to an overall O(BKL) time-complexity algorithm where L is the length of the longest segment.

Let us now consider the case where only the number of breakpoints for each recombinant is known (i.e. without the positions of those last). Then, one can test all the possible positions for each recombinant and apply the previously proposed algorithm. On the whole, since given a recombinant r i with b i breakpoints, there are at most n bi possibilities of placement of the breakpoints, the corresponding algorithm runs in O(n b1 .n b2 . . . . .n bm .BKL) = O(n .BKL).

Finally, we give an algorithm with a lower complexity in the case one only knows the number B of allowed breakpoints. Let first give the maximum number of different strings in D, that is the number m d of different recombinants among the m recombinants of D. In the worst case, there is only one breakpoint on each recombinant, that is B different recombinants. Moreover, there is a maximum number of K different recombinants with no breakpoint in D (they are copies of the founders). Thus, we can give the following upperbound for the number of different recombinants in D : m d ≤ K + B. We also have K ≤ m d ≤ m.

Let now decide which recombinants in D will have some breakpoints. Therefore, we have to choose B recombinants among the m d recombinants with at least one breakpoint. We will try all different possibilities. There are m d B = K+B B = (K +B) B different configurations. For each one, run the exact O(nK 2m ) algorithm of Ukkonen to find the best possible founders. This algorithm is ran among a set of only B chosen recombinants. Thus, the running time of the exact algorithm is O(nK 2B ). On the whole, the total running time of the algorithm when only the number of breakpoints is known is O((K + B) B × nK 2B ). This last result demonstrates that once the number of breakpoints is bounded the problem becomes polynomial.

Open problems

In this article, we showed that the Minimum Mosaic problem given haplotypes with no missing values is hard for an unbounded number K of founders. When K is bounded but the number m of recombinants is not, the problem is still widely open. Indeed, there is an ocean between the linear complexity of the problem when m (and thus K) is bounded, the polynomial-time complexity when K = 2, and the NP-Completeness when K is unbounded. It is also widely open whether the problem admits some PTAS.

the

  SNPs are shown) . . . A . . . C . . . C . . . T . . . G . . . T . . . . . . A . . . C . . . A . . . G . . . C . . . T . . . The corresponding genotype is . . . A/A . . . C/C . . . A/C . . . G/T . . . C/G . . . T/T . . .

Figure 1 :

 1 Figure 1: Recombination of S 1 and S 2 leading to T .

Figure 2 :

 2 Figure 2: A Mosaic issued from a set of founders

W not on the same founder a) ..

Exact algorithms for Minimum Mosaic

In this section, we will provide exact algorithms considering a variant of the Minimum Mosaic/Segmentation problems where extra informations (e.g. position, number, . . .) on the breakpoints are given.

Let us first give a complexity study of the general solution (Section 4) proposed by Ukkonen in [START_REF] Ukkonen | Finding founder sequences from a set of recombinants[END_REF]. Indeed, the corresponding solution shows that, once the number of founders (i.e. K) is bounded, then the minimum mosaic problem is FPT. The main idea of the dynamic programming solution of Ukkonen is to compute all the partition of size K of the m input recombinants for each column. The central recurrence is based on the fact that the "best" K-partitions for the i first columns can be computed using one of all the "best" K-partitions of the i -1 first columns and any K-partition of the i th column; the parsimony criterion being the number of breakpoints induced. One has to store in a dynamic table, the minimal cost of any K-partition of the i first columns for all the K-partitions of the i th column. On the whole, since (1) there are K m Kpartitions of any column i and (2) one has to compute the minimum number of breakpoints considering all the K-partitions of column i -1 for each K-partition of column i, the time complexity is O(nK 2m which leads to a polynomial solution when the number of recombinants (i.e. m) is bounded. It also provides a practical solution when the number of founders (i.e. K) is bounded.

Let us now design exact algorithms for a variant of the Minimum Mosaic/Segmentation problems that considers that extra informations on the breakpoints are given. Let us first consider a kind of reverse problem of Minimum Segmentation problem where given a set of m recombinants D = {r 1 , r 2 , . . . , r m } ⊂ {0, 1} n and a set of B identified breakpoints on Di.e. an overall cost B segmentation of the recombinants S = {I 1 1 , I 1 2 . . . I 1 k1 , . . . I m 1 , I m 2 . . . I m km } such that ( m i=1 k i ) -m = B and I x y is the y th segment of recombinant r x -find a set of K founders F = {f 1 , f 2 , . . . , f K } ⊂ {0, 1} n such that the B cost segmentation can be derived from F . We propose a polynomial-time algorithm (Algorithm 2) that solves this problem. For any I x y , let L x y (resp. R x y ) denote the leftmost (resp. rightmost) position of I x y in [1.

.n].

Algorithm 2 Find K founders according to a segmentation S 1: Let L be a sorted list -according to the leftmost position -of elements of S 2: F = {F 1 , F 2 . . . F K } such that F i = [ "", "", . . . , ""], 1 ≤ i ≤ K 3: while L is not empty do for each F i in F do 7:

Z ← the leftmost empty position in F i greater than L x y 8:

.Z] and I x y share a common prefix then 9:

if LCP < Z -L x y then 10: Let us now prove that this algorithm indeed find a solution if one exists. Roughly, the algorithm tries to reassemble the segments in order to produce at most K founders. To do so, the algorithm computes the founders from left to right using the available segments in L. For each available segment I x y , the algorithm tries first to detect if there is a founder that has the longest common prefix with I x y at the given position L x y . If there exist one then the corresponding founder is "merged" with I x y for the positions from L x y to R x y . Otherwise, the algorithm tries to find a founder such that the position L x y is empty and then merge the founder and I x y . Let us prove that, if there exist a solution then the algorithm will find it. First, notice that if there exist such a solution then any segments can be placed entirely in one of the K founders. Thus when considering all the segments starting at a given position, one should be able to find a distribution of the seg-