Run-Length FM-Index
はじめに
FM-index は検索対象のテキストから予め索引を構築しておくことで，テキストに含まれる任意のパターン文字列の個数を数えるクエリ$ \mathrm{count} を，テキストのサイズに依らず高速に実行できるデータ構造です．加えて，suffix array や inverse suffix array の一部を追加で保持しておくことで，パターンの位置の列挙 $ \mathrm{locate}やテキストの復元 $ \mathrm{extract}といったクエリを高速に実行することができます（自己索引）． FM-Index に関しては，高速文字列解析の世界 や 簡潔データ構造 といった第一人者による書籍に紹介があり，またインターネット上にも id:echizen_tm さんのブログ を初めとする素晴らしい記事があります．なので日本語で読める情報には事欠かきません．しかし，FM-Index に連長圧縮を導入して圧縮効率の改善を図った run-length FM-index については，知る限りでは日本語の解説記事がありませんでした．本記事ではこのデータ構造について紹介します． 本記事で紹介するデータ構造は以下を基にしています：
Veli Mäkinen and Gonzalo Navarro. 2005. Succinct suffix arrays based on run-length encoding. Nordic J. of Computing 12, 1 (March 2005), 40-66.
高速文字列本に倣い，本記事における文字列および配列の添字は 0-based であるものとし，$ \mathrm{rank} や$ \mathrm{select} を以下のように定義します： $ \mathrm{rank}_c(X, i) :列$ X のスライス$ X[0..i] 中の$ cの出現回数を返す
$ \mathrm{select}_c(X, i):列$ Xにおいて$ i + 1番目に出現する$ cの位置を返す
FM-Index
FM-index の詳細については，最初に紹介した文献，もしくは拙著のメモをご覧ください． FM-index は以下2つのデータ構造から成ります．
$ \mathrm{rank}_c(L, i): $ Lのスライス$ L[0..i)中に出現する文字$ cの個数を数える
$ \mathrm{access}_c(L, i) : $ L[i] を得る
$ C: テキストに現れる各アルファベット$ cについて，テキスト中の$ cよりも小さい文字の出現回数を与える
テキスト$ T のBurrows-Wheeler 変換$ L は，$ T の suffix array$ SA をもとに
$ L[p] = \mathrm{if}\;SA[p]>0\;\mathrm{then}\;T[SA[p]-1]\;\mathrm{else}\; T[n-1]
とすれば求めることができますが，$ SAを生成せずに直接計算する手法も提案されています (Okanohara & Sadakane, 2009)．
また，FM-index 中には直接保存しませんが，今後の議論のために配列$ F を$ F[p] = T[SA[p]] として定義します．Suffix array の定義から，$ Fは$ Tを文字の大きさについて（安定）ソートしたものと一致します．下図は$ T = \mathrm{mississippi\$}について$ Lと$ Fを示したものです（太字が$ Tの suffix）．
https://gyazo.com/2a158e1b34e49a2f5d202a6eaade8abc
さて，$ Lの格納に wavelet tree もしくはその亜種，$ Cにアルファベットサイズ分の配列を用いることで，テキスト中のパターン文字列の出現回数を数えるクエリ$ \mathrm{count}を，パターンのサイズ$ mとアルファベットのサイズ$ \sigmaについて$ O(m\log\sigma)時間で計算することができます．アルゴリズムを以下に示します： code:count.rs
let mut s = 0;
let mut e = n;
for i in (0..n).rev() {
if s >= e {
break;
}
}
return e - s
繰り返しごとに suffix の prefix としてパターンが現れる範囲$ [s, e)が狭まっていくようなイメージです．
$ \mathrm{count} の計算に現れる式$ \mathrm{rank}_c(L, p) + C[c] は， $ L[p] = c のとき$ L上の LF-mapping$ LF(p) と等しくなります．LF-mapping は $ L におけるある文字が$ F のどの位置にあるかを与える関数です．例えば上図において$ T[5] = \mathrm{s} は$ L[8] 及び$ F[10] にあるので，$ LF(8) = 10 です．またその逆関数を本記事では inverse LF-mapping と呼びます．Inverse LF-mapping も FM-index によって$ \mathrm{select}_c(L, q - C[c]) とすることで求めることができます．
Run-Length FM-Index
概要
FM-index を構成する$ Lは，使用するデータ構造の種類によりますが，おおよそテキストの長さ$ nに比例する容量を必要とします．しかし Burrows-Wheeler 変換された文字列は比較的圧縮しやすいことが知られており，特にテキスト中に同じパターンが頻出していると$ Lには同じ文字が連続して出現しやすくなります．そこで，$ Lを連長圧縮 (run-length encoding) することで索引サイズを削減しつつ，クエリのオーバーヘッドを定数倍に抑えたものが run-length FM-index です．
構成
$ L中の連続する同一文字の列を run と呼びます．例えば$ L = \mathrm{xxxzyy\$xx}のとき，$ Lは$ \mathrm{xxx}, \mathrm{z}, \mathrm{yy}, \mathrm\$, \mathrm{xx}の5つの run からなります．
Run-length FM-index では，$ L の代わりに$ L の各 run の文字 (run head) ，および各 run の長さ (run length) を順番に並べたものを，それぞれ$ S と$ B として保存します．例えば run head の列は $ L = \mathrm{xxxzyy\$xx} の例だと$ S = \mathrm{xzy\$x} になります．$ S でも$ \mathrm{rank} や$ \mathrm{select} といったクエリが必要になるため，それらをサポートする構造を用いて保存します．論文では各文字ごとに出現位置を記録したビットベクトルを用いていますが，FM-index における$ L と同様に wavelet tree 等を用いる方法についても言及されています．
一方， run length を並べたものは同じ例だと $ [3, 1, 2, 1, 2] になります．これをそのまま格納する代わりに，各要素$ l について1つの$ 1 と$ l - 1 個の$ 0 をマップして連結した列を$ B とします．例えば$ L = \mathrm{xxxzyy\$xx} の例だと$ B = 100110110になります．ここで，$ Bの長さは$ Lの長さに等しく，$ Bに現れる$ 1の数は$ L中の run の数に等しいのがポイントです．$ Bは簡潔ビットベクトルを用いて格納します．
さらに， run length を並べた配列を，対応する run head の文字の大小関係をキーとして安定ソートしたものを考えます．同様の例だと$ [1, 3, 2, 2, 1] になります．これに$ Bと同様のマッピングを施し，$ B'として簡潔ビットベクトルを用いて格納します．例だと$ B' = 110010101になります．
https://gyazo.com/2abdd939fdc621c9c9c676085e2ed005
最後に，$ S に現れる各アルファベット$ c について，$ S 中の$ c より小さい文字の出現回数を数えたものを$ C_S とします．例では$ C_S[\$]=0,\ C_S[\mathrm x]=1,\ C_S[\mathrm y]=6,\ C_S[\mathrm z]=8 です．
FM-index で登場した$ L, Cの代わりに$ S, B, B', C_Sを保存したものが run-length FM-index です．
LF-Mapping
FM-index において$ \mathrm{count} の計算に用いた式 $ \mathrm{rank}_c(L, p) + C[c] を，run-length FM-index のために変形します．
まず， $ pが$ Lにおけるある run の先頭を指している状況を考えます（run の文字は任意）．このとき，任意の文字$ cについて以下の式が成り立ちます：
$ \mathrm{rank}_c(L, p) + C[c] = \mathrm{select}_1(B', C_S[c] + \mathrm{rank}_c(S, \mathrm{rank}_1(B, p)))
右辺の部分式はそれぞれ以下に対応します：
$ p: ある run （下図では$ \mathrm x） の$ Lにおける開始位置（run head の位置）
$ \mathrm{rank}_1(B, p): その run の run head の$ Sにおける位置
$ \mathrm{rank}_c(S, \mathrm{rank}_1(B, p)): $ Lにおいてその run より前に出現する，$ cでできた run の数
https://gyazo.com/8364c8f775ab5d20803fe6ee9aa14cc5
ここで左辺について，$ \mathrm{rank}_c(L, p) は$ L においてある run より前に現れる$ c の数，また$ C[c] は$ L （または$ F ）に現れる$ c より小さい文字の出現回数を表しています．よってそれらの和$ \mathrm{rank}_c(L, p) + C[c] は，その run の直前にある$ cの，$ Fにおける位置 + 1 を示しています（下図の☆部分）．
一方$ B' において，☆で示した位置の手前にある1の個数を数えると，$ c より小さい文字に対応する範囲には$ C_S[c] 個，文字$ c に対応する範囲には$ \mathrm{rank}_c(S, \mathrm{rank}_1(B, p)) 個（＝ ある run より前に出現する$ c でできた run の数）の 1 が出現することがわかります．したがって，それらの和$ C_S[c] + \mathrm{rank}_c(S, \mathrm{rank}_1(B, p) で$ \mathrm{select}を取ることで同じ☆の位置を得ることができます．
https://gyazo.com/c23535d70aa5d1e36a16606dce07627a
この結果を利用して，任意の$ p \in \{0, \ldots, n\} と$ c から$ \mathrm{rank}_c(L, p) + C[c] が計算できます．$ L[p] の値による場合分けになります．
$ L[p] \neq c のとき：
$ \mathrm{rank}_c(L, p) + C[c] = \mathrm{select}_1(B', C_S[c] + \mathrm{rank}_c(S, \mathrm{rank}_1(B, p)))
$ L[p] = c のとき：
$ \mathrm{rank}_c(L, p) + C[c] = \mathrm{select}_1(B', C_S[c] + \mathrm{rank}_c(S, \mathrm{rank}_1(B, p))) + p - \mathrm{select}_1(B, \mathrm{rank}_1(B, p))
まず，$ L[p] = c を判定するために$ L[p] の値を知る必要があります．ここで$ L[p] = S[ \mathrm{rank}_1(B, p + 1) - 1] が成り立つので，$ Sの要素が取得できれば判定が可能です．
$ L[p] \neq c のときは，最初の状況（$ p がある run の先頭位置）と同じに式になっています．$ L において位置$ p を通る run の先頭位置を$ p'とすると，$ \mathrm{rank}_c(L, p)=\mathrm{rank}_c(L, p'), $ \mathrm{rank}_1(B, p) = \mathrm{rank}_1(B, p') が成り立つので，
$ \mathrm{rank}_c(L, p) + C[c]
$ = \mathrm{rank}_c(L, p') + C[c]
$ = \mathrm{select}_1(B', C_S[c] + \mathrm{rank}_c(S, \mathrm{rank}_1(B, p')))
$ = \mathrm{select}_1(B', C_S[c] + \mathrm{rank}_c(S, \mathrm{rank}_1(B, p)))
となります．
また$ L[p] = c のときは，最初の式$ \mathrm{select}_1(B', C_S[c] + \mathrm{rank}_c(S, \mathrm{rank}_1(B, p))) が，$ pをまたぐ$ cの run の前にある$ cの$ Fにおける位置 + 1 ，すなわち$ pを通る run の run head （下図において$ L上の位置$ p'にある$ c）の$ Fにおける位置（$ q'）を指しています．
$ pを通る run の run head の位置は$ p' = \mathrm{select}_1(B, \mathrm{rank}_1(B, p))なので，$ p' - pを上式に足してやれば $ pの LF-mapping になります．
https://gyazo.com/4e8c349edab377d1ffa6bb99f4f2ee12
Inverse LF-Mapping
inverse LF-mapping $ \mathrm{select}_c(L, q - C[c]) も$ S, B, B', C_S から求めることができます．これは$ \mathrm{count} クエリの実行には必要ありませんが，自己索引として索引からテキストを正順に復元する際に便利です．以下は元論文の内容ではなくこの記事の著者による一つのアイディアです．もっとうまい方法があるかもしれません．
$ \mathrm{count} における LF-mapping と異なり，$ F[q] \neq c のケースを考慮する必要がないので$ F[q] = c の場合のみを考えます．
任意の$ q \in \{0, \ldots, n-1\} について，
$ \mathrm{select}_c(L, q - C[c]) = \mathrm{select}_1(B, \mathrm{select}_c(S, \mathrm{rank_1}(B', q + 1) - 1 - C_S[c])) + q - \mathrm{select}_1(B', \mathrm{rank_1}(B', q + 1) - 1 )
ただし$ c = F[q] ．
https://gyazo.com/cf342c54da76dc7499057cc6f6456320
$ F において$ q をまたぐ run を考えます．上図において，$ \mathrm{select}_c(L, q - C[c]) = p ，また$ F と$ L における run の先頭位置をそれぞれ$ q', p' とします．図のように，$ B'[0..q] の中には 1 が$ \mathrm{rank}_1(B', q+1) 個，うち$ c より小さい文字に対応する 1 が$ C_S[c] 個存在するので，その run の$ L における位置は$ p' = \mathrm{select}_1(B, \mathrm{select}_c(S, \mathrm{rank_1}(B', q + 1) - 1 - C_S[c])) になります（上式の第一項）．あとは run の先頭からの offset，すなわち$ q - q' = q - \mathrm{select}_1(B', \mathrm{rank_1}(B', q + 1) - 1 )を足せば終わりです．
$ c\ (= F[q]) の値は，FM-index では $ C[c] \leq q < C[c + 1] であるような$ c を二分探索することで求めることができました．Run-length FM-index においては，$ \mathrm{rank}_1 の単調性と$ C[c] = \mathrm{rank}_1(B', C[c] + 1) - 1 の関係から，$ C_S[c] \leq \mathrm{rank}_1(B', q + 1) - 1 < C_S[c+1] を満たすような$ cを二分探索することで求めることができます．
実装と評価
Rust で FM-index と run-length FM-index を実装しました．
$ Lや$ Sの保存には wavelet matrix (Claude & Navarro, 2012)，簡潔ビットベクトルは rsdic でも実装されている Navarro and Providel (2012) の手法を用いています． ここまで空間計算量に関する議論を避けていました．結果が用いる下部構造に依存するのと，この記事の筆者がその手の計算に不慣れなためです．$ Lが連長圧縮された一方で，新しいデータ$ B, B'が追加されているため常に圧縮が上がるとは限らなさそうです．そこで実際に索引を構築してサイズを比較してみます．
かなり artificial な例ですが，以下のような確率 $ pまたは$ 1 - pで次の状態に遷移して状態ラベルの文字を吐く遷移系によって出力されたテキストを用います．$ pを 1 に近づけることで同じ文字列$ 1234567が出現しやすくなります．
https://gyazo.com/521546fc7a70346770a33966af5ca5a2
テキストを 10,000,000 文字とし， $ pを変化させたときの FM-index と run-length FM-index のサイズは以下のようになりました：
table:result
p 索引サイズ (FM-index) (byte) 索引サイズ (run-length FM-index) (byte) run の数 / テキストの長さ
0.5 2182592 3007432 0.3359
0.9 1473872 1585144 0.1303
0.99 1010064 827208 0.0167
0.999 946088 650976 0.0017
Run-length FM-index が有利になるためには，繰り返して出現する文字列がそれなりに多い必要があるようです．
筆者の環境では，run-length FM-index は FM-index に比べて 3 ~ 4 倍程度のオーバーヘッドがある様子が見られます．それでも1,000,000 文字のテキスト中に現れるパターンの総数，パターン1文字あたり数 μs 程度のオーダーで数えることができています．ただベンチマークのアルファベットサイズは2文字なので，現実的にはもっと時間がかかります（計算量はアルファベットサイズ$ \sigmaに対して$ \log\sigmaに比例する）．
https://gyazo.com/5812eb24449c82c99b56826a821ff895
https://gyazo.com/922cf4582dd9ab734a28ffe198d0d789
発展
Run-length FM-index 中の$ Bや$ B'は， run の数が少ない理想的な状況下だとスパースな配列になるので，それに特化したビットベクトルを用いることでさらなる改善ができるかもしれません．また Gagie et al (2019) は，run-length FM-index の改良の一つとして，$ Bの表現に predecessor search に特化した構造を用いています．
さらに，著者らはパターンの位置を列挙する $ \mathrm{locate}クエリに必要な索引サイズをも run の数 $ rについて$ O(r)の空間計算量で抑える手法を同論文で提案しています（普通は suffix array を適当にサンプルして保存するので，時間計算量を抑えるためには$ O(n)の領域が必要）．実装は https://github.com/nicolaprezza/r-index に公開されています．本当はこちらを advent calendar で紹介する予定でしたが時間の都合上叶いませんでした．無念． まとめ
FM-index に連長圧縮を組み合わせた run-length FM-index を紹介しました．パターンの繰り返しが非常に多いテキストについて全文検索を行いたいときには採用を検討してみるといいかもしれません．
記事でおかしい点を発見された際は twitter 等でご連絡いただければ幸いです．PR も大歓迎です． 修正履歴
参考文献
Veli Mäkinen and Gonzalo Navarro. 2005. Succinct suffix arrays based on run-length encoding. Nordic J. of Computing 12, 1 (March 2005), 40-66.
Daisuke Okanohara and Kunihiko Sadakane. 2009. A Linear-Time Burrows-Wheeler Transform Using Induced Sorting. In Proceedings of the 16th International Symposium on String Processing and Information Retrieval (SPIRE '09), 90-101.
Francisco Claude and Gonzalo Navarro. 2012. The wavelet matrix. In Proceedings of the 19th international conference on String Processing and Information Retrieval (SPIRE'12), 167-179.
Travis Gagie, Gonzalo Navarro, and Nicola Prezza. 2018. Optimal-time text indexing in BWT-runs bounded space. In Proceedings of the Twenty-Ninth Annual ACM-SIAM Symposium on Discrete Algorithms (SODA '18), 1459-1477.