Maxime Crochemore

Gad M Landau
email: landau@poly.edu

Baruch Schieber

Michal Ziv-Ukelson
email: michalz@cs.technion.ac.il

Re-Use Dynamic Programming for Sequence Alignment: An Algorithmic Toolkit *

The problem of comparing two sequences S and T to determine their similarity is one of the fundamental problems in pattern matching. In this manuscript we will be primarily concerned with sequences as our objects and with various string comparison metrics. Our goal is to survey a methodology for utilizing repetitions in sequences in order to speed up the comparison process. Within this framework we consider various methods of parsing the sequences in order to frame their repetitions, and present a toolkit of various solutions whose time complexity depends both on the chosen parsing method as well as on the string-comparison metric used for the alignment.

Preliminaries

Sequence alignment is traditionally based on the transformation of one sequence into the other via operations of substitutions, insertions, and deletions (indels). The costs of the operations are specified in a given scoring matrix δ. An optimal alignment is an alignment that yields the best similarity score -a value computed as the sum of the costs of the operations applied in the transformation. The alignment of two sequences A and B can classically be solved in O(n 2 ) time [START_REF] Levenshtein | Binary Codes Capable of Correcting Deletions, Insertions and Reversals[END_REF][START_REF] Needleman | A general method applicable to the search for similarities in the amino acid sequence of of two proteins[END_REF][START_REF] Smith | Identification of common molecular subsequences[END_REF] and O(n) space [START_REF] Hirschberg | A Linear Space Algorithm for Computing Maximal Common Subsequences[END_REF] by dynamic programming. The dynamic programming solution to the string comparison computation problem can be represented in terms of a weighted alignment graph [START_REF] Gusfield | Algorithms on Strings, Trees, and Sequences[END_REF]. (See Figure 3. Note that all figures can be found in the Appendix.) An alignment graph for A and B is a directed, acyclic, weighted graph containing (|A| + 1)(|B| + 1) nodes, each labelled with a distinct pair (x, w) (0 ≤ x ≤ |A|, 0 ≤ w ≤ |B|). The nodes are organized in a 1 matrix of (|A| + 1) rows and (|B| + 1) columns. The alignment graph contains a directed edge with a weight of δ(-, b w+1 ) from each node (x, w) to (x, w + 1), and a weight of δ(a x+1 , -) from each node (x, w) to (x + 1, w). Node (x, w) has also a diagonal edge with a weight of δ(a x+1 , b w+1 ) to node (x + 1, w + 1). The weight of a given edge can be specified directly on the grid graph, or as is frequently the case in biological applications, is given by a scoring matrix which specifies the substitution cost for each pair of characters and the deletion cost for each character from the alphabet. Best-scoring paths in the alignment graph represent optimal alignments of A and B. The dynamic programming algorithm will set the value of vertex (x, w) in the graph to the total weight of the highest scoring path which originates in vertex [0, 0] and ends in vertex [x, w], a weight which is identical to the alignment score between the first x characters of A and the first w characters of B.

The two widely used classes of scoring schemes are distance scoring, in which the objective is to minimize the total alignment score, and similarity scoring, in which the objective is to maximize the total alignment score. Within these classes, scoring schemes are further characterized by the treatment of gap costs. A gap is the result of the deletion of one or more consecutive characters in one of the sequences. Additive gap costs assign a constant weight to each of the consecutive characters. For other gap functions which have been found useful for biological sequences, see [START_REF] Gusfield | Algorithms on Strings, Trees, and Sequences[END_REF].

Note that in sections 3.1 and 3.2 we will formalize the alignment problems at hand in terms of distance minimization. (For the sake of simplicity, we will restrict our examples to the Edit Distance measure [START_REF] Levenshtein | Binary Codes Capable of Correcting Deletions, Insertions and Reversals[END_REF].) In this respect the first and second subsections differ in style from subsection 3.3 and section 4, where the alignment problems will be formalized in terms of similarity maximization. However, the distance and similarity perspectives are complementary, and any distance problem can be translated into a similarity problem. Correspondingly, the distance-based solutions in this manuscript can easily be translated to score maximization solutions, in order to apply to string comparison metrics which measure similarity, rather than distance. All solutions described in this manuscript assume a scoring scheme with additive gap costs.

A High-Level Overview of the Toolkit

Our goal is to develop a methodology for utilizing repetitions in sequences in order to speed up the comparison process. Within this framework we consider various methods of parsing the sequences in order to frame their repetitions (see Figures 5 to 7).

Stage One: Utilizing Repetitions in the Source Strings.

In Section 3 of this manuscript we present solutions to the Common Substring Alignment Problem, (CSA for short), which is defined as follows. We are given a set of one or more strings S 1 , S 2 . . . S c and a target string T and Y is a common substring of all source strings S i , that is S i = B i Y F i , and Y may be repeated several times in any of the source strings (see Figures 1,3 and 5). The goal is to compute the similarity of all strings S i with T , without computing the part of Y over and over again.

More generally, the common substring Y could be shared by different source strings competing over similarity with a common target, or could appear repeatedly in the same source string. Also, in a given application, we could of course be dealing with more than one repeated or shared subcomponent. It is assumed that the locations of the common subsequence Y in each source sequence S i are known. However, the part of the target T with which Y aligns, may vary according to each B i and F i combination.

We refer the interested reader to [START_REF] Landau | On the Common Substring Alignment Problem[END_REF] for a description of various applications which can be cast as CSA problems. In this manuscript, we will survey the CSA algorithms as macros in a general toolkit for utilizing various forms of repetitions in the graphs used for aligning two sequences.

Let n be the size of the compared sequences, Y denote the common substring, and let ℓ denote the size of Y . Using the classical dynamic programming tables, each appearance of Y in a source string would require the computation of all the values in a dynamic programming table of size O(nℓ).

The three solutions surveyed here are each composed of an encoding stage and an alignment stage. During the first stage, a data structure is constructed which encodes the comparison of Y with T . Then, during the alignment stage, for each comparison of a source S i with T , the pre-compiled data structure is used to speed up the part of Y .

A clear distinction should be made between the off-line pre-processing work and the online encoding stage. In some of the applications for which the CSA algorithms are intended, the source sequence database is prepared off-line, while the target can be viewed as an "unknown" sequence which is received online. The source strings can be pre-processed off-line and parsed into their optimal common substring representation.(The CSA algorithms of Stage One assume that we arrive to the problem after the strings have already been parsed, as opposed to the Stage Two algorithms which will handle the parsing job as well.) Therefore, the CSA algorithms know well beforehand where, in each S i , Y begins and ends. However, the comparison of Y and T cannot be computed until the target is received. Therefore, the encoding stage, as well as the alignment stage -are both online stages, and the tradeoff between the two must be cleverly minimized in order to maximize the efficiency gain by the suggested two-stage scenario. Note that even though both stages are online, they do not bear an equal weight on the time complexity of the algorithm. The efficiency gain is based on the fact that the first stage is executed only once per target, and then the encoding results are used, during the second stage, to speed up the alignment of each appearance of the common subcomponent in any of the source strings. In Section 3 the solutions will be analyzed in terms of encoding stage complexity and alignment stage complexity for a single common substring Y . Clearly, the additional running time depends only on the length of the parts of the strings that are not in any common substring.

The results presented in Section 3 include the following three solutions to the Common Substring Alignment Problem, which appear in [START_REF] Landau | On the Common Substring Alignment Problem[END_REF] and [START_REF] Landau | Sparse LCS Common Substring Alignment[END_REF].

1. Algorithm CSA1: the basic CSA algorithm for unrestricted scoring matrices.

This solution, which is presented in Section 3.1, reduces the O(nℓ) alignment work, for each appearance of the common substring Y in a source string, to O(n) -at the cost of an O(n 2 ) time encoding stage, which is executed only once. This algorithm is general enough to support scoring schemes which use a scoring matrix whose values are real numbers.

2. Algorithm CSA2: a more efficient CSA algorithm for discrete scoring schemes. This solution, which described in Section 3.2, consists of an O(n) work alignment stage, which is run for each appearance of the common substring Y in a source string, at the cost of O(nℓ) encoding work, which is executed only once. This algorithm is intended for all scoring schemes which use a scoring matrix whose values are rational numbers.

3. Algorithm CSA3: an even more efficient CSA algorithm for the LCS metric.

Let L y be the length of the LCS of T and Y , denoted |LCS[T, Y ]|, and L be max{|LCS[T, S i ]|}. This solution consists of an O(nL y ) time encoding stage that is executed once per common substring, and an O(L) time alignment stage that is executed once for each appearance of the common substring in each source string. This algorithm, which is given in Section 3.3, is intended for the LCS scoring scheme.

Stage

Two: Utilizing Repetitions in Both the Source Strings and the Target Strings.

In Section 4 of the manuscript we further advance by utilizing repetitions in both the source and the target strings. Here, two parsing methods, which are associated with known string compression algorithms, are used for defining the repetitions in the compared sequences. One is the parsing associated with Lempel-Ziv Text Compression, and the other is Run-Length Encoding. Here, the compression parsing is employed, not for the traditional objective of saving storage space or bandwidth, but rather in order to speed up the alignment process. The advantage of compressing both strings is that now the two stages (encoding and alignment) can be interleaved, as the blocks of the compression-partitioned dynamic programming graph are traversed in left-to-right, bottom down order and the output border values of each block are computed from the input border values. Due to the incremental properties of the LZ-partitioned blocks and the simplicity of the run-length blocks, the encoding time bottleneck is eliminated, and the solutions suggested in Stage 1 can be extended to efficiently apply to more advanced scoring schemes.

This leads to the following results, which appear in [START_REF] Crochemore | A Sub-quadratic Sequence Alignment Algorithm for Generalized Cost Metrics[END_REF] and are surveyed in Section 4.1.

A Sub Quadratic String Comparison algorithm for Unrestricted Scoring Matrices.

The classical algorithm for computing the similarity between two sequences, under string comparison metrics which use a scoring matrix of real number values [START_REF] Sankoff | Time Warps, String Edits, and Macromolecules: the Theory and Practice of Sequence Comparison[END_REF][START_REF] Smith | Identification of common molecular subsequences[END_REF], uses a dynamic programming matrix, and compares two strings of size n in O(n 2 ) time. The only previously known sub-quadratic global alignment string comparison algorithm, by Masek and Paterson [START_REF] Masek | A faster algorithm for computing string edit distances[END_REF], is based on the Four Russians paradigm. The "Four Russians" algorithm divides the dynamic programming table into uniform sized (log n by log n) blocks, and uses table lookup to obtain an O(n 2 / log n) time complexity string comparison algorithm, based on two assumptions. One is that the sequence elements come from a constant alphabet. The other, which they denote the "discreteness" condition, is that the weights (of substitutions and indels) are all rational numbers.

