お気楽 Common Lisp プログラミング入門

ヒープとハフマン符号

今回は「ヒープ (heap)」というデータ構造を Common Lisp で作ってみましょう。そして、ヒープの応用例として「ハフマン符号」という古典的なデータ圧縮アルゴリズムを取り上げます。

●ヒープとは？

「ヒープ (heap)」は「半順序木 (partial ordered tree)」をベクタで実現したデータ構造です。一般的な二分木では、親よりも左側の子のほうが小さく、親よりも右側の子が大きい、という関係を満たすように作ります。

「半順序木」の場合、親は子より小さいか等しい、という関係を満たすように作ります。したがって、木の根 (ベクタの添字 0) には、必ず最小値のデータが格納されます。下図にヒープとベクタの関係を示します。

            0  1  2  3  4  5  6
    TABLE [10 20 30 40 50 60 70]

         (root)
           10 (0)
         ／   ＼            親の添字を k とすると
       ／       ＼          その子は 2*k+1, 2*k+2 になる。
     20 (1)       30 (2)    子の添字を k とすると
   ／  ＼       ／  ＼      その親は (k - 1) / 2 になる。
 40     50   60      70     親の値 <= 子の値 の関係を満たす。
 (3)    (4)  (5)     (6)

    図 : ヒープとベクタの対応関係

ヒープを利用すると、最小値をすぐに見つけることができ、新しくデータを挿入する場合も、高々要素の個数 (n) の対数 (log₂ n) に比例する程度の時間で済みます。

●ヒープの仕様

今回のプログラムで作成するヒープの操作関数を表に示します。

表 : ヒープのメソッド
関数名	機能
make-heap	ヒープを生成する
heap-push h x	ヒープ h にデータ x を追加する
heap-pop h	ヒープ h からデータを取り出す
heap-peek h	ヒープ h の先頭データを参照する
heap-length h	ヒープ h に格納されている要素数を返す
heap-clear h	ヒープ h を空にする
heap-emptyp h	ヒープ h が空ならば t を返す

●構造体の定義

それではプログラムを作りましょう。最初に構造体 heap を定義します。次のリストを見てください。

リスト : 構造体の定義

