rlevo_evolution/algorithms/

gp_cgp.rs

//! Cartesian Genetic Programming.
//!
//! CGP encodes a directed acyclic computation graph on a fixed
//! `rows × cols` grid. Each node stores `(function_id, input_0, input_1)`,
//! plus the final output gene picks which node produces the output.
//! The genotype is a fixed-length integer vector, so populations are
//! `Tensor<B, 2, Int>` and fit the tensor abstraction cleanly.
//!
//! # Evolutionary engine
//!
//! Canonical CGP uses a `(1 + λ)` Evolution Strategy with point
//! mutation and no crossover. This module re-implements just that
//! engine directly — not via [`crate::algorithms::es_classical`] — so
//! the mutation logic can be specialized to the CGP genome semantics
//! (constrained feed-forward connections, `function_id` range, …).
//!
//! # Function set
//!
//! The v1 function set is the shared [`ArithmeticFunctionSet`]:
//!
//! | id | op | arity | formula |
//! |---|---|---|---|
//! | 0 | add | 2 | `a + b` |
//! | 1 | sub | 2 | `a − b` |
//! | 2 | mul | 2 | `a · b` |
//! | 3 | `protected_div` | 2 | `a / b` (or `a` if `|b| < ε`) |
//! | 4 | sin | 1 | `sin(a)` |
//! | 5 | cos | 1 | `cos(a)` |
//! | 6 | tanh | 1 | `tanh(a)` |
//! | 7 | const 1.0 | 0 | `1.0` |
//!
//! Opcode evaluation is delegated to the [`FunctionSet`] trait: [`evaluate_cgp`]
//! is a thin wrapper over the generic [`evaluate_cgp_with`], which threads a
//! concrete `&F` through the per-node loop so the opcode dispatch inlines. The
//! [`FUNCTION_ARITIES`] / [`NUM_FUNCTIONS`] constants are retained for the
//! mutation logic, which samples function ids in `0..NUM_FUNCTIONS`.
//!
//! # Phenotype evaluation
//!
//! Evaluation runs on the host because the per-node dispatch is not a
//! good fit for dense tensor ops; node values are computed in
//! topological order (left-to-right across the grid columns).
//! Genotype storage stays on-device to match the other strategies.
//!
//! # Reference
//!
//! - Miller (2011), *Cartesian Genetic Programming* (Natural Computing
//!   Series).

use std::marker::PhantomData;

use burn::tensor::{Int, Tensor, TensorData, backend::Backend};
use rand::{Rng, RngExt};

use rlevo_core::config::{self, ConfigError, Validate};
use rlevo_core::probability::Probability;

use crate::function_set::{ArithmeticFunctionSet, FunctionSet, Symbol};
use crate::rng::{SeedPurpose, seed_stream};
use crate::strategy::{Strategy, StrategyMetrics};

/// Fixed v1 function set: arity of each opcode.
pub const FUNCTION_ARITIES: [usize; 8] = [2, 2, 2, 2, 1, 1, 1, 0];
/// Number of opcodes in the v1 function set.
pub const NUM_FUNCTIONS: usize = FUNCTION_ARITIES.len();

/// Static configuration for a [`CartesianGeneticProgramming`] run.
#[derive(Debug, Clone)]
pub struct CgpConfig {
    /// Number of offspring per generation (λ in `(1 + λ)`).
    pub lambda: usize,
    /// Number of inputs (independent variables) the program sees.
    pub n_inputs: usize,
    /// Number of grid rows.
    pub rows: usize,
    /// Number of grid columns.
    pub cols: usize,
    /// Mutation rate applied to each gene of the integer genome. Valid by
    /// construction (`[0, 1]`).
    pub mutation_rate: Probability,
    /// Levels-back parameter: how many previous columns a node can
    /// connect to. `usize::MAX` means "any previous column".
    pub levels_back: usize,
}

impl CgpConfig {
    /// Sensible defaults: 1-output, 1-row, 30-column grid, mutation
    /// rate tuned to flip ~3 genes per genome.
    #[must_use]
    pub fn default_for(n_inputs: usize) -> Self {
        let rows = 1;
        let cols = 30;
        let genes_per_node = 3; // (function, input_0, input_1)
        let output_genes = 1;
        let total_genes = rows * cols * genes_per_node + output_genes;
        #[allow(clippy::cast_precision_loss)]
        let mutation_rate = Probability::new(3.0 / total_genes as f32);
        Self {
            lambda: 4,
            n_inputs,
            rows,
            cols,
            mutation_rate,
            levels_back: usize::MAX,
        }
    }

    /// Genes per node in the genotype layout: `(function_id, input_0, input_1)`.
    pub const GENES_PER_NODE: usize = 3;
    /// Number of output genes per program (one index pointing to the node
    /// whose value is taken as the program output).
    pub const OUTPUT_GENES: usize = 1;

