【Python】DEAPを用いた遺伝的アルゴリズムの実装

## DEAPの概要と基本トピック
### DEAPの概要
DEAP(Distributed Evolutionary Algorithms in Python)は`Python`で実装された進化的計算(Evolutionary Computation)のライブラリです。

https://en.wikipedia.org/wiki/DEAP_(software)


### DEAPのインストールとサンプルコードの実行
DEAPは下記のコマンドを実行することでPyPIからダウンロードすることができます。

```
$ pip install deap
```

https://deap.readthedocs.io/en/master/installation.html

以下、上記に基づいてDEAPのサンプルコードの実行と確認を行います。

https://deap.readthedocs.io/en/master/overview.html

```python
from deap import base, creator
import random
from deap import tools


creator.create("FitnessMin", base.Fitness, weights=(-1.0,))
creator.create("Individual", list, fitness=creator.FitnessMin)

IND_SIZE = 10
toolbox = base.Toolbox()

toolbox = base.Toolbox()
toolbox.register("attribute", random.random)
toolbox.register("individual", tools.initRepeat, creator.Individual,
                 toolbox.attribute, n=IND_SIZE)
toolbox.register("population", tools.initRepeat, list, toolbox.individual)


def evaluate(individual):
    return sum(individual),

toolbox.register("mate", tools.cxTwoPoint)
toolbox.register("mutate", tools.mutGaussian, mu=0, sigma=1, indpb=0.1)
toolbox.register("select", tools.selTournament, tournsize=10)
toolbox.register("evaluate", evaluate)


def main():
    pop = toolbox.population(n=2)
    CXPB, MUTPB, NGEN = 0.5, 0.5, 10

    # Evaluate the entire population
    fitnesses = map(toolbox.evaluate, pop)
    for ind, fit in zip(pop, fitnesses):
        ind.fitness.values = fit

    for g in range(NGEN):
        # Select the next generation individuals
        offspring = toolbox.select(pop, len(pop))
        print(evaluate(offspring[0]), evaluate(offspring[1]))
        # Clone the selected individuals
        offspring = list(map(toolbox.clone, offspring))

        # Apply crossover and mutation on the offspring
        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

        for mutant in offspring:
            if random.random() < MUTPB:
                toolbox.mutate(mutant)
                del mutant.fitness.values

        # Evaluate the individuals with an invalid fitness
        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

        # The population is entirely replaced by the offspring
        pop[:] = offspring
        print(evaluate(pop[0]), evaluate(pop[1]))
        print("===")

    return pop


if __name__ == "__main__":
    random.seed(0)
    results = main()
    print(len(results))
    print(results[0])
```

・実行結果
```
(5.345163753129269,) (5.345163753129269,)
(4.8377211732155745,) (5.345163753129269,)
===
(4.8377211732155745,) (4.8377211732155745,)
(4.8377211732155745,) (4.8377211732155745,)
===
(4.8377211732155745,) (4.8377211732155745,)
(6.373778345485883,) (4.754760051599324,)
===
(4.754760051599324,) (4.754760051599324,)
(4.754760051599324,) (4.754760051599324,)
===
(4.754760051599324,) (4.754760051599324,)
(4.754760051599324,) (3.6226515831862756,)
===
(3.6226515831862756,) (3.6226515831862756,)
(3.063431427375879,) (3.6226515831862756,)
===
(3.063431427375879,) (3.063431427375879,)
(3.063431427375879,) (3.063431427375879,)
===
(3.063431427375879,) (3.063431427375879,)
(3.063431427375879,) (2.286771005230591,)
===
(2.286771005230591,) (2.286771005230591,)
(3.598490829953151,) (2.286771005230591,)
===
(2.286771005230591,) (2.286771005230591,)
(2.286771005230591,) (2.286771005230591,)
===
2
[1.2998733298460659, 0.7579544029403025, 0.420571580830845, 0.25891675029296335, 1.0522301608532234, -1.2026840377582917, 0.7837985890347726, 0.30331272607892745, -0.6555115142606929, -0.7316909826275244]
```

プログラムの`main`関数の主な流れは下記を元に把握すると良いです。

1) 個体の選択(Selection)
→ `toolbox.select`を用いて個体の選択を行う。選択のアルゴリズムは`toolbox.register("select", tools.selTournament, tournsize=10)`で指定しており、`selTournament`の場合は`tournsize=10`だけ個体の重複サンプリングを行い、その中で最適なサンプルを選択する。

2) 交叉(Crossover)
→ `toolbox.mate`を用いて交叉の処理を行う。アルゴリズムの詳細は`toolbox.register("mate", tools.cxTwoPoint)`で指定する。

3) 突然変異(Mutation)
→ `toolbox.mutate`を用いて交叉の処理を行う。アルゴリズムの詳細は`toolbox.register("mutate", tools.mutGaussian, mu=0, sigma=1, indpb=0.1)`で指定する。


また、`select`で選択を行うにあたっての基準は`fitnesses = map(toolbox.evaluate, invalid_ind)`によって計算されます。詳細は`toolbox.register("evaluate", evaluate)`で指定されており、`evaluate`関数を参照します。


