2011-01-20

点線で楕円を描く・まとめ

math

【追記・2011/01/22】最初に書いた内容に大幅に誤りがあったので, 再度編集しています. おっちょこちょいですみません.

【追記・2011/02/03】修正完了. ソースコードを色々変えました.

昨日、一昨日の記事には数式に間違いがあったので訂正したまとめ記事を書きます。

【修正・2013/03/12】数式に対するコメントの通り, 修正しました.

楕円の弧長

まず、 $x = a \cos \varphi,\quad y = b \sin \varphi$ とパラメータ表示して、点 (a, 0) から左回りに弧長 $u(\varphi)$ を測っていくことにします。

$\begin{align*} u &= \int_0^{\varphi} \sqrt{dx^2 + dy ^2} \\ &= \int_0^{\varphi} \sqrt{a^2 \sin^2 \varphi + b^2 \cos^2 \varphi} d\varphi \end{align*}$

となり、この逆関数 $\varphi (u)$ の微分を考えます。(ここの時点で式の形に間違いがありましたね.)

$\frac{d \varphi}{du} = \frac{1}{\frac{du}{d \varphi}} = \frac{1}{\sqrt{a^2 \sin^2 \varphi + b^2 \cos^2 \varphi}}$

これを差分形にして、

$\begin{align*} \Delta \varphi &= \frac{\Delta u}{\sqrt{a^2 \sin^2 \varphi + b^2 \cos^2 \varphi}} \\ \varphi_0 &= 0 \\ \varphi_n &= \varphi_{n - 1} + \frac{\Delta u}{\sqrt{a^2 \sin^2 \varphi_{n-1} + b^2 \cos^2 \varphi_{n-1}}} \end{align*}$

$\sin$ と $\cos$ を間違えていましたね。

ここで $a$ と $b$ の大小については何も言っていないことに注意してください。

これを実装します。

Python ソースコード

開発をローカルから Github に移したので, 最新のソースはこちらを見てください.

import sys
import math

DIVISION = 1000.0
CYCLE = 10

def angles(du, a, b):
    phi = 0
    while phi <= 2 * math.pi:
        yield phi
        phi += du / math.sqrt((a * math.sin(phi))**2 + (b * math.cos(phi))**2)

def coordinate(du, a, b):
    for angle in angles(du, a, b):
        yield (a * math.cos(angle), b * math.sin(angle))

def circumference(a, b):
    """Return a length of a circumference of an ellipse.

    @param a, b length of two semi-axes

    reference: http://en.wikipedia.org/wiki/Ellipse#Circumference
    """
    expr = ((a - b)/(a + b))**2
    return math.pi * (a + b) * (1 + (3 * expr)/(10 + math.sqrt(4 - 3 * expr)))

def draw(argv=None):
    if not argv:
        argv = sys.argv
    filename = argv[1]
    a = float(argv[2])
    b = float(argv[3])
    du = circumference(a, b) / DIVISION

    with open(filename, 'w') as out_file:
        print >> out_file, 'plot "-" w l'
        print >> out_file, '#x\ty'
        for i, coord in enumerate(coordinate(du, a, b)):
            i %= CYCLE * 2
            if i < CYCLE:
                print >> out_file, '{0}\t{1}'.format(*coord)
            else:
                print >> out_file, ''
        print >> out_file, 'end'
        print >> out_file, 'pause -1'


if __name__ == '__main__':
    draw(sys.argv)

コマンド

$ python ellipse.py ellipse.dat 1.0 0.5
gnuplot> plot "ellipse.dat" w l

出力結果

以下のようにちゃんと点線が等長、等間隔で描けています。(実際にプログラムで点線の長さと間隔を測ったところ、1 % 程度の誤差はありましたがだいたい同じでした。)

$a$ と $b$ には大小の制限は無いので縦長の楕円も描けます。

あとがき

昨日一昨日の記事を書いてから、色々な方からコメントやツッコミをいただきました。ありがとうございます。

