極大部分文字列
Twitter で「極大部分文字列を求めるいいライブラリないかなー」とつぶやいていたら id:tkng さんに esaxx という岡野原さんのライブラリを教えてもらった。
esaxx というライブラリ名なのに説明が"stxx is ..."で始まったり、説明がところどころおかしい*1のはご愛敬として(最初は Suffix Tree のライブラリになるはずだったんだろうか)、確かにこれは便利そう。
早速、付属の "enumSubString.cpp" というサンプルをコンパイルして使ってみる。文字列はベタに "abracadabra"。
n:12 alpha:256 node:5 0 2 4 abra 1 5 1 a 2 2 3 bra 3 2 2 ra 4 12 0
あれ？
これは極大部分文字列ではなくて、Suffix Tree の内部ノードだ。
"abra"、"bra"、"ra" はそれぞれ同じ回数(2回)出現するので、"bra" と "ra" は極大部分文字列にはならない。
よく説明を読むと、確かに "... internal nodes in a suffix tree are enumerated ..." と書かれている。
岡野原さんの「全ての部分文字列を考慮した文書分類」を参照すると、
極大部分文字列の必要十分条件は，その部分文字列が内部節点に対応し，かつ，その節点に対応するBWT 配列が二種類以上の異なる文字を持つ場合である．
とある。
実際には、rank をきちんと計算する必要はなく、BWT が変化するたびにインクリメントする配列を用意すればいい。
BWT は、接尾辞配列のひとつ前の位置を見る。
概念的には T[SA[i] - 1] だが、マイナスになるのを防ぐために T[(SA[i] + size - 1) % size] のようにする。
BWT の変化を記録する配列は、esaxx トップページにあるサンプルのように Suffix Array を構築した後、次のようにして作る。
int size = T.size(); vector<int> rank(size); int r = 0; for (int i = 0; i < size; i++) { if (i == 0 || T[(SA[i] + size - 1) % size] != T[(SA[i - 1] + size - 1) % size]) { r++; } rank[i] = r; }
i 番目の内部ノードが極大部分文字列かどうかチェックする時には rank[ R[i] - 1 ] - rank[ L[i] ] を見て、0 より大きければ極大部分文字列。
最終的に、極大部分文字列を求めるサンプルコードは次のようになる。*2（esaxxのページのサンプルコードやenumSubString.cppのコードを使っています。読み込むファイルは "abra.txt" とハードコーディングしてあります。）
// License: MIT http://opensource.org/licenses/MIT #include "esa.hxx" #include <cstdio> #include <iostream> #include <vector> using std::cerr; using std::cout; using std::endl; using std::vector; int readFile(const char* fn, vector<char>& T){ FILE* fp = fopen(fn, "rb"); if (fp == NULL){ cerr << "cannot open " << fn << endl; return -1; } if (fseek(fp, 0, SEEK_END) != 0){ cerr << "cannot fseek " << fn << endl; fclose(fp); return -1; } int n = ftell(fp); rewind(fp); if (n < 0){ cerr << "cannot ftell " << fn << endl; fclose(fp); return -1; } T.resize(n); if (fread(&T[0], sizeof(unsigned char), (size_t)n, fp) != (size_t) n){ cerr << "fread error " << fn << endl; fclose(fp); return -1; } fclose(fp); return 0; } int main(void) { vector<char> T; // input text readFile("abra.txt", T); // read file into T int n = T.size(); vector<int> SA(n); // suffix array vector<int> L (n); // left boundaries of internal node vector<int> R (n); // right boundaries of internal node vector<int> D (n); // depths of internal node int alphaSize = 0x100; // This can be very large int nodeNum = 0; if (esaxx(T.begin(), SA.begin(), L.begin(), R.begin(), D.begin(), n, alphaSize, nodeNum) == -1){ return -1; } int size = T.size(); // BWT の変化を記録 vector<int> rank(size); int r = 0; for (int i = 0; i < size; i++) { if (i == 0 || T[(SA[i] + size - 1) % size] != T[(SA[i - 1] + size - 1) % size]) { r++; } rank[i] = r; } // 極大部分文字列を列挙 for (int i = 0; i < nodeNum; ++i) { int len = D[i]; if (len == 0 || (rank[R[i] - 1] - rank[L[i]] == 0)) { continue; } cout << R[i] - L[i] << "\t" << D[i] << "\t"; int begin = SA[L[i]]; for (int j = 0; j < len; ++j) { cout << T[begin + j]; } cout << endl; } }
"abra.txt" に "abracadabra" を保存し、上記のコード（test.cpp として esaxx のディレクトリに保存）をコンパイルして実行すると、次のようになる。
$ g++ test.cpp $ ./a.out 2 4 abra 5 1 a
1 列目は出現回数、2 列目は文字列の長さ。
ここでは 2 回以上出現するものは 1 文字であっても列挙しているが、それをスキップしたい場合は len < 2 で continue すればよい。
