rlevo_evolution/algorithms/
gp_cgp.rs

1//! Cartesian Genetic Programming.
2//!
3//! CGP encodes a directed acyclic computation graph on a fixed
4//! `rows × cols` grid. Each node stores `(function_id, input_0, input_1)`,
5//! plus the final output gene picks which node produces the output.
6//! The genotype is a fixed-length integer vector, so populations are
7//! `Tensor<B, 2, Int>` and fit the tensor abstraction cleanly.
8//!
9//! # Evolutionary engine
10//!
11//! Canonical CGP uses a `(1 + λ)` Evolution Strategy with point
12//! mutation and no crossover. This module re-implements just that
13//! engine directly — not via [`crate::algorithms::es_classical`] — so
14//! the mutation logic can be specialized to the CGP genome semantics
15//! (constrained feed-forward connections, `function_id` range, …).
16//!
17//! # Function set
18//!
19//! The v1 function set is fixed at construction time:
20//!
21//! | id | op | arity | formula |
22//! |---|---|---|---|
23//! | 0 | add | 2 | `a + b` |
24//! | 1 | sub | 2 | `a − b` |
25//! | 2 | mul | 2 | `a · b` |
26//! | 3 | `protected_div` | 2 | `a / b` (or `a` if `|b| < ε`) |
27//! | 4 | sin | 1 | `sin(a)` |
28//! | 5 | cos | 1 | `cos(a)` |
29//! | 6 | tanh | 1 | `tanh(a)` |
30//! | 7 | const 1.0 | 0 | `1.0` |
31//!
32//! # Phenotype evaluation
33//!
34//! Evaluation runs on the host because the per-node dispatch is not a
35//! good fit for dense tensor ops; node values are computed in
36//! topological order (left-to-right across the grid columns).
37//! Genotype storage stays on-device to match the other strategies.
38//!
39//! # Reference
40//!
41//! - Miller (2011), *Cartesian Genetic Programming* (Natural Computing
42//!   Series).
43
44use std::marker::PhantomData;
45
46use burn::tensor::{Int, Tensor, TensorData, backend::Backend};
47use rand::{Rng, RngExt};
48
49use crate::rng::{SeedPurpose, seed_stream};
50use crate::strategy::{Strategy, StrategyMetrics};
51
52/// Fixed v1 function set: arity of each opcode.
53pub const FUNCTION_ARITIES: [usize; 8] = [2, 2, 2, 2, 1, 1, 1, 0];
54/// Number of opcodes in the v1 function set.
55pub const NUM_FUNCTIONS: usize = FUNCTION_ARITIES.len();
56
57/// Static configuration for a [`CartesianGeneticProgramming`] run.
58#[derive(Debug, Clone)]
59pub struct CgpConfig {
60    /// Number of offspring per generation (λ in `(1 + λ)`).
61    pub lambda: usize,
62    /// Number of inputs (independent variables) the program sees.
63    pub n_inputs: usize,
64    /// Number of grid rows.
65    pub rows: usize,
66    /// Number of grid columns.
67    pub cols: usize,
68    /// Mutation rate applied to each gene of the integer genome.
69    pub mutation_rate: f32,
70    /// Levels-back parameter: how many previous columns a node can
71    /// connect to. `usize::MAX` means "any previous column".
72    pub levels_back: usize,
73}
74
75impl CgpConfig {
76    /// Sensible defaults: 1-output, 1-row, 30-column grid, mutation
77    /// rate tuned to flip ~3 genes per genome.
78    #[must_use]
79    pub fn default_for(n_inputs: usize) -> Self {
80        let rows = 1;
81        let cols = 30;
82        let genes_per_node = 3; // (function, input_0, input_1)
83        let output_genes = 1;
84        let total_genes = rows * cols * genes_per_node + output_genes;
85        #[allow(clippy::cast_precision_loss)]
86        let mutation_rate = 3.0 / total_genes as f32;
87        Self {
88            lambda: 4,
89            n_inputs,
90            rows,
91            cols,
92            mutation_rate,
93            levels_back: usize::MAX,
94        }
95    }
96
97    /// Genes per node in the genotype layout: `(function_id, input_0, input_1)`.
98    pub const GENES_PER_NODE: usize = 3;
99    /// Number of output genes per program (one index pointing to the node
100    /// whose value is taken as the program output).
101    pub const OUTPUT_GENES: usize = 1;
102
103    /// Total genome length (nodes × 3 + outputs).
104    #[must_use]
105    pub fn genome_len(&self) -> usize {
106        self.rows * self.cols * Self::GENES_PER_NODE + Self::OUTPUT_GENES
107    }
108}
109
110/// Generation state for [`CartesianGeneticProgramming`].
111#[derive(Debug, Clone)]
112pub struct CgpState<B: Backend> {
113    /// Parent genotype, shape `(1, genome_len)`.
114    pub parent: Tensor<B, 2, Int>,
115    /// Parent fitness (host-side scalar cache).
116    pub parent_fitness: f32,
117    /// Best-so-far genotype.
118    pub best_genome: Option<Tensor<B, 2, Int>>,
119    /// Best-so-far fitness.
120    pub best_fitness: f32,
121    /// Generation counter.
122    pub generation: usize,
123}
124
125/// Classical Cartesian GP with `(1 + λ)` ES.
126///
127/// # Example
128///
129/// ```no_run
130/// use burn::backend::Flex;
131/// use rlevo_evolution::algorithms::gp_cgp::{CartesianGeneticProgramming, CgpConfig};
132///
133/// let strategy = CartesianGeneticProgramming::<Flex>::new();
134/// let params = CgpConfig::default_for(1);
135/// assert!(params.genome_len() > 0);
136/// let _ = strategy;
137/// ```
138#[derive(Debug, Clone, Copy, Default)]
139pub struct CartesianGeneticProgramming<B: Backend> {
140    _backend: PhantomData<fn() -> B>,
141}
142
143impl<B: Backend> CartesianGeneticProgramming<B> {
144    /// Builds a new (stateless) strategy object.
145    #[must_use]
146    pub fn new() -> Self {
147        Self {
148            _backend: PhantomData,
149        }
150    }
151
152    fn sample_initial_genome(params: &CgpConfig, rng: &mut dyn Rng) -> Vec<i64> {
153        let mut genome = Vec::with_capacity(params.genome_len());
154        for col in 0..params.cols {
155            for _row in 0..params.rows {
156                #[allow(clippy::cast_possible_wrap)]
157                let func = rng.random_range(0..NUM_FUNCTIONS as i64);
158                let (inp0, inp1) = sample_input_pair(col, params, rng);
159                genome.push(func);
160                genome.push(inp0);
161                genome.push(inp1);
162            }
163        }
164        // Output gene: any node index or input index.
165        let max_node_idx = params.n_inputs + params.rows * params.cols;
166        #[allow(clippy::cast_possible_wrap)]
167        genome.push(rng.random_range(0..max_node_idx as i64));
168        genome
169    }
170
171    fn genome_to_host(genome: &Tensor<B, 2, Int>) -> Vec<i64> {
172        genome
173            .clone()
174            .into_data()
175            .into_vec::<i32>()
176            .unwrap_or_default()
177            .into_iter()
178            .map(i64::from)
179            .collect()
180    }
181}
182
183fn sample_input_pair(col: usize, params: &CgpConfig, rng: &mut dyn Rng) -> (i64, i64) {
184    let min_col = col.saturating_sub(params.levels_back);
185    let node_indices_start = params.n_inputs + min_col * params.rows;
186    let node_indices_end = params.n_inputs + col * params.rows;
187    let max = node_indices_end.max(params.n_inputs);
188    // Allowed inputs: 0..n_inputs (graph inputs) ∪ previous nodes.
189    let input_count = params.n_inputs
190        + (max - params.n_inputs)
191            .saturating_sub(node_indices_start.saturating_sub(params.n_inputs));
192    let pool: Vec<i64> = (0..params.n_inputs)
193        .chain(node_indices_start..node_indices_end)
194        .map(|i| {
195            #[allow(clippy::cast_possible_wrap)]
196            let v = i as i64;
197            v
198        })
199        .collect();
200    let pool = if pool.is_empty() {
201        #[allow(clippy::cast_possible_wrap)]
202        (0..params.n_inputs as i64).collect()
203    } else {
204        pool
205    };
206    let _ = input_count;
207    let pick = |rng: &mut dyn Rng| -> i64 {
208        let idx = rng.random_range(0..pool.len());
209        pool[idx]
210    };
211    (pick(rng), pick(rng))
212}
213
214fn mutate_genome(genome: &mut [i64], params: &CgpConfig, rng: &mut dyn Rng) {
215    let genes_per_node = CgpConfig::GENES_PER_NODE;
216    let node_genes = params.rows * params.cols * genes_per_node;
217    for (gene_idx, gene) in genome.iter_mut().enumerate() {
218        if rng.random::<f32>() >= params.mutation_rate {
219            continue;
220        }
221        if gene_idx < node_genes {
222            let within = gene_idx % genes_per_node;
223            let node_idx = gene_idx / genes_per_node;
224            let col = node_idx / params.rows;
225            if within == 0 {
226                // function
227                #[allow(clippy::cast_possible_wrap)]
228                {
229                    *gene = rng.random_range(0..NUM_FUNCTIONS as i64);
230                }
231            } else {
232                let (new0, new1) = sample_input_pair(col, params, rng);
233                *gene = if within == 1 { new0 } else { new1 };
234            }
235        } else {
236            // output gene
237            let max_node_idx = params.n_inputs + params.rows * params.cols;
238            #[allow(clippy::cast_possible_wrap)]
239            {
240                *gene = rng.random_range(0..max_node_idx as i64);
241            }
242        }
243    }
244}
245
246/// Evaluates a CGP genotype at a set of input rows.
247///
248/// `genome` is the host-side integer genotype (length `params.genome_len()`).
249/// `inputs` is a slice of `n_samples` rows, each of length `params.n_inputs`.
250/// Returns one `f32` output per input row.
251///
252/// Out-of-range input/node indices in the genome are clamped to the
253/// last buffer slot rather than panicking — this keeps fitness
254/// evaluation robust to mutated-but-unrepaired genotypes. Non-finite
255/// node values (e.g., `inf` from divisions or `tan`) collapse to `0.0`.
256///
257/// # Panics
258///
259/// Panics if `genome` is empty (the last gene is the output index).
260#[must_use]
261pub fn evaluate_cgp(genome: &[i64], params: &CgpConfig, inputs: &[Vec<f32>]) -> Vec<f32> {
262    let node_count = params.rows * params.cols;
263    let n_inputs = params.n_inputs;
264    #[allow(clippy::cast_sign_loss, clippy::cast_possible_truncation)]
265    let output_idx = genome[genome.len() - 1] as usize;
266
267    let mut outputs = Vec::with_capacity(inputs.len());
268    let mut buf = vec![0.0_f32; n_inputs + node_count];
269
270    for sample in inputs {
271        for (i, v) in sample.iter().enumerate() {
272            buf[i] = *v;
273        }
274        for node in 0..node_count {
275            let base = node * 3;
276            #[allow(clippy::cast_sign_loss, clippy::cast_possible_truncation)]
277            let func = genome[base] as usize;
278            #[allow(clippy::cast_sign_loss, clippy::cast_possible_truncation)]
279            let a_idx = genome[base + 1] as usize;
280            #[allow(clippy::cast_sign_loss, clippy::cast_possible_truncation)]
281            let b_idx = genome[base + 2] as usize;
282            let a = buf[a_idx.min(buf.len() - 1)];
283            let b = buf[b_idx.min(buf.len() - 1)];
284            let v = match func {
285                0 => a + b,
286                1 => a - b,
287                2 => a * b,
288                3 => {
289                    if b.abs() < 1e-6 {
290                        a
291                    } else {
292                        a / b
293                    }
294                }
295                4 => a.sin(),
296                5 => a.cos(),
297                6 => a.tanh(),
298                7 => 1.0,
299                _ => 0.0,
300            };
301            buf[n_inputs + node] = if v.is_finite() { v } else { 0.0 };
302        }
303        outputs.push(buf[output_idx.min(buf.len() - 1)]);
304    }
305
306    outputs
307}
308
309impl<B: Backend> Strategy<B> for CartesianGeneticProgramming<B>
310where
311    B::Device: Clone,
312{
313    type Params = CgpConfig;
314    type State = CgpState<B>;
315    type Genome = Tensor<B, 2, Int>;
316
317    /// Samples the initial parent genome by drawing random node functions and
318    /// feed-forward input connections via `rng`, then uploads the genotype as
319    /// a `(1, genome_len)` integer tensor.
320    fn init(&self, params: &CgpConfig, rng: &mut dyn Rng, device: &<B as burn::tensor::backend::BackendTypes>::Device) -> CgpState<B> {
321        let genome_vec = Self::sample_initial_genome(params, rng);
322        let parent = Tensor::<B, 2, Int>::from_data(
323            TensorData::new(genome_vec, [1, params.genome_len()]),
324            device,
325        );
326        CgpState {
327            parent,
328            parent_fitness: f32::INFINITY,
329            best_genome: None,
330            best_fitness: f32::INFINITY,
331            generation: 0,
332        }
333    }
334
335    /// Returns the offspring population for the current generation.
336    ///
337    /// On the first call (parent fitness not yet set), returns the single
338    /// parent genome unchanged for initial fitness evaluation.
339    /// On subsequent calls, produces `params.lambda` children by cloning the
340    /// parent and applying per-gene point mutation, with mutation draws taken
341    /// from a deterministic `seed_stream` (host-RNG convention).
342    fn ask(
343        &self,
344        params: &CgpConfig,
345        state: &CgpState<B>,
346        rng: &mut dyn Rng,
347        device: &<B as burn::tensor::backend::BackendTypes>::Device,
348    ) -> (Tensor<B, 2, Int>, CgpState<B>) {
349        // First call: evaluate the parent as "offspring" of size 1.
350        if !state.parent_fitness.is_finite() {
351            return (state.parent.clone(), state.clone());
352        }
353
354        let mut mut_rng = seed_stream(
355            rng.next_u64(),
356            state.generation as u64,
357            SeedPurpose::Mutation,
358        );
359        let parent_vec = Self::genome_to_host(&state.parent);
360        let mut offspring_genomes: Vec<i64> =
361            Vec::with_capacity(params.lambda * params.genome_len());
362        for _ in 0..params.lambda {
363            let mut child = parent_vec.clone();
364            mutate_genome(&mut child, params, &mut mut_rng);
365            offspring_genomes.extend(child);
366        }
367        #[allow(clippy::cast_possible_truncation)]
368        let offspring_genomes_i32: Vec<i32> =
369            offspring_genomes.into_iter().map(|v| v as i32).collect();
370        let offspring = Tensor::<B, 2, Int>::from_data(
371            TensorData::new(offspring_genomes_i32, [params.lambda, params.genome_len()]),
372            device,
373        );
374        (offspring, state.clone())
375    }
376
377    /// Applies `(1+λ)` selection and returns the updated state.
378    ///
379    /// The canonical CGP tie-breaking rule is used: an offspring replaces the
380    /// parent when its fitness is **less than or equal to** the parent's,
381    /// allowing neutral mutations to accumulate and maintain genetic diversity
382    /// in the inactive (non-coding) portion of the genome.
383    ///
384    /// The first `tell` after `init` bootstraps the parent fitness from the
385    /// initial single-genome evaluation rather than running selection.
386    fn tell(
387        &self,
388        _params: &CgpConfig,
389        offspring: Tensor<B, 2, Int>,
390        fitness: Tensor<B, 1>,
391        mut state: CgpState<B>,
392        _rng: &mut dyn Rng,
393    ) -> (CgpState<B>, StrategyMetrics) {
394        let fitness_host = fitness.into_data().into_vec::<f32>().unwrap_or_default();
395
396        if !state.parent_fitness.is_finite() {
397            // First tell: initial parent fitness.
398            state.parent_fitness = fitness_host[0];
399            state.generation += 1;
400            update_best(&mut state, &offspring, &fitness_host);
401            let m = StrategyMetrics::from_host_fitness(
402                state.generation,
403                &fitness_host,
404                state.best_fitness,
405            );
406            state.best_fitness = m.best_fitness_ever;
407            return (state, m);
408        }
409
410        // (1+λ): parent survives only if NO offspring strictly beats it;
411        // canonical CGP uses `<=` to break ties in favor of offspring
412        // (neutral mutations accumulate).
413        let best_off_idx = fitness_host
414            .iter()
415            .enumerate()
416            .min_by(|(_, a), (_, b)| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal))
417            .map_or(0, |(i, _)| i);
418        let best_off_fit = fitness_host[best_off_idx];
419        if best_off_fit <= state.parent_fitness {
420            let device = offspring.device();
421            #[allow(clippy::cast_possible_wrap, clippy::cast_possible_truncation)]
422            let idx = Tensor::<B, 1, Int>::from_data(
423                TensorData::new(vec![best_off_idx as i32], [1]),
424                &device,
425            );
426            state.parent = offspring.clone().select(0, idx);
427            state.parent_fitness = best_off_fit;
428        }
429
430        state.generation += 1;
431        update_best(&mut state, &offspring, &fitness_host);
432        let m =
433            StrategyMetrics::from_host_fitness(state.generation, &fitness_host, state.best_fitness);
434        state.best_fitness = m.best_fitness_ever;
435        (state, m)
436    }
437
438    /// Returns the best-so-far genome and its fitness, or `None` before the
439    /// first `tell` call.
440    fn best(&self, state: &CgpState<B>) -> Option<(Tensor<B, 2, Int>, f32)> {
441        state
442            .best_genome
443            .as_ref()
444            .map(|g| (g.clone(), state.best_fitness))
445    }
446}
447
448fn update_best<B: Backend>(state: &mut CgpState<B>, pop: &Tensor<B, 2, Int>, fitness: &[f32]) {
449    if fitness.is_empty() {
450        return;
451    }
452    let mut best_idx = 0usize;
453    let mut best_f = fitness[0];
454    for (i, &f) in fitness.iter().enumerate().skip(1) {
455        if f < best_f {
456            best_f = f;
457            best_idx = i;
458        }
459    }
460    if best_f < state.best_fitness {
461        let device = pop.device();
462        #[allow(clippy::cast_possible_wrap, clippy::cast_possible_truncation)]
463        let idx =
464            Tensor::<B, 1, Int>::from_data(TensorData::new(vec![best_idx as i32], [1]), &device);
465        state.best_genome = Some(pop.clone().select(0, idx));
466        state.best_fitness = best_f;
467    }
468}
469
470#[cfg(test)]
471mod tests {
472    use super::*;
473    use crate::fitness::BatchFitnessFn;
474    use crate::strategy::EvolutionaryHarness;
475    use burn::backend::Flex;
476    type TestBackend = Flex;
477
478    /// Symbolic regression on `x² + 1` over 20 evenly spaced x ∈ [−1, 1].
479    struct SymRegression {
480        params: CgpConfig,
481        xs: Vec<f32>,
482        ys: Vec<f32>,
483    }
484
485    impl SymRegression {
486        #[allow(clippy::cast_precision_loss)]
487        fn new(params: CgpConfig) -> Self {
488            let xs: Vec<f32> = (0..20).map(|i| -1.0 + 2.0 * (i as f32) / 19.0).collect();
489            let ys: Vec<f32> = xs.iter().map(|x| x * x + 1.0).collect();
490            Self { params, xs, ys }
491        }
492    }
493
494    impl<B: Backend> BatchFitnessFn<B, Tensor<B, 2, Int>> for SymRegression {
495        #[allow(clippy::cast_precision_loss)]
496        fn evaluate_batch(
497            &mut self,
498            population: &Tensor<B, 2, Int>,
499            device: &<B as burn::tensor::backend::BackendTypes>::Device,
500        ) -> Tensor<B, 1> {
501            let pop_size = population.dims()[0];
502            let data: Vec<i64> = population
503                .clone()
504                .into_data()
505                .into_vec::<i32>()
506                .unwrap()
507                .into_iter()
508                .map(i64::from)
509                .collect();
510            let gl = self.params.genome_len();
511            let inputs: Vec<Vec<f32>> = self.xs.iter().map(|&x| vec![x]).collect();
512            let mut fitness = Vec::with_capacity(pop_size);
513            for row in 0..pop_size {
514                let genome = &data[row * gl..(row + 1) * gl];
515                let preds = evaluate_cgp(genome, &self.params, &inputs);
516                let mse: f32 = preds
517                    .iter()
518                    .zip(self.ys.iter())
519                    .map(|(p, y)| (p - y).powi(2))
520                    .sum::<f32>()
521                    / (self.ys.len() as f32);
522                fitness.push(mse);
523            }
524            Tensor::<B, 1>::from_data(TensorData::new(fitness, [pop_size]), device)
525        }
526    }
527
528    #[test]
529    #[allow(clippy::cast_precision_loss)]
530    fn cgp_reduces_error_on_square_plus_one() {
531        let device = Default::default();
532        let params = CgpConfig::default_for(1);
533        let landscape = SymRegression::new(params.clone());
534        let initial_error = {
535            // Baseline: random genome MSE on a single seed.
536            use rand::SeedableRng;
537            let mut rng = rand::rngs::StdRng::seed_from_u64(123);
538            let genome = CartesianGeneticProgramming::<TestBackend>::sample_initial_genome(
539                &params, &mut rng,
540            );
541            let inputs: Vec<Vec<f32>> = landscape.xs.iter().map(|&x| vec![x]).collect();
542            let preds = evaluate_cgp(&genome, &params, &inputs);
543            preds
544                .iter()
545                .zip(landscape.ys.iter())
546                .map(|(p, y)| (p - y).powi(2))
547                .sum::<f32>()
548                / (landscape.ys.len() as f32)
549        };
550
551        let mut harness = EvolutionaryHarness::<TestBackend, _, _>::new(
552            CartesianGeneticProgramming::<TestBackend>::new(),
553            params,
554            landscape,
555            21,
556            device,
557            2000,
558        );
559        harness.reset();
560        loop {
561            if harness.step(()).done {
562                break;
563            }
564        }
565        let best = harness.latest_metrics().unwrap().best_fitness_ever;
566        // CGP should substantially beat the random-genome baseline.
567        assert!(
568            best < initial_error,
569            "CGP did not improve: best={best} initial={initial_error}"
570        );
571        // Bias check: ought to beat predicting a constant y=1 (mean ~= 1.33).
572        assert!(best < 0.2, "expected MSE < 0.2 but got {best}");
573    }
574}
rlevo_evolution/algorithms/gp_cgp.rs

rlevo_evolution/algorithms/
gp_cgp.rs
