siEDM: An Efficient String Index and Search Algorithm for Edit Distance with Moves

Abstract

Although several self-indexes for highly repetitive text collections exist, developing an index and search algorithm with editing operations remains a challenge. Edit distance with moves (EDM) is a string-to-string distance measure that includes substring moves in addition to ordinal editing operations to turn one string into another. Although the problem of computing EDM is intractable, it has a wide range of potential applications, especially in approximate string retrieval. Despite the importance of computing EDM, there has been no efficient method for indexing and searching large text collections based on the EDM measure. We propose the first algorithm, named string index for edit distance with moves (siEDM), for indexing and searching strings with EDM. The siEDM algorithm builds an index structure by leveraging the idea behind the edit sensitive parsing (ESP), an efficient algorithm enabling approximately computing EDM with guarantees of upper and lower bounds for the exact EDM. siEDM efficiently prunes the space for searching query strings by the proposed method, which enables fast query searches with the same guarantee as ESP. We experimentally tested the ability of siEDM to index and search strings on benchmark datasets, and we showed siEDM’s efficiency.

1. Introduction

Vast amounts of text data are created, replicated, and modified with the increasing use of the internet and advances of data-centric technology. Many of these data contain repetitions of long substrings with slight differences, so called highly repetitive texts, such as Wikipedia and software repositories like GitHub with a large number of revisions. Also, recent biological databases store a large amount of human genomes while the genetic differences among individuals are less than $0.1$ percent, which results in the collections of human genomes being highly repetitive. Therefore, there is a strong need to develop powerful methods for processing highly repetitive text collections on a large scale.

Building indexes is the de facto standard method to search large databases of highly repetitive texts. Several methods have been presented for indexing and searching large-scale and highly repetitive text collections. Examples include the ESP-index [1], SLP-index [2] and LZ77-based index [3]. Recently, Gagie and Puglisi [4] presented a general framework called kernelization for indexing and searching highly repetitive texts. Although these methods enable fast query searches, their applicability is limited to exact match searches.

The edit distance between two strings is the minimum cost of edit operations (insertions, deletions, and replacements of characters) to transform one string to another. It has been proposed for detecting evolutionary changes in biological sequences [5], detecting typing errors in documents [6], and correcting errors on lossy communication channels [7]. To accelerate the quadratic time upper bound on computing the edit distance, Cormode and Muthukrishnan introduced a new technique called edit sensitive parsing (ESP) [8]. This technique allows us to compute a modified edit distance in near linear time by sacrificing accuracy with theoretical bounds. The modified distance is known as edit distance with moves (EDM) [8], which includes substring move operations in addition to insertions and deletions. While the exact computation of EDM is known to be intractable [9], the approximate computation of EDM using ESP achieves a good approximation ratio $O(lgN{lg}^{*}N)$, and runs in almost linear time $O\left(N{lg}^{*}N\right)$ for the string length N, where lg denotes the logarithm of base two.

ESP is extended to various applications for highly repetitive texts. Examples are data compressions called grammar compression [10,11,12,13], indexes for exact matches [1,14,15], an approximated frequent pattern discovery [16] and an online pattern matching for EDM [17]. Despite several attempts to efficiently compute EDM and various extensions of ESP, there is no method for indexing and searching texts with EDM. Such a method is required in bioinformatics where approximated text searches are used to analyze massive genome sequences. Thus, an open challenge is to develop an efficient string index and search algorithm for EDM.

We propose a novel method called siEDM that efficiently indexes massive text, and performs query searches for EDM. As far as we know, siEDM is the first string index for searching queries for EDM. A space-efficient index structure for a string is built by succinctly encoding a parse tree obtained from ESP, and query searches are performed on the encoded index structures. siEDM prunes useless portions of the search space based on the lower bound of EDM without missing any matching patterns, enabling fast query searches. As in existing methods, similarity searches of siEDM are approximate but have the same guarantee of the approximation ratio as in ESP.

Experiments were performed on indexing and searching repetitive texts for EDM on standard benchmark datasets. The performance comparison with an online pattern matching for EDM [17] demonstrates siEDM’s practicality.

2. Preliminaries

2.1. Basic Notations

Let Σ be a finite alphabet, and σ be $|\Sigma |$. All elements in Σ are totally ordered. Let us denote by ${\Sigma }^{*}$ the set of all strings over Σ, and by ${\Sigma }^q$ the set of strings of length q over Σ, i.e., ${\Sigma }^q=\{w\in {\Sigma }^{*}:\left|w\right|=q\}$ and an element in ${\Sigma }^q$ is called a q-gram. The length of a string S is denoted by $\left|S\right|$. The empty string ϵ is a string of length 0, namely $\left|ϵ\right|=0$. For a string $S=\alpha \beta \gamma$, α, β and γ are called the prefix, substring, and suffix of S, respectively. The i-th character of a string S is denoted by $S\left[i\right]$ for $i\in [1,|S\left|\right]$. For a string S and interval $[i,j]$ ($1\le i\le j\le \left|S\right|$), let $S[i,j]$ denote the substring of S that begins at position i and ends at position j, and let $S[i,j]$ be ϵ when $i>j$. For a string S and integer $q\ge 0$, let $\mathit{pre}(S,q)=S[1,q]$ and $\mathit{suf}(S,q)=S\left[\right|S|-q+1,|S\left|\right]$. We assume a recursive enumerable set $\mathcal{X}$ of variables with $\Sigma \cap \mathcal{X}=\varnothing$. All elements in $\Sigma \cup \mathcal{X}$ are totally ordered, where all elements in Σ must be smaller than those in $\mathcal{X}$. In this paper, we call a sequence of symbols from $\Sigma \cup X$ a string. Let us define ${lg}^{\left(1\right)}u=lgu$, and ${lg}^{(i+1)}u=lg({lg}^{\left(i\right)}u)$ for $i\ge 1$. The iterated logarithm of u is denoted by ${lg}^{*}u$, and defined as the number of times the logarithm function must be applied before the result is less than or equal to 1, i.e., ${lg}^{*}u=min\{i:{lg}^{\left(i\right)}u\le 1\}$.

2.2. Straight-Line Program (SLP)

