Algorithms for Jumbled Pattern Matching in Strings

Jumbled pattern matching is a special case of approximate pattern matching. It has been used as a filtering step in approximate pattern matching algorithms [15], but rarely considered in its own right.

The authors of [7] present an algorithm for finding all occurrences of a Parikh vector in a runlength encoded text. The algorithm’s time complexity is $O(n^{\prime}+\sigma)$, where $n^{\prime}$ is the length of the runlength encoding of $s$. Obviously, if the string is not runlength encoded, a preprocessing phase of time $O(n)$ has to be added. However, this may still be feasible if many queries are expected. To the best of our knowledge, this is the only algorithm that has been presented for the problem we treat here.

An efficient algorithm for computing all Parikh fingerprints of substrings of a given string was developed in [1]. Parikh fingerprints are Boolean vectors where the $k$’th entry is $1$ if and only if $a_{k}$ appears in the string. The algorithm involves storing a data point for each Parikh fingerprint, of which there are at most $O(n\sigma)$ many. This approach was adapted in [10] for Parikh vectors and applied to identifying all repeated Parikh vectors within a given length range; using it to search for queries of arbitrary length would imply using $\Omega(P(s))$ space, where $P(s)$ denotes the number of different Parikh vectors of substrings of $s$. This is not desirable, since, for arbitrary alphabets, there are non-trivial strings of any length with quadratic $P(s)$ [8].

In this paper, we present two novel algorithms which perform significantly better than the simple window algorithm, in the case where many queries arrive.

For the binary case, we present an algorithm which answers decision queries in $O(1)$ time, using a data structure of size $O(n)$ (Interval Algorithm, Sect. 3). The data structure is constructed in $\Theta(n^{2})$ time.

For general alphabets, we present an algorithm with expected sublinear runtime which uses $O(n)$ space to answer occurrence queries (Jumping Algorithm, Sect. 4). We present two different variants of the algorithm. The first one uses a very simple data structure (an inverted table) and answers queries in time $O(\sigma J\log(\frac{n}{J}+m)),$ where $J$ denotes the number of iterations of the main loop of the algorithm. We then show that the expected value of $J$ for the case of random strings and patterns is $O(\frac{n}{\sqrt{m}\sqrt{\sigma\log\sigma}})$, yielding an expected runtime of $O(n(\frac{\sigma}{\log\sigma})^{1/2}\frac{\log m}{\sqrt{m}})$, per query

The second variant of the algorithm uses wavelet trees [13] and has query time $O(\sigma J)$, yielding an overall expected runtime of $O(n(\frac{\sigma}{\log\sigma})^{1/2}\frac{1}{\sqrt{m}})$, per query. (Here and in the following, $\log$ stands for logarithm base $2$.)

Our simulations on random strings and real biological strings confirm the sublinear behavior of the algorithms in practice. This is a significant improvement over the simple window algorithm w.r.t. expected runtime, both for a single query and repeated queries over one string.

The Jumping Algorithm is reminiscent of the Boyer-Moore-like approaches to the classical exact string matching problem [6, 2, 14]. This analogy is used both in its presentation and in the analysis of the number of iterations performed by the algorithm.

Let $s=s_{1}\ldots s_{n}\in\Sigma^{*}$ be given. Recall that $prv(i)$ denotes the Parikh vector of the prefix of $s$ of length $i$, for $i=0,\ldots,n$, where $prv(0)=p(\epsilon)=(0,\ldots,0)$. Consider Parikh vector $p\in{\mathbb{N}}^{\sigma}$, $p\neq(0,\ldots,0)$. We make the following simple observations:

The algorithm moves two pointers $L$ and $R$ along the text, pointing at these potential positions $i$ and $j$. Instead of moving linearly, however, the pointers are updated in jumps, alternating between updates of $R$ and $L$, in such a manner that many positions are skipped. Moreover, because of the way we update the pointers, after any update it suffices to check whether $R-L=|q|$ to confirm that an occurrence has been found (cf. Lemma 2 below).

We first need to define a function firstfit, which returns the smallest potential position where an occurence of a Parikh vector can end. Let $p\in{\mathbb{N}}^{\sigma}$, then

Updating $R$: Assume that the left pointer is pointing at position $L$, i.e. no unreported occurrence starts before $L+1$. Notice that, if there is an occurrence of $q$ ending at any position $j>L$, it must hold that $prv(L)+q\leq prv(j)$. In other words, we must fit both $prv(L)$ and $q$ at position $j$, so we update $R$ to

