Knuth Problem (described in section 3.2 of the 3rd edition of the book "Artificial Intelligence: A Modern Approach", p. 73) with Tree-based Genetic Programming

knuth_problem.py

# You need to have graphviz installed 
# if you want to plot the trees (commented for now).
# It doesn't seem to find the right solution. 
# The best solution seems to be: floor(sqrt(float(sqrt(4)))),
# which produces 1, so it has a fitness of 4.

import math
import operator
import random
import sys

# import matplotlib.pyplot as plt
# import networkx as nx
import numpy as np
from deap import algorithms, base, creator, tools
from deap import gp

# from networkx.drawing.nx_agraph import graphviz_layout

NoneType = type(None)


# def plot_expression(expr):
#     nodes, edges, labels = gp.graph(expr)
#
#     g = nx.Graph()
#     g.add_nodes_from(nodes)
#     g.add_edges_from(edges)
#
#     pos = graphviz_layout(g, prog="dot")
#
#     nx.draw_networkx_nodes(g, pos)
#     nx.draw_networkx_edges(g, pos)
#     nx.draw_networkx_labels(g, pos, labels)
#
#     plt.show()


def fact(x):
    return float(math.factorial(x))


def identity(x):
    return x


primitive_set = gp.PrimitiveSetTyped("f", [int], float)
primitive_set.addPrimitive(math.floor, [float], float)
primitive_set.addPrimitive(math.sqrt, [float], float)
primitive_set.addPrimitive(abs, [float], float)
primitive_set.addPrimitive(abs, [int], int)

primitive_set.addPrimitive(fact, [int], float)
primitive_set.addPrimitive(identity, [int], int, name="int")
primitive_set.addPrimitive(identity, [float], float, name="float")

primitive_set.addTerminal(4, int)
primitive_set.addTerminal(4, float)

primitive_set.renameArguments(ARG0='x')

# Define new classes, whose objects can later be initialized and used.
creator.create("FitnessMax", base.Fitness, weights=(1.0,))
creator.create("Individual", gp.PrimitiveTree, fitness=creator.FitnessMax)

# Create the toolbox.
toolbox = base.Toolbox()

# Define function that can later be called during the evolution process.
toolbox.register("tree_initializer", gp.genFull, pset=primitive_set,
                 min_=1, max_=5)

toolbox.register("individual", tools.initIterate, creator.Individual,
                 toolbox.tree_initializer)
toolbox.register("population", tools.initRepeat, list, toolbox.individual)
toolbox.register("compile", gp.compile, pset=primitive_set)


def evaluate(individual):
    func = toolbox.compile(expr=individual)
    print(individual)
    # plot_expression(individual)
    diff = 5 - func(4)
    return diff,


toolbox.register("evaluate", evaluate)
toolbox.register("select", tools.selTournament, tournsize=3)
toolbox.register("mate", gp.cxOnePoint)
toolbox.register("expr_mut", gp.genFull, min_=0, max_=3)
toolbox.register("mutate", gp.mutUniform, expr=toolbox.expr_mut,
                 pset=primitive_set)

toolbox.decorate("mate", gp.staticLimit(key=operator.attrgetter("height"),
                                        max_value=30))
toolbox.decorate("mutate", gp.staticLimit(key=operator.attrgetter("height"),
                                          max_value=30))


def main():
    seed = random.randrange(sys.maxsize)
    random.seed(seed)

    #  2543404381583366948
    print("Seed =", seed)

    pop = toolbox.population(n=100)
    hof = tools.HallOfFame(1)

    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", np.mean)

    pop, log = algorithms.eaSimple(pop, toolbox, 0.8, 0.2, 100, stats=mstats,
                                   halloffame=hof, verbose=True)
    # print log
    return pop, log, hof


if __name__ == "__main__":
    pop, log, hof = main()
    print("Best individual =", hof[0])
    print("Fitness =", hof[0].fitness)