Canonical Huffman Codes

1999年出版と少し古い書籍ですが Managing Gigabytes を読んでいます。理解のために 2.3 で出て来る Canonical Huffman Codes の習作を作りました。

ハフマン符号は情報圧縮で利用される古典的なアルゴリズムで、圧縮対象データに出現するシンボルの出現確率が分かっているときに、その各シンボルに最適な符号長の接頭語符号を求めるものです。

通常のハフマン符号はポインタで結ばれたハフマン木を構築して、ツリーを辿りながら各シンボルに対する接頭語符号を計算します。このハフマン木には曖昧な箇所が残されています。ハフマン木は木の辺を右に辿るか左に辿るかで符号のビットが決まりますが、右が 0 で左が 1 などというのはどちらでも良いという点です。(曖昧だから駄目、という話ではありません。) 従って、ハフマン木から生成される符号は一意には決まりません。

ここで各シンボルに対応する接頭語符号長はそのままに「各符号の辞書式順序が文字の出現確率の大きい順序と一致する」(http://codezine.jp/article/detail/376) という制約を加えます。すると、符号が一意に決まるようになります。これが Canonical Huffman Codes とのこと。

符号が一意に決まるようになると、符号長だけが分かれば符号が計算できるようになります。よって、ハフマン木を持たずに最適な長さの接頭語符号を得ることができるようになります。(Canonical でない)ハフマン符号は符号化、復号化ともにハフマン木を必要としますが、このハフマン木がシンボルの種類が多くなると大きくなってしまいます。Canonical Huffman Codes の場合はハフマン木は不要で、幾つかの表があれば符号/復号が可能です。これらの表は空間的に非常に小さく、また CPU キャッシュを考慮すると時間的にも有利、とのことです。

以下実装です。余裕があればもう少し綺麗にしたいところ。少し整理しました。 (_min_heapfy() や _build_min_heap() はメソッドではなく CHC クラスの中だけで参照する関数にしたいのですが、Python ではどうしたらいいんでしょう? class 内で def で実装するとメソッドとしてしか呼べない?@staticmethod でいけました。id:hiratara さんありがとうございます。)

#!/usr/bin/env python
# -*- coding: utf-8 -*-

class CHC:
    @staticmethod
    def min_heapfy (a, i, h):
        l = 2 * i + 1
        r = 2 * i + 2

        if l < h and a[a[l]] < a[a[i]]:
            mi = l
        else:
            mi = i

        if r < h and a[a[r]] < a[a[mi]]:
            mi = r

        if mi != i:
            tmp   = a[i]
            a[i]  = a[mi]
            a[mi] = tmp
            CHC.min_heapfy(a, mi, h)

    @staticmethod
    def build_min_heap (a, h):
        for i in xrange((h - 1)/2, -1, -1):
            CHC.min_heapfy(a, i, h)

    def __init__(self, count):
        ## Phase One: Create an array A of 2n words
        n = len(count)
        A = [None] * (2 * n)
        for i in xrange(n):
            A[n + i] = count[i]
            A[i] = n + i
        CHC.build_min_heap(A, n)

        ## Phase Two: Building Huffman tree
        h = n
        while h - 1 > 0:
            m1 = A[0]
            A[0] = A[h - 1]
            h -= 1
            CHC.min_heapfy(A, 0, h)

            m2 = A[0]
            A[h] = A[m1] + A[m2]
            A[0] = h
            A[m1] = A[m2] = h
            CHC.min_heapfy(A, 0, h)

        ## Phase Three: Counting of leaf depths
        A[1] = 0
        for i in xrange(3, 2 * n):
            A[i] = A[A[i]] + 1

        ## codelen
        codelen = [None] * n
        for i in xrange(n):
            codelen[i] = A[n + i]

        maxlength = max(codelen)

        ## numl
        numl = [0] * (maxlength + 1)
        for i in xrange(n):
            numl[ codelen[i] ] += 1

        ## firstcode
        firstcode = [0] * (maxlength + 1)
        for l in xrange(maxlength - 1, 0, -1):
            firstcode[l] = (firstcode[l + 1] + numl[l + 1]) / 2

        ## codeword and symbol
        symbol = [None] * (maxlength + 1)
        for l in xrange(maxlength + 1):
            if numl[l] > 0:
                symbol[l] = [None] * numl[l] 

        nextcode = firstcode[:]
        codeword = [None] * n
        for i in xrange(n):
            l = codelen[i]
            codeword[i] = nextcode[l]
            symbol[l][nextcode[l] - firstcode[l]] = i
            nextcode[l] += 1

        self.codelen   = codelen
        self.codeword  = codeword
        self.firstcode = firstcode
        self.symbol    = symbol

    class Encoder:
        def __init__(self, l, w):
            self.codelen  = l
            self.codeword = w

        def encode(self, i):
            s = ""
            v = self.codeword[i]

            ## decimal to binary strings
            for x in xrange(self.codelen[i]):
                s = str(v % 2) + s
                v = v / 2
            return s

    class Decoder:
        def __init__(self, f, s):
            self.firstcode = f
            self.symbol    = s

        @staticmethod
        def nextinputbit(input):
            if len(input) > 0:
                return int(input.pop(0))
            else:
                return

        def decode(self, input):
            v = CHC.Decoder.nextinputbit(input)
            if v is None:
                return
            l = 1
            while v < self.firstcode[l]:
                v  = 2 * v  + CHC.Decoder.nextinputbit(input)
                l += 1
            return self.symbol[l][v - self.firstcode[l]]

#          a   b   c   d   e   f
count = [ 10, 11, 12, 13, 22, 23 ]
# count = [ 50, 11, 10, 13, 22, 23 ]

chc = CHC(count)
# print chc.codelen
# print chc.codeword

enc = CHC.Encoder(chc.codelen, chc.codeword)
for i in xrange(len(count)):
    print "%d: %s" % (i, enc.encode(i))

dec = CHC.Decoder(chc.firstcode, chc.symbol)
code = list("0111011010001") ## 34521
while True:
    i = dec.decode(code)
    if i is None:
        break
    print i,
% python chc.py
0: 000
1: 001
2: 010
3: 011
4: 10
5: 11
3 4 5 2 1

入出力のビット演算を実装する所は手を抜いて、二進の文字列で扱っています。

符号長を求めるアルゴリズムが面白かったです。シンボルの種類の2倍のサイズの領域を用意して、右半分に各シンボルの出現頻度を保存しておき、左半分の領域で右半分の各シンボルの添字を格納した min ヒープを構成します。この min ヒープを使って、通常のハフマン木を構築するのと同様に、出現頻度の小さい二つの葉をどんどんマージしていくと、この配列がハフマン木(相当)を表現することになり、ここから符号長を簡単に求めることができるというものです。

一度各種表が求まってしまえば、Canonical Huffman Codes の持つ性質から符号化/復号化は簡単です。特に復号が面白い。

解説をほとんど端折りましたが、時間に余裕があるときに実装含め改めて書き直したいと思います。

Managing Gigabytes: Compressing and Indexing Documents and Images, Second Edition (The Morgan Kaufmann Series in Multimedia Information and Systems)

参考文献

Canonical Huffman Code の仕組み, 実装は [1] を参考にしました。[2] はより発展的な話題を多く含んでいます。ヒープの実装は [3] を参考にしました。