読者です 読者をやめる 読者になる 読者になる

俺の備忘録

個人的な備忘録です。

Pythonとmatplotlibで粒子群最適化アルゴリズムをアニメーション表示してみた

python

はじめに

オセロのAI作成は一回お休み。前回 => (Pythonの勉強がてらオセロAIを作ってみる③ - 俺の備忘録)
今回は、PythonPythonのグラフ描画ライブラリmatplotlibで、粒子群最適化アルゴリズムの最適化過程をアニメーション表示してみた。

粒子群最適化アルゴリズムとは

粒子群最適化アルゴリズムとは(Particle Swarm Optimization、以下PSO)とは、遺伝的アルゴリズムと同様のメタヒューリティクスな最適化アルゴリズムの一つで、ナップザック問題や巡回セールスマン問題などの最適化問題に用いられる。

このアルゴリズムは、虫や鳥の行動を模して考案されたもので、 ざっくり言うと、誰か1匹が良い解を見つけると、その近辺にいい解があるだろうと思って、他の虫や鳥も集まってくるという虫/鳥の知恵?修正?を真似したものになっている。 例えば、海で普段カモメは気ままにそれぞれ餌を探しているんだけど、 とある一匹がイワシの群れを見つけたら、他のカモメもその近辺に群がってくるって感かな。

図に表すと以下みたいな感じ。 f:id:magayengineer:20160407233227p:plain

f:id:magayengineer:20160407233235p:plain

f:id:magayengineer:20160407233238p:plain

f:id:magayengineer:20160407233241p:plain

f:id:magayengineer:20160407233244p:plain

このアルゴリズムのいいところは、なんといっても”簡単”であること。 遺伝的アルゴリズム等の知識がなくても、2〜3分あれば、きっとほとんどの人が理解できると思う。8年ぶりぐらいに触ってみたけど、すんなり再理解できた。

PSOアルゴリズム

■数式
PSOの数式は以下で表される。 要は、粒子は移動するんだけど、その時に現在の速度による影響と、その粒子の今まで見つけたベストの場所と、全体で今まで見つけたベストの場所を勘案しながら探索する。 こうすることで今まで見つけたベストの場所に全ての粒子が一直線で進まなくなり、 今まで見つけたベストの場所の周辺等をいい感じに探索することができる。


x = x + v


v = wv + c_1 r_1 (pbx-x) + c_2 r_2 (gbx-x)


x:粒子の現在値


v:現在の速度


w:慣性定数。次の速度を計算する際に現在の速度をどれだけ重視するか


r_1:0から1の乱数


r_2:0から1の乱数


c1:定数。パーソナルベストをどれだけ重視するか


c2:定数。グローバルベストをどれだけ重視するか


pbx:パーソナルベストの位置。つまり、その粒子が今までに見つけた最も良い位置。


gbx:グローバルベストの位置。つまり、全体で今までに見つけた最も良い位置。

■計算手順
具体的な計算手順は以下の通り。

  1. 全ての粒子の現在位置xとvをランダム値で初期化する。
  2. pbxは、各粒子の現在位置xで初期化
  3. gbxは、現在位置xの中で最も良い位置で初期化(最も良い位置は評価関数次第)
  4. 以下ループ。適当な回数で打ち切る。
    1. 前述の式で全粒子の現在位置更新。
    2. 評価関数で全粒子の現在位置を評価。
    3. パーソナルベストより評価が良くなった粒子はパーソナルベストを更新。
    4. グローバルベストより評価が良くなった粒子がいたらグローバルベストを更新。
    5. 各粒子の速度更新。

実装

PSOの粒子を表すクラス。 共通ロジック(座標と速度の更新)の実装と、I/Fの定義。 評価関数ごとに、このクラスのサブクラスを作成し、 評価関数ごとの都合によって変わる部分を実装する。

