1 1 7 9 8
圧縮接尾辞配列ライブラリ csalib
索引
索引とは，本の索引と同じ意味で，検索を高速に行うためのデータのことです．
ただし，本の索引では代表的な言葉のみが登録されていますが，このライブラリの索引は
任意の語が検索できるようになっています．

このライブラリの索引は自己索引 (self-index) と呼ばれるもので，索引自体に
元のファイルの情報を全て含んでいます．ですから検索時には元のファイルは
不要ですし，索引から元のファイルを復元することができます．
高速に検索できることを除けば，自己索引は gzip や bzip2 などで圧縮したものと同じ
役割を持ちます．

このライブラリの索引は filename.idx と filename.??? の2つからなります(3つの場合もあります)．
idxの方はファイルの長さ，文字の頻度，語の出現位置などを格納していて，
圧縮法には依存しません（異なる圧縮法を用いる場合でもidxは共通です）．

???の方は圧縮されたファイルを格納しています．拡張子は圧縮法によって異なります．
2つのファイルに分かれる場合もあります．
BW変換
検索したいファイルの索引を作る場合，まずファイルをBW変換する必要があります．
そして，BW変換されたファイル（以後BWT）から索引を作ります．

補足：なぜこのように2段階にしているかというと，圧縮接尾辞配列ではどの圧縮法でも
まずBW変換を行う必要があり，その部分が最も時間がかかるため，複数の圧縮法を
試す場合にはBW変換の部分を分けておいたほうが効率的だからです．また，BW変換は
圧縮されていない接尾辞配列の作成にも使えるため，分けておいたほうが汎用的です．

BW変換を行うプログラムは2種類用意してあります．dbwtとssssというプログラムです．

dbwtの方はメモリの中で動く通常のプログラムです．入力ファイルのサイズを n バイトとすると，
dbwtで必要なメモリは多くの場合 2.5n バイト以下と省メモリであり，かつ高速に動作します．
ただし，非常に大きなファイルではメモリが不足するため動きません．また，実装の都合上
4ギガバイト以上のファイルについては動作しません．
Unixの端末から  dbwt filename　として実行すると，output.bw とoutput.lst という2つのファイルを
出力します．

ssssは実行時にディスクを作業領域に使うプログラムです（external memory algorithmと呼ばれる）．
メモリが少なくても動作しますが，ディスクを約 30n バイト使用します．
Unixの端末から  ssss -M31 filename として実行すると，メモリを2ギガバイト (2^31) 使用し，
その他にディスクも 30n バイト使用します．-M を省略するとメモリを3ギガバイト使用します．
出力ファイルはdbwt同様，output.bw と output.lst です．

BWTは元のファイルの文字の並び替えなのでサイズは同じです．よってこれを圧縮する必要があります．
索引の作り方
BW変換されたファイルから，圧縮された索引を作ります．

BWTの名前を filename.bw と filename.lst とします．

mkcsa filename [-P[n]:{L}:{O}] [-I:{D}:{D2}]

filename は bw や lst を省いて指定します．

-P は圧縮法や圧縮のオプションを指定します．詳しくは圧縮法のところを参照してください．
例: -P12:512

-I は idx を作成する場合に指定します．
例: -I:128:256
2つの数字は，パタンの出現位置を格納する際のパラメタです．それらを D と D2 とします．
D は検索したパタンの出現位置を求める際の時間に関係するパラメタで，D を大きくすると
検索が遅くなる代わりに索引のサイズが小さくなります．
D2 はファイルの部分復元の時間に関係するパラメタで，D2 を大きくすると復元が遅くなる代わりに
索引のサイズが小さくなります．

補足：D は接尾辞配列を間引くためのパラメタで，接尾辞の辞書順が D の倍数のものについて
その出現位置を格納します．D2 は接尾辞配列の逆配列を間引くためのパラメタで，出現位置が
D2 の倍数の接尾辞についてその辞書順を格納します．n < 2^32 とすると，索引のサイズは
およそ 4n/D + 4n/D2 バイトとなります．4 の部分は n を表現できるバイト数なので，n > 2^32 では
5 以上になります．

D や D2 を 0 にすると，パタンの出現位置を求めたりファイルの部分復元をするための索引は作られません．
しかし，パタンの出現個数や辞書順は求めることができます．
また，部分復元はできませんが全体を復元することはできます．
通常の圧縮法との圧縮率の比較のために余計な索引を付けたくない場合にも使えます．

-I を省略し，-P だけにするとファイルの圧縮だけを行います．-P を省略して -I だけということはできません
（できるようにする予定です）．

圧縮法
圧縮接尾辞配列と一言で言っていますが，大きく分けて2種類の圧縮法が存在します．
compressed suffix array (CSA) と FM-index と呼ばれます．
両者は可能な操作は同じですが速度と圧縮率が違います．