今まで何かを作ったり、記事を書いたりして公開したことはあったけれど、今回もらった反応は今までで一番多いものでした。
そして、何かを発表して反応をもらう喜びを強く感じました。

これは病み付きになりそうです。

もしかして、オープンソースで活動してる人の中には、この喜びを原動力にやってる人もいるのかなぁ? と思ったり。

それでは、お休みなさい。

2011-01-19

点線で楕円を描く・その2

math

[追記] この記事には誤りが含まれています。訂正した記事はこちらです。

昨日の記事の続きです。

昨日は計算式までは出したのですが、実際に計算せずにはいられなかったのでやってみました。

Numpy とかあるのは知ってるのですが、標準の Python だけでどこまでできるかやってみました。
描画は gnuplot 使ってやってます。

ソースコード

ちゃちゃっと書いたので、コメントとかは省いてます。

import sys
import math

def delta(phi, du, k2):
    return phi + du / math.sqrt(1 - k2 * math.sin(phi)**2)

def draw(argv=None):
    if not argv:
        argv = sys.argv
    filename = argv[1]
    a = float(argv[2])
    b = float(argv[3])
    du = float(argv[4])

    k2 = 1 - (b/a)**2
    phi = 0

    with open(filename, 'w') as out_file:
        print >> out_file, '#x\ty'
        i = 0
        while phi <= 2 * math.pi:
            i += 1
            i %= 10
            if i < 5:
                print >> out_file, '{0}\t{1}'.format(a * math.cos(phi), b * math.sin(phi))
            else:
                print >> out_file, ''
            phi = delta(phi, du, k2)

                

if __name__ == '__main__':
    draw(sys.argv)

コマンド

$ python ellipse.py ellipse.dat 1.0 0.5 0.01

gnuplot> plot "ellipse.dat" w l

plot 結果

あとがき

思ったより誤差が発散しなかったので楽でした。

今日ばっかりは俺ってすげーって思いました :P

それでは。

2011-01-18

点線で楕円を描く

math

[追記] この記事には誤りが含まれています。訂正した記事はこちらです。

今日は 1 つネタをもらったので, ガチな数学をやろうと思います.

お題は「点線で楕円を描く」です.

円の場合は, 等角度ごとに区切って点線を書いていけば良いので簡単です.
たいていの場合, 三角関数は標準ライブラリにあるのであまり苦労しません.

しかし楕円の場合となると, 円がつぶれた形をしているので一筋縄ではいきません.

話の流れを追いつつ, どう解決していったか見ていきましょう.

発端

募集：PIL で点線の楕円を描画するのが得意な人。つーか、arc (円弧描画)関数だけで点線が書ける人。
http://twitter.com/#!/tk0miya/status/26986376482267136

楕円の場合は角度あたりの弧の長さが均等にならないよね。つまり点線を描画するにはどうすればいいんだ? 平面幾何なんてさっぱり分からないんだから、こんなん解決できないよ！
http://twitter.com/#!/tk0miya/status/26987619082571776

@tk0miya 楕円の弧長の話ですよね? その計算は楕円積分と言って, 数学で有名な積分計算です. 円もそうですが初等関数にならないので, きっちり計算する以上に楽な方法はありません. 計算式書きましょうか?
http://twitter.com/#!/cocoatomo/status/26991426927599616

実は, この楕円の弧長を求めるという問題は数学でも大きな問題でした.
いったいこれはどんな性質を持っているのか? 三角関数のようなものなのか?

数学のはなし

これを研究していた人たちにヤコビやワイエルシュトラスがいます.
彼らは
$\begin{align*}\begin{cases}x = a \cos \varphi \\ y = b \sin \varphi\end{cases}\end{align*}$
という楕円の周を $(a, 0)$ から左回りにたどったときの長さを求める積分
$\begin{align*}E(k,\quad \varphi) = \int_0^{\varphi} \sqrt{1 - k^2 \sin^2 \varphi} d \varphi \\ (k^2 = \frac{a^2 - b^2}{a^2}) \end{align*}$
の逆関数に着目し, その性質を調べ数学的な業績を挙げました.

