Symbolic Regression

Basic Version

import math
import operator

import numpy
from deap_er import Fitness, Toolbox, creator, gp, tools

tools.rng.seed(1234)  # disables randomization


def safe_div(left, right):
    try:
        return left / right
    except ZeroDivisionError:
        return 1


def evaluate(individual, toolbox, points):
    func = toolbox.compile(expr=individual)
    sq_errors = ((func(x) - x**4 - x**3 - x**2 - x) ** 2 for x in points)
    result = math.fsum(sq_errors) / len(points)
    return (result,)  # The comma is essential here.


def setup():
    pset = gp.PrimitiveSet("MAIN", 1)
    pset.add_primitive(operator.add, 2)
    pset.add_primitive(operator.sub, 2)
    pset.add_primitive(operator.mul, 2)
    pset.add_primitive(safe_div, 2)
    pset.add_primitive(operator.neg, 1)
    pset.add_primitive(math.cos, 1)
    pset.add_primitive(math.sin, 1)
    pset.add_ephemeral_constant("rand101", lambda: tools.rng.randint(-1, 1))
    pset.rename_arguments(ARG0="x")

    creator.create_type("FitnessMin", Fitness, weights=(-1.0,))
    creator.create_type("Individual", gp.PrimitiveTree, fitness=creator.FitnessMin)

    toolbox = Toolbox()
    toolbox.register("expr", gp.gen_half_and_half, prim_set=pset, min_depth=1, max_depth=2)
    toolbox.register("individual", tools.init_iterate, creator.Individual, toolbox.expr)
    toolbox.register("population", tools.init_repeat, list, toolbox.individual)
    toolbox.register("compile", gp.compile_tree, prim_set=pset)
    toolbox.register("mate", gp.cx_one_point)
    toolbox.register("expr_mut", gp.gen_full, min_depth=0, max_depth=2)
    toolbox.register("mutate", gp.mut_uniform, expr=toolbox.expr_mut, prim_set=pset)
    toolbox.register("select", tools.sel_tournament, contestants=3)
    toolbox.register(
        "evaluate", evaluate, toolbox=toolbox, points=[x / 10.0 for x in range(-10, 10)]
    )
    toolbox.decorate("mate", gp.static_limit(limiter=operator.attrgetter("height"), max_value=17))
    toolbox.decorate("mutate", gp.static_limit(limiter=operator.attrgetter("height"), max_value=17))

    stats_size = tools.Statistics(len)
    stats_fit = tools.Statistics(lambda ind: ind.fitness.values)
    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)

    return toolbox, mstats


def print_results(best_ind):
    if best_ind.fitness.values >= (1.0e-3,):
        raise RuntimeError("Evolution failed to converge.")
    print("\nEvolution converged correctly.")


def main():
    toolbox, mstats = setup()
    pop = toolbox.population(size=300)
    hof = tools.HallOfFame(1)
    args = {
        "toolbox": toolbox,
        "population": pop,
        "generations": 40,
        "cx_prob": 0.5,
        "mut_prob": 0.1,
        "hof": hof,
        "stats": mstats,
        "verbose": True,  # prints stats
    }
    tools.ea_simple(**args)
    print_results(hof[0])


if __name__ == "__main__":
    main()

Using Numpy

import numpy
from deap_er import Fitness, Toolbox, creator, gp, tools

tools.rng.seed(1234)  # disables randomization


def safe_div(left, right):
    with numpy.errstate(divide="ignore", invalid="ignore"):
        x = numpy.divide(left, right)
        if isinstance(x, numpy.ndarray):
            x[numpy.isinf(x)] = 1
            x[numpy.isnan(x)] = 1
        elif numpy.isinf(x) or numpy.isnan(x):
            x = 1
    return x


def evaluate(individual, toolbox, samples, values):
    func = toolbox.compile(expr=individual)
    diff = numpy.sum((func(samples) - values) ** 2)
    return (diff,)


