(六) GP實(shí)例——函數(shù)發(fā)現(xiàn)
話不多說直接上代碼
import operator, math, random, numpy
from deap import gp, tools, base, creator, algorithms
toolbox = base.Toolbox()
pset = gp.PrimitiveSet("MAIN", 1)
#define protectedDive function which return their division
def protectedDiv(left, right):
try:
return left / right
except ZeroDivisionError:
return 1
#adding other functions
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)
#add random constants which is an int type
pset.addEphemeralConstant("rand1", lambda: random.randint(-1,1))
#rename augument x
pset.renameArguments(ARG0='x')
#creating MinFit
creator.create("FitnessMin", base.Fitness, weights=(-1.0,))
#creating individual
creator.create("Individual", gp.PrimitiveTree, fitness=creator.FitnessMin)
#import toolbox
toolbox = base.Toolbox()
#resigter expr individual population and compile
toolbox.register("expr", gp.genHalfAndHalf, pset=pset, min_=1, max_=2) #0~3
toolbox.register("individual", tools.initIterate, creator.Individual, toolbox.expr)
toolbox.register("population", tools.initRepeat, list, toolbox.individual)
toolbox.register("compile", gp.compile, pset=pset)
#define evaluating function
def evalSymbReg(individual, points):
# Transform the tree expression in a callable function
func = toolbox.compile(expr=individual)
# Evaluate the mean squared error between the expression
# and the real function : x**4 + x**3 + x**2 + x
sqerrors = ((func(x) - x**4 - x**3 - x**2 - x)**2 for x in points) #define their differential
#return the average fitness
return math.fsum(sqerrors) / len(points),
#register genetic operations(evaluate/selection/mutate/crossover/)
toolbox.register("evaluate", evalSymbReg, points=[x/10. for x in range(-10,10)]) #the range of points in evaluate function
toolbox.register("select", tools.selTournament, tournsize=3)
toolbox.register("mate", gp.cxOnePoint)\
#this is expr_mut, if we want to use a GEP, we can mutute the expr at first, then do the expression
toolbox.register("expr_mut", gp.genFull, min_=0, max_=2)
toolbox.register("mutate", gp.mutUniform, expr=toolbox.expr_mut, pset=pset)
#decorating the operator including crossover and mutate, restricting the tree's height and length
toolbox.decorate("mate", gp.staticLimit(key=operator.attrgetter("height"), max_value=17))
toolbox.decorate("mutate", gp.staticLimit(key=operator.attrgetter("height"), max_value=17))
def main():
#Parameter setting
CXPB = 0.1
MUPB = 0.1
GEN = 400
POP_SIZE = 300
#initializing population
pop = toolbox.population(n = POP_SIZE)
'''start evolution'''
print('Start evolution')
#evaluating the fitness
fitnesses = list(map(toolbox.evaluate, pop))
#assign fitness values
for ind, fit in zip(pop, fitnesses):
ind.fitness.values = fit
print('Evaluated %i individuals' %len(pop))
'''The genetic operations'''
for g in range(GEN):
#select
offspring = toolbox.select(pop, len(pop))
offspring = list(map(toolbox.clone, offspring))
#crossover
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
#mutation
for mutant in offspring:
if random.random() < MUPB:
toolbox.mutate(mutant)
del mutant.fitness.values
#evaluate the invalid individuals
invalid_ind = [ind for ind in offspring if not ind.fitness.valid]
fitnesses = list(map(toolbox.evaluate, invalid_ind))
for ind, fit in zip(invalid_ind, fitnesses):
ind.fitness.values = fit
print('Evaluated %i individuals' %len(invalid_ind))
#update the pop
pop[:] = offspring
#statistics
stats = tools.Statistics(key=lambda ind: ind.fitness.values)
record = stats.compile(pop)
stats.register("avg", numpy.mean, axis=0)
stats.register("min", numpy.min, axis=0)
stats.register("max", numpy.max, axis=0)
record = stats.compile(pop)
logbook = tools.Logbook()
logbook.record(gen=g, evals=30, **record)
logbook.header = "gen", "avg", "min", "max"
print(logbook)
print("-- End of (successful) evolution --")
best_ind = tools.selBest(pop, 1)[0]
print("Best individual is %s, Best results are %s" % (best_ind, best_ind.fitness.values))
if __name__ == "__main__":
main()
大部分注釋都有寫三妈，總結(jié)一下這個(gè)例子：
首先，要發(fā)現(xiàn)的函數(shù)是: f(x) = x^4 - x^3 - x^2 - x氢妈， fitness func為 min{sum(func(x) - f(x))^2}/xnum，x in (-1, 1] 間隔0.1蜕着。
接下來首先搭建GP框架
STEP1. 適應(yīng)度函數(shù)的構(gòu)建酬凳，個(gè)體及種群初始化惠况，所用工具包：creator/toolbox
STEP2. 演化操作定義以及樹的size限制
STEP3. 設(shè)定參數(shù)(交叉率、變異率宁仔、迭代次數(shù)以及種群大小)稠屠，初始化種群并計(jì)算適應(yīng)度
STEP4. 按照選擇、交叉翎苫、變異的順序進(jìn)行演化操作权埠，在一次演化操作完成后選出無效個(gè)體(即適應(yīng)性不如先前個(gè)體的個(gè)體，對(duì)這些個(gè)體的適應(yīng)度不予修改)
STEP5. 更新種群煎谍，統(tǒng)計(jì)當(dāng)前種群的適應(yīng)度信息并輸出
STEP6. 當(dāng)循環(huán)結(jié)束攘蔽，輸出全局最優(yōu)解
注意：在GP這種樹結(jié)構(gòu)的算法中，在適應(yīng)度計(jì)算當(dāng)中呐粘，生成的樹結(jié)構(gòu)首先要通過toolbox.compile轉(zhuǎn)化為函數(shù)满俗，再通過函數(shù)進(jìn)行計(jì)算，同時(shí)所有的遺傳操作均是基于樹結(jié)構(gòu)進(jìn)行的作岖，因此在register這些步驟時(shí)候需要使用特定的方法漫雷，如gp.cxOnePoint等。
