GNU Smalltalk チュートリアル（その4）

gst Smalltalk

続きです。

前回までのあらすじ

前回は、N-Queen パズルの実行時間を計測してみました。その結果

  • 風子さんはかわいいコンピュータって案外遅い

という結論に落ち着いたところでの今回です。（えーっ）

サイを投げてみよう

カモノハシ本の8-Queen のアルゴリズムでは、わたしの 風子さん(Atom N270 なマシン)で、28-Queen で 47分、29-Queen で26分、30-Queenで 1359分(22時間)！さすがにこれは待てません。これだったらサイコロを振ってみたほうがマシかも・・いや、冗談ではなく。

そんなわけで、本当にサイコロ任せでいってみましょう。

作戦はこんなです。まず、縦・横に被らないように Queen を配置します。そしてQueenを一組選んで交換していけば、いつかは、正解にたどり着くはず。絵にするとこんな感じ。

で、本当に適当に交換していったのでは、うまく正解にたどり着くにはよっぽどの運が必要になっちゃうので、ランダムに得た解を今の解と比べて良かったら採用！というアルゴリズムでだんだんと正解までたどり着くようしてみます。

では早速コーディング。まずは問題自体の定義。

  • nqueenProblem.st
Object subclass: NQueenProblem [
    | queenCnt rand |

    NQueenProblem class >> new: queenCnt [
        ^super new initialize: queenCnt ; yourself
    ]

    initialize: num  [
        queenCnt := num
    ]

    getInitialState [
        "初期配置を取得"
        | s0 |
        s0 := OrderedCollection new.
        1 to: queenCnt do: [:x|s0 add: x ].
        ^s0
    ]

    cost: s [
        | set cost |
        cost := 0.

        "/"
        set := Set new.
        1 to: s size do: [:row |
            | col |
            col := s at: row.
            (set includes: col + row) ifTrue: [cost := cost + 1].
            set add: col + row].

        "\"
        set := Set new.
        1 to: s size do: [:row |
            | col |
            col := s at: row.
            (set includes: col - row) ifTrue: [cost := cost + 1].
            set add: col - row].

        ^ cost
    ]

    neighbor: s [
        | row1 row2 col1 col2 |
        rand ifNil: [rand := Random new].

        row1 := rand between:1 and: s size.
        row2 := rand between:1 and: s size.
        col1 := s at:row1.
        col2 := s at:row2.
        ^s copy at: row1 put: col2;
                at: row2 put: col1;
                yourself.
    ]

    isEndCondition: s cost: cost [
        ^cost = 0
    ]
]

ではではコードを説明。getInitStatus で、縦横方向は絶対被らないけれど 斜め方向のことは知ったことじゃない初期配置を作ります。

neighbor: s で、現在の状態 s からランダムに選んだクイーン一組だけ交換した近傍状態を作ります。ちなみにSmalltalk の Random クラスは、Stream の一種になっていて、next を呼ぶことで次の乱数を返してくれます。もちろん値は普通に 0〜1の間ですが、between:and: を使えば その間の値の乱数を返してくれるので便利です。

cost: s は 今の状態がどれだけ正解に近いかを算出します。具体的には斜め方向に盤面を舐めて、どれだけのQueen が取られてしまうかを数え上げた値、要するに「角行」方向チェックです。「飛車」方向チェックについては、getInitStatus と cost: で絶対縦横方向に被らない仕組みになってるので、オミットしました。