The algorithms described in [START_REF] Crochemore | A Sub-quadratic Sequence Alignment Algorithm for Generalized Cost Metrics[END_REF] present a new approach and are better than the above algorithm in two aspects. First, these algorithms are faster for compressible sequences. For such sequences, the complexity of the [START_REF] Crochemore | A Sub-quadratic Sequence Alignment Algorithm for Generalized Cost Metrics[END_REF] algorithms is O(hn 2 / log n), where h ≤ 1 is the entropy of the sequence.

Second, [START_REF] Crochemore | A Sub-quadratic Sequence Alignment Algorithm for Generalized Cost Metrics[END_REF] describes an algorithm that is general enough to support scoring schemes with real number weights. For many scoring schemes, the rational number weights supported by Masek and Paterson's algorithm do not suffice. For example, the entries of PAM similarity matrices, as well as BLOSUM evolutionary distance matrices, are defined to be real numbers, computed as log-odds ratios -and therefore could be irrational.

The paper by Masek and Paterson concludes with the following statement: "The most important problem remaining is finding a better algorithm for the finite (in our terms constant) alphabet case without the discreteness condition". In [START_REF] Crochemore | A Sub-quadratic Sequence Alignment Algorithm for Generalized Cost Metrics[END_REF], more than twenty years later, this important open question has finally been answered! These advantages are based on the following facts. First, the new algorithm does not require any pre-computation of lookup-tables, and therefore can afford more flexible weight values. Also, instead of dividing the dynamic programming matrix into uniform-sized blocks as did Masek and Paterson, the algorithms described in [START_REF] Crochemore | A Sub-quadratic Sequence Alignment Algorithm for Generalized Cost Metrics[END_REF] employ a variable-sized block partition, as induced by Lempel-Ziv factorization of both source and target. The common denominator between blocks, maximized by the compression technique, is then re-cycled and used for computing the relevant information for each block in time which is linear with the length of its sides. In this sense, the approach surveyed in section 4.1 can be viewed as another example of speeding up dynamic programming by keeping and computing only a relevant subset of important values, as demonstrated in [START_REF] Eppstein | Sparse Dynamic Programming I: Linear Cost Functions[END_REF], [START_REF] Eppstein | Sparse Dynamic Programming II: Convex and Concave Cost Functions[END_REF], [START_REF] Landau | On the Common Substring Alignment Problem[END_REF] and [START_REF] Schmidt | All Highest Scoring Paths In Weighted Grid Graphs and Their Application To Finding All Approximate Repeats In Strings[END_REF]. A similar unbalanced strategy has been successfully used for square detection in strings [START_REF] Crochemore | Transducers and Repetitions[END_REF] to speed up the original algorithm based on a divide-and-conquer approach [START_REF] Main | An O(n log n) algorithm for finding all repetitions in a string[END_REF].

Note that the notion of "acceleration by text-compression" has been recently applied to a related problem -that of exact string matching [START_REF] Karkkainen | Lempel-Ziv parsing and sublinear-size index structures for string matching[END_REF], [START_REF] Manber | A text compression scheme that allows fast searching directly in the compressed file[END_REF], [START_REF] Shabita | Speeding up pattern matching by text compression[END_REF]. It should also be mentioned that another related problem -that of exact string matching in compressed text without decoding it, which is often referred to as "compressed pattern matching", has been studied extensively [START_REF] Amir | Let sleeping files lie: Pattern matching in Z-compressed files[END_REF], [21] [55]. Along these lines, string search in compressed text was developed for the compression paradigm of LZ78 [START_REF] Ziv | Compression of individual sequences via variable rate coding[END_REF], and its subsequent variant LZW [START_REF] Welch | A Technique for High Performance Data Compression[END_REF], as described in [START_REF] Kida | Shift-And approach to pattern matching in LZW compressed text[END_REF], [START_REF] Navarro | Boyer-Moore string matching over Ziv-Lempel compressed text[END_REF]. A more challenging problem is that of "fully compressed" pattern matching when both the pattern and text strings are compressed [START_REF] Gasieniec | Randomised efficient algorithms for compressed strings: the finger-print approach[END_REF], [START_REF] Gasieniec | Almost optimal fully LZW compressed pattern matching[END_REF]. For the LZ78-LZW paradigm, compressed matching has been extended and generalized to that of approximate pattern matching (finding all occurrences of a short sequence within a long one allowing up to k changes) in [START_REF] Karkkainen | Approximate String Matching over Ziv-Lempel Compressed Text[END_REF], [START_REF] Navarro | Faster Approximate String Matching Over Compressed Text[END_REF].

Here, the adaptable sub-quadratic string comparison algorithms presented in [START_REF] Crochemore | A Sub-quadratic Sequence Alignment Algorithm for Generalized Cost Metrics[END_REF] will be described in the context of employing algorithms from the CSA toolkit given in Section 3, as follows. Adapting Algorithm CSA1 to work with the LZ-partitioned blocks leads to an O(n 2 / log n) time and space algorithm for an input of constant alphabet size. For most texts, the time and space complexity is actually O(hn 2 / log n) where h ≤ 1 is the entropy of the text. This algorithm applies to string comparison metrics that utilize a scoring matrix of real value numbers.

Adapting Algorithm CSA2 to the LZ-parsing block partitioning framework leads to a reduction of the space complexity to O(h 2 n 2 /(log n) 2 ), while preserving the O(hn 2 / log n) time complexity. This algorithm applies to string comparison metrics that utilize a scoring matrix of rational value numbers.

Note that, for the sake of presentation simplicity, we assume, throughout the description of the global alignment and the local alignment solutions (subsections 4.2 and 4.3), that both input strings A and B are of the same size n, and that both sequences share the same entropy h. For the case of comparing string A of size m and entropy 0 < h A ≤ 1 with string B of size n and entropy 0 < h B ≤ 1, the results of subsections 4.2 and 4.3 are as follows.

• O(mn(h A / log m + h B / log n)) time and space complexity for both global alignment and local alignment, replaces the O(hn 2 / log n) result.

• O(h A h B mn/ log m log n) space complexity for global alignment over "discrete" scoring matrices, replaces the O(h 2 n 2 /(log n) 2 ) result.

String Comparison of Run Length Compressed Strings.

A string S is run-length encoded if it is described as an ordered sequence of pairs (σ, i), often denoted "σ i ," each consisting of an alphabet symbol, σ, and an integer, i. Each pair corresponds to a run in S, consisting of i consecutive occurrences of σ. For example, the string aabbbbbccc can be encoded as a 2 b 5 c 3 . Such a run-length encoded string can be significantly shorter than the expanded string representation after efficiently encoding the integers (see [START_REF] Elias | Universal Codeword Sets and Representation of Integers[END_REF] for example). Run-length encoding serves as a popular image compression technique, since many classes of images (e.g., binary images in facsimile transmission or for use in optical character recognition) typically contain large patches of identically-valued pixels.

Let m and n be the lengths of two run-length encoded strings X and Y , of encoded lengths m ′ and n ′ , respectively. Previous algorithms for the problem compared two run-length encoded strings using the Levenshtein Edit Distance [START_REF] Levenshtein | Binary Codes Capable of Correcting Deletions, Insertions and Reversals[END_REF] and the LCS similarity measure [START_REF] Hirshberg | Algorithms for the longest common subsequence problem[END_REF]. For the LCS metric, Bunke and Csirik [START_REF] Bunke | An improved algorithm for computing the edit distance of run length coded strings[END_REF] presented an O(mn ′ + nm ′ ) time algorithm, while Apostolico, Landau, and Skiena [START_REF] Apostolico | Matching for Run Length Encoded Strings[END_REF] described an O(m ′ n ′ log(m ′ n ′ )) time algorithm. Mitchell [START_REF] Mitchell | A Geometric Shortest Path Problem, with Application to Computing a Longest Common Subsequence in Run-Length Encoded Strings[END_REF] has obtained an O((d Mäkinen, Navarro and Ukkonen [START_REF] Mäkinen | Approximate Matching of Run-Length Compressed Strings[END_REF] posed as an open problem the challenge of extending these results to more general scoring schemes, since in those applications which are related to image compression, the change from a pixel value to the next is smooth. This open problem was answered in [START_REF] Crochemore | A Sub-quadratic Sequence Alignment Algorithm for Generalized Cost Metrics[END_REF], where it was shown that adapting Algorithm CSA1 to the block partitioning obtained by run-length encoding yields an O(m ′ n + n ′ m) complexity algorithm for comparing two run-length encoded strings of length m and n, compressed into m ′ and n ′ runs respectively. This result, which is described in section 4.2, extends to all distance or similarity scoring schemes which use an additive gap penalty.

+ m ′ + n ′ ) log(d + m ′ + n ′ ))
We also observed that additional methods to achieve linear-time propagation across run-length compressed blocks can be obtained by employing algorithm CSA2 or alternatively by adapting the queue algorithms described in [START_REF] Eppstein | Speeding Up Dynamic Programming[END_REF][START_REF] Galil | Speeding Up Dynamic Programming with Applications to[END_REF][START_REF] Hirshberg | The Least Weight Subsequence Problem[END_REF]. In [START_REF] Mäkinen | Approximate Matching of Run-Length Compressed Strings[END_REF], Mäkinen et al. obtained the same bounds for weighted edit distance, using different techniques.

3 Utilizing Repetitions in Source Strings: Common Substring Alignment (CSA) Algorithms.

The alignment graph used for computing the distance between a source string S i = B i Y F i and a target string T can be viewed as a concatenation of three sub-graphs, where the first graph represents the distance between B i and T , the second graph represents the distance between Y and T , and the third graph represents the distance between F i and T (see Figure 3). Let ℓ denote the size of the subsequence Y which is common to all S i .

In this partitioned solution, the weights of the vertices in the last row of the first graph serve as input to initialize the weights of the vertices in the first row of the second graph. The weights of the last row of the second graph can be used to initialize the first row of the third graph. Therefore, let G denote the second sub-graph for comparing Y and T . Let I denote the series of weights of the vertices in the first row of G. Let O denote the series of weights of the vertices of the last row of G. The motivation for breaking the solution into three sub-graphs is that the second sub-graph, representing the distance between Y and T , is identical for all alignment graphs comparing any of the strings S i with T (see Figure 5). More specifically, both the structure and the weights of the edges of all alignment sub-graphs comparing Y with T are identical, but the weights to be assigned to the vertices during the distance computation may vary according to the prefix B i which is specific to the source string. Therefore, an initial investment in the learning of this graph as an encoding stage, and in its representation in a more informative data structure, may pay off later on.

The CSA solutions described in this section comply by the following 2-stage approach.

