A Simple Algorithm for Computing the Document Array

TL;DR

This paper introduces a straightforward, linear-time algorithm for computing the document array from a string collection and its suffix array, using minimal extra space for constant alphabet strings.

Contribution

It provides a simple, efficient method for document array computation with linear time complexity and constant space for constant alphabets.

Findings

Runs in linear time

Uses constant additional space for constant alphabets

Efficiently computes document arrays

Abstract

We present a simple algorithm for computing the document array given a string collection and its suffix array as input. Our algorithm runs in linear time using constant additional space for strings from constant alphabets.

Tables

Table 1. Table 1: Datasets. We used 32 32 32 -bits integers to store 𝖲𝖠 ​ [ 1 , N ] 𝖲𝖠 1 𝑁 \mathsf{SA}[1,N] when N < 2 31 𝑁 superscript 2 31 N<2^{31} (2GB), otherwise we used 64 64 64 -bits. The document array 𝖣𝖠 ​ [ 1 , N ] 𝖣𝖠 1 𝑁 \mathsf{DA}[1,N] is stored using 32 32 32 -bits integers, since d 𝑑 d is always smaller than 2 31 superscript 2 31 2^{31} . Each symbol of T c ​ a ​ t superscript 𝑇 𝑐 𝑎 𝑡 T^{cat} uses 1 byte.

Dataset	$\sigma$	$N/2^{30}$	$d$	$N/d$	longest string
revision	203	0.39	20,433	20,527	2,000,452
influenza	15	0.56	394,217	1,516	2,867
\hdashlinereads	4	2.87	32,621,862	94	101
pages	205	3.74	1,000	4,019,585	362,724,758
wikipedia	208	8.32	3,903,703	2,288	224,488
proteins	25	15.77	50,825,784	333	36,805

Table 2. Table 2: Running time and workspace.

Dataset	Time (seconds)			Workspace (KB)
	Alg. 1	bit_plain	bit_sd	Alg. 1	bit_plain	bit_sd
revision	60.88	11.74	20.37	1	64,002	44
influenza	109.13	20.48	41.24	1	91,168	704
\hdashlinereads	931.35	150.40	549.65	2	470,389	38,980
pages	762.91	141.99	141.25	2	613,341	4
wikipedia	2,947.59	450.64	1,054.08	2	1,363,147	7,096
protein	7,007.87	1,211.13	2,899.63	2	2,583,532	69,423

Equations

$$\mathsf{BWT}[i]=T[\mathsf{SA}[i]-1\mathrm{mod}n].$$

$$\mathsf{LF}(i,c)=\mathsf{C}[c]+{\mathsf{rank}}_c(\mathsf{BWT},i).$$

$$\mathsf{LF}[i]=\mathsf{ISA}[(\mathsf{SA}[i]-1)\mathrm{mod}n]$$

$$bit[i]=1\text{mbox}iff{T}^{cat}[i]=\$\text{mbox}andbit[i]=0\text{mbox},otherwise.$$

$$\mathsf{DA}[i]={\mathsf{rank}}_1(bit,\mathsf{SA}[i])+1,$$

Full text

††journal: Journal of LaTeX Templates

1 Introduction

The suffix array ($\mathsf{SA}$) [1] is a fundamental data structure in string processing that is commonly accompanied by the document array [2] when indexing string collections (e.g. [3, 4, 5, 6, 7]). Given a collection of $d$ strings of total length $N$, the document array is an array of integers $\mathsf{DA}[1,N]$ in the range $[1,d+1]$ that gives which document each suffix in the suffix array belongs to.

It is well-known that $\mathsf{DA}$ can be represented in a compact form by using a bitvector $bit[1,N]$ with support to rank operations, requiring $N+o(N)$ bits of space [8]. However, there are applications where $\mathsf{DA}$ must be accessed sequentially (e.g. [9, 10, 11, 12, 13, 14]), and having the array $\mathsf{DA}[1,N]$ computed explicitly is important.