A context-free grammar (CFG) in Chomsky normal form is a quadruple $G=(\Sigma ,V,D,X_s)$, where V is a finite subset of $\mathcal{X}$, D is a finite subset of $V\times {(V\cup \Sigma )}^2$, and $X_s\in V$ is the start symbol. An element in D is called a production rule. Denote $X_{l\left(k\right)}$ (resp. $X_{r\left(k\right)}$) as a left symbol (resp. right symbol) on the right hand side for a production rule with a variable $X_k$ on the left hand side, $i.e.$, $X_k\to X_{l\left(k\right)}{X}_{r\left(k\right)}$. $\mathit{val}\left(X_i\right)$ for variable $X_i\in V$ denotes the string derived from $X_i$. A grammar compression of S is a CFG G that derives S and only S. The size of a CFG is the number of variables, i.e., $\left|V\right|$ and let $n=\left|V\right|$.

The parse tree of G is a rooted ordered binary tree such that (i) internal nodes are labeled by variables in V and (ii) leaves are labeled by symbols in Σ, i.e., the label sequence in leaves is equal to the input string. In a parse tree, any internal node Z corresponds to a production rule $Z\to XY$, and has the left child with label X and the right child with label Y.

Straight-line program (SLP) [18] is defined as a grammar compression over $\Sigma \cup V$, and its production rules are in the form of $X_k\to X_i{X}_j$ where $X_i,X_j\in \Sigma \cup V$ and $1\le i,j<k\le n+\sigma$.

2.3. Rank/Select Dictionaries

A rank/select dictionary for a bit string B [19] supports the following queries: ${\mathit{rank}}_c(B,i)$ returns the number of occurrences of $c\in \{0,1\}$ in $B[0,i]$; ${\mathit{select}}_c(B,i)$ returns the position of the i-th occurrence of $c\in \{0,1\}$ in B; $\mathit{access}(B,i)$ returns the i-th bit in B. Data structures with only the $\left|B\right|+o\left(\right|B\left|\right)$ bits storage to achieve $O\left(1\right)$ time rank and select queries [20] have been presented.

GMR [21] is a rank/select dictionary for large alphabets and supports rank/ select/access queries for strings in ${(\Sigma \cup V)}^{*}$. GMR uses $(n+\sigma )lg(n+\sigma )+o\left(\right(n+\sigma )lg(n+\sigma ))$ bits while computing both rank and access queries in $O(lglg(n+\sigma ))$ times and also computing select queries in $O\left(1\right)$ time.

3. Problem

We first review the notion of EDM. The distance $d(S,Q)$ between two strings S and Q is the minimum number of edit operations to transform S into Q. The edit operations are defined as follows:

• Insertion: A character a is inserted at position i in S, which generates $S[1,i-1]aS[i,|S\left|\right]$,
• Deletion: A character is deleted at position i in S, which generates $S[1,i-1]S[i+1,|S\left|\right]$,
• Replacement: A character is replaced by a at position i in S, which generates $S[1,i-1]aS[i+1,|S\left|\right]$,
• Substring move: A substring $S[i,j]$ is deleted from the position i, and inserted at the position k in S, which generates $S[1,i-1]S[j+1,k-1]S[i,j]S[k,|S\left|\right]$ for $1\le i\le j\le k\le \left|S\right|$, and $S[1,k-1]S[i,j]S[k,i-1]S[j+1,|S\left|\right]$ for $1\le k\le i\le j\le \left|S\right|$.

Problem 1 

(Query search for EDM). For a string $S\in {\Sigma }^{*}$, a query $Q\in {\Sigma }^{*}$ and a distance threshold $\tau \ge 0$, find all $i\in [1,|S\left|\right]$ satisfying $d\left(S\right[i,i+\left|Q\right|-1],Q)\le \tau$.

Shapira and Storer [9] proved the NP-completeness of EDM and proposed a polynomial-time algorithm for a restricted EDM. Cormode and Muthukrishnan [8] presented an approximation algorithm named ESP for computing EDM. We present a string index and search algorithm by leveraging the idea behind ESP for solving Problem 1. Our method consists of two parts: (i) an efficient index structure for a given string S and (ii) a fast algorithm for searching query Q on the index structure of S with respect to EDM. Although our method is also an approximation algorithm, it guarantees upper and lower bounds for the exact EDM. We first review ESP in the next section and then discuss the two parts.

4. Edit Sensitive Parsing (ESP) for Building SLPs

4.1. ESP Revisit

We review the edit sensitive parsing algorithm for building SLPs [10]. This algorithm, referred to as ESP-comp, computes an SLP from an input sting S. The tasks of ESP-comp are to (i) partition S into $s_1{s}_2\cdots s_{\ell }$ such that $2\le |s_i|\le 3$ for each $1\le i\le \ell$, (ii) if $|s_i|=2$, generate the production rule $X\to s_i$ and replace $s_i$ by X (this subtree is referred to as a 2-tree), and if $|s_i|=3$, generate the production rule $Y\to AX$ and $X\to BC$ for $s_i=ABC$, and replace $s_i$ by Y (referred to as a 2-2-tree), (iii) iterate this process until S becomes a symbol. Finally, the ESP-comp builds an SLP representing the string S.

We focus on how to determine the partition $S=s_1{s}_2\cdots s_{\ell }$. A string of the form $a^r$ with $a\in \Sigma \cup V$ and $r\ge 2$ is called a repetition. First, S is uniquely partitioned into the form $w_1{x}_1{w}_2{x}_2\cdots w_k{x}_k{w}_{k+1}$ by its maximal repetitions, where each $x_i$ is a maximal repetition of a symbol in $\Sigma \cup V$, and each $w_i\in {(\Sigma \cup V)}^{*}$ contains no repetition. Then, each $x_i$ is called type1, each $w_i$ of length at least $2{lg}^{*}\left|S\right|$ is type2, and any remaining $w_i$ is type3. If $|w_i|=1$, this symbol is attached to $x_{i-1}$ or $x_i$ with preference $x_{i-1}$ when both cases are possible. Thus, if $\left|S\right|>2$, each $x_i$ and $w_i$ is longer than or equal to two. One of the substrings is referred to as $S_i$.