Encoding Stage: Given Y and T , encode G in a format which can be efficiently used by the alignment stage.

Alignment Stage: Given the output of the encoding stage and input row I -compute the output row O.

During the encoding stage, a data structure is constructed which encodes the comparison of Y with T . Then, during the alignment stage, for each comparison of a source S i with T , the pre-compiled data structure is used to speed up the part of aligning each appearance of the common substring Y .

3.1 Algorithm CSA1: An O(n 2 ) Encoding and O(n) Alignment CSA Algorithm for Unrestricted Scoring Matrices

The algorithm described in this subsection consists of an O(n 2 ) encoding stage and an O(n) alignment stage. Note that the time complexity of comparing a source string of size ℓ and a target string of size n, using the classical dynamic programming algorithm, is O(nℓ). Therefore, there are cases (n = Ω(ℓk) where k denotes the total number of appearances of the common substring Y in the database) where, due to the encoding stage bottleneck, it is not beneficial or even wasteful to employ this algorithm. In any case, this basic algorithm comes as a stepping stone in preparation for the more efficient two dimensional Common Substring Alignment algorithms, to be presented in section 4. In fact, in section 4 we will show that algorithm CSA1 does very well when both source and target strings are compressed, in which case the alignment stage bottleneck associated with this algorithm is eliminated.

The advantage of Algorithm CSA1 is in the generality of the scoring schemes to which it applies: it supports any additive gap scheme which utilizes a scoring matrix of real values.

The Encoding Stage

The following DIST matrix is computed (see Figure 8).

Definition 1 DIST [i, j], for j = 0 . . . n, i = 0 . . . j, stores the weight of the shortest path from the vertex in column i of the first row of the graph G to the vertex in column j of the last row of the graph G.

A similar encoding of an alignment graph has been used in [START_REF] Apostolico | Efficient parallel algorithms for string editing problems[END_REF][START_REF] Aggarawal | Notes on Searching in Multidimensional Monotone Arrays[END_REF][START_REF] Benson | A space efficient algorithm for finding the best nonoverlapping alignment score[END_REF][START_REF] Kannan | An Algorithm For Locating Non-Overlapping Regions of Maximum Alignment Score[END_REF][START_REF] Schmidt | All Highest Scoring Paths In Weighted Grid Graphs and Their Application To Finding All Approximate Repeats In Strings[END_REF].

A DIST matrix for an n × ℓ sized graph G (encodes the comparison of T versus Y ) can be constructed in O(n 2 + nℓ) time by using the algorithm of [START_REF] Schmidt | All Highest Scoring Paths In Weighted Grid Graphs and Their Application To Finding All Approximate Repeats In Strings[END_REF] (note that ℓ could be greater than n). For the LCS and Edit Distance metrics, DIST can also be constructed in O(n 2 + nℓ) time by employing [START_REF] Landau | Incremental String Comparison[END_REF]. Apostolico et al. [START_REF] Apostolico | New Clique and Independent Set Algorithms for Circle Graphs[END_REF] show how to construct DIST in O(n 2 ) time. Alternatively, one could use the algorithm from [START_REF] Apostolico | Efficient parallel algorithms for string editing problems[END_REF] to construct DIST in O(n 2 log n) time.

The Alignment Stage

Given input row I and sub-graph G, the weight of output row vertex O j can be computed as follows.

O j = j min r=0 {I r + DIST [r, j]}
This computation entails selecting a minimum among up to n sums for each of the n output sources.

(Figure 4 demonstrates an example of an output entry computation.) The above formulation was first presented in [START_REF] Kannan | An Algorithm For Locating Non-Overlapping Regions of Maximum Alignment Score[END_REF]. It is, in essence, a static version of the 1D dynamic programming problem [START_REF] Giancarlo | Dynamic Programming: Special Cases, Pattern Matching Algorithms[END_REF]:

E[j] = j min i=0 {D[i] + w(i, j)}
in which all values of D[i] are specified before any value of E[j] is computed, and the values of function w(i, j) are precomputed for all integers j = 0 . . . n, i = 0 . . . j.

Value O j is the minimum in column j of the following OU T matrix, which merges the information from input row I and DIST . (See Figure 8).

Definition 2 OU T [i, j] = I i + DIST [i, j] for j = 0 . . . n, i = 0 . . . j.
Aggarwal and Park [START_REF] Aggarawal | Notes on Searching in Multidimensional Monotone Arrays[END_REF] and Schmidt [START_REF] Schmidt | All Highest Scoring Paths In Weighted Grid Graphs and Their Application To Finding All Approximate Repeats In Strings[END_REF] observed that DIST matrices are Monge arrays [54].

Definition 3 A matrix M [0 . . . m, 0 . . . n]
is Monge if either condition 1 or 2 below holds for all i = 0 . . . m; j = 0 . . . n:

Condition 1: M [a, d] -M [a, c] ≥ M [b, d] -M [b, c] for all a < b and c < d. Condition 2: M [a, d] -M [a, c] ≤ M [b, d] -M [b, c] for all a < b and c < d.
We use the crossing paths contradiction [START_REF] Schmidt | All Highest Scoring Paths In Weighted Grid Graphs and Their Application To Finding All Approximate Repeats In Strings[END_REF] to show that our version of DIST obeys Monge Condition 1. Consider the four optimal (lowest scoring) paths from vertices a, b in the first row of G to vertices c, d in the last row of G, as shown in Figure 2. Note that P ath w and P ath z must intersect at some column of the Dynamic Programming Graph before or at column c. Let C denote the weight of the prefix of P ath w from its origin up to the intersection point, and let D denote the weight of the suffix of P ath w from the intersection point on. Let A denote the weight of the prefix of P ath z from its origin up to the intersection point, and B denote the weight of the suffix of P ath z from the intersection point on. Let |t|, |w|, |y| and |z| denote the weights of paths t, w, y and z, correspondingly.

Since path t is optimal, |t| ≤ C + B ⇒ |w| -|t| ≥ |w| -C -B = D -B. Since path y is optimal, |y| ≤ A + D ⇒ |y| -|z| ≤ A + D -|z| = D -B.
Therefore,

|w| -|t| = DIST [a, d] -DIST [a, c] ≥ DIST [b, d] -DIST [b, c] = |y| -|z|.
Since DIST is Monge by Condition 1, so is OU T , which is a DIST with constants added to its rows. An important property of Monge arrays is that of being totally monotone.

Definition 4 A matrix M [0 . . . m, 0 . . . n] is totally monotone if either condition 1 or 2 below holds for all a, b = 0 . . . m; c, d = 0 . . . n:

1. M [a, c] ≥ M [b, c] =⇒ M [a, d] ≥ M [b, d] for all a < b and c < d. 2. M [a, c] ≤ M [b, c] =⇒ M [a, d] ≤ M [b, d] for all a < b and c < d.
Note that the Monge property implies total monotonicity, but the converse is not true.

The recursive algorithm of Aggarwal et al [START_REF] Aggarwal | Geometric Applications of a Matrix-Searching Algorithm[END_REF] The second CSA algorithm comes as an answer to the following challenge: can the encoding stage complexity be reduced, without increasing the time complexity of the alignment stage? In order to make the CSA approach beneficial in all cases, our objective is to invest no more than O(nℓ) work in encoding an n × ℓ-sized common substring block. We present a more efficient, non-recursive algorithm which applies to all "discrete" scoring schemes, i.e. any additive gap scoring schemes which use a scoring matrix whose entries are rational numbers. Algorithm CSA2 obtains its strength advantage from the fact that it utilizes the Monge property, while the basic Algorithm CSA1 which was described in the previous subsection utilizes only the weaker Total Monotonicity property. This, combined with the discreteness properties of the allowed scoring schemes, yields a reduction of the encoding stage complexity from O(n 2 ) down to O(nℓ). Since O(nℓ) is the default time complexity for applying classical dynamic programming to any O(nℓ)-sized instance on the common sub-graph, this algorithm is a winner under any common substring parsing configuration.

The Encoding Stage

The encoding stage of Algorithm CSA2 gains its efficiency by utilizing the fact that the number of relevant changes, from one column of both OU T and DIST to the next, is constant [START_REF] Landau | Two Algorithms for LCS Consecutive Suffix Alignment[END_REF][START_REF] Landau | On the Common Substring Alignment Problem[END_REF][START_REF] Schmidt | All Highest Scoring Paths In Weighted Grid Graphs and Their Application To Finding All Approximate Repeats In Strings[END_REF]]. This property allows for a representation of DIST via an O(n) number of "relevant" points.

The DELT A matrix is defined as follows.

Definition 5 DELT A[i, j] = OU T [i, j] -OU T [i, j -1] = DIST [i, j] -DIST [i, j -1] for j = 1 . . . n, i = 0 . . . j -1.
The range of possible values for DELT A[i, j] depends on the scoring scheme which is used for the string comparison, and is actually the upper bound for the value difference between two consecutive elements in the dynamic programming table. (For an example of a DELT A matrix, see Figure 9.)

We will use the term ψ to denote the range bound for DELT A[i, j] values. As an example, if the similarity metric used is LCS, the only possible values for DELT A will be either 1 or 0, and ψ assumes a value of 1. For the Edit Distance metric, on the other hand, ψ is 2, since DELT A can only assume one of the 3 values: -1, 0, 1 [START_REF] Ukkonen | Finding Approximate Patterns in Strings[END_REF]. Our algorithm applies to all scoring scheme metrics for which ψ is a constant. It can be shown [START_REF] Masek | A faster algorithm for computing string edit distances[END_REF] that ψ is a constant for all scoring schemes which use a scoring matrix whose values are rational numbers.

Note that the following two observations apply to any column in DELT A.

Observation 1

Since DIST is a Monge array, each column in DELT A is a series of monotonically non-increasing values.

Observation 2

Since the range of distinct values which DELT A may assume is bounded by a constant (ψ), the number of "steps" (row indices in which the series of column entries increases in value) in each column of DELT A is constant.

As a result, DIST can be represented via an O(ψn) size set of relevant "step" points collected from all columns of DELT A (see Figure 9). Clearly, DELT A has up to ψn Borderline Points.

Observation 3

For any two rows x 1 , x 2 where x 1 ≤ Borderline[α, j] < x 2 for any α = 0 . . . ψ, for some column

j of DELT A, OU T [x 2 , j] -OU T [x 2 , j -1] < OU T [x 1 , j] -OU T [x 1 , j -1]
During the encoding stage, the O(nℓ) time complexity algorithm of [ [START_REF] Schmidt | All Highest Scoring Paths In Weighted Grid Graphs and Their Application To Finding All Approximate Repeats In Strings[END_REF], section 6] is used to compute the Borderline Points of DELT A. For the LCS metric, DELT A can also be constructed in O(nL y ) time, where L y = |LCS[T, Y ]| ≤ ℓ, by employing the algorithms of [START_REF] Landau | Two Algorithms for LCS Consecutive Suffix Alignment[END_REF].