In this paper we show how to compute $\mathsf{DA}$ given a string collection and its suffix array as input in $O(N)$ time.

Our algorithm reuses the space of $\mathsf{SA}$ to store auxiliary arrays used to compute $\mathsf{DA}$. $\mathsf{SA}$ is reconstructed in its original space during $\mathsf{DA}$ computation.

The workspace of our algorithm, that is, the extra space used in addition to the input and output, is $O(\sigma\lg N)$ bits, where $\sigma$ is the alphabet size. Therefore, for constant alphabets, the workspace of our algorithm is constant111We assume a computer word size of $\lg N$ bits..

2 Background

Let $T$ be a string of length $n$, over an alphabet $\Sigma$ of size $\sigma$, such that T[n]=\$$ is an end-marker symbol that does not occur elsewhere in T$and precedes every symbol of$\Sigma$.$T[i,j]$denotes the substring from$T[i]$to$T[j]$inclusive, for$1\leq i\leq j\leq n$. A suffix of$T$is a substring$T[i,n]$. We define$\mathsf{rank}_{c}(T,i)$as the number of occurrences of symbol$c$in$T[1,i]$. The string$T$is stored in$n\lg\sigma$ bits of space.

The suffix array ($\mathsf{SA}$) [1] for $T$ is an array of integers in the interval $[1,n]$ that provides the lexicographical order of all suffixes of $T$. The inverted permutation of $\mathsf{SA}$, denoted as $\mathsf{ISA}$, is defined as $\mathsf{ISA}[\mathsf{SA}[i]]=i$. $\mathsf{SA}$ can be computed in $O(n)$ time using $O(\sigma\lg n)$ bits of workspace [15]. Since, $\mathsf{SA}$ and $\mathsf{ISA}$ are permutations of $[1,n]$, the arrays $\mathsf{SA}$ and $\mathsf{ISA}$ use $n\lg n$ bits of space each one.

The Burrows-Wheeler transform ($\mathsf{BWT}$) [16] of $T$ is obtained by sorting all $n$ rotations of $T$ in a conceptual matrix $\mathcal{M}$, and taking the last column $L$ as the $\mathsf{BWT}$. It can also be defined through the relation

[TABLE]

The $\mathsf{BWT}$ can be computed directly (without $\mathsf{SA}$) in $O(n)$ time using $O(n\lg\sigma)$ bits of workspace [17, 18], or alternatively in $O(n^{2})$ time in-place [19].

The Last-to-First ($\mathsf{LF}$) mapping states that the $k^{th}$ occurrence of a symbol $c$ in column $L$ of $\mathcal{M}$ and the $k^{th}$ occurrence of $c$ in the first column $F$ correspond to the same symbol in $T$. Let $\mathsf{C}[c]$ be the number of symbols $c^{\prime}<c$ in $T[1,n]$. We define

[TABLE]

We use shorthand $\mathsf{LF}(i)$ for $\mathsf{LF}(i,\mathsf{BWT}[i])$. $\mathsf{LF}(i)$ may be computed on-the-fly in $O(\lg\sigma)$ time querying a wavelet tree [20] for the rank queries on Equation 2. The wavelet tree requires additional $(n~{}+~{}o(n))\lg\sigma$ bits of space.

$\mathsf{LF}$-mapping allows us to navigate $T$ right-to-left, given $T[k]=\mathsf{BWT}[i]$, then $T[k-1]=\mathsf{BWT}[\mathsf{LF}(i)]$. $T$ can be reconstructed backwards from $\mathsf{BWT}$ starting with T[n]=\mathsf{BWT}[1]=\$$, and repeatedly applying \mathsf{LF}$for$n$ steps.

2.1 String collections