def setup():
    samples = numpy.linspace(-1, 1, 10000)
    values = samples**4 + samples**3 + samples**2 + samples

    pset = gp.PrimitiveSet("MAIN", 1)
    pset.add_primitive(numpy.add, 2, name="vadd")
    pset.add_primitive(numpy.subtract, 2, name="vsub")
    pset.add_primitive(numpy.multiply, 2, name="vmul")
    pset.add_primitive(safe_div, 2)
    pset.add_primitive(numpy.negative, 1, name="vneg")
    pset.add_primitive(numpy.cos, 1, name="vcos")
    pset.add_primitive(numpy.sin, 1, name="vsin")
    pset.add_ephemeral_constant("rand101", lambda: tools.rng.randint(-1, 1))
    pset.rename_arguments(ARG0="x")

    creator.create_type("FitnessMin", Fitness, weights=(-1.0,))
    creator.create_type("Individual", gp.PrimitiveTree, fitness=creator.FitnessMin)

    toolbox = Toolbox()
    toolbox.register("expr", gp.gen_half_and_half, prim_set=pset, min_depth=1, max_depth=2)
    toolbox.register("individual", tools.init_iterate, creator.Individual, toolbox.expr)
    toolbox.register("population", tools.init_repeat, list, toolbox.individual)
    toolbox.register("compile", gp.compile_tree, prim_set=pset)
    toolbox.register("evaluate", evaluate, toolbox=toolbox, samples=samples, values=values)
    toolbox.register("select", tools.sel_tournament, contestants=3)
    toolbox.register("mate", gp.cx_one_point)
    toolbox.register("expr_mut", gp.gen_full, min_depth=0, max_depth=2)
    toolbox.register("mutate", gp.mut_uniform, expr=toolbox.expr_mut, prim_set=pset)

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

    return toolbox, stats


def print_results(best_ind):
    if best_ind.fitness.values >= (1.0e-3,):
        raise RuntimeError("Evolution failed to converge.")
    print("\nEvolution converged correctly.")


def main():
    toolbox, stats = setup()
    pop = toolbox.population(size=300)
    hof = tools.HallOfFame(1)
    args = {
        "toolbox": toolbox,
        "population": pop,
        "generations": 40,
        "cx_prob": 0.5,
        "mut_prob": 0.1,
        "hof": hof,
        "stats": stats,
        "verbose": True,  # prints stats
    }
    tools.ea_simple(**args)
    print_results(hof[0])


if __name__ == "__main__":
    main()

Using ADF-s

import math
import operator

import numpy
from deap_er import Fitness, Toolbox, creator, gp, tools

tools.rng.seed(1234)  # disables randomization

CX_PROB = 0.5
MUT_PROB = 0.2
GENS = 40


def safe_div(left, right):
    try:
        return left / right
    except ZeroDivisionError:
        return 1


def evaluate(individual, toolbox, points):
    func = toolbox.compile(expr=individual)
    sq_errors = ((func(x) - x**4 - x**3 - x**2 - x) ** 2 for x in points)
    result = math.fsum(sq_errors) / len(points)
    return (result,)  # The comma is essential here.


def add_primitives(prim_set):
    prim_set.add_primitive(operator.add, 2)
    prim_set.add_primitive(operator.sub, 2)
    prim_set.add_primitive(operator.mul, 2)
    prim_set.add_primitive(safe_div, 2)
    prim_set.add_primitive(operator.neg, 1)
    prim_set.add_primitive(math.cos, 1)
    prim_set.add_primitive(math.sin, 1)