Next, ESP-comp parses each $S_i$ depending on the type. For type1 and type3 substrings, the algorithm performs the left aligned parsing as follows. If $|S_i|$ is even, the algorithm builds 2-tree from $S_i[2j-1,2j]$ for each $j\in \{1,2,\dots ,|S_i|/2\}$; otherwise, the algorithm builds a 2-tree from $S_i[2j-1,2j]$ for each $j\in \{1,2,\dots ,\lfloor (|S_i|-3)/2\rfloor \}$ and builds a 2-2-tree from the last trigram $S_i\left[\right|S_i|-2,|S_i\left|\right]$. For type2 $S_i$, the algorithm further partitions it into short substrings of length two or three by alphabet reduction [8].

Alphabet reduction: Given a type2 string S, consider $S\left[i\right]$ and $S[i-1]$ as binary integers. Let p be the position of the least significant bit, in which $S\left[i\right]\ne S[i-1]$, and let $bit(p,S[i\left]\right)$ be the bit of $S\left[i\right]$ at the p-th position. Then, $L\left[i\right]=2p+bit(p,S[i\left]\right)$ is defined for any $i\ge 2$. Because S is repetition-free (i.e., type2), the label string $L\left(S\right)=L\left[2\right]L\left[3\right]\cdots L\left[\right|S\left|\right]$ is also type2. If the number of different symbols in S is n (denoted by $\left[S\right]$), then $\left[L\right(S\left)\right]=O(lgn)$. For the $L\left(S\right)$, the next label string is iteratively computed until the final $L^{*}\left(S\right)$ satisfying $\left[L^{*}\left(S\right)\right]\le {lg}^{*}\left|S\right|$ is obtained. $S\left[i\right]$ is called the landmark if $L\left[i\right]>max\left\{L\right[i-1],L[i+1\left]\right\}$.

The alphabet reduction transforms S into $L^{*}\left(S\right)$ such that any substring of $L^{*}\left(S\right)$ of length at least $2{lg}^{*}\left|S\right|$ contains at least one landmark because $L^{*}\left(S\right)$ is also type2. Using this characteristic, the algorithm ESP-comp determines the bigrams $S[i,i+1]$ to be replaced for any landmark $S\left[i\right]$, where any two landmarks are not adjacent, and then the replacement is deterministic. After replacing all landmarks, any remaining maximal substring s is replaced by the left aligned parsing, where if $\left|s\right|$ = 1, it is attached to its left or right block.

We give an example of the edit sensitive parsing of an input string in Figure 1-(i) and (ii). The input string S is divided into three maximal substrings depending on the types. The label string L is computed for the type2 string. Originally, L is iteratively computed until $\left[L\right]\le {lg}^{*}\left|S\right|$. This case shows that a single iteration satisfies this condition. After the alphabet reduction, three landmarks $S\left[i\right]$ are found, and then each $S[i,i+1]$ is parsed. Any other remaining substrings including type1 and type3 are parsed by the left aligned parsing shown in Figure 1-(ii). In this example, a dashed node denotes that it is an intermediate node in a 2-2-tree. Originally, an ESP tree is a ternary tree in which each node has at most three children. The intermediate node is introduced to represent ESP tree as a binary tree.

As shown in [8], the alphabet reduction approximates the minimum CFG as follows. Let S be a type2 string containing a substring α at least twice. When α is sufficiently long (e.g., $\left|\alpha \right|\ge 2{lg}^{*}\left|S\right|$), there is a partition $\alpha ={\alpha }_1{\alpha }_2$ such that $|{\alpha }_1|=O({lg}^{*}\left|S\right|)$ and each landmark of ${\alpha }_2$ within α is decided by only ${\alpha }_1$. This means the long prefix ${\alpha }_2$ of α is replaced by the same variables, independent of the occurrence of α.

ESP-comp generates a new shorter string $S^{\prime }$ of length from $\left|S\right|/3$ to $\left|S\right|/2$, and it parses $S^{\prime }$ iteratively. Given a string S, ESP builds the ESP-tree of height $O(lg|S\left|\right)$ in $O\left(\right|S\left|{lg}^{*}\right|S\left|\right)$ time and in $O\left(\right|\Sigma \cup V|lg|\Sigma \cup V|)$ space. The approximation ratio of the smallest grammar by ESP is $O({lg}^{*}\left|S\right|lg\left|S\right|)$ [10].

4.2. Approximate Computations of EDM from ESP-Trees

ESP-trees enable us to approximately compute EDM for two strings. After constructing ESP-trees for two strings, their characteristic vectors are defined as follows. Let $T\left(S\right)$ be the ESP-tree for string S. We define that an integer vector $F\left(S\right)$ to be the characteristic vector if $F\left(S\right)\left(X\right)$ represents the number of times the variable X appears in $T\left(S\right)$ as the root of a 2-tree. For a string S, $T\left(S\right)$ and its characteristic vector are illustrated in Figure 2. The EDM between two strings S and Q can be approximated by $L_1$-distance between two characteristic vectors $F\left(S\right)$ and $F\left(Q\right)$ as follows:

$${\parallel F\left(S\right)-F\left(Q\right)\parallel }_1=\sum _{e\in V\left(S\right)\cup V\left(Q\right)}|F\left(S\right)\left(e\right)-F\left(Q\right)\left(e\right)|$$

Cormode and Muthukrishnan showed the upper and lower bounds on the $L_1$-distance between characteristic vectors for the exact EDM.

Theorem 1 

(Upper and lower bounds of the approximated EDM [8]) . For $N=max\left(\right|S|,|Q\left|\right)$,

$$d(S,Q)\le {2\parallel F\left(S\right)-F\left(Q\right)\parallel }_1=O(lgN{lg}^{*}N)d(S,Q)$$

5. Index Structure for ESP-Trees

5.1. Efficient Encoding Scheme

siEDM encodes an ESP-tree built from a string for fast query searches. This encoding scheme sorts the production rules in an ESP-tree such that the left symbols on the right hand side of the production rules are in monotonically increasing order, which enables encoding of these production rules efficiently and supporting fast operations for ESP-trees. The encoding scheme is performed from the first and second levels to the top level (i.e., root) in an ESP-tree.