The Alignment Stage

We now present an algorithm which uses the pre-compiled Borderline Points and input row I to compute output row O. The new algorithm computes the column minima of OU T , in left-to-right order, by executing n iterations. It employs the candidate list concept ( [START_REF] Eppstein | Speeding Up Dynamic Programming[END_REF], [START_REF] Galil | Speeding Up Dynamic Programming with Applications to[END_REF], [START_REF] Hirshberg | The Least Weight Subsequence Problem[END_REF], [START_REF] Miller | Sequence Comparison with Concave Weighting Functions[END_REF]), as follows.

Definition 7 The candidate list, at iteration j of the algorithm, is a subset of the rows of column j of OU T that includes only those rows which are candidates to contain future OU T column minima.

The candidate list is updated from one iteration of the alignment stage algorithm to the next. It is sorted by increasing OU T value and increasing row index. At iteration j of the alignment stage algorithm, the candidate of smallest row index in the list bears the minimal value at column j of OU T .

A high-level outline of the alignment stage algorithm is given below.

Procedure Alignment Stage input:

The set of Borderline Points for the DELT A matrix, and input row I. output: The output row O.

for j := 0 to |T | do 1 Append candidate j, which represents to row j of OU T , to the candidate list. 2 Update the contents of the candidate list. 3 Report output entry O j .

The list contents are updated at each iteration by removing rows which are no longer candidates to produce future column minima. We will denote such rows as extinct, according to total monotonicity condition 1.

Definition 8 A row x 1 is extinct at iteration j of the candidate list algorithm, for 0 < j ≤ n, if {∃ x 2 > x 1 | OU T [x 2 , j] ≤ OU T [x 1 , j]}.
Hence, for any distinct output value achieved during the computation of O j , we only need to keep one representative -that of the greatest candidate to achieve this value. All other candidates which generate the same output value in iteration j become extinct.

More specifically, a candidate becomes extinct, by Definition 8, as a result of two possible events. Event 1. A Bordeline[α, j]. Let x 1 , x 2 denote two rows which appear sequentially on the candidate list at iteration j, where x 1 ≤ Bordeline[α, j] < x 2 . For any two rows x 1 , x 2 where

x 1 ≤ Borderline[α, j] < x 2 . OU T [x 2 , j] -OU T [x 2 , j -1] < OU T [x 1 , j] -OU T [x 1 , j -1]
The candidate list is sorted by increasing OU T value and increasing row index. Therefore, by Observation 3, a Bordeline[α, j] could result in x 2 reaching an identical value to x 1 at column j of OU T . The difference in the OU T value between any two active candidates of index > Bordeline[α, j] is not affected by the Bordeline[α, j] event. The difference in the OU T value between any two active candidates of index ≤ Bordeline[α, j] is also unaffected by the Bordeline[α, j] event. The only two candidates which may now generate the same value at iteration j, as a direct result of the Bordeline[α, j] event, are the candidate x 1 which is of highest index below Bordeline[α, j], and the candidate x 2 which is of lowest row index above Bordeline[α, j], since the difference in the OU T values between these two candidates decreases by one. As a result, row x 1 will be removed from the list (row 5 in column 8 of the OU T matrix of Figure 10).

Event 2. At iteration j, row j is appended to the list. At this point, j is the highest row index in the list, and therefore all the elements with an OU T value which is higher than or equal to that of candidate j will be removed from the list (row 2 in column 3 and rows 1 and 3 in column 4 of the OU T matrix of Figure 10).

Note that two technical challenges need to be met, in order to implement a list manipulation engine, which updates the contents of the candidate list in linear time.

Computing the OU T value of a candidate.

The values of OU T are needed for two purposes. One is the comparison of two adjacent candidates. The other is for reporting the OU T value O j . Since OU T is quadratic in size, it is not possible to compute the OU T value for each candidate on the list in each column of OU T in linear time. Instead, Algorithm CSA2 keeps track of the differences only in OU T values between the members of the candidate list. Column j of the candidate list is computed incrementally from column j -1 of the candidate list by applying the relevant borderline points. The differences in OU T values between two candidates will only change in the positions updated by the borderline points. Since the number of borderline points per column is constant, in [START_REF] Landau | On the Common Substring Alignment Problem[END_REF] we are able to show how to compute all the necessary OU T values in linear time.

2. Efficiently accessing the rows to be removed from the candidate list.

Since not all rows of OU T appear in the candidate list, finding the row to be removed as a result of a Borderline Point is not trivial (see Event 1). However, this can be solved if the candidate list is implemented by employing the incremental tree set union algorithm described in [ [START_REF] Gabow | A Linear Time Algorithm for a Special Case of Disjoint Set Union[END_REF], pp. 216], for the special case in which the union tree is a path.

Time and Complexity Analysis of Algorithm CSA2.

In the encoding stage the Borderline Points are computed in O(nℓ) time and O(n) space using [START_REF] Schmidt | All Highest Scoring Paths In Weighted Grid Graphs and Their Application To Finding All Approximate Repeats In Strings[END_REF].

In [START_REF] Landau | On the Common Substring Alignment Problem[END_REF] it is proven that the alignment stage of Algorithm CSA2 computes output row O, from input row I and the Borderline Points for the comparison of Y with T , in O(n) time and space.

Note that during the encoding stage, the borderline points for the comparison of the prefix Y j 1 with T , can be incrementally computed in O(n) time from the borderline points for the comparison of Y j-1 1 with T , using [START_REF] Schmidt | All Highest Scoring Paths In Weighted Grid Graphs and Their Application To Finding All Approximate Repeats In Strings[END_REF]. Hence, for source sequences with two or more common factors, the time complexity of the encoding stage is further reduced to O(nD), where D is the number of nodes in the dictionary trie for the common factors.

Algorithm CSA3: An O(nL) Encoding and O(L) Alignment CSA Algorithm

for the LCS metric.

The third CSA algorithm comes as an answer to the following challenge: can a more efficient common substring alignment algorithm, which exploits the sparsity inherent to the LCS metric, be designed for the LCS metric? CSA3 exploits sparsity, based on the following two observations.

1. The LCS I and O vectors are both monotone staircases with at most L + 1 unit steps (see Figure 11).

Each row in the LCS

DIST is a monotone staircase with with at most L + 1 unit steps (see Figure 12).

In [START_REF] Landau | Sparse LCS Common Substring Alignment[END_REF] it is proven that the O(L) steps of I are sufficient for the computation of the O(L) steps of O. This allows the two vectors to be compressed into their "step" representation.

The Encoding Stage

In [START_REF] Landau | Sparse LCS Common Substring Alignment[END_REF] we show how the DIST matrix, which is traditionally used to encode the comparison of Y and T , can be replaced with a smaller matrix (see Figure 12), which is in essence the compressed, "step" representation of the rows of the original DIST . This "step" compressed DIST can be computed in O(nL y ) time using the Consecutive Suffix Alignments algorithms described in [START_REF] Landau | Two Algorithms for LCS Consecutive Suffix Alignment[END_REF].

The Alignment Stage

Algorithm CSA3 propagates the steps of I across the "step" compressed DIST matrix. In [START_REF] Landau | Sparse LCS Common Substring Alignment[END_REF] it is shown that the OU T matrix which can be obtained from combining the steps of I with the steps of DIST preserves the Total Monotonicity condition, thus enabling an efficient adaptation of the matrix searching algorithm of Aggarwal et al. 4 Utilizing Repetitions in Source Strings as well as Target Strings.

In this section we expand re-use Dynamic Programming to handle repetitions in both source and target strings. The repetitions are defined via two known parsing algorithms: one is the parsing associated with LZ78 compression. The other is run-length encoding. The results are obtained by adapting and expanding CSA algorithms from the toolkit given in the previous section.

In order to avoid confusion, we remind the reader that throughout this section the alignment problem will be re-formalized in terms of similarity maximization rather than distance minimization. (This, in spite of the fact that in the previous Section 3 algorithms CSA1 and CSA2 were presented using an edit-distance minimization formalism.) Therefore, in both subsections 4.1 and 4.2, we use the term "optimal path" to denote a highest scoring path.

A Block Partitioning of the Alignment Graph based on LZ78 Parsing

The classical algorithm for computing the similarity between two sequences [START_REF] Sankoff | Time Warps, String Edits, and Macromolecules: the Theory and Practice of Sequence Comparison[END_REF][START_REF] Smith | Identification of common molecular subsequences[END_REF] uses a dynamic programming matrix, and compares two strings of size in O(n 2 ) time. (For the sake of presentation simplicity, we assume, throughout this section, that both input strings A and B are of the same size n, and that both sequences share the same entropy h.) The speed-up is achieved by dividing the dynamic programming matrix into variable sized blocks, as induced by Lempel-Ziv parsing of both strings, and utilizing the inherent periodic nature of both strings. The compression mechanism enables the algorithm to adapt to the data and to utilize its repetitions. In this section we will describe an O(n 2 / log n) algorithm for computing the optimal global alignment value of two strings, of size n each, over a constant alphabet. The algorithm is even faster when the sequences are compressible. In fact, for most texts, the complexity of our algorithm is actually O(hn 2 / log n).

For global alignment over "discrete" scoring matrices, we explain how the space complexity can be reduced to O(h 2 n 2 /(log n) 2 ), without impairing the O(hn 2 / log n) time complexity.

The LZ compression methods are based on the idea of self reference: while the text file is scanned, substrings or phrases are identified and stored in a dictionary, and whenever, later in the process, a phrase or concatenation of phrases is encountered again, this is compactly encoded by suitable pointers [START_REF] Lempel | On the complexity of finite sequences[END_REF], [START_REF] Ziv | A Universal Algorithm for Sequential Data Compression[END_REF], [START_REF] Ziv | Compression of individual sequences via variable rate coding[END_REF]. Of the several existing versions of the method, we will use the ones which are denoted LZ78 family [START_REF] Welch | A Technique for High Performance Data Compression[END_REF], [START_REF] Ziv | Compression of individual sequences via variable rate coding[END_REF]. The main feature which distinguishes LZ78 factorization from previous LZ compression algorithms is in the choice of codewords. Instead of allowing pointers to reference any string that has appeared previously, the text seen so far is parsed into phrases, where each phrase is the longest matching phrase seen previously plus one character. For example, the string "B = aacgacgat" is into five phrases: a, ac, g, acg, at. Each phrase is encoded as an index to its prefix, plus the extra character (see Figure 14). The new phrase is then added to the list of phrases that may be referenced.

Since each phrase is distinct from others, the following upper bound applies to the possible number of phrases obtained by LZ78 factorization.