def setup():
    adf_set_2 = gp.PrimitiveSet("ADF2", 2)
    add_primitives(adf_set_2)

    adf_set_1 = gp.PrimitiveSet("ADF1", 2)
    add_primitives(adf_set_1)
    adf_set_1.add_adf(adf_set_2)

    adf_set_0 = gp.PrimitiveSet("ADF0", 2)
    add_primitives(adf_set_0)
    adf_set_0.add_adf(adf_set_1)
    adf_set_0.add_adf(adf_set_2)

    pset = gp.PrimitiveSet("MAIN", 1)
    add_primitives(pset)
    pset.add_ephemeral_constant("rand101", lambda: tools.rng.randint(-1, 1))
    pset.add_adf(adf_set_0)
    pset.add_adf(adf_set_1)
    pset.add_adf(adf_set_2)
    pset.rename_arguments(ARG0="x")

    prim_sets = (pset, adf_set_0, adf_set_1, adf_set_2)

    creator.create_type("FitnessMin", Fitness, weights=(-1.0,))
    creator.create_type("Tree", gp.PrimitiveTree)
    creator.create_type("Individual", list, fitness=creator.FitnessMin)

    toolbox = Toolbox()
    toolbox.register("adf_expr_0", gp.gen_full, prim_set=adf_set_0, min_depth=1, max_depth=2)
    toolbox.register("adf_expr_1", gp.gen_full, prim_set=adf_set_1, min_depth=1, max_depth=2)
    toolbox.register("adf_expr_2", gp.gen_full, prim_set=adf_set_2, min_depth=1, max_depth=2)
    toolbox.register("main_expr", gp.gen_half_and_half, prim_set=pset, min_depth=1, max_depth=2)

    toolbox.register("ADF0", tools.init_iterate, creator.Tree, toolbox.adf_expr_0)
    toolbox.register("ADF1", tools.init_iterate, creator.Tree, toolbox.adf_expr_1)
    toolbox.register("ADF2", tools.init_iterate, creator.Tree, toolbox.adf_expr_2)
    toolbox.register("MAIN", tools.init_iterate, creator.Tree, toolbox.main_expr)

    func_cycle = [toolbox.MAIN, toolbox.ADF0, toolbox.ADF1, toolbox.ADF2]

    toolbox.register("individual", tools.init_cycle, creator.Individual, func_cycle)
    toolbox.register("population", tools.init_repeat, list, toolbox.individual)
    toolbox.register("compile", gp.compile_adf_tree, prim_sets=prim_sets)
    toolbox.register("mate", gp.cx_one_point)
    toolbox.register("expr", gp.gen_full, min_depth=1, max_depth=2)
    toolbox.register("mutate", gp.mut_uniform, expr=toolbox.expr)
    toolbox.register("select", tools.sel_tournament, contestants=3)
    toolbox.register(
        "evaluate", evaluate, toolbox=toolbox, points=[x / 10.0 for x in range(-10, 10)]
    )
    toolbox.decorate("mate", gp.static_limit(limiter=operator.attrgetter("height"), max_value=17))
    toolbox.decorate("mutate", gp.static_limit(limiter=operator.attrgetter("height"), max_value=17))

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

    logbook = tools.Logbook()
    logbook.header = "gen", "evals", "std", "min", "avg", "max"

    return toolbox, stats, logbook, prim_sets


def print_results(best_ind):
    if best_ind.fitness.values >= (0.5,):
        raise RuntimeError("Evolution failed to converge.")
    print("\nEvolution converged correctly.")


def _mate_trees(toolbox, offspring):
    for ind1, ind2 in zip(offspring[::2], offspring[1::2], strict=False):
        for tree1, tree2 in zip(ind1, ind2, strict=False):
            if tools.rng.random() < CX_PROB:
                toolbox.mate(tree1, tree2)
                del ind1.fitness.values
                del ind2.fitness.values


def _mutate_trees(toolbox, offspring, psets):
    for ind in offspring:
        for tree, pset in zip(ind, psets, strict=False):
            if tools.rng.random() < MUT_PROB:
                toolbox.mutate(individual=tree, prim_set=pset)
                del ind.fitness.values


def main():
    toolbox, stats, logbook, psets = setup()
    hof = tools.HallOfFame(1)

    def log_stats(ngen=0):
        hof.update(pop)
        record = stats.compile(pop)
        logbook.record(gen=ngen, evals=len(pop), **record)
        print(logbook.stream)

    pop = toolbox.population(size=100)
    tools.evaluate_invalid(toolbox, pop)

    log_stats()

    for generations in range(1, GENS):
        offspring = toolbox.select(pop, len(pop))
        offspring = [toolbox.clone(ind) for ind in offspring]
        _mate_trees(toolbox, offspring)
        _mutate_trees(toolbox, offspring, psets)
        tools.evaluate_invalid(toolbox, offspring)
        pop = offspring
        log_stats(generations)

    print_results(hof[0])


if __name__ == "__main__":
    main()

Using HARM-GP

import math
import operator

import numpy
from deap_er import Fitness, Toolbox, creator, gp, tools

tools.rng.seed(1234)  # disables randomization


def safe_div(left, right):
    try:
        return left / right
    except ZeroDivisionError:
        return 1


def evaluate(individual, points, toolbox):
    func = toolbox.compile(expr=individual)
    sq_errors = ((func(x) - x**4 - x**3 - x**2 - x) ** 2 for x in points)
    result = math.fsum(sq_errors) / len(points)
    return (result,)  # The comma is essential here.


