Symbolic Regression

Basic Version

 1from deap_er import creator
 2from deap_er import tools
 3from deap_er import base
 4from deap_er import gp
 5import operator
 6import random
 7import numpy
 8import math
 9
10
11random.seed(1234)  # disables randomization
12
13
14def safe_div(left, right):
15    try:
16        return left / right
17    except ZeroDivisionError:
18        return 1
19
20
21def evaluate(individual, toolbox, points):
22    func = toolbox.compile(expr=individual)
23    sq_errors = ((func(x) - x**4 - x**3 - x**2 - x)**2 for x in points)
24    result = math.fsum(sq_errors) / len(points)
25    return result,  # The comma is essential here.
26
27
28def setup():
29    pset = gp.PrimitiveSet("MAIN", 1)
30    pset.add_primitive(operator.add, 2)
31    pset.add_primitive(operator.sub, 2)
32    pset.add_primitive(operator.mul, 2)
33    pset.add_primitive(safe_div, 2)
34    pset.add_primitive(operator.neg, 1)
35    pset.add_primitive(math.cos, 1)
36    pset.add_primitive(math.sin, 1)
37    pset.add_ephemeral_constant("rand101", lambda: random.randint(-1, 1))
38    pset.rename_arguments(ARG0='x')
39
40    creator.create("FitnessMin", base.Fitness, weights=(-1.0,))
41    creator.create("Individual", gp.PrimitiveTree, fitness=creator.FitnessMin)
42
43    toolbox = base.Toolbox()
44    toolbox.register("expr", gp.gen_half_and_half, prim_set=pset, min_depth=1, max_depth=2)
45    toolbox.register("individual", tools.init_iterate, creator.Individual, toolbox.expr)
46    toolbox.register("population", tools.init_repeat, list, toolbox.individual)
47    toolbox.register("compile", gp.compile_tree, prim_set=pset)
48    toolbox.register("mate", gp.cx_one_point)
49    toolbox.register("expr_mut", gp.gen_full, min_depth=0, max_depth=2)
50    toolbox.register("mutate", gp.mut_uniform, expr=toolbox.expr_mut, prim_set=pset)
51    toolbox.register("select", tools.sel_tournament, contestants=3)
52    toolbox.register("evaluate", evaluate, toolbox=toolbox, points=[x / 10. for x in range(-10, 10)])
53    toolbox.decorate("mate", gp.static_limit(limiter=operator.attrgetter("height"), max_value=17))
54    toolbox.decorate("mutate", gp.static_limit(limiter=operator.attrgetter("height"), max_value=17))
55
56    stats_size = tools.Statistics(len)
57    stats_fit = tools.Statistics(lambda ind: ind.fitness.values)
58    mstats = tools.MultiStatistics(fitness=stats_fit, size=stats_size)
59    mstats.register("avg", numpy.mean)
60    mstats.register("std", numpy.std)
61    mstats.register("min", numpy.min)
62    mstats.register("max", numpy.max)
63
64    return toolbox, mstats
65
66
67def print_results(best_ind):
68    if not best_ind.fitness.values < (1.0e-3,):
69        raise RuntimeError('Evolution failed to converge.')
70    print(f'\nEvolution converged correctly.')
71
72
73def main():
74    toolbox, mstats = setup()
75    pop = toolbox.population(size=300)
76    hof = tools.HallOfFame(1)
77    args = dict(
78        toolbox=toolbox,
79        population=pop,
80        generations=40,
81        cx_prob=0.5,
82        mut_prob=0.1,
83        hof=hof,
84        stats=mstats,
85        verbose=True  # prints stats
86    )
87    tools.ea_simple(**args)
88    print_results(hof[0])
89
90
91if __name__ == "__main__":
92    main()

Using Numpy

 1from deap_er import creator
 2from deap_er import tools
 3from deap_er import base
 4from deap_er import gp
 5import random
 6import numpy
 7
 8
 9random.seed(1234)  # disables randomization