Theorem 4.1.1 (Lempel and Ziv 1976 [42]) Given a sequence S of size n over a constant alphabet. The maximal number of distinct phrases in S is O( n log n ).

Even though the upper bound above applies to any possible sequence over a constant alphabet, it has been shown that in many cases we can do better than that. Intuitively, the LZ78 algorithm compresses the sequence because it is able to discover some repeated patterns. Therefore, in order to compute a tighter upper bound on the number of phrases obtained by LZ78 factorization for "compressible" sequences, the repetitive nature of the sequence needs to be quantified. One of the fundamental ideas in information theory is that of entropy, denoted by the real number h, 0 < h ≤ 1, which measures the amount of disorder or randomness, or inversely, the amount of order or redundancy in a sequence. Entropy is small when there is a lot of order (i.e, the sequence is repetitive) and large when there is a lot of disorder. The entropy of a sequence should ideally reflect the ratio between the size of the sequence after it has been compressed, and the length of the uncompressed sequence. The number of distinct phrases obtained by LZ78 factorization has been shown to be O(hn/ log n) for most texts [START_REF] Bell | Text Compression[END_REF], [START_REF] Crochemore | Text Algorithms[END_REF], [START_REF] Lempel | On the complexity of finite sequences[END_REF], [START_REF] Ziv | Compression of individual sequences via variable rate coding[END_REF]. Note that for any text over a constant alphabet, the upper bound above still applies by setting h to 1.

The alignment graph will be partitioned as follows. Strings A and B will be parsed separately using LZ78 factorization. This induces a partition of the alignment graph for comparing A with B into variable-sized blocks (see Figure 7). Each block will correspond to a comparison of an LZ phrase of A with an LZ phrase of B.

Let xa denote a phrase in A obtained by extending a previous phrase x of A with character a, and yb denote a phrase in B, obtained by extending a previous phrase of B with character b.

From now on we will focus on the computations necessary for a single block of the alignment graph. We define the following three pref ix blocks of G.

1. The lef t pref ix of G denotes the block comparing phrase xa of A and phrase y of B.

2. The diagonal pref ix of G denotes the block comparing phrase x of A and phrase y of B.

3. The top pref ix of G denotes the block comparing phrase x of A and phrase yb of B.

Observation 4

When traversing the blocks of an LZ78 parsed alignment graph in a left-to-right, top-to-bottom order, the blocks for the left prefix, diagonal prefix and top prefix of G are encountered prior to block G (see Figure 7). Note that the graph for the left prefix of G is identical to the subgraph of G containing all columns but the last one. More specifically, both the structure and the weights of edges of these two graphs are identical, but the weights to be assigned to the vertices during the similarity computation may vary according to the input border values. Similarly, for the top prefix and diagonal prefix graphs. The only new cell in G, which does not appear in any of its prefix block graphs, is the cell in the bottom-rightmost corner of G, which is the cell for comparing a and b.

4.2 Combining Algorithm CSA1 with LZ78 Block Partitioning.

Similarly to Algorithm CSA1, the work for each block consists of two stages.

The Encoding Stage

The structure of G is encoded by computing a DIST matrix containing weights of optimal paths connecting each entry of its input border with each entry of its output border. Here we will use the incremental properties of LZ78 parsing to implement the encoding for each block in time linear with the size of its borders. Assuming that a Prefix Block can be accessed in constant time, the new DIST information for G can be computed in O(t) time, by using the DIST representations previously computed for its prefix blocks, plus the information of its new cell, as follows (see Figure 13). Only one new DIST column, column ℓ c , needs to be computed for the encoding of G, and this new column contains the weights of optimal paths from each input entry to the single new vertex of G, which is vertex (ℓ r , ℓ c ) in the lowest, rightmost corner of G. All the previous DIST columns are obtained by pointers to left prefix blocks of G. All the following DIST columns are obtained by pointers to top prefix blocks of G.

Furthermore, the new column of DIST is computed in O(t) time by using the 3 new edges in the new cell of G and the 3 DIST columns computed for the top, diagonal and left prefix blocks of G as follows. Entry [i] in column ℓ c of the DIST for G contains the weight of the optimal path from entry i in the input border of G to vertex (ℓ r , ℓ c ). This path must go through one of the three vertices (ℓ r -1, ℓ c ), (ℓ r -1, ℓ c -1) or (ℓ r , ℓ c -1). Therefore, the weight of the optimal (highest scoring) path from entry i in the input border of G to (ℓ r , ℓ c ) is equal to the maximum among the following three values:

1 Entry [i] of column ℓ c -1 of the DIST for the left prefix of G, plus the weight of the horizontal edge leading into (ℓ r , ℓ c ).

2 Entry [i] of column ℓ c -1 of the DIST for the diagonal prefix of G, plus the weight of the diagonal edge leading into (ℓ r , ℓ c ).

3 Entry [i] of column ℓ c of the DIST for the top prefix of G, plus the weight of the vertical edge leading into (ℓ r , ℓ c ).

In order to demonstrate how a prefix block can be accessed in constant time, we refer the reader to Figure 7. Here we have the tries for the phrases of A and B, as formed during the LZ78 parsing of the strings. Each of the partitioned blocks in the alignment graph for comparing strings A and B corresponds to a unique node in the trie for A and a unique node in the trie for B. For example, block 5/4 corresponds to node 5 in the trie for A (phrase "ag") and to node 4 in the trie for B (phrase "acg"). The tries can be used to access direct prefix blocks in constant time. For example, to obtain the diagonal prefix block of G, we just have to backtrack up one edge in the trie for A (this leads to node 3), and up one edge in the trie for B (this leads to node 2), and this way we compute in constant time the id (3,2) of the diagonal prefix block of G. We refer the reader to [START_REF] Crochemore | A Sub-quadratic Sequence Alignment Algorithm for Generalized Cost Metrics[END_REF] for the engineering details of the encoding algorithm.

The Alignment Stage

Given I and the DIST for G which was computed in the encoding stage, the alignment stage can be implemented in O(t) time, by using a variation of the alignment stage of Algorithm CSA1, extended to support full propagation from the leftmost and upper boundaries to the bottom and rightmost boundaries of G [START_REF] Crochemore | A Sub-quadratic Sequence Alignment Algorithm for Generalized Cost Metrics[END_REF]. The number of blocks in the LZ78-partitioned alignment graph is equal to the number of phrases in A times number of phrases in B, and is therefore O(h 2 n 2 /(log n) 2 ), assuming that both sequences share the same entropy h. When computing the optimal global alignment value with scoring matrices that follow the "discreteness" condition, algorithm CSA2 from section 3.2 can be extended to support full propagation from the leftmost and upper boundaries to the bottom and right most boundaries of G. This extended propagation algorithm can then be used to compute the values of the global alignment O for G, given the I for G and a minimal encoding of the DIST for G. The advantage of this minimal encoding of DIST is that rather than saving an O(t) sized DIST column per block, we only need to save a constant number of values per block. The encoding for the new DIST column of each block can be computed and stored in constant time and space from the information stored for the left, diagonal and top prefix blocks of G, using the technique described in Section 6 of [START_REF] Schmidt | All Highest Scoring Paths In Weighted Grid Graphs and Their Application To Finding All Approximate Repeats In Strings[END_REF]. This reduces the space complexity to O(h 2 n 2 /(log n) 2 ), while preserving the O(hn 2 / log n) time complexity.

Combining Algorithm CSA1 with Run Length Encoding Block Partitioning

In this solution, the alignment graph is also partitioned into blocks, where each block here consists of two runs -one of X and one of Y . This results in the partition of the alignment graph into m ′ n ′ blocks (see Figure 6.) The suggested algorithm also propagates accumulated scores from the left and upper boundaries of each block, to its bottom and right boundaries.

The Encoding Stage

We claim that this algorithm does not require any encoding stage, as any value of the DIST used for encoding a run-length block can be computed in constant time during the propagation stage, as follows. Consider the block R for comparing the run α i of X with the run β j of Y . An edge in R could be assigned one of three possible weight values: D(diagonal), H(horizontal) and V (vertical).

Let ∆ h and ∆ w denote the difference in row index values and column index values respectively, between entry i on the input border of R, and entry j on the output border of R.

We show how to compute DIST[i, j] (which is the cost of the best scoring path from entry i in the input border of the block, to entry j in the output border of the block) in constant time, given ∆ h and ∆ w for the input and output entries, and the values D, H and V .

• H + V ≤ D. Clearly, an optimal (highest scoring) path from i to j can use all possible diagonal edges and only then the minimal number of remaining H and V edges necessary to reach j.

Therefore, DIST[i, j] obtains one of three values:

1. If ∆ w = ∆ h , then DIST[i, j] = D × ∆ h . 2. If ∆ w > ∆ h , then DIST[i, j] = D × ∆ h + H × (∆ w -∆ h ). 3. If ∆ w < ∆ h , then DIST[i, j] = D × ∆ w + V × (∆ h -∆ w ).
• H + V > D. In this case, an optimal path never uses any diagonal edge. The path includes only the minimal number of H edges, and the minimal number of V edges necessary to reach j from i. In this case, DIST

[i, j] = H × ∆ w + V × ∆ h .
Therefore, DIST[i, j] can be easily computed in constant time for any additive-gap scoring scheme which uses a scoring matrix of real numbers.

The Alignment Stage

Given I and the DIST for G which was computed in the encoding stage, the alignment stage can be implemented for each block in time linear with the size of its I/O borders, by using a variation of the alignment stage of Algorithm CSA1, extended to support full propagation from the leftmost and upper boundaries to the bottom and rightmost boundaries of R [START_REF] Crochemore | A Sub-quadratic Sequence Alignment Algorithm for Generalized Cost Metrics[END_REF].

Time and Space Analysis of the CSA1/LZ78 Combination.

The O vector for each block is computed using SMAWK. Vector I for block R can be easily obtained from the O vectors for the block above R and the block to its left, in time linear with the sides of R. Therefore, any value OUT[i, j] = I[i] + DIST[i, j] can be computed in constant time.

Since the work for each block is linear with the size of its I/O borders, the total time complexity is linear with the total size of the borders of the blocks, which is O(m ′ n + n ′ m).

Since all the relevant DIST entry values are computed "on the fly" and do not need to be stored, Hirshberg's method [START_REF] Hirschberg | A Linear Space Algorithm for Computing Maximal Common Subsequences[END_REF] can be applied to achieve an algorithm with a space complexity which is linear with the size of the uncompressed strings.