    /// Total genome length (nodes × 3 + outputs).
    #[must_use]
    pub fn genome_len(&self) -> usize {
        self.rows * self.cols * Self::GENES_PER_NODE + Self::OUTPUT_GENES
    }
}

impl Validate for CgpConfig {
    fn validate(&self) -> Result<(), ConfigError> {
        const C: &str = "CgpConfig";
        config::at_least(C, "lambda", self.lambda, 1)?;
        config::at_least(C, "n_inputs", self.n_inputs, 1)?;
        config::at_least(C, "rows", self.rows, 1)?;
        config::at_least(C, "cols", self.cols, 1)?;
        // `mutation_rate` is a `Probability`: valid by construction, so no
        // `in_range` check here — see ADR 0031.
        config::at_least(C, "levels_back", self.levels_back, 1)?;
        Ok(())
    }
}

/// Generation state for [`CartesianGeneticProgramming`].
#[derive(Debug, Clone)]
pub struct CgpState<B: Backend> {
    /// Parent genotype, shape `(1, genome_len)`.
    pub parent: Tensor<B, 2, Int>,
    /// Parent fitness (host-side scalar cache).
    ///
    /// `None` until the first `tell` bootstraps it. Using `Option` — rather
    /// than a `f32::NEG_INFINITY` "unset" sentinel — keeps "uninitialised"
    /// distinct from a legitimately sanitized `−∞` fitness. A `Minimize`
    /// landscape whose natural cost is `+∞` canonicalizes to `−∞` (ADR 0034);
    /// with the old sentinel that value re-triggered the bootstrap branch on
    /// the next `ask`, collapsing the `(1+λ)` loop to a single parent forever.
    pub parent_fitness: Option<f32>,
    /// Best-so-far genotype.
    pub best_genome: Option<Tensor<B, 2, Int>>,
    /// Best-so-far fitness.
    pub best_fitness: f32,
    /// Generation counter.
    pub generation: usize,
}

/// Classical Cartesian GP with `(1 + λ)` ES.
///
/// # Example
///
/// ```no_run
/// use burn::backend::Flex;
/// use rlevo_evolution::algorithms::gp_cgp::{CartesianGeneticProgramming, CgpConfig};
///
/// let strategy = CartesianGeneticProgramming::<Flex>::new();
/// let params = CgpConfig::default_for(1);
/// assert!(params.genome_len() > 0);
/// let _ = strategy;
/// ```
#[derive(Debug, Clone, Copy, Default)]
pub struct CartesianGeneticProgramming<B: Backend> {
    _backend: PhantomData<fn() -> B>,
}

impl<B: Backend> CartesianGeneticProgramming<B> {
    /// Builds a new (stateless) strategy object.
    #[must_use]
    pub fn new() -> Self {
        Self {
            _backend: PhantomData,
        }
    }

    fn sample_initial_genome(params: &CgpConfig, rng: &mut dyn Rng) -> Vec<i64> {
        let mut genome = Vec::with_capacity(params.genome_len());
        for col in 0..params.cols {
            for _row in 0..params.rows {
                #[allow(clippy::cast_possible_wrap)]
                let func = rng.random_range(0..NUM_FUNCTIONS as i64);
                let (inp0, inp1) = sample_input_pair(col, params, rng);
                genome.push(func);
                genome.push(inp0);
                genome.push(inp1);
            }
        }
        // Output gene: any node index or input index.
        let max_node_idx = params.n_inputs + params.rows * params.cols;
        #[allow(clippy::cast_possible_wrap)]
        genome.push(rng.random_range(0..max_node_idx as i64));
        genome
    }

    fn genome_to_host(genome: &Tensor<B, 2, Int>) -> Vec<i64> {
        genome
            .clone()
            .into_data()
            .into_vec::<i32>()
            .expect("genome tensor must be readable as i32")
            .into_iter()
            .map(i64::from)
            .collect()
    }
}

fn sample_input_pair(col: usize, params: &CgpConfig, rng: &mut dyn Rng) -> (i64, i64) {
    let min_col = col.saturating_sub(params.levels_back);
    let node_indices_start = params.n_inputs + min_col * params.rows;
    let node_indices_end = params.n_inputs + col * params.rows;
    let max = node_indices_end.max(params.n_inputs);
    // Allowed inputs: 0..n_inputs (graph inputs) ∪ previous nodes.
    let input_count = params.n_inputs
        + (max - params.n_inputs)
            .saturating_sub(node_indices_start.saturating_sub(params.n_inputs));
    let pool: Vec<i64> = (0..params.n_inputs)
        .chain(node_indices_start..node_indices_end)
        .map(|i| {
            #[allow(clippy::cast_possible_wrap)]
            let v = i as i64;
            v
        })
        .collect();
    let pool = if pool.is_empty() {
        #[allow(clippy::cast_possible_wrap)]
        (0..params.n_inputs as i64).collect()
    } else {
        pool
    };
    let _ = input_count;
    if pool.is_empty() {
        // No legal input candidates (only reachable when n_inputs == 0, which
        // CgpConfig::validate rejects). Defensive: emit a benign (0, 0) so the
        // primitive is total for direct callers that bypass the validating harness.
        return (0, 0);
    }
    let pick = |rng: &mut dyn Rng| -> i64 {
        let idx = rng.random_range(0..pool.len());
        pool[idx]
    };
    (pick(rng), pick(rng))
}