### DEAPのオペレータ(Operators)
![DEAP1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3936087/9088134c-33d6-48f3-b88f-9da819afb84e.png)
[DEAPのドキュメント](https://deap.readthedocs.io/en/master/api/tools.html) より

DEAPにおける処理は上記のオペレータを指定することで実行が可能です。上記には前項で用いた`cxTwoPoint`(Crossover)や`mutGaussian`(Mutation)、`selTournament`(Selection)などが含まれています。



## DEAPの応用例
### 巡回セールスマン問題(TSP)
https://github.com/DEAP/deap/blob/master/examples/ga/tsp.py

当項では以下、上記を参考にDEAPにおける巡回セールスマン問題の最適解の近似について確認します。

```python
import array
import random
import json

import numpy as np

from deap import algorithms
from deap import base
from deap import creator
from deap import tools

n_pop, tourn_size, n_iter = 500, 20, 100


# gr*.json contains the distance map in list of list style in JSON format
# Optimal solutions are : gr17 = 2085, gr24 = 1272, gr120 = 6942
with open("tsp/gr120.json", "r") as tsp_data:
    tsp = json.load(tsp_data)

distance_map = tsp["DistanceMatrix"]
IND_SIZE = tsp["TourSize"]

creator.create("FitnessMin", base.Fitness, weights=(-1.0,))
creator.create("Individual", array.array, typecode='i', fitness=creator.FitnessMin)

toolbox = base.Toolbox()

# Attribute generator
toolbox.register("indices", random.sample, range(IND_SIZE), IND_SIZE)

# Structure initializers
toolbox.register("individual", tools.initIterate, creator.Individual, toolbox.indices)
toolbox.register("population", tools.initRepeat, list, toolbox.individual)

def evalTSP(individual):
    distance = distance_map[individual[-1]][individual[0]]
    for gene1, gene2 in zip(individual[0:-1], individual[1:]):
        distance += distance_map[gene1][gene2]
    return distance,

toolbox.register("mate", tools.cxPartialyMatched)
toolbox.register("mutate", tools.mutShuffleIndexes, indpb=0.05)
toolbox.register("select", tools.selTournament, tournsize=tourn_size)
toolbox.register("evaluate", evalTSP)

def main():
    random.seed(169)

    pop = toolbox.population(n=n_pop)

    hof = tools.HallOfFame(1)
    stats = tools.Statistics(lambda ind: ind.fitness.values)
    stats.register("avg", np.mean)
    stats.register("std", np.std)
    stats.register("min", np.min)
    stats.register("max", np.max)

    algorithms.eaSimple(pop, toolbox, 0.7, 0.2, n_iter, stats=stats, 
                        halloffame=hof)

    return pop, stats, hof

if __name__ == "__main__":
    pop, stats, hof = main()
    
    optim_idx = 0
    for i in range(1, n_pop):
        if evalTSP(pop[i])[0] > evalTSP(pop[optim_idx])[0]:
            optim_idx = i
            
    print("Optimal Path: {}, Total Distance: {}".format(np.array(pop[i]), evalTSP(pop[optim_idx])[0]))
```

・実行結果
```
gen	nevals	avg    	std    	min  	max  
0  	500   	52236.7	2308.02	45169	58507
1  	367   	49288  	2461.32	42273	57108
2  	386   	47777.5	3046.69	41320	56607
3  	391   	45884.8	3328.6 	41032	58617
4  	411   	42648.9	1603.21	39629	50939
5  	400   	41620.5	1660.14	38553	51171
6  	366   	40106.3	1391.94	37953	47475
7  	370   	39153.7	1413.43	37538	45846
8  	387   	38534  	1299.62	36377	45816
9  	395   	37826.6	1288.27	36042	45754
10 	410   	37156.1	1623.56	34851	45847
11 	387   	36370.5	1432.12	33891	44189
12 	363   	35747.3	1736.18	33579	43873
13 	384   	34848.5	1756.38	32581	42646
14 	386   	34132.9	1621.75	32188	41444
15 	378   	33351  	1732.46	32022	41444
16 	386   	32906.9	1609.93	31896	41367
17 	380   	32785.7	1582.83	31666	41058
18 	380   	32770.9	1959.11	31666	41664
19 	390   	32545.4	1901.79	31308	41322
20 	376   	32364.4	1763.74	31308	42296
Optimal Path: [ 96  87   1  91  93  73   8 102  35  83  14 117  67  95  23  92 108  81
  52  76  25  51  57  16  68  17  97  50  22  10   0 104  19  41  65  80
  74 107 118  72 100 101  49  78 116  21  15  32 109  62  38   4  36  66
  99  24  60  86  58  28  29  55  37  39  82 105   9  48  12  63  79  98
 103  31 110  94  11  61  53  27  75   5  56  26 114   6  71  54   3  59
 111   2  33  46  88  43  64  45 106  30  89  70  34  47  69 119  84  42
  18  85 112  77 115  44  90  13  40   7  20 113], Total Distance: 31308
```

実行結果の`min`の列に着目することにより徐々に最適解に近づいていることが確認できます。セールスマン巡回問題における交叉ではPMX(Partially Matched Crossover)という処理が用いられており、プログラムの`toolbox.register("mate", tools.cxPartialyMatched)`対応している点に注目しておくと良いです。