class BaseParticle:
    """
    PSOの基礎クラス
    """

    def __init__(self, w, x, v, c1, c2):
        """
        コンストラクタ
        :param w: 慣性定数
        :param x: 現在地(numpy.ndarray)
        :param v: 速度  (numpy.ndarray)
        :param c1:群のうちで良い位置に向かう粒子の割合
        :param c2:群のうちで良い位置に向かう粒子の割合
        :return 粒子のインスタンス
        """
        self._x = x
        self._v = v
        self._w = w
        self._c1 = c1
        self._c2 = c2

        # パーソナルベストを初期化
        self._p_best_x = x.copy()
        self._p_best_score = self.eval()

    def update(self, g_best_x):
        """
        粒子の更新を行う
        :param g_best_x:グローバルベスト
        :return:void
        """
        # 粒子の位置を更新
        self._x = self._x + self._v

        # 粒子の位置を補正
        self.clip_x()

        # 粒子の速度を更新
        r1 = np.random.rand()
        r2 = np.random.rand()
        # v = wv + c1r1(px - x) * c2r2 (gp - x)
        self._v = self._w * self._v + self._c1 * r1 * (self._p_best_x - self._x) + self._c2 * r2 * (g_best_x - self._x)

        # 速度を補正
        self.clip_v()

        # 評価
        score = self.eval()

        # パーソナルベスト更新
        if not self.is_better_than(score):
            self._p_best_score = score
            self._p_best_x = self._x.copy()

    @property
    def x(self):
        """
        :return: 現在の位置
        """
        return self._x

    @property
    def p_best_x(self):
        """
        :return: パーソナルベストの位置
        """
        return self._p_best_x

    @property
    def p_best_score(self):
        """
        :return:パーソナルベストの値
        """
        return self._p_best_score

    def eval(self):
        """
        評価。ここに評価関数を記述する。
        :return:評価値
        """
        pass

    @classmethod
    def create_particle(cls, w, x, v, c1, c2):
        """
        粒子作成用のファクトリメソッド
        :param w: 慣性定数
        :param x: 現在地(numpy.ndarray)
        :param v: 速度  (numpy.ndarray)
        :param c1:群のうちで良い位置に向かう粒子の割合
        :param c2:群のうちで良い位置に向かう粒子の割合
        :return 粒子のインスタンス
        """
        pass

    def is_better_than(self, score):
        """
        パーソナルベストが、引数scoreより良いか判定
        :param score:比較対象の値
        :return:パーソナルベストのほうが良ければtrue
        """
        pass

    def clip_x(self):
        """
        位置を補正する
        :return: void
        """
        pass

    def clip_v(self):
        """
        速度を補正する
        :return: void
        """
        pass

    @staticmethod
    def show_animation(animation_data):
        """
        最適化の過程をアニメーション表示する
        :param animation_data: アニメーション用データ
        :return: void
        """
        pass

各粒子の初期化、位置・速度・グローバルベストの更新制御処理

def calc(particle_class, times, num_of_particle, show_animation=False):
    """
    計算実行
    :param particle_class: 粒子のクラス
    :param times: 何回試行するかするか
    :param num_of_particle: 粒子数
    :param show_animation:グラフ表示有無
    :return: void
    """

    # グローバルベスト
    global_best_x = None
    global_best_score = None

    # 粒子郡
    particles = []

    if show_animation:
        animation_data = AnimationData()

    # 粒子初期化
    for n in range(num_of_particle):

        particle = particle_class.create_particle()
        particles.append(particle)

        if n == 0:
            global_best_x = particle.p_best_x
            global_best_score = particle.p_best_score
        else:
            if particle.is_better_than(global_best_score):
                global_best_score = particle.p_best_score
                global_best_x = particle.p_best_x

    # アニメーション表示用データ準備
    if show_animation:
        animation_data.add_best_particle_position_data(global_best_x, global_best_score)
        animation_data.add_particle_position_data(particles)

    for t in range(times):
        # 粒子更新
        for n, particle in enumerate(particles):
            particle.update(global_best_x)
            if particle.is_better_than(global_best_score):
                global_best_score = particle.p_best_score
                global_best_x = particle.p_best_x

        # アニメーション表示用データ準備
        if show_animation:
            animation_data.add_best_particle_position_data(global_best_x, global_best_score)
            animation_data.add_particle_position_data(particles)

        print "t =", t, "x =", global_best_x, "score =", global_best_score
    print "t =", t, "x =", global_best_x, "score =", global_best_score

    # graph表示
    if show_animation:
        particle_class.show_animation(animation_data)

最適解を求める対象の関数

Rastrigin Function


Z = 10n + \sum_i^n(x_i^2 - 10cos(2\pi x_i))


探索範囲:−5.12≦x_i≦5.12


最適解:f(0, 0, 0...)=0

n=2のときの外観は以下。 f:id:magayengineer:20160407233249p:plain

Rastrigin関数の最適解を求める粒子クラス