CSAでは，接尾辞配列 SA[i] の代わりに，Ψ[i] = ISA[SA[i]+1] を格納します
（ISAは接尾辞配列の逆配列で，SA[i]= j のとき ISA[j] = i となります）．

FM-indexでは，SA[i] の代わりに BW[i] = T[SA[i]-1] を格納します．
つまりBWTをそのまま（ただし圧縮して）格納します．
BWを格納しておくと，LF[i] = ISA[SA[i]-1] が計算できます．

ΨとLFは互いに逆関数になっています．つまり Ψ[LF[i]] = LF[Ψ[i]] = i です．
また，片方のみを格納しておけば，その逆関数は効率的に計算できます．
このことから，CSAとFM-indexは機能が同じであることがわかります．
ただし速度は少し違います．

csalib の提供する圧縮法の詳細は以下の通りです．
{} に囲まれている部分は省略可能です．
全圧縮法において，L が小さいと検索・復元速度が速いですが索引サイズが大きくなります．

圧縮率と速度のバランスがいいものは
    -P4:512
    -P12:512
    -P3:512
あたりです．比較結果についても書く予定です．

-P1:{L}:{O}
    Ψは長さ L のブロックに分割され，ブロックの先頭ではΨ[i]をそのまま格納．
    先頭以外はΨを隣の値との差分（Ψ[i] - Ψ[i-1]）に変換し，差分をガンマ符号で符号化．
    Lは128から256程度の値が速度と圧縮率のバランスが良いです．
    Oは1にするとブロックの先頭へのポインタと先頭のΨの値も圧縮します(sparsearrayを用いる)．
    索引のファイル名は .psd と .psi．実行時には .psd の方だけを指定すると
    自動的に .psi も読み込みます．

-P2:{L}:{O}
    -P1とほぼ同じですが，差分の値で1が連続する場合にはそれを連長圧縮します．
    圧縮率の高いファイルでは-P1よりもさらに小さくなります．
    索引のファイル名は .prd と .pri．実行時には .prd の方だけを指定すると
    自動的に .pri も読み込みます．

-P3:{L}
    DNA配列専用．入力ファイルは acgt の4種類以外の文字を含んではいけません．
    索引のファイル名は .bwd と .bw2．実行時には .bwd の方だけを指定すると
    自動的に .bw2 も読み込みます．
    .bw2 はBWTを1文字あたり2ビットで表現したものです．DNA配列は圧縮が難しく，
    圧縮しても検索速度が遅くなるだけで無理に圧縮はしていません．そのかわり
    非常に高速に検索を行うことができます．bwa や SOAP2 と同様の手法です．
    L はBWTでのrankを何文字おきに格納するかのパラメタです．256か512がいいでしょう．

-P4:{L}:{O}
    BWTをウェーブレット木で表現したもの．ウェーブレット木を表現するビットベクトルは
    comparrayで圧縮されている．
    L はビットベクトルでのrankを計算するための索引のパラメタで，ビットベクトルを
    長さ L のブロックに分割し，各ブロックの先頭でのrankの値を格納している．
    Lは512程度がいいでしょう．
    Oが1のとき，ブロックへのポインタ等を圧縮します．2のとき，ベクトルを1の連続する数と0の連続する数の
    ガンマ符号で表現します．4のとき，小ブロックの1の数と実際の符号を分離して格納します．
    索引のファイル名は .wtd です．

-P5:{L}:{O}
    -P4とほぼ同じですが，小ブロックの1の数を固定長ではなくハフマン符号で圧縮します．
    -P4よりも圧縮率が良くなります．

-P6
    Ψをsparsearrayで符号化します．復号は速いですが圧縮率は良くありません．

-P7:{L}
    BWTのウェーブレット木で表現し，そのビットベクトルを無圧縮で（densearrayで）格納します．
    L はブロックのサイズです．

-P8
    -P1とほぼ同じですが，ブロックの先頭の値とポインタをsparsearrayで格納します．
    索引のファイル名は .pss と .psi．.psi は -P1 と共通です．
    実行時には .pss の方だけを指定すると自動的に .psi も読み込みます．

-P9
    -P2とほぼ同じですが，ブロックの先頭の値とポインタをsparsearrayで格納します．
    索引のファイル名は .prs と .pri．.psi は -P1 と共通です．
    実行時には .prs の方だけを指定すると自動的に .pri も読み込みます．

-P10
    BWTのウェーブレット木で表現し，そのビットベクトルをsparsearrayで格納します．
    索引のファイル名は .wsa です．
 
-P11:{L}:{O}
    Ψをビットベクトルで表現しているとみなし，それを連続する1の数と0の数のガンマ符号で表現します．
    L はブロックの長さです．O が 1 のとき，ブロックの先頭の値をsparsearrayで圧縮します．
    索引のファイル名は .pxd と .pxi です．
    実行時には .pxd の方だけを指定すると自動的に .pxi も読み込みます．