Let $\mathcal{T}=T_{1},T_{2},\dots,T_{d}$ be a collection of $d$ strings of lengths $n_{1},n_{2},\dots,n_{d}$. The suffix array of $\mathcal{T}$ is the $\mathsf{SA}$ built for the concatenation of all strings $T^{cat}[1,N]=T_{1}T_{2}\dots T_{d}\#$, with total size $N=\sum_{i=1}^{d}(n_{1})+1$ and a new end-marker \#<\$$. \mathsf{SA}$can be computed in$O(N)$time using$O(\sigma\lg N)bits of workspace [[21](#bib.bib21)], such that an end-marker$$from string$T_{i}$is smaller than a$$$from string$T_{j}$iff$i<j$, which is equivalent to using$d$ different end-markers as separators, without increasing the alphabet size.

The $\mathsf{BWT}$ may also be generalized for string collections. $\mathsf{BWT}$ of $\mathcal{T}$ is obtained from $\mathsf{SA}$ and $T^{cat}$ through the Equation 1. However, $\mathsf{LF}$-mapping through Equation 2 does not work for symbols \$$, since the k^{th}$symbol$$$in column$L$does not (necessarily) corresponds to the$k^{th}$symbol$$$in column$F$, in this case$\mathsf{LF}(i,$)$ is undefined [7].

$\mathsf{LF}(i)$ can be pre-computed in an array $\mathsf{LF}[1,n]$ through Equation 3 such that $\mathsf{LF}$ still works for \$$-symbols. Given \mathsf{SA}$and$\mathsf{ISA}$, we have

[TABLE]

The array $\mathsf{LF}$ uses $N\lg N$ bits of space.

2.2 Document array

The document array ($\mathsf{DA}$) is an array of integers in the interval $[1,d+1]$ that tells us which document $j\in\mathcal{T}$ each suffix in the $\mathsf{SA}$ belongs to [2]. We define $\mathsf{DA}[i]=j$ iff suffix $T^{cat}[\mathsf{SA}[i],N]$ came from string $T_{j}\in\mathcal{T}$. $\mathsf{DA}[1]=d+1$ for the last suffix $T^{cat}[N,N]=\#$. $\mathsf{DA}[1,N]$ uses $N\lg(d+1)$ bits of space.

The array $\mathsf{DA}[1,N]$ can also be represented using a wavelet tree [20], within the same $N\lg(d+1)$ bits with additional functionalities [3]. $\mathsf{DA}$ can still be compressed using grammars when the string collection is repetitive [22].

2.3 Related work

Given $T^{cat}$ and $\mathsf{SA}$, the document array $\mathsf{DA}$ can be constructed in $O(N)$ time using $N\lg N$ additional bits to store $\mathsf{ISA}$, such that $\mathsf{DA}[\mathsf{ISA}[i]]=j$ for $i=\ell_{j-1},\dots,\ell_{j}$, with $\ell_{0}=1$ and $\ell_{j}=\sum_{k=1}^{j}n_{k}$, see [23, Alg. 5.29].

$\mathsf{DA}$ can also be computed in the same fashion as the text $T^{cat}$ is reconstructed from its $\mathsf{BWT}$. Given $T^{cat}$ and $\mathsf{SA}$, we can compute $\mathsf{ISA}[1,N]$ and then array $\mathsf{LF}[1,N]$ (Equation 3). $\mathsf{DA}$ is obtained in $O(N)$ time during the $\mathsf{BWT}$ inversion using $N\lg N$ bits of workspace, see [23, Alg. 7.30].

In particular, in Section 3 we show an alternative algorithm that reuses the space of $\mathsf{SA}$ to compute $\mathsf{LF}$ without $\mathsf{ISA}$. Our algorithm uses $O(\sigma\lg N)$ bits of workspace and reconstructs $\mathsf{SA}$ during $\mathsf{DA}$ computation.

Lightweight alternative

$\mathsf{DA}$ can be computed using a compact data structure composed by a bitvector $bit[1,N]$ with rank support operation. $bit$ is built over $T^{cat}[1,N]$, such that

[TABLE]

$\mathsf{DA}[i]$ can be obtained using $bit$ and $\mathsf{SA}$ as follows [23, Alg. 7.29]:

