【BWT】Burrows-Wheeler Transform in C

圧縮に関しては素人ですが, 全文検索やデータ圧縮など様々な場面で使われる BWT/MTF や関連技術に興味があったので調べてみました。

Burrows-Wheeler Transform

BWT (Burrows-Wheeler Transform)はBlock Sortingとも呼ばれます。主にデータ圧縮の前処理として用いられます。

BWTでは, 接尾辞ではなく巡回シフトで得られた巡回行列を prefixで辞書順ソートしていきます。
ソート後の行列の最終列を取り出してます。

BWT後はデータの大きさは変換前とほぼ同じですが, 同じ文字が連続して出現されやすくなることから圧縮効率が高まります。

SA

SA (Suffix Array)は, ある文字列を先頭から一文字ずつ削って得られる接尾辞 (suffix)の文字列を配列 (array)に入れていきます。
効率的に文字列を検索できることから文字列探索, 全文検索などで使われ, namazu作者が書いたSaryでも使われています。

Move to Front

MTF (Move to Front)は先頭移動法という通り文字を先頭に移動していく方法です。これもBWT同様, より圧縮率を高めるような変換を行います。

MTFの歴史は古く, 1963年にTsetlinがロシアの学術雑誌に論文を発表しています。机の上の本の山から必要な本を一冊取り出し, 読んだら一番上に置くという”積ん読”のモデルとして導入されています。
読まれる確率の高い本が上に配置されやすいので, 誰かに勝手に掃除されてごちゃまぜにされると不便になってしまうことを証明したとも言えます。

圧縮技術におけるMTFの例として ‘aaabbccc’ という文字列を変換してみます。この文字列中に存在する文字は {abc} の3文字から成る集合を文字リストとします。
まず, 文字列 ‘aaabbccc’の先頭は a でこれは文字リストの1番目と同じです。結果の文字列に 0 を追加します。次に, 一文字ずつ移動して b の時は文字リスト2番目になります。
ここで, 文字リストの b を先頭に移動します。文字リストは {bac} になります。従って, 同じ文字が続く限り0が連続します。
これを繰り返すことで, 最終的に ‘00010200’ となり, MTF後の文字列は 0,1,2… と小さな数字が登場し易い偏った状態の文字列に変換されます。
この偏りを生むことで, エントロピー符号化による圧縮率が向上します。

BWT/MTFで変換したデータは, ランレングス圧縮やハフマン符号と相性が良いです。

BWTを実装してみた

BWTをCで実装してみました。アルゴリズム確認のため非効率な面があります。

BWT Implementation

まず, qsortでソートするための構造体を作ります。qsortはstdlibの関数です。
compare関数は qsortで使う比較関数です。

typedef struct {
	int  num;
	char str[MAX_LEN];
} Pair;

/*
* compare string
*
*/ 
int compare( const void *x, const void *y){
	return strcmp(((Pair*)x)->str ,((Pair*)y)->str);
}

目印となる特殊文字には$を使います。$でなくても良いですがソート時に一番小さい値にする必要があります。
先頭の文字から一文字ずつずらした文字列(巡回シフト)をPair構造体に格納していきます。
qsort後の Pairから suffixの文字を配列に入れて完成です。

