Pyevolveで学ぶ遺伝的アルゴリズム
Pyevolveとは、Pythonで書かれた遺伝的アルゴリズムのフレームワークです。公式サイトによれば、Pyevolveの方針は、
・ pure python で書く
・ APIを簡単に使えるようにする
・ 進化過程をグラフ等で見れる
・ 拡張性をもたせる
・ パフォーマンスを第一にデザインする
・ 一般的な特徴を持つ
・ 全てのオプションにデフォルト値がある
・ open-source (GPL)
…とのこと。遺伝的アルゴリズムは前から少し興味があったので、これを使って勉強してみたいと思います。
Pyevolveは日本語はおろか英語での情報もあまりありません。以下に参考になるサイトを挙げておきます。
公式サイト
Welcome to Pyevolve documentation ! — Pyevolve v0.5 documentation
日本語でPyevolveについて書かれた唯一(?)のページ。ナップザック問題と巡回セールスマン問題をPyevolveで解いています。
pyevolveによる遺伝的アルゴリズム(1) - Pashango’s Blog
Pyevolveを紹介している講演の映像(このブログはPyevolveの作者の方が書いてます)
http://pyevolve.sourceforge.net/wordpress/?p=1309
遺伝的アゴリズムそのものについて
遺伝的アルゴリズム
村上・泉田研究室 遺伝的アルゴリズム
予め言っておくと、私は遺伝的アルゴリズムについて全く勉強したことがないので、以下の記述に間違いがあるかもしれません^^;
遺伝的アルゴリズム?遺伝的プログラミング?
まぁまずそもそも遺伝的アルゴリズムって何?という話ですが、それについては以下サイトが参考になります。
簡単にアルゴリズムの流れをまとめると、
1. 初期化 (要素Nの集団を生成)
2. 集団の評価 (収束判定)
3. 交叉[組み換え] (子孫を生成)
4. 突然変異 (一定率に応じて個体を変化させる)
5. 評価 (各個体の適応度を評価)
6. 淘汰 (子孫を生成した分、数が増えたので、適応度を目安に数を減らす)
7. 2に戻る
となります(これだけ言われてもなんのことやら…)。手順は入れかわることもあるようです。
で、てっきり自分は遺伝的アルゴリズムと遺伝的プログラムって同じでしょ、と思ってたんですが、実は違います。
遺伝的アルゴリズムのデータ構造が配列であるのに対して、木構造を使ったものが遺伝的プログラムというそうです。
Pyevolveのインストール
まぁ言葉とかいきなりたくさん定義されても訳が分からなくなるので、とりあえずプログラムを作ってみましょう。
Pyevolveは、ここからダウンロードできます。
ダウンロードしたあと、sudo python setup.py install とかすればインストールされます。Python2.5か2.6で動きます。付属のexampleを利用していきます。
Pyevolveを使ってみる
それでは、公式のチュートリアルにしたがってプログラムを作成していきます(英語読める人は普通にチュートリアル読んだほうがいいです^^;)。ちなみにバージョン0.5と0.6rc1とで微妙に変数名が違ったりするみたいです。
今回求めたいことは、ずばり
・要素数20の配列xがあるとき、以下の関数f(x)が最大となるxを求めよ
見てのとおり、上の関数f(x)は、xの要素のうち0であるものの数を求めています(つまりf([0,1,3,4,0,5]) = 2)。
そんなこと言われなくても全部0の配列に決まってんじゃん!と思うかもしれませんが、まぁ最初は簡単な問題からです。
まず始めに、評価関数を作成します。
def eval_func(chromosome): score = 0.0 for value in chromosome: if value == 0: score += 1 return score
ただ上の関数を定義しただけです。scoreというのがf(x)の値ですが、これがようするに評価(適応度)となります。chromosomeは日本語でいうと染色体。いきなり染色体と言われてなんのこっちゃって感じですが、ここではこの染色体は1次元配列のことです。この染色体がそれぞれ情報(遺伝子)を持ってるわけです。
遺伝的アルゴリズムでは、まず複数の個体(染色体)を生成し、それを評価してその中から子孫をいくつか作って…ということを繰替えしていきます。
それでは次に、プログラムで利用する、ゲノムのインスタンスを生成します。ゲノムとは全遺伝情報ってことですが…ここでは個体(染色体)でいいんだと思います。
# ゲノムの生成 genome = G1DList.G1DList(20)
1Dとは1 Demension のことで、つまり1次元配列のことです。遺伝情報を格納する方法はいくつかあって、
・値が0か1のみの1次元配列 (G1DBinaryString)
・1次元配列 (G1DList)
・2次元配列 (G2DList)
・木構造
などです。ここでは1次元の配列の構造を利用します(もちろんgenomeは実際にはクラスのインスタンスです)。
それから、このゲノムに対して設定をしていきます。
# 評価関数をセット genome.evaluator.set(eval_func)
evalutor.set()で先程の評価関数をゲノムに結びつけます。他にもいろいろオプションがありますが、他の値は設定しなければデフォルトの値が利用されます。
たとえば、ここで配列を初期化していませんが(初期化関数を設定することになります)、この場合は、pyevolveで定義される[Consts.CDefRangeMin, Consts.CDefRangeMax]の間の数になります。
ここではこれ以上の設定はしません。
次に、実際に進化をおこなう環境となる遺伝的アルゴリズムエンジン(Genetic Algorithm Engine)を作成します。
# 評価をするGAEngine を作成
ga = GSimpleGA.GSimpleGA(genome)
これで完成です。あとは動かすのみです。
交叉の方法も突然変異の方法も選択の方法も指定してませんがデフォルトの値が利用されます。
# 評価開始 ga.evolve(freq_stats=10) # 10世代おきに途中経過を表示 # 最も評価の高い個体を表示 print ga.bestIndividual()
ga.evolve()で評価開始です。freq_statsを指定すればそれに応じて途中経過が出力されます。
ga.bestIndividual()を出力することで結果を得ます。
以下がソースをまとめたものです。
#!/usr/bin/env python2.6 # -*- coding: utf-8 -*- # f(x) = ��{(x[i] == 0) ? 1 : 0} # の評価関数で最大の評価が得られるものを探す from pyevolve import G1DList from pyevolve import GSimpleGA # 評価関数 def eval_func(chromesome): score = 0.0 for value in chromesome: if value == 0: score += 1 return score # Genome Instance genome = G1DList.G1DList(20) # 評価関数をセット genome.evaluator.set(eval_func) # 評価をするGAEngine を作成 ga = GSimpleGA.GSimpleGA(genome) # 評価開始 ga.evolve(freq_stats=10) # 10世代おきに途中経過を表示 # 最も評価の高い個体を表示 print ga.bestIndividual()
出力結果はこんな感じです。
Gen. 0 (0.00%): Max/Min/Avg Fitness(Raw) [0.25(2.00)/0.21(0.00)/0.21(0.21)]
Gen. 10 (10.00%): Max/Min/Avg Fitness(Raw) [10.80(10.00)/7.20(8.00)/9.00(9.00)]
Gen. 20 (20.00%): Max/Min/Avg Fitness(Raw) [19.00(19.00)/19.00(19.00)/19.00(19.00)]
Gen. 30 (30.00%): Max/Min/Avg Fitness(Raw) [20.00(20.00)/20.00(20.00)/20.00(20.00)]
Gen. 40 (40.00%): Max/Min/Avg Fitness(Raw) [20.00(20.00)/20.00(20.00)/20.00(20.00)]
Gen. 50 (50.00%): Max/Min/Avg Fitness(Raw) [20.00(20.00)/20.00(20.00)/20.00(20.00)]
Gen. 60 (60.00%): Max/Min/Avg Fitness(Raw) [20.00(20.00)/20.00(20.00)/20.00(20.00)]
Gen. 70 (70.00%): Max/Min/Avg Fitness(Raw) [20.00(20.00)/20.00(20.00)/20.00(20.00)]
Gen. 80 (80.00%): Max/Min/Avg Fitness(Raw) [20.00(20.00)/20.00(20.00)/20.00(20.00)]
Gen. 90 (90.00%): Max/Min/Avg Fitness(Raw) [20.00(20.00)/20.00(20.00)/20.00(20.00)]
Gen. 100 (100.00%): Max/Min/Avg Fitness(Raw) [20.00(20.00)/20.00(20.00)/20.00(20.00)]
Total time elapsed: 0.767 seconds.
- GenomeBase
Score: 20.000000
Fitness: 20.000000
Params: {}
Slot [Evaluator] (Count: 1)
Name: eval_func - Weight: 0.50
Slot [Initializator] (Count: 1)
Name: G1DListInitializatorInteger - Weight: 0.50
Doc: Integer initialization function of G1DList
This initializator accepts the *rangemin* and *rangemax* genome parameters.
Slot [Mutator] (Count: 1)
Name: G1DListMutatorSwap - Weight: 0.50
Doc: The mutator of G1DList, Swap Mutator
.. note:: this mutator is :term:`Data Type Independent`
Slot [Crossover] (Count: 1)
Name: G1DListCrossoverSinglePoint - Weight: 0.50
Doc: The crossover of G1DList, Single Point
.. warning:: You can't use this crossover method for lists with just one element.
- G1DList
List size: 20
List: [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]おーちゃんと要素が0だけの配列が表示されてます。
それではここでデフォルトの値と、交叉・選択・突然変異の方法についてちょこっと説明しときます。
配列の要素に使われる値を指定したい場合は、
genome.setParams(rangemin=0, rangemax=10)
みたいにします。
また、集団の数(Population)は80、世代は100まで計算をします。突然変異率は2%、交差率は80%です。
G1DListの場合、交差方法は1点交叉です。これは、染色体の一点を選んでそれをいれかえます。
010101|000 --> 010101|111
111010|111 --> 111010|000
その他の交叉方法として、2点を選んでその間をいれかえる2点交叉、
010|1010|00 --> 010|0101|00
111|0101|11 --> 111|1010|11
マスクを利用して交叉する一様(Uniform)交叉
110100001(マスク)
010101000 --> 011110110 (マスクが1のとき上段の親から、マスクが0のとき下段の親から受け継ぐ)
111010111 --> 110001001 (上の逆)
などがあります。
突然変異の方法は、デフォルトでは遺伝子をランダムに入れ変えます。他の方法としては、binary stringの場合、特定の遺伝子を1にしたり、といったことがあります。
選択の方法は、デフォルトではランキング方式(エリート選択方式)です。つまり、評価の高いものから選んでいくわけです。
他の方法としてはルーレット方式、トーナメント方式があります。ルーレット方式では、適応度に応じて当たりやすくなるようなルーレットを作って、それを元に選択していきます。
…といことで長くなりましたが、Pyevolveを使えば、求めたい問題を適切に遺伝子型に変形してあげて、評価関数を作成すれば遺伝的アルゴリズムが利用できます。
関数の最小値を求める
それでは次に関数の最小値を求めてみたいと思います。example3のpyevolve_ex3_schaffer.py のコードを利用します。
ここでは Schaffer F6 Function という関数の最小値を求めます。初めて聞いた名前ですけど、
という関数らしいです。これをグラフであらわしてみたものはこちらを参照してみてください^^;
グラフを見ると分かるとおり、極値がたくさんあります。関数の最小値を求める方法として、最小点が存在する区間を狭めていく黄金分割法がありますが、その方法では最小を見つけることは難しいです(極値を求める可能性が高い)。そこで、遺伝的アルゴリズムを利用します。
今度もさっきと同じような感じですが、少しデフォルト値を変更してます。該当箇所にコメントつけました。
# -*- coding:utf-8 -*- from pyevolve import G1DList, GSimpleGA, Selectors from pyevolve import Initializators, Mutators, Consts import math # Schaffer F6 Function # 評価関数 def schafferF6(genome): t1 = math.sin(math.sqrt(genome[0]**2 + genome[1]**2)); t2 = 1.0 + 0.001*(genome[0]**2 + genome[1]**2); score = 0.5 + (t1*t1 - 0.5)/(t2*t2) return score def run_main(): # Genome instance genome = G1DList.G1DList(2) # x,yの2つをあらわす # rangemax,rangemin … xとyの範囲 genome.setParams(rangemin=-100.0, rangemax=100.0, bestrawscore=0.0000, rounddecimal=4) # G1DListをrangemin,rangemaxの間で初期化 genome.initializator.set(Initializators.G1DListInitializatorReal) # rangemin,rangemaxの値の間で突然変異する genome.mutator.set(Mutators.G1DListMutatorRealGaussian) #評価関数のセット genome.evaluator.set(schafferF6) ga = GSimpleGA.GSimpleGA(genome) # 選択方法はルーレット方式 ga.selector.set(Selectors.GRouletteWheel) ga.setMinimax(Consts.minimaxType["minimize"]) # 評価関数を最小化するようにする ga.setGenerations(8000) # 8000世代まで計算 ga.setMutationRate(0.05) # 突然変異率 0.05 ga.setPopulationSize(100) # 集団の数は100 ga.terminationCriteria.set(GSimpleGA.RawScoreCriteria) # 個体から得られるbestrawscoreによって進化を打ち切る # 評価開始,250世代ごとに途中経過を表示 ga.evolve(freq_stats=250) # 最も良いものを表示 best = ga.bestIndividual() print best print "Best individual score: %.2f" % best.getRawScore() if __name__ == "__main__": run_main()
bastrawscoreとrounddecimalについて補足しておくと、今回は関数の性質上値は正になるので、評価関数の値が小数点第四位を四捨五入した時に0ならばそこで進化(計算)を打ち切ります(打ち切るように実際に指定してるのはga.terminationCriteria.set())。また、今回は評価関数の値は低いほうがいいので、ga.setMinimax(Consts.minimaxType["minimize"])としています。
自分が実行してみた結果はx = 0.0014013025801391521, y = -0.004028033222658399 で最小値0.00 みたいな感じになりました。実際この関数はx=y=0の時最小値0をとります。関数の最小値を求めるのに遺伝的アルゴリズムはいいみたいですね。
こんな感じでだいたいPyevolveの使い方は分かってきたんじゃないでしょうか?
nクイーン問題
今度はnクイーン問題に挑戦してみます。n個のクイーンをどの駒も他のに取られるような位置におかないようにするっていうあれです。ソースはexample21の pyevolve_ex21_nqueens.py です。
まぁたまには普通に遺伝的アルゴリズムを使わないでやってみます。う〜んこんな感じでしょうか…
#!/usr/bin/python import sys N = 8 def f(queen,i): if i == N: if check(queen): print queen sys.exit(0) else: for j in xrange(N): queen[i] = j f(queen,i+1) def check(queen): for i in xrange(N-1): for j in xrange(i+1,N): if (queen[i] == queen[j] or abs(queen[i]-queen[j]) == j-i): return False return True f([0]*N,0)
なんか予想以上に短く書けてしまいましたが(w)、まぁこの実装はかなり手抜きです。サイズが8ならまだいいですが、大きな数になるといくらまっても計算が終わりません(もちろんもっといい実装方法はありますが)。
Pyevolveでやるとこんな感じになります。
# -*- coding:utf-8 -*- from pyevolve import G1DList from pyevolve import Mutators, Crossovers from pyevolve import Consts, GSimpleGA #from pyevolve import DBAdapters from random import shuffle # 盤のサイズ(64 × 64) BOARD_SIZE = 64 # 評価関数 # collisionsが少ない方がよい def queens_eval(genome): collisions = 0 for i in xrange(0, BOARD_SIZE): # 駒が置かれていない行がある場合 if i not in genome: return 0 for i in xrange(0, BOARD_SIZE): col = False # 他の駒と当たらないかチェック for j in xrange(0, BOARD_SIZE): if (i != j) and (abs(i-j) == abs(genome[j]-genome[i])): col = True if col == True: collisions +=1 return BOARD_SIZE-collisions # 初期化関数 # ランダムにクイーンを並べる def queens_init(genome, **args): genome.genomeList = range(0, BOARD_SIZE) shuffle(genome.genomeList) def run_main(): genome = G1DList.G1DList(BOARD_SIZE) genome.setParams(bestrawscore=BOARD_SIZE, rounddecimal=2) genome.initializator.set(queens_init) genome.mutator.set(Mutators.G1DListMutatorSwap) genome.crossover.set(Crossovers.G1DListCrossoverCutCrossfill) genome.evaluator.set(queens_eval) ga = GSimpleGA.GSimpleGA(genome) ga.terminationCriteria.set(GSimpleGA.RawScoreCriteria) ga.setMinimax(Consts.minimaxType["maximize"]) ga.setPopulationSize(100) ga.setGenerations(250) ga.setMutationRate(0.02) ga.setCrossoverRate(1.0) ga.evolve(freq_stats=10) best = ga.bestIndividual() print best print "Best individual score: %.2f\n" % (best.getRawScore(),) if __name__ == "__main__": run_main()
まず初期化の際に各列各行に1つずつクイーンをランダムに置きます。あとは他の例と同じ感じです。自分の環境の場合だいたい40秒程度で答えがでました(bestrawscoreで計算を打ち切っているのできちんと解がでています)。
ナップサック問題
それでは今度は遺伝的アルゴリズムでよく出てくるナップサック問題を解いてみます。ナップサック問題とこの後の巡回セールスマン問題については初めに紹介したこのページに日本語での解説がありますが、こちらでもやってみます。
ナップサック問題とは、重さと価格が決まっている複数の商品があるとき、ナップサックの重量制限内で商品の合計金額が最大になる組合せを求める問題です。単純に考えると、単位重さあたりの価格が最大のものから順に選んでいけばよさそうですが、必ずしもそうはならないのがこの問題の特徴です。ナップサック問題のデータと評価関数はこのページのものを利用することにします。評価関数は組合せが重量制限内なら各商品の価格の和、それ以外は1です。
↓名前、重さ、価格のcsvデータ(名前は今回は連番です)
1,2,21 2,10,22 3,7,28 4,2,21 5,4,12 6,9,24 7,10,15 8,7,2 9,8,25 10,5,28 11,3,4 12,10,22 13,9,36 14,8,2 15,8,7 16,5,40 17,7,14 18,3,40 19,9,33 20,7,21 21,2,28 22,10,22 23,7,14 24,9,36 25,7,28 26,2,21 27,10,18 28,4,12 29,9,24 30,10,15 31,4,21 32,7,2 33,8,25 34,5,28 35,2,28 36,3,4 37,10,22 38,9,36 39,7,31 40,8,2 41,8,7 42,5,40 43,7,14 44,5,4 45,7,28 46,3,40 47,9,33 48,7,35 49,7,21 50,9,20
#!/usr/bin/python2.6 # -*- coding:utf-8 -*- from pyevolve import G1DBinaryString from pyevolve import GSimpleGA import csv # ナップサックにはいるのは40kgまで MAX_SIZE = 40 data = [] # 評価関数 def eval_func(chromosome): p = w = 0 for i,v in enumerate(chromosome): if v == 1: w += data[i]['weight'] p += data[i]['price'] if w > MAX_SIZE: return 1 else: return p def run_main(): genome = G1DBinaryString.G1DBinaryString(len(data)) genome.evaluator.set(eval_func) ga = GSimpleGA.GSimpleGA(genome) ga.setPopulationSize(200) ga.setGenerations(1000) ga.setMutationRate(0.5) ga.setCrossoverRate(1.0) ga.evolve(freq_stats=100) best = ga.bestIndividual() p = w = 0 for i in xrange(len(best)): if best[i] == 1: print data[i]['name'], print "price:", data[i]['price'], print "size:",data[i]['weight'] p += data[i]['price'] w += data[i]['weight'] print "total price:",p,"total weight:",w if __name__ == "__main__": # データの読み込み for name,weight,price in csv.reader(open("data.csv","rb")): data.append({'name':name,'price':int(price),'weight':int(weight)}) run_main()
今回は個体を表すのにG1DBinaryStringを利用しています。これは値が1か0のみの配列です。data[]の各要素に対して、1ならばそれを利用する、0ならば利用しない、という風に情報をあらわしている訳です。さすがにこんだけの量になると、これだけの計算では最適解はでません。自分が何回かやってみた感じ500〜700ちょいぐらいまでの値が出ました。もっとよい解を見つけるには、単純に計算回数(世代)を増やしてみたり、突然変異率や評価関数を変えてみたり…といったことが必要になります。いかによい設定をするかというのが重要です(どうすればいいんでしょうね…)。
巡回セールスマン問題
最後は、これもよく遺伝的アルゴリムの例にでてくる巡回セールスマン問題(Travelling Salesman Problem : TSP)を解いてみます。ソースはexample12の pyevolve_ex12_tsp.py です。
巡回セールスマン問題とは、複数の都市があるとき、その全ての都市を一度ずつ巡って出発点にもどる経路が最小になるものを求める問題です。
まずは、どのように経路情報を遺伝子として保持するかを決めます。ここでは、単純にG1DListに順番ずつに都市番号を入れていくことにします(最後は出発点に戻ってくる)。
また、交叉方法ですが、各都市には1度しか行かないので、遺伝子に同じ値があってはいけません。そこで、Crossovers.G1DListCrossoverEdgeを使います。これは Edge recombination operator と呼ばれる方法で、各頂点ではなく頂点間の繋りを重視して現在ある経路から、違う経路を作ります。
Wikipediaの例をそのまま使うと、
1: -> [A] <-> [B] <-> [C] <-> [E] <-> [F] <-> [D] <- 2: -> [C] <-> [A] <-> [B] <-> [D] <-> [E] <-> [F] <-
という2つの経路があったとします。これから、各頂点がどの要素と隣りあわせなのか表を作ります(隣接行列)。
1 2 A: B D A: C B B: A C B: A D C: B E C: F A D: F A D: B E E: C F E: D F F: E D F: E C
これらの論理和をとると、
A: B C D = {B,D} ∪ {C,B}
B: A C D = {A,C} ∪ {A,D}
C: A B E F = {B,E} ∪ {F,A}
D: A B E F = {F,A} ∪ {B,E}
E: C D F = {C,F} ∪ {D,F}
F: C D E = {E,D} ∪ {E,C}
となります。これで、新たな隣接行列ができました。これを基に、次のアルゴリズムに従って新たな経路を作ります。
1. ランダムに1つ頂点を選び、経路情報に追加
2. 経路情報が埋まったら終了。
3. 経路情報の最後の要素(頂点)を取り出す
4. その頂点について隣接要素が存在すれば、その各要素の中で隣接要素が最もすくない要素を選ぶ(それか、適当に選ぶ)
そうでなければ、まだ経路情報に追加されていない頂点の中から適当に選び、それを経路情報に追加
5. 2に戻る
う〜ん自分の日本語能力がすごく残念な感じなので分かりにくい…。
例えば、例の2つの経路からABDFCEという経路ができた場合、次のようなことが行なわれています。
A={B,C,D}でAの隣接要素がB,C,Dであることをあらわしています。
1 () -> A. Aをランダムに選択。全ての隣接情報からAを消す。A={B,C,D}で、B={C,D}、C={B,E,F}、D={B,E,F}なので最も隣接要素のすくないBを選ぶ。
2 AB. B={C,D}。C={E,F}、D={E,F}で隣接要素の数が同じなので、適当にDを選ぶ。
3 ABD. D={E,F}。E={C,F}、F={C,E}で隣接要素の数が同じなので、適当にFを選ぶ。
4 ABDF. F={C,E}。C={E}、E={C}なので適当にCを選ぶ。
5 ABDFC. C={E}、E={}なのでEを選ぶ.
6 ABDFCE. 完成。
また、今回は経路を示すリストを受けとっても何がなんだかわからないので、PILを使ってグラフを作成します。グラフを作成する関数はwrite_tour_to_img()です。今回はPILについての説明は避けます。
グラフは、最初と最後、それに途中の適当な所で作成します。途中でグラフを表示させるにはeval_func()が終わったあとに呼びだされるコールバック関数を作成し、ga.stepCallback.set()で設定します。グラフの保存場所として、プログラムのある場所にtspimgという名前のディレクトリを作る必要があります(プログラムで作ってもいいですが^^;)。
今回のプログラムは今までのものよりもちょっと長いですが、実際の処理を見れば特に変わりはありません。
# -*- coding: utf-8 -*- from pyevolve import G1DList from pyevolve import GSimpleGA from pyevolve import Mutators from pyevolve import Crossovers from pyevolve import Consts import sys, random random.seed(1024) from math import sqrt PIL_SUPPORT = None try: import Image, ImageDraw, ImageFont PIL_SUPPORT = True except: PIL_SUPPORT = False # 都市間の距離を格納する配列 cm = [] # 各街の位置(座標)を格納する配列 coords = [] # 街の数 CITIES = 100 # 街の大きさ WIDTH = 1024 HEIGHT = 768 # 最後の成績(コールバック関数で利用) LAST_SCORE = -1 # 全ての都市間の距離を計算 def cartesian_matrix(coords): matrix={} for i,(x1,y1) in enumerate(coords): for j,(x2,y2) in enumerate(coords): dx, dy = x1-x2, y1-y2 dist=sqrt(dx*dx + dy*dy) matrix[i,j] = dist return matrix # 評価関数 # 経路の合計の長さを求める # matrix[(i,j)]で都市iから都市jまでの距離が求まる def tour_length(matrix, tour): total = 0 t = tour.getInternalList() for i in range(CITIES): j = (i+1)%CITIES total += matrix[t[i], t[j]] return total # PILでグラフを作成 def write_tour_to_img(coords, tour, img_file): padding=20 coords=[(x+padding,y+padding) for (x,y) in coords] maxx,maxy=0,0 for x,y in coords: maxx, maxy = max(x,maxx), max(y,maxy) maxx+=padding maxy+=padding img=Image.new("RGB",(int(maxx),int(maxy)),color=(255,255,255)) font=ImageFont.load_default() d=ImageDraw.Draw(img); num_cities=len(tour) for i in range(num_cities): j=(i+1)%num_cities city_i=tour[i] city_j=tour[j] x1,y1=coords[city_i] x2,y2=coords[city_j] d.line((int(x1),int(y1),int(x2),int(y2)),fill=(0,0,0)) d.text((int(x1)+7,int(y1)-5),str(i),font=font,fill=(32,32,32)) for x,y in coords: x,y=int(x),int(y) d.ellipse((x-5,y-5,x+5,y+5),outline=(0,0,0),fill=(196,196,196)) del d img.save(img_file, "PNG") print "The plot was saved into the %s file." % (img_file,) # 初期化関数 # ランダムに初期化 def G1DListTSPInitializator(genome, **args): lst = [i for i in xrange(genome.getListSize())] random.shuffle(lst) genome.setInternalList(lst) # eval_func()の後に呼び出されるコールバック関数 # 世代が100で割りきれて、かつ今までよりも良い結果が得られたとこに # グラフを作成する def evolve_callback(ga_engine): global LAST_SCORE if ga_engine.getCurrentGeneration() % 100 == 0: best = ga_engine.bestIndividual() if LAST_SCORE != best.getRawScore(): write_tour_to_img( coords, best, "tspimg/tsp_result_%d.png" % ga_engine.getCurrentGeneration()) LAST_SCORE = best.getRawScore() return False def main_run(): global cm, coords, WIDTH, HEIGHT # 都市の作成 coords = [(random.randint(0, WIDTH), random.randint(0, HEIGHT)) for i in xrange(CITIES)] # 都市間の距離を求める cm = cartesian_matrix(coords) genome = G1DList.G1DList(len(coords)) genome.evaluator.set(lambda chromosome: tour_length(cm, chromosome)) genome.crossover.set(Crossovers.G1DListCrossoverEdge) genome.initializator.set(G1DListTSPInitializator) ga = GSimpleGA.GSimpleGA(genome) ga.setGenerations(200000) ga.setMinimax(Consts.minimaxType["minimize"]) # 評価(適応度)は低い方がよい ga.setCrossoverRate(1.0) ga.setMutationRate(0.02) ga.setPopulationSize(80) # コールバック関数の設定 ga.stepCallback.set(evolve_callback) ga.evolve(freq_stats=500) best = ga.bestIndividual() if PIL_SUPPORT: write_tour_to_img(coords, best, "tsp_result.png") else: print "No PIL detected, cannot plot the graph !" if __name__ == "__main__": main_run()
20万世代まで計算するので、このプログラムの実行はかなり時間がかかります。自分の環境だと3時間ぐらいかかりました。
↓最初
↓最後
↓せっかくなのでGIFアニメにしてみました(微妙かな…^^;)
う〜んもうちょっといい経路があるような気がしますが…まぁそれなりな経路が見つかってます。
ちなみに今回途中で作成されたグラフは87枚です。200000/500=400なので最大で400枚グラフができるわけですけど、なかなか上手くはいきませんね。
とりあえず以上で今回の話は終了です。遺伝的アルゴリズムについて勉強しようとしたとき、遺伝的アルゴリズムを説明するページはけっこうある気がしますが、じゃぁ実際どうやって作るの?といった情報はあまりないように思います。Pyevolveを使えば結構手軽に遺伝的アルゴリムができるんだなーと思いました。本当はPyevolveにはもっと機能があって、遺伝的プログラミングもできるんですが、いかんせん自分がよく分かってないのでその話はまた次の機会にでも…(^^;。まぁ本当に遺伝的アルゴリズムを勉強したいんであればPyevolveがやっているコアな部分も勉強する必要があるでしょうね。
なんだか思ったより長くなってしまました。最後になりましたが、ここまで読んでくださった方、ありがとうございました。