class RastriginParticle(BaseParticle):
    """
    Rastrigin関数
    """

    w = 0.9
    c1 = 0.9
    c2 = 0.9

    def eval(self):
        x = self._x
        x1 = x[0]
        x2 = x[1]
        return 10 * 2 + (x1 * x1 - 10 * np.cos(2 * np.pi * x1)) + (x2 * x2 - 10 * np.cos(2 * np.pi * x2))

    @classmethod
    def create_particle(cls):
        x = np.array([5.12 * 2 * np.random.rand() - 5.12, 5.12 * 2 * np.random.rand() - 5.12])
        v = np.array([5.12 * 2 * np.random.rand() - 5.12, 5.12 * 2 * np.random.rand() - 5.12])
        return RastriginParticle(cls.w, x, v, cls.c1, cls.c2)

    def is_better_than(self, score):
        if self._p_best_score < score:
            return True
        else:
            return False

    def clip_x(self):
        self._x = np.clip(self._x, -5.12, 5.12)

    def clip_v(self):
        self._v = np.clip(self._v, -5.12, 5.12)

計算結果

メイン関数。 30粒子で30ループさせてみる。

if __name__ == "__main__":
    calc(RastriginParticle, 100, 50, True)

無事に17回目のループで最適解を見つけた模様。

試行回数,   粒子の位置,                粒子の評価値
t = 0 x = [-1.05704654  0.79439026] score = 9.63082872935
t = 1 x = [ 0.86879518  0.85701424] score = 8.47060907598
t = 2 x = [-0.07911788  0.06427271] score = 2.02515850631
t = 3 x = [-0.07911788  0.06427271] score = 2.02515850631
t = 4 x = [-0.07911788  0.06427271] score = 2.02515850631
t = 5 x = [-0.07911788  0.06427271] score = 2.02515850631
t = 6 x = [-0.07911788  0.06427271] score = 2.02515850631
t = 7 x = [-0.07472074 -0.03142844] score = 1.28289421909
t = 8 x = [-0.07472074 -0.03142844] score = 1.28289421909
t = 9 x = [-0.07472074 -0.03142844] score = 1.28289421909
t = 10 x = [-0.03988766 -0.01807089] score = 0.378723517767
t = 11 x = [-0.03988766 -0.01807089] score = 0.378723517767
t = 12 x = [-0.03212093  0.01240108] score = 0.234496246395
t = 13 x = [-0.03212093  0.01240108] score = 0.234496246395
t = 14 x = [-0.03212093  0.01240108] score = 0.234496246395
t = 15 x = [-0.03212093  0.01240108] score = 0.234496246395
t = 16 x = [-0.02446902 -0.01417394] score = 0.158382181234
t = 17 x = [ 0.  0.] score = 0.0
t = 18 x = [ 0.  0.] score = 0.0
t = 19 x = [ 0.  0.] score = 0.0
t = 20 x = [ 0.  0.] score = 0.0
t = 21 x = [ 0.  0.] score = 0.0
t = 22 x = [ 0.  0.] score = 0.0
t = 23 x = [ 0.  0.] score = 0.0
t = 24 x = [ 0.  0.] score = 0.0
t = 25 x = [ 0.  0.] score = 0.0
t = 26 x = [ 0.  0.] score = 0.0
t = 27 x = [ 0.  0.] score = 0.0
t = 28 x = [ 0.  0.] score = 0.0
t = 29 x = [ 0.  0.] score = 0.0
t = 29 x = [ 0.  0.] score = 0.0

アニメーション表示すると

こんな感じ。赤い点がその時点でのグローバルベスト。 各粒子がグローバルベストを中心に探索しているのが見て取れる。

f:id:magayengineer:20160407233313g:plain

ソースコード全文

#!/usr/bin/python
# -*- coding: utf-8 -*-
import numpy as np
import math
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import matplotlib.animation as animation


