(六) GP实例——函数发现
话不多说直接上代码

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()

大部分注释都有写，总结一下这个例子：
首先，要发现的函数是: 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. 适应度函数的构建，个体及种群初始化，所用工具包：creator/toolbox
STEP2. 演化操作定义以及树的size限制
STEP3. 设定参数(交叉率、变异率、迭代次数以及种群大小)，初始化种群并计算适应度
STEP4. 按照选择、交叉、变异的顺序进行演化操作，在一次演化操作完成后选出无效个体(即适应性不如先前个体的个体，对这些个体的适应度不予修改)
STEP5. 更新种群，统计当前种群的适应度信息并输出
STEP6. 当循环结束，输出全局最优解

注意：在GP这种树结构的算法中，在适应度计算当中，生成的树结构首先要通过toolbox.compile转化为函数，再通过函数进行计算，同时所有的遗传操作均是基于树结构进行的，因此在register这些步骤时候需要使用特定的方法，如gp.cxOnePoint等。