First, the set of production rules at the first and second levels in the ESP-tree is sorted in increasing order of the left symbols on the right hand of the production rules, i.e., $X_{l\left(i\right)}$ in the form of $X_i\to X_{l\left(i\right)}{X}_{r\left(i\right)}$, which results in a sorted sequence of these production rules. The variables in the left hand side in the sorted production rules are renamed in the sorted order, generating a set of new production rules that is assigned to the corresponding nodes in the ESP-tree. The same scheme is applied to the next level of the ESP-tree, which iterates until it reaches the root node.

Figure 3 shows an example of the encoding scheme for the ESP-tree built from an input string $S=babababaaba$. At the first and second levels in the ESP-tree, the set of production rules, $\{X_1\to ab,X_2\to b{X}_1,X_3\to aa,X_4\to ba\}$, is sorted in the lexicographic order of the left symbols on right hand sides of production rules, which results in the sequence of production rules, $(X_1\to ab,X_3\to aa,X_2\to b{X}_1,X_4\to ba)$. The variables on the right hand side of the production rules are renamed in the sorted order, resulting in the new sequence $(X_1\to ab,X_2\to aa,X_3\to b{X}_1,X_4\to ba)$, whose production rules are assigned to the corresponding nodes in the ESP-tree. This scheme is repeated until it reaches level 4.

Using the above encoding scheme, we obtain a monotonically increasing sequence of left symbols on the right hand side of the production rules, i.e., $X_{l\left(i\right)}$ in the form of $X_i\to X_{l\left(i\right)}{X}_{r\left(i\right)}$. Let $A_l$ be the increasing sequence; $A_l$ can be efficiently encoded into a bit string by using the gap-encoding and the unary coding. For example, the gap-encoding represents the sequence $(1,1,3,5,8)$ by $(1,0,2,2,3)$, and it is further transformed to the bit string $0^1{10}^0{10}^2{10}^2{10}^{3}1=0110010010001$ by unary coding. Generally, for a sequence $(x_1,x_2,\dots ,x_n)$, its unary code U represents $x_i$ by ${\mathit{rank}}_0(U,{\mathit{select}}_1(U,i))$. Because the number of 0s and the number of 1s is $n+\sigma$ and n, respectively, the size of U is $2n+\sigma$ bits. The bit string is indexed by the rank/select dictionary.

Let $A_r$ be the sequence consisting of the right symbols on the right hand side of the production rules, i.e., $X_r\left(i\right)$ in the form of $X_i\to X_{l\left(i\right)}{X}_{r\left(i\right)}$. $A_r$ is represented using $(n+\sigma )lg(n+\sigma )$ bits. $A_r$ is indexed by GMR [21].

The space for storing $A_l$ and $A_r$ is $(n+\sigma )lg(n+\sigma )+2n+\sigma +o\left(\right(n+\sigma )lg(n+\sigma ))$ bits in total. $A_l$ and $A_r$ enable us to simulate fast queries on encoded ESP-trees, which is presented in the next subsection.

5.2. Query Processing on Tree

The encoded ESP-trees support four tree operations, $\mathit{LeftChild}$, $\mathit{RightChild}$, $\mathit{LeftParents}$ and $\mathit{RightParents}$, which are used in our search algorithm. $\mathit{LeftChild}\left(X_k\right)$ returns the left child $X_{l\left(k\right)}$ of $X_k$ and can be implemented on bit string $A_l$ in $O\left(1\right)$ time as $m={\mathit{select}}_1(A_l,X_k)$ and $\mathit{LeftChild}\left(X_k\right)=m-X_k$. $\mathit{RightChild}\left(X_k\right)$ returns the right child $X_{r\left(k\right)}$ of $X_k$ and can be implemented on array $A_r$ in $O(lglg(n+\sigma ))$ time as $X_j=\mathit{access}(A_r,X_k)$.

$\mathit{LeftParents}\left(X_k\right)$ and $\mathit{RightParents}\left(X_k\right)$ return sets of parents of $X_k$ as left and right children, respectively, i.e., $\mathit{LeftParents}\left(X_k\right)=\{X_i\in V:X_i\to X_k{X}_j{,}^{\forall }{X}_j\in (\Sigma \cup V)\}$ and $\mathit{RightParents}\left(X_k\right)=\{X_i\in V:X_i\to X_j{X}_k{,}^{\forall }{X}_j\in (\Sigma \cup V)\}$.

Because $A_l$ is a monotonic sequence, any $X_k$ appears consecutively in $A_l$. Using the unary encoding of $A_l$, $\mathit{LeftParents}\left(X_k\right)$ is computed by $\{p+i:p={\mathit{select}}_1(A_l,X_k),{\mathit{rank}}_0(A_l,p+i)={\mathit{rank}}_0(A_l,p)\}$ in $O\left(\right|\mathit{LeftParents}\left(X_k\right)\left|\right)$ time. $\mathit{RightParents}\left(X_k\right)$ can be computed by repeatedly applying select operations for $X_k$ on $A_r$ until no more $X_k$ appear, i.e., ${\mathit{select}}_{X_k}(A_r,p)$ for $1\le p\le n$. Thus, $\mathit{RightParents}\left(X_k\right)$ for $X_k\in V$ can be computed in $O\left(\right|\mathit{RightParents}\left(X_k\right)\left|\right)$ time.

5.3. Other Data Structures

As a supplemental data structure, siEDM computes the node characteristic vector, denoted by $F\left(X_i\right)$, for each variable $X_i$: the characteristic vector consisting of the frequency of any variable derived from $X_i$. The space for storing all node characteristic vectors of n variables is at most $n^{2}lg\left|S\right|$ bits. Figure 3-(V) shows an example of the node characteristic vectors for ESP-tree in Figure 3-(III). In addition, let $V\left(X_i\right)$ be a set of $X_i$ and variables appearing in all the descendant nodes under $X_i$, i.e., $V\left(X_i\right)=\{e\in (V\cup \Sigma ):F\left(X_i\right)\left(e\right)\ne 0\}$. Practically, $F\left(X_i\right)$ is represented by a sequence of a pair of $X_j\in V\left(X_i\right)$ and $F\left(X_i\right)\left(X_j\right)$. Additionally, because $F\left(X_i\right)=F\left(\mathit{LeftChild}\left(X_i\right)\right)+F\left(\mathit{RightChild}\left(X_i\right)\right)+(X_i,1)$ ($+(X_i,1)$ represents adding 1 to dimension $X_i$), the characteristic vectors can be stored per level 2 of the ESP-tree. The data structure is represented by a bit array $FB$ indexed by a rank/select dictionary and the characteristic vectors reduced per level 2 of ESP-tree. $FB$ is set to 1 for i-th bit if $F\left(X_i\right)$ is stored, otherwise it is 0. Then, $F\left(X_i\right)$ can be computed by ${\mathit{rank}}_1(FB,i)$-th characteristic vector if the i-th bit of $FB$ is 1; otherwise, $F\left(\mathit{LeftChild}\left(X_i\right)\right)+F\left(\mathit{RightChild}\left(X_i\right)\right)+(X_i,1)$.