class BaseParticle:
    """
    PSOの基礎クラス
    """

    def __init__(self, w, x, v, c1, c2):
        """
        コンストラクタ
        :param w: 慣性定数
        :param x: 現在地(numpy.ndarray)
        :param v: 速度  (numpy.ndarray)
        :param c1:群のうちで良い位置に向かう粒子の割合
        :param c2:群のうちで良い位置に向かう粒子の割合
        :return 粒子のインスタンス
        """
        self._x = x
        self._v = v
        self._w = w
        self._c1 = c1
        self._c2 = c2

        # パーソナルベストを初期化
        self._p_best_x = x.copy()
        self._p_best_score = self.eval()

    def update(self, g_best_x):
        """
        粒子の更新を行う
        :param g_best_x:グローバルベスト
        :return:void
        """
        # 粒子の位置を更新
        self._x = self._x + self._v

        # 粒子の位置を補正
        self.clip_x()

        # 粒子の速度を更新
        r1 = np.random.rand()
        r2 = np.random.rand()
        # v = wv + c1r1(px - x) * c2r2 (gp - x)
        self._v = self._w * self._v + self._c1 * r1 * (self._p_best_x - self._x) + self._c2 * r2 * (g_best_x - self._x)

        # 速度を補正
        self.clip_v()

        # 評価
        score = self.eval()

        # パーソナルベスト更新
        if not self.is_better_than(score):
            self._p_best_score = score
            self._p_best_x = self._x.copy()

    @property
    def x(self):
        """
        :return: 現在の位置
        """
        return self._x

    @property
    def p_best_x(self):
        """
        :return: パーソナルベストの位置
        """
        return self._p_best_x

    @property
    def p_best_score(self):
        """
        :return:パーソナルベストの値
        """
        return self._p_best_score

    def eval(self):
        """
        評価。ここに評価関数を記述する。
        :return:評価値
        """
        pass

    @classmethod
    def create_particle(cls, w, x, v, c1, c2):
        """
        粒子作成用のファクトリメソッド
        :param w: 慣性定数
        :param x: 現在地(numpy.ndarray)
        :param v: 速度  (numpy.ndarray)
        :param c1:群のうちで良い位置に向かう粒子の割合
        :param c2:群のうちで良い位置に向かう粒子の割合
        :return 粒子のインスタンス
        """
        pass

    def is_better_than(self, score):
        """
        パーソナルベストが、引数scoreより良いか判定
        :param score:比較対象の値
        :return:パーソナルベストのほうが良ければtrue
        """
        pass

    def clip_x(self):
        """
        位置を補正する
        :return: void
        """
        pass

    def clip_v(self):
        """
        速度を補正する
        :return: void
        """
        pass

    @staticmethod
    def show_animation(animation_data):
        """
        最適化の過程をアニメーション表示する
        :param animation_data: アニメーション用データ
        :return: void
        """
        pass


class XxYyParticle(BaseParticle):
    """
    Z = X * X + Y * Yの最適化を行うための粒子クラス
    """
    w = 0.9
    c1 = 0.9
    c2 = 0.9

    def eval(self):
        x = self._x
        return x[0] * x[0] + x[1] * x[1]

    @classmethod
    def create_particle(cls):

        x = np.array([200 * np.random.rand() - 100, 200 * np.random.rand() - 100])
        v = np.array([20 * np.random.rand() - 10, 20 * np.random.rand() - 10])
        return XxYyParticle(cls.w, x, v, cls.c1, cls.c2)

    def is_better_than(self, score):
        if self._p_best_score < score:
            return True
        else:
            return False

    def clip_x(self):
        self._x = np.clip(self._x, -100, 100)

    def clip_v(self):
        self._v = np.clip(self._v, -10, 10)

    @staticmethod
    def show_animation(graph_data):
        # 関数の外観を書く ######################################
        x = np.arange(-100, 100, 10)
        y = np.arange(-100, 100, 10)

        mesh_x, mesh_y = np.meshgrid(x, y)
        mesh_z = mesh_x * mesh_x + mesh_y * mesh_y

        show_3d_animation(mesh_x, mesh_y, mesh_z, graph_data)


class MatyasParticle(BaseParticle):
    """
    Matyasの最適化を行うための粒子クラス
    """
    w = 0.9
    c1 = 0.9
    c2 = 0.9

    def eval(self):
        x = self._x
        return 0.26*(x[0] * x[0] + x[1] * x[1])-0.48*x[0]*x[1]

    @classmethod
    def create_particle(cls):

        x = np.array([20 * np.random.rand() - 10, 20 * np.random.rand() - 10])
        v = np.array([2 * np.random.rand() - 1, 2 * np.random.rand() - 1])
        return MatyasParticle(cls.w, x, v, cls.c1, cls.c2)

    def is_better_than(self, score):
        if self._p_best_score < score:
            return True
        else:
            return False

    def clip_x(self):
        self._x = np.clip(self._x, -10, 10)

    def clip_v(self):
        self._v = np.clip(self._v, -10, 10)

    @staticmethod
    def show_animation(graph_data):
        # 関数の外観を書く ######################################
        x = np.arange(-15, 15, 0.5)
        y = np.arange(-15, 15, 0.5)

        mesh_x, mesh_y = np.meshgrid(x, y)
        mesh_z = 0.26 * (mesh_x * mesh_x + mesh_y * mesh_y) - 0.48 * mesh_x * mesh_y
        show_3d_animation(mesh_x, mesh_y, mesh_z, graph_data)