(defstruct heap
  (buff (make-array 8 :fill-pointer 0 :adjustable t))
  (key  #'identity)
  (obj> #'>))

スロット BUFF にはデータを格納するベクタをセットします。Common Lisp の場合、make-array で :adjustable に t を指定すると可変長ベクタとして利用することができます。詳しい説明は拙作のページ「配列」をお読みください。

KEY にはキーを取り出す関数を、OBJ> にはデータを比較する関数をセットします。OBJ> は第 1 引数が第 2 引数よりも大きいとき T を返す述語で、デフォルト値は #'> とします。これで小さいデータから順番に取り出すことができます。

●ヒープの構築 (1)

ヒープは、次の手順で作ることができます。

TABLE [＊ ＊ ＊ ＊ ＊ ＊ ＊ ＊ ＊ ＊]     最初は空

      [80 ＊ ＊ ＊ ＊ ＊ ＊ ＊ ＊ ＊]     最初のデータをセット

      [80 10 ＊ ＊ ＊ ＊ ＊ ＊ ＊ ＊]     次のデータをセットし親と比較
       親 子                              親の位置 0 = (1 - 1)/2

      [10 80 ＊ ＊ ＊ ＊ ＊ ＊ ＊ ＊]     順序が違っていたら交換

      [10 80 60 ＊ ＊ ＊ ＊ ＊ ＊ ＊]     データをセットし比較
       親    子                           親の位置 0 = (2 - 1)/2

      [10 80 60 20 ＊ ＊ ＊ ＊ ＊ ＊]     データをセットし比較
          親    子                        親の位置 1 = (3 - 1)/2

      [10 20 60 80 ＊ ＊ ＊ ＊ ＊ ＊]     交換する

      ・・・・データがなくなるまで繰り返す・・・・

                図 : ヒープの構築 (1)

まず、データを最後尾に追加します。そして、このデータがヒープの条件を満たしているかチェックします。もしも、条件を満たしていなければ、親と子を入れ換えて、次の親をチェックします。これを木のルート方向 (添字 0 の方向) に向かって繰り返します。条件を満たすか木のルート (添字 0) まで到達すれば処理を終了します。これをデータの個数だけ繰り返します。

このアルゴリズムを Common Lisp でプログラムすると、次のようになります。

リスト : ヒープの構築

;;; 要素の比較
(defun obj> (h x y)
  (funcall (heap-obj> h)
           (funcall (heap-key h) (aref (heap-buff h) x))
           (funcall (heap-key h) (aref (heap-buff h) y))))

;;; 要素の交換
(defun swap (buff x y)
  (psetf (aref buff x) (aref buff y)
         (aref buff y) (aref buff x)))

;;; ヒープの構築
(defun upheap (h n)
  (do ((p (floor (1- n) 2) (floor (1- n) 2)))
      ((or (minusp p) (not (obj> h p n))))
      (swap (heap-buff h) p n)
      (setf n p)))

関数 upheap はヒープを満たすように N 番目の要素をルート方向に向かって移動させます。0 から N - 1 番目までの要素はヒープの条件を満たしているものとします。N の親を P とすると、P は (N - 1) / 2 で求めることができます。そして、P が 0 以上で、かつ P の要素が N の要素よりも大きいのであれば、P と N の要素を交換して次の親子関係をチェックします。そうでなければ、ヒープの条件を満たしているので処理を終了します。

●ヒープの再構築

次に、最小値を取り出したあとで新しいデータを追加し、ヒープを再構築する手順を説明します。

TABLE [10 20 30 40 50 60 70 80 90 100]    ヒープを満たしている

      [＊ 20 30 40 50 60 70 80 90 100]    最小値を取り出す

      [66 20 30 40 50 60 70 80 90 100]    新しい値をセット

      [66 20 30 40 50 60 70 80 90 100]    小さい子と比較する
       ^  ^                               (2*0+1) < (2*0+2)
       親 子 子

      [20 66 30 40 50 60 70 80 90 100]    交換して次の子と比較
          ^     ^                         (2*1+1) < (2*1+2)
          親    子 子

      [20 40 30 66 50 60 70 80 90 100]    交換して次の子と比較
                ^        ^                (2*3+1) < (2*3+2)
                親       子 子            親が小さいから終了

                図 : ヒープの再構築

最初に、ヒープの最小値である添字 0 の位置にあるデータを取り出します。次に、その位置に新しいデータをセットし、ヒープの条件を満たしているかチェックします。ヒープの構築とは逆に、葉の方向 (添字の大きい方向) に向かってチェックしていきます。

まず、2 つの子の中で小さい方の子を選び、それと挿入したデータを比較します。もしも、ヒープの条件を満たしていなければ、親と子を交換して、その次の子と比較します。この処理を、ヒープの条件を満たすか子がなくなるまで繰り返します。

このアルゴリズムを Common Lisp でプログラムすると次のようになります。

リスト : ヒープの再構築

(defun downheap (h n nums)
  (do ((c (+ (* n 2) 1) (+ (* n 2) 1)))
      ((>= c nums))
      (when (and (< (1+ c) nums) (obj> h c (1+ c)))
        (incf c))
      (when (not (obj> h n c))
        (return))
      (swap (heap-buff h) n c)
      (setf n c)))

関数 downheap はヒープを満たすように N 番目の要素を葉の方向へ移動させます。N + 1 番目から最後までの要素はヒープの条件を満たしているものとします。最初に、N の子 C を求めます。これが NUM 以上であれば処理を終了します。もう一つの子 (C + 1) がある場合は、値が小さい方を選択します。そして、N の要素が C の要素よりも大きい場合はヒープの条件を満たしていないので、N 番目と C 番目の要素を交換して処理を繰り返します。、

なお、最小値を取り出したあと新しいデータを挿入しない場合は、新しいデータのかわりにベクタ BUFF の最後尾のデータを先頭にセットしてヒープを再構築します。上図の例でいえば、100 を BUFF の 0 番目にセットして、ヒープを再構築すればいいわけです。この場合、ヒープに格納されているデータの個数は一つ減ることになります。

●ヒープの構築 (2)

ところで、N 個のデータをヒープに構築する場合、N - 1 回 upheap を呼び出さなければいけません。ところが、すべてのデータをベクタに格納したあとで、ヒープを構築するうまい方法があります。次の図を見てください。

TABLE [100 90 80 70 60|50 40 30 20 10]    後ろ半分が葉に相当

      [100 90 80 70|60 50 40 30 20 10]    60 を挿入する
                    ^
      [100 90 80 70|60 50 40 30 20 10]    子供と比較する
                    ^              ^       (2*4+1), (2*4+2)
                    親             子

      [100 90 80 70|10 50 40 30 20 60]    交換する

      ・・・ 70 80 90 を順番に挿入し修正する ・・・

      [100|10 40 20 60 50 80 30 70 90]    90 を挿入し修正した

      [100 10 40 20 60 50 80 30 70 90]    100 を挿入、比較
        ^  ^  ^                           (2*0+1), (2*0+2)
        親 子 子

      [10 100 40 20 60 50 80 30 70 90]    小さい子と交換し比較
           ^     ^  ^                     (2*1+1), (2*1+2)
           親    子 子

      [10 20 40 100 60 50 80 30 70 90]    小さい子と交換し比較
                 ^           ^  ^         (2*3+1), (2*3+2)
                 親          子 子

      [10 20 40 30 60 50 80 100 70 90]    交換して終了

                図 : ヒープの構築 (2)

ベクタを前半と後半の 2 つに分けると、後半部分はデータがつながっていない葉の部分になります。つまり、後半部分の要素は互いに関係がなく、前半部分の親にあたる要素と関係しているだけなのです。したがって、後半部分だけを見れば、それはヒープを満たしていると考えることができます。

あとは、前半部分の要素に対して、葉の方向に向かってヒープの関係を満たすよう修正していけば、ベクタ全体がヒープを満たすことになります。興味のある方はプログラムを作ってみてください。

●操作関数の作成

次は操作関数 heap-push, heap-peek, heap-pop を作ります。次のリストを見てください。

リスト : 操作関数の定義

;;; データの追加
(defun heap-push (h x)
  (vector-push-extend x (heap-buff h))
  (upheap h (- (heap-length h) 1)))

;;; 先頭データの参照
(defun heap-peek (h)
  (if (heap-emptyp h)
      (error "heap : heap is empty")
    (aref (heap-buff h) 0)))

;;; データの取り出し
(defun heap-pop (h)
  (prog1
      (heap-peek h)
    (unless (heap-emptyp h)
      (let ((buff (heap-buff h)))
        (setf (aref buff 0) (vector-pop buff))
        (downheap h 0 (fill-pointer buff))))))

heep-push は vector-push-extend で可変長ベクタの末尾にデータ X を追加します。それから、関数 upheap で X をルート方向に移動してヒープを修正します。heap-peek は heap-emptyp を呼び出してヒープが空かチェックします。空の場合はエラーを送出します。データがある場合は可変長ベクタの 0 番目の要素を返します。

heap-pop は heap-peek で可変長ベクタの 0 番目の要素を求め、その後で 0 番目の要素を削除します。vector-pop で最後尾のデータを求め、それを 0 番目にセットします。そして、そのデータを関数 downheap でルートから葉の方向へ移動してヒープを修正します。あとの関数は簡単なので説明は割愛いたします。詳細はプログラムリスト１をお読みください。

●実行例

それでは簡単な実行例を示します。

リスト : 簡単なテスト

(require :heap "heap.lisp")
(use-package :heap)

(defun test0 ()
  (let ((h (make-heap)))
    (print (heap-length h))
    (print (heap-emptyp h))
    (dolist (x '(5 6 4 7 3 8 2 9 1 0))
      (heap-push h x))
    (print (heap-length h))
    (print (heap-emptyp h))
    (loop until (heap-emptyp h)
          do (print (heap-pop h)))))

(defun test1 ()
  (let ((h (make-heap :key #'car :obj> #'<)))
    (dolist (x '((5 a) (6 b) (4 c) (7 d) (0 e)
                 (3 f) (8 g) (2 h) (9 i) (1 j)))
      (heap-push h x))
    (print (heap-emptyp h))
    (print (heap-length h))
    (loop until (heap-emptyp h)
          do (print (heap-pop h)))))

* (test0)

0
T
10
NIL
0
1
2
3
4
5
6
7
8
9
NIL
* (test1)

NIL
10
(9 I)
(8 G)
(7 D)
(6 B)
(5 A)
(4 C)
(3 F)
(2 H)
(1 J)
(0 E)
NIL

このように、ヒープを使うと最小値 (または最大値) のデータを簡単に求めることができます。

●ハフマン符号

それでは、簡単な例題として「ハフマン符号」を取り上げます。ハフマン符号は 1952 年にハフマン (D. Huffman) が考案した、平均符号長を最小にすることができる符号化法です。古典的なデータ圧縮アルゴリズムですが、ほかのアルゴリズムと簡単に組み合わせることができるため、ハフマン符号は今でも現役のアルゴリズムです。最初にハフマン符号のアルゴリズムを簡単に説明します。

●ハフマン符号のアルゴリズム

ハフマン符号の構成は符号木を作ることで行います。ハフマン符号を構成するアルゴリズムを以下に示します。

各記号に対応する葉を作成する。この葉には、記号の出現頻度をあらかじめ格納しておく。
出現頻度の小さい方から 2 つの葉を取り出す。この葉を格納する新しい節を一つ作り、左右の枝に符号 0 と 1 を割り当てる。この節には 2 つの葉の出現頻度を足した値を格納し、新しい葉として追加する。
葉が一つになるまで手順 2 を繰り返すと、二分木を作成することができる。これをハフマン木と呼ぶ。根から記号に達するまでの枝をたどったときに得られる 0 と 1 の系列が、その記号の符号となる。

それでは、記号列 abccddeeeeffffgggggggghhhhhhhh を入力したときの、ハフマン符号化の具体的な構成例を示しましょう。

(８) (８) (４) (４) (２) (２) ─→(１) (１)
 ｈ   ｇ   ｆ   ｅ   ｄ   ｃ       ｂ   ａ

1. ａとｂを取り出す。

                                     N1
(８) (８) (４) (４) (２) (２)       (２)
 ｈ   ｇ   ｆ   ｅ   ｄ   ｃ       ／  ＼
                                (１)    (１)
                                 ｂ      ａ

2. 新しい節 N1 を作りａとｂを格納する。

                         N1
(８) (８) (４) (４)     (２)    (２) (２)
 ｈ   ｇ   ｆ   ｅ     ／  ＼    ｄ   ｃ
                    (１)    (１)
                     ｂ      ａ
3. N1 を登録する。

   図 : ハフマン符号の構成 (その１)

まず、各記号の出現頻度を求めて「節」の集合を構成します。この集合の中から、出現頻度の小さい方から 2 つ取り出して、新しい節に格納します。最初は、a と b を取り出して N1 に格納します。このとき、N1 の出現頻度は a と b を足した値をセットします。そして、この節 N1 を節の集合に登録します。この時点で節の集合は、{ c, d, N1, e, f, g, h } となります。あとは、この操作を節が一つになるまで繰り返します。

              N2                   N1
(８) (８)    (４)    (４) (４)    (２)
 ｈ   ｇ    ／  ＼    ｆ   ｅ    ／  ＼
         (２)    (２)         (１)    (１)
          ｄ      ｃ           ｂ      ａ

4. ｄとｃを取り出して新しい節 N2 を作る。

          N4                    N3
         (８)      (８) (８)   (６)
        ／  ＼      ｈ   ｇ   ／  ＼
     (４)    (４)          (２)    (４)
    ／  ＼    ｆ          ／  ＼    ｅ
 (２)    (２)          (１)    (１)
  ｄ      ｃ            ｂ      ａ

5. N1 とｅを取り出して新しい節 N3 を作る。
6. N2 とｆを取り出して新しい節 N4 を作る。

   図 : ハフマン符号の構成 (その２)

同様に、節の集合の中から d と c を取り出して、新しい節 N2 にセットして集合に登録します。節の集合は {N1, e, f, N2, g, h} となり、この中から頻度 2 の N1 と頻度 4 の e を取り出して N3 を登録します。すると、節の集合は {f, N2, N3, g, h} となり、その中から頻度 4 の N2 と f を取り出して N4 を登録します。

                   ROOT:N7          左の枝を 0 とすると
                  ／  ＼
                ／      ＼              a : 1101
           N6(16)        (14)N5         b : 1100
            ／  ＼      ／  ＼          c : 0001
         (８)    (８)(８)    (６)       d : 0000
        ／  ＼    ｈ  ｇ    ／  ＼      e : 111
     (４)    (４)        (２)    (４)   f : 001
    ／  ＼    ｆ        ／  ＼    ｅ    g : 10
 (２)    (２)        (１)    (１)       h : 01
  ｄ      ｃ          ｂ      ａ        平均符号長 = 80 / 30

 7. N3 とｇを取り出して新しい節 N5 を作る。
 8. N4 とｈを取り出して新しい節 N6 を作る。
 9. N5 と N6 を取り出して新しい節 N7 を作る。
10. 節が N7 の一つしかなくなったので終了。

           図 : ハフマン符号の構成 (その３)

この時点で節の集合は {N3, g, h, N4} の 4 つあります。小さい方から N3 と g を取り出して N5 を登録します。次に、h と N4 を取り出して N6 を登録します。節の集合は {N5, N6} となり、この 2 つを一つにまとめてハフマン木が完成します。

各記号の符号語は、ハフマン木の ROOT から葉に向かってたどっていくことで求めることができます。左右の枝にラベル 0 と 1 を割り当てることにすると、記号 a は「右、右、左、右」と枝をたどって葉に到達するので、符号語は 1101 となります。ほかの記号も同様に求めることができます。

なお、ハフマン符号は「葉」の組み合わせ方によって、異なる符号が得られます。しかしながら、どのハフマン符号でも同一の平均符号長が得られるので、圧縮率は同じになります。

●符号木の定義

それでは、ハフマン符号のプログラムを作りましょう。最初に符号木 (二分木) の節を定義します。

リスト : 符号木の節

(defstruct node
  (sym nil)     ; 記号
  (cnt 0)       ; 出現頻度
  (left nil)    ; 左の子
  (right nil))  ; 右の子

節を表す構造体は NODE としました。スロット SYM に記号を、CNT に出現回数をセットします。LEFT と RIGHT には左右の子を格納します。終端は NIL で表します。ハフマン符号は、ヒープを使うと簡単にプログラムを作ることができます。

●出現頻度表の作成

次は記号の出現頻度表を作成する関数 make-frequency を作ります。

リスト : 出現頻度表の作成

(defun make-frequency (xs &aux a)
  (dolist (x xs a)
    (let ((cell (assoc x a)))
      (if cell
          (incf (cdr cell))
        (push (cons x 1) a)))))

出現頻度表は連想リストで表します。コンスセルの CAR 部に記号を、CDR 部に出現回数を格納します。符号化するデータは引数 XS で受け取ります。dolist で XS から要素をひとつずつ取り出し、assoc で連想リスト A から記号を探索します。連想リスト内に記号 X がある場合、incf で出現回数を +1 します。見つからなかった場合、新しいセル (cons x 1) を生成して、連想リスト A に追加します。

それでは実際に実行してみましょう。

* (make-frequency '(a a a a b c c c d d))

((D . 2) (C . 3) (B . 1) (A . 4))
* (make-frequency '(a b c d a b c d e a b))

((E . 1) (D . 2) (C . 2) (B . 3) (A . 3))

●ハフマン木の生成

次は符号木を作る関数 make-huffman-tree を作ります。

リスト : ハフマン木の生成

(defun make-huffman-tree (ls)
  (if (null (cdr ls))
      (push (cons 'eof 0) ls))
  (let ((hp (make-heap :key #'node-cnt)))
    (dolist (x ls)
      (heap-push hp (make-node :sym (car x) :cnt (cdr x))))
    (do ()
        ((= (heap-length hp) 1) (heap-pop hp))
      (let ((a (heap-pop hp))
            (b (heap-pop hp)))
        (heap-push hp
                   (make-node :cnt (+ (node-cnt a) (node-cnt b))
                              :left a
                              :right b))))))

引数 LS には make-frequency で作成した出現頻度表 (連想リスト) を渡します。LS に要素がひとつしかないとハフマン木を構成できないので、ダミーのデータ (eof . 0) を追加します。次に、make-heap でヒープを生成して変数 HP にセットします。このとき、キーは記号の出現回数になるので、節 NODE から出現回数を求める関数 #'node-cnt を :key に指定します。それから、関数 make-node で節を生成し、heap-push でヒープに追加ます。

次に、ヒープからデータを取り出して、ハフマン木を構成します。ヒープにデータがひとつしかない場合、それがハフマン木のルートになります。heap-pop で節を取り出して返します。そうでなければ、ヒープから節を 2 つ取り出して変数 A と B にセットします。そして、新しい節を生成してヒープに追加します。このとき、A と B を左右の子にセットし、その節の出現回数は (+ (node-cnt a) (node-cnt b)) となります。これで、ハフマン木を構成することができます。

それでは、ここでハフマン木を表示する関数 print-huffman-tree を作成し、簡単なテストを行ってみましょう。

リスト : ハフマン木の表示

(defun print-huffman-tree (node &optional (n 0))
  (when node
    (print-huffman-tree (node-left node) (+ n 1))
    (dotimes (x n) (princ "    "))
    (princ (node-sym node))
    (terpri)
    (print-huffman-tree (node-right node) (+ n 1))))

* (print-huffman-tree (make-huffman-tree (make-frequency '(a a b a b c a b c d))))
    A
NIL
        B
    NIL
            D
        NIL
            C
NIL

このように、ハフマン符号では出現回数が多い記号ほど経路長 (符号語長) が短くなります。

●符号化と復号

最後に、符号化と復号を行う関数 huffman-encode と huffman-decode を作ります。符号化を行う関数 huffman-encode は次のようになります。

リスト : 符号化処理

;;; ハフマン符号を求める
(defun make-huffman-code (node cs code)
  (if (leafp node)
      (cons (cons (node-sym node) (reverse cs)) code)
    (make-huffman-code (node-right node)
                       (cons 1 cs)
                       (make-huffman-code (node-left node)
                                          (cons 0 cs)
                                          code))))

;;; 符号化
(defun huffman-encode (ls)
  (let* ((tree (make-huffman-tree (make-frequency ls)))
         (code (make-huffman-code tree '() '())))
    (values tree
            (loop for c in ls append (cdr (assoc c code))))))

関数 make-huffman-code は符号木を巡回して、記号と符号語を連想リストに格納して返します。CODE が連想リストで、CS が記号の符号語を表します。NODE が葉の場合、記号と符号語を cons でセルにまとめて CODE に追加します。そうでなければ、make-huffman-code を再帰呼び出しして符号木をたどります。左の枝をたどるときは CS に 0 を追加し、右の枝をたどるときは CS に 1 を追加します。

符号化を行う huffman-encode は簡単です。変数 TREE にハフマン木を、CODE にハフマン符号をセットします。あとは loop マクロで記号を符号語に変換するだけです。最後に values でハフマン木と符号を返します。

復号を行う関数 huffman-decode は次のようになります。

リスト : 復号

(defun decode-sub (tree node ls &optional a)
  (cond
   ((leafp node)
    (decode-sub tree tree ls (cons (node-sym node) a)))
   ((null ls)
    (reverse a))
   ((zerop (car ls))
    (decode-sub tree (node-left node) (cdr ls) a))
   (t
    (decode-sub tree (node-right node) (cdr ls) a))))

(defun huffman-decode (tree ls)
 (decode-sub tree tree ls))

引数 TREE がハフマン木で LS が符号を格納したリストです。局所関数 decode-sub でハフマン木をたどり、NODE が葉に到達したら記号を node-sym で取り出して累積変数 A にセットします。LS が空リストになれば復号は終了です。累積変数 A を reverse で反転して返します。あとは、符号 0 の場合は左の部分木をたどり、1 の場合は右の部分木をたどります。