Another data structure that siEDM uses is a non-negative integer vector named length vector, each dimension of which is the length of the substring derived from the corresponding variable (See Figure 3-(VI)). The space for storing length vectors of n variables is $nlg\left|S\right|$ bits.

From the above argument, the space of the siEDM’s index structure for n variables is $n(n+1)lg\left|S\right|+(n+\sigma )lg(n+\sigma )+2n+\sigma +o\left(\right(n+\sigma )lg(n+\sigma ))$ bits in total.

6. Search Algorithm

6.1. Baseline Algorithm

Given a $T\left(S\right)$, the maximal subtree decomposition of $S[i,j]$ is a sequence $(X_1,X_2,\dots ,X_m)$ of variable in $T\left(S\right)$ defined recursively as follows. $X_1$ is the variable of the root of the maximal subtree satisfying that $S\left[i\right]$ is its leftmost leaf and $\left|val\right(X_1\left)\right|\le j-i$. If $val\left(X_1\right)=S[i,j]$, then $\left(X_1\right)$ is the maximal subtree decomposition of $S[i,j]$. Otherwise, let $X_1,X_2,\dots ,X_m$ be already determined and $|val\left(X_1\right)val\left(X_2\right)\cdots val\left(X_m\right)|=k<j-i$. Then, let $X_{m+1}$ be the variable of the root of the maximal subtree satisfying that $S[i+k+1]$ is its leftmost leaf and $\left|val\right(X_{m+1}\left)\right|\le j-i-k$. Repeating this process until $val\left(X_1\right)val\left(X_2\right)\cdots val\left(X_m\right)=S[i,j]$, the maximal subtree decomposition is determined.

Based on the maximal subtree decomposition, we explain the outline of the baseline algorithm, called online ESP [17], for computing an approximation of EDM between two strings. $T\left(S\right)$ is constructed beforehand. Given a pattern Q, the online ESP computes $T\left(Q\right)$, and for each substring $S[i,j]$ of length $\left|Q\right|$, it computes the approximate EDM as follows. It computes the maximal subtree decomposition $(X_1,X_2,\dots ,X_m)$ of $S[i,j]$. Then, the distance ${\parallel F\left(Q\right)-F\left(S[i,j]\right)\parallel }_1$ is approximated by $\parallel F\left(Q\right)-{\sum }_{k=1}^{m}F\left(X_k\right){\parallel }_1$ because ESP-tree is balanced and then $\parallel F\left(S[i,j]\right)-{\sum }_{k=1}^{m}F\left(X_k\right){\parallel }_1=O(lgm)$. This baseline algorithm is, however, required to compute the characteristic vector of $S[i,j]$ at each position i. Next, we improve the time and space of the online ESP by finding those $\left|Q\right|$-grams for each variable X in $V\left(S\right)$ instead of each position i.

6.2. Improvement

The siEDM approximately solves Problem 1 with the same guarantees presented in Theorem 1. Let $X_i\in V\left(S\right)$ such that $\left|\mathit{val}\right(X_i\left)\right|>\left|Q\right|$. There are $\left|Q\right|$-grams formed by the string $\mathit{suf}(\mathit{val}\left(X_{l\left(i\right)}\right),\left|Q\right|-k)\mathit{pre}(\mathit{val}\left(X_{r\left(i\right)}\right),k)$ with $k=1,2,\dots ,\left(\right|Q|-1)$. Then, the variable $X_i$ is said to stab the $\left|Q\right|$-grams. The set of the $\left|Q\right|$-grams stabbed by $X_i$ is denoted by $\mathit{itv}\left(X_i\right)$. Let $\mathit{itv}\left(S\right)$ be the set of $\mathit{itv}\left(X_i\right)$ for all $X_i$ appearing in $T\left(S\right)$. An important fact is that $\mathit{itv}\left(S\right)$ includes any $\left|Q\right|$-gram in S. Using this characteristic, we can reduce the search space .

If a $\left|Q\right|$-gram R is in $\mathit{itv}\left(X_i\right)$, there exists a maximal subtree decomposition $X_{i_1},X_{i_2},\dots ,X_{i_m}$. Then, the $L_1$-distance of $F\left(Q\right)$ and ${\sum }_{j=1}^{m}F\left(X_{i_j}\right)$ guarantees the same upper bounds in the original ESP as follows.

Theorem 2. 

Let $R\in \mathit{itv}\left(X_i\right)$ be a $\left|Q\right|$-gram on S and $X_{i_1},X_{i_2},\dots ,X_{i_m}$ be its maximal subtree decomposition in the tree $T\left(X_i\right)$. Then, it holds that

$$\parallel F\left(Q\right)-\sum _{j=1}^{m}F\left(X_{i_j}\right){\parallel }_1=O(lg\left|Q\right|{lg}^{*}\left|S\right|)d(Q,R)$$

Proof. 

By Theorem 1, ${\parallel F\left(Q\right)-F\left(R\right)\parallel }_1=O(lg\left|Q\right|{lg}^{*}\left|S\right|)d(Q,R)$. On the other hand, for an occurrence of R in S, let $T\left(X_i\right)$ be the smallest subtree in $T\left(S\right)$ containing the occurrence of R, i.e., $R\in \mathit{itv}\left(X_i\right)$. For $T\left(R\right)$ and $T\left(X_i\right)$, let $s\left(R\right)$ and $s\left(X_i\right)$ be the sequences of the level 2 symbols in $T\left(R\right)$ and $T\left(X_i\right)$, respectively. By the definition of the ESP, it holds that $s\left(R\right)=\alpha \beta \gamma$ and $s\left(X_i\right)={\alpha }^{\prime }\beta {\gamma }^{\prime }$ for some strings satisfying $|\alpha {\alpha }^{\prime }\gamma {\gamma }^{\prime }|=O({lg}^{*}\left|S\right|\left)\right|$, and this is true for the remaining string β iteratively. Thus, $\parallel F\left(R\right)-F\left(X_i\right){\parallel }_1=O(lg|R|{lg}^{*}\left|S\right|)$ since the trees are balanced. Hence, by the equation