[TABLE]

$bit[1,N]$ can be pre-processed in $O(N)$ time so that rank queries are supported in $O(1)$ time using additional $o(N)$ bits [24]. This procedure computes $\mathsf{DA}$ in $O(N)$ time using $N+o(N)$ bits of workspace.

3 Computing $\mathsf{DA}$

In this section we show how to compute $\mathsf{DA}$ given $\mathsf{SA}$ built for $T^{cat}$. Our algorithm runs in $O(N)$ time using $O(\sigma\lg N)$ bits of workspace, which is constant when $\sigma=O(1)$.

At a glance, we reuse the space of $\mathsf{SA}$ to compute the $\mathsf{LF}$-array, which is used to traverse $T^{cat}[1,N]$ from right-to-left applying the $\mathsf{LF}$-mapping $N$ times. We compute $\mathsf{DA}[1],\mathsf{DA}[\mathsf{LF}^{1}(1)],\mathsf{DA}[\mathsf{LF}^{2}(1)],\dots,\mathsf{DA}[\mathsf{LF}^{N-1}(1)]$. Starting with $doc=d+1$, each $\mathsf{DA}[i]$ receives $doc$, and whenever T^{cat}[\mathsf{SA}[i]-1]=\mathsf{BWT}[i]=\$$, that is, \mathsf{LF}[i]\in[2,d+1]$,$doc$ is decremented by one.