-P12:{L}:{O}
    BWTのウェーブレット木で表現し，そのビットベクトルを連続する1の数と0の数のガンマ符号で表現します．
    L と O は -P4 と同じです．
    索引のファイル名は .wxd です．

-P15:{L}
    ビットベクトル専用．入力ファイルは 0 と 1 以外の文字を含んではいけません．
    索引のファイル名は .bwb と .bw1．実行時には .bwb の方だけを指定すると
    自動的に .bw1 も読み込みます．
    .bw1 はBWTを無圧縮のビットベクトルで表現したものです．
    L はBWTでのrankを何文字おきに格納するかのパラメタです．
コンパイルの仕方
単に make するだけですが，いくつか注意点があります．

4ギガバイト以上のファイルを扱うため，ライブラリ内ではほとんどの値は64ビットで表現されています．
そのため，OSも64ビット板の方がいいでしょう．32ビット板では動作確認はしていません．
また，ファイルポインタを64ビットにするオプションもつけています．

CPUのエンディアンには依存しないはずですが，リトルエンディアンでしか動作確認していません．

-P3 と -P15 では，インテルのSSE4.2命令（popcount）を使用して高速化しています．
コンパイルするには，gcc-4.3 以降が必要です．MakefileのCCの項でコンパイラを指定してください．
gcc-4.3以降が使えない場合にはSSE4.2命令は使えません．また，実行したいCPUがSSE4.2命令を
サポートしていない場合には不正命令となり実行できないので，古いgccでコンパイルする必要があります．
実行時にSSE4.2命令が使用可能か判定して，使える場合のみ使うようにする予定です．

コンパイルして作成されるファイルは次のものです．
    csa.a     圧縮接尾辞配列ライブラリ本体．自作のCプログラムではこれをリンクして使います．
    mkcsa   BWTから索引を作成するプログラム．「索引の作り方」を参照してください．
    unbwt    索引から元のファイルを復元するプログラム．出力ファイル名は output.dec です．

    csa        索引を使った検索のサンプルプログラム
    mkcst    圧縮接尾辞木を作成するサンプルプログラム
    approx   近似文字列照合を行うサンプルプログラム
C言語でのライブラリの利用
自作のCプログラムからライブラリを利用する場合，csa.h をインクルードし，csa.a をリンクする必要があります．