Updating $L$: Assume that $R$ has just been updated. Thus, $prv(R)-prv(L)\geq q$ by definition of firstfit. If equality holds, then we have found an occurrence of $q$ in position $(L+1,R)$, and $L$ can be incremented by $1$. Otherwise $prv(R)-prv(L)>q$, which implies that, interspersed between the characters that belong to $q$, there are some “superfluous” characters. Now the first position where an occurrence of $q$ can start is at the beginning of a contiguous sequence of characters ending in $R$ which all belong to $q$. In other words, we need the beginning of the longest suffix of $s[L+1,R]$ with Parikh vector $\leq q$, i.e. the smallest position $i$ such that $prv(R)-prv(i)\leq q$, or, equivalently, $prv(i)\geq prv(R)-q$. Thus we update $L$ to

Finally, in order to check whether we have found an occurrence of query $q$, after each update of $R$ or $L$, we check whether $R-L=|q|$. In Figure 2, we give the pseudocode of the algorithm.

It remains to see how to compute the firstfit and $prv$ functions. We first prove that the algorithm is correct. For this, we will need the following lemma.

We illustrate the proof in Fig. 3.

Storing all prefix vectors of $s$ would require $O(\sigma n)$ storage space, which may be too much. Instead, we construct an “inverted prefix vector table” $I$ containing the increment positions of the prefix vectors: for each character $a_{k}\in\Sigma$, and each value $j$ up to $p(s)_{k}$, the position in $s$ of the $j$’th occurrence of character $a_{k}$. Formally, $I[k][j]=\min\{i\mid prv(i)_{k}\geq j\}$ for $j\geq 1$, and $I[k][0]=0$. Then we have

We can also compute the prefix vectors $prv(i)$ from table $I$: For $k=1,\ldots,\sigma$,

The obvious way to find these values is to do binary search for $j$ in each row of $I$. However, this would take time $\Theta(\sigma\log n)$; a better way is to use information already acquired during the run of the algorithm. By Lemma 2, it always holds that $L\leq R$. Thus, for computing $prv(R)_{k}$, it suffices to search for $R$ between $prv(L)_{k}$ and $prv(L)_{k}+(R-L)$. This search takes time proportional to $\log(R-L)$. Moreover, after each update of $L$, we have $L\geq R-m$, so when computing $prv(L)_{k}$, we can restrict the search for $L$ to between $prv(R)_{k}-m$ and $prv(R)_{k}$, in time $O(\log m)$. For more details, see Section 4.4.

Table $I$ can be computed in one pass over $s$ (where we take the liberty of identifying character $a_{k}\in\Sigma$ with its index $k$). The variables $c_{k}$ count the number of occurrences of character $a_{k}$ seen so far, and are initialized to $0$.

• Algorithm Construct $I$
• 1.

  for $i=1$ to $n$

• 2.

  $c_{s_{i}}=c_{s_{i}}+1$;

• 3.

  $I[s_{i}][c_{s_{i}}]=i;$

Table $I$ requires $O(n)$ storage space (with constant 1). Moreover, the string $s$ can be discarded, so we have zero additional storage. (Access to $s_{i},1\leq i\leq n,$ is still possible, at cost $O(\sigma\log n)$.)

A wavelet tree on $s\in\Sigma^{*}$ allows rank, select, and access queries in time $O(\log\sigma)$. For $a_{k}\in\Sigma$, ${\textit{rank}}_{k}(s,i)=|\{j\mid s_{j}=a_{k},j\leq i\}|$, the number of occurrences of character $a_{k}$ up to and including position $i$, while ${\textit{select}}_{k}(s,i)=\min\{j\mid{\textit{rank}}_{k}(s,j)\geq i\}$, the position of the $i$’th occurrence of character $a_{k}$. When the string is clear, we just use ${\textit{rank}}_{k}(i)$ and ${\textit{select}}_{k}(i)$. Notice that

• •

  $prv(j)=({\textit{rank}}_{1}(j),\ldots,{\textit{rank}}_{\sigma}(j))$, and

• •

  for a Parikh vector $p=(p_{1},\ldots,p_{\sigma})$, $\textsc{firstfit}(p)=\max_{k=1,\ldots,\sigma}\{{\textit{select}}_{k}(p_{k})\}$.

So we can use a wavelet tree of string $s$ to implement those two functions. We give a brief recap of wavelet trees, and then explain how to implement the two functions above in $O(\sigma)$ time each.

