FM-Adaptive: A Practical Data-Aware FM-Index

The FM-index [10, 11, 13, 12, 15, 24, 18, 29, 31] and the compressed suffix array ($\mathrm{𝙲𝚂𝙰}$) [22, 21, 23, 52, 53, 15, 42, 32, 19, 17] are space-efficient text indexes whose query times are proportional to the query pattern size plus the product of the output size and a small polylog function of $n$. The former maintains a succinct representation of the $\mathrm{𝙱𝚆𝚃}$ and the latter maintains a succinct representation of the neighbor function $\mathrm{\Psi }$ [22, 23]. Both are self-indexes in that they represent the original text, and thus it can be discarded.

Grossi and Vitter [22, 23] and Sadakane [52] introduced the compressed suffix array and the neighbor function $\mathrm{\Psi }$. Ferragina and Manzini [10, 11] designed the original entropy-compressed FM-index based upon the Burrows-Wheeler transform ($\mathrm{𝙱𝚆𝚃}$) [3, 44] and the mapping function $\mathrm{𝐿𝐹}$. The mapping function $\mathrm{𝐿𝐹}$ and the neighbor function $\mathrm{\Psi }$ are inverses of one other.

Grossi, Gupta and Vitter [21] introduced the elegant data structure known as the wavelet tree, which has since become ubiquitous in text indexing. Using the wavelet tree, they were the first to establish a $\mathrm{𝙲𝚂𝙰}$ self-index that provably achieves the asymptotically optimal space bound with leading coefficient 1 (i.e., $\sim {\mathscr{H}}_k(𝒯)$ bits, the $k$th order entropy of $𝒯$, as defined subsequently in Definition 2). Their analysis of the wavelet tree also applied to the FM-index and was the first to show asymptotic optimality as well for the FM-index [21, 12, 13, 40, 41]. The original space bound derived for the BWT had a leading term of $5n{\mathscr{H}}_k(𝒯)$ [10], which was improved to the asymptotically optimal $n{\mathscr{H}}_k(𝒯)$ [21, 12] using wavelet trees.

Subsequently to the early theoretical breakthroughs [22, 10, 21, 11, 23], there have been numerous improvements. Simpler implementations for the $\mathrm{𝙲𝚂𝙰}$ and the FM-index achieve high-order compression without explicit partitioning into separate $k$-contexts and thus using a single wavelet tree for the entire text $𝒯$ [15, 40, 41, 14]. (Here “$k$-context” denotes a length-$k$ prefix in $𝒯$ of a symbol $𝒯[\mathrm{𝑆𝐴}[i]]$, for some $k\ge 0$, where $\mathrm{𝑆𝐴}$ is the suffix array of $𝒯$.) We refer the reader to the nice survey of Navarro and Mäkinen [49], Navarro’s book [48] and other references in the literature [40, 41, 18, 25, 20, 37, 28, 27, 33, 50, 46, 16, 47, 5, 38, 31, 17, 32, 1, 34, 2, 42, 31].

For example, Foschini et al. [15] encoded the single wavelet tree of the $\mathrm{𝙱𝚆𝚃}$ sequence using run-length encoding and Elias $\gamma$ code to encode the runs; if the $k$-contexts are hypothetically overlaid onto the $\mathrm{𝙱𝚆𝚃}$ sequence, the encoding of each run length adapts implicitly to the frequency statistics of the current $k$-context(s), thus achieving zero-order compression, and thus by Definition 2 of $k$th-order entropy, they encode the overall node in $k$th-order entropy space. We use that same idea in this paper. Mäkinnen and Navarro [40, 41] did a careful analysis to show the surprising result that applying RRR [51] to a single Wavelet Tree of the entire $\mathrm{𝙱𝚆𝚃}$ sequence without any partitioning achieves the $n{\mathscr{H}}_k(𝒯)$ bits leading space term, for a similar reason as in [15]; the sublocks formed by [51] implicitly encode each $k$th-order context in roughly zero-order entropy space. Interesting experiments on the practical performance of these RLE and RRR approaches and related methods appear in [24]. Gog et al. [18] used a fixed-block boosting technique and the RRR method [51] to implement the FM-index in asymptotically optimal entropy-compressed space with an extra cost of $o(n\log \sigma )$ bits. Mäkinen et al. [42] proposed a $\mathrm{𝙲𝚂𝙰}$-based index for highly repetitive data collections.