まず，索引を読み込み，それを構造体で表します．
int main(int argc, char *argv[]) {
CSA SA; // 読み込んだ索引を表す構造体
csa_read(&SA,argc-1, argv+1); // 索引の読み込み
csa_readの2番目の引数は読み込む索引の数（通常は2），3番目の引数は索引のファイル名の配列です．
この例ではプログラムの実行時の引数をそのまま渡しています．ただし，argv[0] はプログラム名が入っていて
実際の引数は argv[1] からなので，argv を1増やしています．同じ理由で argc は1減らしています．

パタンの検索は次のように行います．
i64 s, l, r;
uchar *key;
s = SA.search(key,strlen(key),&SA,&l,&r);
検索したいパタン key と，その長さを指定します．
結果は関数の返り値 s と，l と r に入ります．
s と key の長さが等しいならば，key が存在することになります．その際，key から始まる接尾辞の辞書順の
最小値と最大値が l と r に格納されています．つまり，SA[l], SA[l+1],..., SA[r] の接尾辞は，先頭が key になっています．また，key の出現回数は r-l+1 であることがわかります．
s が key の長さより小さい場合，key の後ろから s 文字のパタンが見つかったことを表します．
つまりkeyの最長一致接尾辞を求めていることになります．この場合の l と r は，見つかった s 文字のパタンの
辞書順を表しています．完全一致を見つける場合は，必ず返り値の値を確認してください．

索引に格納されているファイルの長さを n とすると，l と r は1 から n の範囲です．

関数は構造体の中に格納されています．これは，索引ごとに実際の関数は異なるため，そうしないと
索引ごとに異なるプログラムが必要になってしまうからです．

見つかった key の出現位置を求める場合は次のようにします．
for (i=l; i <= r; i++) {
  j = SA.lookup(&SA,i);
  printf("SA[%ld] = %ld \n", i, j);
}
引数 i は 0 から n まで，返り値は 0 から n までです．
なお，常に SA[0] = n となっています．ファイルは T[0..n-1] ですが，その後の T[n] に終端記号があり，
それはどの文字よりも小さいという仮定をしています．

ファイルの一部分を復元するときは次のようにします．
uchar buf[100];
SA.text(buf,&SA,s,t);
これでファイルの s 文字目から t 文字目までが buf に格納されます．
buf は当然 t-s+1 バイト以上必要です．C言語の文字列とは違い，最後に 0 を書かないので注意してください．
s と t は 0 から n-1 の範囲です．

接尾辞配列の置き換えとして使うならばこれだけでも可能ですが，もっと複雑な問い合わせもできます．
i = SA.inverse(&SA,j);
SAの逆関数 i = ISA[j] を求めます．これは接尾辞 T[j..n] の辞書順が i であることを表します．

c = SA.T(&SA, i);
辞書順が i の接尾辞の先頭の文字 c を求めます．つまり c = T[SA[i]] です．
i は 0 から n までですが，i=0 のときは終端記号なので c = -1 となり，通常の文字ではないので注意してください．
なお，T[SA[i]] は i が増えると単調増加します．これは接尾辞が辞書順にソートされているからです．

c = SA.BW(&SA, i);
BWTの i 番目の文字を求めます．つまり c = BW[i] = T[SA[i]-1] です．
BW[0] はファイルの一番最後の文字になります．
また，一番長い接尾辞 T[0..n] (ファイル全体) の辞書順を x とすると，BW[x] = -1 です．

j = SA.psi(&SA,i);
j = Ψ[i] = ISA[SA[i]+1] を計算します．SA[i] = s とすると，SA[j] = s+1 になっています．
つまり，接尾辞 T[s..n] の辞書順から1つ短い接尾辞 T[s+1..n] の辞書順を計算できます．
このとき s の値は必要ありません．
i の範囲は 0 から n．x = Ψ[0] とすると，x は T[0..n] の辞書順になります．つまり BW[x] = -1．

j = SA.LF(&SA,i);
j = LF[i] = ISA[SA[i]-1] を計算します．

なお，Ψを用いた索引，BWTを用いた索引どちらの場合でも，psi, LF, BW は全て計算できます．
ただしΨを用いた索引でLFやBWを計算すると，アルファベットサイズに比例した時間がかかります．
BWTを用いた索引でpsiを計算すると，log n に比例した時間がかかります．

なお，ファイルの部分復元は，inverse, T, psi （またはLF）で実現されています．

SA.searchsub(c, &SA, &l, &r);
パタン P の辞書順 [l..r] がすでに求まっているとします．このとき，P の直前に文字 c をつけた
cP というパタンの辞書順を求めます．結果は l と r に上書きされます．
もし l > r であればパタン cP は存在しないことを意味します．

k = csa_search_r(len, c, &SA,&l, &r);
パタン P の辞書順 [l..r] がすでに求まっているとします．またパタンの長さを len とします．
このとき，Pの右に文字 c をつけた Pc というパタンを検索します．
返り値 k は，Pc の接尾辞でファイル中に存在するものの長さの最大値です．
k = len+1 のとき，Pc が存在することになります．l と r は Pc の辞書順を表します．


uchar head[SIGMA];
i64 L[SIGMA], R[SIGMA];
k = SA.child_l(&SA,l,r,head, L, R);
パタン P の辞書順 [l..r] がすでに求まっているとします．このとき，cP というパタンでファイル中に存在するもの
全てに対し，そのときの c と cP の辞書順を配列 head, L, R に格納します．配列のサイズは SIGMA = 256
以上ある必要があります．関数の返り値 k は存在する cP の種類数となります．これは 0 から 255 までの
全ての c に対して searchsub を実行した場合と出力は同じですが，索引がBWTを用いたものの場合には
k に比例した時間で実行できるため高速です．

uchar head[SIGMA];
i64 L[SIGMA], R[SIGMA];
k = SA.child_r(&SA,len,l,r,head, L, R);
パタン P の辞書順 [l..r] がすでに求まっているとします．またパタンの長さを len とします．
このとき，Pの右に文字 c をつけた Pc というパタンでファイル中に存在するもの全てを求めます．
出力は child_l と同じです．k に比例した時間で実行できます．
Rubyでのライブラリの利用
Rubyの拡張ライブラリとしても利用できます．

ソースコードのrubyディレクトリでコンパイルします．
Pizza&Chili Corpus API
Pizza&Chili CorpusのAPIにも準拠しています．

ソースコードのpizzachiliディレクトリでコンパイルします．
その際，Pizza&Chili Corpusのサイトの interface.h と run_queries.c が必要です．
コンパイルすると csa_run というプログラムが作成されます．これはPizza&Chiliでの
検索速度のベンチマーク用の標準的なプログラムです．
サンプルプログラム

サンプルプログラム

{{cabinetFile.CabinetFile.filename}} >  サンプルプログラム
名前 更新日
leftfrequent_c.txt 605
左極大頻出パタン列挙
2010/08/06
maximalfrequent_c.txt 1283
極大頻出パタン列挙
2010/08/12
rightfrequent_c.txt 593
2010/08/10
substringfreq_c.txt 1078
2010/08/20