An Appendix Including All Figures %& 9$ %" 2 % % #" %$ '& #@ '& 9$ %" %" %" %$ '& %( " %@ %@ '& 9$ %$ %$ '& 9$ '& $ %@ %@ %@ '& %& 9$ %$ %" %$ & 9@ %@ %@ %@ '& 9$ %" 2 #" ( 9@ %@ %@ %@ %@ '& 9$ %" %$ ) 9@ %@ %@ %@ %@ %@ '& 9$ %" 0 9@ %@ %@ %@ %@ %@ %@ '& 9$ 1 %@ %@ %@ %@ %@ %@ %@ %@ '& OUT[x,j] = Ix + DIST[x,j] '0 9) %( %& %& %& %( %) '0 #@ #0 9) %( %( %( %) '0 21 " %@ %@ #0 9) %) %) '0 9) '0 $ %@ %@ %@ ') %) %( %( %& %( & 9@ %@ %@ %@ '( %& 9$ %" %$ ( 9@ %@ %@ %@ %@ ') %( %& %( ) 9@ %@ %@ %@ %@ %@ #0 9) %( 0 2@ %@ %@ %@ %@ %@ %@ %1 #0 1 %@ %@ %@ %@ %@ %@ %@ %@ 'A B ¤¥ 4 3 ¤¥ ¤4 C5 ¡ ¨7 B 0 9) %( %& %& %& 9$ %" %$

T = "BCBADBDCD" Y = "BCBD" S1 = "BC BCBD C" B1 = "BC" F1 = "C" S2 = "E BCBD DBCBD A" B2a = "E" F2a = "DBCBD A" B 2b = "EBCBD D" F 2b = "A"

  time algorithm for a more general string matching problem in run-length encoded strings, where d is the number of matches of compressed characters. Both Arbell et al.[START_REF] Arbell | Edit distance of run-length encoded strings[END_REF] and Mäkinen et al.[START_REF] Mäkinen | Approximate Matching of Run-Length Compressed Strings[END_REF] independently obtained an O(m ′ n+n ′ m) time algorithm for computing the edit distance between two run-length encoded strings for the Levenshtein distance metric.

P

  ath t -corresponds to DIST [a, c]. P ath w -corresponds to DIST [a, d]. P ath z -corresponds to DIST [b, c]. P ath y -corresponds to DIST [b, d].

Definition 6

 6 Let Borderline[α, j], for α = 0 . . . ψ and j = 0 . . . n, denote a row index of a "step" of size 1 in the series of monotonically non-increasing values of column j of DELT A. (A step of size k is represented by k different Borderline Points).

  This algorithm consists of an O(nL y ) encoding stage, and an O(L) alignment stage, where L y denotes the length of the LCS of T and Y , denoted by |LCS[T, Y ]|, and L denotes max{|LCS[T, S i ]|}.When the problem is sparse (L y << |Y |, L << n), the time bounds of Algorithm CSA3 are better than those of Algorithms CSA1 and CSA2. Even when the data is dense, the time bounds of this algorithm are not worse than the previous two algorithms.

  [START_REF] Aggarwal | Geometric Applications of a Matrix-Searching Algorithm[END_REF]. This adaptation is then used to compute the O(L) steps of O from the O(L) steps of I in O(L) time.

  Consider the block G which corresponds to the comparison of xa and yb. Extending the definitions of section 3.1, we define input border I as the left and top borders of G, and output border O as the bottom and right borders of G. (The node entries on the input border are numbered in a clockwise direction, and the node entries on the output border are numbered in a counter-clockwise direction.) Let ℓ r denote the number of rows in G, ℓ r = |xa|. Let ℓ c denote the number of columns in G, ℓ c = |yb|. Let t = ℓ r + ℓ c . Clearly, |I| = |O| = t.

4. 3

 3 Combining Algorithm CSA2 with LZ78 Block Partitioning.

Figure 1 :Figure 2 :Figure 3 : 6 Figure 4 :

 12364 Figure 1: An example of two different source strings S 1 , S 2 sharing a common substring Y , and a target string T .

5 YFigure 5 :Figure 6 :Figure 7 :

 5567 Figure 5: Examples of two dynamic programming graphs for the alignment of two source sequences S 3 = "BCBCBD" and S 4 = "CABCBDCBCBD", sharing the common substring Y = "BCBD", with target string T = "BCBADBDCD". Note that Y is repeated twice in S 4 . The target sequence is T = "BCBADBDCD".

Figure 8 :

 8 Figure 8: The DIST and OU T matrices which correspond to the sequences T = "DCBADBDC", Y = "DCBD" and the input row I for Bi = "CBA". This figure continues the example of Figures 3 and 4.

Figure 9 :Figure 10 :

 910 Figure 9: An example of a DELT A matrix and its Borderline Points. Note that for the Edit Distance metric, which is used in this example, ψ = 2, and therefore the number of Borderline Points for DELT A is bounded by 2n. This figure continues the example of Figures 3, 4 and 8.¢¡ ¤£ ¥£ ¢¦§ ©¢ ! !" !# !$ &% $ !# !" ! !" !# &$ ' $ !# !" " !" !# !$ &% ' !' !$ # # !# !$ !# &$ ' !' !' # # !" !" ! &" ' !' !' !' " ( ) " ' !' !' !' ' # !" " # ' !' !' !' ' !' $ # " $ ' !' !' !' ' !' !' % $ % ' !' !' !' ' !' !' ' % 0 21 23 54 ¤6 87 ¤9 2@ A0 ! ( ! B" C ¤9 ED EF E ¢D 5ËG 21 5 IH E9 5G 21 E6 5P ' !' !' !' !' ! Q" ' !' !' !' ' !' R TS U 5¦ ¢1 5F 5G 2F E¦ E6 E WV ¤G EP ¢6 YX D E9 2@ `G 1 5F E Ea ¤b c !c !c !c c !d !d !d &d e !e !e fd !S !S !S &g h i g !g &p p &q ¢4 ¤6 3 54 56 W7 ¤9 2@ ¥ p !g !S !d d !d i (h )i

Figure 11 :Figure 12 :

 1112 Figure 11: Due to the unit-step properties of LCS, both I and O are monotonically nondecreasing series, and their values go up by unit steps. Therefore, the compressed Isteps and Osteps vectors are constructed which store only the steps of I and O, correspondingly.

Figure 13 :Figure 14 :

 1314 Figure 13: [13] Computing a single entry in the new DIST column for the new cell of G in constant time. This figure continues Figures 7 and 14.

  , nicknamed SM AW K in the literature, can be used to compute in O(n) time all the column minima of an n × n totally monotone matrix, by querying only O(n) elements of the array. Hence, one could use SM AW K to compute the output row O by querying only O(n) elements of OU T . Clearly, if both the full DIST and all entries of I are available, then accessing an element of OU T is O(1) work.One obstacle which comes up during this implementation is that DIST is not rectangular. Only the values in the upper triangle are defined. However, this can be resolved by setting the undefined values in the lower triangle to ∞.

Time and Space Complexity Analysis of Algorithm CSA1. DIST is constructed during the encoding stage in O(n 2 ) time and space. The alignment stage is done in O(n) time and space, for each appearance of the common substring Y in a source sequence.

3.2 Algorithm CSA2: An O(nℓ) Encoding and O(n) Alignment CSA Algorithm for Discrete Scoring Schemes

  Time and Space Analysis of the CSA1/LZ78 Combination Algorithm.Assuming sequence size n and sequence entropy h ≤ 1. The LZ78 factorization algorithm parses the strings and constructs the tries for A and B in O(n) time. The resulting number of phrases in both A and B is O(hn/ log n). Since the work and space for each block (both encoding and propagation) is linear with the size of its I/O borders, the total time and space complexity is linear with the total size of the borders of the blocks. The block borders form O(hn/ log n) rows of size |B| each, and O(hn/ log n) columns of size |A| each, in the alignment graph. Therefore, the total time and space complexity is O(hn 2 / log n).

Acknowledgements.

Many thanks to Amihood Amir, Alberto Apostolico and Ilan Newman for their helpful comments.

; partially supported by NSF grant CCR-0104307, by the Israel Science Foundation grant 282/01, by the FIRST Foundation of the Israel Academy of Science and Humanities, by IBM Faculty Partnership Award, and by the NATO Science Programme. ; partially supported by the Aly Kaufman Post Doctoral Fellowship, by the Bar-Nir Bergreen Software Technology Center of Excellence, and by the NATO Science Programme.


Maxime Crochemore

Renaud V Erin

Direct construction of Compact Directed Acyclic Word Graphs

Keywords: pattern matching algorithm, su x automaton, DAWG, Compact DAWG, su x tree, index on text

The Directed Acyclic Word Graph (DAWG) is an e cient data structure to treat and analyze repetitions in a text, especially in DNA genomic sequences. Here, we consider the Compact Directed Acyclic Word Graph of a word. We give the rst direct algorithm to construct it. It runs in time linear in the length of the string on a xed alphabet. Our implementation requires half the memory space used by DAWGs.

Introduction

In the classical string-matching problem for a word w and a text T , we want to know if w occurs in T , i.e., if w is a factor of T . In many applications, the same text is queried several times. So, e cient solutions are based on data structures built on the text that serve as an index to look for any word w in T . The typical running of various implementations of the search is O(jwj) (on a xed alphabet).

Among the implementations, the su x tree [START_REF] Mccreight | A space-economical su x tree construction algorithm[END_REF]) is the most popular. Its size and construction time are linear in the length of the text. It has been studied and used extensively. Apostolico 2] lists over 40 references on it, and Manber and Myers 12] mention several others. Many variants have been developed, like su x arrays 12], PESTry 11], su x cactus 10], or su x binary search trees 9]. Besides, the su x trie, the non-compact version of the su x tree, has been re ned to the su x automaton (Directed Acyclic Word Graph, DAWG). This automaton is a good alternative to represent the whole set of factors of a text. It is the minimal automaton accepting this set. It has been fully exposed by Blumer 3] and Crochemore 7]. As for the su x tree, its construction and size is linear in the length of the text.

In the genome research eld, DNA sequences can be viewed as words over the alphabet fa; c; g; tg. They become subjects for linguistic and statistic analysis.

For this purpose, su x automata are useful data structures. Indeed, the structure is fast to compute and easy to use. Meanwhile, the length of sequences in databases grows rapidly and the bottleneck to using the above data structures is their size. Keeping the index in main memory is more and more di cult for large sequences. So, having a structure using as little space as possible is appreciable for its construction as well as for its utilization. Compression methods are of no use to reduce the memory space of such indexes because they eliminate the direct access to substrings. On the contrary, the Compact Directed Acyclic Word Graph (CDAWG) keeps the direct access while requiring less memory space. The structure has been introduced by Blumer et al. [START_REF] Blumer | Complete inverted les for e cient text retrieval and analysis[END_REF][START_REF] Blumer | Average sizes of su x trees and dawgs[END_REF]). The automaton is based on the concatenation of factors issued from a same context. This concatenation induces the deletion of all states of outdegree one and of their corresponding transitions, excepting terminal states. This saves 50% of memory space. At the same time, the reduction of the number of states (2=3 less) and transitions (about half less) makes the applications run faster. Both time and space are saved.

