2009年7月9日木曜日

Python: 階層的クラスタリングのデンドログラム描画と閾値による区分け

このブログ記事をはてなブックマークに追加

前回、「C++: マルチコアCPUを利用した並列化による高速な階層的クラスタリング」でクラスタリングを行ったのだが、ここではその出力データをPython+PILにより、デンドログラムを描画したり、指定した閾値で区分けを行ってみる。

まず、クラスタリングの出力データは前回のtest.outを使用する。内容は以下の通りだ。

0.0833487 4 14 0.11483 3 15 0.123895 0 5 0.126783 7 11 0.144271 16 -4 0.14854 8 9 0.253505 -5 -6 0.264889 -1 -3 0.301108 2 10 0.366858 6 -2 0.382649 13 -7 0.439469 17 19 0.588505 18 -8 0.648837 -10 -11 0.80762 -9 -13 1.03717 1 12 1.16488 -14 -15 1.46078 -12 -16 2.92199 -17 -18

これを下記のPythonのコードで処理すると、以下のようなデンドログラムが作成される。


使い方は、

draw_clusters.py test.out test.jpg

または、

draw_clusters.py test.out test.jpg labels.txt

とする。test.outは入力するクラスタデータ、test.jpgは出力されるデンドログラムの画像ファイルだ。また、二つ目に示した使い方では、labels.txtにデータ数分の文字列を入れておくことで、デンドログラムの数字を任意のラベルにすることができる。因みにこのデンドログラム表示プログラムについても、「集合知プログラミング」を参考にさせてもらった。

draw_clusters.py

#!/usr/bin/env python """ Draw dendrogram for hierarchical clustering by nox, 2009.07.02 """ import sys, os from PIL import Image, ImageDraw HEAD_LENGTH = 50 TAIL_LENGTH = 150 WIDTH = 1200 HEIGHT = 20 class bicluster: def __init__(self, left=None, right=None, distance=0.0, id=None): self.left = left self.right = right self.id = id self.distance = distance def find_cluster(clust, id): for i, c in enumerate(clust): if c.id == id: return i return None def read_clusters(fname, clust, labels): has_labels = False if labels: has_labels = True lines = file(fname).readlines() num = len(lines) + 1 for i in range(num): clust.append(bicluster(id=i)) if not has_labels: labels.append(i) for i, l in enumerate(lines): i = -(i + 1) data = l.strip().split() left_id = find_cluster(clust, int(data[1])) right_id = find_cluster(clust, int(data[2])) clust.append(bicluster(id=i, left=clust[left_id], right=clust[right_id], distance=float(data[0]))) del clust[right_id] del clust[left_id] def get_height(clust): if clust.left == None and clust.right == None: return 1 return get_height(clust.left) + get_height(clust.right) def get_depth(clust): if clust.left == None and clust.right == None: return 0 return max(get_depth(clust.left), get_depth(clust.right)) + clust.distance def draw_dendrogram(clust, labels, image_file): h = get_height(clust) * HEIGHT depth = clust.distance scaling = float(WIDTH - (HEAD_LENGTH + TAIL_LENGTH)) / depth img = Image.new("RGB", (WIDTH, h), (255, 255, 255)) draw = ImageDraw.Draw(img) draw.line((0, h / 2, HEAD_LENGTH, h / 2), fill=(255, 0, 0)) draw_node(draw, clust, HEAD_LENGTH, (h / 2), scaling, labels) img.save(image_file) def draw_node(draw, clust, x, y, scaling, labels): if clust.id < 0: h1 = get_height(clust.left) * HEIGHT h2 = get_height(clust.right) * HEIGHT top = y - (h1 + h2) / 2 bottom = y + (h1 + h2) / 2 ll = clust.distance * scaling lll = ll - clust.left.distance * scaling llr = ll - clust.right.distance * scaling draw.line((x, top + h1 / 2, x, bottom - h2 / 2), fill=(255, 0, 0)) draw.line((x, top + h1 / 2, x + lll, top + h1 / 2), fill=(255, 0, 0)) draw.line((x, bottom - h2 / 2, x + llr, bottom - h2 / 2), fill=(255, 0, 0)) draw_node(draw, clust.left, x + lll, top + h1 / 2, scaling, labels) draw_node(draw, clust.right, x + llr, bottom - h2 / 2, scaling, labels) else: draw.line((x, y, WIDTH - TAIL_LENGTH, y), fill=(255, 0, 0)) draw.text((WIDTH - TAIL_LENGTH + 5, y - 7), str(labels[clust.id]), (0, 0, 0)) def main(args): if len(args) < 2: print >>sys.stderr, "Draw dendrogram for hierarchical clustering." print >>sys.stderr print >>sys.stderr, " Usage: %s cluster_data_file [image_file=clusters.jpg label_data_file=[]]" % os.path.basename(args[0]) sys.exit(1) labels = [] if len(args) > 2: image_file = args[2] else: image_file = "clusters.jpg" if len(args) > 3: for l in file(args[3]): labels.append(l.strip()) clust = [] read_clusters(args[1], clust, labels) draw_dendrogram(clust[0], labels, image_file) if __name__ == "__main__": main(sys.argv)