def setup():
    pset = gp.PrimitiveSet("MAIN", 1)
    pset.add_primitive(operator.add, 2)
    pset.add_primitive(operator.sub, 2)
    pset.add_primitive(operator.mul, 2)
    pset.add_primitive(safe_div, 2)
    pset.add_primitive(operator.neg, 1)
    pset.add_primitive(math.cos, 1)
    pset.add_primitive(math.sin, 1)
    pset.add_ephemeral_constant("rand101", lambda: tools.rng.randint(-1, 1))
    pset.rename_arguments(ARG0="x")

    creator.create_type("FitnessMin", Fitness, weights=(-1.0,))
    creator.create_type("Individual", gp.PrimitiveTree, fitness=creator.FitnessMin)

    toolbox = Toolbox()
    toolbox.register("expr", gp.gen_half_and_half, prim_set=pset, min_depth=1, max_depth=2)
    toolbox.register("individual", tools.init_iterate, creator.Individual, toolbox.expr)
    toolbox.register("population", tools.init_repeat, list, toolbox.individual)
    toolbox.register("compile", gp.compile_tree, prim_set=pset)
    toolbox.register("mate", gp.cx_one_point)
    toolbox.register("expr_mut", gp.gen_full, min_depth=0, max_depth=2)
    toolbox.register("mutate", gp.mut_uniform, expr=toolbox.expr_mut, prim_set=pset)
    toolbox.register("select", tools.sel_tournament, contestants=3)
    toolbox.register(
        "evaluate", evaluate, points=[x / 10.0 for x in range(-10, 10)], toolbox=toolbox
    )
    toolbox.decorate("mate", gp.static_limit(limiter=operator.attrgetter("height"), max_value=17))
    toolbox.decorate("mutate", gp.static_limit(limiter=operator.attrgetter("height"), max_value=17))

    stats_size = tools.Statistics(len)
    stats_fit = tools.Statistics(lambda ind: ind.fitness.values)
    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)

    return toolbox, mstats


def print_results(best_ind):
    if best_ind.fitness.values >= (1.0e-3,):
        raise RuntimeError("Evolution failed to converge.")
    print("\nEvolution converged correctly.")


def main():
    toolbox, mstats = setup()
    pop = toolbox.population(size=300)
    hof = tools.HallOfFame(1)
    args = {
        "toolbox": toolbox,
        "population": pop,
        "generations": 40,
        "cx_prob": 0.5,
        "mut_prob": 0.1,
        "hof": hof,
        "stats": mstats,
        "verbose": True,  # prints stats
    }
    gp.harm(**args)
    print_results(hof[0])


if __name__ == "__main__":
    main()

Using SLIM semantic GP

Non-bloating geometric semantic variation with HARM size control:

import math
import operator

import numpy
from deap_er import Fitness, Toolbox, creator, gp, tools

tools.rng.seed(1234)  # disables randomization


def lf(x):
    if x > 0:
        return 1 / (1 + math.exp(-x))
    exp_x = math.exp(x)
    return exp_x / (1 + exp_x)


def safe_div(left, right):
    try:
        return left / right
    except ZeroDivisionError:
        return 1


def evaluate(individual, points, toolbox):
    func = toolbox.compile(individual)
    sq_errors = ((func(x) - x**4 - x**3 - x**2 - x) ** 2 for x in points)
    result = math.fsum(sq_errors) / len(points)
    return (result,)


def setup():
    pset = gp.PrimitiveSet("MAIN", 1)
    pset.add_primitive(operator.add, 2)
    pset.add_primitive(operator.sub, 2)
    pset.add_primitive(operator.mul, 2)
    pset.add_primitive(safe_div, 2)
    pset.add_primitive(operator.neg, 1)
    pset.add_primitive(lf, 1, name="lf")
    pset.add_ephemeral_constant("rand101", lambda: tools.rng.randint(-1, 1))
    pset.rename_arguments(ARG0="x")

    creator.create_type("FitnessMin", Fitness, weights=(-1.0,))
    creator.create_type(
        "Individual",
        gp.SlimTree,
        fitness=creator.FitnessMin,
    )

    toolbox = Toolbox()
    toolbox.register("expr", gp.gen_half_and_half, prim_set=pset, min_depth=1, max_depth=2)
    toolbox.register("individual", tools.init_iterate, creator.Individual, toolbox.expr)
    toolbox.register("population", tools.init_repeat, list, toolbox.individual)
    toolbox.register("compile", gp.compile_slim_tree, prim_set=pset)
    toolbox.register("mate", gp.cx_slim_donor, prim_set=pset)
    toolbox.register("mutate", gp.mut_slim, prim_set=pset, inflate_prob=0.3)
    toolbox.register("select", tools.sel_tournament, contestants=3)
    toolbox.register(
        "evaluate", evaluate, points=[x / 10.0 for x in range(-10, 10)], toolbox=toolbox
    )
    toolbox.register("clone", tools.clone_individual)

    stats_size = tools.Statistics(len)
    stats_fit = tools.Statistics(lambda ind: ind.fitness.values)
    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)

    return toolbox, mstats


