Example for evolutionary regression with parametrized nodes

Example demonstrating the use of Cartesian genetic programming for a regression task that requires fine tuning of constants in parametrized nodes. This is achieved by introducing a new node, “ParametrizedAdd” which produces a scaled and shifted version of the sum of its inputs.

# The docopt str is added explicitly to ensure compatibility with
# sphinx-gallery.
docopt_str = """
  Usage:
    example_parametrized_nodes.py [--max-generations=<N>]

  Options:
    -h --help
    --max-generations=<N>  Maximum number of generations [default: 500]
"""

import functools
import math

import matplotlib.pyplot as plt
import numpy as np
import scipy.constants
import torch
from docopt import docopt

import cgp

args = docopt(docopt_str)

We first define a new node that adds the values of its two inputs then scales and finally shifts the result. The scale (“w”) and shift factors (“b”) are parameters that are adapted by local search. We need to define the arity of the node, callables for the initial values for the parameters and the operation of the node as a string. In this string parameters are enclosed in angle brackets, inputs are denoted by “x_i” with i representing their corresponding index.

class ParametrizedAdd(cgp.OperatorNode):
    """A node that adds its two inputs.

    The result of addition is scaled by w and shifted by b. Both these
    parameters can be adapted via local search are passed on from
    parents to their offspring.

    """

    _arity = 2
    _initial_values = {"<w>": lambda: 1.0, "<b>": lambda: 0.0}
    _def_output = "<w> * (x_0 + x_1) + <b>"

We define a target function which contains numerical constants that are not available as constants for the search and need to be found by local search on parameterized nodes.

def f_target(x):
    return math.pi * (x[:, 0] + x[:, 1]) + math.e

Then we define a differentiable(!) inner objective function for the evolution. This function accepts a torch class as a parameter. It returns the mean-squared error between the output of the forward pass of this class and the target function evaluated on a set of random points. This inner objective is then used by actual objective function to determine the fitness of the individual.

def inner_objective(f, seed):
    torch.manual_seed(seed)
    batch_size = 500
    x = torch.DoubleTensor(batch_size, 2).uniform_(-5, 5)
    y = f(x)
    return torch.nn.MSELoss()(f_target(x), y[:, 0])


def objective(individual, seed):
    if not individual.fitness_is_None():
        return individual

    f = individual.to_torch()
    loss = inner_objective(f, seed)
    individual.fitness = -loss.item()

    return individual

Next, we define the parameters for the population, the genome of individuals, the evolutionary algorithm, and the local search. Note that we add the custom node defined above as a primitive.

population_params = {"n_parents": 1, "seed": 818821}

genome_params = {
    "n_inputs": 2,
    "n_outputs": 1,
    "n_columns": 5,
    "n_rows": 1,
    "levels_back": None,
    "primitives": (ParametrizedAdd, cgp.Add, cgp.Sub, cgp.Mul),
}

ea_params = {"n_offsprings": 4, "tournament_size": 1, "mutation_rate": 0.04, "n_processes": 2}

evolve_params = {"max_generations": int(args["--max-generations"]), "termination_fitness": 0.0}

local_search_params = {"lr": 1e-3, "gradient_steps": 9}

We then create a Population instance and instantiate the local search and evolutionary algorithm.

pop = cgp.Population(**population_params, genome_params=genome_params)

local_search = functools.partial(
    cgp.local_search.gradient_based,
    objective=functools.partial(inner_objective, seed=population_params["seed"]),
    **local_search_params,
)

ea = cgp.ea.MuPlusLambda(**ea_params, local_search=local_search)

We define a recording callback closure for bookkeeping of the progression of the evolution.

history = {}
history["fitness_champion"] = []
history["expr_champion"] = []


def recording_callback(pop):
    history["fitness_champion"].append(pop.champion.fitness)
    history["expr_champion"].append(pop.champion.to_sympy())

We fix the seed for the objective function to make sure results are comparable across individuals and, finally, we call the evolve method to perform the evolutionary search.

obj = functools.partial(objective, seed=population_params["seed"])

cgp.evolve(pop, obj, ea, **evolve_params, print_progress=True, callback=recording_callback)

Out:

