Pythonと機械学習

Pythonも機械学習も初心者ですが、頑張ってこのブログで勉強してこうと思います。

トップ > Python > DEAP

DEAP

Python 最適化

目次

はじめに

Pythonで遺伝的アルゴリズムのライブラリーが無いか検索してみたところ、どうもDEAPというライブラリーが良いらしいので使い方を調べてみることにしました。

ソースコードはGitHub上にあげられています。

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個で行いましたが、さらに遺伝子の個数が増えた場合や、取り得る範囲をもっと広くとった場合は結構時間がかかるかもですね。