巡回セールスマン問題といえば, xy空間上などに適当に散布している各点を通るような最も最短の経路を求める問題です. 下図の実行結果の各ノードを結ぶ線が巡回経路になります. 組み合わせが爆発するので, (点の数が極端に小さい場合を除き)総当り法では解けず, 一般には何らかの最適化アルゴリズムを使って近似解を導くわけですが, そのうちの一つが遺伝的アルゴリズムというわけです.
SGA(単純GA)
SGAの仕組みは(Wikipediaに詳しい説明がありますが)基本的に,個体の初期集団の生成 → 選択 → 交叉(交配) → 突然変異 → 選択(2回目) → 交叉(2回目) → 突然変異(2回目) →...
と続いていく処理です.
「個体の初期集団」とは解の集合で, これは最初はランダムにいくつか生成.
選択のフェーズでは, 生成された解から次世代を担う解を選択. スコアの高いものが残る確率が高くなるように選びます. 初期集団の個体と同じ数を揃えます(選択の重複や選択したもののコピーも有). 選択方法には, ルーレット選択や, ランキング選択, トーナメント選択などいくつか種類があります.
交叉のフェーズでは, 2つの個体(解)を選び, その解のそれぞれの部分列をお互いに入れ替え,新しい解を生成します. これも, 一点交叉, 二点交叉などいくつか種類があります.
突然変異のフェーズでは, 一定確率で, 個体の集合の一部を適当に書き換えます. 書き換える確率は問題によって異なります.
上記の3ステップで, 解の2世代目が生成されます. これをn回繰り返して, 1, 2, ... ... ,n世代目と解の集合がより洗練されたものになっていくという仕組み.
繰り返しの上限などを決めておき, 適当なところで止めて, その時の世代(解の集合)の中でもっともスコアの高い解が最終的な計算結果です.
初期集団をどれくらい作るのか, 何回繰り返すのか, 選択方式など各種パラメータは, 計算量や生成される最適解とにらめっこしながら実装者が試行錯誤しながら決めていくもののようです.
特に優れた解は, 選択/交叉/突然変異のプロセスを飛ばして, 次世代に残したり(エリート抽出)しますが, 今回の実装ではしてません.
Grefenstetteらの手法
今回は巡回セールスマン問題を扱うので, 解は各ノードのリストとなります. 最終的な解は以下のようなデータ構造です.[:a :f :e :g :d :h :c :b]
交叉とは例えば, 二つの個体, 01000111 と 11101101について, [010][00111] と[111][01101]に分割して, [010][01101] と [111][00111] の二つの新しい個体を生成する作業なのです. 末尾の5bitを入れ替えました. 何bitのところで交叉させるかは, ランダムに決めます.
これで新しい解が生成されるわけですが, これを上記のノードのリストに直接適用しようとすると不都合が生じかねません. 解がbit列で表現されているときは, まだいいのですが. 上記のような重複しない要素の順列が解になる場合 abdc と cdab について解の中央で交叉させると, [ab][ab]と[cd][dc]という全く意味を成さない解が生成されてしまいます.
これについて, いろいろ試行錯誤ができるのですが, Grefenstetteらの手法(特に呼び名がないのでこう呼びますが)では, ノードのリストを次のような数値列にエンコードします.
エンコードされた解が(x0, x1, ..., xn)のとき x0は, 全ノードのソートされた列P(例えば (a, b, c, ...))の x0番目の要素を表し, x1は, x0番目の要素が取り除かれた後のP'のx1番目の要素を表します.
xiは, x0からx(i-1)を取り除く作業が終わった後の, 全ノードの(整列された)列P''のxi番目の要素を表すことになります.
例えば, P = (a, b, c, d)があったときに, (3, 1, 2, 1)は, (c, a, d, b)を表すわけですが, 必ずソートされたすべての要素を含む列Pを参照するので, 要素が重複することはありません.
先頭の3は, Pの3番目cを表し, その次に来る1は, P' = (a, b, d)の1番目aを表し, 2は, P''=(b, d)の2番目dを表し, 最後の1は, P'''=(b)の1番目bを表します.
したがって, 交叉しても, 解として意味を成さない解は生成されません (この手法はこの手法でGA的に別の問題がありますが, ここでは割愛します). これにより, 交叉可能なデータへエンコードできるわけです.
この手法は, (著) 棟朝 雅晴, 遺伝的アルゴリズム -その理論と先端的手法-, 森北出版, 2008に説明があります.
これで新しい解が生成されるわけですが, これを上記のノードのリストに直接適用しようとすると不都合が生じかねません. 解がbit列で表現されているときは, まだいいのですが. 上記のような重複しない要素の順列が解になる場合 abdc と cdab について解の中央で交叉させると, [ab][ab]と[cd][dc]という全く意味を成さない解が生成されてしまいます.
これについて, いろいろ試行錯誤ができるのですが, Grefenstetteらの手法(特に呼び名がないのでこう呼びますが)では, ノードのリストを次のような数値列にエンコードします.
エンコードされた解が(x0, x1, ..., xn)のとき x0は, 全ノードのソートされた列P(例えば (a, b, c, ...))の x0番目の要素を表し, x1は, x0番目の要素が取り除かれた後のP'のx1番目の要素を表します.
xiは, x0からx(i-1)を取り除く作業が終わった後の, 全ノードの(整列された)列P''のxi番目の要素を表すことになります.
例えば, P = (a, b, c, d)があったときに, (3, 1, 2, 1)は, (c, a, d, b)を表すわけですが, 必ずソートされたすべての要素を含む列Pを参照するので, 要素が重複することはありません.
先頭の3は, Pの3番目cを表し, その次に来る1は, P' = (a, b, d)の1番目aを表し, 2は, P''=(b, d)の2番目dを表し, 最後の1は, P'''=(b)の1番目bを表します.
したがって, 交叉しても, 解として意味を成さない解は生成されません (この手法はこの手法でGA的に別の問題がありますが, ここでは割愛します). これにより, 交叉可能なデータへエンコードできるわけです.
この手法は, (著) 棟朝 雅晴, 遺伝的アルゴリズム -その理論と先端的手法-, 森北出版, 2008に説明があります.
Clojureのコード
今回はClojureでの実装が簡単そうな, 一点交叉, トーナメント選択で, 突然変異する確率は, 0.2%に設定してあります.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| (use '[clojure.core.match :only (match)] | |
| '[clojure.math.numeric-tower :as math]) | |
| ;;--------------------------------------------------------------- | |
| ;; profiler | |
| (def profiling-points (ref {})) | |
| (def temp-value (ref [])) | |
| (def profile-program true) | |
| (defmacro prof [dist-name s-exp] | |
| (if profile-program | |
| `(dosync | |
| (let [elapsed-time-message# | |
| (with-out-str (time (ref-set temp-value ~s-exp))) | |
| elapsed-time# | |
| (Double/parseDouble | |
| (nth (clojure.string/split elapsed-time-message# #" ") 2))] | |
| (if (contains? @profiling-points ~dist-name) | |
| (do | |
| (ref-set profiling-points | |
| (update-in @profiling-points [~dist-name :call-count] inc)) | |
| (ref-set profiling-points | |
| (update-in @profiling-points [~dist-name :elapsed-time] | |
| (fn [x#] (+ x# elapsed-time#))))) | |
| (ref-set profiling-points | |
| (assoc @profiling-points | |
| ~dist-name | |
| {:call-count 1 | |
| :elapsed-time elapsed-time#})))) | |
| @temp-value) | |
| s-exp)) | |
| (defn refresh-all [] | |
| (dosync (ref-set profiling-points {}))) | |
| (defn print-exec-time [] | |
| (loop [points (keys @profiling-points)] | |
| (if-not (empty? points) | |
| (let [point (first points) state (point @profiling-points)] | |
| (println point ":" | |
| "call=" (:call-count state) | |
| ",elapsed-time=" (:elapsed-time state)) | |
| (recur (rest points)))))) | |
| ;;--------------------------------------------------------------- | |
| (def sample-points-with-coord | |
| {:a [10 100] :b [30 20] :c [80 60] :d [60 100] :e [70 10] | |
| :f [5 100] :g [80 80] :h [30 30] :i [110 20] :j [0 120]}) | |
| ;;--------------------------------------------------------------- | |
| ;; find opt path | |
| ;;---------------------------- | |
| ;; Grefenstette encoding/decoding | |
| (defn symbolic->grefnum | |
| "grefenstetteらによる Symbols->Numbersの変換" | |
| [points given-route] | |
| (loop [s [] c points p given-route] | |
| (if (empty? c) (reverse s) | |
| (let [target (first p)] | |
| (recur | |
| (cons (inc (.indexOf c target)) s) | |
| (remove #(= target %) c) | |
| (rest p)))))) | |
| (defn grefnum->symbolic | |
| "symbolic->grefnumの逆変換 Numbers->Symbolsの変換" | |
| [points encoded-route] | |
| (loop [decoded-route [] c points e encoded-route] | |
| (if (empty? e) (reverse decoded-route) | |
| (let [target (nth c (dec (first e)))] | |
| (recur | |
| (cons target decoded-route) | |
| (remove #(= target %) c) | |
| (rest e)))))) | |
| ;;---------------------------- | |
| ;; selection | |
| (defn get-winner [a b] | |
| (if (< (second a) (second b)) a b)) | |
| (defn tournament-selection [scored-population] | |
| (->> (partition 2 scored-population) | |
| (map #(first (get-winner (first %) (second %)))) ;; fetch winner | |
| (repeat 2) ;; copy | |
| (apply concat))) | |
| (defn selection | |
| [scored-population] | |
| (prof | |
| :selection | |
| (doall | |
| (tournament-selection scored-population)))) | |
| ;;---------------------------- | |
| ;; crossover | |
| (defn one-point | |
| [at genes] | |
| (let [a (split-at at (first genes)) | |
| b (split-at at (second genes))] | |
| [(concat (first a) (second b)) (concat (first b) (second a))])) | |
| (defn crossover [population gene-size] | |
| (prof | |
| :crossover | |
| (doall | |
| (->> (partition 2 population) | |
| (map (partial one-point (rand-int gene-size))) | |
| (reduce concat))))) | |
| ;;---------------------------- | |
| ;; mutation | |
| (defn mutation | |
| [population gene-size] | |
| (letfn [(endamage | |
| [gene] | |
| (loop [i (dec gene-size) gene gene] | |
| (cond | |
| (= 0 i) gene | |
| (= 0 (rand-int 500)) ;; mutation (rate = 0.2%) | |
| (recur | |
| (dec i) | |
| (let [tail (drop (inc i) gene)] | |
| (concat (take i gene) [(inc (rand-int (count tail)))] tail))) | |
| :else | |
| (recur (dec i) gene))))] | |
| (prof | |
| :mutation | |
| (doall (map endamage population))))) | |
| ;;---------------------------- | |
| ;; find-opt-path | |
| (defn mean-square [x1 x2 y1 y2] | |
| (math/sqrt (+ (* (- x1 x2) (- x1 x2)) (* (- y1 y2) (- y1 y2))))) | |
| (defn find-opt-path | |
| [points-with-coord loop-limit pop-size] | |
| (let [points (sort (keys points-with-coord)) | |
| gene-size (count points) | |
| pop-size (if (odd? pop-size) (inc pop-size) pop-size) | |
| init-population | |
| (->> #(->> points shuffle (symbolic->grefnum points)) | |
| repeatedly (take pop-size))] | |
| (letfn [(get-distance ;; 二乗距離 | |
| [point-A point-B] | |
| (let [a (get points-with-coord point-A) | |
| b (get points-with-coord point-B)] | |
| (mean-square (first a) (first b) (second a) (second b)))) | |
| (evaluate ;; 適応度(=移動距離)の評価 | |
| [route] | |
| (->> (cons (last route) (drop-last route)) | |
| (map get-distance route) | |
| (reduce +))) | |
| (get-scored-routes ;; [ルートの表現, 適藤度]のペアのリスト | |
| [population] | |
| (map #(list %1 (evaluate (grefnum->symbolic points %1))) population))] | |
| ;; main routine | |
| (loop [population init-population rest-loops loop-limit] | |
| (let [scored-routes (get-scored-routes (shuffle population))] | |
| (if (< rest-loops 1) | |
| ;; get result | |
| (->> scored-routes (reduce get-winner) first (grefnum->symbolic points)) | |
| ;; selection -> mutation -> crossover repeatedly | |
| (recur | |
| (-> scored-routes selection (mutation gene-size) (crossover gene-size)) | |
| (dec rest-loops)))))))) | |
| ;;--------------------------------------------------------------- | |
| ;; graph plot | |
| (import (javax.swing JFrame) | |
| (java.awt Color) | |
| (java.awt Graphics)) | |
| (def radius 8) | |
| (defn x-extend [x] | |
| (+ (* (- x radius) 3) 50)) | |
| (defn y-extend [y] | |
| (+ (- (* (- y radius) 3)) 400)) | |
| (defn plot-graph [points-with-coord route] | |
| (def frame (JFrame. "Route")) | |
| (doto frame | |
| (.setSize 450 450) | |
| (.setVisible true) | |
| (.setResizable false)) | |
| (Thread/sleep 100) ;; wait for the generateion of Window | |
| (def graphics (.. frame (getGraphics))) | |
| (defn plot-point [name] | |
| (let [coord (name points-with-coord) | |
| r (* radius 2) | |
| x-pos (x-extend (first coord)) | |
| y-pos (y-extend (second coord))] | |
| (doto graphics | |
| (.setColor (Color. 255 100 100)) | |
| (.fillOval x-pos y-pos r r) | |
| (.setColor (Color. 100 190 190)) | |
| (.drawOval x-pos y-pos r r) | |
| (.setColor (Color. 0 0 0)) | |
| (.drawString (str name) x-pos (+ y-pos 25))))) | |
| (defn draw-line [point-a point-b] | |
| (let [xy-a (point-a points-with-coord) | |
| xy-b (point-b points-with-coord) | |
| get-x #(+ radius (x-extend (first %1))) | |
| get-y #(+ radius (y-extend (second %1)))] | |
| (doto graphics | |
| (.setColor (Color. 0 0 0)) | |
| (.drawLine (get-x xy-a) (get-y xy-a) (get-x xy-b) (get-y xy-b))))) | |
| ;; draw points | |
| (doall (map plot-point (keys points-with-coord))) | |
| ;; draw route | |
| (doall (map draw-line route (cons (last route) (drop-last route)))) | |
| (println "done!")) |
上半分は簡単なプロファイラ. 適当な使いやすい実行時間計測のためのプロファイラがなかったので自作. Clojureは, こういうのが簡単に作れるところがいいですね.
中央部がSGAの実装.
find-opt-pathがSGAの本体です. メインルーチンは, find-opt-path関数の下から数行の部分です. 選択(selection), 交叉(crossover), 突然変異(mutation)の実装はそれぞれ, 見やすいので, map/reduceで楽して書いています. 実行時間に影響がありそうですが.
実行(計算)結果の可視化は, Clojureなので, Awtを使えば簡単に表示できます. Y軸を反転させて, 座標系を拡大し, 原点を下にもってきて, 線を描くだけです.
適当にreplにロードすれば使えます.
user> (load-file "find-opt-path.clj")find-opt-path関数に, ターゲットの座標の情報と, ループ数, 個体群のサイズを指定して実行.
user> (find-opt-path sample-points-with-coord 200 200) (:a :f :j :d :g :i :c :e :b :h)(最適化されているはずの)記号列を出力します. もちろん, 実行毎によって結果は変化しますが, 大体似たような結果になることが多いです.
実行結果の可視化は, plot-graphとfind-opt-pathを組み合わせて実行します.
user> (plot-graph sample-points-with-coord (find-opt-path sample-points-with-coord 200 200))
表示結果はこんな感じ.
気になる実行時間ですが,
user> (time (find-opt-path sample-points-with-coord 200 200)) "Elapsed time: 1849.978839 msecs" user> (print-exec-time) :crossover : call= 200 ,elapsed-time= 1651.9376810000003 :mutation : call= 200 ,elapsed-time= 0.5620599999999997 :selection : call= 200 ,elapsed-time= 3.0204269999999958 nil
書いている途中で気がついたのですが, 上記のプログラムは大体, パラメータさえ決まってしまえば, 使うメモリ空間は一定でほとんどin-placeアルゴリズムなわけです. ところが, 関数型言語風に書いたため, ひたすらconsingばかりしています. ちょっと効率悪いかなと.
例えば, 突然変異は, 個体を表す数値を少し書き換えればよいのですが, ばらして, (変異させるbitを)選んで, 再構成して, という処理で, 無駄な部分がいくつかあります. 特に実行時間の長かった交叉(crossover)は, リストの再構成にconcatを多用していて, 破壊的な処理に書き直したくなりました.
というわけでClojureは, こういうガリガリ回す数値計算に向いていないのか, と思ったのですが, これはどうやら間違いでした.
make-arrayを使えばJavaの配列が使えて, それなりに早いようです.
Java の配列を利用してメモリと時間の壁を突破する． (tnoda-clojure)
というわけでやり直し.