A wavelet tree is a complete binary tree with $\sigma=|\Sigma|$ many leaves. To each inner node, a bitstring is associated which is defined recursively, starting from the root, in the following way. If $|\Sigma|=1$, then there is nothing to do (in this case, we have reached a leaf). Else split the alphabet into two roughly equal parts, $\Sigma_{\rm left}$ and $\Sigma_{\rm right}$. Now construct a bitstring of length $n$ from $s$ by replacing each occurrence of a character $a$ by $0$ if $a\in\Sigma_{\rm left}$, and by $1$ if $a\in\Sigma_{\rm right}$. Let $s_{\rm left}$ be the subsequence of $s$ consisting only of characters from $\Sigma_{\rm left}$, and $s_{\rm right}$ that consisting only of characters from $\Sigma_{\rm right}$. Now recurse on the left child with string $s_{\rm left}$ and alphabet $\Sigma_{\rm left}$, and on the right child with $s_{\rm right}$ and $\Sigma_{\rm right}$. An illustration is given in Fig. 4. At each inner node, in addition to the bitstring $B$, we have a data structure of size $o(|B|)$, which allows to perform rank and select queries on bit vectors in constant time ([20, 9, 21]).

Now, using the wavelet tree of $s$, any rank or select operation on $s$ takes time $O(\log\sigma)$, which would yield $O(\sigma\log\sigma)$ time for both $prv(j)$ and $\textsc{firstfit}(p)$. However, we can implement both in a way that they need only $O(\sigma)$ time: In order to compute ${\textit{rank}}_{k}(j)$, the wavelet tree, which has $\log\sigma$ levels, has to be descended from the root to leaf $k$. Since for $prv(j)$, we need all values ${\textit{rank}}_{1}(j),\ldots,{\textit{rank}}_{\sigma}(j)$ simultaneously, we traverse the complete tree in $O(\sigma)$ time.

For computing $\textsc{firstfit}(p)$, we need $\max_{k}\{{\textit{select}}_{k}(p_{k})\}$, which can be computed bottom-up in the following way. We define a value $x_{u}$ for each node $u$. If $u$ is a leaf, then $u$ corresponds to some character $a_{k}\in\Sigma$; set $x_{u}=p_{k}$. For an inner node $u$, let $B_{u}$ be the bitstring at $u$. We define $x_{u}$ by $x_{u}=\max\{{\textit{select}}_{0}(B_{u},x_{\rm left}),$ ${\textit{select}}_{1}(B_{u},x_{\rm right})\}$, where $x_{\rm left}$ and $x_{\rm right}$ are the values already computed for the left resp. right child of $u$. The desired value is equal to $x_{\rm root}$.

Let $\mathbb{A}_{1}(s,q)$ denote the running time of the Jumping Algorithm using inverted tables over a text $s$ and a Parikh vector $q$, and $\mathbb{A}_{2}(s,q)$ that of the Jumping Algorithm using a wavelet tree. Further, let $J=J(s,q)$ be the number of iterations performed in the while loop in line 2, i.e., the number of jumps performed by the algorithm on the input $q.$

The time spent in each iteration depends on how the functions firstfit and $prv$ are implemented (lines 3 and 7). In the wavelet tree implementation, as we saw before, both take time $O(\sigma)$, so the overall runtime of the algorithm is

For the inverted table implementation, it is easy to see that computing firstfit takes $O(\sigma)$ time. Now denote, for each $i=1,\dots,J,$ by $\hat{L}_{i},\,\hat{R}_{i}$ the value of $L$ and $R$ after the $i$’th execution of line 3 of the algorithm, respectively.¹¹1The $\hat{L}_{i}$ and $\hat{R}_{i}$ coincide with the $L_{k}$ and $R_{k}$ from the proof of Theorem 4.1 almost but not completely: When an occurrence is found after the update of $L$, then the corresponding pair $L_{k},R_{k}$ is skipped here. The reason is that now we are only considering those updates that carry a computational cost. The computation of $prv(\hat{L_{i}})$ in line 3 takes $O(\sigma\log m)$: For each $k=1,\dots,\sigma,$ the component $prv(\hat{L_{i}})_{k}$ can be determined by binary search over the list $I[k][prv(\hat{R}_{i-1})_{k}-m],I[k][prv(\hat{R}_{i-1})_{k}-m+1],\dots,I[k][prv(\hat{R}_{i-1})_{k}].$ By $\hat{L_{i}}\geq\hat{R}_{i-1}-m,$ the claim follows.

The computation of $prv(\hat{R_{i}})$ in line 7 takes $O(\sigma\log(\hat{R_{i}}-\hat{R}_{i-1}+m)).$ Simply observe that in the prefix ending at position $\hat{R_{i}}$ there can be at most $\hat{R_{i}}-\hat{L_{i}}$ more occurrences of the $k$’th character than there are in the prefix ending at position $\hat{L_{i}}.$ Therefore, as before, we can determine $prv(\hat{R_{i}})_{k}$ by binary search over the list $I[k][prv(\hat{L}_{i})_{k}],I[k][prv(\hat{L}_{i})_{k}+1],\dots,I[k][prv(\hat{L}_{i})_{k}+\hat{R_{i}}-\hat{L_{i}}].$ Using the fact that $\hat{L_{i}}\geq\hat{R}_{i-1}-m,$ the desired bound follows.