[2/500] max fitness: -127.01931908693929
[3/500] max fitness: -127.01931908693929
[4/500] max fitness: -127.01931908693929
[5/500] max fitness: -5.7621478742860965
[6/500] max fitness: -4.400266790086014
[7/500] max fitness: -3.3439171938076777
[8/500] max fitness: -2.5221616375748344
[9/500] max fitness: -1.886291288388642
[10/500] max fitness: -1.3981594801407393
[11/500] max fitness: -1.027024955717336
[12/500] max fitness: -0.7477863327318205
[13/500] max fitness: -0.5399321058490655
[14/500] max fitness: -0.3868341292604577
[15/500] max fitness: -0.2751900492412008
[16/500] max fitness: -0.19452652473250845
[17/500] max fitness: -0.13673494099465527
[18/500] max fitness: -0.09563997095003499
[19/500] max fitness: -0.06661059705457541
[20/500] max fitness: -0.046222253258311974
[21/500] max fitness: -0.031973909639731866
[22/500] max fitness: -0.02205885115609153
[23/500] max fitness: -0.01518420142810349
[24/500] max fitness: -0.010432247679957566
[25/500] max fitness: -0.007156051160756987
[26/500] max fitness: -0.004902206168242835
[27/500] max fitness: -0.003354500929279119
[28/500] max fitness: -0.0022933114582781846
[29/500] max fitness: -0.0015666233077337203
[30/500] max fitness: -0.0010695205965242655
[31/500] max fitness: -0.0007297666386346675
[32/500] max fitness: -0.0004977239202683268
[33/500] max fitness: -0.00033934028308000566
[34/500] max fitness: -0.00023128736715539105
[35/500] max fitness: -0.0001576015752678627
[36/500] max fitness: -0.00010736927896605433
[37/500] max fitness: -7.313509670770486e-05
[38/500] max fitness: -4.980933852908367e-05
[39/500] max fitness: -3.391918703676957e-05
[40/500] max fitness: -2.3096096228535896e-05
[41/500] max fitness: -1.5725246069023054e-05
[42/500] max fitness: -1.0706019354335675e-05
[43/500] max fitness: -7.288451067003709e-06
[44/500] max fitness: -4.961616392292634e-06
[45/500] max fitness: -3.3774988989541173e-06
[46/500] max fitness: -2.2990801926213286e-06
[47/500] max fitness: -1.564956210506372e-06
[48/500] max fitness: -1.0652249266356427e-06
[49/500] max fitness: -7.25058544542861e-07
[50/500] max fitness: -4.935131658081641e-07
[51/500] max fitness: -3.359072595462748e-07
[52/500] max fitness: -2.2863140927698259e-07
[53/500] max fitness: -1.556141119044658e-07
[54/500] max fitness: -1.059154323848727e-07
[55/500] max fitness: -7.208869798177893e-08
[56/500] max fitness: -4.906515782366229e-08
[57/500] max fitness: -3.33947068058455e-08
[58/500] max fitness: -2.2729022729656698e-08
[59/500] max fitness: -1.5469732821737942e-08
[60/500] max fitness: -1.0528923586553802e-08
[61/500] max fitness: -7.166125354812942e-09
[62/500] max fitness: -4.877353378342053e-09
[63/500] max fitness: -3.31958310056805e-09
[64/500] max fitness: -2.2593444953407397e-09
[65/500] max fitness: -1.5377333234503386e-09
[66/500] max fitness: -1.0465965912643985e-09
[67/500] max fitness: -7.123236680097046e-10
[68/500] max fitness: -4.848141005422031e-10
[69/500] max fitness: -3.299688575234535e-10
[70/500] max fitness: -2.2457971871821608e-10
[71/500] max fitness: -1.5285090235237916e-10
[72/500] max fitness: -1.0403162728212182e-10
[73/500] max fitness: -7.080480056918192e-11

After the evolutionary search has ended, we print the expression with the highest fitness and plot the progression of the search.

print(f"Final expression {pop.champion.to_sympy()} with fitness {pop.champion.fitness}")

print("Best performing expression per generation (for fitness increase > 0.5):")
old_fitness = -np.inf
for i, (fitness, expr) in enumerate(zip(history["fitness_champion"], history["expr_champion"])):
    delta_fitness = fitness - old_fitness
    if delta_fitness > 0.5:
        print(f"{i:3d}: {fitness}, {expr}")
        old_fitness = fitness
print(f"{i:3d}: {fitness}, {expr}")

width = 9.0

fig = plt.figure(figsize=(width, width / scipy.constants.golden))

ax_fitness = fig.add_subplot(111)
ax_fitness.set_xlabel("Generation")
ax_fitness.set_ylabel("Fitness")
ax_fitness.set_yscale("symlog")

ax_fitness.axhline(0.0, color="k")
ax_fitness.plot(history["fitness_champion"], lw=2)

plt.savefig("example_parametrized_nodes.pdf", dpi=300)

Out:

Final expression 3.1415928340353736*x_0 + 3.1415928340353736*x_1 + 2.7182734383513227 with fitness -7.080480056918192e-11
Best performing expression per generation (for fitness increase > 0.5):
  0: -127.01931908693929, x_1
  4: -5.7621478742860965, 3.1337934090450837*x_0 + 3.1337934090450837*x_1 + 0.31836625198458554
  5: -4.400266790086014, 3.1565329832083608*x_0 + 3.1565329832083608*x_1 + 0.62085799830068183
  6: -3.3439171938076777, 3.15892017999446*x_0 + 3.15892017999446*x_1 + 0.89027846015470884
  7: -2.5221616375748344, 3.1597540007392134*x_0 + 3.1597540007392134*x_1 + 1.1311091306962316
  8: -1.886291288388642, 3.1595918340334699*x_0 + 3.1595918340334699*x_1 + 1.3460568187315501
 10: -1.027024955717336, 3.1575106310584499*x_0 + 3.1575106310584499*x_1 + 1.7062610297769132
 13: -0.3868341292604577, 3.1529224955139751*x_0 + 3.1529224955139751*x_1 + 2.0975544857381176
 72: -7.080480056918192e-11, 3.1415928340353736*x_0 + 3.1415928340353736*x_1 + 2.7182734383513227