DEAP
目次
はじめに
Pythonで遺伝的アルゴリズムのライブラリーが無いか検索してみたところ、どうもDEAPというライブラリーが良いらしいので使い方を調べてみることにしました。
PyPIにも登録されているので、pipでインストールできます。
pip install -U deap
DEAPは遺伝的アルゴリズムだけでなく、他にも色々な進化的計算のアルゴリズムが実装されているみたいです。
遺伝的プログラミング?。というのもできるらしく、よくわからないですがなんか凄そうです。
とりあえず今回の目標は遺伝的アルゴリズムの使い方を覚えることとして、また後で気が向いたら他にも調べてみようと思います。
使い方の解説
DEAPのドキュメントでOneMax問題(一つある目的関数を最大化する問題)のチュートリアルがあるので、使い方を解説していこうと思います。
設計変数である遺伝子は全部0か1のビット値を取り、目的関数は個体内の遺伝子の総和としています。
実は、前回の遺伝的アルゴリズムのサンプルコードはまさにこのチュートリアルコードの内容と一緒です。
チュートリアルコードは、GitHub上のexamples/ga/onemax.pyにあります。
下記のコードは、チュートリアルコードのコメントを全部取っ払って(表示すると長くなるので)私のGistにコピーしたものです。
GAの実装に関して必要な変数と関数を定義するのに、以下2つの重要な関数がDEAPで用意されています。
- creator.create()
- base.Toolbox.register()
この2つの関数の使い方についてみていきましょう。
creator.create()関数
遺伝子(設計変数)のセットを表す個体を定義するのに、creator
というモジュールが用意されており、複雑な構造の個体を表現できるみたいです。
例えば、遺伝子のセットを保存したい場合、list
を取りあえず使えばいいと思いますが、遺伝子が保存されたlist
は個体を表しているので、適応度も一緒に保存出来たらいいな。と思ってしまいます。
creator.create()
関数を使うと、list
に適応度(fitness)を保存するメンバ変数を、簡単に追加することができます。
creator.create()
関数は、あるクラスにメンバ変数を追加して子クラスを新たに作成する関数です。
creator.create()
関数を使う時は、少なくとも2この引数を指定する必要があります。
第2引数で指定したクラスを継承して、第1引数で指定する名前のクラスをcreator
モジュール内に新しく生成します。
第3引数以降は、生成する子クラスに追加するメンバ変数を指定します。
適応度の定義
7行目の記述は、base.Fitness
クラスを継承して、weights=(1.0,)
というメンバ変数を追加した適応度を表すFitnessMax
というクラスをcreator
モジュール内に作成しています。
creator.create("FitnessMax", base.Fitness, weights=(1.0,))
weights=(1.0,)
は、評価する目的関数が1つの場合で、個体の適応度を最大化したい場合にこう書きます。1.0
の後に,
がついていることに注意してください。
適応度を最小化したい場合は、マイナスを付けてweights=(-1.0,)
と書きます。
また、評価したい目的関数が2つある場合は、1つの個体につき2つの適応度を評価する必要があります。
例えば2つある適応度の1つを最大化、もう1つを最小化する場合は、weights=(1.0,-1.0)
と書いてやります。
個体の定義
8行目の記述は、list
クラスを継承して、fitness=creator.FitnessMax
というメンバ変数を追加したIndividual
クラスを作成しています。
creator.create("Individual", list, fitness=creator.FitnessMax)
遺伝子を保存する個体list
に、適応度fitness
というメンバ変数を追加してIndividual
としているわけです。
base.Toolbox.register()関数
base.Toolbox.register()
関数を使うと、引数のデフォルト値が無い関数に、デフォルト値を設定することができます。
第2引数で指定する関数に、第3引数以降で指定するデフォルト値を設定して、第1引数で指定する名前でbase.Toolbox
内に新しく関数を作成します。
例えば以下の様に、デフォルト値がない引数a
、b
、c
をもつ関数test
を、引数なしで実行するとエラーになってしまいますが、
def test(a,b,c): print a,b,c test()
以下の様に、base.Toolbox.register()関数を使って引数のデフォルト値を指定してやると、引数を指定せずに実行できる関数を作成することができます。
toolbox = base.Toolbox() toolbox.register("test_wda", test, 1,2,3) toolbox.test_wda()
チュートリアルコードのtoolbox.register()
関数を使っている部分を一つずつ見ていきましょう。
遺伝子を作成する関数
12行目の記述は、ランダムに整数を生成する関数random.randint()
のデフォルト引数に0
と1
を指定してattr_bool
という名前の関数を作成しています。
toolbox.register("attr_bool", random.randint, 0, 1)
デフォルト引数に設定した0
と1
ですが、ランダムに生成する整数の範囲(minとmax値)を表しており、0
か1
をランダムで返す引数なしで実行できる関数attr_bool
が出来上がります。
個体を作成する関数
13行目の記述では、よくわからない関数tools.initRepeat
をindividual
によみかえていますね。
toolbox.register("individual", tools.initRepeat, creator.Individual, toolbox.attr_bool, 100)
tools.initRepeat
の定義は、deap/tools/init.py
に書かれています。
def initRepeat(container, func, n): return container(func() for _ in xrange(n))
tools.initRepeat
は、3この引数container
、func
、n
を持っており、n
回func
を実行して、その値をcontainer
に格納して返す関数です。
つまり、100
回toolbox.attr_bool
を実行し、その値をcreator.Individual
に格納して返す関数individual
を作成しています。引数なしで個体を生成できるindividual
関数ができたわけです。
世代を作成する関数
14行目の記述も、tools.initRepeat
を読みかえていますね。
toolbox.register("population", tools.initRepeat, list, toolbox.individual)
個体をtoolbox.individual
で作成してlist
に格納し、世代を生成するpopulation
関数を作っています。
世代内の個体数あたるn
は、後でmain
関数内で指定するので、デフォルト値は設定してないでおきます。
attr_bool
で遺伝子を生成し、individual
で個体を生成、population
で世代を生成するというわけです。
19~22行目は、目的関数、また交差・突然変異・選択を求める関数を読みかえています。
toolbox.register("evaluate", evalOneMax) toolbox.register("mate", tools.cxTwoPoint) toolbox.register("mutate", tools.mutFlipBit, indpb=0.05) toolbox.register("select", tools.selTournament, tournsize=3)
それぞれの関数の定義を見ていきましょう。
目的関数
16~17行目では、目的関数を指定していますね。
def evalOneMax(individual): return sum(individual),
個体individual
の遺伝子の総和をかえしています。return sum(individual),
ですが、最後に,
を付けることに注意してください。
交叉関数
DEAPのソースコードのdeap/tools/crossover.py
に定義があります。
def cxTwoPoint(ind1, ind2): size = min(len(ind1), len(ind2)) cxpoint1 = random.randint(1, size) cxpoint2 = random.randint(1, size - 1) if cxpoint2 >= cxpoint1: cxpoint2 += 1 else: # Swap the two cx points cxpoint1, cxpoint2 = cxpoint2, cxpoint1 ind1[cxpoint1:cxpoint2], ind2[cxpoint1:cxpoint2] \ = ind2[cxpoint1:cxpoint2], ind1[cxpoint1:cxpoint2] return ind1, ind2
二つの個体ind1
, ind2
を引数として、固体内の遺伝子を2点交叉で入れ替えて返す関数です。
突然変異関数
deap/tools/mutation.py
に定義があります。
def mutFlipBit(individual, indpb): for i in xrange(len(individual)): if random.random() < indpb: individual[i] = type(individual[i])(not individual[i]) return individual,
一つの個体individual
と遺伝子突然変異率indpb
を引数にして、固体内の遺伝子を突然変異させて返していますね。
遺伝子が0
と1
のビット値なので、ビット値を入れ替えることで突然変異を表しています。
21行目の読み替え記述では、indpb=0.05
というデフォルト値を与えています。遺伝子突然変異率は5%で決め打ちという事です。
選択関数
deap/tools/selection.py
に定義があります。
def selTournament(individuals, k, tournsize, fit_attr="fitness"): chosen = [] for i in xrange(k): aspirants = selRandom(individuals, tournsize) chosen.append(max(aspirants, key=attrgetter(fit_attr))) return chosen
トーナメント選択を実施する関数です。引数はindividuals
, k
, tournsize
, fit_attr="fitness"
の4つですね。
individuals
は、個体individual
の後ろにs
がついており複数形になっています。紛らわしいですが、複数ある個体は世代を表すので、この引数には世代(population
)を入れてやります。
k
は、世代individuals
中の個体の数です。tournsize
はトーナメントサイズです。22行目の読み替え時は、tournsize=3
で決め打ちしています。
fit_attr="fitness"
は、固体のメンバ変数になっている適応度の変数名を指定しています。fit_attr
はfitness_attribute
(適応度_属性)の略だと思います。
デフォルト値がfitness
になってますね。8行目の個体定義時に、fitness
という名前以外で適応度を登録した場合、ここで変えないとエラーになるわけですね。
ちなみに最近知りましたが、Pythonではメンバ変数の事を属性(attribute)とも呼ぶらしいです。
注意点
GA実装時は、目的関数は自分で作成する必要があり、また交叉、突然変異、選択の関数は、DEAPに適当なものがなければ、自分で作成する必要があります。
目的関数と、突然変異関数の戻り値の最後に,
を付ける必要があるので注意してください。この,
ですが、戻り値をtupleにして返す場合に,
を付けるみたいです。
main関数内の処理
続いてmain関数内の処理を見ていきましょう。
28行目で初期世代を作成しています。ここでn=300
で世代内の個体数を指定しています。
pop = toolbox.population(n=300)
toolbox.population
を実行するだけで、遺伝子、固体、世代を一気に作れるというわけです。
29行目の記述は、交叉率、個体突然変異率、ループを回す世代数を指定しています。
CXPB, MUTPB, NGEN = 0.5, 0.2, 40
33~35行目の記述は、初期世代の個体の適応度を計算し、メンバ変数fitness
に登録しています。個体内の適応度にはfitness.values
でアクセスします。
fitnesses = list(map(toolbox.evaluate, pop)) for ind, fit in zip(pop, fitnesses): ind.fitness.values = fit
42、43行目は、選択の操作ですね。42行目で選択した個体をoffspring
に格納してから、43行目でクローンを作ってまたoffspring
に入れてますね。
offspring = toolbox.select(pop, len(pop)) offspring = list(map(toolbox.clone, offspring))
toolbox.clone
の実体はcopyモジュールのdeepcopy
です。これ以降の記述で、個体の適応度fitness
の値をすべて参照渡ししているために、ここで参照を切っているんだと思いますが、ちょっと複雑で具体的な挙動は追えませんでした。43行目をコメントアウトすると、最適化がすすまなくなってしまうので、まあ気にせず入れておきましょう。
45~50行目は、交叉の処理ですね。
for child1, child2 in zip(offspring[::2], offspring[1::2]): if random.random() < CXPB: toolbox.mate(child1, child2) del child1.fitness.values del child2.fitness.values
offspring
の偶数インデックスoffspring[::2]
と奇数インデックスoffspring[1::2]
でそれぞれの個体を交叉させてます。
zip
ってあまり使ったことがないですが、こうしてみると結構使えますね。
交叉で遺伝子を取り換えた後は、適応度を再計算する必要があるためdel
で削除してます。
52~56行目は、突然変異です。
for mutant in offspring: if random.random() < MUTPB: toolbox.mutate(mutant) del mutant.fitness.values
個体突然変異率MUTPB
で突然変異する個体をoffspring
の中から選び、あとは、toolbox.mutate
関数に突っ込んで、固体内の遺伝子を突然変異させてます。
こちらも突然変異が発生した個体の適応度を再計算する必要があるためdel
で削除してます。
58~61行目は、交叉と突然変異でdel
した適応度を再計算しています。
invalid_ind = [ind for ind in offspring if not ind.fitness.valid] fitnesses = map(toolbox.evaluate, invalid_ind) for ind, fit in zip(invalid_ind, fitnesses): ind.fitness.values = fit
リスト内包表記の仕方がうまいですね。if not ind.fitness.valid
を選んで適応度を再計算させると。こんな書き方できるんですね。
ここで一連のステップが完了するので、これ以降の記述は次世代offspring
を現世代pop
にコピー(65行目のpop[:] = offspring
)してまたループするという流れになってます。
個体にnumpyのndarrayを使う場合
上記のチュートリアルでは、個体の遺伝子保存コンテナにlist
を使っており、その後の処理もlist
を想定して書かれています。
list
ではなく、numpy
のndarray
を使う方が目的関数内で複雑な計算する場合は計算が楽になるので、できればndarray
を使いたいなぁ。と思ってたらしっかりチュートリアルが用意されていました。
こちらドキュメントとチュートリアルコードのリンクです。
下記のコードはチュートリアルコードのコメントを削除して、私のGistにコピーしたものです。
10行目の個体Individual
の定義時に、継承元がlist
ではなくnumpy.ndarray
になってます。
creator.create("Individual", numpy.ndarray, fitness=creator.FitnessMax)
個体にnumpyのndarrayを使う場合は、交叉関数を少し書き換える必要があります。
30行目で個体内の遺伝子を入れ替える時に、.copy()
を使っていることに注意してください。
ind1[cxpoint1:cxpoint2], ind2[cxpoint1:cxpoint2] = ind2[cxpoint1:cxpoint2].copy(), ind1[cxpoint1:cxpoint2].copy()
こちらのチュートリアルコードでは、main関数内の処理が随分すっきりしちゃいましたね。世代ループをalgorithms.eaSimple
という関数内に突っ込んじゃってます。(53行目です。)
algorithms.eaSimple(pop, toolbox, cxpb=0.5, mutpb=0.2, ngen=40, stats=stats,halloffame=hof)
algorithms.eaSimple
の定義は、deap/algorithms.py内にあります。
このalgorithms.eaSimple
を使って世代ループを回す際に、渡してやる引数halloffame=hof
も、numpy用に書き換える必要があります。(45行目の記述です。)
hof = tools.HallOfFame(1, similar=numpy.array_equal)
tools.HallOfFame
は、世代内のベスト解(適応度が最大の個体)を保存する為のクラスで、deap/tools/support.pyに定義があります。
ちなみにこのHall of fame
をGoogle翻訳すると栄誉殿堂だそうです。ベスト解を殿堂入りさせているんですね。
遺伝子が連続値の場合
チュートリアルでは設計変数である遺伝子を0と1のビット値としていましたが、実際は実数の連続値になるケースが多いと思います。
以下遺伝子が連続値になる様にnumpy版のチュートリアルコードを書き換えてしてみました。
個体内の遺伝子の個数は100個とし、遺伝子は-1.0~1.0の範囲を取る実数値です。
目的関数は遺伝子のsum値で、個体内の遺伝子が全て1.0になる場合に、最適解100が得られます。
変更点
14行目~16行目で、遺伝子の取り得る範囲を指定しています。個体内の遺伝子の個数がn_gene = 100
で、取り得る範囲のmin値、max値をmin_ind
、max_ind
で設定しています。
n_gene = 100 min_ind = numpy.ones(n_gene) * -1.0 max_ind = numpy.ones(n_gene) * 1.0
18行目~22行目で、個体を生成する関数を定義しています。
def create_ind_uniform(min_ind, max_ind): ind = [] for min, max in zip(min_ind, max_ind): ind.append(random.uniform(min, max)) return ind
引数に遺伝子の取り得る範囲min_ind
、max_ind
を渡してやり、個体をlist
で返します。遺伝子の生成は、random.uniform
で一様乱数を生成しています。
24行目の記述は、本来は遺伝子を作成する関数を設定しますが、引数なしでlist
で個体を返す関数にしています。
toolbox.register("create_ind", create_ind_uniform, min_ind, max_ind)
というのも、取り得る範囲を各遺伝子で設定したい場合、どうしても範囲を設定したい遺伝子のインデックスが必要になるので、デフォルト引数を各遺伝子で設定するのが難しいです。引数なしで遺伝子作成ができる関数を作れません。
こういう場合は、引数なしでlist
で個体を作成する関数を作っておき、適応度をメンバ変数に持つIndividual
に渡してやります。
25行目の記述でそれを行っています。
toolbox.register("individual", tools.initIterate, creator.Individual, toolbox.create_ind)
tools.initIterate
関数ですが、定義はdeap/tools/init.pyにあります。
def initIterate(container, generator): return container(generator())
個体を生成するgenerator
関数を、個体の保存先container
につっこんでいます。
また、突然変異関数も変更する必要があります。44~49行目に定義しています。
突然変異する遺伝子の取り得る範囲を指定して一様乱数で変異させます。
def mutUniformDbl(individual, min_ind, max_ind, indpb): size = len(individual) for i, min, max in zip(xrange(size), min_ind, max_ind): if random.random() < indpb: individual[i] = random.uniform(min, max) return individual,
実行結果
以下、1000回の世代ループを回した結果です。
gen nevals avg std min max 0 300 0.0468393 5.64126 -12.315 17.335 1 233 5.20779 5.1087 -10.5526 20.6311 2 219 9.60382 4.31149 -1.569 20.6311 3 216 13.0049 3.59862 -1.23289 24.8766 ... ... ... 997 225 93.8435 3.21631 83.95 96.8 998 232 93.6538 3.32555 81.7815 96.8 999 224 93.6433 3.1658 83.4545 96.8 1000 222 93.4084 3.3193 83.9801 96.8
max
が世代内のベスト解です。17.335から始まって、1000回目のループで、96.8まで最適化できました。
やはり実数となると遺伝子の取り得る値が拡張されるので、かなりのループを回さないと中々最適解が得られないですね。
最後の方になるとベスト解が更新されなくなってきます。交叉率や突然変異率を調整する必要があるみたいです。
遺伝子の個数を100個で行いましたが、さらに遺伝子の個数が増えた場合や、取り得る範囲をもっと広くとった場合は結構時間がかかるかもですね。