fn mutate_genome(genome: &mut [i64], params: &CgpConfig, rng: &mut dyn Rng) {
    let genes_per_node = CgpConfig::GENES_PER_NODE;
    let node_genes = params.rows * params.cols * genes_per_node;
    for (gene_idx, gene) in genome.iter_mut().enumerate() {
        if rng.random::<f32>() >= params.mutation_rate.get() {
            continue;
        }
        if gene_idx < node_genes {
            let within = gene_idx % genes_per_node;
            let node_idx = gene_idx / genes_per_node;
            let col = node_idx / params.rows;
            if within == 0 {
                // function
                #[allow(clippy::cast_possible_wrap)]
                {
                    *gene = rng.random_range(0..NUM_FUNCTIONS as i64);
                }
            } else {
                let (new0, new1) = sample_input_pair(col, params, rng);
                *gene = if within == 1 { new0 } else { new1 };
            }
        } else {
            // output gene
            let max_node_idx = params.n_inputs + params.rows * params.cols;
            #[allow(clippy::cast_possible_wrap)]
            {
                *gene = rng.random_range(0..max_node_idx as i64);
            }
        }
    }
}

/// Evaluates a CGP genotype at a set of input rows.
///
/// `genome` is the host-side integer genotype (length `params.genome_len()`).
/// `inputs` is a slice of `n_samples` rows, each of length `params.n_inputs`.
/// Returns one `f32` output per input row.
///
/// Out-of-range input/node indices in the genome are clamped to the
/// last buffer slot rather than panicking — this keeps fitness
/// evaluation robust to mutated-but-unrepaired genotypes. Non-finite
/// node values (e.g., `inf` from divisions or `tan`) collapse to `0.0`.
///
/// # Panics
///
/// Panics if `genome` is empty (the last gene is the output index).
#[must_use]
pub fn evaluate_cgp(genome: &[i64], params: &CgpConfig, inputs: &[Vec<f32>]) -> Vec<f32> {
    evaluate_cgp_with(genome, params, inputs, &ArithmeticFunctionSet)
}

/// Evaluates a CGP genotype against an arbitrary [`FunctionSet`].
///
/// This is the generic core of [`evaluate_cgp`]: the opcode at each node is
/// dispatched through `fs.apply` rather than a hard-coded match, so the same
/// CGP engine can run any function set. [`evaluate_cgp`] calls this with the
/// default [`ArithmeticFunctionSet`].
///
/// `fs` is taken as a concrete monomorphized `&F` (never `&dyn FunctionSet`)
/// so the `apply` dispatch inlines in the per-node × per-sample inner loop.
///
/// Each node supplies up to two argument slots (`input_0`, `input_1`). The
/// opcode's [`arity`](FunctionSet::arity) selects how many of them reach
/// `apply`: arity-2 ops receive both, arity-1 ops receive only the first, and
/// zero-arity ops (constants) receive an empty slice. Out-of-range
/// input/node/opcode ids are clamped or treated as inert (arity 0) rather than
/// panicking, keeping evaluation robust to mutated-but-unrepaired genotypes.
/// Non-finite node values collapse to `0.0`.
///
/// # Panics
///
/// Panics if `genome` is empty (the last gene is the output index).
#[must_use]
pub fn evaluate_cgp_with<F: FunctionSet>(
    genome: &[i64],
    params: &CgpConfig,
    inputs: &[Vec<f32>],
    fs: &F,
) -> Vec<f32> {
    let node_count = params.rows * params.cols;
    let n_inputs = params.n_inputs;
    #[allow(clippy::cast_sign_loss, clippy::cast_possible_truncation)]
    let output_idx = genome[genome.len() - 1] as usize;

    let mut outputs = Vec::with_capacity(inputs.len());
    let mut buf = vec![0.0_f32; n_inputs + node_count];

    for sample in inputs {
        for (i, v) in sample.iter().enumerate() {
            buf[i] = *v;
        }
        for node in 0..node_count {
            let base = node * 3;
            #[allow(clippy::cast_sign_loss, clippy::cast_possible_truncation)]
            let func = genome[base] as i32;
            #[allow(clippy::cast_sign_loss, clippy::cast_possible_truncation)]
            let a_idx = genome[base + 1] as usize;
            #[allow(clippy::cast_sign_loss, clippy::cast_possible_truncation)]
            let b_idx = genome[base + 2] as usize;
            let a = buf[a_idx.min(buf.len() - 1)];
            let b = buf[b_idx.min(buf.len() - 1)];
            let sym = Symbol::from_raw(func);
            let arity = fs.arity(sym);
            let arg_buf = [a, b];
            let v = fs.apply(sym, &arg_buf[..arity.min(arg_buf.len())]);
            buf[n_inputs + node] = crate::function_set::finite_or_zero(v);
        }
        outputs.push(buf[output_idx.min(buf.len() - 1)]);
    }

    outputs
}

