|x| by strongly typed GP
The problem
Symbolic regression again, as in Koza's quartic: given points (x, y), find a formula for y. Here y = |x|, at 20 values of x drawn at random from [−1, 1], and the formula is found when it's recovered exactly: an error at the level of rounding on the 20 training points and on 101 test points spread evenly over [−1, 1].
The absolute value has a kink at 0, and no polynomial has one. It takes a decision: −x where x is
negative, else x. So the formulas are built from two types of values, real numbers and Booleans:
less compares two reals and returns a Boolean, if takes a Boolean and two reals, and and,
or and not combine Booleans. This is strongly typed genetic programming (Montana 1995): every
function says what types its arguments and its result have, and the search only makes programs in
which they fit.
The Python version declares the same typed set with gx.gp.PrimitiveSetBuilder and evaluates
each tree with numpy, tree.evaluate({"x": xs}, functions): one call per node on the columns of
all the points (np.add, np.less, np.where for if), whose IEEE arithmetic gives the bits that
Rust computes point by point. It prints the same output and writes the same trace.
What makes it hard
Without types, a function set has to make every function accept every value: Koza (1992) avoided predicates such as a comparison that returns a Boolean, as Montana notes. Types let the set say what a program means: a comparison goes where a condition goes, and a number never does. With types, generation, crossover and mutation have to keep every program well typed: a subtree is only ever replaced by one of its own type.
The problem itself is small. Close fits without a decision miss by far more than rounding (x² has an RMSE of 0.18 over [−1, 1]), and recovering |x| exactly needs the comparison at 0 and the two branches right.
Representation
A tree of a genetic program (gp::Tree) over a gp::PrimitiveSet of two types, real and
bool: add, sub and mul of two reals; less(a, b), true when a < b; and, or and not of Booleans;
if(c, a, b), a when c is true, else b; the variable x; and ephemeral random constants, the integers
from −2 to 2, which a real can be. Trees return a real. gp::Gp sets Koza's limits and
initialization: depth at most 17 and at most 1024 nodes, ramped half-and-half with depths 2 to 6,
the initial population divided evenly without duplicates (Gp::ramped_half_and_half).
The fitness is the root mean squared error on the 20 points, minimized. Each tree is evaluated at
each point with Tree::evaluate, on a stack of values of an enum with a real and a Boolean
variant: the set's types guarantee which one each function gets.
Algorithm
A genetic algorithm of 1000 trees, with double tournaments (fitness tournaments of 7, the smaller of their two winners chosen with probability 0.7) against bloat, subtree crossover at a rate of 0.9, and at a rate of 0.1 one of four mutations: subtree (weight 0.5), point (0.3), hoist (0.1) and shrink (0.1). Every operator keeps the trees typed. The run stops at an RMSE of 10⁻¹⁰ times the standard deviation of the 20 values, or after 100 generations.
Output
The first two lines give the problem and the setting. Then come how the run stopped, after how many generations and evaluations, the RMSE of the best tree on the training and the test points, and the tree itself, with its size and depth.
The project page plays this run back.
Good results
The optimum is |x| itself, recovered exactly. The run of output.txt found it after 2 generations
and 2600 evaluations, as
if(less(x, 0.0), mul(-1.0, x), sub(x, 0.0)),
which is −x where x < 0, else x: |x|, with an error of 0 on the training and the test points.
Over seeds 1 to 100 (the seed of the algorithm and of the initial population, with the same data), all 100 runs recovered |x|: after 3 generations and 3285 evaluations in the median, and 12 generations at most. The expressions found had 20 nodes in the median and 110 at most; some reach the same function by longer routes, such as 2 · if(x < 0, −x, 0) + x.
Reference: Montana, D. J. (1995). Strongly typed genetic programming. Evolutionary Computation 3(2): 199-230.
Known optimum: |x| exactly (an RMSE of 0, up to rounding)
Source: examples/abs_typed
Interactive run: tachsin.gr/projects/genoxide/examples/abs-typed
cargo run --release --example abs_typed
//! |x| by strongly typed genetic programming: find the absolute value from 20 points of it, with
//! a comparison that returns a Boolean and a conditional that takes one.
//!
//! Two types, real numbers and Booleans (Montana 1995): `less` compares two reals and returns a
//! Boolean, `if` takes a Boolean and two reals, and every tree genoxide makes puts a Boolean
//! where a Boolean goes. The run stops at exact recovery: an error at the level of rounding, on
//! the 20 training points and on 101 test points across [−1, 1], with the expression printed.
//!
//! With `GENOXIDE_TRACE=<file>`, it also writes a trace of its run for the plot on the example's
//! page, with `trace.rs`.
//!
//! ```text
//! cargo run --release --example abs_typed
//! ```
mod trace;
use genoxide::gp::{Constants, Gp, Mutations, PrimitiveSet, SubtreeCrossover, Tree};
use genoxide::prelude::*;
use rand::RngExt;
// the primitives: arithmetic on reals, a comparison, Boolean functions and a conditional
#[derive(Clone, Copy, Debug, PartialEq)]
pub enum Op {
Add,
Sub,
Mul,
Less,
And,
Or,
Not,
If,
X,
}
pub fn primitives() -> Result<PrimitiveSet<Op>> {
let mut set = PrimitiveSet::builder();
let real = set.new_type("real");
let boolean = set.new_type("bool");
set.function("add", Op::Add, [real, real], real)
.function("sub", Op::Sub, [real, real], real)
.function("mul", Op::Mul, [real, real], real)
.function("less", Op::Less, [real, real], boolean)
.function("and", Op::And, [boolean, boolean], boolean)
.function("or", Op::Or, [boolean, boolean], boolean)
.function("not", Op::Not, [boolean], boolean)
.function("if", Op::If, [boolean, real, real], real)
.terminal("x", Op::X, real)
// ephemeral random constants: the integers -2 to 2
.constants(real, Constants::integers(-2..=2)?);
set.build(real)
}
// a value of either type; the set's types guarantee which one each primitive gets
#[derive(Clone, Copy, Debug)]
enum Value {
Real(f64),
Bool(bool),
}
impl Value {
fn real(self) -> f64 {
match self {
Value::Real(value) => value,
Value::Bool(_) => unreachable!("a real, by the set's types"),
}
}
fn bool(self) -> bool {
match self {
Value::Bool(value) => value,
Value::Real(_) => unreachable!("a Boolean, by the set's types"),
}
}
}
// the value of `op` from its arguments' values, at x
fn apply(op: Op, args: &[Value], x: f64) -> Value {
let real = |i: usize| args[i].real();
let bool = |i: usize| args[i].bool();
match op {
Op::Add => Value::Real(real(0) + real(1)),
Op::Sub => Value::Real(real(0) - real(1)),
Op::Mul => Value::Real(real(0) * real(1)),
Op::Less => Value::Bool(real(0) < real(1)),
Op::And => Value::Bool(bool(0) && bool(1)),
Op::Or => Value::Bool(bool(0) || bool(1)),
Op::Not => Value::Bool(!bool(0)),
Op::If => args[if bool(0) { 1 } else { 2 }],
Op::X => Value::Real(x),
}
}
// points and |x| at them
pub struct Data {
pub xs: Vec<f64>,
pub ys: Vec<f64>,
}
impl Data {
fn new(xs: Vec<f64>) -> Self {
let ys = xs.iter().map(|x| x.abs()).collect();
Self { xs, ys }
}
// the tree's values at the points
pub fn predict(&self, set: &PrimitiveSet<Op>, tree: &Tree) -> Vec<f64> {
let mut stack = Vec::new();
self.xs
.iter()
.map(|&x| {
let value = tree.evaluate(
set,
&mut stack,
|op, args| apply(op, args, x),
|_, constant| Value::Real(constant), // only reals have constants
);
value.real()
})
.collect()
}
// the root mean squared error of the tree
pub fn rmse(&self, set: &PrimitiveSet<Op>, tree: &Tree) -> f64 {
let predictions = self.predict(set, tree);
let mut sum = 0.0;
for (prediction, y) in predictions.iter().zip(&self.ys) {
let error = prediction - y;
sum += error * error;
}
(sum / self.ys.len() as f64).sqrt()
}
// the standard deviation of the values
fn deviation(&self) -> f64 {
let n = self.ys.len() as f64;
let mean = self.ys.iter().sum::<f64>() / n;
let squares: f64 = self.ys.iter().map(|y| (y - mean) * (y - mean)).sum();
(squares / n).sqrt()
}
}
const POPULATION: usize = 1000;
fn main() -> Result<()> {
// 20 training points uniform in [-1, 1], from a fixed seed, and 101 test points evenly spaced
let mut rng = StreamRng::seed_from_u64(20);
let training = Data::new((0..20).map(|_| rng.random_range(-1.0..=1.0)).collect());
let test = Data::new((0..=100).map(|i| f64::from(i) / 50.0 - 1.0).collect());
// exact recovery: an error of at most 1e-10 of the values' standard deviation
let tolerance = 1e-10 * training.deviation();
println!("|x| from 20 points in [-1, 1], with real and Boolean types");
println!("{POPULATION} trees, until the RMSE is at most {tolerance:.2e}\n");
let gp = Gp::builder(primitives()?).build()?;
let set = gp.primitives().clone();
let initial = gp.ramped_half_and_half(POPULATION, &mut StreamRng::seed_from_u64(1))?;
let ga = Ga::builder(gp)
.population_size(POPULATION)
.initial_genomes(initial)
.select(DoubleTournament::new(7, 1.4)?)
.crossover(SubtreeCrossover::new())
.mutate(
Mutations::builder()
.subtree(0.5)
.point(0.3)
.hoist(0.1)
.shrink(0.1)
.build()?,
)
.crossover_rate(0.9)
.mutation_rate(0.1)
.minimize()
.seed(1)
.build()?;
let mut trace = trace::Trace::from_env(&training, &set);
let outcome = Engine::new(ga, |tree: &Tree| training.rmse(&set, tree))
.stop_when(Stop::target(tolerance).or(Stop::generations(100)))
.on_generation(|snapshot| trace.record(snapshot))
.run()?;
let best = outcome.best_genome();
println!(
"{:?} after {} generations and {} evaluations",
outcome.stop_reason(),
outcome.generations(),
outcome.evaluations()
);
println!(
"RMSE on the 20 training points: {:.2e}",
training.rmse(&set, best)
);
println!(
"RMSE on the 101 test points: {:.2e}",
test.rmse(&set, best)
);
println!(
"\nthe expression, {} nodes of depth {}:\n{}",
best.len(),
best.depth(&set),
best.display(&set)
);
trace.write();
Ok(())
}
python examples/abs_typed/main.py
"""|x| by strongly typed genetic programming: find the absolute value from 20 points of it, with a
comparison that returns a Boolean and a conditional that takes one.
Two types, real numbers and Booleans (Montana 1995): ``less`` compares two reals and returns a
Boolean, ``if`` takes a Boolean and two reals, and every tree genoxide makes puts a Boolean where
a Boolean goes. The run stops at exact recovery: an error at the level of rounding, on the 20
training points and on 101 test points across [-1, 1], with the expression printed.
The primitives are the program's own (``gx.gp.PrimitiveSetBuilder``), evaluated by numpy: each
function is called once per node, on the columns of all the points. Their arithmetic is IEEE's,
the same bits as Rust's point by point, so the run is the Rust example's, to the bit.
With ``GENOXIDE_TRACE=<file>``, it also writes a trace of its run for the plot on the example's
page, with trace.py.
python examples/abs_typed/main.py
"""
import math
import numpy as np
import genoxide as gx
from trace import Trace
POPULATION = 1000
def primitives():
"""The primitives: arithmetic on reals, a comparison, Boolean functions and a conditional."""
builder = gx.gp.PrimitiveSetBuilder()
real = builder.new_type("real")
boolean = builder.new_type("bool")
builder.function("add", [real, real], real)
builder.function("sub", [real, real], real)
builder.function("mul", [real, real], real)
builder.function("less", [real, real], boolean)
builder.function("and", [boolean, boolean], boolean)
builder.function("or", [boolean, boolean], boolean)
builder.function("not", [boolean], boolean)
builder.function("if", [boolean, real, real], real)
builder.terminal("x", real)
# ephemeral random constants: the integers -2 to 2
builder.constants(real, gx.gp.Constants.integers(-2, 2))
return builder.build(real)
# what the functions mean, on numpy columns: each is called once per node, on all the points
FUNCTIONS = {
"add": np.add,
"sub": np.subtract,
"mul": np.multiply,
"less": np.less,
"and": np.logical_and,
"or": np.logical_or,
"not": np.logical_not,
"if": np.where,
}
class Data:
"""Points and |x| at them."""
def __init__(self, xs):
self.xs = np.asarray(xs, dtype=float)
self.ys = np.abs(self.xs)
def predict(self, tree):
"""The tree's values at the points (a constant broadcast to all of them)."""
values = tree.evaluate({"x": self.xs}, FUNCTIONS)
return np.broadcast_to(np.asarray(values, dtype=float), self.xs.shape)
def rmse(self, tree):
"""The root mean squared error of the tree, summed in order as in Rust."""
total = 0.0
for error in (self.predict(tree) - self.ys).tolist():
total += error * error
return math.sqrt(total / len(self.ys))
def deviation(self):
"""The standard deviation of the values."""
ys = self.ys.tolist()
mean = 0.0
for y in ys:
mean += y
mean /= len(ys)
squares = 0.0
for y in ys:
squares += (y - mean) * (y - mean)
return math.sqrt(squares / len(ys))
def scientific(value):
"""``value`` with 2 decimals and an exponent, as Rust's ``{:.2e}`` writes it: 1.27e-10,
0.00e0."""
if value != value:
return "NaN"
mantissa, exponent = f"{value:.2e}".split("e")
return f"{mantissa}e{int(exponent)}"
# 20 training points uniform in [-1, 1], from a fixed seed (Koza's sampling, the points of his
# quartic), and 101 test points evenly spaced
training = Data(gx.gp.regression.problems.Koza1().dataset().training.x[:, 0])
test = Data([i / 50.0 - 1.0 for i in range(101)])
# exact recovery: an error of at most 1e-10 of the values' standard deviation
tolerance = 1e-10 * training.deviation()
print("|x| from 20 points in [-1, 1], with real and Boolean types")
print(f"{POPULATION} trees, until the RMSE is at most {scientific(tolerance)}")
print()
gp = gx.gp.Gp(primitives())
ga = gx.Ga(
gp,
population_size=POPULATION,
initial_genomes=gp.ramped_half_and_half(POPULATION, 1),
select=gx.DoubleTournament(7, 1.4),
crossover=gx.gp.SubtreeCrossover(),
mutation=gx.gp.Mutations(
[
(0.5, gx.gp.SubtreeMutation()),
(0.3, gx.gp.PointMutation(count=1)),
(0.1, gx.gp.HoistMutation()),
(0.1, gx.gp.ShrinkMutation()),
]
),
crossover_rate=0.9,
mutation_rate=0.1,
objective="minimize",
seed=1,
)
trace = Trace(training)
result = ga.run(
training.rmse, target=tolerance, generations=100, on_generation=trace.on_generation
)
best = result.best_genome
print(
f"{result.stop_reason.capitalize()} after {result.generations} generations and "
f"{result.evaluations} evaluations"
)
print(f"RMSE on the 20 training points: {scientific(training.rmse(best))}")
print(f"RMSE on the 101 test points: {scientific(test.rmse(best))}")
print()
print(f"the expression, {len(best)} nodes of depth {best.depth}:")
print(best)
trace.write()
What it prints, from a seeded run:
|x| from 20 points in [-1, 1], with real and Boolean types
1000 trees, until the RMSE is at most 3.39e-11
Target after 2 generations and 2600 evaluations
RMSE on the 20 training points: 0.00e0
RMSE on the 101 test points: 0.00e0
the expression, 10 nodes of depth 2:
if(less(x, 0.0), mul(-1.0, x), sub(x, 0.0))