symclaw/crates/symclaw-bio/tests/integration.rs

//! Integration tests for symclaw-bio.

use std::collections::HashMap;
use symclaw_bio::genome::{SignedPerm, genetic_code};
use symclaw_bio::ode_model::OdeModel;
use symclaw_bio::population::{HardyWeinberg, JC69};
use symclaw_bio::reactions::{Reaction, ReactionNetwork};

// ── ODE Model ────────────────────────────────────────────────────

#[test]
fn ode_model_state_count() {
    let m = OdeModel::new(vec!["x", "y"], vec!["k"]);
    assert_eq!(m.n_states(), 2);
}

#[test]
fn ode_model_parameter_count() {
    let m = OdeModel::new(vec!["x"], vec!["k1", "k2", "k3"]);
    assert_eq!(m.n_params(), 3);
}

#[test]
fn ode_model_add_ode_and_output() {
    let mut m = OdeModel::new(vec!["x"], vec!["k"]);
    m.add_ode("x", "k * x");
    m.add_output("obs", "x");
    assert_eq!(m.n_outputs(), 1);
    assert!(m.ode_expr("x").is_ok(), "ODE expression for x should parse");
}

// ── Reaction Network ─────────────────────────────────────────────

#[test]
fn reaction_network_stoichiometry_shape() {
    // A → B: stoichiometry matrix should be 2 rows × 1 column
    let mut net = ReactionNetwork::new();
    net.add_species("A");
    net.add_species("B");
    let mut reactants = HashMap::new();
    reactants.insert("A".to_owned(), 1.0);
    let mut products = HashMap::new();
    products.insert("B".to_owned(), 1.0);
    net.add_reaction(Reaction {
        name: "r1".into(),
        reactants,
        products,
        rate: "k1".into(),
    });
    let s = net.stoichiometry_matrix();
    assert_eq!(s.len(), 2, "2 rows (species)");
    assert_eq!(s[0].len(), 1, "1 column (reaction)");
}

#[test]
fn reaction_network_deficiency_reversible() {
    // A ⇌ B: deficiency = 0
    let mut net = ReactionNetwork::new();
    net.add_species("A");
    net.add_species("B");
    let mut r1 = HashMap::new();
    r1.insert("A".to_owned(), 1.0);
    let mut p1 = HashMap::new();
    p1.insert("B".to_owned(), 1.0);
    let mut r2 = HashMap::new();
    r2.insert("B".to_owned(), 1.0);
    let mut p2 = HashMap::new();
    p2.insert("A".to_owned(), 1.0);
    net.add_reaction(Reaction {
        name: "r1".into(),
        reactants: r1,
        products: p1,
        rate: "k1".into(),
    });
    net.add_reaction(Reaction {
        name: "r2".into(),
        reactants: r2,
        products: p2,
        rate: "k2".into(),
    });
    let d = net.deficiency();
    assert_eq!(d, 0, "A⇌B deficiency = 0, got {d}");
}

// ── Population Genetics ───────────────────────────────────────────

#[test]
fn hardy_weinberg_frequencies_sum_to_one() {
    let hw = HardyWeinberg::new(0.3);
    let (aa, ab, bb) = hw.genotype_freqs();
    let total = aa + ab + bb;
    assert!((total - 1.0).abs() < 1e-10, "HW freqs sum to 1: {total}");
}

#[test]
fn hardy_weinberg_max_het_at_p_half() {
    let hw = HardyWeinberg::new(0.5);
    let (_, ab, _) = hw.genotype_freqs();
    assert!((ab - 0.5).abs() < 1e-10, "2pq at p=0.5 = 0.5, got {ab}");
}

#[test]
fn jc69_rate_matrix_symmetry() {
    let jc = JC69::new(0.01);
    let r = jc.rate_matrix();
    for (i, row_i) in r.iter().enumerate() {
        for (j, &val) in row_i.iter().enumerate() {
            if i != j {
                assert!(
                    (val - r[j][i]).abs() < 1e-12,
                    "JC69 rate[{i}][{j}]={} ≠ [{j}][{i}]={}",
                    val,
                    r[j][i]
                );
            }
        }
    }
}

#[test]
fn jc69_rate_rows_sum_to_zero() {
    let jc = JC69::new(0.02);
    let r = jc.rate_matrix();
    for (i, row) in r.iter().enumerate() {
        let sum: f64 = row.iter().sum();
        assert!(sum.abs() < 1e-10, "JC69 row {i} = {sum}, expected 0");
    }
}

// ── Genome Algebra ────────────────────────────────────────────────

#[test]
fn genetic_code_has_64_codons() {
    assert_eq!(genetic_code().len(), 64);
}

#[test]
fn genetic_code_atg_is_met() {
    // ATG → Met (either "M" or "Met" depending on encoding)
    let aa = genetic_code().get("ATG").copied().unwrap_or("");
    assert!(aa == "M" || aa == "Met", "ATG should be Met/M, got: {aa}");
}

#[test]
fn genetic_code_stop_codons() {
    let code = genetic_code();
    for stop in &["TAA", "TAG", "TGA"] {
        let aa = code.get(stop).copied().unwrap_or("");
        assert!(
            aa == "Stop" || aa == "*" || aa == "Ter",
            "{stop} should be a stop codon, got: {aa}"
        );
    }
}

#[test]
fn signed_perm_double_reversal_identity() {
    let p = SignedPerm::new(vec![1, 2, 3, 4]);
    let rev = p.reverse_segment(1, 3);
    let back = rev.reverse_segment(1, 3);
    assert_eq!(back.to_string(), p.to_string(), "rev(rev(p)) = p");
}

#[test]
fn signed_perm_identity_is_identity() {
    assert!(SignedPerm::identity(4).is_identity());
}

#[test]
fn signed_perm_breakpoints_sorted() {
    assert_eq!(SignedPerm::identity(4).breakpoints(), 0);
}