impl<B: Backend> Strategy<B> for CartesianGeneticProgramming<B>
where
    B::Device: Clone,
{
    type Params = CgpConfig;
    type State = CgpState<B>;
    type Genome = Tensor<B, 2, Int>;

    /// Samples the initial parent genome by drawing random node functions and
    /// feed-forward input connections via `rng`, then uploads the genotype as
    /// a `(1, genome_len)` integer tensor.
    fn init(
        &self,
        params: &CgpConfig,
        rng: &mut dyn Rng,
        device: &<B as burn::tensor::backend::BackendTypes>::Device,
    ) -> CgpState<B> {
        debug_assert!(
            params.validate().is_ok(),
            "invalid CgpConfig reached init: {params:?}"
        );
        let genome_vec = Self::sample_initial_genome(params, rng);
        let parent = Tensor::<B, 2, Int>::from_data(
            TensorData::new(genome_vec, [1, params.genome_len()]),
            device,
        );
        CgpState {
            parent,
            parent_fitness: None,
            best_genome: None,
            best_fitness: f32::NEG_INFINITY,
            generation: 0,
        }
    }

    /// Returns the offspring population for the current generation.
    ///
    /// On the first call (parent fitness not yet set), returns the single
    /// parent genome unchanged for initial fitness evaluation.
    /// On subsequent calls, produces `params.lambda` children by cloning the
    /// parent and applying per-gene point mutation, with mutation draws taken
    /// from a deterministic `seed_stream` (host-RNG convention).
    fn ask(
        &self,
        params: &CgpConfig,
        state: &CgpState<B>,
        rng: &mut dyn Rng,
        device: &<B as burn::tensor::backend::BackendTypes>::Device,
    ) -> (Tensor<B, 2, Int>, CgpState<B>) {
        // First call: evaluate the parent as "offspring" of size 1.
        if state.parent_fitness.is_none() {
            return (state.parent.clone(), state.clone());
        }

        let mut mut_rng = seed_stream(
            rng.next_u64(),
            state.generation as u64,
            SeedPurpose::Mutation,
        );
        let parent_vec = Self::genome_to_host(&state.parent);
        let mut offspring_genomes: Vec<i64> =
            Vec::with_capacity(params.lambda * params.genome_len());
        for _ in 0..params.lambda {
            let mut child = parent_vec.clone();
            mutate_genome(&mut child, params, &mut mut_rng);
            offspring_genomes.extend(child);
        }
        #[allow(clippy::cast_possible_truncation)]
        let offspring_genomes_i32: Vec<i32> =
            offspring_genomes.into_iter().map(|v| v as i32).collect();
        let offspring = Tensor::<B, 2, Int>::from_data(
            TensorData::new(offspring_genomes_i32, [params.lambda, params.genome_len()]),
            device,
        );
        (offspring, state.clone())
    }

    /// Applies `(1+λ)` selection and returns the updated state.
    ///
    /// The canonical CGP tie-breaking rule is used: an offspring replaces the
    /// parent when its fitness is **less than or equal to** the parent's,
    /// allowing neutral mutations to accumulate and maintain genetic diversity
    /// in the inactive (non-coding) portion of the genome.
    ///
    /// The first `tell` after `init` bootstraps the parent fitness from the
    /// initial single-genome evaluation rather than running selection.
    fn tell(
        &self,
        _params: &CgpConfig,
        offspring: Tensor<B, 2, Int>,
        fitness: Tensor<B, 1>,
        mut state: CgpState<B>,
        _rng: &mut dyn Rng,
    ) -> (CgpState<B>, StrategyMetrics) {
        let fitness_host = fitness
            .into_data()
            .into_vec::<f32>()
            .expect("fitness tensor must be readable as f32");

        if fitness_host.is_empty() {
            // A generation with no offspring (only reachable when lambda == 0, which
            // CgpConfig::validate rejects). Defensive: advance the counter and emit a
            // worst-case metric without touching selection, so tell is total for
            // direct callers that bypass the validating harness.
            state.generation += 1;
            let m = StrategyMetrics::from_host_fitness(
                state.generation,
                &[f32::NEG_INFINITY],
                state.best_fitness,
            );
            state.best_fitness = m.best_fitness_ever();
            return (state, m);
        }

        if state.parent_fitness.is_none() {
            // First tell: initial parent fitness. Sanitize so a NaN seed cannot
            // masquerade as a finite parent in the later `>=` comparison.
            state.parent_fitness = Some(crate::fitness::sanitize_fitness(fitness_host[0]));
            state.generation += 1;
            update_best(&mut state, &offspring, &fitness_host);
            let m = StrategyMetrics::from_host_fitness(
                state.generation,
                &fitness_host,
                state.best_fitness,
            );
            state.best_fitness = m.best_fitness_ever();
            return (state, m);
        }

        // (1+λ): parent survives only if NO offspring strictly beats it;
        // canonical CGP uses `>=` (under the maximise convention) to break
        // ties in favor of offspring (neutral mutations accumulate).
        // Sanitize NaN → −inf (worst) so a NaN offspring can never be picked as
        // best; the raw `best_off_fit >= parent` check below then rejects it.
        let best_off_idx = fitness_host
            .iter()
            .map(|&f| crate::fitness::sanitize_fitness(f))
            .enumerate()
            .max_by(|(_, a), (_, b)| a.total_cmp(b))
            .map_or(0, |(i, _)| i);
        let best_off_fit = crate::fitness::sanitize_fitness(fitness_host[best_off_idx]);
        let parent_fit = state
            .parent_fitness
            .expect("parent_fitness is Some after the bootstrap tell");
        // `total_cmp` (not `>=`) so the comparison is well-defined even at the
        // `−∞` worst-sentinel; ties still favour the offspring (neutral drift).
        if best_off_fit.total_cmp(&parent_fit) != std::cmp::Ordering::Less {
            let device = offspring.device();
            #[allow(clippy::cast_possible_wrap, clippy::cast_possible_truncation)]
            let idx = Tensor::<B, 1, Int>::from_data(
                TensorData::new(vec![best_off_idx as i32], [1]),
                &device,
            );
            state.parent = offspring.clone().select(0, idx);
            state.parent_fitness = Some(best_off_fit);
        }

        state.generation += 1;
        update_best(&mut state, &offspring, &fitness_host);
        let m =
            StrategyMetrics::from_host_fitness(state.generation, &fitness_host, state.best_fitness);
        state.best_fitness = m.best_fitness_ever();
        (state, m)
    }

    /// Returns the best-so-far genome and its fitness, or `None` before the
    /// first `tell` call.
    fn best(&self, state: &CgpState<B>) -> Option<(Tensor<B, 2, Int>, f32)> {
        state
            .best_genome
            .as_ref()
            .map(|g| (g.clone(), state.best_fitness))
    }
}