$$\begin{array}{cc}\hfill \parallel F\left(Q\right)-\sum _{j=1}^{m}F\left(X_{i_j}\right){\parallel }_1& =O(lg\left|Q\right|{lg}^{*}\left|S\right|)d(Q,R)+O(lg\left|Q\right|{lg}^{*}\left|S\right|)\hfill \\ & =O(lg\left|Q\right|{lg}^{*}\left|S\right|)d(Q,R)\hfill \end{array}$$

we obtain the approximation ratio. ☐

To further enhance the search efficiency, we present a lower bound of the $L_1$-distance between characteristic vectors, which can be used for reducing the search space.

Theorem 3 

(A lower bound μ). For any $X_i\in V\left(S\right)\cup V\left(Q\right)$, the inequality ${\parallel F\left(S\right)-F\left(Q\right)\parallel }_1\ge \mu \left(X_i\right)$ holds where

$$\mu \left(X_i\right)=\sum _{e\in V\left(S\right)⨁V\left(Q\right)}F\left(X_i\right)\left(e\right)$$

Proof. 

The $L_1$ distance between $F\left(S\right)$ and $F\left(Q\right)$ is divided into four classes of terms: (i) both members in $F\left(S\right)$ and $F\left(Q\right)$ are non-zero, (ii) both members in $F\left(S\right)$ and $F\left(Q\right)$ are zero, (iii) the members in $F\left(S\right)$ and $F\left(Q\right)$ are zero and non-zero, (iv) the members in $F\left(S\right)$ and $F\left(Q\right)$ are non-zero and zero, respectively. Terms consisting of class (iii) and (iv) can be written as ${\sum }_{e\in V\left(S\right)⨁V\left(Q\right)}F\left(S\right)\left(e\right)$, which is a lower bound of the $L_1$-distance. Thus, ${\parallel F\left(S\right)-F\left(Q\right)\parallel }_1\ge {\sum }_{e\in V\left(S\right)⨁V\left(Q\right)}F\left(S\right)\left(e\right)$. ☐

Theorem 4 

(Monotonicity of μ). If a variable $X_i$ derives $X_k$, the inequality $\mu \left(X_i\right)\ge \mu \left(X_k\right)$ holds.

Proof. 

Every entry in $F\left(X_k\right)$ is less than or equal to the corresponding entry in $F\left(X_i\right)$. Thus, the inequality holds. ☐

6.3. Candidate Finding

By Theorems 2, 3 and 4, the task of the algorithm is reduced to finding a maximal subtree decomposition $(X_{i_1},X_{i_2},\dots ,X_{i_m})$ within $X_i$. Given a threshold $\tau \ge 0$, for each $\left|Q\right|$-gram in $\mathit{itv}\left(S\right)$, the algorithm finds the candidate: the maximal subtree decomposition $(X_{i_1},X_{i_2},\dots ,X_{i_m})$ satisfying $\mu \left(X_{i_1}\right)+\mu \left(X_{i_2}\right)+\cdots +\mu \left(X_{i_m}\right)\le \tau$.

For an $X_i$ and an occurrence of some $\left|Q\right|$-gram in $\mathit{itv}\left(X_i\right)$, the $\left|Q\right|$-gram is formed by the expression $\mathit{suf}(\mathit{val}\left(X_{l\left(i\right)}\right),\left|Q\right|-k)\mathit{pre}(\mathit{val}\left(X_{r\left(i\right)}\right),k)$ for a k $(1\le k\le |Q|-1)$. The algorithm computes the maximal subtree decompositions $(x_1,x_2,\dots ,x_p)$ covering $\mathit{suf}(\mathit{val}\left(X_{l\left(i\right)}\right),|Q|-k)$ and $(y_1,y_2,\dots ,y_q)$ covering $\mathit{pre}(\mathit{val}\left(X_{r\left(i\right)}\right),k)$, and outputs $(x_1,\dots ,x_p,y_1,\dots ,y_q)$ covering the $\left|Q\right|$-gram when ${\sum }_{1\le i\le p}\mu \left(x_i\right)+{\sum }_{1\le i\le q}\mu \left(y_i\right)\le \tau$. We illustrate the computation of candidates satisfying $\mu \left(X_{i_1}\right)+\mu \left(X_{i_2}\right)+\cdots +\mu \left(X_{i_m}\right)\le \tau$ in Figure 4 and show the pseudo-code in Algorithm 1.

Applying all variables to Algorithm 1 enables us to find the candidates covering all solutions. There are no possibilities for missing any $\left|Q\right|$-grams in $\mathit{itv}\left(S\right)$ such that the $L_1$-distances between their characteristic vectors and $F\left(Q\right)$ are at most τ, i.e., false negatives. The set may include a false positive, i.e., the solution set encodes a $\left|Q\right|$-gram such that the $L_1$-distance between its characteristic vector and $F\left(Q\right)$ is more than τ. However, false positives are efficiently removed by computing the $L_1$-distance $\parallel F\left(Q\right)-{\sum }_{j=1}^{m}F\left(X_{i_j}\right){\parallel }_1$ as a post-processing.

Theorem 5. 

The computation time of FindCandidates is $O\left(n\right|Q|lglg(n+\sigma )(lg\left|S\right|+lg\left|Q\right|\left)\right)$.

Proof. 

Because the height of the ESP-tree is $O(lg\left|S\right|)$, for each variable X, the number of visited nodes is $O(lg\left|Q\right|+lg|S\left|\right)$. The computation time of $\mathit{LeftChild}\left(X\right)$ and $\mathit{RightChild}\left(X\right)$ is $O(lglg(n+\sigma ))$, and the time of FindLeft and FindRight is $O\left(\right|Q|lglg(n+\sigma )(lg\left|S\right|+lg\left|Q\right|\left)\right)$. Thus, for n iterations of the functions, the total computation time is $O\left(n\right|Q|lglg(n+\sigma )(lg\left|S\right|+lg\left|Q\right|\left)\right)$. ☐

