2014年12月7日日曜日

巡回セールスマン問題にアニメーションを付けてみる

前回pythonで巡回セールスマン問題をとく遺伝的アルゴリズムを作成したのですが、
最終結果だけの巡回順路を与えられても、うまいこと進化しているのかがわかりにくいです。
そこで、進化の過程をアニメーションで表示させてみようと思います。
どのようなアニメーションを作りたいかというと、
一つ目は、各世代の最適解において、巡回経路がどのように進化しているのかを見れるアニメーションと、
2つ目はその最適解のスコア(適応値)が世代でどのような数値をとっているかを追っていけるアニメーションです。

アニメーションを作るには2つのやり方があります。
こちらの一つ目は実行中に完全データを作ってしまい、そのデータを使ってアニメーションを作る方法と、
2つめは実行しながら値を更新していきアニメーションを作成する方法です。

最終的に実行結果はどちらも同じようになると思いますが、
1つ目の方ではメインループで結果を配列に格納していき、ループが終了してからその配列を元にアニメーションを作成するので、作成したコードに追加するだけでできます。
それに対して、2つ目はメインループ自体をアニメーション関数に置き換えるので、僕が前回作ったようなコードでは、関数を抜本的に書き換えなければならないと思われます。
ですので今回は1つ目にやり方でやることにしました。

最適解のスコアと世代数のグラフなんかは値を更新しながら作りたかったのですが、
かたっぽを計算が終わってから、もうかたっぽを計算中にとするには力量が足らなかったので、今回はまとめて計算を行い計算結果からアニメーションを作っていきます。

では実際のコードです。
# -+- coding: utf-8 -*-
import random #ランダムモジュール
import math
import copy
import numpy as np
import pylab
import matplotlib.pyplot as plt
import matplotlib.animation as animation

def city_init():
    #都市座標
    citylist_x = [2,10,7,32,35,9,28,3,8,5,13,45,63,11,17,34,55,38,41,15,6,53,46,18,23,16,19,27,42,12]
    citylist_y = [36,22,47,41,32,24,15,2,25,28,48,5,11,16,27,8,46,35,58,14,26,21,49,34,1,7,19,43,29,39]

    return [citylist_x,citylist_y]

##評価関数
def eval_func(gean):
    x,y =city_init()
    x_cp = list(x)
    y_cp = list(y)
##    print x
    #a,b,c,d,e,f,のうち行った都市は消していく
    #c,e,a,f,d,bの順で回る場合は
    #[2,3,0,2,1,0]となる
    ##gean = [2, 3, 1, 1, 0]
    ##gean = [2,4,1,3,0]
    gean_decode = []
    geannum = 0
    while geannum<=len(geana)-1:
        num = x.index(x_cp[gean[geannum]])
        gean_decode.append(num)
        del x_cp[gean[geannum]]
        geannum +=1
    #評価の基準は距離が短いことである。
    dist_stok = []
    geannum = 0
    while geannum<=len(gean)-2:
        dist = (x[gean_decode[geannum+1]]-x[gean_decode[geannum]])**2 + (y[gean_decode[geannum+1]]-y[gean_decode[geannum]])**2
        dist_stok.append(dist)
        geannum +=1
    alldistance = sum(dist_stok)

    last_dist = (x[gean_decode[len(gean)-1]]-x[gean_decode[0]])**2 + (y[gean_decode[len(gean)-1]]-y[gean_decode[0]])**2
    alldistance = alldistance + last_dist
##    print alldistance
    return alldistance


def geneticoptimize(maxiter = 100,maximize = False,popsize = 50,popnum = 30,elite = 0.99,mutprob =0.2):
    """
    maxiter = 1, 繰り返し数
    maximize = True,    スコアを最大化
    popsize = 50,   個体数
    popnum = 10,    遺伝子数(長さ)
    elite = 0.2,    生き残る遺伝子の割合
    mutprob =0.2    突然変異のおこる確立
    """
    # 突然変異
    def mutate(vec):
        i = random.SystemRandom().randint(0,popnum-1)
        after_mut = random.randint(0,popnum-i-1)
        return vec[:i] + [after_mut]+vec[i+1:]

     # 1点交叉 非推奨
    def one_point_crossover(r1,r2):
        i = random.SystemRandom().randint(1,popnum-2)

        return random.SystemRandom().choice((r1[0:i] + r2[i:], r2[0:i] + r1[i:]))

    # 2点交叉
    def two_point_crossover(r1,r2):
        i, j = sorted(random.SystemRandom().sample(range(popnum),2))
        return random.SystemRandom().choice((r1[0:i] + r2[i:j] + r1[j:] , r2[0:i] + r1[i:j] + r2[j:]))

    # 一様交叉
    def uniform_crossover(r1, r2):
        q1 = copy.copy(r1)
        q2 = copy.copy(r2)
        for i in range(len(r1)):
            if random.SystemRandom().random() < 0.5:
                q1[i], q2[i] = q2[i], q1[i]

        return random.SystemRandom().choice([q1,q2])
    def anime_flam(gean):
        copy_x = list(citylist_x)
        copy_y = list(citylist_y)
        ordercity_x =[]
        ordercity_y =[]
        for i in gean[0][1]:
            ordercity_x.append(copy_x[i])
            ordercity_y.append(copy_y[i])
            del copy_x[i],copy_y[i]
        return ordercity_x,ordercity_y
    #遺伝子の初期化
    pop = []
    vec = []
    for i in range(popsize):
        selected_city = 1
        while selected_city <= popnum:
            vec_vec = random.SystemRandom().randint(0,popnum-selected_city)
            vec.append(vec_vec)
            selected_city +=1
        pop.append(list(vec))   #listはコピーを渡さないと変更元まで変更される
        del vec[:]


    #遺伝子の初期化

    # 交叉アルゴリズムの選択
    crossover = two_point_crossover
    #crossover = uniform_crossover

    #メインループ
    topelite = int(elite * popsize)
    score_stoc = []
    citylist_x,citylist_y = city_init()
    fig = plt.figure()
    ax1=fig.add_subplot(1,2,1)
    ax2=fig.add_subplot(1,2,2)
    ims = []
    im2_pointset =[]
    for i in range(maxiter):
        scores=[(eval_func(v),v) for v in pop]
        scores.sort()
        if maximize:
            scores.reverse()
        ranked = [v for (s,v) in scores]
        # 弱い遺伝子は淘汰される