その逆関数を $\varphi(k,\quad u)$ とするとその微分は以下のように求まります. (みなさん高校でやりましたね!?)
$\begin{align*}\frac{d \varphi}{du}(k,\quad u) &= \frac{1}{E'(k,\quad \varphi)} \\ &= \frac{1}{\sqrt{1 - k^2 \sin^2 \varphi}} \\ (u = E(k,\quad \varphi)) \end{align*}$

これを微分の形から差分の形に直して
$\begin{align*} \Delta \varphi &= \frac{\Delta u}{\sqrt{1 - k^2 \sin^2 \varphi}}\end{align*}$
これをさらに総和の形に直すと
$\begin{align*}\varphi_0 &= 0 \\ \varphi_n &= \varphi_{n-1} + \Delta \varphi_{n-1} \\ &= \varphi_{n-1} + \frac{\Delta u}{\sqrt{1 - k^2 \sin^2 \varphi_{n-1}}} \\ (\varphi_n \approx \varphi(n * \Delta u) \) \end{align*}$
となる.

なので, $\varphi (k,\quad u)$ の値を求めるときは $n$ を十分大きく (つまり $\Delta u = u/n$ 十分小さく) とれば $\varphi$ の値が差分の総和を計算することで求まり, $x = a \cos \varphi,\quad y = b \sin \varphi$ という計算で座標の値に戻せる.

誤差評価を行ってないので収束あたりがまだあやしいけれど, 一応目処は立った.
これからは誤差が出てしまったときの対処が残っているかな.

あとがき

こんなんでどうでしょう, tk0miya さん.

2011-01-17

グラフのはなし・その11

math ComputerScience order graph Python

さて前回 mro のトポロジカルソートを使った実装までできたのですが, mro の失敗が検知できない不十分なアルゴリズムで実装してしまいました.

今回はその不十分なところを修正するために, 今まで 2 つ (「未探索」と「探索済」) しかなかったノードの状態を 3 つ (「未探索」と「探索中」と「探索済」) に増やします.

ソースコード

改良したソースコードはこんなふうになります.

変わっている部分は PyClass オブジェクトの checked フィールドが color フィールドに変わり, その取る値が「真偽値」から 'WHITE', 'GRAY', 'BLACK' の 3 色になったところです.

それに合わせて, 探索の状態に応じて色を変えるようにし, 探索済をチェックするコードも色で判定するように変えました.
まずは実行してみてください.

# -*- coding: utf-8 -*-
class PyClass(object):
    def __init__(self, name=''):
        self.name = name
        self.color = 'WHITE'  # 最初はみんな白
        self.parents = []

    def __str__(self):
        return self.name

    def __repr__(self):
        return str(self)

# クラスの準備
D = PyClass('D')
C1 = PyClass('C1')
C2 = PyClass('C2')
B1 = PyClass('B1')
B2 = PyClass('B2')
A = PyClass('A')

# 矢印の準備
D.parents = [C2, C1]
C1.parents = [B2, B1]
C2.parents = [B2, B1]
B1.parents = [A]
B2.parents = [A]

# 探索完了の記録用リスト
searched = []

# クラスの継承関係を探索する
def dfs(klass):
    print('-> {0.name}'.format(klass))
    if klass.color == 'BLACK' or klass.color == 'GRAY': # 色のチェック
        print(u'{0.name} はチェック済み. 色: {0.color}'.format(klass))
        print('<- {0.name}'.format(klass))
        return

    klass.color = 'GRAY' # 探索中の色に塗る
    for parent in klass.parents:
        dfs(parent) # さらに継承関係を探索する

    print(u'{0.name} を黒に塗って引き返す'.format(klass)) # 祖先クラスの探索が終わった引き返す
    klass.color = 'BLACK' # 探索完了の色に塗る
    global searched
    searched.append(klass) # 探索が完了したら記録
    print('{1} を追加: {0}'.format(searched, klass))
    print('<- {0.name}'.format(klass))


if __name__ == '__main__':
    dfs(D)
    searched.reverse()
    print(searched)

きっとこんな結果になったはずです.

-> D
-> C2
-> B2
-> A
A を黒に塗って引き返す
A を追加: [A]
<- A
B2 を黒に塗って引き返す
B2 を追加: [A, B2]
<- B2
-> B1
-> A
A はチェック済み. 色: BLACK
<- A
B1 を黒に塗って引き返す
B1 を追加: [A, B2, B1]
<- B1
C2 を黒に塗って引き返す
C2 を追加: [A, B2, B1, C2]
<- C2
-> C1
-> B2
B2 はチェック済み. 色: BLACK
<- B2
-> B1
B1 はチェック済み. 色: BLACK
<- B1
C1 を黒に塗って引き返す
C1 を追加: [A, B2, B1, C2, C1]
<- C1
D を黒に塗って引き返す
D を追加: [A, B2, B1, C2, C1, D]
<- D
[D, C1, C2, B1, B2, A]

引き返すときに色を黒く塗っているのが分かりますね.

またごちゃごちゃするので出力はしていませんが, あるクラスの探索を始めるときにはそのクラスを灰色に塗っています.

あとがき

今日のところはここまでにします.

次回はこの色付けがどのように役立つのか, やや理論寄りな話をしたいと思います.
理論と言ってもできるだけ分かりやすく書くつもりなので, 紙とペン, 鉛筆を用意しておいてください.

それでは.

2011-01-16

グラフのはなし・その10

math ComputerScience order graph Python

さて今回は C3 アルゴリズムをトポロジカルソートで書けることを示してみましょう.

戦略

さて C3 アルゴリズムの肝は何だったかと言うと, それぞれの親クラスにとっての祖先クラスの順序と, 親クラスどうしの順序が食い違わないように, 祖先クラスを並べることでした. その操作を $merge$ 関数で行っていたのですね.

トポロジカルソートでもこれを真似してみましょう.
注意する点としては, このブログで出したトポロジカルソートのアルゴリズムでは, リストに点を追加していって最後に引っ繰り返していることです.

なので, あるクラスの親クラスリストでは順序を逆に並べます.

実装して実行してみる

こんな例を使いましょう. この継承関係は食い違いは起こしていませんね? 確認してみてください.

クラス: 親クラス
B1: A
B2: A
C1: B1, B2
C2: B1, B2
D: C1, C2

このような継承関係を Python スクリプトで書くとこうなります.

# クラスの準備
D = PyClass('D')
C1 = PyClass('C1')
C2 = PyClass('C2')
B1 = PyClass('B1')
B2 = PyClass('B2')
A = PyClass('A')

# 矢印の準備
D.parents = [C2, C1]
C1.parents = [B2, B1]
C2.parents = [B2, B1]
B1.parents = [A]
B2.parents = [A]

優先度とは逆に親クラスが並べてあるのが分かりますね.

さて先祖クラスの順序を求める以下の Python スクリプトを動かしてみましょう.
コメントや文字列が少し変更されてますが, スクリプトの中身はほとんど昨日使っていたものと同じです.
(Python2.7 で動くことを確認しています. Python2.6 や Python3.1 でも動くはずです.)

# -*- coding: utf-8 -*-
class PyClass(object):
    def __init__(self, name=''):
        self.name = name
        self.checked = False
        self.parents = []

    def __str__(self):
        return self.name

    def __repr__(self):
        return str(self)

# クラスの準備
D = PyClass('D')
C1 = PyClass('C1')
C2 = PyClass('C2')
B1 = PyClass('B1')
B2 = PyClass('B2')
A = PyClass('A')

# 矢印の準備
D.parents = [C2, C1]
C1.parents = [B2, B1]
C2.parents = [B2, B1]
B1.parents = [A]
B2.parents = [A]

# 探索完了の記録用リスト
searched = []

# クラスの継承関係を探索する
def dfs(klass):
    print('-> {0.name}'.format(klass))
    if klass.checked:
        print(u'{0.name} はチェック済み'.format(klass))
        print('<- {0.name}'.format(klass))
        return

    klass.checked = True
    for parent in klass.parents:
        dfs(parent) # さらに継承関係を探索する
    else:
        print(u'{0.name} から引き返す'.format(klass)) # 祖先クラスの探索が終わった引き返す

    global searched
    searched.append(klass) # 探索が完了したら記録
    print('{1} を追加: {0}'.format(searched, klass))
    print('<- {0.name}'.format(klass))


if __name__ == '__main__':
    dfs(D)
    searched.reverse()
    print(searched)

さて, みなさんの環境で動かすことはできたでしょうか?
スクリプトの動きについてはコメントや出力メッセージでなんとなく分かるように作ってあります. (分かりづらければコメントくださいな.)

ちゃんと動けばこんな出力になっているはずです.

-> D
-> C2
-> B2
-> A
A から引き返す
A を追加: [A]
<- A
B2 から引き返す
B2 を追加: [A, B2]
<- B2
-> B1
-> A
A はチェック済み
<- A
B1 から引き返す
B1 を追加: [A, B2, B1]
<- B1
C2 から引き返す
C2 を追加: [A, B2, B1, C2]
<- C2
-> C1
-> B2
B2 はチェック済み
<- B2
-> B1
B1 はチェック済み
<- B1
C1 から引き返す
C1 を追加: [A, B2, B1, C2, C1]
<- C1
D から引き返す
D を追加: [A, B2, B1, C2, C1, D]
<- D
[D, C1, C2, B1, B2, A]

探索が完了すると無事に正しい先祖クラスの順序になっていますね.

問題点

さてこのスクリプトには 1 つ問題があります. それは何でしょうか?

そうです. 継承関係の食い違いがあっても実行できて, 何らかの結果が出てきてしまうことです.
これはどこに問題があったのでしょうか?

PyClass というクラスのオブジェクトには checked というフィールドがあり, ここの値に「あるクラスが既に探索されたか?」という情報 (真偽値) を保存しておき二重探索を防いでいました.

しかし実は, 継承関係の食い違いを見付け出すためには真と偽の 2 つの値では不足していて, 「未探索」「探索中」「探索済」の 3 種類の値が無いといけないのです.

深さ優先探索のアルゴリズムから修正しないといけないので, また明日以降に今のスクリプトを改造していきましょう.

あとがき

なんだか最初の順序のはなしからだいぶ離れてしまったように感じますか?
これでも当初の目論見どおりなのです.

ここ最近, グラフやグラフアルゴリズムの勉強をしていたのですが, そこでふと「グラフアルゴリズムって, 点の整列を行っている」ことに気付いたのです.

「もしかして,『順序』という視点でグラフの話ができないか? それに最愛の言語 Python の mro (method resolution order, メソッド解決順序) もグラフアルゴリズムだ. そして Hadoop の MapReduce のジョブフローもグラフと見れ, その上のアルゴリズムは重要になるだろう. よし単純な順序の話を切り口に, 主にグラフの話しをしよう.」

こう思ってこの連載はスタートしました.
1ヶ月ネタが続くか不安でしたが, なんとか折り返し地点まで来れました.
これも毎日読んでいてくださるみなさまのおかげです.

後半も気を引き締めて, 回を落とさないように頑張ります.
よろしくです.

追伸

グラフアルゴリズムを解説するときは, グラフィカルな表示も欲しいですよね.
自分も欲しいです.

そう思ったので現在鋭意製作中です.
今は PyOpenGL で作っています.

完成をお楽しみに.

2011-01-15

グラフのはなし・その9

math ComputerScience order graph Python

今日は Perl や Python のクラスの多重継承で使われているアルゴリズム C3 を紹介します.

多重継承でメソッドを呼び出した場合, 複数の親クラスを持つので今呼び出しているメソッドがどの親に属しているのかを調べなければなりません.
そしてその順序も決めておき, 必ず決まったメソッドが呼ばれるようにしなくてはなりません.

それを実現する C3 アルゴリズムとはどんなものなのでしょうか?

C3 アルゴリズム概要

このページを参考にして解説を行います.
http://www.python.org/download/releases/2.3/mro/

数式で書いてしまえば,

$\mathrm{mro}(C) = [C] + \mathrm{merge}(\mathrm{mro}(P_1),\quad \mathrm{mro}(P_2),\quad \dots,\quad \mathrm{mro}(P_n),\quad [P_1,\quad P_2,\dots,\quad P_n])$

([P1, P2,..., Pn] は C の親クラス)

となります.
[C] は C 1 つからなるリストで, その後ろの "+" はリストの結合です.

右辺で再度 mro 関数が呼び出されていますが, 親クラスを持たないクラス C に対して $\mathrm{mro}(C) = C$ と定義されるので, この計算はどこかで止まります.

merge 関数では引数に入ってきたリストの順序を崩さずに, 新たなリストを作っているので, まさにトポロジカルソートが行われています.

C3 アルゴリズムの動き

具体的な merge 関数の動きを見ていきましょう.

まず, 先頭の引数 $\mathrm{mro}(P_1)$ (実体はリスト) から先頭の元 $P$ を取ってきます.
$P$ が他の引数 $\mathrm{mro}(P_2),\quad\mathrm{mro}(P_3),\dots,\quad\mathrm{mro}(P_n),\quad[P_1,\quad P_2,\dots,\quad P_n]$ のそれぞれの先頭にあるか, それとも含まれないか, をチェックします.
そのチェックが通った場合は $P$ を括り出します. (例はまた後で)
もしそのチェックに通らなければ, $\mathrm{mro}(P_2)$ に対して同じチェックを行います.
全ての引数でチェックが通らなければ $\mathrm{mro}$ 失敗です.

例

言葉だけの説明では分かりづらいので, 具体例を扱っていきましょう. (何度も言ってますが, 具体例を考えるのはすごく重要です.)

クラス: 親クラス
C: B

単純な継承です. すごく簡単だし, わざわざ計算しなくてもすぐに分かりますが, 練習として解いてみましょう.

$\begin{align*}\mathrm{mro}(C) &= [C] + \mathrm{merge}(\mathrm{mro}(B),\quad [B]) \\ &= [C] + \mathrm{merge}([B],\quad [B]) \\ &= [C] + [B] \\ &= [C,\quad B]\end{align*}$

確かに望んだ結果が出てきましたね.

では次はもう少し出てくるクラスを増やしましょう.

クラス: 親クラス
C: B1, B2

単純な多重継承です.

$\begin{align*}\mathrm{mro}(C) &= [C] + \mathrm{merge}(\mathrm{mro}(B_1),\quad \mathrm{mro}(B_2) [B_1,\quad B_2]) \\ &= [C] + \mathrm{merge}([B_1],\quad [B_2],\quad [B_1,\quad B_2]) \\ &= [C] + [B_1] + \mathrm{merge}([],\quad [B_2],\quad [B_2]) \\ &= [C] + [B_1] + [B_2] \\ &= [C] + [B_1,\quad B_2] \\ &= [C,\quad B_1,\quad B_2]\end{align*}$

こんな継承関係を考えてみましょう. これは有名なダイヤモンド継承ですね.

クラス: 親クラス
C: B1, B2
B1: A
B2: A

$\begin{align*}\mathrm{mro}(C) &= [C] + \mathrm{merge}(\mathrm{mro}(B_1),\quad\mathrm{mro}(B_2),\quad[B_1,\quad B_2]) \\ &= [C] + \mathrm{merge}([B_1,\quad A],\quad [B_2,\quad A] ,\quad[B_1,\quad B_2])\\ &= [C] + [B_1] + \mathrm{merge}([A],\quad [B_2,\quad A],\quad [B_2] ) \\ &= [C] + [B_1] + [B_2] + \mathrm{merge}([A],\quad [A],\quad []) \\ &= [C] + [B_1] + [B_2] + [A] \\ &= [C,\quad B_1,\quad B_2,\quad A] \end{align}$

今度はこんな計算になり, 上手く多重継承の先祖クラスの順序を扱えていますね.

4 段目で $A$ ではなく $B_2$ が括り出されているのは, $\mathrm{merge}$ 関数の第 2 引数に $[B_2,\quad A]$ というリストがあり, この順序を守るためです.

では最後に mro に失敗する例を考えてみましょう.

クラス: 親クラス
B1: A
B2: A
C1: B1, B2
C2: B2, B1
D: C1, C2

この継承関係では $C_1$ の親クラスの順序と $C_2$ の親クラスの順序が食い違っているので, きちんと先祖クラスの順序が付けられません.

実際に Python でクラス D を作ろうとすると, 以下のようにエラーが出て作ることができません.

>>> A = object
>>> class B1(A): pass
... 
>>> class B2(A): pass
... 
>>> class C1(B1, B2): pass
... 
>>> class C2(B2, B1): pass
... 
>>> class D(C1, C2): pass
... 
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
TypeError: Error when calling the metaclass bases
    Cannot create a consistent method resolution
order (MRO) for bases B2, B1

さて $\mathrm{mro}$ で計算してみましょう.

$\begin{align*}\mathrm{mro}(D) &= [D] + \mathrm{merge}(\mathrm{mro}(C_1),\quad \mathrm{mro}(C_2),\quad [C_1,\quad C_2]) \\ &= [D] + \mathrm{merge}([C_1,\quad B_1,\quad B_2,\quad A],\quad [C_2,\quad B_2,\quad B_1,\quad A],\quad [C_1,\quad C_2]) \\ &= [D] + [C_1,\quad C_2] + \mathrm{merge}([B_1,\quad B_2,\quad A],\quad [B_2,\quad B_1,\quad A],\quad [])\end{align*}$

とここまできて困ってしまいました. B1 と B2 の順序をどちらを先としてもうまく $\mathrm{merge}$ できません.
なので $\mathrm{mro}$ の計算は失敗し, ここで終了です.

あとがき

今日はトポロジカルソートが使われている一例として, Python の mro を扱いました.
昨日のトポロジカルソートとのつながりが話せていないので, 次回は Python の mro を実装し, 今までのトポロジカルソートと比べてみます.

2011-01-14

グラフのはなし・その8

math ComputerScience order graph Python

さて今日は目標の1つであったトポロジカルソートのプログラムを見てみましょう.
基本的な構造は昨日の深さ優先探索 (DFS) と同じです.

今日トポロジカルソートを行う DAG (輪っか無しの矢印グラフ) はこれです.

0 -> 1 -> 3
 |     `-> 4
 |         V
  `-> 2 -> 5
       `-> 6

昨日のグラフとだいたい同じですが,「4 → 5」の矢印が増えたところが違いますね. ここは後でポイントになりますので, 覚えておいてください.

トポロジカルソートを行う Python スクリプト

# -*- coding: utf-8 -*-                                                                                                                                         
class Node(object):
    def __init__(self, value=0):
        self.value = value
        self.checked = False
        self.children = []

    def __str__(self):
        return 'node{0}'.format(self.value)

    def __repr__(self):
        return str(self)

# 点の準備                                                                                                                                                      
node0 = Node(0)
node1 = Node(1)
node2 = Node(2)
node3 = Node(3)
node4 = Node(4)
node5 = Node(5)
node6 = Node(6)

# 矢印の準備                                                                                                                                                    
node0.children = [node1, node2]
node1.children = [node3, node4]
node2.children = [node5, node6]
node4.children = [node5]

# 探索完了の記録用リスト
searched = []

# 分かれ道にきたらこの手順を踏む
def dfs(node):
    print('-> {0.value}'.format(node))
    if node.checked:
	print(u'チェック済み')
	print('<- {0.value}'.format(node))
	return

    node.checked = True
    for child in node.children:
	dfs(child) # さらに分かれ道にきたら同じ手順を繰り返す
    else:
        print(u'引き返す')

    global searched
    searched.append(node) # 探索が完了したら記録
    print(searched)
    print('<- {0.value}'.format(node))


if __name__ == '__main__':
    dfs(node0)
    searched.reverse()
    print(searched)

昨日のソースコードと違っているところはどこでしょうか?

searched というリストが増えた.
dfs というメソッドの終わりの方に searched に node を追加する操作が増えた.

この 2 点だけです.

実行結果は以下のようになりました. 最後の行がトポロジカルソートを行った結果です.
昨日の結果と比べて違っているところはどこでしょうか??

-> 0
-> 1
-> 3
引き返す
[node3]
<- 3
-> 4
-> 5
引き返す
[node3, node5]
<- 5
引き返す
[node3, node5, node4]
<- 4
引き返す
[node3, node5, node4, node1]
<- 1
-> 2
-> 5
チェック済み
<- 5
-> 6
引き返す
[node3, node5, node4, node1, node6]
<- 6
引き返す
[node3, node5, node4, node1, node6, node2]
<- 2
引き返す
[node3, node5, node4, node1, node6, node2, node0]
<- 0
[node0, node2, node6, node1, node4, node5, node3]

もちろんアルゴリズムは変更していないので, 昨日のグラフと今日のグラフの違うところが影響しているわけですね.

今日のグラフは「4 → 5」の矢印が増えているので, 4 の場所まできたときにさらに 5 の場所まで探索できます. それが出力の↓この部分に出ています.

-> 4
-> 5
引き返す
[node3, node5]
<- 5

考察

さて, 昨日のソースコードには無かった処理「searched.append(node) # 探索が完了したら記録」の部分について考えてみましょう.

当然この部分がトポロジカルソートの肝なのですが,「探索が完了」するとはどういうことでしょうか?
深さ優先探索とは,「探索していないところがあったらどんどん先に進んで, 行き止まりに到達したら引き返す」というものでしたね.
つまりここで言っている「探索が完了」というのは,「ここから行けるところはもう全て探索したよ」という意味なのです.

なので, searched (探索完了) という名前のリスト中のある点 A を見たとき, その点 A から到達できる点は, 必ずそれより前の方にあることになります.
ここがポイントなのですが分かりましたか?
自分でグラフを書いて, 上の処理を 1 つ 1 つ追っていってみてください.

こういうときに重要なのは具体例に対して自分の手を動かすこと.
頭だけで考えただけだと大抵はあやふやな理解しかできていません.
はっきり言いましょう. 手を動かさないで理解した気になっているのは勘違いです!! 手を動かして実験してみてください.

最後に reverse では結果が左右逆になっていて見づらいので, 引っくり返しているのです.

あとがき

昨日の深さ優先探索だったスクリプトが突然トポロジカルソートのスクリプトに変わって驚いたでしょうか?
実は, 私も深さ優先探索とトポロジカルソートのこの関係を知って驚きました.
しかし, よくよく考えて頭で味わってみると, なるほどよくできている関係でトポロジカルソートの性質をすっきり言い切っているなぁ, と感心しました.

こういう 1 つの概念やイメージが頭の中にスッと入ってくる瞬間が楽しくて私は数学をやっています.
みなさんもこの連載でそのような瞬間が味わえると望外の幸せです.

それでは.