class LeviFunctionParticle(BaseParticle):
    """
    LeviFunction関数
    """
    w = 0.9
    c1 = 0.9
    c2 = 0.9

    def eval(self):
        x = self._x
        x1 = x[0]
        x2 = x[1]
        a = math.pow(np.sin(3 * np.pi * x1), 2)
        b = math.pow(x1 - 1, 2) * (1 + math.pow(np.sin(3 * np.pi * x2), 2))
        c = math.pow(x2 - 1, 2) * (1 + math.pow(np.sin(2 * np.pi * x2), 2))
        return a + b + c

    @classmethod
    def create_particle(cls):
        x = np.array([20 * np.random.rand() - 10, 20 * np.random.rand() - 10])
        v = np.array([20 * np.random.rand() - 10, 20 * np.random.rand() - 10])
        return LeviFunctionParticle(cls.w, x, v, cls.c1, cls.c2)

    def is_better_than(self, score):
        if self._p_best_score < score:
            return True
        else:
            return False

    def clip_x(self):
        self._x = np.clip(self._x, -10, 10)

    def clip_v(self):
        self._v = np.clip(self._v, -10, 10)

    @staticmethod
    def show_animation(graph_data):
        # 関数の外観を書く ######################################
        x = np.arange(-10, 10, 0.25)
        y = np.arange(-10, 10, 0.25)

        mesh_x, mesh_y = np.meshgrid(x, y)
        a = np.power(np.sin(3 * np.pi * mesh_x), 2)
        b = np.power(mesh_x - 1, 2) * (1 + np.power(np.sin(3*np.pi*mesh_y), 2))
        c = np.power(mesh_y - 1, 2) * (1 + np.power(np.sin(2*np.pi*mesh_x), 2))
        mesh_z = a + b + c

        show_3d_animation(mesh_x, mesh_y, mesh_z, graph_data)


class RastriginParticle(BaseParticle):
    """
    Rastrigin関数
    """

    w = 0.9
    c1 = 0.9
    c2 = 0.9

    def eval(self):
        x = self._x
        x1 = x[0]
        x2 = x[1]
        return 10 * 2 + (x1 * x1 - 10 * np.cos(2 * np.pi * x1)) + (x2 * x2 - 10 * np.cos(2 * np.pi * x2))

    @classmethod
    def create_particle(cls):
        x = np.array([5.12 * 2 * np.random.rand() - 5.12, 5.12 * 2 * np.random.rand() - 5.12])
        v = np.array([5.12 * 2 * np.random.rand() - 5.12, 5.12 * 2 * np.random.rand() - 5.12])
        return RastriginParticle(cls.w, x, v, cls.c1, cls.c2)

    def is_better_than(self, score):
        if self._p_best_score < score:
            return True
        else:
            return False

    def clip_x(self):
        self._x = np.clip(self._x, -5.12, 5.12)

    def clip_v(self):
        self._v = np.clip(self._v, -5.12, 5.12)

    @staticmethod
    def show_animation(graph_data):
        # 関数の外観を書く ######################################
        x = np.arange(-5.12, 5.12, 0.1)
        y = np.arange(-5.12, 5.12, 0.1)

        mesh_x, mesh_y = np.meshgrid(x, y)
        a = mesh_x * mesh_x - 10 * np.cos(2 * np.pi * mesh_x)
        b = mesh_y * mesh_y - 10 * np.cos(2 * np.pi * mesh_y)
        mesh_z = 10 * 2 + a + b

        show_3d_animation(mesh_x, mesh_y, mesh_z, graph_data)