Algorithm 1 to output the candidate $R\subseteq V\left(S\right)$ for $X\in V\left(S\right)$, a query pattern Q and a distance threshold τ.

1: function FindCandidates(X,Q,τ) 2:     for $j=1,2,\dots ,\left|Q\right|$ do 3:         $R\leftarrow \varphi$                              ▹ Initialize solution set  4:         $q_1\leftarrow$FindLeft($\mathit{LeftChild}\left(X\right),\left|Q\right|-j,0,R$)                ▹ for left child  5:         $q_2\leftarrow$FindRight($\mathit{RightChild}\left(X\right),j,0,R$)                 ▹ for right child  6:         if $q_1=0$ and $q_2=0$ then  7:            Output R  8: function FindLeft($X,q,d,R$)  9:     if $d>\tau$ then  10:         return ∞  11:     else if $q=0$ then  12:         return 0  13:     else if $X\in \Sigma$ then  14:         $d\leftarrow d+1$ if $X\notin V\left(Q\right)$  15:         $R\leftarrow R\cup \left\{X\right\}$  16:         return $q-1$  17:     else if $\left|\mathit{val}\right(X\left)\right|\le q$ then  18:         $d\leftarrow d+\mu \left(X\right)$  19:         $R\leftarrow R\cup \left\{X\right\}$  20:         return $q-\left|\mathit{val}\right(X\left)\right|$  21:     $q^{\prime }\leftarrow$FindLeft($\mathit{RightChild}\left(X\right),q,d,R$)  22:     if $q^{\prime }\ne 0$ then  23:         return FindLeft$(\mathit{LeftChild}\left(X\right),q^{\prime },d,R)$  24: function FindRight($X,q,d,R$)  25:     if $d>\tau$ then  26:         return ∞  27:     else if $q=0$ then  28:         return 0  29:     else if $X\in \Sigma$ then  30:         $d\leftarrow d+1$ if $X\notin V\left(Q\right)$  31:         $R\leftarrow R\cup \left\{X\right\}$  32:         return $q-1$  33:     else if $\left|\mathit{val}\right(X\left)\right|\le q$ then  34:         $d\leftarrow d+\mu \left(X\right)$  35:         $R\leftarrow R\cup \left\{X\right\}$  36:         return $q-\left|\mathit{val}\right(X\left)\right|$  37:     $q^{\prime }\leftarrow$FindRight($\mathit{LeftChild}\left(X\right),q,d,R$)  38:     if $q^{\prime }\ne 0$ then  39:         return FindRight$(\mathit{RightChild}\left(X\right),q^{\prime },d,R)$

6.4. Computing Positions

The algorithm also computes all the positions of $\mathit{val}\left(X_i\right)$, denoted by $P\left(X_i\right)=\{p\in \{1,2,\dots ,\left|S\right|\}:S[p,p+|\mathit{val}\left(X_i\right)|-1]=\left(X_i\right)\}$. Starting from $X_i$, the algorithm goes up to the root in the ESP-tree built from S. p is initialized to 0 at $X_i$. If $X_k$ through the pass from $X_i$ to the root is the parent with the right child $X_{r\left(k\right)}$ on the pass, non-negative integer $\left(\right|\mathit{val}\left(X_k\right)|-|\mathit{val}\left(X_{r\left(k\right)}\right)\left|\right)$ is added to p. Otherwise, nothing is added to p. When the algorithm reaches the root, p represents a start position of $\mathit{val}\left(X_i\right)$ on S, i.e., $\mathit{val}\left(X_i\right)=S[p,p+|\mathit{val}\left(X_i\right)|-1]$. To compute the set $P\left(X_i\right)$, the algorithm starts from $X_i$ and goes up to the root for each parent in $\mathit{RightParents}\left(X_i\right)$ and $\mathit{LeftParents}\left(X_i\right)$, which return sets of parents for $X_i$. Algorithm 2 shows the pseudo-code.

Algorithm 2 to compute the set P of all occurrence of $\mathit{val}\left(X\right)$ on S for $X\in V\left(S\right)$.

1: function ComputePosition(X)  2:     $P\leftarrow \left\{\varnothing \right\}$                    ▹ Initialize solution set  3:     Recursion(X, 1)  4: function Recursion(X,p)  5:     if X is the root node then  6:         $P\leftarrow P\cup \left\{p\right\}$  7:         return 8:     for each $X_p\in \mathit{RightParents}\left(X\right)$do      ▹X is the right child of $X_p$  9:         Recursion($X_p$,$p+|\mathit{val}\left(X_p\right)|-|\mathit{val}\left(X\right)|$)  10:     for each $X_p\in \mathit{LeftParents}\left(X\right)$do       ▹X is the left child of $X_p$  11:         Recursion($X_p$,p)

Theorem 6. 

The computation time of $P\left(X\right)$ is $O(\mathit{occ}lg\left|S\right|)$, where $\mathit{occ}$ is the number of occurrences of X in $T\left(S\right)$.

Proof. 

Using the index structures of $\mathit{RightParents}\left(X\right)$ and $\mathit{LeftParents}\left(X\right)$, we can traverse the path from any node with label $\left(X\right)$ to the root of $T\left(S\right)$ counting the position. The length of the path is $O(lg|S\left|\right)$. ☐

Theorem 7. 

The search time is $O\left(n\right|Q|lglg(n+\sigma )(lg\left|S\right|+lg\left|Q\right|)+occlg\left|S\right|)$ using the data structure of size $n(n+1)lg\left|S\right|+(n+\sigma )lg(n+\sigma )+2n+\sigma +o\left(\right(n+\sigma )lg(n+\sigma ))$ bits.

Proof. 

The time for computing $T\left(Q\right)$ and $F\left(Q\right)$ is $t_1=O\left(\right|Q|{lg}^{*}\left|S\right|)$. The time for finding candidates and computing $\parallel F\left(Q\right)-{\sum }_{j=1}^m{X}_{i_j}{\parallel }_1$ is $t_2=O\left(n\right|Q|lglg(n+\sigma )(lg\left|S\right|+lg\left|Q\right|))$ by Theorem 5. The time for computing positions is $O(occlg\left|S\right|)$ by Theorem 6. Thus, the total time for a query search is $t_1+t_2+t_3=O\left(n\right|Q|lglg(n+\sigma )(lg\left|S\right|+lg\left|Q\right|)+occlg\left|S\right|)$. The size of the data structure is derived by the results in Section 5. ☐