##        pop = ranked[0:topelite]
        eliminate_cycle = 0

        while eliminate_cycle < popsize:
            if eliminate_cycle == 0:
                pop = ranked[:1]
            else:
                if random.SystemRandom().random() < elite**eliminate_cycle:
                    pop.append(ranked[eliminate_cycle])
            eliminate_cycle +=1
        topelite = len(pop)
        # 生き残った遺伝子同士で交叉したり突然変異したり
        while len(pop) < popsize:
            if random.SystemRandom().random() < mutprob:
                # 突然変異
                c = random.SystemRandom().randint(0,topelite)
                pop.append(mutate(ranked[c]))
                if random.SystemRandom().random() < mutprob:
                    selected_city = 1
                    while selected_city <= popnum:
                        vec_vec = random.SystemRandom().randint(0,popnum-selected_city)
                        vec.append(vec_vec)
                        selected_city +=1
                    pop.append(list(vec))
                    del vec[:]
            else:
                # 交叉
                c1 = random.SystemRandom().randint(0,topelite)
                c2 = random.SystemRandom().randint(0,topelite)
                pop.append(crossover(ranked[c1],ranked[c2]))
        # 暫定の値を出力
        
        score_stoc.append(scores[0][0])
        #scores[0]はスコアと遺伝子(トップの)
        #scores[0][0]はスコア
        #scores[0][1]は遺伝子
        x,y = anime_flam(scores)
        #ordercityは暫定での訪問順

        im, = ax1.plot(x,y,'o-k')
        im2, = ax2.plot(score_stoc[0:i],'o-k')

        #plotの場合は変数にカンマを付ける
        ims.append([im,im2])
    ani = animation.ArtistAnimation(fig, ims, interval=1,blit=True, repeat_delay=0)
    plt.show()
    return scores[0]



def main():
    x,y =city_init()
    ans = geneticoptimize()
    ans_x = []
    ans_y = []
    print("Ans:",ans)
    for i in ans[1]:
        ans_x.append(x[i])
        ans_y.append(y[i])
        del x[i],y[i]


if __name__ == '__main__':
    main()

特別なことは特にやっていませんが、特に苦戦したところは、
グラフの方式がplotの時はそれを代入する変数に,(カンマ)を付けなければならないようです。
この詳細は不明なのですが、とりあえず今の段階ではカンマを付けると覚えておきます。

こんなグラフを作りました。
アニメーションの保存方法を模索中なのでまだ動きません。

参考
Animation using matplotlib with subplots and ArtistAnimation
http://stackoverflow.com/questions/17853680/animation-using-matplotlib-with-subplots-and-artistanimation
LineCollection animation when updating the frames in a loop
http://stackoverflow.com/questions/20250154/linecollection-animation-when-updating-the-frames-in-a-loop
pythonで散布図アニメーションを試してみた
http://cflat-inc.hatenablog.com/entry/2014/03/17/214719
matplotlibによるサイクロイドのアニメーションにて
http://matsulib.hatenablog.jp/entry/2013/07/15/203419
matplotlibでランダムウォークをアニメーション
http://d.hatena.ne.jp/xef/20120629/p2
Python:matplotlibでScatterを使った2Dアニメーション
http://venuschjp.blogspot.jp/2013/03/pythonmatplotlibscatter2d.html
Gifアニメーション
http://ipythonnb4jpnexp.blogspot.jp/2013/09/gif-matplotlib1.html
matplotlibでアニメーション
http://blog.monophile.net/2013/12/matplotlib.html
matplotlib/plot編集する
http://seesaawiki.jp/met-python/d/matplotlib/plot

0 件のコメント:

コメントを投稿