class AnimationData:
    """
    Animationデータの保持クラス
    """
    def __init__(self):
        """
        コンストラクタ
        :return:
        """
        # 粒子の時系列ごとの座標データ
        self._particles_position_t = []
        # グローバルベストの時系列座標データ
        self._best_particle_position_t = None

    def add_particle_position_data(self, particles):
        """
        時刻tの粒子群の座標をデータに追加する
        :param particles: 時刻tの粒子データ
        :return: void
        """
        ptx = None
        for n, particle in enumerate(particles):
            if n == 0:
                ptx = np.array([particle.x.copy()])
                ptx = np.append(ptx, particle.eval())
            else:
                tmp = np.array([particle.x.copy()])
                tmp = np.append(tmp, particle.eval())
                ptx = np.vstack((ptx, tmp))
        self._particles_position_t.append(ptx)

    def add_best_particle_position_data(self, global_best_position, global_best_score):
        """
        時刻tの粒子群の座標をデータに追加する
        :param global_best_position:時刻tのグローバルベストの座標
        :param global_best_score:時刻tのグローバルベストの値
        :return:void
        """
        if self._best_particle_position_t is None:
            self._best_particle_position_t = np.array([global_best_position.copy()])
            self._best_particle_position_t = np.append(self._best_particle_position_t, global_best_score)
        else:
            tmp = np.array([global_best_position.copy()])
            tmp = np.append(tmp, global_best_score)
            self._best_particle_position_t = np.vstack((self._best_particle_position_t, tmp))

    @property
    def global_best_position(self):
        return self._best_particle_position_t

    @property
    def particle_position_t(self):
        return self._particles_position_t


def show_3d_animation(mesh_x, mesh_y, mesh_z, animation_data):
    """
    最適化の過程をアニメ表示
    :param mesh_x: グラフ描画用のメッシュデータ
    :param mesh_y: グラフ描画用のメッシュデータ
    :param mesh_z: グラフ描画用のメッシュデータ
    :param animation_data:粒子の座標データ(アニメーション用)
    :return:void
    """
    fig = plt.figure()
    ax = Axes3D(fig)
    ax.plot_wireframe(mesh_x, mesh_y, mesh_z, color='#B8F28C')
    ax.view_init(elev=75, azim=-50) #カメラ位置調整

    # アニメーションデータ作成

    ims = []
    gbp = animation_data.global_best_position
    for i in range(len(animation_data.particle_position_t)):
        pts = animation_data.particle_position_t[i]
        tim = ax.scatter3D(pts[:, 0], pts[:, 1], pts[:, 2],  c='black')
        gim = ax.scatter3D(gbp[i, 0], gbp[i, 1], gbp[i, 2], s=100, c='red')
        ims.append([tim, gim])

    # アニメーション作成 ※なぜか左辺を作らなきゃ動かないので注意
    ani = animation.ArtistAnimation(fig, ims, interval=1000, repeat_delay=10000)

    # グラフ表示
    plt.show()

    # アニメーション保存
    #ani.save("test.mp4")


def calc(particle_class, times, num_of_particle, show_animation=False):
    """
    計算実行
    :param particle_class: 粒子のクラス
    :param times: 何回試行するかするか
    :param num_of_particle: 粒子数
    :param show_animation:グラフ表示有無
    :return: void
    """

    # グローバルベスト
    global_best_x = None
    global_best_score = None

    # 粒子郡
    particles = []

    if show_animation:
        animation_data = AnimationData()

    # 粒子初期化
    for n in range(num_of_particle):

        particle = particle_class.create_particle()
        particles.append(particle)

        if n == 0:
            global_best_x = particle.p_best_x
            global_best_score = particle.p_best_score
        else:
            if particle.is_better_than(global_best_score):
                global_best_score = particle.p_best_score
                global_best_x = particle.p_best_x

    # アニメーション表示用データ準備
    if show_animation:
        animation_data.add_best_particle_position_data(global_best_x, global_best_score)
        animation_data.add_particle_position_data(particles)

    for t in range(times):
        # 粒子更新
        for n, particle in enumerate(particles):
            particle.update(global_best_x)
            if particle.is_better_than(global_best_score):
                global_best_score = particle.p_best_score
                global_best_x = particle.p_best_x

        # アニメーション表示用データ準備
        if show_animation:
            animation_data.add_best_particle_position_data(global_best_x, global_best_score)
            animation_data.add_particle_position_data(particles)

        print "t =", t, "x =", global_best_x, "score =", global_best_score
    print "t =", t, "x =", global_best_x, "score =", global_best_score

    # graph表示
    if show_animation:
        particle_class.show_animation(animation_data)

if __name__ == "__main__":
    #calc(XxYyParticle, 100, 30, True)
    #calc(MatyasParticle, 10, 30, True)
    #calc(LeviFunctionParticle, 30, 10, True)
    calc(RastriginParticle, 30, 30, True)

参考サイト

https://ja.wikipedia.org/wiki/%E7%B2%92%E5%AD%90%E7%BE%A4%E6%9C%80%E9%81%A9%E5%8C%96 http://qiita.com/tomitomi3/items/d4318bf7afbc1c835dda