These references also mention many practical applications. An example is a special-purpose $\mathrm{𝙱𝚆𝚃}$-based compressed index for biological FASTQ data that exploits specific characteristics of next-generation sequence data [31].

Developing space-efficient entropy-compressed self-indexes that achieve fast query performance both in theory and practice has been a challenging problem. In this paper, we propose FM-Adaptive, a fast data-aware FM-index applicable for a wide range of text strings with different alphabet sizes. The key is a new representation of the wavelet tree of the $\mathrm{𝙱𝚆𝚃}$, with new tradeoffs between space occupancy and search time. We propose several auxiliary data structures to support fast access to $\mathrm{𝙱𝚆𝚃}$.

Using gap $\gamma$ code instead of Elias-Fano (EF) code [9, 7] used in [31] for DNA sequences, we deduce a new space bound for general texts while the space bound derived by using the EF code in [31] holds only when the entropy is $o(1)$. In addition, the compressed index in [31] is geared towards FASTQ data, which have a specific format, and its queries are not identical to the queries on the general-purpose compressed indexes in this paper. Our index has an additional efficiency in that it allows the gap $\gamma$ code and the run-length $\gamma$ code to share the same lookup table as in [28], though the tables for the gap $\gamma$ code and for the run-length $\gamma$ code which we introduce in Section 3.1 for encoding methods have different meanings, thus avoiding the need to store the lookup table $R_{\mathrm{𝐸𝐹}}$ [31] for Elias-Fano code.

Ferragina, Giancarlo, and Manzini [14] discussed the run-length code and the gap code for the implementations of the Wavelet Tree. For a bit string $ℬ$ of length $n$, they can represent $ℬ$ in $|\mathrm{𝑅𝑙𝑒}(ℬ)|\le \max (2a,b+2a/3)|ℬ|{\mathscr{H}}_0(ℬ)+3b+1=4n{\mathscr{H}}_0(ℬ)+4$ bits (Lemma 3.3) using the run-length $\gamma$ code and in $|\mathrm{𝐺𝑒}(ℬ)|\le \max (a,b)|ℬ|{\mathscr{H}}_0(ℬ)+a\log |ℬ|+b+2=2n{\mathscr{H}}_0(ℬ)+2\log n+3$ bits (Lemma 3.4) using the gap $\gamma$ code for $a=2$ and $b=1$. Instead, we can represent $ℬ$ in $2n{\mathscr{H}}_0(ℬ)+2\log n+2$ bits [24, 28] using the run-length $\gamma$ code and in $2n{\mathscr{H}}_0(ℬ)$ bits in Lemma 3 using the gap $\gamma$ code, which avoids any second order terms.

Ferragina et al. [14] addressed the problem of compressing the $\mathrm{𝙱𝚆𝚃}$ $L$, but not the indexing mechanics of how to provide fast random access to $L$ how to decode $L$. In this paper, we address compression as well as the problem of how to index the data. We designed several efficient auxiliary data structures and the lookup table for fast access and decoding to support efficient pattern matching in a text. We also present a new locating algorithm that substantially improves the locating time in practice.

For input text $𝒯$ of $n$ symbols, we summarize below the novel contributions we make in this paper:

1. (1)

   We present a new compressed representation for the wavelet tree of the $\mathrm{𝙱𝚆𝚃}$, called FM-Adaptive, that incorporates the gap $\gamma$-encoding. Using this, we can implement FM-Adaptive in $2n{\mathscr{H}}_k+o(n\log \sigma )$ bits, for $k\le \alpha {\log }_{\sigma }n-1$ and any fixed constant , where ${\mathscr{H}}_k$ denotes the $k$th-order empirical entropy of $𝒯$. We can construct FM-Adaptive in $𝒪(n\log \sigma )$ time using $𝒪(n\log n)$ bits of space. In addition, for a bit string $ℬ$ of length $n$, we can represent $ℬ$ in $2n{\mathscr{H}}_0(ℬ)$ bits using the gap $\gamma$ code, which avoids and improves any second order terms in space bounds in [14, 24], where ${\mathscr{H}}_0(ℬ)$ is zero-order empirical entropy of $ℬ$, defined in Definition 2.

2. (2)

   We present an improved algorithm that provides substantially faster performance in practice for locating patterns in a text, based upon memoization, which takes advantage of the overlapping-subproblems property.

3. (3)

   We design table $R_{\gamma }$ (see Section 3.4) for accelerating the decoding of the gap $\gamma$-encoded blocks to support fast pattern matching in a text.