最後、isEndCondition:cost: は、終了条件チェックなのですが、もちろん cost が 0 (Queen がどれも取られない解）が見つかったらおしまいです。


次は正解への近づき方。

  • randomSolver.st
Object subclass: RandomSolver [
    | problem curS curCost bestS bestCost time |

    setProblem: aProblem [
        problem := aProblem
    ]

    problem [
        ^problem
    ]

    time [
        ^time
    ]

    curS [
        ^curS
    ]

    curCost [
        ^curCost
    ]


    solveProblem: s0 maxTime: maxTime [
        "
         @param s0  初期状態
         @param maxTime 最大時間
         "
        curS  := s0.
        bestS := curS.

        curCost  := problem cost: curS.
        bestCost := problem cost: bestS.
        time  := 0.

        [time < maxTime] whileTrue: [
            self tryAnswer.
            (problem isEndCondition: curS cost: curCost) ifTrue: [^curS].
            time := time +1
        ].
        ^curS
    ]

    tryAnswer [
        | newS newCost dCost |

        newS    := problem neighbor: curS.
        newCost := problem cost: newS.
        dCost   := newCost - curCost.

        (dCost < 0) ifTrue: [self updateS: newS cost: newCost]
    ]


    updateS: newS cost: newCost [
        curS    := newS.
        curCost := problem cost: curS.

        (newCost < bestCost) ifTrue: [
            bestS    := newS.
            bestCost := problem cost: bestS ].

        self changed
    ]
]

といってもぐるぐるループを回して、ランダムに得られた次の状態が今の状態よりよかったら、解を上書きしていく・・という単純なものです。問題が解けなくても無限ループにならないように maxTime で最大ループ回数を制限しています。

で、これらを使って問題を解いてみます。

#!/usr/bin/env gst

Object subclass: SolutionPrinter [
    update: solver [
        self print: solver
    ]

    print: solver [
        Transcript 
            show: '(';
            show: (solver curS collect:[:x| x asString, ' ']) join;
            show: ')';
            show: '(time:';
            show: solver time printString;
            show: ' cost:';
            show: solver curCost printString;
            show: ')';
            cr.
    ]
]

Eval [
    | problem s0 sa ret printer |

    FileStream fileIn: 'nqueenProblem.st'.
    FileStream fileIn: 'randomSolver.st'.
    FileStream fileIn: 'nqueenResultPrinter.st'.

    problem := NQueenProblem new: 8.
    s0      := problem getInitialState.

    solver := RandomSolver new.
    solver setProblem: problem.

    printer := SolutionPrinter new.
    solver addDependent: printer.

    ret := solver solveProblem: s0 maxTime: 3000.

    '=========================' printNl.
    printer print: solver.
    NQueenResultPrinter new print: ret
]

早速 8-Queen で試してみました。

んー。確かに解けるときには解けるのですが（そして早い）、結構な頻度で解けなかったりもします。まさに運まかせです。さすがにサイコロを振っているのでその辺は仕方がないというところでしょうか(^^;

ところがこの解けなかった方の答えをよぉーく見てみると、どの Queen を互いに交換しても、つまりこの先どーんなに Queen を交換しつづけようと正解にたどり着けない状態になっています。ちょうどカモノハシ本の解き方で、ある列のQueenが全部の行に進んでみてもどれもNGだった様なケース(隣の部分解を崩してやり直しするケース）みたいな状態。

$ gst run.st 
OrderedCollection (1 2 3 4 5 6 7 8 )
(1 2 3 4 6 5 7 8 )(time:2 cost:6)
(1 6 3 4 2 5 7 8 )(time:3 cost:5)
(1 6 3 4 7 5 2 8 )(time:5 cost:4)
(1 8 3 4 7 5 2 6 )(time:10 cost:2)
(4 8 3 1 7 5 2 6 )(time:65 cost:1)
"Global garbage collection... done"
'========================='
(4 8 3 1 7 5 2 6 )(time:3000 cost:1)
+---+---+---+---+---+---+---+---+
|   |   |   | Q |   |   |   |   |
+---+---+---+---+---+---+---+---+
|   |   |   |   |   |   | Q |   |
+---+---+---+---+---+---+---+---+
|   |   | Q |   |   |   |   |   |
+---+---+---+---+---+---+---+---+
| Q |   |   |   |   |   |   |   |
+---+---+---+---+---+---+---+---+
|   |   |   |   |   | Q |   |   |
+---+---+---+---+---+---+---+---+
|   |   |   |   |   |   |   | Q |
+---+---+---+---+---+---+---+---+
|   |   |   |   | Q |   |   |   |
+---+---+---+---+---+---+---+---+
|   | Q |   |   |   |   |   |   |
+---+---+---+---+---+---+---+---+

・・サイコロが振り足りないわけじゃないのね。

なるほど、これが有名な局所的最適解(local minimum)にはまり込んだというやつですねっ！！

・・・・白々しくってすみません。

Simulated Annealing

この手の組合せ最適化問題の解き方は色々あります。とっても名前がメジャーな「遺伝的アルゴリズム(Genetic Algorithm)」はじめ、「タブーサーチ(Tabu Search)」、「蟻コロニー最適化(Ant Colony Optimization)」、必殺技のようなカッコいい「シミュレーティド・エボリューション(Simulated Evolution)」！などなど。

そんな、名前だけでもオラ何だかワクワクしてしきてしまうヤツらの中から、今回はあまり普通の人の間で有名でない「シミュレーティド・アニーリング(Simulated Annealing)」 をチョイスします。


* * *


Simulated Annealing は、焼き鈍し(Annealing)を模したアルゴリズムです。

焼き鈍しとは、高温に熱した金属をゆっくりと冷やすことて鈍すこと。急に温度を下げると綺麗な結晶にならないですが、こんなふうにゆっくり冷やすと秩序のあるきれいな結晶構造、つまり系の最小エネルギー状態に戻るわけで、例えば加工硬化しちゃった金属なんかを再び柔らかくしたりします。

Simurated Annealing の特徴は、現状より悪化する状態にも時として遷移すること。これにより先ほどの局所的な最適解の谷間にはまり込んでもそこから抜け出すことが出きるようになります。

ただ、いつでも単純にサイコロ次第で現状より悪化してもOK では、なかなか解にたどり着けません。Simulated Annealing は焼き鈍しを模すこと、つまり温度が高いほどコストを無視した状態により遷移しやすく、冷えてきたらあまりコストに逆行する状態に遷移しにくくなることで、単なるランダムよりも効率よく、もっともコストが低い状態（答え）にたどり着きやすくなるわけです。

・・・と、教科書丸写しの薀蓄は置いといて、早速コードを書いてみましょう。Simulated Annealing は、その特徴としてアルゴリズムが簡単で、プログラムにするのが楽ちんであることも挙げられていて、そんな意味でもワナビーなわたしにはまさにぴったりです。

問題側 NQueenProblem クラスはそのままに、先ほどの RandomSolver だけ新しいクラス SA に置き換えます。（ハイ、最初からそのつもりでした）

  • sa.st
Object subclass: SA [
    | problem curS curCost bestS bestCost time temp rand |

    setProblem: aProblem [
        problem := aProblem
    ]

    problem [
        ^problem
    ]

    time [
        ^time
    ]

    curS [
        ^curS
    ]

    curCost [
        ^curCost
    ]

    temp [
        ^temp
    ]

    simulatedAnnealing: s0 t: t0 alpha: alpha beta: beta maxTime: maxTime [
        "
         @param s0  初期状態
         @param t0  初期温度
         @param alpha 冷却率 (0 < alpha < 1)
         @param beta  アニーリングに費やされる時間の増加率 (bata >= 1)
         @param maxTime 最大時間

         m 次パラメータ更新までの時間
         "
        | m |
        temp  := t0.
        curS  := s0.
        bestS := curS.
        m := 10.

        rand := Random new.

        curCost  := problem cost: curS.
        bestCost := problem cost: bestS.
        time  := 0.

        [time < maxTime] whileTrue: [
            self metropolis: temp m: m asInteger.

            time := time  + m asInteger.
            temp := alpha * temp.
            m    := beta  * m.

            (problem isEndCondition: curS cost: curCost) ifTrue: [^curS]
        ].
        ^curS
    ]

    metropolis:temp m: m [
        | newS newCost dCost |
        m timesRepeat: [
            newS    := problem neighbor: curS.
            newCost := problem cost: newS.
            dCost   := newCost - curCost.

            (dCost < 0)
                ifTrue: [self updateS: newS cost: newCost]
                ifFalse: [
                    (rand next < (Float e raisedTo: -1 * (dCost / temp)))
                        ifTrue: [self updateS: newS cost: newCost]
                ]
        ]
    ]

    updateS: newS cost: newCost [
        curS    := newS.
        curCost := problem cost: curS.

        (newCost < bestCost) ifTrue: [
            bestS    := newS.
            bestCost := problem cost: bestS ].

        self changed
    ]
]

こんな感じ。これもまた教科書のコードそのまんま(^^;

組合せ最適化アルゴリズムの最新手法―基礎から工学応用まで

組合せ最適化アルゴリズムの最新手法―基礎から工学応用まで

(↑今回の教科書）

プログラム化が平易といわれるだけあって、パッとみればすぐ分かる簡単なコードになっていますが、一応分かりにくそうなところを解説。なんだか名前から何者かが分からない alpha(=α)とbate(=β)はアニーリングの冷却具合の調整用。一回の温度下降あたりどれくらいの割合で低下するかの係数がαで、同一温度で何回試行するかの伸び率がβ。Simulated Annealing の効率は冷却パラメータでずいぶんとかわるので、重要です。

悪化する解を選ぶ確率は P(RANDOM < e^-Δcost/T) で与えられる・・とあるので、そのまんまコーディング。てっ、そう言えば Smalltalk のべき乗って、#raisedTo: が基本なのですが、VisualWorks では #** も使えたりします*1。で、GNU Smalltalk ではどうなのかというと・・・、残念ながら #** は使えません。せっかくなので定義してみましょう。

余談: GNU Smalltalk で既存クラスにメソッドを足す

Smalltalk では出来合いの既存クラスに メソッドを足すなんてことは日常茶飯事なのですが、システムブラウザから直接変える 普通のSmalltalk とちがって、ファイルにスクリプトを書く GNU Smalltalk ではどう書くばよいのかな？・・と調べてみました。

Number extend [
    ** num [
        ^self raisedTo: num
    ]
]

意外に簡単。これで先の SA >> metropolis:m: メソッドの

(rand next < (Float e raisedTo: -1 * (dCost / temp)))

をすっきりわかりやすく

(rand next < (Float e ** (-1 * (dCost / temp))))

風に書き換えられます。（って、2項メッセージになっちゃったのでカッコが一個増えちゃって、反って読みにくくなったかも(^^; )

Simulated Annealing すごい！

というわけで、早速 これを使って解いてみます。

  • run.st
#!/usr/bin/env gst

Object subclass: SAPrinter [
    update: sa [
        self print: sa
    ]

    print: sa [
        Transcript 
            show: '(';
            show: (sa curS collect:[:x| x asString, ' ']) join;
            show: ')';
            show: '(time:';
            show: sa time printString;
            show: ' temp:'; 
            show: sa temp printString; 
            show: ' cost:';
            show: sa curCost printString;
            show: ')';
            cr.
    ]
]


Eval [
    | problem s0 sa ret printer |

    FileStream fileIn: 'nqueenProblem.st'.
    FileStream fileIn: 'sa.st'.
    FileStream fileIn: 'nqueenResultPrinter.st'.

    problem := NQueenProblem new: 50.
    s0      := problem getInitialState.
    s0 printNl.

    sa := SA new.
    sa setProblem: problem.

    printer := SAPrinter new.
    sa addDependent: printer.

    ret := sa simulatedAnnealing: s0 
              t: 2 alpha: 0.9 beta: 1.0 maxTime: 3000.

    '=========================' printNl.
    printer print: sa.
    NQueenResultPrinter new print: ret
]

すごい！この Simulated Annealing 凄いよ！さすがターンAのお兄さん！！本当にサクサク解けます。Queen の 20個や30個は何のその。50-Queenだってホラ、

この通り。まさに人類の英知の結晶ですね！

って、「早い早い」って口で言うだけじゃわかりませんね。とりあえず前回風子さんが苦戦した30-Queen で計測してみます。

実行時間 (msec) 解 time (step) コスト
1回目 4425 17 6 21 9 30 4 27 15 12 29 26 20 2 24 10 13 28 18 3 8 23 5 19 25 1 7 14 22 11 16 2130 0
2回目 4684 17 19 24 11 7 25 18 8 15 1 6 27 22 26 5 3 16 4 13 30 9 23 28 2 21 29 10 20 14 12 2210 0
3回目 6343 13 15 19 5 23 4 12 27 22 6 25 20 7 16 26 2 9 21 8 29 3 1 30 14 18 11 24 28 17 10 3000 1
4回目 1291 14 19 23 25 27 4 30 12 9 24 3 21 10 29 1 6 2 20 8 11 15 18 7 17 26 13 22 16 28 5 590 0
5回目 3015 3 12 27 15 19 29 11 5 2 30 17 4 20 22 28 6 8 14 18 7 24 16 23 10 26 21 25 1 13 9 1350 0
6回目 3966 15 17 19 7 28 25 11 9 5 3 16 18 2 24 13 23 4 27 30 6 21 12 1 26 22 10 8 20 14 29 1910 0
7回目 2842 20 17 11 18 8 12 4 29 25 15 22 16 13 1 28 26 10 2 5 27 30 23 21 14 6 3 9 7 24 19 1390 0
8回目 1037 21 18 6 3 11 19 30 9 5 22 20 1 27 16 23 8 17 28 7 24 14 29 10 15 4 2 25 13 26 12 460 0
9回目 1569 23 11 30 10 16 20 5 7 27 15 26 19 6 4 2 13 3 22 29 8 21 17 12 25 28 1 18 24 9 14 730 0
10回目 1946 7 22 12 18 28 8 2 13 17 9 29 25 30 21 3 26 14 16 4 24 15 6 23 1 5 27 20 10 19 11 910 0

おぉっ！。たしかにたまに失敗はするけれど、この頻度でこれくらいの時間なら「失敗したらもう一度」で全然OKですね。（てかちゃんと冷却スケジュールを調整すればこんなに失敗しないのだけれども）

ちなみに 50-Queen だとこんな感じ。

実行時間 (msec) 解 time (step) コスト
1回目 6296 35 25 20 46 16 34 47 11 2 23 42 19 31 12 3 6 32 44 33 9 50 17 1 18 41 8 48 15 30 7 27 29 43 40 37 10 49 22 39 45 24 28 4 13 36 38 26 21 14 5 1520 0
2回目 12172 7 28 16 24 2 6 34 19 30 43 40 36 47 17 10 38 9 26 50 1 15 41 14 11 8 49 23 44 3 32 48 20 37 39 21 29 12 4 18 5 42 13 25 22 35 31 27 46 33 45 3000 1
3回目 11254 27 17 37 29 21 8 20 49 42 7 14 26 36 38 15 44 3 47 50 25 22 39 43 13 9 5 48 19 33 2 24 10 41 1 4 31 11 32 46 16 23 34 28 6 18 12 45 35 30 40 2760 0
4回目 12121 43 37 34 23 29 27 13 33 19 4 2 26 10 18 49 46 1 11 5 25 28 48 45 24 15 39 16 30 38 6 47 3 21 35 12 17 36 50 7 20 44 40 8 22 41 31 14 9 32 42 2900 0
5回目 5354 18 27 24 34 46 11 16 50 33 4 48 35 32 6 41 47 14 7 29 1 5 38 43 22 10 44 28 21 3 20 9 40 36 30 26 8 17 42 39 13 2 25 19 31 12 45 49 23 37 15 1290 0
6回目 12239 5 27 22 36 28 4 19 45 18 34 21 11 35 50 37 44 24 33 1 15 47 6 16 10 25 46 48 41 2 32 7 29 14 42 38 30 12 17 40 23 26 3 31 49 9 13 43 8 39 20 3000 1
7回目 11990 47 33 48 5 41 39 11 16 25 17 31 27 8 42 15 12 43 20 18 45 40 7 21 46 37 49 6 38 2 44 1 26 36 9 24 28 34 3 50 14 32 35 10 13 23 4 29 19 30 22 2940 0
8回目 9556 42 7 22 31 18 9 30 36 39 16 6 34 27 8 15 50 24 2 44 32 10 48 4 23 46 43 35 1 25 12 14 41 20 38 26 3 47 49 11 28 17 40 13 37 19 21 33 45 29 5 2260 0
9回目 3529 22 24 11 14 3 31 50 42 40 33 18 21 25 47 15 48 30 37 43 34 1 49 28 44 7 41 8 16 23 29 9 2 12 5 46 20 4 10 26 6 32 27 17 19 13 39 36 38 35 45 830 0
10回目 7804 29 50 25 16 22 6 27 17 28 41 32 19 11 49 38 3 18 20 37 30 24 33 47 40 7 10 42 1 15 36 2 8 14 23 48 35 45 21 34 9 5 46 31 4 39 43 13 44 26 12 1880 0

コスト計算とかが効いてくるので、Queen が増えれば増えるほど 重くなっていくので、ちょっと遅くなっています。


* * *


というわけで、全4回かけて GNU Smalltalk で N-Queen パズルをSimulated Annealing で解くところまでをお送りしました、が、・・・・・・・一体誰に向けて何のチュートリアルを書いているんでしょうかね、わたしは...orz

ただ、こう言うお題って、プログラミングしていてとても愉しいことだと思うんです。

以前 月が地球に落ちるお話で、Etoys で引力シミュレータを書いたのですが、引力もなんで公転するのかも、話で聞いて何となく分かった気になるのと、プログラムで実際に動かして理解するのは、本当に雲泥の差です。たとえばケプラーの第二法則があるだけで、簡単に惑星が恒星のまわりを安定して回るようになるところとか、本当にアハッとする経験です。

ちゃんと考えて論理的に筋道通した解き方よりも、盲目的にサイコロを降りつつ解を求めた方が、よっぽど早く、確実に答えが出ちゃうのは、たいていの直感に反します（奇跡や偶然という概念から）。でもこうしてやってみれば、あまりの効率のよさになるほど納得してしまうのです。

うん、やっぱり教育にプログラミングは必要ですよ！


* * *


そんなわけで、このネタを思いついたときには「シンプルでかつプログラミングの面白みを感じるネタだし、すごくいいんじゃない！」「コンピュータって意外に遅いんだよ、と言うのが実感できるのも愉しいし！」「Simulated Annealing は小学生にもかけるし、焼き鈍しの模倣が効率の良いアルゴリズムになるってあたりも自然科学チックで楽しそうだし！(Etoysでやって見るのも面白そう)」みたいに、私的にすごく盛り上がっていたのですが、出来てみればこの有様。

うーん、、、、、、、、大失敗でし。なんでかなー。しくしく。

*1:どうでもいい べき乗 の話 - みねこあ