10
11
12def safe_div(left, right):
13    with numpy.errstate(divide='ignore', invalid='ignore'):
14        x = numpy.divide(left, right)
15        if isinstance(x, numpy.ndarray):
16            x[numpy.isinf(x)] = 1
17            x[numpy.isnan(x)] = 1
18        elif numpy.isinf(x) or numpy.isnan(x):
19            x = 1
20    return x
21
22
23def evaluate(individual, toolbox, samples, values):
24    func = toolbox.compile(expr=individual)
25    diff = numpy.sum((func(samples) - values)**2)
26    return diff,
27
28
29def setup():
30    samples = numpy.linspace(-1, 1, 10000)
31    values = samples ** 4 + samples ** 3 + samples ** 2 + samples
32
33    pset = gp.PrimitiveSet("MAIN", 1)
34    pset.add_primitive(numpy.add, 2, name="vadd")
35    pset.add_primitive(numpy.subtract, 2, name="vsub")
36    pset.add_primitive(numpy.multiply, 2, name="vmul")
37    pset.add_primitive(safe_div, 2)
38    pset.add_primitive(numpy.negative, 1, name="vneg")
39    pset.add_primitive(numpy.cos, 1, name="vcos")
40    pset.add_primitive(numpy.sin, 1, name="vsin")
41    pset.add_ephemeral_constant("rand101", lambda: random.randint(-1, 1))
42    pset.rename_arguments(ARG0='x')
43
44    creator.create("FitnessMin", base.Fitness, weights=(-1.0,))
45    creator.create("Individual", gp.PrimitiveTree, fitness=creator.FitnessMin)
46
47    toolbox = base.Toolbox()
48    toolbox.register("expr", gp.gen_half_and_half, prim_set=pset, min_depth=1, max_depth=2)
49    toolbox.register("individual", tools.init_iterate, creator.Individual, toolbox.expr)
50    toolbox.register("population", tools.init_repeat, list, toolbox.individual)
51    toolbox.register("compile", gp.compile_tree, prim_set=pset)
52    toolbox.register("evaluate", evaluate, toolbox=toolbox, samples=samples, values=values)
53    toolbox.register("select", tools.sel_tournament, contestants=3)
54    toolbox.register("mate", gp.cx_one_point)
55    toolbox.register("expr_mut", gp.gen_full, min_depth=0, max_depth=2)
56    toolbox.register("mutate", gp.mut_uniform, expr=toolbox.expr_mut, prim_set=pset)
57
58    stats = tools.Statistics(lambda ind: ind.fitness.values)
59    stats.register("avg", numpy.mean)
60    stats.register("std", numpy.std)
61    stats.register("min", numpy.min)
62    stats.register("max", numpy.max)
63
64    return toolbox, stats
65
66
67def print_results(best_ind):
68    if not best_ind.fitness.values < (1.0e-3,):
69        raise RuntimeError('Evolution failed to converge.')
70    print(f'\nEvolution converged correctly.')
71
72
73def main():
74    toolbox, stats = setup()
75    pop = toolbox.population(size=300)
76    hof = tools.HallOfFame(1)
77    args = dict(
78        toolbox=toolbox,
79        population=pop,
80        generations=40,
81        cx_prob=0.5,
82        mut_prob=0.1,
83        hof=hof,
84        stats=stats,
85        verbose=True  # prints stats
86    )
87    tools.ea_simple(**args)
88    print_results(hof[0])
89
90
91if __name__ == "__main__":
92    main()