4. (4)

   Given any pattern $𝒫$ of $m$ symbols encoded in $m\log \sigma$ bits, using FM-Adaptive, we can count the number of occurrences of $𝒫$ in $𝒯$ in $𝒪(m\log \sigma )$ time, we can locate all $\mathrm{𝑜𝑐𝑐}$ occurrences in $𝒪(\mathrm{𝑜𝑐𝑐}{\log }^{1+ϵ}n)$ additional time, and we can retrieve a text substring $𝒯[\mathrm{𝑠𝑡𝑎𝑟𝑡},\mathrm{𝑠𝑡𝑎𝑟𝑡}+\mathrm{\ell }-1]$ in $𝒪(\mathrm{\ell }\log \sigma +{\log }^{1+ϵ}n)$ time.

5. (5)

   Extensive experiments demonstrate that FM-Adaptive generally provides faster query performance, often by a considerable amount, and comparable or better compression than other state-of-the-art FM-index methods. The source code is available online [30].

In Section 3 we introduce a new compression structure for wavelet trees. The search and locate functions are described and analyzed in Section 4. We report the experimental results in Section 5 and give final comments in Section 6.

Let $𝒯=𝒯[0,n-1]$ be a text string consisting of $n$ symbols from alphabet $\mathrm{\Sigma }$ of size $\sigma$, and let $\mathrm{𝑆𝐴}$ denote the suffix array of $𝒯$. The $\mathrm{𝙱𝚆𝚃}$ $L$ of $𝒯$ is defined as $L[i]=𝒯[\mathrm{𝑆𝐴}[i]-1modn]$. As shown in Table 1, the $\mathrm{𝙱𝚆𝚃}$ of $𝒯$ is the string $L$ formed by sorting the suffixes of $𝒯$ in lexicographical order, choosing the preceding symbol for each suffix, and concatenating those preceding symbols. (When the suffix is the entire string $𝒯$, we use the symbol $\mathrm{\#}$ as its preceding symbol.)

The Burrows-Wheeler transform is invertible. We can retrieve $𝒯$ from its $\mathrm{𝙱𝚆𝚃}$ $L$ by backward search [10] using the mapping function $\mathrm{𝐿𝐹}$: The value $\mathrm{𝐿𝐹}(i)$ is the lexicographical rank of the suffix with prefix $L[i]$, namely, $F[LF[i]]=L[i]$, where $F$ is the first string of the $\mathrm{𝙱𝚆𝚃}$.

We can compute $\mathrm{𝐿𝐹}(i)$ by $\mathrm{𝐿𝐹}(i)=𝒞(L[i])+\mathrm{𝑂𝑐𝑐}(L[i],i)$, where $\mathrm{𝑂𝑐𝑐}(L[i],i)$ is the number of occurrences of symbol $L[i]$ in $L[0,i]$, and $𝒞[L[i]]$ is the number of occurrences of symbols in $𝒯$ that are smaller than $L[i]$. Table 1 shows the $\mathrm{𝑆𝐴}$ and $\mathrm{𝙱𝚆𝚃}L$ of the string $𝒯=\text{tcaaaatatatgcaacatatagtattagattgtat\#}$ in which for each $i$ we show the first four symbols (starting with $F[i]$) and the last symbol (namely, $L[i]$) of conceptual suffixes of $𝒯$ in lexicographical order. The mapping function $\mathrm{𝐿𝐹}$ and the neighbor function $\mathrm{\Psi }$ are inverses of one another: $\mathrm{𝐿𝐹}(\mathrm{\Psi }(i))=\mathrm{\Psi }(\mathrm{𝐿𝐹}(i))$, as shown in Table 1.

Let $S$ denote a text string of length $n$ over an alphabet $\mathrm{\Sigma }=\{0,1,\mathrm{\dots },\sigma -1\}$. We define the wavelet tree ($\mathrm{𝚆𝚃}$) [21] of $S$ with $\sigma$ leaves labelled by the symbols of the alphabet $\mathrm{\Sigma }$ as follows: For each node $v$ in $\mathrm{𝚆𝚃}$, let ${\mathrm{\Sigma }}_v$ be the subset of symbols in the subtree rooted at $v$, and let $S_v$ be the substring of $S$ consisting of all the symbols of ${\mathrm{\Sigma }}_v$. For each internal node $v$ in $\mathrm{𝚆𝚃}$, let $B_v$ be a bit string with the same length as $S_v$ defined as $B_v[i]=0$ if $S_v[i]$ is in the left subtree of $v$; $B_v[i]=1$ if $S_v[i]$ is in the right subtree of $v$. Such a bit string representation happens recursively at each internal node; the substring at the internal node consists of the characters dispatched from the parent node, with their order in the text string preserved. The collective size of the bit strings at any given level of the tree is bounded by $n$. With proper encoding, the wavelet tree representation of $S$ achieves the zero-order entropy space bound [21].

The powerful wavelet tree data structure [21, 15, 24] reduces the problem of compressing a string over a finite alphabet to the problem of compressing a set of bit (i.e., binary) strings. It is an elegant and versatile data structure that allows efficient access, rank, and select queries on a text string $𝒯$ of $n$ symbols from a multisymbol alphabet:

• $\mathrm{■}$

  access$(𝒯,i)$: returns the symbol $𝒯[i]$.

• $\mathrm{■}$

  rank${}_c{}^{}(𝒯,i)$: returns the number of occurrences of symbol $c$ in $𝒯[0,i]$ for any $0\le i\le n-1$ and $c\in \mathrm{\Sigma }$.

• $\mathrm{■}$

  select${}_c{}^{}(𝒯,j)$: returns the position in $L$ of the $j$th occurrence of symbol $c$ for any $1\le j\le n$.

A key result of Grossi et al. [21] is that if each bit string ${ℬ}_i$ of the internal nodes of the wavelet tree is compressed to zero-order entropy, then the cumulative encoding is a zero-order entropy encoding of $𝒯$:

With the extra cost for auxiliary data structures, the total space occupied by a wavelet tree is $n{\mathscr{H}}_0(𝒯)+o(n\log \sigma )$ bits of space [21, 15, 24].

We use a single wavelet tree to represent the $\mathrm{𝙱𝚆𝚃}$ $L$, as first done in Foschini et al. [15]. Mäkinen and Navarro [39] used a single wavelet tree to represent the $S$ string of the run-length encoding of the FM-index (i.e., the run heads). Figure 1 shows the balanced wavelet tree for the example $L$ = tcacaattttcatttgtgaattaatagaaag#ataa from Table 1. There are various ways [15] to form the wavelet tree, such as using a Huffman criterion [26]. Each leaf of the wavelet tree corresponds to a distinct symbol. The column labeled $ℬ$ in Table 1 is the root bit string of the wavelet tree in Figure 1. A value of 0 (resp., 1) in $ℬ$ means that the corresponding symbol in $L$ lies in the left (resp., right) subtree. The bit string for the roots of each subtree are defined recursively.

Using the wavelet tree, we can turn each $\mathrm{𝑂𝑐𝑐}(L[i],i)$ computation into one ${\mathrm{𝚛𝚊𝚗𝚔}}_c()$ computation on the wavelet tree.

We let $ℬ$ denote the bit string of length $n$. Let $r$ be the number of occurrences of the least frequent bit in $ℬ$ such that $1\le r\le n/2$. Let denote the increasing sequence of positions of the least frequent bit in $ℬ$. We assume for convenience $p_0=0$. Let $g_j=p_j-p_{j-1}$ denote the gap sequence of pairwise differences of neighboring values of the position increasing sequence of the least frequent bit in $ℬ$ for $j=1,2,\mathrm{\dots },r$. We represent each gap using Elias $\gamma$ code [8]. The resulting gap $\gamma$-encoding of $ℬ$ is the bit string ${\gamma }_{\mathrm{𝑔𝑎𝑝}}(ℬ)=\gamma (g_1)\gamma (g_2)\mathrm{\cdots }\gamma (g_r)$. Obviously, ${\sum }_{j=1}^r{g}_i=p_r\le n$.

On the other hand, we can view $ℬ$ as a sequence of maximal runs of identical bits $ℬ=b_1^{{\mathrm{\ell }}_1},b_2^{{\mathrm{\ell }}_2}\mathrm{\cdots }{b}_s^{{\mathrm{\ell }}_s}$ for some $s\le n$, where $b_i$ and $b_{i+1}\in \{0,1\}$ and $b_i\ne b_{i+1}$ for . We can represent the length of each run using Elias $\gamma$ code [8]. The resulting run-length $\gamma$ encoding of $ℬ$ is the bit string ${\gamma }_{\mathrm{𝑟𝑙}}(ℬ)=\gamma ({\mathrm{\ell }}_1)\gamma ({\mathrm{\ell }}_2)\mathrm{\cdots }\gamma ({\mathrm{\ell }}_s)$.

The starting point of the structure of the FM-Adaptive algorithm follows from [31], which we incorporate into the following discussion. We start the summary by introducing the compression structure to represent bit strings of the nodes of the wavelet tree. We then give an example to show how the newly introduced gap $\gamma$ code works, which is the key to improve general-purpose performance.

To enable effective compression on the bit strings of the wavelet tree, we partition the bit strings of the wavelet tree nodes into blocks and categorize the blocks into three basic types:

1. 1.

   blocks consisting only of all 0s or of all 1s;

2. 2.

   blocks having relatively long runs of 0s or of 1s; and

3. 3.

   blocks having a random-like sequence of 0s and 1s.

For each block, we choose one of the following four compression methods to minimize its coding length: $\mathrm{All0}$/$\mathrm{All1}$, $\mathrm{GapG0}$/$\mathrm{GapG1}$, $\mathrm{RLG0}$/$\mathrm{RLG1}$, and $\mathrm{Plain}$.

Specifically, for the block of type 1 consisting only of all 0s or all 1s, we use the $\mathrm{All0}$/$\mathrm{All1}$ code to encode it; that is, no additional bits are needed to store it. For the block of type 2 having relatively long 0-runs or 1-runs, we use either the $\mathrm{GapG}$ code or the $\mathrm{RLG}$ code to encode it. Here $\mathrm{GapG}$ means that we apply the gap $\gamma$ code (see Section 3.1) to encode the block. $\mathrm{RLG}$ means that we apply the run-length $\gamma$ code (see Section 3.1) to encode the block. For the blocks of type 3 having a random-like sequence of 0s and 1s, we keep the block unchanged, denoted as $\mathrm{Plain}$. This partitioning is combined with the mixed encoding to represent the wavelet tree nodes. The compressed wavelet trees ($\mathrm{𝙲𝚆𝚃}$) are the ones whose bit strings of nodes are partition-based and mixed-encoded.

Let $ℬ$ denote a bit string of an internal node of a wavelet tree. We use following three steps to obtain a succinct representation of $ℬ$:

1. 1.

   Partition $ℬ$ into blocks ${ℬ}^j$ of size $b$, except possibly the last one.

2. 2.

   Combine $a/b$ contiguous blocks to form a superblock of size $a$.

3. 3.

   Apply the mixed encoding to encode each block. We build the encoded sequence $𝒮$ by concatenating the encoded blocks.

We also maintain five extra structures to support fast access to $ℬ$: $\mathrm{𝑆𝐵𝑟𝑎𝑛𝑘}$, $\mathrm{𝐵𝑟𝑎𝑛𝑘}$, $\mathrm{𝑆𝐵}$, $B$, and $M$, where $\mathrm{𝑆𝐵𝑟𝑎𝑛𝑘}$ stores the number of 1s in $ℬ$ preceding the current superblock; $\mathrm{𝐵𝑟𝑎𝑛𝑘}$ stores the number of 1s in $ℬ$ preceding the current block relative to the beginning of its enclosing superblock; $\mathrm{𝑆𝐵}$ stores the number of bits in $𝒮$ preceding the current superblock; $B$ stores the number of bits in $𝒮$ preceding the current block relative to the beginning of its enclosing superblock; and $M$ indicates the encoding method used in a block.

Table 2 shows the compression structure for the root bit string $ℬ=$ 100000111100111111001100101000100100 of the wavelet tree of the $\mathrm{𝙱𝚆𝚃}$ for the example string $𝒯$ from Table 1, where $n=36$, $b=6$, and $a=18$. The figure shows four encodings: $\mathrm{All1}$, $\mathrm{GapG}$, $\mathrm{RLG}$, and $\mathrm{Plain}$. For purposes of illustrating all four methods, the $\mathrm{RLG}$ code is used to encode block 4. But in the actual algorithm, $\mathrm{GapG}$ would instead be selected to encode block 4, since it would result in a shorter encoding.

Using encoding sequence $𝒮$ of $ℬ$ and auxiliary structures $\mathrm{𝑆𝐵𝑟𝑎𝑛𝑘}$, $\mathrm{𝐵𝑟𝑎𝑛𝑘}$, $\mathrm{𝑆𝐵}$, $B$, and $M$, we can compute ${\mathrm{𝚛𝚊𝚗𝚔}}_1(ℬ,i)$ by Equation (1). The $\mathrm{𝑙𝑟𝑎𝑛𝑘}$ operation performs $\mathrm{All0}$/$\mathrm{All1}$, $\mathrm{GapG0}$/$\mathrm{GapG1}$, $\mathrm{RLG0}$/$\mathrm{RLG1}$, or $\mathrm{Plain}$ decodings depending upon $M[\lfloor i/b\rfloor ]$ to return the number of 1s up to $\mathrm{𝑜𝑓𝑓𝑠𝑒𝑡}$ within block $\lfloor i/b\rfloor$ and $\mathrm{𝑜𝑓𝑓𝑠𝑒𝑡}=imodb$. The starting decoding position on $𝒮$ is $\mathrm{𝑆𝐵}[\lfloor i/a\rfloor ]+B[\lfloor i/b\rfloor ]$. We can access $ℬ[i]$ in a way similar to ${\mathrm{𝚛𝚊𝚗𝚔}}_1(ℬ,i)$.

In this section, we describe the construction of the compressed wavelet tree $\mathrm{𝙲𝚆𝚃}$. The $\mathrm{𝙲𝚆𝚃}$ consists of compressed bit strings of the wavelet tree nodes. We describe the compression in Algorithm 3 in Appendix A, consisting of the encoded sequence $𝒮$ and five auxiliary structures $\mathrm{𝑆𝐵𝑟𝑎𝑛𝑘}$, $\mathrm{𝐵𝑟𝑎𝑛𝑘}$, $\mathrm{𝑆𝐵}$, $B$, and $M$. It is simple to see we can construct the compressed wavelet tree in $𝒪(n\log \sigma )$ time.

In order to accelerate the ${\mathrm{𝑟𝑎𝑛𝑘}}_1(ℬ,i)$ computation for a gap $\gamma$-encoded block (GapG), we design table $R_{\gamma }$ for the gap $\gamma$-encoding for every possible chunk $w$ of $(\log n)/2$ bits. Using $R_{\gamma }$, we can process each $\mathrm{\Theta }(\log n)$ bits of a $\gamma$-encoded sequence in constant time. Every $\mathrm{𝑖𝑛𝑑𝑒𝑥}$ of $R_{\gamma }$, corresponding to a bit string of $w$ bits, contains three components: $R_1$, $R_2$, and $R_3$:

• $\mathrm{■}$

  $R_1$ stores the total number of decoded entries in the $\mathrm{𝑖𝑛𝑑𝑒𝑥}$. It is the number of decoded gaps.

• $\mathrm{■}$

  $R_2$ stores the cumulative sum of decoded digits (gap values) in the $\mathrm{𝑖𝑛𝑑𝑒𝑥}$, which corresponds to the maximum value of decoded positions.

• $\mathrm{■}$

  $R_3$ stores the total number of decoded bits in the $\mathrm{𝑖𝑛𝑑𝑒𝑥}$.

The function $\mathit{GapG1rank}()$ in Algorithm 4 in Appendix A implements the $\mathrm{𝑙𝑟𝑎𝑛𝑘}()$ operation using $R_{\gamma }$ for the GapG1-encoded block. The bit operator $\ll$ shifts the bits in the argument to the left by the designated number of positions.

We let ${ℬ}_i$ of length $n_i$ denote the bit string of internal node $i$ of the wavelet tree of the $\mathrm{𝙱𝚆𝚃}$ $L$ for $i=1,2,\mathrm{\dots },t=\sigma -1$, and let $n_i$ denote the number of bits contained in ${ℬ}_i$. We let $r$ be the number of occurrences of the least frequent bit in ${ℬ}_i$ such that $1\le r\le n_i/2$. Applying the gap $\gamma$ code to ${ℬ}_i$, we get the resulting bit string ${\gamma }_{\mathrm{𝑔𝑎𝑝}}({ℬ}_i)=\gamma (g_1)\gamma (g_2)\mathrm{\cdots }\gamma (g_r)$, whose space bound is given in Lemma 3; it represents an improved space bound over the ones in [14, 24]:

We sample suffix array $\mathrm{𝑆𝐴}$ in the same manner as in Huo et al. [32], and we denote by ${\mathrm{𝑆𝐴}}_s$ the $\mathrm{𝑆𝐴}$ sampling, where $s$ is the $\mathrm{𝑆𝐴}$ sampling step size. We use the conceptual bit array $D$ of length $n$ to record the indices corresponding to the sampled values in ${\mathrm{𝑆𝐴}}_s$. The locate algorithm without memoization is given in Appendix A.

In this section, we consider how to retrieve text substrings $𝒯[\mathrm{𝑠𝑡𝑎𝑟𝑡},\mathrm{𝑠𝑡𝑎𝑟𝑡}+\mathrm{\ell }-1]$ using the $\mathrm{𝙲𝚆𝚃}$ of $L$ and ${\mathrm{𝑆𝐴}}_s^{-1}$. We sample inverse suffix array ${\mathrm{𝑆𝐴}}^{-1}$ in the same manner as in Huo et al. [32]. We denote by ${\mathrm{𝑆𝐴}}_s^{-1}$ the ${\mathrm{𝑆𝐴}}^{-1}$ sampling.

The workflow to retrieve a substring using the $\mathrm{𝙲𝚆𝚃}$ of $L$ and ${\mathrm{𝑆𝐴}}_s^{-1}$ is that we first transform a given position into its rank $i$ in $\mathrm{𝑆𝐴}$ and then retrieve $𝒯[\mathrm{𝑠𝑡𝑎𝑟𝑡},\mathrm{𝑠𝑡𝑎𝑟𝑡}+\mathrm{\ell }-1]$ by walking on $\mathrm{𝐿𝐹}$. The first step is done by the $\mathrm{𝑡𝑟𝑎𝑛𝑠𝑓𝑜𝑟𝑚}()$ procedure and the second step by the $\mathrm{𝑟𝑒𝑡𝑟𝑖𝑒𝑣𝑒}()$ procedure. Both procedures are given in Algorithm 2.

In this section we describe experiments performed using the environment described in [32]. We used C++ to implement the algorithms and constructed the suffix array with parameter 64 using Mori’s fast lightweight suffix-sorting library.³³3github.com/y-256/libdivsufsort/

Table 3 summarizes some statistical characteristics of the data sets we used for the experiments. The data sets consist of four datasets from the ⁴⁴4pizzachili.dcc.uchile.cl/texts.html (datasets 1–4), the highly repetitive data sets from the $\mathrm{𝑅𝑒𝑝𝑒𝑡𝑖𝑡𝑖𝑣𝑒}\mathrm{𝐶𝑜𝑟𝑝𝑢𝑠}$⁵⁵5pizzachili.dcc.uchile.cl/repcorpus.html (datasets 5–8), the human genome (called hg38 at UCSC) from UCSC⁶⁶6hgdownload.cse.ucsc.edu/goldenPath/hg38 (dataset 9), and NA12877R10 (dataset 10), formed by extracting 10 gigabytes of reads from NA12877_1, downloaded from EBI.⁷⁷7www.ebi.ac.uk/ena/browser/view/ERR194146 Datasets 1–4 and 9–10 are nonrepetitive data sets, and datasets 5–8 are highly repetitive data sets. The expression $n/r$ denotes the average run length in the $\mathrm{𝙱𝚆𝚃}L$, and $\sigma$ is the alphabet size of the input data. For hg38, we exclude the query on pattern NN…N.

In our experiments, we examine the following state-of-the-art algorithms for their space usage and query time:

1. 1.

   AF-Index: the implementation⁸⁸8pizzachili.dcc.uchile.cl/indexes of alphabet-friendly FM-index [13, 12], which combines an existing compression boosting technique with the wavelet tree data structure. We used its latest version af-index_v2.1.

2. 2.

   FMI-Hybrid: the recent implementation⁹⁹9github.com/simongog/sdsl-lite of the FM-index [18], which uses the fixed-block boosting technique [36], in which bitvectors by default are implemented by hybrid encoding [37].

3. 3.

   RL-CSA: the implementation¹⁰¹⁰10jltsiren.kapsi.fi/rlcsa of the compressed suffix array [42] that uses run-length encoding and has been optimized for highly repetitive data. We used its current version of May 2016.

4. 4.

   GeCSA: the recent implementation¹¹¹¹11codeocean.com/capsule/3554560/tree/v1 of the compressed suffix array [32], in which a new run-length $\delta$ code is introduced for the gap sequence of $\mathrm{\Psi }$ [22, 23].

5. 5.

   FM-Adaptive: our method described in this paper.

In this section, we make an experimental comparison of the locate query before and after the improvement. The key point of the memoization improvement is to remember the computed suffix positions for some $j\in [l,r]$ and then use these computed positions and other auxiliary information to determine the suffix positions of the remaining $j$ so as to reduce the number of accesses to $\mathrm{𝐿𝐹}$ and thus speed up the locate query.

In the experiments, we perform 500 locate queries and calculate the average number of access to $\mathrm{𝐿𝐹}$ for each locate query. Figure 3 shows the average number of access to $\mathrm{𝐿𝐹}$ for each locate query before (in blue) and after (in orange) the improvement by memoization.

As can be seen in Figure 3, the memoization of some computed suffix positions generally improves the locate query time substantially. The average number of accesses to $\mathrm{𝐿𝐹}$ is decreased by $88.69\text{<em class="ltx\_emph ltx\_font\_italic">–</em>}99.11$ percent on english, sources, para, w.leaders, and kernel, by $44.88\text{<em class="ltx\_emph ltx\_font\_italic">–</em>}47.31$ percent on hg38, by $6.29\text{<em class="ltx\_emph ltx\_font\_italic">–</em>}17.82$ percent on NA12877R10, by $0.05\text{<em class="ltx\_emph ltx\_font\_italic">–</em>}2.2$ percent on proteins, by $0.01\text{<em class="ltx\_emph ltx\_font\_italic">–</em>}0.08$ percent on DNA, and by $0.0\text{<em class="ltx\_emph ltx\_font\_italic">–</em>}0.01$ percent on influenza for all $s$.

Figure 3 shows substantial reductions in $\mathrm{𝐿𝐹}$ invocations compared with the locate algorithm without memoization shown in Appendix A for most tested data. This memoization did not work on the datasets DNA, proteins, and influenza because there are few short patterns occurring frequently, which make $i={\mathrm{𝐿𝐹}}^k(j)\in [l,r]$ rarely satisfied.

To show that the improvement of the locate query is largely due to the memoization rather than the suffix array sampling method, we made additional experiments on three locate algorithms on the , as well as the $\mathrm{𝑅𝑒𝑝𝑒𝑡𝑖𝑡𝑖𝑣𝑒}\mathrm{𝐶𝑜𝑟𝑝𝑢𝑠}$, which we show in Table 4.

We denote the locate algorithm with position sampling in suffix array in [28] as $\mathit{locate15}$, the locate algorithm without memoization with suffix array value sampling (Algorithm 5 in Appendix A) as $\mathrm{𝑙𝑜𝑐𝑎𝑡𝑒}$, and the locate algorithm with memoization with suffix array value sampling (Algorithm 1) as $\mathrm{𝑓𝑎𝑠𝑡𝐿𝑜𝑐𝑎𝑡𝑒}$. We randomly select $1000$ patterns with length $20$ from the input text, and average the total time over the $1000$ locate queries on each of the three locate algorithms, denoted as $T_{\mathrm{𝑝𝑠}}$ for $\mathit{locate15}$, $T_{\mathrm{𝑣𝑠}}$ for $\mathrm{𝑙𝑜𝑐𝑎𝑡𝑒}$, and $T_{\mathrm{𝑓𝑎𝑠𝑡}}$ for $\mathrm{𝑓𝑎𝑠𝑡𝐿𝑜𝑐𝑎𝑡𝑒}$ in milliseconds, shown in Table 4. We take block size $b=256$ for the three locate algorithms. It can be seen that $\mathrm{𝑓𝑎𝑡𝑠𝐿𝑜𝑐𝑎𝑡𝑒}$ provides generally several times to orders of magnitude faster than $\mathit{locate15}$ and $\mathrm{𝑙𝑜𝑐𝑎𝑡𝑒}$ for different suffix array sampling step size $s=32,64$ and $128$.

In the following sections we compare the performance of our proposed algorithms with the state-of-the-art indexing methods described above, based upon the performance criteria of compression ratio, locate query time, and extract query time. We use the following parameters for each index:

• $\mathrm{■}$

  AF-Index: $b=32\ast 8$, and $s$ = 64, 128, and 256;

• $\mathrm{■}$

  FMI-Hybrid: $b=256$, and $s$ = 32, 64, and 128;

• $\mathrm{■}$

  RL-CSA: $b=64$, $s$ = 32, 64, and 128;

• $\mathrm{■}$

  GeCSA: $b=256$, and $s$ = 128, 256, and 512;

• $\mathrm{■}$

  FM-Adaptive: $b=256$, and $s$ = 64, 128, and 256.

The other parameters used are the default parameters. The tested method used is the same as those in [32]. The compression ratio is defined as the ratio of the index size with the original size of the input text. The parameters in our index are $a=16b$, $d=s$, and $w=16$.