In this paper, we give an algorithm to build compact DAWGs. This direct construction avoids constructing the DAWG rst, which makes it suitable for the actual DNA sequences (more than 1:5 million nucleotides for some of them). The compact DAWG allows to apply standard treatment on sequences twice as long in reasonable time (a few minutes).

In Section 2 we recall the basic notions on DAWGs. Section 3 introduces the compact DAWG, also called compact su x automaton, with the bounds on its size. We show in Section 4 how to build the CDAWG from the DAWG in time linear in the size of this latter structure. The direct construction algorithm for the CDAWG is given in Section 5. A conclusion follows.

2 De nitions Let be a nonempty alphabet and the set of words over , with " as the empty word. If w is a word in , jwj denotes its length, w i its i th letter, and w i::j its factor (subword) w i w i+1 : : :w j . If w = xyz with x; y; z 2 , then x, y, and z denote some factors or subwords of w, x is a pre x of w, and z is a su x of w. S(x) denotes the set of all su xes of x and F(x) the set of its factors.

For an automaton, the tuple (p; a; q) denotes a transition of label a starting at p and ending at q. A roman letter is used for mono-letter transitions, a greek letter for multi-letter transitions. Moreover, (p; ] denotes a transition from p for which is a pre x of its label.

Here, we recall the de nition of the DAWG, and a theorem about its implementation and its size proved in 3] and 7].

De nition1. The Su x Automaton of a word x, denoted DAWG(x), is the minimal deterministic automaton (not necessarily complete) that accepts S(x), the ( nite) set of su xes of x.

For example, Figure 1 shows the DAWG of the word gtagtaaac. States which are double circled are terminal states.

Theorem 2. The size of the DAWG of a word x is O(jxj) and the automaton can be computed in time O(jxj). The maximum number of states of the automaton is 2jxj 1, and the maximum number of edges is 3jxj 4. Recall that the right context of a factor u of x is u 1 S(x). The syntactic congruence, denoted by S(x) , associated with S(x) is de ned, for x; u; v 2 , by: u S(x) v () u 1 S(x) = v 1 S(x). We call classes of factors the congruence classes of the relation S(x) . The longest word of a class of factors is called the representative of the class. States of DAWG(x) are exactly the classes of the relation S(x) . Since this automaton is not required to be complete, the class of words not occurring in x, corresponding to the empty right context, is not a state of DAWG(x).

Moreover, we induce a selection among the congruence classes that we call strict classes of factors of S(x) and that are de ned as follows: De nition3. Let u be a word of C, a class of factors of S(x) . If at least two letters a and b of exist such that ua and ub are factors of x, then we say that C is a strict class of factors of S(x) .

We also introduce the function endpos x : F(x) ! N, de ned, for every word u, by: endpos x (u) = minfjwj j w pre x of x and u su x of wg and the function length x de ned on states of DAWG(x) by: length x (p) = juj; with u representative of p: The word u also corresponds to the concatenated labels of transitions of the longest path from the initial state to p in DAWG(x). The transitions that belong to the spanning tree of longest paths from the initial state are called solid transitions. Equivalently, for each transition (p; a; q) we have the property: (p; a; q) is solid () length x (q) = length x (p) + 1: The function length x works as well for multi-letter transitions, just replacing 1 in the above equivalence by the length of the label of the transition. This extends the notion of solid transitions to multi-letter transitions: (p; ; q) is solid () length x (q) = length x (p) + j j: In addititon, we de ne the su x link for a state of DAWG(x) by: De nition4. Let p be a state of DAWG(x), di erent from the initial state, and let u a word of the equivalence class p. The su x link of p, denoted by s x (p), is the state q which representative v is the longest su x z of u such that u 6 S(x) z. Note that, consequently to this de nition, we have length x (q) < length x (p). Then, by iteration, su x links induce su x paths in DAWG(x), which is an important notion used by the construction algorithm. Indeed, as a consequence of the above inequality, the sequence (p; s x (p); s 2 x (p); :::) is nite and ends at the initial state of DAWG(x). This sequence is called the su x path of p.

3 Compact Directed Acyclic Word Graphs

De nition

The compression of DAWGs is based on the deletion of some states and their corresponding transitions. This is possible using multi-letter transitions and the selection of strict classes of factors de ned in the previous section (De nition 3). Thus, we de ne the Compact DAWG as follows.

De nition5. The Compact Directed Acyclic Word Graph of a word x, denoted by CDAWG(x), is the compaction of DAWG(x) obtained by keeping only states that are either terminal states or strict classes of factors according to S(x) , and by labeling transitions accordingly. Consequently to De nition 3, the strict classes of factors correspond to the states that have an outdegree greater than one. So, we can delete every state having outdegree one exactly, except terminal states. Note that initial and nal states are terminal states too, so they are not deleted. The construction of the DAWG of a word including some repetitions shows that many states have outdegree one only. For example, in Figure 1, the DAWG of the word gtagtaaac has 12 states, 7 of which have outdegree one; it has 18 transitions. Figure 2 displays the result after the deletion of these states, using multi-letter transitions. The resulting automaton has only 5 states and 11 edges.

According to experiments to construct DAWGs of biological DNA sequences, considering them as words over the alphabet = fa; c; g; tg, we got that more than 60% of states have an outdegree one. So, the deletion of these states is worth, it provides an important saving. The average analysis of the number of states and edges is done in 5] in a Bernouilly model of probability.

When a state p is deleted, the deletion of outgoing edges is realized by adding the label of the outgoing edge of the deleted state to the labels of its incoming edges. For example, let r, p and q be states linked by transitions (r; b; p) and (p; a; q). We replace the edges (r; b; p) and (p; a; q) by the edge (r; ba; q). By recursion, we extend this method to every multi-letter transition (r; ; p).

In the example (Figure 1), one can note that, inside the word gtagtaaac, occurrences of g are followed by ta, and those of t and gt by a. So, gta is the representative of state 3 and it is not necessary to create states for g and (gt or t). Then, we directly connect state I to state 3 with edges (I,gta,3) and (I,ta,3). States 1 and 2 are so deleted.

The su x links de ned on states of DAWGs remain valid when we reduce them to CDAWGs because of the next lemma. Lemma 6. If p is a state of CDAWG(x), then s x (p) is a state of CDAWG(x).

Size bounds

By Theorem 2 DAWG(x) is linear in jxj. As we shall see below (Section 3.3), labels of multi-letter transitions are implemented in constant space. So, the size of CDAWG(x) is also O(jxj). Meanwhile, as we delete many states and edges, we review the exact bounds on the number of states and edges of CDAWG(x).

They are respectively denoted by States(x) and Edges(x). Corollary7. Given x 2 , if jxj = 0, then States(x) = 1; if jxj = 1, then States(x) = 2; else jxj 2, then 2 States(x) jxj+ 1 and the upper bound is reached when x is in the form a jxj , where a 2 . Corollary8. Given x 2 , if jxj = 0, Edges(x) = 0; if jxj = 1, Edges(x) = 1; else jxj 2, then Edges(x) 2jxj 2 and this upper bound is reached when x is in the form a jxj 1 c, where a and c are two di erent letters of .

Implementation and Results

Transition matrices and adjacency lists are the classical implementations of automata. Their principal di erence lies in the implementation of transitions. The rst one gives a direct access to transitions, but requires O(States(x) card( )).

The second one stores only the exact number of transitions in memory, but needs O(log card( )) time to access them. When the size of the alphabet is big and the transition matrix is sparse, adjacency lists are preferable. Otherwise, like for genomic sequences, transition matrix is a better choice, as shown by the experiments below. So, we only consider here transition matrices to implement CDAWGs.

We now describe the exact implementation of states and edges. We do this on a four-letter alphabet, so characters take 0:25 byte. We use integers encoded with 4 bytes. For each state, to encode the target state of outgoing edges, transitions matrices need a vector of 4 integers. Adjacency lists need, for each edge, 2 integers, one for the target state and another one for the pointer to the next edge.

The basic information required to construct the DAWG is composed of a table to implement the function s x and one boolean value (0:125 byte) for each edge to know if it is solid or not. For the CDAWG, in order to implement multiletter transitions, we need one integer for the endpos x value of each state, and another integer for the label length of each edge. And that is all.

Indeed, we can nd the label of a transition by cutting o the length of this transition from the endpos x value of its ending state. Then, we got the position of the label in the source and its length. Keeping the source in memory is negligible considering the global size of the automaton (0:25 byte by character). This is quite a convenient solution also used for su x trees. Figure 3 Then, respectively for transitions matrices and adjacency lists, each state requires 20:5 and 17:13 bytes for the DAWG, and 40:5 and 41:21 bytes for the CDAWG. As a reference, su x trees, as implemented by McCreight 13], need 28:25 and 20:25 bytes per state. Moreover, for CDAWG and su x trees the source has to be stored in main memory. Theoretical average numbers of states, calculated by Blumer et al. [START_REF] Blumer | Average sizes of su x trees and dawgs[END_REF]), are 0:54n for CDAWG, 1:62n for DAWG, and 1; 62n for su x trees, when n is the length of x. This gives respective sizes in bytes per character of the source: 45:68 and 32:70 for su x trees, 33:26 and 27:80 for DAWGs, and 22:40 and 22:78 for CDAWGs.

Considering the complete data structures required for applications, the function endpos x has to be added for the DAWG and the su x tree. In addition, the occurrence number of each factor has to be stored in each state for all the structures. Therefore, the respective sizes in bytes per character of the source become : 58:66 and 45:68 for su x trees, 46:24 and 40:78 for DAWGs, and 24:26 and 24:72 for CDAWGs. Indeed, we just need to apply the de nition of the CDAWG recursively. This is computed by the function Reduction, given below. Observe that, in this function, state(p; a] denotes the state pointed to by the transition (p; a]. The computation is done with a depth-rst traversal of the automaton, and runs in time linear in the number of transitions of DAWG(x). Then, by theorem 2, the computation also runs in time linear in the length of the text.

However, this method needs to construct the DAWG rst, which spends time and memory space proportional to DAWG(x), though CDAWG(x) is signi cantly smaller. So, it is better to construct the CDAWG directly.

Reduction ( Return (E,1);

Direct Construction of CDAWG

In this section, we give the direct construction of CDAWGs and show that the running time is linear in the size of the input word x on a xed alphabet.

Algorithm

Since the CDAWG of x is a minimization of its su x tree, it is rather natural to base the direct construction on McCreight's algorithm 13]. Meanwhile, properties of the DAWG construction are also used, especially su x links (notion that is di erent from the su x links of McCreight's algorithm), lengths, and positions, as explained in the previous section.