fn update_best<B: Backend>(state: &mut CgpState<B>, pop: &Tensor<B, 2, Int>, fitness: &[f32]) {
    if fitness.is_empty() {
        return;
    }
    // Sanitize (NaN → −∞) then order with `total_cmp`: the §3 correctness floor
    // for a direct (non-harness) caller. `best_fitness` seeds at `−∞`, so a
    // legitimately sanitized `−∞` fitness is treated as the worst, not skipped.
    let sane: Vec<f32> = fitness
        .iter()
        .map(|&f| crate::fitness::sanitize_fitness(f))
        .collect();
    let mut best_idx = 0usize;
    let mut best_f = sane[0];
    for (i, &f) in sane.iter().enumerate().skip(1) {
        if f.total_cmp(&best_f) == std::cmp::Ordering::Greater {
            best_f = f;
            best_idx = i;
        }
    }
    if best_f.total_cmp(&state.best_fitness) == std::cmp::Ordering::Greater {
        let device = pop.device();
        #[allow(clippy::cast_possible_wrap, clippy::cast_possible_truncation)]
        let idx =
            Tensor::<B, 1, Int>::from_data(TensorData::new(vec![best_idx as i32], [1]), &device);
        state.best_genome = Some(pop.clone().select(0, idx));
        state.best_fitness = best_f;
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::fitness::BatchFitnessFn;
    use crate::strategy::EvolutionaryHarness;
    use burn::backend::Flex;
    type TestBackend = Flex;

    #[test]
    fn default_config_validates() {
        assert!(CgpConfig::default_for(1).validate().is_ok());
    }

    #[test]
    fn rejects_zero_rows() {
        let mut cfg = CgpConfig::default_for(1);
        cfg.rows = 0;
        assert_eq!(cfg.validate().unwrap_err().field, "rows");
    }

    /// Defensive guard for the empty input pool (issue #154, `gp_cgp` §5.2).
    /// With `n_inputs == 0` and `col == 0` both the built pool and the fallback
    /// are empty, so the raw `random_range(0..0)` would panic. A direct caller
    /// bypassing the validating harness must get a benign `(0, 0)` instead.
    #[test]
    fn sample_input_pair_empty_pool_returns_benign_zero() {
        use rand::SeedableRng;
        use rand::rngs::StdRng;

        // Construct an invalid config (n_inputs == 0) by mutating a default,
        // bypassing `validate` which would reject it.
        let mut cfg = CgpConfig::default_for(1);
        cfg.n_inputs = 0;
        let mut rng = StdRng::seed_from_u64(7);
        assert_eq!(super::sample_input_pair(0, &cfg, &mut rng), (0, 0));
    }

    /// Defensive guard for the empty `fitness_host` (issue #154, `gp_cgp` §5.3).
    /// A generation with no offspring (`lambda == 0`) makes both the bootstrap
    /// and selection paths index `fitness_host[0]` and panic. A direct caller
    /// bypassing the validating harness must instead advance the generation
    /// counter and return without touching selection or the parent fitness.
    #[test]
    fn tell_empty_fitness_does_not_panic() {
        use rand::SeedableRng;
        use rand::rngs::StdRng;

        let device = Default::default();
        let strategy = CartesianGeneticProgramming::<TestBackend>::new();

        // Invalid config (lambda == 0) via mutation, bypassing `validate`.
        let mut params = CgpConfig::default_for(1);
        params.lambda = 0;
        let mut rng = StdRng::seed_from_u64(11);
        // Build the bootstrap state directly: `init` debug-asserts the config,
        // which lambda == 0 would trip, and `tell` ignores `params` anyway.
        let parent = Tensor::<TestBackend, 2, Int>::zeros([1, params.genome_len()], &device);
        let state: CgpState<TestBackend> = CgpState {
            parent,
            parent_fitness: None,
            best_genome: None,
            best_fitness: f32::NEG_INFINITY,
            generation: 0,
        };

        let empty_fitness =
            Tensor::<TestBackend, 1>::from_data(TensorData::new(Vec::<f32>::new(), [0]), &device);
        let empty_offspring = Tensor::<TestBackend, 2, Int>::from_data(
            TensorData::new(Vec::<i32>::new(), [0, params.genome_len()]),
            &device,
        );

        let (state1, _) = strategy.tell(&params, empty_offspring, empty_fitness, state, &mut rng);
        assert_eq!(
            state1.generation, 1,
            "empty generation must advance the counter"
        );
        assert_eq!(
            state1.parent_fitness, None,
            "empty generation must not bootstrap or mutate parent fitness"
        );
    }

    /// Regression for the `is_finite()` bootstrap sentinel vs the sanitize-to-`−∞`
    /// convention (issue #132, `gp_cgp` §1.1 / ADR 0034). A canonical `−∞` parent
    /// fitness (a `Minimize` `+∞` cost canonicalizes to `−∞`) must not re-trigger
    /// the bootstrap branch: the `(1+λ)` loop has to keep emitting λ offspring.
    #[test]
    fn neg_inf_parent_fitness_does_not_collapse_lambda_loop() {
        use rand::SeedableRng;
        use rand::rngs::StdRng;

        let device = Default::default();
        let strategy = CartesianGeneticProgramming::<TestBackend>::new();
        let params = CgpConfig::default_for(1);
        let mut rng = StdRng::seed_from_u64(3);
        let state = strategy.init(&params, &mut rng, &device);

        // Bootstrap ask returns the single parent for initial evaluation.
        let (boot, next) = strategy.ask(&params, &state, &mut rng, &device);
        assert_eq!(boot.dims()[0], 1, "bootstrap ask returns the single parent");

        // Bootstrap tell with a canonical −∞ fitness. Under the old
        // `is_finite()` sentinel this left `parent_fitness` non-finite and the
        // next `ask` collapsed back to a single genome.
        let neg_inf = Tensor::<TestBackend, 1>::from_data(
            TensorData::new(vec![f32::NEG_INFINITY], [1]),
            &device,
        );
        let (state1, _) = strategy.tell(&params, boot, neg_inf, next, &mut rng);
        assert_eq!(
            state1.parent_fitness,
            Some(f32::NEG_INFINITY),
            "bootstrap must store the sanitized −∞ parent fitness, not re-arm the sentinel"
        );

        // Next ask must produce a full λ offspring population.
        let (offspring, _) = strategy.ask(&params, &state1, &mut rng, &device);
        assert_eq!(
            offspring.dims()[0],
            params.lambda,
            "post-bootstrap ask must emit λ offspring even with a −∞ parent fitness"
        );
    }

    /// `update_best` treats a sanitized `−∞` fitness as the worst value (never a
    /// champion) yet still promotes a finite winner in the same generation
    /// (issue #132, `gp_cgp` §1.1).
    #[test]
    fn update_best_treats_neg_inf_as_worst() {
        let device = Default::default();
        let parent = Tensor::<TestBackend, 2, Int>::zeros([3, 4], &device);
        let mut state: CgpState<TestBackend> = CgpState {
            parent: parent.clone(),
            parent_fitness: None,
            best_genome: None,
            best_fitness: f32::NEG_INFINITY,
            generation: 0,
        };

        // An all-`−∞` generation must not promote any champion.
        update_best(&mut state, &parent, &[f32::NEG_INFINITY; 3]);
        assert!(
            state.best_genome.is_none(),
            "an all −∞ generation must not promote a champion"
        );

        // A finite winner (index 1) is recorded despite the `−∞` neighbours.
        update_best(
            &mut state,
            &parent,
            &[f32::NEG_INFINITY, 2.5, f32::NEG_INFINITY],
        );
        approx::assert_relative_eq!(state.best_fitness, 2.5, epsilon = 1e-6);
        assert!(
            state.best_genome.is_some(),
            "a finite winner must be recorded even beside −∞ members"
        );
    }

    /// `mutate_genome` preserves genome length and, at `mutation_rate == 0.0`,
    /// is a no-op (`rng.random() >= 0.0` always holds, so every gene is skipped).
    #[test]
    fn mutate_genome_zero_rate_is_length_preserving_noop() {
        use rand::SeedableRng;
        use rand::rngs::StdRng;

        let mut params = CgpConfig::default_for(2);
        let mut rng = StdRng::seed_from_u64(5);
        let mut genome =
            CartesianGeneticProgramming::<TestBackend>::sample_initial_genome(&params, &mut rng);
        let before = genome.clone();

        params.mutation_rate = Probability::new(0.0);
        mutate_genome(&mut genome, &params, &mut rng);
        assert_eq!(genome.len(), before.len(), "length must be preserved");
        assert_eq!(genome, before, "rate 0.0 must leave the genome untouched");
    }

    /// Feed-forward / `levels_back` invariant: for a fresh genome and for a
    /// fully mutated one (`rate == 1.0`), every node input gene references only
    /// an earlier node or a graph input, and node references stay within the
    /// `levels_back` column window. Verified on a single-row grid so the node
    /// global index is `n_inputs + col`.
    #[test]
    fn genome_respects_feedforward_and_levels_back() {
        use rand::SeedableRng;
        use rand::rngs::StdRng;

        let n_inputs = 2usize;
        let cols = 8usize;
        let levels_back = 3usize;
        let mut params = CgpConfig::default_for(n_inputs);
        params.rows = 1;
        params.cols = cols;
        params.levels_back = levels_back;

        // Assert the invariant for an arbitrary host genome layout.
        let check = |genome: &[i64]| {
            for col in 0..cols {
                let base = col * CgpConfig::GENES_PER_NODE;
                let func = genome[base];
                #[allow(clippy::cast_possible_wrap)]
                let num_funcs = NUM_FUNCTIONS as i64;
                assert!(
                    (0..num_funcs).contains(&func),
                    "function id {func} out of range at col {col}"
                );
                #[allow(clippy::cast_possible_wrap)]
                let node_global = (n_inputs + col) as i64;
                #[allow(clippy::cast_possible_wrap)]
                let n_in = n_inputs as i64;
                #[allow(clippy::cast_possible_wrap)]
                let min_node_ref = (n_inputs + col.saturating_sub(levels_back)) as i64;
                for &inp in &[genome[base + 1], genome[base + 2]] {
                    assert!(
                        inp < node_global,
                        "col {col}: input {inp} must be strictly earlier than node {node_global}"
                    );
                    assert!(
                        inp < n_in || inp >= min_node_ref,
                        "col {col}: node input {inp} outside the levels_back window \
                         [{min_node_ref}, {node_global})"
                    );
                }
            }
        };

        let mut rng = StdRng::seed_from_u64(19);
        let genome =
            CartesianGeneticProgramming::<TestBackend>::sample_initial_genome(&params, &mut rng);
        check(&genome);

        // Fully mutate, then re-check: mutation must also honor the invariant.
        let mut mutated = genome.clone();
        params.mutation_rate = Probability::new(1.0);
        mutate_genome(&mut mutated, &params, &mut rng);
        assert_eq!(mutated.len(), genome.len(), "mutation preserves length");
        check(&mutated);
    }

    /// Out-of-range input and output gene indices are clamped to the last
    /// buffer slot rather than panicking, keeping evaluation robust to
    /// mutated-but-unrepaired genotypes.
    #[test]
    fn evaluate_cgp_clamps_out_of_range_indices() {
        let params = CgpConfig {
            lambda: 1,
            n_inputs: 1,
            rows: 1,
            cols: 1,
            mutation_rate: Probability::new(0.1),
            levels_back: usize::MAX,
        };
        // [func=add, in0, in1, output] with every index far out of range.
        let genome: Vec<i64> = vec![0, 999, 999, 999];
        let inputs: Vec<Vec<f32>> = vec![vec![2.0], vec![3.0]];
        let out = evaluate_cgp(&genome, &params, &inputs);
        assert_eq!(out.len(), inputs.len(), "one output per input row");
        assert!(
            out.iter().all(|v| v.is_finite()),
            "clamped evaluation must stay finite, got {out:?}"
        );
    }

    /// `best()` returns `None` before the first `tell` (no champion recorded
    /// yet), matching the documented contract.
    #[test]
    fn best_is_none_before_first_tell() {
        use rand::SeedableRng;
        use rand::rngs::StdRng;

        let device = Default::default();
        let strategy = CartesianGeneticProgramming::<TestBackend>::new();
        let params = CgpConfig::default_for(1);
        let mut rng = StdRng::seed_from_u64(0);
        let state = strategy.init(&params, &mut rng, &device);
        assert!(
            strategy.best(&state).is_none(),
            "best must be None before the first tell"
        );
    }

    /// Symbolic regression on `x² + 1` over 20 evenly spaced x ∈ [−1, 1].
    struct SymRegression {
        params: CgpConfig,
        xs: Vec<f32>,
        ys: Vec<f32>,
    }

    impl SymRegression {
        #[allow(clippy::cast_precision_loss)]
        fn new(params: CgpConfig) -> Self {
            let xs: Vec<f32> = (0..20).map(|i| -1.0 + 2.0 * (i as f32) / 19.0).collect();
            let ys: Vec<f32> = xs.iter().map(|x| x * x + 1.0).collect();
            Self { params, xs, ys }
        }
    }

    impl<B: Backend> BatchFitnessFn<B, Tensor<B, 2, Int>> for SymRegression {
        #[allow(clippy::cast_precision_loss)]
        fn evaluate_batch(
            &mut self,
            population: &Tensor<B, 2, Int>,
            device: &<B as burn::tensor::backend::BackendTypes>::Device,
        ) -> Tensor<B, 1> {
            let pop_size = population.dims()[0];
            let data: Vec<i64> = population
                .clone()
                .into_data()
                .into_vec::<i32>()
                .expect("genome host-read of a tensor this test just built")
                .into_iter()
                .map(i64::from)
                .collect();
            let gl = self.params.genome_len();
            let inputs: Vec<Vec<f32>> = self.xs.iter().map(|&x| vec![x]).collect();
            let mut fitness = Vec::with_capacity(pop_size);
            for row in 0..pop_size {
                let genome = &data[row * gl..(row + 1) * gl];
                let preds = evaluate_cgp(genome, &self.params, &inputs);
                let mse: f32 = preds
                    .iter()
                    .zip(self.ys.iter())
                    .map(|(p, y)| (p - y).powi(2))
                    .sum::<f32>()
                    / (self.ys.len() as f32);
                fitness.push(mse);
            }
            Tensor::<B, 1>::from_data(TensorData::new(fitness, [pop_size]), device)
        }

        fn sense(&self) -> rlevo_core::objective::ObjectiveSense {
            rlevo_core::objective::ObjectiveSense::Minimize
        }
    }

    #[test]
    #[allow(clippy::cast_precision_loss)]
    fn cgp_reduces_error_on_square_plus_one() {
        let device = Default::default();
        let params = CgpConfig::default_for(1);
        let landscape = SymRegression::new(params.clone());
        let initial_error = {
            // Baseline: random genome MSE on a single seed.
            use rand::SeedableRng;
            let mut rng = rand::rngs::StdRng::seed_from_u64(123);
            let genome = CartesianGeneticProgramming::<TestBackend>::sample_initial_genome(
                &params, &mut rng,
            );
            let inputs: Vec<Vec<f32>> = landscape.xs.iter().map(|&x| vec![x]).collect();
            let preds = evaluate_cgp(&genome, &params, &inputs);
            preds
                .iter()
                .zip(landscape.ys.iter())
                .map(|(p, y)| (p - y).powi(2))
                .sum::<f32>()
                / (landscape.ys.len() as f32)
        };

        let mut harness = EvolutionaryHarness::<TestBackend, _, _>::new(
            CartesianGeneticProgramming::<TestBackend>::new(),
            params,
            landscape,
            21,
            device,
            2000,
        )
        .expect("valid params");
        harness.reset();
        loop {
            if harness.step(()).done {
                break;
            }
        }
        let best = harness.latest_metrics().unwrap().best_fitness_ever();
        // CGP should substantially beat the random-genome baseline.
        assert!(
            best < initial_error,
            "CGP did not improve: best={best} initial={initial_error}"
        );
        // Bias check: ought to beat predicting a constant y=1 (mean ~= 1.33).
        assert!(best < 0.2, "expected MSE < 0.2 but got {best}");
    }
}