In Theorem 7, n and $occ$ are incomparable because $occ>n$ is possible for a highly repetitive string.

7. Experiments

We evaluated the performance of siEDM on one core of a quad-core Intel Xeon Processor E5540 (2.53GHz) machine with 144GB memory. We implemented siEDM using the rank/select dictionary and GMR in libcds (https://github.com/fclaude/libcds). We used two standard benchmark datasets of einstein and cere from repetitive text collections in the pizza and chili corpus (http://pizzachili.dcc.uchile.cl/repcorpus.html), which is detailed in

Table 1. Summary of datasets.

Dataset	Length	$\mathbf{|}\mathbf{\Sigma }\mathbf{|}$	Size (MB)
einstein	$467,626,544$	139	446
cere	$461,286,644$	5	440

. As a comparison method, we used the online pattern matching for EDM called online ESP (baseline) [17] that approximates EDM between a query Q and substrings of the length of $\left|Q\right|$ of each position of an input text. We randomly selected $S[i,j]$ as the query pattern Q for each $\left|Q\right|=50,100,500,1000$ and examined the performance.

Table 2. Comparison of the memory consumption for the query search.

Dataset	Einstein	Cere
siEDM (MB)	$17.12$	$254.75$
baseline (MB)	$6.98$	$10.95$

shows the memory consumption in the search of the siEDM and baseline. The memory consumption of siEDM was larger than the baseline for both texts because the baseline does not have characteristic vectors of each node and length vector.

Table 3. Comparison of the index size and construction time.

Dataset		Einstein	Cere
	Encoded ESP-tree (MB)	$1.18$	$19.92$
Index Size	Characteristic vector F (MB)	$15.35$	$227.34$
	Length vector L (MB)	$0.59$	$7.49$
Construction time (sec)		$117.65$	$472.21$

shows the size for each component of the index structure and the time for building the index structure on einstein and cere datasets. Most of the size of the index structure was consumed by the characteristic vector F. The index size of cere was much larger than that of einstein. The index sizes of cere and einstein were approximately 16 megabytes and 256 megabytes, respectively, because the number of variables generated from cere was much larger than that generated from einstein. The number of variables generated from einstein was $305,098$ and the number of variables generated from cere was $4,512,406$. The construction times of the index structures were 118 s for einstein and 472 s for cere. The results for constructing the index structures demonstrate the applicability of siEDM to moderately large, repetitive texts.

Figure 5 shows the total search time (sec.) of siEDM and the baseline for einstein and cere in distance thresholds τ from 10 to 60. In addition, this result does not contain the case $\tau <10$ because siEDM found no candidate under the condition. The query length is one of $\{50,100,500,1000\}$. Because the search time of baseline is linear in $\left|S\right|+\left|Q\right|$, we show only the fastest case: $q=\left|Q\right|=50$. The search time of siEDM was faster than baseline in most cases.

Figure 6 shows the detailed search time in second. CF is the time for finding candidates of Q in $T\left(S\right)$, DIST is the time for computing approximated $L_1$ distance by characteristic vectors, and PC is the time for determining the positions of all $\left|Q\right|$-grams within the threshold τ.

Figure 7 shows the number of nodes $T\left(S\right)$ visited by the algorithm, #TN, the number of candidate $\left|Q\right|$-grams computed by $FindCandidates$, #CAND, the number of true positives among candidate $\left|Q\right|$-grams, #TP, and the number of occurrences, #OCC. The most time-consuming task is the candidate finding.

By the monotonicity of characteristic vectors, pruning the search space for small distance thresholds and long query length is more efficient. Thus, it is expected that siEDM is faster for smaller distance thresholds and longer query lengths and the experimental results support this. The search time on cere is much slower than that on einstein because the number of generated production rules from cere is much larger than that from einstein, and a large number of iterations of FindCandidates is executed. In addition, the comparison of #CAND and #TP validates the efficiency of siEDM for candidate finding with the proposed pruning method.

In Figure 7, the algorithm failed to find a candidate. Such a phenomenon often appears when the required threshold τ is too small, because the ESP-tree $T\left(Q\right)$ is not necessarily identical to $T\left(S\right[i,j\left]\right)$ even if $Q=S[i,j]$. Generally, the parsing of $T\left(S\right[i,j\left]\right)$ is affected by a suffix of $S[1,i-1]$ and a prefix of $S[j+1,|S\left|\right]$ of length at most ${lg}^{*}\left|S\right|$.

As shown in Table 3 and Figure 5, the search time of siEDM depends on the size of encoded ESP-tree for the input. Finally, we confirm this feature by an additional experiment for other repetitive texts.

Table 4. Summary of additional datasets.

Dataset	Length	$\mathbf{|}\mathbf{\Sigma }\mathbf{|}$	Size (MB)
influenza	154808555	15	$147.64$
Escherichia_Coli	112689515	15	$107.47$

,

Table 5. Comparison of the memory consumption for the query search.

Dataset	Influenza	Escherichia_Coli
siEDM (MB)	$164.87$	$262.01$
baseline (MB)	$53.01$	$100.81$

and

Table 6. Comparison of the index size and construction time for additional datasets.

Dataset		Influenza	Escherichia_Coli
	Encoded ESP-tree (MB)	$9.92$	$20.21$
Index Size	Characteristic vector F (MB)	$150.87$	$234.91$
	Length vector L (MB)	$4.08$	$6.88$
Construction time (sec)		$290.33$	$420.43$

are the description of several datasets from the pizza & chili corpus. Figure 8 shows the search time of siEDM and baseline. This result supports our claim that siEDM is suitable for computing EDM of repetitive texts.

8. Conclusions

We have proposed siEDM, an efficient string index for computing approximate searching based on EDM. Experimental results demonstrated the applicability of siEDM to real-world repetitive text collections as well as a longer pattern search. Future work will make the search algorithm in siEDM faster, which would be beneficial for users performing query searches for EDM.