Using ADF-s

  1from deap_er import creator
  2from deap_er import tools
  3from deap_er import base
  4from deap_er import gp
  5import operator
  6import random
  7import numpy
  8import math
  9
 10
 11random.seed(1234)  # disables randomization
 12
 13CX_PROB = 0.5
 14MUT_PROB = 0.2
 15GENS = 40
 16
 17
 18def safe_div(left, right):
 19    try:
 20        return left / right
 21    except ZeroDivisionError:
 22        return 1
 23
 24
 25def evaluate(individual, toolbox, points):
 26    func = toolbox.compile(expr=individual)
 27    sq_errors = ((func(x) - x**4 - x**3 - x**2 - x)**2 for x in points)
 28    result = math.fsum(sq_errors) / len(points)
 29    return result,  # The comma is essential here.
 30
 31
 32def add_primitives(prim_set):
 33    prim_set.add_primitive(operator.add, 2)
 34    prim_set.add_primitive(operator.sub, 2)
 35    prim_set.add_primitive(operator.mul, 2)
 36    prim_set.add_primitive(safe_div, 2)
 37    prim_set.add_primitive(operator.neg, 1)
 38    prim_set.add_primitive(math.cos, 1)
 39    prim_set.add_primitive(math.sin, 1)
 40
 41
 42def setup():
 43    adf_set_2 = gp.PrimitiveSet("ADF2", 2)
 44    add_primitives(adf_set_2)
 45
 46    adf_set_1 = gp.PrimitiveSet("ADF1", 2)
 47    add_primitives(adf_set_1)
 48    adf_set_1.add_adf(adf_set_2)
 49
 50    adf_set_0 = gp.PrimitiveSet("ADF0", 2)
 51    add_primitives(adf_set_0)
 52    adf_set_0.add_adf(adf_set_1)
 53    adf_set_0.add_adf(adf_set_2)
 54
 55    pset = gp.PrimitiveSet("MAIN", 1)
 56    add_primitives(pset)
 57    pset.add_ephemeral_constant("rand101", lambda: random.randint(-1, 1))
 58    pset.add_adf(adf_set_0)
 59    pset.add_adf(adf_set_1)
 60    pset.add_adf(adf_set_2)
 61    pset.rename_arguments(ARG0='x')
 62
 63    prim_sets = (pset, adf_set_0, adf_set_1, adf_set_2)
 64
 65    creator.create("FitnessMin", base.Fitness, weights=(-1.0,))
 66    creator.create("Tree", gp.PrimitiveTree)
 67    creator.create("Individual", list, fitness=creator.FitnessMin)
 68
 69    toolbox = base.Toolbox()
 70    toolbox.register('adf_expr_0', gp.gen_full, prim_set=adf_set_0, min_depth=1, max_depth=2)
 71    toolbox.register('adf_expr_1', gp.gen_full, prim_set=adf_set_1, min_depth=1, max_depth=2)
 72    toolbox.register('adf_expr_2', gp.gen_full, prim_set=adf_set_2, min_depth=1, max_depth=2)
 73    toolbox.register('main_expr', gp.gen_half_and_half, prim_set=pset, min_depth=1, max_depth=2)
 74
 75    toolbox.register('ADF0', tools.init_iterate, creator.Tree, toolbox.adf_expr_0)
 76    toolbox.register('ADF1', tools.init_iterate, creator.Tree, toolbox.adf_expr_1)
 77    toolbox.register('ADF2', tools.init_iterate, creator.Tree, toolbox.adf_expr_2)
 78    toolbox.register('MAIN', tools.init_iterate, creator.Tree, toolbox.main_expr)
 79
 80    func_cycle = [toolbox.MAIN, toolbox.ADF0, toolbox.ADF1, toolbox.ADF2]
 81
 82    toolbox.register('individual', tools.init_cycle, creator.Individual, func_cycle)
 83    toolbox.register('population', tools.init_repeat, list, toolbox.individual)
 84    toolbox.register('compile', gp.compile_adf_tree, prim_sets=prim_sets)
 85    toolbox.register('mate', gp.cx_one_point)
 86    toolbox.register('expr', gp.gen_full, min_depth=1, max_depth=2)
 87    toolbox.register('mutate', gp.mut_uniform, expr=toolbox.expr)
 88    toolbox.register('select', tools.sel_tournament, contestants=3)
 89    toolbox.register('evaluate', evaluate, toolbox=toolbox, points=[x / 10. for x in range(-10, 10)])
 90    toolbox.decorate("mate", gp.static_limit(limiter=operator.attrgetter("height"), max_value=17))
 91    toolbox.decorate("mutate", gp.static_limit(limiter=operator.attrgetter("height"), max_value=17))
 92
 93    stats = tools.Statistics(lambda ind: ind.fitness.values)
 94    stats.register("avg", numpy.mean)
 95    stats.register("std", numpy.std)
 96    stats.register("min", numpy.min)
 97    stats.register("max", numpy.max)
 98
 99    logbook = tools.Logbook()
100    logbook.header = "gen", "evals", "std", "min", "avg", "max"
101
102    return toolbox, stats, logbook, prim_sets
103
104
105def print_results(best_ind):
106    if not best_ind.fitness.values < (0.5,):
107        raise RuntimeError('Evolution failed to converge.')
108    print(f'\nEvolution converged correctly.')
109
110
111def main():
112    toolbox, stats, logbook, psets = setup()
113    hof = tools.HallOfFame(1)
114
115    def log_stats(ngen=0):
116        hof.update(pop)
117        record = stats.compile(pop)
118        logbook.record(gen=ngen, evals=len(pop), **record)
119        print(logbook.stream)
120
121    pop = toolbox.population(size=100)
122    for ind in pop:
123        ind.fitness.values = toolbox.evaluate(ind)
124
125    log_stats()
126
127    for generations in range(1, GENS):
128        offspring = toolbox.select(pop, len(pop))
129        offspring = [toolbox.clone(ind) for ind in offspring]
130
131        for ind1, ind2 in zip(offspring[::2], offspring[1::2]):
132            for tree1, tree2 in zip(ind1, ind2):
133                if random.random() < CX_PROB:
134                    toolbox.mate(tree1, tree2)
135                    del ind1.fitness.values
136                    del ind2.fitness.values
137
138        for ind in offspring:
139            for tree, pset in zip(ind, psets):
140                if random.random() < MUT_PROB:
141                    toolbox.mutate(individual=tree, prim_set=pset)
142                    del ind.fitness.values
143
144        invalids = [ind for ind in offspring if not ind.fitness.is_valid()]
145        for ind in invalids:
146            ind.fitness.values = toolbox.evaluate(ind)
147        pop = offspring
148
149        log_stats(generations)
150
151    print_results(hof[0])
152
153
154if __name__ == "__main__":
155    main()

Using HARM-GP

 1from deap_er import creator
 2from deap_er import tools
 3from deap_er import base
 4from deap_er import gp
 5import operator
 6import random
 7import numpy
 8import math
 9
10
11random.seed(1234)  # disables randomization
12
13
14def safe_div(left, right):
15    try:
16        return left / right
17    except ZeroDivisionError:
18        return 1
19
20
21def evaluate(individual, points, toolbox):
22    func = toolbox.compile(expr=individual)
23    sq_errors = ((func(x) - x**4 - x**3 - x**2 - x)**2 for x in points)
24    result = math.fsum(sq_errors) / len(points)
25    return result,  # The comma is essential here.
26
27
28def setup():
29    pset = gp.PrimitiveSet("MAIN", 1)
30    pset.add_primitive(operator.add, 2)
31    pset.add_primitive(operator.sub, 2)
32    pset.add_primitive(operator.mul, 2)
33    pset.add_primitive(safe_div, 2)
34    pset.add_primitive(operator.neg, 1)
35    pset.add_primitive(math.cos, 1)
36    pset.add_primitive(math.sin, 1)
37    pset.add_ephemeral_constant("rand101", lambda: random.randint(-1, 1))
38    pset.rename_arguments(ARG0='x')
39
40    creator.create("FitnessMin", base.Fitness, weights=(-1.0,))
41    creator.create("Individual", gp.PrimitiveTree, fitness=creator.FitnessMin)
42
43    toolbox = base.Toolbox()
44    toolbox.register("expr", gp.gen_half_and_half, prim_set=pset, min_depth=1, max_depth=2)
45    toolbox.register("individual", tools.init_iterate, creator.Individual, toolbox.expr)
46    toolbox.register("population", tools.init_repeat, list, toolbox.individual)
47    toolbox.register("compile", gp.compile_tree, prim_set=pset)
48    toolbox.register("mate", gp.cx_one_point)
49    toolbox.register("expr_mut", gp.gen_full, min_depth=0, max_depth=2)
50    toolbox.register("mutate", gp.mut_uniform, expr=toolbox.expr_mut, prim_set=pset)
51    toolbox.register("select", tools.sel_tournament, contestants=3)
52    toolbox.register("evaluate", evaluate, points=[x / 10. for x in range(-10, 10)], toolbox=toolbox)
53    toolbox.decorate("mate", gp.static_limit(limiter=operator.attrgetter("height"), max_value=17))
54    toolbox.decorate("mutate", gp.static_limit(limiter=operator.attrgetter("height"), max_value=17))
55
56    stats_size = tools.Statistics(len)
57    stats_fit = tools.Statistics(lambda ind: ind.fitness.values)
58    mstats = tools.MultiStatistics(fitness=stats_fit, size=stats_size)
59    mstats.register("avg", numpy.mean)
60    mstats.register("std", numpy.std)
61    mstats.register("min", numpy.min)
62    mstats.register("max", numpy.max)
63
64    return toolbox, mstats
65
66
67def print_results(best_ind):
68    if not best_ind.fitness.values < (1.0e-3,):
69        raise RuntimeError('Evolution failed to converge.')
70    print(f'\nEvolution converged correctly.')
71
72
73def main():
74    toolbox, mstats = setup()
75    pop = toolbox.population(size=300)
76    hof = tools.HallOfFame(1)
77    args = dict(
78        toolbox=toolbox,
79        population=pop,
80        generations=40,
81        cx_prob=0.5,
82        mut_prob=0.1,
83        hof=hof,
84        stats=mstats,
85        verbose=True  # prints stats
86    )
87    gp.harm(**args)
88    print_results(hof[0])
89
90
91if __name__ == "__main__":
92    main()

Using Epsilon Lexicase

 1from deap_er import creator
 2from deap_er import tools
 3from deap_er import base
 4from deap_er import gp
 5import operator
 6import random
 7import numpy
 8import math
 9
10
11random.seed(1234)  # disables randomization
12
13
14def safe_div(left, right):
15    try:
16        return left / right
17    except ZeroDivisionError:
18        return 1
19
20
21def evaluate(individual, toolbox, points):
22    func = toolbox.compile(expr=individual)
23    sq_errors = ((func(x) - x**4 - x**3 - x**2 - x)**2 for x in points)
24    result = math.fsum(sq_errors) / len(points)
25    return result,  # The comma is essential here.
26
27
28def setup():
29    pset = gp.PrimitiveSet("MAIN", 1)
30    pset.add_primitive(operator.add, 2)
31    pset.add_primitive(operator.sub, 2)
32    pset.add_primitive(operator.mul, 2)
33    pset.add_primitive(safe_div, 2)
34    pset.add_primitive(operator.neg, 1)
35    pset.add_primitive(math.cos, 1)
36    pset.add_primitive(math.sin, 1)
37    pset.add_ephemeral_constant("rand101", lambda: random.randint(-1, 1))
38    pset.rename_arguments(ARG0='x')
39
40    creator.create("FitnessMin", base.Fitness, weights=(-1.0,))
41    creator.create("Individual", gp.PrimitiveTree, fitness=creator.FitnessMin)
42
43    toolbox = base.Toolbox()
44    toolbox.register("expr", gp.gen_half_and_half, prim_set=pset, min_depth=1, max_depth=2)
45    toolbox.register("individual", tools.init_iterate, creator.Individual, toolbox.expr)
46    toolbox.register("population", tools.init_repeat, list, toolbox.individual)
47    toolbox.register("compile", gp.compile_tree, prim_set=pset)
48    toolbox.register("mate", gp.cx_one_point)
49    toolbox.register("expr_mut", gp.gen_full, min_depth=0, max_depth=2)
50    toolbox.register("mutate", gp.mut_uniform, expr=toolbox.expr_mut, prim_set=pset)
51    toolbox.register("select", tools.sel_epsilon_lexicase)
52    toolbox.register("evaluate", evaluate, toolbox=toolbox, points=[x / 10. for x in range(-10, 10)])
53    toolbox.decorate("mate", gp.static_limit(limiter=operator.attrgetter("height"), max_value=17))
54    toolbox.decorate("mutate", gp.static_limit(limiter=operator.attrgetter("height"), max_value=17))
55
56    stats_size = tools.Statistics(len)
57    stats_fit = tools.Statistics(lambda ind: ind.fitness.values)
58    mstats = tools.MultiStatistics(fitness=stats_fit, size=stats_size)
59    mstats.register("avg", numpy.mean)
60    mstats.register("std", numpy.std)
61    mstats.register("min", numpy.min)
62    mstats.register("max", numpy.max)
63
64    return toolbox, mstats
65
66
67def print_results(best_ind):
68    if not best_ind.fitness.values < (1.0e-3,):
69        raise RuntimeError('Evolution failed to converge.')
70    print(f'\nEvolution converged correctly.')
71
72
73def main():
74    toolbox, mstats = setup()
75    pop = toolbox.population(size=300)
76    hof = tools.HallOfFame(1)
77    args = dict(
78        toolbox=toolbox,
79        population=pop,
80        generations=40,
81        cx_prob=0.5,
82        mut_prob=0.1,
83        hof=hof,
84        stats=mstats,
85        verbose=True  # prints stats
86    )
87    tools.ea_simple(**args)
88    print_results(hof[0])
89
90
91if __name__ == "__main__":
92    main()