Recall that \mathsf{LF}(i,\)$is undefined then we cannot traverse backwards$T^{cat}[1,N]$with the$\mathsf{LF}-mapping given by Equation [2](#S2.E2). Alternatively, given the \mathsf{BWT}$of$T^{cat}$and an auxiliary array$\mathsf{C}[1,\sigma]$initialized with$\mathsf{C}[c]$equal to the number of symbols$c^{\prime}<c$in$T^{cat}[1,N]$, we can pre-compute correct$\mathsf{LF}$entries for every position with a corresponding$\mathsf{BWT}$symbol$c\neq$$. For$i=1,\dots,N$,$\mathsf{LF}[i]=\mathsf{C}[\mathsf{BWT}[i]]$, and$\mathsf{C}[\mathsf{BWT}[i]]$is incremented by one. The resulting (incorrect)$\mathsf{LF}$-positions, corresponding to$\mathsf{BWT}[i]=$$, will be in the interval$[2,d+1]. These values will be computed correctly by Algorithm [1](#algorithm1) on-the-fly during the right-to-left T^{cat}[1,N]$ traversal.

Algorithm 1

The algorithm starts with $\mathsf{SA}$ stored in $\mathsf{A}[1,N]$. We use $N\lg N$ bits to store $\mathsf{A}[1,N]$, and $N\lg(d+1)$ bits to store $\mathsf{DA}[1,N]$. First, we overwrite $\mathsf{SA}$ with the $\mathsf{BWT}$ in $\mathsf{A}[1,N]$ (Lines 1-3). Then, we overwrite the $\mathsf{BWT}$ with the $\mathsf{LF}$-array computed as described above (Lines 4-6). Recall that positions with $\mathsf{A}[i]\in[2,d+1]$ are not correct. In the sequel, $\mathsf{DA}[1,N]$ is computed while $\mathsf{SA}$ is reconstructed in the space of $\mathsf{A}[1,N]$ as follows. Initially, $pos=1$ and $doc=d+1$ (Lines 7-8). At each step $i=N,\dots,1$ (Lines 9-18), the value in $\mathsf{A}[pos]$ (corresponding to $\mathsf{LF}(pos)$) is stored in a temporary variable (Line 10) and replaced by $\mathsf{SA}[pos]=i$ (Line 11), then $\mathsf{DA}[pos]=doc$ (Line 12).

Whenever $\mathsf{LF}[pos]=tmp\in[2,d+1]$, $\mathsf{BWT}(pos)$ is a \$$-symbol and we have to compute correctly its \mathsf{LF}$-mapping.

In particular, when we reach the first $tmp\in[2,d+1]$, we reach the BWT position corresponding to the $d^{th}$ \$$-symbol in T^{cat}$(the last one), because we traverse$T^{cat}[1,N]$right-to-left, and its correct$\mathsf{LF}$-mapping is$tmp=d+1$. The next iteration we reach$tmp\in[2,d+1]$we are at the BWT position corresponding to the$(d-1)^{th}$$$$-symbol in$T^{cat}$, and$tmp=d$, and so on. Therefore, whenever$tmp\in[2,d+1]$we update$tmp$with the correct$\mathsf{LF}$-mapping value stored in$doc$(Line 14), and$doc$is decremented by one for the next iterations (Line 15). The next step will visit position$pos=tmp=\mathsf{LF}(pos)$(Line 17). At the end,$\mathsf{DA}[1,N]$is completely computed and$\mathsf{SA}$is reconstructed in the same space of$A[1,N]$.

Theoretical costs

The number of steps is $N$ and only array $\mathsf{C}[1,\sigma]$ was needed in addition to the input and output. Therefore, the algorithm runs in $O(N)$ time, using $O(\sigma\lg N)$ bits of workspace.

Discussion

We remark that one can use a standard suffix sorting algorithm (e.g. [15]) to compute the suffix array for $T^{cat}$, such that \$$-symbols are not considered different symbols in T^{cat}(see Section [2.1](#S2.SS1)), then\mathsf{LF}(i,$)$ is well-defined and Algorithm 1 can be applied with line 14 commented. Notice that, in this case, during suffix sorting unnecessary comparisons may be performed, depending on the order of the strings in the collection, what may deteriorate the practical performance of suffix sorting (see [21] for details).

4 Experimental results

We compared our algorithm with the lightweight alternative described in Section 2.3. We evaluated two versions of this procedure, using compressed (bit_sd) and plain bitvectors (bit_plain). We used C++ and SDSL library [25] version 2.0. The algorithms receive as input the concatenated string ($T^{cat}$) and its suffix array ($\mathsf{SA}$), which was computed using gSACA-K [21]. Our algorithm was implemented in ANSI C. The source codes are available at https://github.com/felipelouza/document-array/.

The experiments were conducted on a machine with Debian GNU/Linux 8 64 bits OS (kernel 3.16.0-4) with processor Intel Xeon E5-2630 v3 20M Cache $2.40$-GHz, $386$ GB of RAM and a $13$ TB SATA disk. We used real data collections described in Table 1.

Table 2 shows the running time (in seconds) and workspace (in KB) of each algorithm. The workspace is the peak space used subtracted by the space used for the input, $T^{cat}[1,N]$ and $\mathsf{SA}[1,N]$, and for the output, $\mathsf{DA}[1,N]$.

Results

bit_plain was the fastest algorithm in all tests. bit_plain was $2.19$ times faster than bit_sd, and $5.73$ times faster than Alg. 1, on the average. bit_sd was still $3$ times faster than Alg. 1, which shows that Alg. 1 is not competitive in practice. On the other hand, Alg. 1 was the only algorithm that kept the workspace constant, namely $1$ KB for inputs smaller than $2^{31}$ (2 GB) and $2$ KB otherwise, which correspond to the space used by the auxiliary array $C[1,\sigma]$ used to compute $\mathsf{LF}$. The workspace of bit_plain and bit_sd were much larger, bit_plain spent $0.16\times N$ bytes, whereas bit_sd spent $0.003\times N$ bytes, on the average.

Competing interests

The author declare that there is no competing interests.

Acknowledgments

We thank the anonymous reviewers for comments that improved the manuscript. We thank Giovanni Manzini, Travis Gagie and Nicola Prezza for helpful discussions.

Funding

F.A.L. was partially supported by the grants $\#$2017/09105-0 and $\#$2018/21509-2 from the São Paulo Research Foundation (FAPESP).