コード内の変数、WIDTHを変更することで画像ファイルの幅を、HEIGHTを変更することで1データの高さを指定できる。初期値では幅1,200ピクセル、1データの高さを20ピクセルとしている。

次に、任意の閾値でクラスタデータを分けてみる。使い方は、

group_clusters.py test.out test.grp 1.0

のようにする。test.outを入力とし、出力ファイルにtest.grp、閾値を1.0としている。実行後の画面には以下のように表示され、

Input cluster file: test.out Output group file : test.grp Threshold : 1.000000 Number of elements: 20 Number of clusters: 5

出力されるtest.grpは、以下の内容となる。

001: 1 002: 12 003: 6 3 15 13 16 7 11 8 9 004: 2 10 18 4 14 0 5 005: 17 19

ソースコードを以下に示しておく

group_clusters.py

#!/usr/bin/env python """ Grouping clustering data by nox, 2009.07.02 """ import sys, os, math, re def expand_group(clust, n): group = [] for x in clust[n][1:]: if x >= 0: group.append(x) else: group.extend(expand_group(clust, -(x + 1))) return group def arrange_cluster(clust, threshold): tt = [[i, data[1:]] for i, data in enumerate(clust) if data[0] > threshold] t = [] for d in tt: for x in d[1]: if -x < tt[0][0] + 1: t.append(x) groups = [] for x in t: g = [] if x >= 0: g = [x] else: g = expand_group(clust, -(x + 1)) groups.append(g) if not groups: groups = [range(1, len(clust) + 2)] return groups def read_clust(fname, clust_data): for l in file(fname): data = l.strip().split() clust_data.append((float(data[0]), int(data[1]), int(data[2]))) def output_group(fname, groups): f = file(fname, "w") count = 1 for c in groups: f.write("%03d:" % count) for d in c: f.write(" %d" % d) f.write("\n") count += 1 f.close() def main(args): if len(args) < 4: print >>sys.stderr, "Grouping clustering data." print >>sys.stderr print >>sys.stderr, " Usage: %s cluster_data_file output_group_file threshold" % os.path.basename(args[0]) sys.exit(1) clust_file = args[1] group_file = args[2] threshold = float(args[3]) print "Input cluster file: %s" % clust_file print "Output group file : %s" % group_file print "Threshold : %f" % threshold clust_data = [] read_clust(clust_file, clust_data) groups = arrange_cluster(clust_data, threshold) print "Number of elements: %d" % sum([len(x) for x in groups]) print "Number of clusters: %d" % len(groups) output_group(group_file, groups) if __name__ == "__main__": main(sys.argv)

0 コメント: