【遗传规划/计算智能】 彻底学会 DEAP 框架，从零搭建 GP

1. 应用方向

DEAP是一个新颖的进化计算框架，用于快速原型设计和测试。它旨在使算法清楚和数据结构透明。 它可以在并行机制之间完美协调，例如多处理和SCOOP（简单并发性面向对象程序设计）。

DEAP包括以下特征：

• （1）在遗传算法中可以使用任何你能想到的表示，例如：列表 List，数组Array，集合 Set，字典 Dictionary，树 Tree，Numpy Array 等等。
• （2）基因编程使用前缀树：Loosely typed, Strongly typed（理解为稀疏类型和紧凑类型）；自定义函数；
• （3）进化策略（包括CMA-ES）：Covariance Matrix Adaptation Evolutionary Strategies的缩写，中文名称是协方差矩阵自适应进化策略
• （4）多目标优化：包括 NSGA-II, SPEA2, MO-CMA-ES；
• （5）多种群（多人口）共同进化；
• （6）评估的并行化（和更多）；
• （7）名人堂，居住在人群中最好的人物；
• （8）检查点定期拍摄系统的快照；
• （9）基准模块包含最常用的测试功能；
• （10）进化的家谱（与NetworkX兼容）
• （11）其他算法的例子：PSO，差分演化，分布算法估计；

Core 核心模块：

• base： 基本结构。包含Toolbox(存储自定的EA运行所需的对象和操作)， Fitness(个体的适应度的基类)等。
• creator： 允许通过动态添加属性或函数来创建符合问题需求的类，常用来创建个体。
• tools：包含多种选择(selection)操作函数，遗传算子操作函数(多种crossover, mutation)等。

2. 安装

pip install deap

3. DEAP 优化问题的定义

单目标优化： creator.create('FitnessMin', base.Fitness, weights=(-1.0, ))

• 在创建单目标优化问题时，weights 用来指示最大化和最小化。此处 -1.0 即代表问题是一个最小化问题，对于最大化，应将 weights 改为正数，如 1.0。

• 另外即使是单目标优化，weights 也需要是一个 tuple，以保证单目标和多目标优化时数据结构的统一。

• 对于单目标优化问题，weights 的绝对值没有意义，只要符号选择正确即可。

• 多目标优化： creator.create('FitnessMulti', base.Fitness, weights=(-1.0, 1.0))

• 对于多目标优化问题，weights 用来指示多个优化目标之间的相对重要程度以及最大化最小化。如示例中给出的 (-1.0, 1.0) 代表对第一个目标函数取最小值，对第二个目标函数取最大值。

4. DEAP 个体编码

4.1 实数编码(Value encoding)：

直接用实数对变量进行编码。优点是不用解码，基因表达非常简洁，而且能对应连续区间。但是实数编码后搜索区间连续，因此容易陷入局部最优

from deap import base, creator, tools
import random

IND_SIZE = 5
creator.create('FitnessMin', base.Fitness, weights=(-1.0,)) # 优化目标：单变量，求最小值
creator.create('Individual', list, fitness = creator.FitnessMin) # 创建Individual 类，继承 list

toolbox = base.Toolbox()

# random.random 进行实数编码
toolbox.register('Attr_float', random.random)
toolbox.register('Individual', tools.initRepeat, creator.Individual, toolbox.Attr_float, n=IND_SIZE)

ind1 = toolbox.Individual()
print(ind1)

# 结果：[0.8579615693371493, 0.05774821674048369, 0.8812411734389638, 0.5854279538236896, 0.12908399219828248]

random.random 来产生随机数，如果喜欢 numpy 的，也可以用 np.random.rand 代替，DEAP 的灵活性还是很强的。

4.2 二进制编码(Binary encoding)：

以随机生成一个长度为 10 的二进制编码为例，本身 DEAP 库中没有内置的 Binary encoding，可以借助 Scipy 模块中的伯努利分布来生成一个二进制序列。

from deap import base, creator, tools
from scipy.stats import bernoulli

creator.create('FitnessMin', base.Fitness, weights=(-1.0,)) #优化目标：单变量，求最小值
creator.create('Individual', list, fitness = creator.FitnessMin) #创建Individual类，继承list

GENE_LENGTH = 10

toolbox = base.Toolbox()

# 实数编码
toolbox.register('Binary', bernoulli.rvs, 0.5) #注册一个Binary的alias，指向scipy.stats中的bernoulli.rvs，概率为0.5
toolbox.register('Individual', tools.initRepeat, creator.Individual, toolbox.Binary, n = GENE_LENGTH) #用tools.initRepeat生成长度为GENE_LENGTH的Individual

ind1 = toolbox.Individual()
print(ind1)

# 结果：[1, 0, 0, 0, 0, 1, 0, 1, 1, 0]

4.3 序列编码(Permutation encoding)：

通常在求解顺序问题时用到，例如TSP问题。序列编码中的每个染色体都是一个序列。

from deap import base, creator, tools
import random
creator.create("FitnessMin", base.Fitness, weights=(-1.0,))
creator.create("Individual", list, fitness=creator.FitnessMin)

IND_SIZE=10

toolbox = base.Toolbox()

# 序列编码
toolbox.register("Indices", random.sample, range(IND_SIZE), IND_SIZE)
toolbox.register("Individual", tools.initIterate, creator.Individual,toolbox.Indices)
ind1 = toolbox.Individual()
print(ind1)

#结果：[0, 1, 5, 8, 2, 3, 6, 7, 9, 4]

5. 初始种群的创建

一般种群的 DEAP 实现： toolbox.register('population', tools.initRepeat, list, toolbox.individual)

以二进制编码为例，以下代码可以生成由10个长度为5的随机二进制编码个体组成的一般种群：

from deap import base, creator, tools
from scipy.stats import bernoulli

# 定义问题
creator.create('FitnessMin', base.Fitness, weights=(-1.0,)) # 单目标，最小化
creator.create('Individual', list, fitness = creator.FitnessMin)

# 生成个体
GENE_LENGTH = 5
toolbox = base.Toolbox() #实例化一个Toolbox
toolbox.register('Binary', bernoulli.rvs, 0.5)
toolbox.register('Individual', tools.initRepeat, creator.Individual, toolbox.Binary, n=GENE_LENGTH)

# 生成初始族群
N_POP = 10
toolbox.register('Population', tools.initRepeat, list, toolbox.Individual)
toolbox.Population(n = N_POP)

# 结果：
# [[1, 0, 1, 1, 0],
# [0, 1, 1, 0, 0],
# [0, 1, 0, 0, 0],
# [1, 1, 0, 1, 0],
# [0, 1, 1, 1, 1],
# [0, 1, 1, 1, 1],
# [1, 0, 0, 0, 1],
# [1, 1, 0, 1, 0],
# [0, 1, 1, 0, 1],
# [1, 0, 0, 0, 0]]

6. DEAP 中评价定义（fitness）

评价部分是根据任务的特性高度定制的，DEAP 库中并没有预置的评价函数模版。

❗️ 在使用 DEAP 时，需要注意的是，无论是单目标还是多目标优化，评价函数的返回值必须是一个 tuple 类型。

from deap import base, creator, tools
import numpy as np
# 定义问题
creator.create('FitnessMin', base.Fitness, weights=(-1.0,)) #优化目标：单变量，求最小值
creator.create('Individual', list, fitness = creator.FitnessMin) #创建Individual类，继承 list

# 生成个体
IND_SIZE = 5
toolbox = base.Toolbox()

# 实数编码
toolbox.register('Attr_float', np.random.rand)
toolbox.register('Individual', tools.initRepeat, creator.Individual, toolbox.Attr_float, n=IND_SIZE)

# 生成初始族群
N_POP = 5
toolbox.register('Population', tools.initRepeat, list, toolbox.Individual)
pop = toolbox.Population(n = N_POP) # 产生 N_POP 个个体，每个个体都是由 list 的实数组成

# 定义评价函数
def evaluate(individual):
  return sum(individual), # 注意这个逗号，即使是单变量优化问题，也需要返回tuple

# 评价初始族群
toolbox.register('Evaluate', evaluate)
fitnesses = map(toolbox.Evaluate, pop) # 表示 pop 中每个个体都利用函数 toolbox.Evaluate 进行一次计算

for i in pop:
    print(i) # 打印每个个体的实数编码结果

for fit in fitnesses:
    print(fit) # 打印每个个体的 fitness 值

# 将个体 population 与对应的 fitness 绑定
for ind, fit in zip(pop, fitnesses):
    ind.fitness.values = fit
    print(ind.fitness.values) # 打印每个个体的 fitness 值

# 结果：
# (2.593989197511478,)
# (1.1287944225903104,)
# (2.6030877077096717,)
# (3.304964061515382,)
# (2.534627558467466,)

zip 函数详解： https://blog.csdn.net/PaulZhn/article/details/104391756

7. DEAP 内置的选择操作

函数	简介
selTournament()	锦标赛选择
selRoulette()	轮盘赌选择（不能用于最小化或者适应度会小于等于0的问题）
selNSGA2()	NSGA-II 选择，适用于多目标遗传算法
selSPEA2()	SPEA2 选择，该函数实现有误，没有为个体分配距离，不建议使用
selRandom()	有放回的随机选择
selBest()	选择最佳
selWorst()	选择最差
selTournamentDCD()	Dominance/Crowding distance锦标赛选择，目前版本的实现也有些问题
selDoubleTournament()	Size+Fitness双锦标赛选择
selStochasticUniversalSampling()	Size+Fitness双锦标赛选择
selStochasticUniversalSampling()	词典选择
selEpsilonLexicase()
selAutomaticEpsilonLexicase()

from deap import base, creator, tools
import numpy as np

# 定义问题
creator.create('FitnessMin', base.Fitness, weights=(-1.0,))  # 优化目标：单变量，求最小值
creator.create('Individual', list, fitness=creator.FitnessMin)  # 创建Individual类，继承list

# 生成个体
IND_SIZE = 5
toolbox = base.Toolbox()
toolbox.register('Attr_float', np.random.rand) # 实数编码
toolbox.register('Individual', tools.initRepeat, creator.Individual, toolbox.Attr_float, n=IND_SIZE)

# 初始化种群
N_POP = 5
toolbox.register('Population', tools.initRepeat, list, toolbox.Individual)
pop = toolbox.Population(n=N_POP)


# 定义评价函数
def evaluate(individual):
    return sum(individual),  # 注意这个逗号，即使是单变量优化问题，也需要返回tuple


# 评价初始种群
toolbox.register('Evaluate', evaluate)
fitnesses = map(toolbox.Evaluate, pop)
for ind, fit in zip(pop, fitnesses):
    ind.fitness.values = fit

# 选择方式1：锦标赛选择
toolbox.register('TourSel', tools.selTournament, tournsize=2)  # 注册 Tournsize 为 2 的锦标赛选择
selectedTour = toolbox.TourSel(pop, 5)  # 选择 5 个个体
print('锦标赛选择结果：')
for ind in selectedTour:
    print(ind)
    print(ind.fitness.values)

# 选择方式2: 轮盘赌选择
toolbox.register('RoulSel', tools.selRoulette)
selectedRoul = toolbox.RoulSel(pop, 5) # 选择 5 个个体
print('轮盘赌选择结果：')
for ind in selectedRoul:
    print(ind)
    print(ind.fitness.values)

# 选择方式3: 随机普遍抽样选择
toolbox.register('StoSel', tools.selStochasticUniversalSampling)
selectedSto = toolbox.StoSel(pop, 5) # 选择 5 个个体
print('随机普遍抽样选择结果：')
for ind in selectedSto:
    print(ind)
    print(ind.fitness.values)

#结果：
#锦标赛选择结果：
#[0.2673058115582905, 0.8131397980144155, 0.13627430737326807, 0.10792026110464248, 0.4166962522797264]
#(1.741336430330343,)
#[0.5448284697291571, 0.9702727117158071, 0.03349947770537576, 0.7018813286570782, 0.3244029157717422]
#(2.5748849035791603,)
#[0.8525836387058023, 0.28064482205939634, 0.9235436615033125, 0.6429467684175085, 0.5965523553349544]
#(3.296271246020974,)
#[0.5243293164960845, 0.37883291328325286, 0.28423194217619596, 0.5005947374376103, 0.3017896612109636]
#(1.9897785706041071,)
#[0.4038211036464676, 0.841374996509095, 0.3555644512425019, 0.5849111474726337, 0.058759891556433574]
#(2.2444315904271317,)
#轮盘赌选择结果：
#[0.42469039733882064, 0.8411201950346711, 0.6322812691061555, 0.7566549973076343, 0.9352307652371067]
#(3.5899776240243884,)
#[0.42469039733882064, 0.8411201950346711, 0.6322812691061555, 0.7566549973076343, 0.9352307652371067]
#(3.5899776240243884,)
#[0.5448284697291571, 0.9702727117158071, 0.03349947770537576, 0.7018813286570782, 0.3244029157717422]
#(2.5748849035791603,)
#[0.630305953330188, 0.09565983206218687, 0.890691659939096, 0.8706091807317707, 0.19708949882847437]
#(2.684356124891716,)
#[0.5961060867498598, 0.4300051776616509, 0.4512760237511251, 0.047731561819711055, 0.009892120639829804]
#(1.5350109706221766,)
#随机普遍抽样选择结果：
#[0.2673058115582905, 0.8131397980144155, 0.13627430737326807, 0.10792026110464248, 0.4166962522797264]
#(1.741336430330343,)
#[0.4038211036464676, 0.841374996509095, 0.3555644512425019, 0.5849111474726337, 0.058759891556433574]
#(2.2444315904271317,)
#[0.630305953330188, 0.09565983206218687, 0.890691659939096, 0.8706091807317707, 0.19708949882847437]
#(2.684356124891716,)
#[0.40659881466060876, 0.8387139101647804, 0.28504735705240236, 0.46171554118627334, 0.7843353275244066]
#(2.7764109505884718,)
#[0.42469039733882064, 0.8411201950346711, 0.6322812691061555, 0.7566549973076343, 0.9352307652371067]
#(3.5899776240243884,)

8. DEAP 内置的交叉(Crossover)操作

函数	简介	适用编码方式
cxOnePoint()	单点交叉	实数、二进制
cxTwoPoint()	两点交叉	实数、二进制
cxUniform()	均匀交叉	实数、二进制
cxPartialyMatched()	部分匹配交叉PMX	序列
cxUniformPartialyMatched()	PMX变种，改两点为均匀交叉	序列
cxOrdered()	有序交叉	序列
cxBlend()	混合交叉	实数
cxESBlend()	带进化策略的混合交叉
cxESTwoPoint()	带进化策略的两点交叉
cxSimulatedBinary()	模拟二值交叉	实数
cxSimulatedBinaryBounded()	有界模拟二值交叉	实数
cxMessyOnePoint()	混乱单点交叉	实数、二进制

from deap import base, creator, tools
import random

random.seed(42) # 保证结果可复现
IND_SIZE = 8
creator.create('FitnessMin', base.Fitness, weights=(-1.0, ))  # 优化目标：单变量，求最小值
creator.create('Individual', list, fitness = creator.FitnessMin)

# 创建两个序列编码个体
toolbox = base.Toolbox()
toolbox.register('Indices', random.sample, range(IND_SIZE), IND_SIZE)
toolbox.register('Individual', tools.initIterate, creator.Individual, toolbox.Indices)

ind1, ind2 = [toolbox.Individual() for _ in range(2)]
print(ind1, '\n', ind2)
# 结果：[1, 0, 5, 2, 7, 6, 4, 3]
# [1, 4, 3, 0, 6, 5, 2, 7]

# 单点交叉
child1, child2 = [toolbox.clone(ind) for ind in (ind1, ind2)]
tools.cxOnePoint(child1, child2)
print(child1, '\n', child2)
#结果：[1, 4, 3, 0, 6, 5, 2, 7]
# [1, 0, 5, 2, 7, 6, 4, 3]
# 可以看到从第四位开始被切开并交换了

# 两点交叉
child1, child2 = [toolbox.clone(ind) for ind in (ind1, ind2)]
tools.cxTwoPoint(child1, child2)
print(child1, '\n', child2)
# 结果：[1, 0, 5, 2, 6, 5, 2, 3]
# [1, 4, 3, 0, 7, 6, 4, 7]
# 基因段[6, 5, 2]与[7, 6, 4]互换了

# 均匀交叉
child1, child2 = [toolbox.clone(ind) for ind in (ind1, ind2)]
tools.cxUniform(child1, child2, 0.5)
print(child1, '\n', child2)
# 结果：[1, 0, 3, 2, 7, 5, 4, 3]
# [1, 4, 5, 0, 6, 6, 2, 7]

# 部分匹配交叉
child1, child2 = [toolbox.clone(ind) for ind in (ind1, ind2)]
tools.cxPartialyMatched(child1, child2)
print(child1, '\n', child2)
# 结果：[1, 0, 5, 2, 6, 7, 4, 3]
# [1, 4, 3, 0, 7, 5, 2, 6]
# 可以看到与之前交叉算子的明显不同，这里的每个序列都没有冲突

# 有序交叉
child1, child2 = [toolbox.clone(ind) for ind in (ind1, ind2)]
tools.cxOrdered(child1, child2)
print(child1, '\n', child2)
# 结果：[5, 4, 3, 2, 7, 6, 1, 0]
# [3, 0, 5, 6, 2, 7, 1, 4]

# 混乱单点交叉
child1, child2 = [toolbox.clone(ind) for ind in (ind1, ind2)]
tools.cxMessyOnePoint(child1, child2)
print(child1, '\n', child2)
# 结果：[1, 0, 5, 2, 7, 4, 3, 0, 6, 5, 2, 7]
# [1, 6, 4, 3]
# 注意个体序列长度的改变

9. DEAP 内置的突变(Mutation)操作

函数	简介	适用编码方式
mutGaussian()	高斯突变	实数
mutShuffleIndexes()	乱序突变	实数、二进制、序列
mutFlipBit()	位翻转突变	二进制
mutPolynomialBounded()	有界多项式突变	实数
mutUniformInt()	均匀整数突变	实数、序列
mutESLogNormal()

from deap import base, creator, tools
import random

random.seed(42) # 保证结果可复现
IND_SIZE = 5
creator.create('FitnessMin', base.Fitness, weights=(-1.0, ))
creator.create('Individual', list, fitness = creator.FitnessMin)

# 创建一个实数编码个体
toolbox = base.Toolbox()
toolbox.register('Attr_float', random.random)
toolbox.register('Individual', tools.initRepeat, creator.Individual, toolbox.Attr_float, IND_SIZE)

ind1 = toolbox.Individual()
print(ind1)
# 结果：[0.6394267984578837, 0.025010755222666936, 0.27502931836911926, 0.22321073814882275, 0.7364712141640124]

# 高斯突变
mutant = toolbox.clone(ind1)
tools.mutGaussian(mutant, 3, 0.1, 1)
print(mutant)
# 结果：[3.672658632864655, 2.99827700737295, 3.2982590920597916, 3.339566606808737, 3.6626390539295306]
# 可以看到当均值给到3之后，变异形成的个体均值从0.5也增大到了3附近

# 乱序突变
mutant = toolbox.clone(ind1)
tools.mutShuffleIndexes(mutant, 0.5)
print(mutant)
# 结果：[0.22321073814882275, 0.7364712141640124, 0.025010755222666936, 0.6394267984578837, 0.27502931836911926]

# 有界多项式突变
mutant = toolbox.clone(ind1)
tools.mutPolynomialBounded(mutant, 20, 0, 1, 0.5)
print(mutant)
# 结果：[0.674443861742489, 0.020055418656044655, 0.2573977358171454, 0.11555018832942898, 0.6725269223692601]

# 均匀整数突变
mutant = toolbox.clone(ind1)
tools.mutUniformInt(mutant, 1, 5, 0.5)
print(mutant)
# 结果：[0.6394267984578837, 3, 0.27502931836911926, 0.22321073814882275, 0.7364712141640124]
# 可以看到在第二个位置生成了整数3

10. DEAP 环境选择

环境选择也就是重插入(Reinsertion)，在选择、交叉和突变之后，得到的育种后代族群规模与父代相比可能增加或减少。为保持族群规模，需要将育种后代插入到父代中，替换父代种群的一部分个体，或者丢弃一部分育种个体。

• 重插入分为全局重插入(Global reinsertion)和本地重插入(Local reinsertion)两种，后者只有在使用含有本地邻域的算法时使用。常见的全局重插入操作有以下四种：
  • 完全重插入(Pure reinsertion)： 产生与父代个体数量相当的配种个体，直接用配种个体生成新一代族群。
  • 均匀重插入(Uniform reinsertion)： 产生比父代个体少的配种个体，用配种个体随机均匀地替换父代个体。
  • 精英重插入(Elitist reinsertion)： 产生比父代个体少的配种个体，选取配种后代中适应度最好的一些个体，插入父代中，取代适应度较低的父代个体。
  • 精英保留重插入(Fitness-based reinsertion)： 产生比父代个体多的配种个体，选取其中适应度最大的配种个体形成新一代族群。

通常来说后两种方式由于精英保留的缘故，收敛速度较快，因此比较推荐。

DEAP 中没有设定专门的 reinsertion 操作，可以用选择操作选择中的 selBest, selWorst, selRandom 来对育种族群和父代族群进行操作。

11. 数据统计

🌈 11.1 Computing Statistics

通常情况下，想要在优化过程中编辑数据。Statistic 模块 可以在任何设计好的目标上改变一些本不可改变的数据。为了达到这个目的，需要使用与工具箱中完全相同的语法在静态数据中注册统计函数。

states = tools.Statistics(key = lambda ind : ind.fitness.values)

• 使用 key 的第一个参数作为统计对象。这个 key 必须支持一个可以在之后被应用到数据上的函数从而得到统计结果。上面的例子使用了 ind.fitness.values中每一个元素的属性。

states.register('avg', numpy.mean)
states.register('std', numpy.std)
states.register('min', numpy.min)
states.register('max', numpy.max)

• 这些统计函数现在被注册了。register 函数希望 alias(‘名字’) 作为第一个属性以及一个在向量上操作的函数 (numpy.mean) 作为第二个属性。在调用时，任何在后面的元素都会传到函数上。统计目标的创建已经完成。

🎉 11.2 Predefined Algorithms

当使用一个预定义的算法时，例如 esSimple()、eaMuPlusLambada()、eaMuCommaLambda()、eaGenerateUpdata()，之前创建的统计目标可以作为算法的属性。

pop, logbook = algorithms.esSimple(pop, toolbox, cxpb = 0.5, mutpb = 0.2, ngen = 0, stats = stats, verbose = True)

统计将会在每一次迭代中自动的进行计算。详细参数在优化过程中会打印在屏幕上。一旦算法返回，最终的种群和一个 logbook 将会返回。在下面可以看到更详细的信息。

💛11.3 Writing Your Own Algorithm

• 当编写自己的算法时，包含统计是十分简单的。只用在需要的目标上编写统计。例如，在一个给定种群 pop 上编写统计需要调用 compile() 方法完成。

record = stats.compile(pop)

• 用于编辑函数的属性必须在一个迭代元素中，这样这些 key 才会被调用。

• 这里，种群(pop)包含了许多个体。统计目标将会在每一个个体上调用 key 函数获取 fitness.values 属性的值。这个结果数组的值最终会给到每一个统计函数并且将结果输入到 record 字典中，每一个 key 都会与相应的函数相关联。

• 在 main 函数中把这个命令放在更新种群之后就可以得到一系列的统计值：标准差-最大值-平均值-最小值。

💖11.4 Multi-objective Statistics

正如统计可以通过 numpy 函数直接进行计算，所有的目标将会通过默认 numpy 的属性联合在一起。接下来，一个需要明确的事情是每一个 axis 的操作。这会通过给予 axis 一个额外的属性作为注册函数达成。

stats = tools.Statistics(key=lambda ind: ind.fitness.values)
stats.register("avg", numpy.mean, axis=0)
stats.register("std", numpy.std, axis=0)
stats.register("min", numpy.min, axis=0)
stats.register("max", numpy.max, axis=0)

✨11.5 Multiple Statistics

计算种群个体的不同属性也是可以的。例如，在遗传程序设计(GP)中除了对适应度统计之外，对表达式树的高度统计也是常见的。一个可以使用多元 Statistics 目标的函数是 MultiStatistics。

stats_fit = tools.Statistics(ket= lambda ind: ind.fitness.values)
stats_size = tools.Statistics(key = len)
mstats = tools.MultiStatistics(fitness = stats_fit, size = stats_size)

两个统计目标使用与之前相同的方式被创建。第二个目标将会通过调用 len() 获取每一个个体的长度。一旦创建完成，统计目标将会被给到 MultiStatistics 函数中，这里地元素都是使用 keywords 来定义的。这些 keywords 将会提供不同统计量的定义。这些统计函数仅仅可以在 multi-statistics 被注册一次

mstats.register("avg", numpy.mean)
mstats.register("std", numpy.std)
mstats.register("min", numpy.min)
mstats.register("max", numpy.max)

多元统计目标可以给到一个算法当中，或者他们可以以相似的形式运用到单一统计中

record = mstats.compile(pop)

在这种情况下，record 就是一个字典的字典类型 (嵌套)。第一个等级包括一些统计信息的关键字，第二个等级包括一些预制件相似的简单统计目标

>>> print(record)
{'fitness': {'std': 1.64, 'max': 6.86, 'avg': 1.71, 'min': 0.166},
'size': {'std': 1.89, 'max': 7, 'avg': 4.54, 'min': 3}}

⭐️11.6 Logging Data

一旦数据通过统计产生，可以使用 Logbook 对它进行存储。Logbook 是用来
按时间顺序排列的条目（如字典）。它会直接兼容数据类型并且返回统计目标，但是不会被数据限制。实际上，任何东西都可以包含在日志的条目中。

logbook = tools.Logbook()
logbook.record(gen=0, evals=30, **record)

• record() 方法采用可变数字作为参数，每一个参数作为数据都会被记录。
• 在最后一个例子中，保存了一代，结果的数量和任何包含在 record() 中的东西都会通过这个方法产生统计学目标。所有的记录都会被保存在 Logbook 中，直到它被销毁。在一些列的记录之后，可能会想去调用 logbook 中的信息

gen, avg = logbook.select("gen", "avg")

select() 方法提供了一种方法去调用所有的与 record 中 keyword 相关的信息。这个方法采取了可变数量的字符串元素，就是在 record 或者 statistics 中的 keywords。这里，我们调用代数和平均适应度，使用一个单独的 select 进行调用。一个 logbook 是一个可选择的目标(只要插入数据是可选择的)，它可以提供一种很好的方法去存储演化过程中的统计参数。

import pickle
pickle.dump(logbook, lb_file)

• 注意： 每一个算法返回的 logbook 包含着每一代的信息和适应度在整个进化过程中的数目。

⛄️11.7 Printing to Screen

一个 logbook 可以打印在屏幕或者文档中。它的 str() 方法在第一个 key 中返回每一个 key 的头，同时使用这些 keys 完成 logbook。按照行插入的时间顺序排列，而列将处于未定义的顺序。指定顺序的最简单方法是将 header 属性设置为指定列顺序的字符串列表。

logbook.header = 'gen', 'avg', 'spam'

结果为：

>>>print(logbook)
gen avg spam
0     [ 50, 2 ]

0 [ 50, 2 ] 每一个列名包含了没有明确记录的列名，它将会被空着，就像 spam 一样。

一个 logbook 同样包含尚未打印的条目的流属性。

>>> print(logbook.stream)
gen   avg      spam
0     [ 50.2]
>>> logbook.record(gen=1, evals=15, **record)
>>> print(logbook.stream)
1     [ 50.2]

🐾11.8 Dealing with Multi-statistics

Logbook 可以与字典类型合作，返回 MultiStatistics 目标。事实上，它将在 chapters 中为每一个包含在字典中的子字典记录数据，接下来，一个 multi record 可以被看做一个 record

logbook = tools.Logbook()
logbook.record(gen=0, evals=30, **record)

与列排序不同的一点是，当我们明确 chapters 的顺序时，它们的内容如下

logbook.header = "gen", "evals", "fitness", "size"
logbook.chapters["fitness"].header = "min", "avg", "max"
logbook.chapters["size"].header = "min", "avg", "max"

(看起来就是 fitness 包含 min、avg、max；size 也包含这些，然后它们的输出就都分别包含这些)

>>> print(logbook)
                     fitness                 size
            -------------------------  ---------------
gen   evals min      avg      max      min   avg   max
0     30    0.165572 1.71136  6.85956  3     4.54  7

这些数据的调用同样也可以通过 chapter 进行调用

gen = logbook.select("gen")
fit_mins = logbook.chapters["fitness"].select("min")
size_avgs = logbook.chapters["size"].select("avg")

• 迭代次数、最小适应度和平均尺寸就可以按照时间顺序获得了。如果一些数据不可用，一个 None 会出现在向量中。

💗11.9 Some Plotting Sugar

在优化过程中最常用的操作就是在图中显示进化过程。Logbook 可以有效地执行这一操作。使用选择方法，我们可以调用需要的数据并且使用 matplotlib 去绘制图形。

gen = logbook.select("gen")
fit_mins = logbook.chapters["fitness"].select("min")
size_avgs = logbook.chapters["size"].select("avg")

import matplotlib.pyplot as plt

fig, ax1 = plt.subplots()
line1 = ax1.plot(gen, fit_mins, "b-", label="Minimum Fitness")
ax1.set_xlabel("Generation")
ax1.set_ylabel("Fitness", color="b")
for tl in ax1.get_yticklabels():
    tl.set_color("b")

ax2 = ax1.twinx()
line2 = ax2.plot(gen, size_avgs, "r-", label="Average Size")
ax2.set_ylabel("Size", color="r")
for tl in ax2.get_yticklabels():
    tl.set_color("r")

lns = line1 + line2
labs = [l.get_label() for l in lns]
ax1.legend(lns, labs, loc="center right")

plt.show()

11. 基于 DEAP 的遗传编程实现

DEAP 官方文档给出的一个简单GP例子 (符号回归问题)：Symbolic Regression Problem: Introduction to GP

令期望曲线为：
$$x^4+x^3+x^2+x$$

并在 [-1, 1] 区间设置 20 个等分点来评估适应度(fitness)，来寻找一个 GP 生成的最拟合的表达式。

以下是对它的讲解，代码部分的 comment 当中也加入了对代码过程和操作的解读。

1. 导入模块

import math
import random
import operator  # python内部操作符，包含对象比较、逻辑比较、算术运算和序列操作
import numpy
# deap的常用五个模块齐了
from deap import base    # core
from deap import creator # core
from deap import tools	 # core
from deap import gp 
from deap import algorithms

2. Primitive Set 的创建

• 注： Primitive set 为 function set 和 terminal set 的并集。primitives 不是指 functions + terminals，primitives 仅仅指 functions

def protectedDiv(left, right):  # 保护性除法，防崩
    try:
        return left / right
    except ZeroDivisionError:
        return 1

pset = gp.PrimitiveSet("MAIN", 1)  # "MAIN"是名称， 1为程序的inputs数量

# primitive set有了，可以向里面添加function(primitive)和terminal了
pset.addPrimitive(operator.add, 2) 
pset.addPrimitive(operator.sub, 2)
pset.addPrimitive(operator.mul, 2)
pset.addPrimitive(protectedDiv, 2)
pset.addPrimitive(operator.neg, 1)
pset.addPrimitive(math.cos, 1)
pset.addPrimitive(math.sin, 1)

# Ephemeral Constant: 没有固定值的"常量"
# 当这个节点被选择加到树里面时，所包含的函数被执行,函数输出被作为一个constant加入树中。
# 本质上它还是一个终止符，但这里不是用 addTerminal 方法生成的
pset.addEphemeralConstant("rand101", lambda: random.randint(-1, 1))  # (随机-1, 0, 1)

pset.renameArguments(ARG0='x')  # 程序inputs名称默认为ARG0,ARG1,ARG2,...,可自行修改
								# 上面设置了程序inputs数量为1， 所以这里可以改一个ARG0

这部分主要使用 gp 模块创建 primitive set 并向里面添加函数(function/primitive) 和终止符 (terminal)。
gp 模块两种 primitive set 的创建：

1. Loosely Typed GP： 不指定节点的输入和输出类型(type)。常用于较为单一处理类型的问题，比如我们现在看的这个只用到浮点数 float 的符号回归。

   创建 primitive set 使用 gp.PrimitiveSet(名称, inputs数量, prefix=‘ARG’)
   添加函数节点：gp.addPrimitive(函数, 参数数量, name=None)
   添加终止符: gp.addTerminal(终止符, name=None)

2. Strongly Typed GP：指定每一个节点输入和输出的类型 (它真的很严格 ! )。在Srong Typed GP的随机过程中，如果一个下层节点的输出类型和它连接的一个上层节点输入类型不一致，那么存在这种连接的树会被丢弃。

   创建 primitive set 使用 gp.PrimitiveSetTyped(名称, 输入types, 返回/输出type,prefix=‘ARG’)
   添加函数节点: addPrimitive(函数, 输入types, 返回/输出type, name=None)
   添加终止符: addTerminal(终止符, 返回/输出type, name=None)

gp 模块方法官方文档：https://deap.readthedocs.io/en/master/api/gp.html

2.1 Primitive Set 参数说明

pset = gp.PrimitiveSet("MAIN", 13)

• pset.arguments：自变量列表
• pset.primitives：存储树中所有的函数节点，node
• pset.terminals：存储树中所有的终端节点
• node.ret：节点类型
• node.args：节点参数

3. creator 创建个体类

# 使用 creator.create 创建基因型(genotype)的个体和适应度(fitness)
# weights*必须*用以元组表达，这里我们要一个minimizing fitness，故传入-1.0
creator.create("FitnessMin", base.Fitness, weights=(-1.0,))
creator.create("Individual", gp.PrimitiveTree, fitness=creator.FitnessMin)

creator 允许我们动态添加函数、属性来创建自己的自定类型。

• evolutionary program 中，至少需要用 creator.create 创建两个类型：个体类和它对应的适应度的类。个体类表现为一个 gp.PrimitiveTree，适应度表现为 base 模块中的 Fitness 基类。
• 我们一会用到的 fitness function 使用均方误差来评估个体适应度，越小的值表明个体越优良，所以我们想要一个 minimizing fitness。

creator 模块方法官方文档： https://deap.readthedocs.io/en/master/api/creator.html

base 模块官方文档： https://deap.readthedocs.io/en/master/api/base.html

4. Toolbox

“注册一些工具”：使用 register 将自定函数填充到工具箱 (base.Toolbox()) 当中，在之后的算法部分可以通过 toolbox.name 调用 ☆☆☆。算法部分，暨遗传迭代的实行部分，可以自定、也可使用 algorithms 模块中提供的算法方法。使用 algorithms 模块中的方法时，工具箱注册的函数名称固定，如适应度评估必须注册为 “evaluate”, 遗传操作交叉注册为 “mate”，变异注册为 “mute” 等。

# 过程中所需参数动态绑定
toolbox = base.Toolbox()
toolbox.register("expr", gp.genHalfAndHalf, pset=pset, min_=1, max_=2) # 利用 half and half 的形式产生树结构，树形结果的最小深度为 2， 最大深度为 6
toolbox.register("individual", tools.initIterate, creator.Individual, toolbox.expr)

# 注册 toolbox。population， 并注册 toolbox.compile 对个体树形表达式转化为可编译函数
toolbox.register("population", tools.initRepeat, list, toolbox.individual)
toolbox.register("compile", gp.compile, pset=pset)

# 评估个体的函数，即 fitness
def evalSymbReg(individual, points):
    # compile: 将表达式转换为一个可调用的函数
    # 这里将individual的树形表现形式转换为executable形式, 得以输入参数得到输出来和期望输出进行比较
    func = toolbox.compile(expr=individual)
    
    # 均方误差评估个体适应度
    sqerrors = ((func(x) - x**4 - x**3 - x**2 - x)**2 for x in points)
    return math.fsum(sqerrors) / len(points),  # 注意这里 return一个tuple 
    										   # (*DEAP存储 fitness为可迭代类型)

toolbox.register("evaluate", evalSymbReg, points=[x/10. for x in range(-10,10)])
toolbox.register("select", tools.selTournament, tournsize=3) # 选择
toolbox.register("mate", gp.cxOnePoint) # 交叉
toolbox.register("expr_mut", gp.genFull, min_=0, max_=2)
toolbox.register("mutate", gp.mutUniform, expr=toolbox.expr_mut, pset=pset) # 变异

# 装饰器: 一个包装函数，在被包装函数被调用前和调用后执行一些操作 (initializaiton and termination)
# 这里是分别对之前注册的mate函数和mutate函数进行包装，限制树形最大深度为17
toolbox.decorate("mate", gp.staticLimit(key=operator.attrgetter("height"), max_value=17))
toolbox.decorate("mutate", gp.staticLimit(key=operator.attrgetter("height"), max_value=17))

对于 register 的过程理解，我们来体会一下在代码段中的两个注册例子：

Example 1 :

  toolbox.register("mate", gp.cxOnePoint)

过程解读： 将 gp.cxOnePoint 函数命名为 “mate”，放入工具箱中。可以看成一个 重命名的过程。

执行时，调用工具箱的 “mate” 即通过调用 toolbox.mate，便会直接调用 gp.cxOnePoint ，进行交叉。

注意： 现在这句里只有两个参数，名称和被注册函数，但并不代表不需要给 gp.cxOnePoint 这个被注册函数某些参数。register 中包含一个偏函数partial，可以让我们在之后再以字典方式传入参数。下面的例句同样验证了这一点。

Example 2 :

toolbox.register("population", tools.initRepeat, list, toolbox.individual)

过程解读： 这是种群列表的注册，注册名称 “population”，tools.initRepeat 为被注册函数，后面紧跟的 list, toolbox.individual 为 tools.initRepeat 的参数。initRepeat 还应该接收一个数量 n 为参数才能正常执行，可以看到这行语句没有任何参数 n，但在 6. Launching 部分(本文后面) 的语句 pop = toolbox.population(n=300) 中，参数 n 以字典形式传入了注册函数 population中。 执行 toolbox.population(n=300) 时，tools.initRepeat(list, toolbox.individual, n=300) 被调用，它会进行 300 次调用 toolbox.individual 并将返回的结果填入 list 结构的操作 (toolbox.population() 返回一个list)。同理，toolbox.individual 调用 tools.initIterate(creator.Individual, toolbox.expr) ，将 toolbox.expr 返回值填充入 creator.Individual 类型。

后面使用的 evolutionary algorithms 会抽调这些注册的操作。

tools模块方法官方文档： https://deap.readthedocs.io/en/master/api/tools.html

5. 统计数据

    stats_fit = tools.Statistics(lambda ind: ind.fitness.values)
    stats_size = tools.Statistics(len)
    # 创建一个有着两个统计目标的统计对象
    mstats = tools.MultiStatistics(fitness=stats_fit, size=stats_size)
    # 注册统计目标的统计项
    mstats.register("avg", numpy.mean)
    mstats.register("std", numpy.std)
    mstats.register("min", numpy.min)
    mstats.register("max", numpy.max)

设置统计项，注册到统计容器 mstats 中，在后面能够通过 log 观察运行过程中每一代的数据。

6. Launching

    pop = toolbox.population(n=300)  # 设置种群数量
    hof = tools.HallOfFame(1)		 # 名人堂最多存储一个个体，即最后产生的最好个体数量
    pop, log = algorithms.eaSimple(pop, toolbox, 0.5, 0.1, 40, stats=mstats, halloffame=hof, verbose=True)

初始化、运行。

HallOfFame 存储最适合的个体，传入参数为 1 则其包含一个最适合的个体。

这里使用 algoritms 模块的 eaSimple，一个最基础的算法。传入参数分别为种群对象，toolbox 工具箱，crossover 概率，mutate 概率，演化多少代，统计容器，名人堂，verbose=true 设置显示 statistics。

eaSimple 算法会自动抽调传入的 toolbox 里注册过的工具，toolbox 中必须按特定名称注册对应操作的工具，如 “mate”, “mutate”, “select”, “evaluate” 这些一定要有，不能随意起名字…

可以在运行结束后通过查看名人堂 hof 来找到 GP 生成的表达式。如下直接打印出来：

print(str(hof.items[0]))
print(hof.keys[0])
print(hof.items[0].fitness)

# 输出:
# add(mul(x, sub(x, neg(mul(x, sub(x, neg(mul(x, x))))))), x)
# (5.1229736520700476e-33,)
# (5.1229736520700476e-33,)

之前设置名人堂数量为 1，所以直接找 hof.items[0] 就可以。将打印出的语法树转为表达式为 $x+x^2+x^3+x^4$ ，GP 找到的最优解即为我们的期望曲线。
algorithms 模块方法官方文档：https://deap.readthedocs.io/en/master/api/algo.html

7. 用法技巧总结

spt = random.randint(0, 3)
node = par[spt]  # par 为通过 Primitive Set 产生的树结构，在 par 树中随机选择一个节点

[p for p in pset.primitives[node.ret] if p.args == node.args] # 提取同 node.ret 同类型,且参数相同个数的节点

pset.add_index=[primitive_index for primitive_index in range(len(pset.primitives[pset.ret])) if pset.primitives[pset.ret][primitive_index].name=="add"][0]

12. 利用DEAP库创建一个GP表达式

（1）首先创建 primitive set：函数集 + 终点集

import operator
from deap import gp
import math
import random

# 创建一个 primitive set: primitiveset（名称，变量个数）  
pset = gp.PrimitiveSet('main111',2)
# 变量默认叫ARG0,ARG1...可以自己改名字,这里的输入参数可以作为terminal
pset.renameArguments(ARG0 = 'x')
pset.renameArguments(ARG1 = 'y')

# 定义function set：addPrimitive(函数名,参数个数)  不同的函数需要的参数个数不一样
# 函数可以是已有的
pset.addPrimitive(operator.add,2)
pset.addPrimitive(operator.mul,2)
pset.addPrimitive(operator.neg,1)  # 取负值
pset.addPrimitive(math.sin,1)
# 函数也可以自己定义
def if_then_else(constraint, output1,output2):
    if constraint:
        return output1
    else:
        return output2
pset.addPrimitive(if_then_else,3)

# 定义terminal set
pset.addTerminal(3)  # 常数
pset.addEphemeralConstant('num',lambda:random.randint(-5,5))  # 随机数：(name,随机数函数)

（2）根据上述 primitive set 创建表达式

# 三种生成方法,full,grow,half(primitive set,min_depth,max_depth)
# python中树的最大深度是91,最小深度是0（只有一个叶子节点）
expr1 = gp.genFull(pset,0,2)    # list
expr2 = gp.genGrow(pset,1,3)
expr3 = gp.genHalfAndHalf(pset,1,3)
tree = gp.PrimitiveTree(expr1)  # 将表达式形式转换成树形结构


# 编译表达式：必须经过编译才能进行运算
function1 = gp.compile(tree,pset)  
# function2 = gp.compile(expr2,pset)  
result = function1(1,2)   # 输入变量值能够得到该表达式的结果
print('result:',result)

• 注意：步骤（2）只是创建一个gp表达式、编译表达式的示例，实际使用deap库时，直接进行步骤（3）即可。

（3）创建个体（individual）

from deap import creator,tools,base

# create(类名，继承的类，参数)
creator.create('FitnessMin',base.Fitness,weights=(-1.0,))
creator.create('Individual',gp.PrimitiveTree,fitness = creator.FitnessMin,pset = pset) 
# Individual继承的类是"gp.PrimitiveTree",生成tree-based individual

# 向toolbox注册函数
toolbox = base.Toolbox()   # 实例
# register(函数别名，函数，函数需要的参数)
toolbox.register('expr',gp.genHalfAndHalf,pset = pset,min_ = 3,max_ = 3)  
# 也可以根据需要用genFull、genGrow
toolbox.register('individual',tools.initIterate,creator.Individual,toolbox.expr)

# 生成一个个体
one_individual = toolbox.individual()
# print('individual:',one_individual)

13. DEAP 内置算法（Algorithms）

算法	简介
eaSimple()
eaMuPlusLambda()
eaMuCommaLambda()
eaGenerateUpdate()
StrategyOnePlusLambda()
StrategyMultiObjective()

14. DEAP 实现 NSGAII 算法

其中个体编码方式前三个变量用标签编码，后面两个使用实数编码，因此，在交叉变异方法，前面使用的是单点交叉和变异，后面使用的是模拟二进制交叉和变异。

import pandas as pd
import os

# 查看当前工作路径
print("当前工作路径:", os.getcwd())

os.chdir("D:/08Code/C题")


data=pd.read_csv('data.csv')
# 查看更改后的工作路径
print("更改后的工作路径:", os.getcwd())

data['励磁波形']=data['励磁波形'].map({'正弦波': 0, '三角波': 1, '梯形波': 2})
data['磁芯材料']=data['磁芯材料'].map({'材料1': 0, '材料2': 1, '材料3': 2, '材料4':3})
data['温度，oC']=data['温度，oC'].map({25: 0, 50: 1, 70: 2, 90:3})

subset = data.iloc[:, 5:]
# 计算每行的最大值并形成新的一列
data['磁通密度峰值'] = subset.max(axis=1)

X=pd.concat([data[['温度，oC','磁芯材料' ,'励磁波形']].reset_index(drop=True),data[['频率，Hz','磁通密度峰值']].reset_index(drop=True)],axis=1)
Y = data['磁芯损耗，w/m3'].copy()

from sklearn.ensemble import RandomForestClassifier, AdaBoostClassifier, GradientBoostingClassifier, \
    RandomForestRegressor, AdaBoostRegressor, GradientBoostingRegressor

import xgboost as xgb
model = xgb.XGBRegressor(random_state=42)
model.fit(X, Y)

# 定义变量上下界
temperatures = [25, 50, 70, 90]
materials = ['材料1', '材料2', '材料3', '材料4']
waveforms = ['正弦波', '三角波', '梯形波']
freq_min, freq_max = 49990, 501180
flux_density_min, flux_density_max = 0.00963815, 0.313284469

def flatten_individual(ind):
    # 检查最后一个元素是否是列表
    if isinstance(ind[-1], list):
        # 将列表中的元素展平并替换
        ind[:] = ind[:-1] + ind[-1]
    return ind

import numpy as np
from deap import base, tools, creator, algorithms
import random
import matplotlib.pyplot as plt

#定义问题
creator.create('FitnessMinMax',base.Fitness,weights=(-1.0, 1.0)) # 两个目标，最小化磁芯损耗 P, 最大化传输磁能f·B
creator.create('Individual',list,fitness = creator.FitnessMinMax) #创建individual类

toolbox  = base.Toolbox()

# 属性生成器
toolbox.register("attr_temperature", random.randint, 0, len(temperatures)-1)
toolbox.register("attr_material", random.randint, 0, len(materials)-1)
toolbox.register("attr_waveform", random.randint, 0, len(waveforms)-1)
toolbox.register("attr_frequency", random.uniform, freq_min, freq_max)
toolbox.register("attr_flux_density", random.uniform, flux_density_min, flux_density_max)