First, we introduce the notions used by the algorithm, some of them are taken from 13]. The algorithm constructs the CDAWG of the word x of length n, noted x 0::n 1 . The automaton is de ned by a set of states and transitions, especially with I and F, the initial and nal states. A partial path represents a connected sequence of edges between two states of the automaton. A path is a partial path that begins at I. The label of a path is the concatenation of the labels of corresponding edges.

The locus, or exact locus, of a string is the end of the path labeled by the string. The contracted locus of a string is the locus of the longest pre x of whose locus is de ned.

Preliminary Algorithm Basically, the algorithm to build CDAWG inserts the paths corresponding to all the su xes of x from the longest to the shortest. We de ne suf i as the su x x i::n 1 of x. We denote by A i the automaton constructed after the insertion of all the suf j for 0 j i. which are used subsequently. We initialize the automaton A " with states I and F. At step i (i > 0), the algorithm inserts a path corresponding to suf i in A i 1 and produces A i . The algorithm satis es the following invariant properties: P1: at the beginning of step i, all su xes suf j , 0 j < i, are paths in A i 1 . P2: at the beginning of step i, the states of A i 1 are in one-to-one correspondence with the longest common pre xes of pairs of su xes longer than suf j . We de ne head i as the longest pre x of suf i which is also a pre x of suf j for some j < i. Equivalently, head i is the longest pre x of suf i which is also a path of A i 1 . We de ne tail i as head 1 i suf i . At step i, the preliminary algorithm has to insert tail i from the locus of head i in A i 1 (see Figure 5).

To do so, the contracted locus of head i in A i 1 is found with the help of function SlowFind that compares letter-to-letter the right path of A i 1 to suf i . This is similar to the corresponding McCreight's procedure, except on what is explained below. Then, if necessary, a new state is created to split the last encountered edge, state that is the locus of head i . The automaton B of Figure 4, displays the creation of state 1 during the insertion of suf 1 =abbabbc. Note that, if an already existing state matches the strict class of factor of head i , the last encountered edge is split in the same way, but it is redirected to this state. Such an example appears in the same example (case D): the insertion of suf 5 =bbc induces the redirection of the edge (2,babbc,F) that becomes (2,b,3). Then, an edge labeled by tail i is created from the locus of head i to F. We can write the preliminary algorithm as follows:

Preliminary Algorithm 1. For all suf i (i 2 0..n-1]) Do 2.

(q; ) SlowFind(I); 3.

If ( = ") Then 4.
insert (q,taili,F);

create v locus of headi splitting (q; ] and insert (v,taili,F); or redirect (q; ] onto v, the last created state;

7. End For all;

8. mark terminal states;

Note rst that SlowFind returns the last encountered state. This keeps accessible the transition (q; ] that can be split if this state is not an exact locus.

Second, as in the DAWG construction, if a non-solid edge is encountered during SlowFind, its target state has to be duplicated in a clone and the nonsolid edge is redirected to this clone. But, if the clone has just been created at the previous step, the edge is redirected to this state. Note that, in the two cases, the redirected transition becomes solid.

Finally, when tail i = " at the end of the construction, terminal states are marked along the su x path of F.

From the above discussion, a proof of the invariance of properties P1 and P2 can be derived. Thus, at the end of the algorithm all subwords of x and only these words are labels of paths in the automaton (property P1). By property P2, states correspond to strict classes of factors (when the longest common pre x of a pair of su xes is not equal to any of them) or to terminal states (when the contrary holds). This gives a sketch of the correctness of the algorithm.

The running time of the preliminary algorithm is O(jxj 2 ) (with an implementation by transition matrix), like is the sum of lengths of all su xes of the word x.

Linear Algorithm To get a linear-time algorithm, we use together properties of DAWGs construction and of su x trees construction. The main feature is the notion of su x links. They are de ned as for DAWGs in Section 2. They are the clue for the linear-running-time of the algorithm.

Three elements have to be pointed out about su x links in the CDAWG. First, we do not need to initialize su x links. Indeed, when suf 0 is inserted, x 0 is obviously a new letter, which directly induces s x (F)=I. Note that s x (I) is never used, and so never de ned. Second, traveling along the su x path of a state p does not necessarily end at state I. Indeed, with multi-letter transitions, if s x (p)=I we have to treat the su x a 1 (a 2 ) where is the representative of p. And third, su x links induce the following invariant property satis ed at step i: P3: at the beginning of step i, the su x links are de ned for each state of A i 1 according to De nition 4. The next remark allows redirections without having to search with SlowFind for existing states belonging to a same class of factors.

The algorithm has to deal with su x links each time a state is created. This happens when a state is duplicated, and when a state is created after the execution of SlowFind.

In the duplication, su x links are updated as follows. Let w be the clone of q. In regard to strict classes of factors and De nition 4, the class of w is inserted between the ones of q and s x (q). So, we update su x links by setting s x (w)=s x (q) and s x (q)=w.

Moreover, the duplication has the same properties as in the DAWG construction. Let (p; ; q) be the transition redirected during the duplication of q. We can redirect all non-solid edges that end the partial path and that start from a state of the su x path of p. This is done until the rst edge that is solid. We are helped in this operation by the function FastFind, similar to the one used in McCreight's algorithm 13], that goes through transitions just comparing the rst letters of their labels. This function returns the last encountered state and edge. Note that it is not necessary to nd each time the partial path from a su x of p, we just need to take the su x link of the last encountered state and the label of the previous redirected transition.

Let # be the representative of a state of the su x path of p. Observe that the corresponding redirection is equivalent to insert suf i+j j j#j . Indeed, all operations done after this redirection will be the same as for the insertion of suf i , since they go through the same path. After the execution of SlowFind, if state v is created, we have to compute its su x link. Let be the label of the transition starting at q and ending at v. To compute the su x link, the algorithm goes through the path having label from the su x link of q, s = s x (q). The operation is repeated if necessary. Figure 6 displays a scheme of this search. The thick dashed edges represent paths in the automaton, and the thin dashed edge represents the su x link of q. This search will allow to insert, as for the duplication, the su xes suf j , for i < j < i+jhead i j.

To travel along the path, we use again the function FastFind. Let r and (r; ] be the last state and transition encountered by FastFind. If r is the exact locus of , it is the wanted state, and we set then s x (v) = r. Else, if (r; ] is a solid edge, then we have to create a new node w. The edge (r; ] is split, it becomes (r; ; w), and we insert the transition (w,tail i ,F). Else, (r; ] is non-solid. Then, it is split and becomes (r; ; v). In the two last cases, since s x (v) is not found, we run FastFind again with s x (r) and , and this goes on until s x (v) is eventually found, that is, when = ".

The discussion shows how su x links are updated to insure that property P3 is satis ed. The operations do not in uence the correctness of the algorithm, sketched in the last section, but yield the following linear-time algorithm. Its time complexity is discussed in the next section. Linear Algorithm 1. p I; i 0; 2. While not end of x Do 3.

(q; ) SlowFind(p); 4.

If ( = ") Then 5.

insert (q,taili,F); 6. sx(F) q; 7.

If (q 6 = I) Then p sx(q) Else p I;

Sketch of the proof

It can be proved that each step of the algorithm leads to increase strictly variables j or k in the generic situation displayed in Figure 7. These variables respectively represent the index of the current su x being inserted, and a pointer on the text. These variables never decrease. Therefore, the total running time of the algorithm is linear in the length of x.

Conclusion

We have considered the Compact Direct Acyclic Word Graph, which is an ecient compact data structure to represent all su xes of a word. There are many data structures representing this set. But, this one allows an interesting space gain compared to the well-known DAWG, which is a reference. Indeed, on the one hand, the upper bounds are of jxj + 1 states and 2jxj 2 transitions. This saves jxj states and jxj transitions of the DAWG, which leads to faster utilisation.

On the other hand, experiments on genomic DNA sequences and random strings display a memory space gain of 50% according to the DAWG. Moreover, when the size of the alphabet is small, transition matrices do not take more space than adjacency lists, keeping direct access to transitions. Thus, we can construct the data structure of twice larger strings, keeping them in main memory, which is actually important to get e cient treatments. This work shows that the CDAWG can be constructed directly. The algorithm is linear in the length of the text. Of course, it is easier to compute, by reduction, the CDAWG from the DAWG. On the contrary, our algorithm saves time and space simultaneously.

Fig. 1 .

 1 Fig. 1. DAWG(gtagtaaac)

Fig. 2 .

 2 Fig. 2. CDAWG(gtagtaaac)

Fig. 3 .

 3 Fig. 3. Data Structure of CDAWG(gtagtaaac)

Fig. 4 .Figure 4

 44 abbc c

Fig. 5 .

 5 Fig. 5. Scheme of the insertion of a suf i in Ai 1.

Fig. 6 .

 6 Fig. 6. Scheme of the search using su x links

  The algorithm that builds the CDAWG of a word x of can be implemented in time O(jxj) and in space O(jxj card( )) with a transition matrix, or in time O(jxj log card( )) and in space O(jxj) with adjacency lists.

	5.2 Complexity
	x Theorem9. I	i	head i j q	suf i	v	k	tail i
		I	s	r
		8. 9. 10. 11. 12. 13. 14. 15. End While; Else create v locus of headi splitting (q; ]; insert (v,taili,F); sx(F) v; nd r = sx(v) with FastFind; p r; update i; 16. mark terminal states;

Fig. 7. Positions of labels when suf i is inserted

Remark. Let have locus p and assume that q = s x (p) is the locus of . Then, p is the locus of su xes of whose lengths are greater than j j.

Moreover, Table 1 compares sizes of DAWG and CDAWG meant for applications to DNA sequences. Sizes for random words of di erent lengths and j j = 4 are also given. DNA sequences are Saccharomyces cerevisiae yeast chromosome II (chro II), a contig of Escherichia Coli DNA sequence (coli), and contigs 1 and 115 of Bacillus Subtilis DNA sequence (bs). Number of states and edges according to the length of the source and the memory space gain are displayed. Theoretical average ratios are given, calculated from Blumer et al. [START_REF] Blumer | Average sizes of su x trees and dawgs[END_REF]). First, we observe there are 2=3 less states in the CDAWG, and near of half edges. Second, the memory space saving is about 50%. Third, the number of edges by state is going up to 2:66. With a four-letter alphabet, this is interesting because the transition matrix becomes smaller than adjacency lists. At the same time, we keep a direct access to transitions.

Constructing CDAWG from DAWG

The DAWG construction is fully exposed and demonstrated in 3] and 7]. As we show in this section, the CDAWG is easily derived from the DAWG.