/*
* create_suffix_array
*
*/ 
const char *create_suffix_array (char *str, int len){
     // add $ -> len+1
     Pair p[len+1];
     memset(p, 0, sizeof(p));
     char  strs[len+1];
     sprintf(strs,"%s%s",str,"

quot;); char tmp; // create_char_pair printf("-- Before Sort --\n"); for (int i = 0; i &lt; len+1; i++) { p[i].num = i; const char * ctmp = strs; strcpy(p[i].str,ctmp); // shift tmp=strs[len]; int j; for(j=len; j>0; j--){ strs[j]=strs[j-1]; } strs[j] = tmp; printf("%d : [%s]\n", p[i].num, p[i].str); } printf("-- Quick Sort --\n"); // quick_sort qsort(p, len+1, sizeof( Pair), compare); // create suffix array printf("-- After Sort --\n"); char SA[len]; memset(SA, 0, sizeof(SA)); for(int i = 0; i<len+1; i++){ const char *ptr = &p[i].str[len]; const char *preptr; if (i!=0){ preptr=&p[i-1].str[len]; } else{ preptr = ""; } // suffix add if (SA == NULL){ sprintf(SA,"%c", p[i].str[len]); }else{ sprintf(SA,"%s%c",SA, p[i].str[len]); } } const char *s = &SA[0]; return s; }

実際に使ってみます。
入力文字列は英語版wikiの LinuxKernelの冒頭部分です。

The Linux kernel is a Unix-like operating system kernel used by a variety of operating systems based on it, which are usually in the form of Linux distributions.The Linux kernel is a prominent example of free and open source software.

BWTの結果, 文字列全体の長さは 234文字($除く)で変わりませんが, 連続した文字が出現し易い傾向がわかります。

feeassyehsdxtefyllnxxmmeyddefaenggnlea,.txes   $.    ux  wvbrr i rieen hhehlkcrrssrnnnttpnppkkki ooo  onncTTtwrhrl mttLLLt  d ni   eeep-aleroaeoeiarrriiiUoeiii   sr i   fsooom eeuaaftaoeeepiimnua  yyiu  niss aausfeso  bnnn t uuuiebltss

Inverse BWT Implementation

続いて逆変換です。BWTには LF-mappingという性質があります。
LF-mappingは, 元の文字列の1文字前の位置を返すことでBWTから復元するために用いられます。
(Lは一番右の列, Fは一番左の列でLFの関係をmappingする)

文字列を一文字づつ Pair構造体に入れていきます。qsortでソートしたら特殊文字$を探し, そこから辿って元の文字列を復元していきます。

/*
* create_inverse_suffix_array
*
*/ 
const char *create_inverse_suffix_array (char *str, int len){
  printf("-- starting --\n");
  Pair p[len];
  memset(p, 0, sizeof(p));
  char  strs[len];
  sprintf(strs,"%s",str);
  char tmp;

  // create_char_pair
  printf("-- Before Sort --\n");
  for (int i = 0; i <len; i++)
  {
     p[i].num = i;
     char a = strs[i];
     const char * s = &a;
     strcpy(p[i].str,s);
     
     printf("{num: %d,str: %c}\n", p[i].num, a);
  }

  printf("-- Quick Sort --\n");
  qsort(p, len, sizeof( Pair), compare);

  // create suffix array
  printf("-- After Sort --\n");
  char SA[len];
  memset(SA, 0, sizeof(SA));
  
  char ary[len];
  int index,first=0;

  // get index
  for (int i = 0;i < len; i++)
  {
    char a = p[i].str[0];    
    if (a == '

#039;){ first = i; break; } } for(int i=0, index=first; i<len ; i++){ ary[i] = p[index].str[0]; printf("{num: %d,str: %c}\n",p[index].num,p[index].str[0]); index = p[index].num; } const char *s = &ary[0]; return s; }

MTFを実装してみた

コードは GitHub です。
以下の文字列を入力します。

aaabdlgagggassssggggggddddffffggggggcccccckgggsjskkkkkkkxxxxccppppoooooozzzkkjjjjjoooooo

MTFで変換した結果, 0の出現率が高くなっています。

0 0 0 1 3 11 7 4 1 0 0 1 18 0 0 0 2 0 0 0 0 0 4 0 0 0 8 0 0 0 2 0 0 0 0 0 7 0 0 0 0 0 12 2 0 0 5 12 1 3 0 0 0 0 0 0 23 0 0 0 5 0 17 0 0 0 17 0 0 0 0 0 25 0 0 5 0 7 0 0 0 0 3 0 0 0 0 0

圧縮アルゴリズムにかけてみた

圧縮前の前処理をしたので, 圧縮にかけてみます。
圧縮にも様々なアルゴリズムがあると思いますが, 比較的理解し易いのは以下だと思います。

ランレングス圧縮

‘aaabbccc’ の入力に対して ‘a3b2c3’ と連続した回数を数字で表現します。
連続した文字列になり易いBWTや, データがバイナリの場合だと回数だけになり圧縮効率が高いですが, データが連続していない場合, 元の文字列よりサイズが大きくなる欠点があります。

ランレングス圧縮も簡単なアルゴリズムなので, 実装して今回のBWT文字列を圧縮したところ元データの文字列に対して28%圧縮率が高くなりました。

エントロピー符号化

可変長符号化と呼ばれます。文字の出現率が異なるのに8bitで等しく1文字表現するのはもったいないことから, よく出現する文字には短い bit-codeを割当てることで効率よく表現します。
代表的な可変長符号化はハフマン符号化です。

LZ圧縮

別名, 増分分解法や動的辞書法とも呼ばれます。
繰り返し同じ文字列が出現する場合に有効な圧縮方法です。復元が早いのも特徴です。
gzip(GnuZip)では LZ77とハフマン符号化が使われています。

gzipでBWT前の元データを圧縮すると 169Byte, BWT後は 194Byte, BWT/MTF後は 174Byteと完敗でした…
こうなったらgzipのソースコード見て勉強します。

参考情報

今回のコードは GitHub にあります。
効率良くSAを構築できる Induced Sortingやハフマン符号化も書いてみたいです。今回, 参考にさせて頂いた本は『高速文字列解析の世界――データ圧縮・全文検索・テキストマイニング』です。

[1] 高速文字列解析の世界(SupportPage)
[2] ウェーブレット木解説
[3] 文法圧縮入門：超高速テキスト処理のためのデータ圧縮（NLP2014チュートリアル）