# 结构初始化器
toolbox.register("individual", tools.initCycle, creator.Individual,
                 (toolbox.attr_temperature, toolbox.attr_material, toolbox.attr_waveform,
                  toolbox.attr_frequency, toolbox.attr_flux_density), n=1) # 注册个体生成工具
toolbox.register("population", tools.initRepeat, list, toolbox.individual)

def fitnessfunction(ind):
    # ind = [int(round(val)) if idx < 3 else float(val) for idx, val in enumerate(ind)]

    ind[:] = flatten_individual(ind[:])

    t = ind[:]
    t = np.reshape(t, (1, -1))
    # print(t)
    t = t.tolist()
    f1 = model.predict(t)

    # 如果 f1 <= 0，返回一个非常大的值
    if f1[0] <= 0:
        f1[0] = float('inf')

    f2 = ind[3] * ind[4]

    return f1[0], f2

toolbox.register('evaluate', fitnessfunction)#注册评价函数

N_POP = 300 # 种群内个体数量，参数过小，搜索速度过慢
toolbox.register('population',tools.initRepeat,list,toolbox.individual) # 注 册种群生成工具
pop = toolbox.population(n=N_POP) # 建立种群pop
# for ind in pop:#打印一个种群检查
#   print(ind)

toolbox.register('selectGen1', tools.selTournament, tournsize=2)
toolbox.register('select', tools.emo.selTournamentDCD) # 该函数是 binary tournament，不需要 tournsize

# 自定义交叉函数：分类标签用单点交叉，实数用模拟二进制交叉
def custom_mate(ind1, ind2):
    # 进行分类标签的单点交叉
    ind1[:3], ind2[:3] = tools.cxOnePoint(ind1[:3], ind2[:3])

    # 对实数部分进行模拟二进制交叉
    ind1[3:], ind2[3:] = tools.cxSimulatedBinaryBounded(ind1[3:], ind2[3:], eta=20.0, low=[freq_min, flux_density_min],
                                   up=[freq_max, flux_density_max])

    # 强制前三个为整数
    ind1[0], ind1[1], ind1[2] = int(round(ind1[0])), int(round(ind1[1])), int(round(ind1[2]))
    ind2[0], ind2[1], ind2[2] = int(round(ind2[0])), int(round(ind2[1])), int(round(ind2[2]))

    return ind1, ind2

toolbox.register('mate', custom_mate)

# 自定义变异函数：分类标签随机变异，实数用多项式变异
def custom_mutate(ind):
    # 对分类标签部分进行随机变异
    if random.random() < 0.5:
        ind[0] = random.randint(0, len(temperatures) - 1)
    if random.random() < 0.5:
        ind[1] = random.randint(0, len(materials) - 1)
    if random.random() < 0.5:
        ind[2] = random.randint(0, len(waveforms) - 1)

    # 对实数部分进行多项式变异
    ind[3:] = tools.mutPolynomialBounded(ind[3:], eta=20.0, low=[freq_min, flux_density_min],
                               up=[freq_max, flux_density_max], indpb=1.0 / len(ind[3:]))

    # 强制前三个为整数
    ind[0], ind[1], ind[2] = int(round(ind[0])), int(round(ind[1])), int(round(ind[2]))

    return ind,

toolbox.register('mutate', custom_mutate)

NGEN = 500 # 迭代步数，参数过小，在收敛之前就结束搜索
CXPB = 0.8 # 交叉概率，参数过小，族群不能有效更新
MUTPB = 0.2 # 突变概率，参数过小，容易陷入局部最优

#开始迭代
#第一代与第二代之后的代数操作不同
fitnesses = toolbox.map(toolbox.evaluate, pop)#为初始种群进行评价，得到适应度
for ind,fitness in zip(pop,fitnesses):
    ind.fitness.values = fitness # 这里先有一个适应度才能进行快速非支配排序


fronts = tools.emo.sortNondominated(pop, k=N_POP, first_front_only=False) # 快速非支配排序，得到不同前沿的 pareto 层集合fronts
# print(fronts )#打印一个pareto层集合fronts检查，每层前沿都是一个列表

for idx,front in enumerate(fronts): # 使用枚举循环得到各层的标号与 pareto 解
    # print(idx,front) # 打印检查前沿标号与每层的pareto解
    for ind in front:
        ind.fitness.values = (idx+1), # 将个体的适应度设定为 pareto 解的前沿次序


offspring = toolbox.selectGen1(pop, N_POP)  #进行选择

# for i in range(len(offspring) - 1):
#     print(offspring[i])

offspring = algorithms.varAnd(offspring, toolbox, CXPB, MUTPB) # 只做一次交叉与变异操作



for i in range(len(offspring)):
    offspring[i][:] = flatten_individual(offspring[i][:])
    # print(offspring[i])

#从第二代开始循环
for gen in range(1, NGEN):
    # print("第：" + str(gen)+ + "开始")
    combinedPop = pop + offspring#将子代与父代结合成一个大种群

    fitnesses = toolbox.map(toolbox.evaluate, combinedPop) # 对该大种群进行适应度计算
    for ind,fitness in zip(combinedPop,fitnesses):
        ind.fitness.values = fitness

    fronts = tools.emo.sortNondominated(combinedPop,k=N_POP,first_front_only=False)#对该大种群进行快速非支配排序

    for front in fronts:
        tools.emo.assignCrowdingDist(front)#计算拥挤距离
    pop = []
    for front in fronts:
        pop += front
    pop = toolbox.clone(pop)
    pop = tools.selNSGA2(pop,k=N_POP,nd='standard')#基于拥挤度实现精英保存策略

    offspring = toolbox.select(pop,N_POP)#选择
    offspring = toolbox.clone(offspring)
    offspring = algorithms.varAnd(offspring,toolbox,CXPB,MUTPB)#交叉变异

bestInd = tools.selBest(pop,1)[0]#选择出种群中最优个体
bestFit = bestInd.fitness.values
print('best solution:',bestInd)
print('best fitness:',bestFit)

front = tools.emo.sortNondominated(pop,len(pop))[0] # 返回的不同前沿的pareto层集合fronts中第一个front为当前最优解集
# # 图形化显示
# for ind in front:
#     plt.plot(ind.fitness.values[0],ind.fitness.values[1],'r.',ms=2)
# plt.xlabel('f1')
# plt.ylabel('f2')
# plt.tight_layout()
# plt.show()

Reference

[1] 【遗传编程/基因规划】python DEAP框架学习笔记

github 源码： https://gitcode.net/mirrors/DEAP/deap?utm_source=csdn_github_accelerator

DEAP 文档： http://deap.gel.ulaval.ca/doc/dev/index.html

DEAP 中文介绍： https://www.jianshu.com/p/8fa044ed9267

Poli, Riccardo, et al. A Field Guide to Genetic Programming. Lulu Press, 2008.