def print_results(best_ind, initial_fitness):
    if best_ind.fitness.values[0] >= initial_fitness:
        raise RuntimeError("Evolution failed to improve the population.")
    if best_ind.fitness.values[0] >= 0.1:
        raise RuntimeError("Evolution failed to reach the demo fitness target.")
    if len(best_ind) > 250:
        raise RuntimeError("SLIM elite grew too large for a readable model.")
    print("\nEvolution converged correctly.")


def main():
    toolbox, mstats = setup()
    pop = toolbox.population(size=300)
    hof = tools.HallOfFame(1)
    for ind in pop:
        ind.fitness.values = toolbox.evaluate(ind)
    initial_fitness = min(ind.fitness.values[0] for ind in pop)
    args = {
        "toolbox": toolbox,
        "population": pop,
        "generations": 80,
        "cx_prob": 0.5,
        "mut_prob": 0.2,
        "hof": hof,
        "stats": mstats,
        "verbose": True,
    }
    gp.harm(**args)
    print_results(hof[0], initial_fitness)


if __name__ == "__main__":
    main()

Using Epsilon Lexicase

import math
import operator

import numpy
from deap_er import Fitness, Toolbox, creator, gp, tools

tools.rng.seed(1234)  # disables randomization


def safe_div(left, right):
    try:
        return left / right
    except ZeroDivisionError:
        return 1


def evaluate(individual, toolbox, points):
    func = toolbox.compile(expr=individual)
    sq_errors = ((func(x) - x**4 - x**3 - x**2 - x) ** 2 for x in points)
    result = math.fsum(sq_errors) / len(points)
    return (result,)  # The comma is essential here.


def setup():
    pset = gp.PrimitiveSet("MAIN", 1)
    pset.add_primitive(operator.add, 2)
    pset.add_primitive(operator.sub, 2)
    pset.add_primitive(operator.mul, 2)
    pset.add_primitive(safe_div, 2)
    pset.add_primitive(operator.neg, 1)
    pset.add_primitive(math.cos, 1)
    pset.add_primitive(math.sin, 1)
    pset.add_ephemeral_constant("rand101", lambda: tools.rng.randint(-1, 1))
    pset.rename_arguments(ARG0="x")

    creator.create_type("FitnessMin", Fitness, weights=(-1.0,))
    creator.create_type("Individual", gp.PrimitiveTree, fitness=creator.FitnessMin)

    toolbox = Toolbox()
    toolbox.register("expr", gp.gen_half_and_half, prim_set=pset, min_depth=1, max_depth=2)
    toolbox.register("individual", tools.init_iterate, creator.Individual, toolbox.expr)
    toolbox.register("population", tools.init_repeat, list, toolbox.individual)
    toolbox.register("compile", gp.compile_tree, prim_set=pset)
    toolbox.register("mate", gp.cx_one_point)
    toolbox.register("expr_mut", gp.gen_full, min_depth=0, max_depth=2)
    toolbox.register("mutate", gp.mut_uniform, expr=toolbox.expr_mut, prim_set=pset)
    toolbox.register("select", tools.sel_epsilon_lexicase)
    toolbox.register(
        "evaluate", evaluate, toolbox=toolbox, points=[x / 10.0 for x in range(-10, 10)]
    )
    toolbox.decorate("mate", gp.static_limit(limiter=operator.attrgetter("height"), max_value=17))
    toolbox.decorate("mutate", gp.static_limit(limiter=operator.attrgetter("height"), max_value=17))

    stats_size = tools.Statistics(len)
    stats_fit = tools.Statistics(lambda ind: ind.fitness.values)
    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)

    return toolbox, mstats


def print_results(best_ind):
    if best_ind.fitness.values >= (1.0e-3,):
        raise RuntimeError("Evolution failed to converge.")
    print("\nEvolution converged correctly.")


def main():
    toolbox, mstats = setup()
    pop = toolbox.population(size=300)
    hof = tools.HallOfFame(1)
    args = {
        "toolbox": toolbox,
        "population": pop,
        "generations": 40,
        "cx_prob": 0.5,
        "mut_prob": 0.1,
        "hof": hof,
        "stats": mstats,
        "verbose": True,  # prints stats
    }
    tools.ea_simple(**args)
    print_results(hof[0])


if __name__ == "__main__":
    main()