北海道大学で行った集中講義の資料を公開します．
http://researchmap.jp/sada/succinct/

部分文字列の頻度を求めるプログラムを作成しました．
サンプルプログラムのところからダウンロードできます．

文字列 T (長さ n) と文字列 P (長さ m) が与えられたとき，P の各部分文字列の T の中での
出現頻度を求めます．

実行は次のように行います．

./a.out -f1 pattern index.idx index.wxd

patternは P に対応するファイルで，index.idx と index.wxd は T の索引です．
-f1 は出力するパタンを指定するフラグで， -f0 から -f2 まであります．省略すると -f0 と同じです．

-f0 のとき，P の各位置 i に対し，そこから始まる部分文字列 P[i..j] (i <= j < n) に対し，
その出現頻度を表示します．ある i に対し，j を順に増やしていき頻度を求め，出現頻度が
0 になったら i を増やします．このときすでに頻度を求めた部分文字列に含まれるものについては
頻度を表示しません．例えば，P[0..0]，P[0..1]，P[0..2]，...，P[0..9] までの部分文字列が存在し，
P[0..10] が存在しない場合，P[1..9]，P[2..9]，...，P[9..9] は存在することは明らかなので飛ばし，
P[x..10] が存在するような最小の x を次の i とします．

-f1 のときは，-f0 の出力のうちで極大のものだけ表示します．つまり，上の例では P[0..9] のみ表示します．
-f2 のときは，出現する全ての部分文字列を表示します．

アルゴリズムの計算量は，O(md) です（圧縮接尾辞配列の各関数が定数時間としたとき）．ここで
d は P の部分文字列で T に存在するものの長さの最大値です．極端な場合，P と T が同じ文字列ならば
m = d = n なので，計算量は O(n^2) になってしまいます．なお，-f0 から -f2 の違いは結果を表示するか
しないかだけなので，計算量は同じです．

極大頻出パタン，つまり左極大かつ右極大な頻出パタンを求めるプログラムを作りました．
サンプルプログラムのところからダウンロードできます．

このプログラムでは，まず左極大パタンを全て求め，配列に格納します．
1つのパタンは (s, t, len) で表されます．s と t はパタンの辞書順の左端と右端，len はパタンの長さです．
左極大なパタン P が右極大でないということは，Pc というパタンが存在するということです．
前者の辞書順を [s1,t1]，後者の辞書順を [s2,t2] とすると，s1 <= s2 <= t2 <= t1 が成り立ちます．
つまり，入れ子になっています．また，共通の接頭辞を持たない2つのパタンに対しては，それらの
辞書順の区間は共通部分を持ちません．

以上のことから，パタンの辞書順をソートすると，共通の接頭辞を持つパタンを集めることができます．
ただし，ソートの基準は，(1) まずパタンの辞書順の左端の小さい順にソート，(2) 左端が同じ場合は
右端の大きい順にソート，(3) 両端が同じ場合は短い順にソート，とします．
ソートされたパタンを先頭から順に見ていくと，パタンの集合を表現する木構造の深さ優先探索を
シミュレートできます．パタンが入れ子になっていればそれらをスタックに積んでいきます．入れ子
でなければスタックの一番上のパタンが右極大となります．そして，新しいパタンと入れ子にならない
ものをスタックから捨てます．以上を繰り返すことで，右極大（かつ左極大）なパタンを列挙できます．

このアルゴリズムはデータを圧縮していない通常のアルゴリズムで，左極大なパタン数に比例する
時間で動きます．問題は左極大なパタンが何個あるかです．この数が多いと，全てをメモリに格納
することができなくなります．

パタンが頻出であることの定義を，頻度が n * t 以上（n は文字列長，t は 0 以上 1 以下の実数）
とします．すると，長さが k である頻出パタンは 1/t 個以下です．なぜならば，長さが同じである
2つの異なるパタンの辞書順の範囲は重ならないので，1つのパタンの範囲の幅が n*t 以上ならば，
そのようなパタンは 1/t 個以下しかないからです．長さが異なる2つのパタンの辞書順の範囲は
重なることがあるので，パタン長の最大値を M とすると，頻出パタンの数は M/t 個以下となります．
M と t を定数とすると，頻出パタンの数も定数となり，全てをメモリ内に格納することができます．

前回とほぼ同じですが，右極大頻出パタンを列挙するサンプルプログラムを作成しました．

左極大との違いは，基本的には child_l の代わりに child_r を用いたことだけです．

パタン P が，それを接頭辞に持つ接尾辞の辞書順 [l,r] で表現されているとします．
このとき，child_r 関数を用いると，Pc というパタンで存在するもの全ての辞書順を求めることができます．
ただし，child_l との違いは，child_r ではパタンの長さ m も指定する必要がある点です．

child_l と child_r は計算量も違います．出力サイズ (c の数) を k とすると，child_l の時間は k に比例しますが，
child_r の時間は k と m の積に比例します（その他に圧縮法によってはさらに時間がかかります）．

また，child_r を計算するにはΨを用いるので，索引がBWTに基づく場合にはΨをシミュレートする必要があります．

2種類の索引で左極大・右極大頻出パタンを求めてみました．データは前回と同じ NTCIR-4 PATENT で，
索引はBWTを用いたもの (-P12) とΨを用いたもの (-P11) です．出現頻度の閾値は文字列長の 0.0001 倍です．
用いた計算機はメモリが96GBあります．ただし実際に使われるのは索引サイズと同じ量です．

圧縮法	左極大	右極大
BWT(P12)	4.25	258.99
Ψ (P11)	48.12	48.87

