読者です 読者をやめる 読者になる 読者になる

BWT : bzip2

common lisp algorithm

gzipに続いて、bzip2圧縮(もしかしたら解凍も)common lisp(sbcl)で実装してみようと思う。

まだ詳しく調べてないので不確かだが、bzip2では「入力データ -> BWT*1 -> MTF変換 -> ハフマン符号 -> 圧縮データ」のようなプロセスで圧縮を行っているようだ。

今回は、この内のBWTを扱う。

BWT

BWTは圧縮そのものと云うよりは、入力データを圧縮しやすいものに変換する技法らしい。
詳細は、Wikipediaその他を見てもらえれば分かると思うので飛ばして、実装だけを載せる。

;; BWT前処理: 引数の文字列末尾にヌル文字を追加する
(defun append-null (str)
  (concatenate 'string str (list #\Null)))

;; 引数の文字列の先頭の文字を、末尾に移動する
(defun rotate (str)
  (let ((len (length str)))
    (concatenate 'string
                 (subseq str 1 len)
                 (list (char str 0)))))
> (rotate "引数の文字列")
--> "数の文字列引"

;; 文字列の循環列を生成する
(defun gen-all-rotate-pattern (str)
  (loop REPEAT (length str)
        FOR r = str THEN (rotate r)
        COLLECT r))
> (format t "~{~S~%~}" (gen-all-rotate-pattern "引数の文字列"))
"引数の文字列"
"数の文字列引"
"の文字列引数"
"文字列引数の"
"字列引数の文"
"列引数の文字"
--> NIL

;; bwt関数
;;  処理内容は、以下のようなもの (関数定義そのままだが...)
;;  1] 引数の文字列に終端(ヌル)文字を追加する
;;  2] 1の文字列の循環列を作成
;;  3] 2の文字列群を、string<でソート
;;     --> この時点で、先頭部分が等しい文字列は隣接することになる。
;;  4] 3でソートした文字列群の末尾の文字を集める
;;     --> ・末尾の文字 == 各文字列の先頭部分の文字列の一つ前の文字
;;         ・3で、先頭部分が等しい文字列は隣接するようになったので、
;;           この処理は、後ろの部分文字列が等しい文字を近くにまとめて取り出すことになる
;;         ・一般的に、後半部分が似ている文字列は、前半も似ている傾向があるので、
;;           この処理で同じ文字をだいたい隣接させられるようになる、らしい。
;;  5] 終了
(defun bwt (str)
  (flet ((last-char (s) (char s (1- (length s)))))

    (let ((rotated-strs (gen-all-rotate-pattern
                           (append-null str))))
       (map 'string #'last-char
                    (sort rotated-strs #'string<)))))

;; 例
;; 上のコメント(の一部)から改行・空白を除いたものにbwtを適用してみる。
> (defvar *text* "・末尾の文字 == 各文字列の先頭部分の文字列の一つ前の文字・3で、先頭部分が等しい文字列は隣接するようになったので、この処理は、後ろの部分文字列が等しい文字を近くにまとめて取り出すことになる・一般的に、後半部分が似ている文字列は、前半も似ている傾向があるので、この処理で同じ文字を隣接させられるようになる、らしい。")

;; 変換後の文字列では、同じ文字が隣接する傾向があることが分かる
> (bwt *text*)
--> "。字=・= ででるでははにいがしててししよよ向分列分近す、、接ら等等同出接さっ
    な一似似めのの3理こまににに的ううとくるた列列ここ尾前分ろ列理列にと半るる、
    せ取なあすれないいら後字字字る\Nullの・もがる、のののり部部部部字字字字字つ、
    前後て で傾文文文文文文文文文末、、隣隣のいじの分の各るい・処処般がが一を半
    頭頭のをは先先"

;; bwt関数定義内で、(sort ...)が返す文字列リストの中身は次のような感じ
( ;; ... 略 ... 
  ;; 以下は、先頭が"字"で始まる文字列群
  ;; 末尾の(= "字"の一つ前に位置していた)文字は、すべて"文" ※ この文字が、この後の処理で取り出される
  "字 ==    ... 文"
  "字を近く ... 文"
  "字を隣接 ... 文"
  "字・3で  ... 文"
  "字列が等 ... 文"
  "字列の一 ... 文"
  "字列の先 ... 文"
  "字列は、 ... 文"
  "字列は隣 ... 文"
  "尾の文字 ... 末"
  ;; ... 略 ...
)

bwt関数を『こころ』に適用してみた例。
参照: read-file

> (subseq (bwt (subseq (read-file "/path/to/kokoro") 5000 10000)) 200 400)
--> "にはででと時で時ではにかねてをは人もをで末にてのばて時らてとははとはやてがに
     がにはもたたただたたたたたたたたたたたるたたたたたたたたたたかかるたたにたた
     たたたよいたたたたたたたたたたたたたたたたたたたるたたたかたたたたたたたねる
     るたるるたたたたたたたたたたたたたたたたにたるたたたがたたたたたたたたたたた
     たたたたたかたたたたたたたたたたたたたたたたた種何折若何時色時段或或拶拶的後
     後暗好好入傷傷杏杏生"


次は、bwt関数が適用された文字列を、もとに戻す関数の定義。

;; 補助関数
;; 列の各要素に0から始まるインデックスを付与する関数を生成する
(defun indexing ()
  (let ((i -1))
    (lambda (x)
      (cons (incf i) x))))
> (mapcar (indexing) '(a b c d))
--> ((0 . A) (1 . B) (2 . C) (3 . D))

;; bwtの逆変換を行う関数
;;  1] bwt-strは、循環列の末尾の文字を並べたものなので、tailsとする
;;  2] bwt-strを、文字順でソートしたものは、元の循環列の先頭文字を並べたものと等しいので、headsとする
;;       --> この時点で、元々の循環列の先頭文字と末尾文字の情報が揃う
;;  4] 開始位置をtailsの#\Null文字の位置として、headsとtailsを終端文字に達するまで
;;     順番に辿っていけば、元の文字列が復元できる
;;       --> tails[n]の次の文字は、heads[n]のcdr部。heads[n]の次の文字は、tails[(car heads[n])]。その次の文字は...
(defun bwt-inverse (bwt-str)
  (let* ((tails bwt-str)      
         (heads (stable-sort   ; tailsとの対応(添字)を保存して、文字列を(安定)ソートする
                 (map 'vector (indexing) tails)
                 #'char< :key #'cdr))
         (start-idx (position #\Null tails)))
    (coerce
     (loop FOR (idx . ch) = (aref heads start-idx) THEN (aref heads idx)
           UNTIL (char= ch #\Null)
           COLLECT ch)
     'string)))

逆変換例。

;; bwt -> bwt-inverse
> (subseq (bwt-inverse (bwt (subseq (read-file "/path/to/kokoro") 5000 10000))) 200 400)
-->"のは不思議だといったりした。私は最後に先生に向かって、どこかで先生を見たように思うけれ
    ども、どうしても思い出せないといった。若い私はその時|暗《あん》に相手も私と同じような
    感じを持っていはしまいかと疑った。そうして腹の中で先生の返事を予期してかかった。ところ
    が先生はしばらく沈吟《ちんぎん》したあとで、「どうも君の顔には見覚《みおぼ》えがありま
    せんね。人違いじゃないですか」といったので私は変に一種の"

結構面白い。
分かってみると納得だが、bwt関数適用後の文字列から最初の文字列を復元できるのが、初めは不思議だった。

ちなみに、上の実装は返ってくる結果の分かりやすさを優先しており、恐ろしく非効率*2

*1:Burrow Wheeler's Transform

*2:上の『こころ』での例を何度か試していたら、ヒープ領域が足りなくなってしまった。この実装だと文字列の長さ^2のメモリが必要。