表は実行時間（秒）を表します．BWTでの左極大が一番速いことがわかります．Ψでの左極大が遅い理由は，
child_lの実行時間がアルファベットサイズに必ず比例するからです．BWTの右極大が遅い理由は，child_rでは
パタン長に比例する時間がかかることと，BWTを用いてΨをシミュレートするときに，BWTでのselectを計算する
必要がありますが，それをrankを用いた2分探索で行っていることです．selectの索引を追加すれば高速になる
可能性があります．Ψでの右極大は左極大と同じような時間になっていますが，これは偶然だと思います．
右極大を求める場合はアルファベットサイズには依存しないので早いですが，パタン長に比例する時間はかかって
しまいます．それがたまたま釣り合っているのだと思います．

前回のΨの左極大はBWTの30倍ほど時間がかかっていましたが，それはデバッグ用の遅いルーチンを用いて
いたからです．高速化しましたがそれでも10倍以上時間がかかっています．

索引サイズはBWTが19.4GB，Ψが21.2GBです．このプログラムではパタンの出現位置は使用しないので
サンプルしたSAは不要です．

圧縮接尾辞配列ライブラリcsalibを更新しました．

変更点
・Ψを用いて検索するアルゴリズムのバグを修正．具体的には，csa.c の csa_psi_pred_naive と csa_psi_succ_naive．
・Ψを用いて検索するアルゴリズムの高速化．psi1.c の psi1_pred と psi1_succ．
・サンプルプログラム（左/右極大頻出パタン列挙）を追加．

圧縮接尾辞配列ライブラリのサンプルプログラムを作りました．
左極大頻出パタンを求めるプログラムです．

圧縮接尾辞配列ライブラリのサンプルプログラムのところにあります．
システムの都合でファイル名の拡張子に .c が使えないので，.txt にしてあります．

頻出パタンとは，テキスト中に指定された頻度以上出現するパタンのことです．
左極大頻出パタン P とは，P が頻出で，かつ全ての文字 c に対して cP が頻出ではないもののことです．
プログラムでは，そのようなパタンで，かつ長さが 100 以内のものを全て出力します．

頻出パタンには単調性があります．つまり，あるパタン P が頻出のとき，その任意の部分文字列も頻出です．
全ての頻出パタンを列挙すると数が多くなり，何が重要なのかわからなくなりますが，極大な頻出パタンだけ
求めることにすると，数が少なくなりわかりやすくなります．また，必要ならば極大な頻出パタンから全ての頻出
パタンを求めることができます．

なお，このプログラムで求めるのは左極大な頻出パタンです．右極大な頻出パタンも同様に定義でき，極大な
頻出パタンは左極大かつ右極大なものですが，右極大を効率よく求めるには接尾辞木のようなデータ構造が
必要なため，求めていません（左極大なパタンを全部メモリに記憶する方法ならばできます）．

アルゴリズムでは，頻出パタンの単調性を利用して，短いパタンから徐々に長いものを見つけていきます．
単調性より，P が頻出でなければ，どんな cP も頻出ではないため，そこで探索を打ち切ることができます．

パタン P は，それを接頭辞に持つ接尾辞の辞書順 [l,r] で表現されています．
ライブラリの child_l 関数を用いると，cP というパタンで存在するもの全ての辞書順を求めることができます．
その中で頻度が閾値以上のものに対して，再帰的に頻出パタンを求めていきます．
もし cP というパタンが頻出であるならば，P は頻出だが左極大ではないので P は出力しません．

P が頻出かつ左極大であったとき，それを表示しますが，入力ファイルがUnicodeであると仮定して，
パタンがUnicodeとして正しいものであったときのみそれを表示します．

コンパイルは次のように行います．

gcc leftfrequent.c csa.a

実行は次のようにおこないます．

./a.out 0.0001 filename.idx filename.wxd

引数の 0.0001 は頻出かどうかの閾値を表します．値が 1 より小さい場合には，検索している文字列
の長さ n に対する割合を表します．つまり，n * (割合) 回以上出現するパタンを求めます．
値が 1 以上の場合には頻度そのものを表します．

試しに NTCIR-4 PATENT 日本語データで実験してみました．これは1993-1997の日本の特許広報全文
データで，Unicodeで約113GBあります．これに対して -P12 の圧縮法で索引を作ると約21.6GBになります．
ちなみに同じテキストを bzip2 で圧縮すると約15.2GBです．

この索引を使って閾値 0.0001 で頻出パタンを求めると，約 900 個見つかります．
メモリ12GBのマシンで3秒以内です．なお，1回目の実行では索引をディスクから読み込むので
時間がかかりますが，2回目以降の実行では索引がメモリにキャッシュされるのでこの時間で終わります．

以下は見つかった頻出パタンの例です．
特許の文章に特徴的なパタンが見つかっています．

なお，child_l 関数を使っているので，Ψに基づく索引よりもBWTに基づく索引の方が高速です．
Ψだと実行時間は約30倍になります．
（追記：こんなに遅い原因は，Ψを計算するルーチンがデバッグ用の遅いものになっているから
だと思います．後日修正します．）

圧縮接尾辞配列ライブラリ csalib の使い方を書きました．

サンプルプログラムも今後公開する予定です．

dbwtにバグがありました．Ubuntuでは動かない可能性が高いです．
dbwt 2010-07-30を使ってください．

文字列を圧縮したまま検索するためのライブラリを公開しました．
詳しい説明はまだ用意できてないですが，ご自由に使ってください．
http://researchmap.jp/sada/csalib/

情報工学概論(アルゴリズムとデータ構造)の講義資料第13回(7月7日分)と

サンプルプログラム(Zアルゴリズム，接尾辞配列，C言語)を公開しました．