RustedSciThe 0.4.20

//! # BVP_sci Story Tests
//!
//! End-to-end correctness and performance story tests for the `BVP_sci` solver.
//!
//! ## Design Principles
//!
//! 1. **Extensibility via `BvpSciStatistics` HashMaps** — `RaceRow` has minimal hardcoded
//!    fields (source, matrix, variant, bootstrap_hint, total_ms, max_abs_solution,
//!    solve_diff, rel_x_diff, status). All solver-specific metrics are read from
//!    `BvpSciStatistics.counters`, `.timers`, `.diagnostics` at summarization time.
//!    Adding a new backend (Banded, Dense) = adding a `RaceVariant` entry, no struct changes.
//!
//! 2. **Key comparison axes** — Sparse (only backend available now) vs future Banded/Dense.
//!    AtomView vs ExprLegacy. Toolchain matrix: Lambdify (baseline), AOT/gcc, AOT/tcc,
//!    AOT/zig, AOT/rust.
//!
//! 3. **Statistics source** — `BvpSciStatistics` (the internal instrument) exclusively.
//!    No handrolled wrappers. The `stats_timer_ms()`, `stats_diagnostic_*()`, `stats_count()`
//!    helpers mirror `BVP_Damp_tests4.rs` exactly.
//!
//! 4. **Future-proof** — `RaceVariant` uses `BvpSciGeneratedBackendConfig` directly.
//!    When Banded backend is added to BVP_sci, new variants with `matrix: "Banded"` are
//!    added to the variant generator. No structural changes needed.
//!
//! 5. **Per-stage information** — `BvpSciStatistics.diagnostics` records `n_variables`,
//!    `n_mesh_points`, `matrix_nnz`, `solver_strategy`. The `timers` HashMap records
//!    per-stage timers (symbolic, residual, jacobian, linear system, grid refinement).
//!    The `counters` HashMap records iterations, linear solves, jacobian recalculations,
//!    backtracking steps.
//!
//! ## Test Index
//!
//! - `combustion_200_lambdify_baseline_story` — ExprLegacy lambdify baseline, 3 reps
//! - `combustion_200_aot_correctness_story` — All 4 AOT toolchains (gcc/tcc/zig/rust), 3 reps
//! - `combustion_1000_release_matrix_story` — Full matrix (all variants), 5 reps
//! - `combustion_1000_sparse_aot_whole_vs_chunk4_story` — Sparse AOT whole vs chunk4, 5 reps
//! - `combustion_1000_sparse_residual_aot_whole_vs_chunk4_story` — Sparse residual whole vs chunk4, 5 reps
//! - `combustion_3000_sparse_aot_chunking_story` — Heavy sparse AOT chunking matrix, 5 reps
//! - `combustion_3000_sparse_residual_aot_chunking_story` — Heavy sparse residual chunking matrix, 5 reps
//! - `combustion_200_exprlegacy_stability_story` — ExprLegacy stability, 5 reps
//! - `combustion_3000_sparse_isolated_stress_story` — Process-isolated cold stress, 3000 mesh, 2 reps

#[cfg(test)]
mod tests {
    use crate::numerical::BVP_Damp::NR_Damp_solver_damped::{
        DampedBvpStatistics, DampedSolverOptions, NRBVP as DampNRBVP, SolverParams,
    };
    use crate::numerical::BVP_Damp::generated_solver_handoff::{
        AotBuildPolicy as DampAotBuildPolicy, AotExecutionPolicy as DampAotExecutionPolicy,
        GeneratedBackendConfig as DampGeneratedBackendConfig,
    };
    use crate::numerical::BVP_sci::BVP_sci_aot::{
        BvpSciGeneratedBackendConfig, BvpSciGeneratedMatrixBackend,
    };
    use crate::numerical::BVP_sci::BVP_sci_faer::BvpSciLinearSolvePolicy;
    use crate::numerical::BVP_sci::BVP_sci_symb::{BVPwrap, BvpSciSolverOptions, BvpSciStatistics};
    use crate::symbolic::codegen::codegen_backend_selection::BackendSelectionPolicy;
    use crate::symbolic::codegen::codegen_runtime_api::ResidualChunkingStrategy;
    use crate::symbolic::codegen::codegen_tasks::{BandedChunkingStrategy, SparseChunkingStrategy};
    use crate::symbolic::symbolic_engine::Expr;
    use crate::symbolic::symbolic_functions_BVP::BvpSymbolicAssemblyBackend;
    use nalgebra::DMatrix;
    use std::collections::HashMap;
    use std::fs;
    use std::io::{self, Write};
    use std::panic::{AssertUnwindSafe, catch_unwind};
    use std::path::PathBuf;
    use std::process::Command;
    use std::thread;
    use std::time::{Duration, Instant};

    // ============================================================
    // Section A: Data structures (extensible via BvpSciStatistics)
    // ============================================================

    const CHUNKING_STORY_SCHEMA: &str = "bvp-sci-chunking-callback-stages-v2";

    #[derive(Clone)]
    struct RaceVariant {
        source: &'static str,
        matrix: &'static str,
        variant: &'static str,
        bootstrap_hint: &'static str,
        config: BvpSciGeneratedBackendConfig,
    }

    /// Minimal hardcoded fields. All solver-specific metrics are read from
    /// `BvpSciStatistics` HashMaps at summarization time via `RaceSample`.
    #[derive(Clone, Debug)]
    struct RaceRow {
        source: &'static str,
        matrix: &'static str,
        variant: &'static str,
        bootstrap_hint: &'static str,
        total_ms: f64,
        max_abs_solution: f64,
        solve_diff: f64,
        rel_x_diff: f64,
        status: String,
    }

    /// Wraps a `RaceRow` with the full `BvpSciStatistics` so that
    /// `summarize_variant` can read any diagnostic/timer/counter from the HashMaps.
    #[derive(Clone)]
    struct RaceSample {
        row: RaceRow,
        statistics: BvpSciStatistics,
    }

    #[derive(Clone, Copy, Debug)]
    struct Aggregate {
        mean: f64,
        stddev: f64,
        min: f64,
        max: f64,
    }

    #[derive(Clone, Debug)]
    struct RaceSummaryRow {
        source: &'static str,
        matrix: &'static str,
        variant: &'static str,
        bootstrap_hint: &'static str,
        runs: usize,
        ok_runs: usize,
        total_ms: Aggregate,
        max_abs_solution: Aggregate,
        solve_diff: Aggregate,
        rel_x_diff: Aggregate,
        status: String,
    }

    #[derive(Clone, Debug)]
    struct CrossSolverMemoryRow {
        solver: &'static str,
        dense_kib: f64,
        dense_mib: f64,
        sparse_kib: Option<f64>,
        sparse_mib: Option<f64>,
        nnz: usize,
        status: &'static str,
    }

    // ============================================================
    // Section B: Statistics helpers (mirror BVP_Damp_tests4.rs)
    // ============================================================

    fn stats_count(stats: &BvpSciStatistics, key: &str) -> usize {
        stats.counters.get(key).copied().unwrap_or(0)
    }

    fn stats_timer_ms(stats: &BvpSciStatistics, prefix: &str) -> f64 {
        stats
            .timers
            .iter()
            .find(|(key, _)| key.starts_with(prefix))
            .and_then(|(key, value)| timer_value_to_ms(key, value))
            .unwrap_or(f64::NAN)
    }

    fn stats_diagnostic_usize(stats: &BvpSciStatistics, key: &str) -> f64 {
        stats
            .diagnostics
            .get(key)
            .and_then(|value| value.parse::<usize>().ok())
            .or_else(|| stats.counters.get(key).copied())
            .map(|value| value as f64)
            .unwrap_or(f64::NAN)
    }

    fn stats_counter_us_as_ms(stats: &BvpSciStatistics, key: &str) -> f64 {
        stats_count(stats, key) as f64 / 1_000.0
    }
    #[allow(dead_code)]
    fn stats_diagnostic_ms(stats: &BvpSciStatistics, key: &str) -> f64 {
        stats
            .diagnostics
            .get(key)
            .and_then(|value| value.parse::<f64>().ok())
            .unwrap_or(f64::NAN)
    }

    fn stats_diagnostic_string(stats: &BvpSciStatistics, key: &str) -> String {
        stats
            .diagnostics
            .get(key)
            .cloned()
            .unwrap_or_else(|| "-".to_string())
    }

    fn damp_dense_memory_kib(stats: &DampedBvpStatistics) -> f64 {
        stats
            .counters
            .get("jacobian memory, MB")
            .copied()
            .unwrap_or(0) as f64
            * 1024.0
    }

    fn damp_jacobian_elements(stats: &DampedBvpStatistics) -> usize {
        stats
            .counters
            .get("number of jacobian elements")
            .copied()
            .unwrap_or(0)
    }

    fn timer_value_to_ms(key: &str, value: &str) -> Option<f64> {
        let duration = value
            .split(',')
            .next_back()
            .and_then(|part| part.trim().parse::<f64>().ok())?;
        let multiplier = if key.contains("ms") {
            1.0
        } else if key.contains("min") {
            60_000.0
        } else if key.contains('h') {
            3_600_000.0
        } else {
            1_000.0
        };
        Some(duration * multiplier)
    }

    // ============================================================
    // Section C: Solution helpers
    // ============================================================

    fn solution_max_abs(solution: &DMatrix<f64>) -> f64 {
        solution
            .iter()
            .copied()
            .map(f64::abs)
            .fold(0.0_f64, f64::max)
    }

    fn solution_linf_diff(lhs: &DMatrix<f64>, rhs: &DMatrix<f64>) -> f64 {
        lhs.iter()
            .zip(rhs.iter())
            .map(|(l, r)| (l - r).abs())
            .fold(0.0_f64, f64::max)
    }

    fn solution_rel_diff(lhs: &DMatrix<f64>, rhs: &DMatrix<f64>) -> f64 {
        solution_linf_diff(lhs, rhs) / solution_max_abs(rhs).max(1.0)
    }

    fn mean(values: &[f64]) -> f64 {
        if values.is_empty() {
            return f64::NAN;
        }
        values.iter().copied().sum::<f64>() / values.len() as f64
    }

    fn uniform_initial_guess(variable_count: usize, n_steps: usize, value: f64) -> DMatrix<f64> {
        DMatrix::from_element(variable_count, n_steps, value)
    }

    // ============================================================
    // Section D: Combustion problem setup
    // ============================================================

    fn combustion_problem_options(n_steps: usize) -> BvpSciSolverOptions {
        let values = vec![
            "Teta".to_string(),
            "q".to_string(),
            "C0".to_string(),
            "J0".to_string(),
            "C1".to_string(),
            "J1".to_string(),
        ];
        let teta = Expr::parse_expression("Teta");
        let q = Expr::parse_expression("q");
        let c0 = Expr::parse_expression("C0");
        let j0 = Expr::parse_expression("J0");
        let j1 = Expr::parse_expression("J1");

        let q_heat = 3000.0 * 1e3 * 0.034;
        let dt = 600.0;
        let t_scale = 600.0;
        let l: f64 = 3e-4;
        let m0 = 34.2 / 1000.0;
        let lambda = 0.07;
        let p = 2e6;
        let tm = 1500.0;
        let c1_0 = 1.0;
        let t_initial = 1000.0;
        let pe_q = 0.0090168;
        let d_ro = 2.88e-4;
        let pe_d = 1.50e-3;
        let ro_m_ = m0 * p / (8.314 * tm);

        let dt_sym = Expr::Const(dt);
        let t_scale_sym = Expr::Const(t_scale);
        let lambda_sym = Expr::Const(lambda);
        let q_heat = Expr::Const(q_heat);
        let a = Expr::Const(1.3e5);
        let e = Expr::Const(5000.0 * 4.184);
        let m = Expr::Const(m0);
        let r_g = Expr::Const(8.314);
        let ro_m = Expr::Const(ro_m_);
        let qm = Expr::Const(l.powf(2.0) / t_scale);
        let qs = Expr::Const(l.powf(2.0));
        let pe_q_sym = Expr::Const(pe_q);
        let ro_d = [Expr::Const(d_ro), Expr::Const(d_ro)];
        let pe_d = [Expr::Const(pe_d), Expr::Const(pe_d)];
        let minus = Expr::Const(-1.0);
        let m_reag = Expr::Const(0.342);

        let rate = a
            * Expr::exp(-e / (r_g * (teta.clone() * t_scale_sym + dt_sym)))
            * c0.clone()
            * (ro_m.clone() / m_reag.clone());
        let eq_system = vec![
            q.clone() / lambda_sym,
            q * pe_q_sym - q_heat * rate.clone() * qm,
            j0.clone() / ro_d[0].clone(),
            j0 * pe_d[0].clone()
                - (m.clone() * minus * rate.clone() * ro_m.clone() / m.clone()) * qs.clone(),
            j1.clone() / ro_d[1].clone(),
            j1 * pe_d[1].clone() - (m.clone() * rate * ro_m / m) * qs,
        ];

        let boundary_conditions = HashMap::from([
            ("Teta".to_string(), vec![(0, (t_initial - dt) / t_scale)]),
            ("q".to_string(), vec![(1, 1e-10)]),
            ("C0".to_string(), vec![(0, c1_0)]),
            ("J0".to_string(), vec![(1, 1e-7)]),
            ("C1".to_string(), vec![(0, 1e-3)]),
            ("J1".to_string(), vec![(1, 1e-7)]),
        ]);
        let bounds = HashMap::from([
            ("Teta".to_string(), vec![(0, 0.0), (1, 10.0)]),
            ("q".to_string(), vec![(0, -1e20), (1, 1e20)]),
            ("C0".to_string(), vec![(0, 0.0), (1, 1.5)]),
            ("J0".to_string(), vec![(0, -1e2), (1, 1e2)]),
            ("C1".to_string(), vec![(0, 0.0), (1, 1.5)]),
            ("J1".to_string(), vec![(0, -1e2), (1, 1e2)]),
        ]);

        BvpSciSolverOptions::new(
            None,
            Some(0.0),
            Some(1.0),
            Some(n_steps),
            eq_system,
            values,
            vec![],
            None,
            boundary_conditions,
            "x".to_string(),
            1e-6,
            n_steps * 2,
            uniform_initial_guess(6, n_steps, 0.99),
        )
        .with_bounds(Some(bounds))
        .with_loglevel(Some("off".to_string()))
    }

    fn exponential_endpoint_problem_options(n_steps: usize) -> BvpSciSolverOptions {
        let eq_system = Expr::parse_vector_expression(vec!["z", "-(2.0/4.0)*(1+2.0*ln((y)))*y"]);
        let values = vec!["y".to_string(), "z".to_string()];
        let a = 4.0_f64;
        let bc_val = (-1.0 / a).exp();
        let boundary_conditions =
            HashMap::from([("y".to_string(), vec![(0usize, bc_val), (1usize, bc_val)])]);
        let bounds = HashMap::from([("y".to_string(), vec![(0usize, 1e-10)])]);

        BvpSciSolverOptions::new(
            None,
            Some(-1.0),
            Some(1.0),
            Some(n_steps),
            eq_system,
            values,
            vec![],
            None,
            boundary_conditions,
            "x".to_string(),
            1e-6,
            (n_steps * 2).max(2048),
            DMatrix::zeros(2, n_steps),
        )
        .with_bounds(Some(bounds))
        .with_loglevel(Some("off".to_string()))
    }

    fn make_combustion_solver(n_steps: usize, config: BvpSciGeneratedBackendConfig) -> BVPwrap {
        let use_native_banded = config.matrix_backend == BvpSciGeneratedMatrixBackend::Banded;
        let mut options = combustion_problem_options(n_steps).with_generated_backend_config(config);
        if use_native_banded {
            options = options.with_linear_solve_policy(BvpSciLinearSolvePolicy::AutoBanded);
        }
        BVPwrap::new_with_options(options)
    }

    fn make_solver_with_options_and_linear_policy(
        options: BvpSciSolverOptions,
        config: BvpSciGeneratedBackendConfig,
        policy: BvpSciLinearSolvePolicy,
    ) -> BVPwrap {
        BVPwrap::new_with_options(
            options
                .with_generated_backend_config(config)
                .with_linear_solve_policy(policy),
        )
    }

    fn damp_combustion_problem_options(n_steps: usize) -> DampedSolverOptions {
        let bounds = HashMap::from([
            ("Teta".to_string(), (0.0, 10.0)),
            ("q".to_string(), (-1e20, 1e20)),
            ("C0".to_string(), (0.0, 1.5)),
            ("J0".to_string(), (-1e2, 1e2)),
            ("C1".to_string(), (0.0, 1.5)),
            ("J1".to_string(), (-1e2, 1e2)),
        ]);
        let rel_tolerance = HashMap::from([
            ("Teta".to_string(), 1e-5),
            ("q".to_string(), 1e-5),
            ("C0".to_string(), 1e-5),
            ("J0".to_string(), 1e-5),
            ("C1".to_string(), 1e-5),
            ("J1".to_string(), 1e-5),
        ]);
        let strategy_params = SolverParams {
            max_jac: Some(6),
            max_damp_iter: Some(6),
            damp_factor: Some(0.5),
            adaptive: None,
        };

        DampedSolverOptions::sparse_damped()
            .with_strategy_params(Some(strategy_params))
            .with_abs_tolerance(1e-6)
            .with_rel_tolerance(rel_tolerance)
            .with_max_iterations(n_steps * 2)
            .with_bounds(bounds)
            .with_loglevel(Some("off".to_string()))
    }

    fn damp_combustion_lambdify_config() -> DampGeneratedBackendConfig {
        DampGeneratedBackendConfig::default()
            .with_backend_policy_override(Some(BackendSelectionPolicy::LambdifyOnly))
            .with_symbolic_assembly_backend(BvpSymbolicAssemblyBackend::ExprLegacy)
            .with_aot_execution_policy(DampAotExecutionPolicy::SequentialOnly)
            .with_aot_build_policy(DampAotBuildPolicy::RebuildAlways {
                profile:
                    crate::numerical::BVP_Damp::generated_solver_handoff::AotBuildProfile::Release,
            })
    }

    fn make_damp_combustion_solver(
        n_steps: usize,
        config: DampGeneratedBackendConfig,
    ) -> DampNRBVP {
        let unknowns_str: Vec<&str> = vec!["Teta", "q", "C0", "J0", "C1", "J1"];
        let unknowns: Vec<Expr> = Expr::parse_vector_expression(unknowns_str.clone());
        let teta = unknowns[0].clone();
        let q = unknowns[1].clone();
        let c0 = unknowns[2].clone();
        let j0 = unknowns[3].clone();
        let j1 = unknowns[5].clone();

        let q_heat = 3000.0 * 1e3 * 0.034;
        let dt = 600.0;
        let t_scale = 600.0;
        let l: f64 = 3e-4;
        let m0 = 34.2 / 1000.0;
        let lambda = 0.07;
        let p = 2e6;
        let tm = 1500.0;
        let c1_0 = 1.0;
        let t_initial = 1000.0;
        let pe_q = 0.0090168;
        let d_ro = 2.88e-4;
        let pe_d = 1.50e-3;
        let ro_m_ = m0 * p / (8.314 * tm);

        let dt_sym = Expr::Const(dt);
        let t_scale_sym = Expr::Const(t_scale);
        let lambda_sym = Expr::Const(lambda);
        let q_heat = Expr::Const(q_heat);
        let a = Expr::Const(1.3e5);
        let e = Expr::Const(5000.0 * 4.184);
        let m = Expr::Const(m0);
        let r_g = Expr::Const(8.314);
        let ro_m = Expr::Const(ro_m_);
        let qm = Expr::Const(l.powf(2.0) / t_scale);
        let qs = Expr::Const(l.powf(2.0));
        let pe_q_sym = Expr::Const(pe_q);
        let ro_d = vec![Expr::Const(d_ro), Expr::Const(d_ro)];
        let pe_d = vec![Expr::Const(pe_d), Expr::Const(pe_d)];
        let minus = Expr::Const(-1.0);
        let m_reag = Expr::Const(0.342);

        let rate = a
            * Expr::exp(-e / (r_g * (teta.clone() * t_scale_sym + dt_sym)))
            * c0.clone()
            * (ro_m.clone() / m_reag.clone());
        let eq_system = vec![
            q.clone() / lambda_sym,
            q * pe_q_sym - q_heat * rate.clone() * qm,
            j0.clone() / ro_d[0].clone(),
            j0 * pe_d[0].clone()
                - (m.clone() * minus * rate.clone() * ro_m.clone() / m.clone()) * qs.clone(),
            j1.clone() / ro_d[1].clone(),
            j1 * pe_d[1].clone() - (m.clone() * rate * ro_m / m) * qs,
        ];

        let boundary_conditions = HashMap::from([
            ("Teta".to_string(), vec![(0, (t_initial - dt) / t_scale)]),
            ("q".to_string(), vec![(1, 1e-10)]),
            ("C0".to_string(), vec![(0, c1_0)]),
            ("J0".to_string(), vec![(1, 1e-7)]),
            ("C1".to_string(), vec![(0, 1e-3)]),
            ("J1".to_string(), vec![(1, 1e-7)]),
        ]);
        DampNRBVP::new_with_options(
            eq_system,
            DMatrix::from_element(6, n_steps, 0.99),
            unknowns_str.iter().map(|value| value.to_string()).collect(),
            "x".to_string(),
            boundary_conditions,
            0.0,
            1.0,
            n_steps,
            damp_combustion_problem_options(n_steps).with_generated_backend_config(config),
        )
    }

    // ============================================================
    // Section E: Runner functions
    // ============================================================

    fn run_race_variant(
        n_steps: usize,
        source: &'static str,
        matrix: &'static str,
        variant: &'static str,
        bootstrap_hint: &'static str,
        config: BvpSciGeneratedBackendConfig,
    ) -> (RaceSample, Option<DMatrix<f64>>) {
        let total_begin = Instant::now();
        let mut solver = make_combustion_solver(n_steps, config);
        let solve_status = catch_unwind(AssertUnwindSafe(|| solver.try_solve()));
        let total_ms = total_begin.elapsed().as_secs_f64() * 1_000.0;
        let statistics = solver.get_statistics();

        let (row, solution) = match solve_status {
            Ok(Ok(())) => match solver.get_result() {
                Some(solution) => {
                    let max_abs_solution = solution_max_abs(&solution);
                    (
                        RaceRow {
                            source,
                            matrix,
                            variant,
                            bootstrap_hint,
                            total_ms,
                            max_abs_solution,
                            solve_diff: 0.0,
                            rel_x_diff: 0.0,
                            status: "ok".to_string(),
                        },
                        Some(solution),
                    )
                }
                None => (
                    RaceRow {
                        source,
                        matrix,
                        variant,
                        bootstrap_hint,
                        total_ms,
                        max_abs_solution: f64::NAN,
                        solve_diff: f64::NAN,
                        rel_x_diff: f64::NAN,
                        status: "no_result".to_string(),
                    },
                    None,
                ),
            },
            Ok(Err(err)) => (
                RaceRow {
                    source,
                    matrix,
                    variant,
                    bootstrap_hint,
                    total_ms,
                    max_abs_solution: f64::NAN,
                    solve_diff: f64::NAN,
                    rel_x_diff: f64::NAN,
                    status: format!("err:{}", err),
                },
                None,
            ),
            Err(panic_err) => (
                RaceRow {
                    source,
                    matrix,
                    variant,
                    bootstrap_hint,
                    total_ms,
                    max_abs_solution: f64::NAN,
                    solve_diff: f64::NAN,
                    rel_x_diff: f64::NAN,
                    status: format!("panic:{:?}", panic_err),
                },
                None,
            ),
        };

        (RaceSample { row, statistics }, solution)
    }

    fn run_race_variant_with_problem_options_and_linear_policy(
        options: BvpSciSolverOptions,
        source: &'static str,
        matrix: &'static str,
        variant: &'static str,
        bootstrap_hint: &'static str,
        config: BvpSciGeneratedBackendConfig,
        policy: BvpSciLinearSolvePolicy,
    ) -> (RaceSample, Option<DMatrix<f64>>) {
        let total_begin = Instant::now();
        let mut solver = make_solver_with_options_and_linear_policy(options, config, policy);
        let solve_status = catch_unwind(AssertUnwindSafe(|| solver.try_solve()));
        let total_ms = total_begin.elapsed().as_secs_f64() * 1_000.0;
        let statistics = solver.get_statistics();

        let (row, solution) = match solve_status {
            Ok(Ok(())) => match solver.get_result() {
                Some(solution) => {
                    let max_abs_solution = solution_max_abs(&solution);
                    (
                        RaceRow {
                            source,
                            matrix,
                            variant,
                            bootstrap_hint,
                            total_ms,
                            max_abs_solution,
                            solve_diff: 0.0,
                            rel_x_diff: 0.0,
                            status: "ok".to_string(),
                        },
                        Some(solution),
                    )
                }
                None => (
                    RaceRow {
                        source,
                        matrix,
                        variant,
                        bootstrap_hint,
                        total_ms,
                        max_abs_solution: f64::NAN,
                        solve_diff: f64::NAN,
                        rel_x_diff: f64::NAN,
                        status: "no_result".to_string(),
                    },
                    None,
                ),
            },
            Ok(Err(err)) => (
                RaceRow {
                    source,
                    matrix,
                    variant,
                    bootstrap_hint,
                    total_ms,
                    max_abs_solution: f64::NAN,
                    solve_diff: f64::NAN,
                    rel_x_diff: f64::NAN,
                    status: format!("err:{}", err),
                },
                None,
            ),
            Err(panic_err) => (
                RaceRow {
                    source,
                    matrix,
                    variant,
                    bootstrap_hint,
                    total_ms,
                    max_abs_solution: f64::NAN,
                    solve_diff: f64::NAN,
                    rel_x_diff: f64::NAN,
                    status: format!("panic:{:?}", panic_err),
                },
                None,
            ),
        };

        (RaceSample { row, statistics }, solution)
    }

    fn run_race_samples(
        variants: &[RaceVariant],
        n_steps: usize,
        repetitions: usize,
    ) -> Vec<RaceSample> {
        let mut samples = Vec::with_capacity(variants.len() * repetitions);
        for repetition in 0..repetitions {
            println!(
                "[BVP_sci story] starting repetition {}/{}",
                repetition + 1,
                repetitions
            );
            let mut rows = Vec::with_capacity(variants.len());
            let mut solutions = Vec::with_capacity(variants.len());
            for variant in variants {
                println!(
                    "[BVP_sci story] running source={} matrix={} variant={} bootstrap_hint={}",
                    variant.source, variant.matrix, variant.variant, variant.bootstrap_hint
                );
                let _ = io::stdout().flush();
                let (sample, solution) = run_race_variant(
                    n_steps,
                    variant.source,
                    variant.matrix,
                    variant.variant,
                    variant.bootstrap_hint,
                    variant.config.clone(),
                );
                println!(
                    "[BVP_sci story] finished source={} matrix={} variant={} status={}",
                    sample.row.source, sample.row.matrix, sample.row.variant, sample.row.status
                );
                let _ = io::stdout().flush();
                rows.push(sample);
                solutions.push(solution);
            }
            fill_solution_diffs(&mut rows, &solutions);
            samples.extend(rows);
        }
        samples
    }

    fn fill_solution_diffs(samples: &mut [RaceSample], solutions: &[Option<DMatrix<f64>>]) {
        let baseline = solutions.iter().find_map(|solution| solution.as_ref());
        if let Some(baseline) = baseline {
            for (sample, solution) in samples.iter_mut().zip(solutions.iter()) {
                if let Some(solution) = solution {
                    sample.row.solve_diff = solution_linf_diff(solution, baseline);
                    sample.row.rel_x_diff = solution_rel_diff(solution, baseline);
                }
            }
        }
    }

    // ============================================================
    // Section F: Aggregation and summarization
    // ============================================================

    fn aggregate(values: impl IntoIterator<Item = f64>) -> Aggregate {
        let values: Vec<f64> = values.into_iter().collect();
        let n = values.len() as f64;
        if n == 0.0 {
            return Aggregate {
                mean: f64::NAN,
                stddev: f64::NAN,
                min: f64::NAN,
                max: f64::NAN,
            };
        }
        let mean = values.iter().sum::<f64>() / n;
        let variance = values.iter().map(|v| (v - mean).powi(2)).sum::<f64>() / n;
        let stddev = variance.sqrt();
        let min = values.iter().copied().fold(f64::INFINITY, f64::min);
        let max = values.iter().copied().fold(f64::NEG_INFINITY, f64::max);
        Aggregate {
            mean,
            stddev,
            min,
            max,
        }
    }
    #[allow(dead_code)]
    fn summarize_reason<'a>(values: impl IntoIterator<Item = &'a str>) -> String {
        let reasons: Vec<&str> = values.into_iter().collect();
        if reasons.is_empty() {
            return "-".to_string();
        }
        let first = reasons[0];
        if reasons.iter().all(|r| *r == first) {
            first.to_string()
        } else {
            let unique: std::collections::HashSet<&str> = reasons.into_iter().collect();
            let mut sorted: Vec<&&str> = unique.iter().collect();
            sorted.sort();
            sorted
                .iter()
                .map(|r| r.to_string())
                .collect::<Vec<_>>()
                .join("|")
        }
    }

    fn summarize_variant(variant: &RaceVariant, samples: &[RaceSample]) -> RaceSummaryRow {
        let rows: Vec<&RaceSample> = samples
            .iter()
            .filter(|s| {
                s.row.source == variant.source
                    && s.row.matrix == variant.matrix
                    && s.row.variant == variant.variant
                    && s.row.bootstrap_hint == variant.bootstrap_hint
            })
            .collect();
        let ok_runs = rows.iter().filter(|s| s.row.status == "ok").count();
        let status = if ok_runs == rows.len() {
            format!("ok {ok_runs}/{}", rows.len())
        } else {
            let first_failure = rows
                .iter()
                .find(|s| s.row.status != "ok")
                .map(|s| s.row.status.as_str())
                .unwrap_or("unknown");
            format!("ok {ok_runs}/{}, first_failure={first_failure}", rows.len())
        };

        RaceSummaryRow {
            source: variant.source,
            matrix: variant.matrix,
            variant: variant.variant,
            bootstrap_hint: variant.bootstrap_hint,
            runs: rows.len(),
            ok_runs,
            total_ms: aggregate(rows.iter().map(|s| s.row.total_ms)),
            max_abs_solution: aggregate(rows.iter().map(|s| s.row.max_abs_solution)),
            solve_diff: aggregate(rows.iter().map(|s| s.row.solve_diff)),
            rel_x_diff: aggregate(rows.iter().map(|s| s.row.rel_x_diff)),
            status,
        }
    }

    fn summarize_samples(variants: &[RaceVariant], samples: &[RaceSample]) -> Vec<RaceSummaryRow> {
        variants
            .iter()
            .map(|variant| summarize_variant(variant, samples))
            .collect()
    }

    // ============================================================
    // Section G: Table printers
    // ============================================================

    fn fmt_agg(value: Aggregate) -> String {
        if value.mean.is_finite() {
            format!(
                "{:.3} +/- {:.3} [{:.3}, {:.3}]",
                value.mean, value.stddev, value.min, value.max
            )
        } else {
            "-".to_string()
        }
    }

    fn fmt_agg_short(value: Aggregate) -> String {
        if value.mean.is_finite() {
            format!("{:.3} +/- {:.3}", value.mean, value.stddev)
        } else {
            "-".to_string()
        }
    }

    fn fmt_agg_exp(value: Aggregate) -> String {
        if value.mean.is_finite() {
            format!("{:.3e} +/- {:.1e}", value.mean, value.stddev)
        } else {
            "-".to_string()
        }
    }

    fn print_race_summary_table(title: &str, rows: &[RaceSummaryRow]) {
        println!("{title}");
        println!("[BVP_sci story] summary table: all time columns are milliseconds.");
        println!(
            "source   | matrix | variant | runs | total_ms mean+/-std [min,max] | solve_diff | rel_x_diff | max_abs_sol | status"
        );
        println!("{}", "-".repeat(190));
        for row in rows {
            println!(
                "{:<8} | {:<6} | {:<8} | {:>2}/{:<2} | {:<32} | {:<18} | {:<18} | {:<18} | {}",
                row.source,
                row.matrix,
                row.variant,
                row.ok_runs,
                row.runs,
                fmt_agg(row.total_ms),
                fmt_agg_exp(row.solve_diff),
                fmt_agg_exp(row.rel_x_diff),
                fmt_agg_exp(row.max_abs_solution),
                row.status
            );
        }
        println!();
    }

    fn print_e2e_correctness_table(title: &str, rows: &[RaceSummaryRow]) {
        println!("{title}");
        println!(
            "[BVP_sci e2e] correctness table: all solution diffs are against the first successful Lambdify baseline in each repetition."
        );
        println!(
            "source   | matrix | variant    | bootstrap_hint  | ok/runs | solve_diff mean+/-std | rel_x_diff mean+/-std | max_abs_sol mean+/-std | status"
        );
        println!("{}", "-".repeat(190));
        for row in rows {
            println!(
                "{:<8} | {:<6} | {:<10} | {:<15} | {:>2}/{:<2}  | {:<20} | {:<21} | {:<22} | {}",
                row.source,
                row.matrix,
                row.variant,
                row.bootstrap_hint,
                row.ok_runs,
                row.runs,
                fmt_agg_exp(row.solve_diff),
                fmt_agg_exp(row.rel_x_diff),
                fmt_agg_exp(row.max_abs_solution),
                row.status
            );
        }
        println!();
    }

    fn print_e2e_performance_table(title: &str, rows: &[RaceSummaryRow]) {
        println!("{title}");
        println!(
            "[BVP_sci e2e] timing/counter table: all time columns are milliseconds; counters are counts."
        );
        println!(
            "source   | matrix | variant    | bootstrap_hint  | total_ms mean+/-std [min,max] | status"
        );
        println!("{}", "-".repeat(140));
        for row in rows {
            println!(
                "{:<8} | {:<6} | {:<10} | {:<15} | {:<32} | {}",
                row.source,
                row.matrix,
                row.variant,
                row.bootstrap_hint,
                fmt_agg(row.total_ms),
                row.status
            );
        }
        println!();
    }

    fn print_e2e_stage_breakdown_table(title: &str, samples: &[RaceSample]) {
        println!("{title}");
        println!(
            "[BVP_sci e2e] stage breakdown table: symbolic/prep, residual, Jacobian, linear solve, and grid refinement totals are all milliseconds."
        );
        println!(
            "source   | matrix | variant    | bootstrap_hint  | symbolic_ms | residual_ms | jacobian_ms | linear_ms | grid_refine_ms | niter | linsys | jac_rebuilds | status"
        );
        println!("{}", "-".repeat(190));
        for sample in samples {
            println!(
                "{:<8} | {:<6} | {:<10} | {:<15} | {:<11.3} | {:<11.3} | {:<11.3} | {:<9.3} | {:<14.3} | {:>5} | {:>6} | {:>12} | {}",
                sample.row.source,
                sample.row.matrix,
                sample.row.variant,
                sample.row.bootstrap_hint,
                sample.statistics.symbolic_prepare_ms_total,
                sample.statistics.residual_ms_total,
                sample.statistics.jacobian_ms_total,
                sample.statistics.linear_system_ms_total,
                sample.statistics.grid_refinement_ms_total,
                sample.statistics.number_of_iterations,
                sample.statistics.number_of_linear_solves,
                sample.statistics.number_of_jacobian_recalculations,
                sample.row.status
            );
        }
        println!();
    }

    fn print_generated_callback_stage_table(title: &str, samples: &[RaceSample]) {
        println!("{title}");
        println!(
            "[BVP_sci AOT callback stages] wall-clock substage timers inside linked generated callbacks. These split broad solver residual_ms/jacobian_ms into argument packing, generated values, and matrix assembly."
        );
        println!(
            "source   | matrix | variant    | bootstrap_hint  | res_args | res_values | res_assembly | jac_args | jac_values | jac_assembly | param_jac | status"
        );
        println!("{}", "-".repeat(205));
        for sample in samples {
            let sparse_assembly =
                stats_timer_ms(&sample.statistics, "Callback Jacobian Matrix Assembly");
            let banded_assembly =
                stats_timer_ms(&sample.statistics, "Callback Jacobian Banded Assembly");
            let jacobian_assembly = if sparse_assembly.is_finite() {
                sparse_assembly
            } else {
                banded_assembly
            };
            println!(
                "{:<8} | {:<6} | {:<10} | {:<15} | {:>8.3} | {:>10.3} | {:>12.3} | {:>8.3} | {:>10.3} | {:>12.3} | {:>9.3} | {}",
                sample.row.source,
                sample.row.matrix,
                sample.row.variant,
                sample.row.bootstrap_hint,
                stats_timer_ms(&sample.statistics, "Callback Residual Args"),
                stats_timer_ms(&sample.statistics, "Callback Residual Values"),
                stats_timer_ms(&sample.statistics, "Callback Residual Matrix Assembly"),
                stats_timer_ms(&sample.statistics, "Callback Jacobian Args"),
                stats_timer_ms(&sample.statistics, "Callback Jacobian Values"),
                jacobian_assembly,
                stats_timer_ms(&sample.statistics, "Callback Parameter Jacobian Values"),
                sample.row.status
            );
        }
        println!();
    }

    fn print_generated_runtime_diagnostics_table(title: &str, samples: &[RaceSample]) {
        println!("{title}");
        println!(
            "[BVP_sci AOT runtime] linked generated callback facts. `actual_jobs` is the effective mesh-parallel worker count; `chunk_count` is the number of linked generated chunk symbols available for that stage."
        );
        println!(
            "source   | matrix | variant    | bootstrap_hint  | res_jobs | jac_jobs | res_chunks | jac_chunks | res_work/job | jac_work/job | res_mesh_par | jac_mesh_par | res_fallback | jac_fallback | status"
        );
        println!("{}", "-".repeat(245));
        for sample in samples {
            println!(
                "{:<8} | {:<6} | {:<10} | {:<15} | {:>8.0} | {:>8.0} | {:>10.0} | {:>10.0} | {:>12.0} | {:>12.0} | {:<12} | {:<12} | {:<12} | {:<12} | {}",
                sample.row.source,
                sample.row.matrix,
                sample.row.variant,
                sample.row.bootstrap_hint,
                stats_diagnostic_usize(&sample.statistics, "aot.runtime.residual.actual_jobs"),
                stats_diagnostic_usize(
                    &sample.statistics,
                    "aot.runtime.sparse_jacobian.actual_jobs"
                ),
                stats_diagnostic_usize(&sample.statistics, "aot.runtime.residual.chunk_count"),
                stats_diagnostic_usize(
                    &sample.statistics,
                    "aot.runtime.sparse_jacobian.chunk_count"
                ),
                stats_diagnostic_usize(&sample.statistics, "aot.runtime.residual.work_per_job"),
                stats_diagnostic_usize(
                    &sample.statistics,
                    "aot.runtime.sparse_jacobian.work_per_job"
                ),
                stats_diagnostic_string(&sample.statistics, "aot.runtime.residual.mesh_parallel"),
                stats_diagnostic_string(
                    &sample.statistics,
                    "aot.runtime.sparse_jacobian.mesh_parallel"
                ),
                stats_diagnostic_string(&sample.statistics, "aot.runtime.residual.fallback_reason"),
                stats_diagnostic_string(
                    &sample.statistics,
                    "aot.runtime.sparse_jacobian.fallback_reason"
                ),
                sample.row.status
            );
        }
        println!();
    }

    fn print_true_banded_route_table(title: &str, samples: &[RaceSample]) {
        println!("{title}");
        println!(
            "source   | matrix | variant    | direct_assembly | sparse_bypass | bordered_factor | sparse_fallback | jac_chunks | jac_jobs | jac_fallback"
        );
        println!("{}", "-".repeat(155));
        for sample in samples {
            let jacobian_prefix = if sample.row.matrix == "Banded" {
                "aot.runtime.banded_jacobian"
            } else {
                "aot.runtime.sparse_jacobian"
            };
            println!(
                "{:<8} | {:<6} | {:<10} | {:>15} | {:>13} | {:>15} | {:>15} | {:>10} | {:>8} | {}",
                sample.row.source,
                sample.row.matrix,
                sample.row.variant,
                stats_count(&sample.statistics, "bvp sci direct banded assembly calls"),
                stats_count(
                    &sample.statistics,
                    "bvp sci global sparse jacobian bypasses"
                ),
                stats_count(&sample.statistics, "bvp sci bordered factorization calls"),
                stats_count(
                    &sample.statistics,
                    "bvp sci linear backend sparse fallback solves"
                ),
                stats_diagnostic_usize(
                    &sample.statistics,
                    &format!("{jacobian_prefix}.chunk_count")
                ),
                stats_diagnostic_usize(
                    &sample.statistics,
                    &format!("{jacobian_prefix}.actual_jobs")
                ),
                stats_diagnostic_string(
                    &sample.statistics,
                    &format!("{jacobian_prefix}.fallback_reason")
                ),
            );
        }
        println!();
    }

    fn print_e2e_lifecycle_table(title: &str, rows: &[RaceSummaryRow]) {
        println!("{title}");
        println!(
            "[BVP_sci e2e] lifecycle table: source | matrix | variant | bootstrap_hint | refinements | final_grid_points | status"
        );
        println!("{}", "-".repeat(120));
        for row in rows {
            println!(
                "{:<8} | {:<6} | {:<10} | {:<15} | {}",
                row.source, row.matrix, row.variant, row.bootstrap_hint, row.status
            );
        }
        println!();
    }

    fn print_generated_backend_lifecycle_samples(title: &str, samples: &[RaceSample]) {
        println!("{title}");
        println!(
            "[BVP_sci AOT lifecycle] per-run generated backend diagnostics; RequirePrebuilt rows must report reused_linked."
        );
        println!(
            "source   | matrix | variant | bootstrap_hint              | action          | policy          | toolchain | problem_key       | status"
        );
        println!("{}", "-".repeat(165));
        for sample in samples {
            println!(
                "{:<8} | {:<6} | {:<7} | {:<27} | {:<15} | {:<15} | {:<9} | {:<17} | {}",
                sample.row.source,
                sample.row.matrix,
                sample.row.variant,
                sample.row.bootstrap_hint,
                stats_diagnostic_string(&sample.statistics, "generated_backend_action"),
                stats_diagnostic_string(&sample.statistics, "generated_backend_policy"),
                stats_diagnostic_string(&sample.statistics, "generated_backend_toolchain"),
                stats_diagnostic_string(&sample.statistics, "generated_backend_problem_key"),
                sample.row.status
            );
        }
        println!();
    }

    fn print_linear_policy_route_table(title: &str, samples: &[RaceSample]) {
        println!("{title}");
        println!(
            "[BVP_sci linear policy] route table: counters are accumulated solver statistics. AutoBanded promotes supported parameter-free endpoint-BC matrices to the bordered solver; ExperimentalBorderedBanded is the strict no-fallback diagnostic policy."
        );
        println!(
            "source   | matrix       | variant    | total_ms | sparse | sparse_fb | full_banded | bordered | extract_ms | factor_ms | solve_ms | factor_calls | solve_calls | reuse | ls | dense_kib | sparse_kib | dense/sparse | route_full | route_bordered | route_sparse | status"
        );
        println!("{}", "-".repeat(315));
        for sample in samples {
            let dense_kib =
                stats_diagnostic_usize(&sample.statistics, "global jacobian dense equivalent kib");
            let sparse_kib =
                stats_diagnostic_usize(&sample.statistics, "global jacobian sparse storage kib");
            let dense_to_sparse = stats_diagnostic_usize(
                &sample.statistics,
                "global jacobian dense to sparse permille",
            ) / 1000.0;
            println!(
                "{:<8} | {:<12} | {:<10} | {:>8.3} | {:>6} | {:>9} | {:>11} | {:>8} | {:>10.3} | {:>9.3} | {:>8.3} | {:>12} | {:>11} | {:>5} | {:>2} | {:>9.0} | {:>10.0} | {:>12.3} | {:>10} | {:>14} | {:>12} | {}",
                sample.row.source,
                sample.row.matrix,
                sample.row.variant,
                sample.row.total_ms,
                stats_count(&sample.statistics, "bvp sci linear backend sparse solves"),
                stats_count(
                    &sample.statistics,
                    "bvp sci linear backend sparse fallback solves"
                ),
                stats_count(
                    &sample.statistics,
                    "bvp sci linear backend full banded solves"
                ),
                stats_count(
                    &sample.statistics,
                    "bvp sci linear backend bordered structured solves"
                ),
                stats_counter_us_as_ms(&sample.statistics, "bvp sci bordered extraction us"),
                stats_counter_us_as_ms(&sample.statistics, "bvp sci bordered factorization us"),
                stats_counter_us_as_ms(&sample.statistics, "bvp sci bordered structured solve us"),
                stats_count(&sample.statistics, "bvp sci bordered factorization calls"),
                stats_count(
                    &sample.statistics,
                    "bvp sci bordered structured solve calls"
                ),
                stats_count(&sample.statistics, "bvp sci bordered reuse solve calls"),
                stats_count(
                    &sample.statistics,
                    "bvp sci bordered line search solve calls"
                ),
                dense_kib,
                sparse_kib,
                dense_to_sparse,
                stats_count(&sample.statistics, "bvp sci route full scalar banded"),
                stats_count(
                    &sample.statistics,
                    "bvp sci route bordered banded candidate"
                ),
                stats_count(&sample.statistics, "bvp sci route sparse fallback"),
                sample.row.status
            );
        }
        println!();
    }

    fn print_jacobian_memory_table(title: &str, samples: &[RaceSample]) {
        println!("{title}");
        println!(
            "[BVP_sci memory] Jacobian footprint is reported via dense-equivalent and sparse CSC diagnostics. This keeps the metric family aligned with BVP_Damp story conclusions."
        );
        println!(
            "source   | matrix       | variant    | dense_kib | sparse_kib | dense/sparse_ratio | nnz | status"
        );
        println!("{}", "-".repeat(110));
        for sample in samples {
            let dense_to_sparse = stats_diagnostic_usize(
                &sample.statistics,
                "global jacobian dense to sparse permille",
            ) / 1000.0;
            println!(
                "{:<8} | {:<12} | {:<10} | {:>9} | {:>10} | {:>18.3} | {:>3} | {}",
                sample.row.source,
                sample.row.matrix,
                sample.row.variant,
                stats_diagnostic_usize(&sample.statistics, "global jacobian dense equivalent kib"),
                stats_diagnostic_usize(&sample.statistics, "global jacobian sparse storage kib"),
                dense_to_sparse,
                stats_diagnostic_usize(&sample.statistics, "global jacobian nnz"),
                sample.row.status
            );
        }
        println!();
    }

    fn print_cross_solver_memory_comparison(title: &str, rows: &[CrossSolverMemoryRow]) {
        println!("{title}");
        println!(
            "[BVP_sci vs BVP_Damp memory] compare dense-equivalent Jacobian footprint on the same combustion setup; BVP_sci also reports sparse CSC storage."
        );
        println!("solver   | dense_kib | dense_mib | sparse_kib | sparse_mib | nnz     | status");
        println!("{}", "-".repeat(98));
        for row in rows {
            let sparse_kib = row
                .sparse_kib
                .map(|value| format!("{value:>10.0}"))
                .unwrap_or_else(|| format!("{:>10}", "-"));
            let sparse_mib = row
                .sparse_mib
                .map(|value| format!("{value:>10.0}"))
                .unwrap_or_else(|| format!("{:>10}", "-"));
            println!(
                "{:<8} | {:>9.0} | {:>9.0} | {} | {} | {:>7} | {}",
                row.solver,
                row.dense_kib,
                row.dense_mib,
                sparse_kib,
                sparse_mib,
                row.nnz,
                row.status
            );
        }
        println!();
    }

    // ============================================================
    // Section H: Variant generators
    // ============================================================

    fn combustion_sparse_variants() -> Vec<RaceVariant> {
        vec![
            // Lambdify baseline (ExprLegacy)
            RaceVariant {
                source: "Lambdify",
                matrix: "Sparse",
                variant: "ExprLegacy",
                bootstrap_hint: "baseline",
                config: BvpSciGeneratedBackendConfig::default_lambdify(),
            },
            // AOT toolchains (AtomView)
            RaceVariant {
                source: "AOT",
                matrix: "Sparse",
                variant: "gcc",
                bootstrap_hint: "build_if_missing",
                config: BvpSciGeneratedBackendConfig::sparse_atomview_build_if_missing_release_gcc(
                    "output_bvp_sci",
                ),
            },
            RaceVariant {
                source: "AOT",
                matrix: "Sparse",
                variant: "tcc",
                bootstrap_hint: "build_if_missing",
                config: BvpSciGeneratedBackendConfig::sparse_atomview_build_if_missing_release_tcc(
                    "output_bvp_sci",
                ),
            },
            RaceVariant {
                source: "AOT",
                matrix: "Sparse",
                variant: "zig",
                bootstrap_hint: "build_if_missing",
                config: BvpSciGeneratedBackendConfig::sparse_atomview_build_if_missing_release_zig(
                    "output_bvp_sci",
                ),
            },
            RaceVariant {
                source: "AOT",
                matrix: "Sparse",
                variant: "rust",
                bootstrap_hint: "build_if_missing",
                config: BvpSciGeneratedBackendConfig::sparse_atomview_build_if_missing_release_rust(
                    "output_bvp_sci",
                ),
            },
        ]
    }

    fn combustion_sparse_whole_vs_chunk4_variants() -> Vec<RaceVariant> {
        vec![
            RaceVariant {
                source: "Lambdify",
                matrix: "Sparse",
                variant: "ExprLegacy",
                bootstrap_hint: "baseline",
                config: BvpSciGeneratedBackendConfig::default_lambdify(),
            },
            RaceVariant {
                source: "AOT",
                matrix: "Sparse",
                variant: "tcc/whole",
                bootstrap_hint: "build_if_missing",
                config: BvpSciGeneratedBackendConfig::sparse_atomview_build_if_missing_release_tcc(
                    "output_bvp_sci_sparse_chunking",
                )
                .with_residual_chunking_strategy(ResidualChunkingStrategy::Whole)
                .with_sparse_jacobian_chunking_strategy(SparseChunkingStrategy::Whole),
            },
            RaceVariant {
                source: "AOT",
                matrix: "Sparse",
                variant: "tcc/chunk4",
                bootstrap_hint: "build_if_missing",
                config: BvpSciGeneratedBackendConfig::sparse_atomview_build_if_missing_release_tcc(
                    "output_bvp_sci_sparse_chunking",
                )
                .with_residual_chunking_strategy(ResidualChunkingStrategy::ByTargetChunkCount {
                    target_chunks: 4,
                })
                .with_sparse_jacobian_chunking_strategy(
                    SparseChunkingStrategy::ByTargetChunkCount { target_chunks: 4 },
                ),
            },
        ]
    }

    fn combustion_sparse_heavy_chunking_variants() -> Vec<RaceVariant> {
        vec![
            RaceVariant {
                source: "Lambdify",
                matrix: "Sparse",
                variant: "ExprLegacy",
                bootstrap_hint: "baseline",
                config: BvpSciGeneratedBackendConfig::default_lambdify(),
            },
            RaceVariant {
                source: "AOT",
                matrix: "Sparse",
                variant: "tcc/whole",
                bootstrap_hint: "build_if_missing",
                config: BvpSciGeneratedBackendConfig::sparse_atomview_build_if_missing_release_tcc(
                    "output_bvp_sci_sparse_chunking_heavy",
                )
                .with_residual_chunking_strategy(ResidualChunkingStrategy::Whole)
                .with_sparse_jacobian_chunking_strategy(SparseChunkingStrategy::Whole),
            },
            RaceVariant {
                source: "AOT",
                matrix: "Sparse",
                variant: "tcc/chunk4",
                bootstrap_hint: "build_if_missing",
                config: BvpSciGeneratedBackendConfig::sparse_atomview_build_if_missing_release_tcc(
                    "output_bvp_sci_sparse_chunking_heavy",
                )
                .with_residual_chunking_strategy(ResidualChunkingStrategy::ByTargetChunkCount {
                    target_chunks: 4,
                })
                .with_sparse_jacobian_chunking_strategy(
                    SparseChunkingStrategy::ByTargetChunkCount { target_chunks: 4 },
                ),
            },
            RaceVariant {
                source: "AOT",
                matrix: "Sparse",
                variant: "tcc/chunk8",
                bootstrap_hint: "build_if_missing",
                config: BvpSciGeneratedBackendConfig::sparse_atomview_build_if_missing_release_tcc(
                    "output_bvp_sci_sparse_chunking_heavy",
                )
                .with_residual_chunking_strategy(ResidualChunkingStrategy::ByTargetChunkCount {
                    target_chunks: 8,
                })
                .with_sparse_jacobian_chunking_strategy(
                    SparseChunkingStrategy::ByTargetChunkCount { target_chunks: 8 },
                ),
            },
            RaceVariant {
                source: "AOT",
                matrix: "Sparse",
                variant: "tcc/chunk12",
                bootstrap_hint: "build_if_missing",
                config: BvpSciGeneratedBackendConfig::sparse_atomview_build_if_missing_release_tcc(
                    "output_bvp_sci_sparse_chunking_heavy",
                )
                .with_residual_chunking_strategy(ResidualChunkingStrategy::ByTargetChunkCount {
                    target_chunks: 12,
                })
                .with_sparse_jacobian_chunking_strategy(
                    SparseChunkingStrategy::ByTargetChunkCount { target_chunks: 12 },
                ),
            },
        ]
    }

    fn combustion_true_banded_aot_variants() -> Vec<RaceVariant> {
        let output_dir = "output_bvp_sci_true_banded_aot";
        vec![
            RaceVariant {
                source: "Lambdify",
                matrix: "Sparse",
                variant: "ExprLegacy",
                bootstrap_hint: "baseline",
                config: BvpSciGeneratedBackendConfig::default_lambdify(),
            },
            RaceVariant {
                source: "AOT",
                matrix: "Sparse",
                variant: "tcc/whole",
                bootstrap_hint: "global_sparse",
                config: BvpSciGeneratedBackendConfig::sparse_atomview_build_if_missing_release_tcc(
                    output_dir,
                )
                .with_residual_chunking_strategy(ResidualChunkingStrategy::Whole)
                .with_sparse_jacobian_chunking_strategy(SparseChunkingStrategy::Whole),
            },
            RaceVariant {
                source: "AOT",
                matrix: "Banded",
                variant: "tcc/whole",
                bootstrap_hint: "native_bordered",
                config: BvpSciGeneratedBackendConfig::banded_atomview_build_if_missing_release_tcc(
                    output_dir,
                )
                .with_residual_chunking_strategy(ResidualChunkingStrategy::Whole)
                .with_banded_jacobian_chunking_strategy(BandedChunkingStrategy::Whole),
            },
            RaceVariant {
                source: "AOT",
                matrix: "Banded",
                variant: "tcc/chunk4",
                bootstrap_hint: "native_bordered",
                config: BvpSciGeneratedBackendConfig::banded_atomview_build_if_missing_release_tcc(
                    output_dir,
                )
                .with_residual_chunking_strategy(ResidualChunkingStrategy::ByTargetChunkCount {
                    target_chunks: 4,
                })
                .with_banded_jacobian_chunking_strategy(
                    BandedChunkingStrategy::ByTargetChunkCount { target_chunks: 4 },
                ),
            },
        ]
    }

    fn combustion_sparse_residual_whole_vs_chunk4_variants() -> Vec<RaceVariant> {
        vec![
            RaceVariant {
                source: "Lambdify",
                matrix: "Sparse",
                variant: "ExprLegacy",
                bootstrap_hint: "baseline",
                config: BvpSciGeneratedBackendConfig::default_lambdify(),
            },
            RaceVariant {
                source: "AOT",
                matrix: "Sparse",
                variant: "tcc/residual_whole",
                bootstrap_hint: "build_if_missing",
                config: BvpSciGeneratedBackendConfig::sparse_atomview_build_if_missing_release_tcc(
                    "output_bvp_sci_sparse_residual_chunking",
                )
                .with_sparse_jacobian_chunking_strategy(SparseChunkingStrategy::Whole)
                .with_residual_chunking_strategy(ResidualChunkingStrategy::Whole),
            },
            RaceVariant {
                source: "AOT",
                matrix: "Sparse",
                variant: "tcc/residual_chunk4",
                bootstrap_hint: "build_if_missing",
                config: BvpSciGeneratedBackendConfig::sparse_atomview_build_if_missing_release_tcc(
                    "output_bvp_sci_sparse_residual_chunking",
                )
                .with_sparse_jacobian_chunking_strategy(SparseChunkingStrategy::Whole)
                .with_residual_chunking_strategy(
                    ResidualChunkingStrategy::ByTargetChunkCount { target_chunks: 4 },
                ),
            },
        ]
    }

    fn combustion_sparse_residual_heavy_chunking_variants() -> Vec<RaceVariant> {
        vec![
            RaceVariant {
                source: "Lambdify",
                matrix: "Sparse",
                variant: "ExprLegacy",
                bootstrap_hint: "baseline",
                config: BvpSciGeneratedBackendConfig::default_lambdify(),
            },
            RaceVariant {
                source: "AOT",
                matrix: "Sparse",
                variant: "tcc/residual_whole",
                bootstrap_hint: "build_if_missing",
                config: BvpSciGeneratedBackendConfig::sparse_atomview_build_if_missing_release_tcc(
                    "output_bvp_sci_sparse_residual_chunking_heavy",
                )
                .with_sparse_jacobian_chunking_strategy(SparseChunkingStrategy::Whole)
                .with_residual_chunking_strategy(ResidualChunkingStrategy::Whole),
            },
            RaceVariant {
                source: "AOT",
                matrix: "Sparse",
                variant: "tcc/residual_chunk4",
                bootstrap_hint: "build_if_missing",
                config: BvpSciGeneratedBackendConfig::sparse_atomview_build_if_missing_release_tcc(
                    "output_bvp_sci_sparse_residual_chunking_heavy",
                )
                .with_sparse_jacobian_chunking_strategy(SparseChunkingStrategy::Whole)
                .with_residual_chunking_strategy(
                    ResidualChunkingStrategy::ByTargetChunkCount { target_chunks: 4 },
                ),
            },
            RaceVariant {
                source: "AOT",
                matrix: "Sparse",
                variant: "tcc/residual_chunk8",
                bootstrap_hint: "build_if_missing",
                config: BvpSciGeneratedBackendConfig::sparse_atomview_build_if_missing_release_tcc(
                    "output_bvp_sci_sparse_residual_chunking_heavy",
                )
                .with_sparse_jacobian_chunking_strategy(SparseChunkingStrategy::Whole)
                .with_residual_chunking_strategy(
                    ResidualChunkingStrategy::ByTargetChunkCount { target_chunks: 8 },
                ),
            },
        ]
    }

    fn combustion_tcc_lifecycle_variants() -> Vec<RaceVariant> {
        let output_dir = "output_bvp_sci_lifecycle";
        vec![
            RaceVariant {
                source: "Lambdify",
                matrix: "Sparse",
                variant: "ExprLegacy",
                bootstrap_hint: "baseline",
                config: BvpSciGeneratedBackendConfig::default_lambdify(),
            },
            RaceVariant {
                source: "AOT",
                matrix: "Sparse",
                variant: "tcc",
                bootstrap_hint: "build_if_missing_or_reuse",
                config: BvpSciGeneratedBackendConfig::sparse_atomview_build_if_missing_release_tcc(
                    output_dir,
                ),
            },
            RaceVariant {
                source: "AOT",
                matrix: "Sparse",
                variant: "tcc",
                bootstrap_hint: "require_prebuilt_in_process",
                config: BvpSciGeneratedBackendConfig::sparse_atomview_require_prebuilt()
                    .with_output_parent_dir(output_dir),
            },
        ]
    }

    fn combustion_linear_policy_variants() -> Vec<(RaceVariant, BvpSciLinearSolvePolicy)> {
        vec![
            (
                RaceVariant {
                    source: "Lambdify",
                    matrix: "Sparse",
                    variant: "ExprLegacy",
                    bootstrap_hint: "linear_policy_baseline",
                    config: BvpSciGeneratedBackendConfig::default_lambdify(),
                },
                BvpSciLinearSolvePolicy::Sparse,
            ),
            (
                RaceVariant {
                    source: "Lambdify",
                    matrix: "AutoBanded",
                    variant: "ExprLegacy",
                    bootstrap_hint: "auto_banded_policy",
                    config: BvpSciGeneratedBackendConfig::default_lambdify(),
                },
                BvpSciLinearSolvePolicy::AutoBanded,
            ),
            (
                RaceVariant {
                    source: "Lambdify",
                    matrix: "ExperimentalBordered",
                    variant: "ExprLegacy",
                    bootstrap_hint: "experimental_bordered_policy",
                    config: BvpSciGeneratedBackendConfig::default_lambdify(),
                },
                BvpSciLinearSolvePolicy::ExperimentalBorderedBanded,
            ),
        ]
    }

    fn env_usize_or(name: &str, default: usize) -> usize {
        std::env::var(name)
            .ok()
            .and_then(|value| value.parse::<usize>().ok())
            .filter(|value| *value > 0)
            .unwrap_or(default)
    }

    fn run_combustion_linear_policy_samples(
        n_steps: usize,
        repetitions: usize,
    ) -> (Vec<RaceVariant>, Vec<RaceSample>) {
        run_linear_policy_samples_with_options(|| combustion_problem_options(n_steps), repetitions)
    }

    fn run_linear_policy_samples_with_options(
        mut options_factory: impl FnMut() -> BvpSciSolverOptions,
        repetitions: usize,
    ) -> (Vec<RaceVariant>, Vec<RaceSample>) {
        let variant_specs = combustion_linear_policy_variants();
        let variants: Vec<RaceVariant> = variant_specs
            .iter()
            .map(|(variant, _)| variant.clone())
            .collect();
        let mut samples = Vec::with_capacity(variant_specs.len() * repetitions);

        for _rep in 0..repetitions {
            let mut rep_samples = Vec::with_capacity(variant_specs.len());
            let mut rep_solutions = Vec::with_capacity(variant_specs.len());
            for (variant, policy) in variant_specs.iter().cloned() {
                let (sample, solution) = run_race_variant_with_problem_options_and_linear_policy(
                    options_factory(),
                    variant.source,
                    variant.matrix,
                    variant.variant,
                    variant.bootstrap_hint,
                    variant.config,
                    policy,
                );
                rep_samples.push(sample);
                rep_solutions.push(solution);
            }
            fill_solution_diffs(&mut rep_samples, &rep_solutions);
            samples.extend(rep_samples);
        }

        (variants, samples)
    }

    fn assert_linear_policy_contract(summary: &[RaceSummaryRow], samples: &[RaceSample]) {
        for row in summary {
            assert!(
                row.ok_runs == row.runs,
                "{} {} {}: expected every linear-policy run to pass, got {}",
                row.source,
                row.matrix,
                row.variant,
                row.status
            );
            assert!(
                row.solve_diff.mean < 1e-8,
                "{} {} {}: solve_diff too large: {:e}",
                row.source,
                row.matrix,
                row.variant,
                row.solve_diff.mean
            );
        }

        let auto_samples: Vec<&RaceSample> = samples
            .iter()
            .filter(|sample| sample.row.matrix == "AutoBanded")
            .collect();
        assert!(
            auto_samples.iter().all(|sample| stats_count(
                &sample.statistics,
                "bvp sci route bordered banded candidate"
            ) > 0),
            "AutoBanded should detect bordered-banded endpoint-BC structure"
        );
        assert!(
            auto_samples.iter().all(|sample| stats_count(
                &sample.statistics,
                "bvp sci linear backend full banded solves"
            ) == 0),
            "AutoBanded must not force full scalar banded on endpoint-BC matrices"
        );
        assert!(
            auto_samples.iter().all(|sample| stats_count(
                &sample.statistics,
                "bvp sci linear backend bordered structured solves"
            ) > 0),
            "AutoBanded should promote parameter-free endpoint-BC systems to bordered solves"
        );
        assert!(
            auto_samples.iter().all(|sample| stats_count(
                &sample.statistics,
                "bvp sci linear backend sparse fallback solves"
            ) == 0),
            "successful AutoBanded bordered solves should not fall back to Sparse"
        );

        let experimental_samples: Vec<&RaceSample> = samples
            .iter()
            .filter(|sample| sample.row.matrix == "ExperimentalBordered")
            .collect();
        assert!(
            experimental_samples.iter().all(|sample| stats_count(
                &sample.statistics,
                "bvp sci route bordered banded candidate"
            ) > 0),
            "ExperimentalBordered should detect bordered-banded endpoint-BC structure"
        );
        assert!(
            experimental_samples.iter().all(|sample| stats_count(
                &sample.statistics,
                "bvp sci linear backend bordered structured solves"
            ) > 0),
            "ExperimentalBordered must use the structured bordered solver"
        );
        assert!(
            experimental_samples.iter().all(|sample| stats_count(
                &sample.statistics,
                "bvp sci bordered extraction calls"
            ) > 0),
            "ExperimentalBordered must report bordered extraction diagnostics"
        );
        assert!(
            experimental_samples.iter().all(|sample| stats_count(
                &sample.statistics,
                "bvp sci bordered structured solve calls"
            ) > 0),
            "ExperimentalBordered must report bordered structured solve diagnostics"
        );
        assert!(
            experimental_samples.iter().all(|sample| stats_count(
                &sample.statistics,
                "bvp sci bordered factorization calls"
            ) > 0),
            "ExperimentalBordered must factor the bordered structured system"
        );
        assert!(
            experimental_samples.iter().all(|sample| stats_count(
                &sample.statistics,
                "bvp sci bordered structured solve us"
            ) > 0),
            "ExperimentalBordered must report non-zero bordered solve timing"
        );
        assert!(
            experimental_samples.iter().all(|sample| stats_count(
                &sample.statistics,
                "bvp sci linear backend sparse fallback solves"
            ) == 0),
            "ExperimentalBordered must not silently fall back to Sparse"
        );
    }

    // ============================================================
    // Section I: Story tests
    // ============================================================

    /// I.1: Combustion_200 — Lambdify baseline (ExprLegacy)
    ///
    /// Verifies that the ExprLegacy lambdify workflow converges on the combustion
    /// problem with 200 mesh points. Reports correctness and performance tables.
    /// This serves as the baseline for all AOT comparisons.
    #[test]
    fn combustion_200_lambdify_baseline_story() {
        let variants: Vec<RaceVariant> = combustion_sparse_variants()
            .into_iter()
            .filter(|v| v.source == "Lambdify")
            .collect();
        let samples = run_race_samples(&variants, 200, 3);
        let summary = summarize_samples(&variants, &samples);
        print_e2e_correctness_table("Combustion 200: Lambdify baseline (ExprLegacy)", &summary);
        print_e2e_performance_table("Combustion 200: Lambdify performance", &summary);
        // Assert: Lambdify baseline converges
        for row in &summary {
            assert!(
                row.ok_runs > 0,
                "{} {} {}: all runs failed",
                row.source,
                row.matrix,
                row.variant
            );
        }
    }

    /// I.2: Combustion_200 — AOT correctness (gcc, tcc, zig, rust)
    ///
    /// Verifies that all AOT toolchains produce correct solutions on the
    /// combustion problem with 200 mesh points. Reports correctness,
    /// performance, and lifecycle tables.
    #[test]
    fn combustion_200_aot_correctness_story() {
        let variants: Vec<RaceVariant> = combustion_sparse_variants()
            .into_iter()
            .filter(|v| v.source == "AOT")
            .collect();
        let samples = run_race_samples(&variants, 200, 3);
        let summary = summarize_samples(&variants, &samples);
        print_e2e_correctness_table("Combustion 200: AOT correctness (all toolchains)", &summary);
        print_e2e_performance_table("Combustion 200: AOT performance", &summary);
        print_e2e_lifecycle_table("Combustion 200: AOT lifecycle", &summary);
        // Assert: all AOT variants converge
        for row in &summary {
            assert!(
                row.ok_runs > 0,
                "{} {} {}: all runs failed",
                row.source,
                row.matrix,
                row.variant
            );
            if row.solve_diff.mean.is_finite() {
                assert!(
                    row.solve_diff.mean < 1e-6,
                    "{} {} {}: solve_diff too large: {:e}",
                    row.source,
                    row.matrix,
                    row.variant,
                    row.solve_diff.mean
                );
            }
        }
    }

    /// I.3: Combustion_1000 — Full release matrix (all variants, 5 reps)
    ///
    /// Runs all variants (Lambdify + AOT) on the combustion problem with
    /// 1000 mesh points and 5 repetitions. Reports summary, correctness,
    /// and performance tables.
    #[test]
    fn combustion_1000_release_matrix_story() {
        let variants = combustion_sparse_variants();
        let samples = run_race_samples(&variants, 1000, 5);
        let summary = summarize_samples(&variants, &samples);
        print_race_summary_table("Combustion 1000: Full release matrix (5 reps)", &summary);
        print_e2e_correctness_table("Combustion 1000: Correctness", &summary);
        print_e2e_performance_table("Combustion 1000: Performance", &summary);
        // Assert: all variants converge
        for row in &summary {
            assert!(
                row.ok_runs > 0,
                "{} {} {}: all runs failed",
                row.source,
                row.matrix,
                row.variant
            );
        }
    }

    /// I.3b: Combustion_1000 — Sparse AOT whole vs chunk4 release comparison
    ///
    /// Compares a Lambdify baseline with two sparse AOT variants using the same
    /// toolchain but different sparse Jacobian chunking strategies. This is the
    /// first release story for sparse chunking on BVP_sci.
    #[test]
    #[ignore]
    fn combustion_1000_sparse_aot_whole_vs_chunk4_story() {
        let variants = combustion_sparse_whole_vs_chunk4_variants();
        let samples = run_race_samples(&variants, 1000, 5);
        let summary = summarize_samples(&variants, &samples);

        print_race_summary_table(
            "Combustion 1000: Sparse AOT whole vs chunk4 (5 reps)",
            &summary,
        );
        print_e2e_correctness_table(
            "Combustion 1000: Sparse AOT whole vs chunk4 correctness",
            &summary,
        );
        print_e2e_performance_table(
            "Combustion 1000: Sparse AOT whole vs chunk4 timing",
            &summary,
        );
        print_generated_callback_stage_table(
            "Combustion 1000: Sparse AOT whole vs chunk4 callback stages",
            &samples,
        );
        print_generated_runtime_diagnostics_table(
            "Combustion 1000: Sparse AOT whole vs chunk4 runtime diagnostics",
            &samples,
        );
        print_generated_backend_lifecycle_samples(
            "Combustion 1000: Sparse AOT whole vs chunk4 generated backend actions",
            &samples,
        );

        for row in &summary {
            assert!(
                row.ok_runs > 0,
                "{} {} {}: all runs failed",
                row.source,
                row.matrix,
                row.variant
            );
        }
        assert!(
            samples.iter().any(|sample| {
                sample.row.source == "AOT"
                    && sample.row.variant.contains("chunk4")
                    && stats_diagnostic_usize(
                        &sample.statistics,
                        "aot.runtime.residual.actual_jobs",
                    ) > 1.0
                    && stats_diagnostic_string(
                        &sample.statistics,
                        "aot.runtime.residual.fallback_reason",
                    ) == "none"
                    && stats_diagnostic_usize(
                        &sample.statistics,
                        "aot.runtime.sparse_jacobian.actual_jobs",
                    ) > 1.0
                    && stats_diagnostic_string(
                        &sample.statistics,
                        "aot.runtime.sparse_jacobian.fallback_reason",
                    ) == "none"
            }),
            "chunk4 BVP_sci AOT story must prove real residual and sparse-Jacobian runtime parallelism"
        );
    }

    /// I.3b1: Combustion_1000 — Sparse residual AOT whole vs chunk4 release comparison
    ///
    /// Mirrors the sparse-Jacobian story, but isolates residual chunking so we
    /// can prove the residual runtime path itself is mesh-parallel.
    #[test]
    #[ignore]
    fn combustion_1000_sparse_residual_aot_whole_vs_chunk4_story() {
        let variants = combustion_sparse_residual_whole_vs_chunk4_variants();
        let samples = run_race_samples(&variants, 1000, 5);
        let summary = summarize_samples(&variants, &samples);

        print_race_summary_table(
            "Combustion 1000: Sparse residual AOT whole vs chunk4 (5 reps)",
            &summary,
        );
        print_e2e_correctness_table(
            "Combustion 1000: Sparse residual AOT whole vs chunk4 correctness",
            &summary,
        );
        print_e2e_performance_table(
            "Combustion 1000: Sparse residual AOT whole vs chunk4 timing",
            &summary,
        );
        print_generated_callback_stage_table(
            "Combustion 1000: Sparse residual AOT whole vs chunk4 callback stages",
            &samples,
        );
        print_generated_runtime_diagnostics_table(
            "Combustion 1000: Sparse residual AOT whole vs chunk4 runtime diagnostics",
            &samples,
        );
        print_generated_backend_lifecycle_samples(
            "Combustion 1000: Sparse residual AOT whole vs chunk4 generated backend actions",
            &samples,
        );

        for row in &summary {
            assert!(
                row.ok_runs > 0,
                "{} {} {}: all runs failed",
                row.source,
                row.matrix,
                row.variant
            );
        }
        assert!(
            samples.iter().any(|sample| {
                sample.row.source == "AOT"
                    && sample.row.variant.contains("chunk4")
                    && stats_diagnostic_usize(
                        &sample.statistics,
                        "aot.runtime.residual.actual_jobs",
                    ) > 1.0
                    && stats_diagnostic_string(
                        &sample.statistics,
                        "aot.runtime.residual.fallback_reason",
                    ) == "none"
                    && stats_diagnostic_usize(
                        &sample.statistics,
                        "aot.runtime.sparse_jacobian.actual_jobs",
                    ) == 1.0
                    && stats_diagnostic_string(
                        &sample.statistics,
                        "aot.runtime.sparse_jacobian.fallback_reason",
                    ) == "single_chunk"
            }),
            "chunk4 BVP_sci AOT story must prove real residual runtime parallelism"
        );
    }

    /// I.3c: Combustion_3000 — Sparse AOT chunking matrix on the 12-core source of truth
    ///
    /// This is the heavier sparse chunking story.  It keeps the Lambdify
    /// baseline, then compares `tcc/whole` against several chunk counts that are
    /// plausible on a 12-core machine.  The intent is to see whether the larger
    /// combustion solve reaches a real break-even point for more aggressive
    /// chunking.
    #[test]
    #[ignore]
    fn combustion_3000_sparse_aot_chunking_story() {
        println!("[BVP_sci story schema] {CHUNKING_STORY_SCHEMA}");
        let variants = combustion_sparse_heavy_chunking_variants();
        let samples = run_race_samples(&variants, 3000, 5);
        let summary = summarize_samples(&variants, &samples);

        print_race_summary_table(
            "Combustion 3000: Sparse AOT chunking matrix (5 reps)",
            &summary,
        );
        print_e2e_correctness_table("Combustion 3000: Sparse AOT chunking correctness", &summary);
        print_e2e_performance_table("Combustion 3000: Sparse AOT chunking timing", &summary);
        print_e2e_stage_breakdown_table(
            "Combustion 3000: Sparse AOT chunking stage breakdown",
            &samples,
        );
        print_generated_callback_stage_table(
            "Combustion 3000: Sparse AOT chunking callback stages",
            &samples,
        );
        print_generated_runtime_diagnostics_table(
            "Combustion 3000: Sparse AOT chunking runtime diagnostics",
            &samples,
        );
        print_generated_backend_lifecycle_samples(
            "Combustion 3000: Sparse AOT chunking generated backend actions",
            &samples,
        );

        for row in &summary {
            assert!(
                row.ok_runs > 0,
                "{} {} {}: all runs failed",
                row.source,
                row.matrix,
                row.variant
            );
        }
        assert!(
            samples.iter().any(|sample| {
                sample.row.source == "AOT"
                    && stats_diagnostic_usize(
                        &sample.statistics,
                        "aot.runtime.residual.actual_jobs",
                    ) > 1.0
                    && stats_diagnostic_string(
                        &sample.statistics,
                        "aot.runtime.residual.fallback_reason",
                    ) == "none"
                    && stats_diagnostic_usize(
                        &sample.statistics,
                        "aot.runtime.sparse_jacobian.actual_jobs",
                    ) > 1.0
                    && stats_diagnostic_string(
                        &sample.statistics,
                        "aot.runtime.sparse_jacobian.fallback_reason",
                    ) == "none"
            }),
            "chunked BVP_sci AOT story must prove real residual and sparse-Jacobian runtime parallelism"
        );
    }

    /// I.3c2: Combustion_1000 — true Banded AOT versus the global Sparse route.
    ///
    /// Unlike the older AutoBanded stories, the Banded rows compile a native
    /// diagonal-major pointwise Jacobian and assemble bordered blocks directly.
    /// The route table is a hard proof: Banded rows must report direct assembly,
    /// global sparse bypasses and bordered factorization without Sparse fallback.
    #[test]
    #[ignore]
    fn combustion_1000_true_banded_aot_vs_sparse_story() {
        println!("[BVP_sci story schema] true-banded-aot-v1");
        let variants = combustion_true_banded_aot_variants();
        let samples = run_race_samples(&variants, 1_000, 3);
        let summary = summarize_samples(&variants, &samples);

        print_race_summary_table(
            "Combustion 1000: Lambdify / Sparse AOT / true Banded AOT (3 reps)",
            &summary,
        );
        print_e2e_correctness_table("Combustion 1000: true Banded AOT correctness", &summary);
        print_e2e_performance_table(
            "Combustion 1000: true Banded AOT wall-clock and solver stages",
            &summary,
        );
        print_e2e_stage_breakdown_table(
            "Combustion 1000: true Banded AOT stage breakdown",
            &samples,
        );
        print_generated_callback_stage_table(
            "Combustion 1000: true Banded AOT callback stages",
            &samples,
        );
        print_true_banded_route_table(
            "Combustion 1000: proof of native Banded AOT route",
            &samples,
        );
        print_generated_backend_lifecycle_samples(
            "Combustion 1000: true Banded AOT artifact lifecycle",
            &samples,
        );

        for row in &summary {
            assert_eq!(
                row.ok_runs, row.runs,
                "{} {} {}: every story repetition must solve",
                row.source, row.matrix, row.variant
            );
            assert!(
                row.solve_diff.mean <= 1e-7,
                "{} {} {}: solution diff is too large: {:.3e}",
                row.source,
                row.matrix,
                row.variant,
                row.solve_diff.mean
            );
        }

        let banded_samples = samples
            .iter()
            .filter(|sample| sample.row.matrix == "Banded")
            .collect::<Vec<_>>();
        assert!(!banded_samples.is_empty());
        for sample in banded_samples {
            assert!(
                stats_count(&sample.statistics, "bvp sci direct banded assembly calls") > 0,
                "{} must assemble bordered blocks directly",
                sample.row.variant
            );
            assert!(
                stats_count(
                    &sample.statistics,
                    "bvp sci global sparse jacobian bypasses"
                ) > 0,
                "{} must bypass the global Sparse Jacobian",
                sample.row.variant
            );
            assert!(
                stats_count(&sample.statistics, "bvp sci bordered factorization calls") > 0,
                "{} must use native bordered factorization",
                sample.row.variant
            );
            assert_eq!(
                stats_count(
                    &sample.statistics,
                    "bvp sci linear backend sparse fallback solves"
                ),
                0,
                "{} must not silently fall back to Sparse LU",
                sample.row.variant
            );
        }
        assert!(
            samples.iter().any(|sample| {
                sample.row.matrix == "Banded"
                    && sample.row.variant == "tcc/chunk4"
                    && stats_diagnostic_usize(
                        &sample.statistics,
                        "aot.runtime.banded_jacobian.chunk_count",
                    ) > 1.0
                    && stats_diagnostic_usize(
                        &sample.statistics,
                        "aot.runtime.banded_jacobian.actual_jobs",
                    ) > 1.0
                    && stats_diagnostic_string(
                        &sample.statistics,
                        "aot.runtime.banded_jacobian.fallback_reason",
                    ) == "none"
            }),
            "chunk4 row must prove real native Banded Jacobian parallelism"
        );
    }

    /// I.3c3: Combustion_3000 — true Banded AOT versus the global Sparse route.
    ///
    /// This is the heavy release-oriented follow-up to PS.8a3.  It uses the same
    /// true banded route proof, but on a much larger mesh so we can see whether
    /// the native bordered assembly keeps its advantage when the problem size is
    /// large enough for Jacobian work to matter more.
    #[test]
    #[ignore]
    fn combustion_3000_true_banded_aot_vs_sparse_story() {
        println!("[BVP_sci story schema] true-banded-aot-v1");
        let variants = combustion_true_banded_aot_variants();
        let samples = run_race_samples(&variants, 3_000, 3);
        let summary = summarize_samples(&variants, &samples);

        print_race_summary_table(
            "Combustion 3000: Lambdify / Sparse AOT / true Banded AOT (3 reps)",
            &summary,
        );
        print_e2e_correctness_table("Combustion 3000: true Banded AOT correctness", &summary);
        print_e2e_performance_table(
            "Combustion 3000: true Banded AOT wall-clock and solver stages",
            &summary,
        );
        print_e2e_stage_breakdown_table(
            "Combustion 3000: true Banded AOT stage breakdown",
            &samples,
        );
        print_generated_callback_stage_table(
            "Combustion 3000: true Banded AOT callback stages",
            &samples,
        );
        print_true_banded_route_table(
            "Combustion 3000: proof of native Banded AOT route",
            &samples,
        );
        print_generated_backend_lifecycle_samples(
            "Combustion 3000: true Banded AOT artifact lifecycle",
            &samples,
        );

        for row in &summary {
            assert_eq!(
                row.ok_runs, row.runs,
                "{} {} {}: every story repetition must solve",
                row.source, row.matrix, row.variant
            );
            assert!(
                row.solve_diff.mean <= 1e-7,
                "{} {} {}: solution diff is too large: {:.3e}",
                row.source,
                row.matrix,
                row.variant,
                row.solve_diff.mean
            );
        }
        assert!(
            samples.iter().any(|sample| {
                sample.row.source == "AOT"
                    && sample.row.matrix == "Banded"
                    && sample.row.variant.contains("chunk4")
                    && stats_diagnostic_usize(
                        &sample.statistics,
                        "bvp sci direct banded assembly calls",
                    ) > 0.0
                    && stats_diagnostic_usize(
                        &sample.statistics,
                        "bvp sci global sparse jacobian bypasses",
                    ) > 0.0
                    && stats_diagnostic_usize(
                        &sample.statistics,
                        "bvp sci bordered factorization calls",
                    ) > 0.0
            }),
            "true Banded AOT chunk4 must report direct assembly, sparse bypass, and bordered factorization"
        );
    }

    /// I.3c1: Combustion_3000 — Sparse residual AOT chunking matrix on the 12-core source of truth
    ///
    /// This mirrors the sparse-Jacobian heavy story, but keeps the Jacobian
    /// whole and varies residual chunking only.  It proves that the residual
    /// runtime path itself can go mesh-parallel on the larger combustion case.
    #[test]
    #[ignore]
    fn combustion_3000_sparse_residual_aot_chunking_story() {
        println!("[BVP_sci story schema] {CHUNKING_STORY_SCHEMA}");
        let variants = combustion_sparse_residual_heavy_chunking_variants();
        let samples = run_race_samples(&variants, 3000, 5);
        let summary = summarize_samples(&variants, &samples);

        print_race_summary_table(
            "Combustion 3000: Sparse residual AOT chunking matrix (5 reps)",
            &summary,
        );
        print_e2e_correctness_table(
            "Combustion 3000: Sparse residual AOT chunking correctness",
            &summary,
        );
        print_e2e_performance_table(
            "Combustion 3000: Sparse residual AOT chunking timing",
            &summary,
        );
        print_e2e_stage_breakdown_table(
            "Combustion 3000: Sparse residual AOT chunking stage breakdown",
            &samples,
        );
        print_generated_callback_stage_table(
            "Combustion 3000: Sparse residual AOT chunking callback stages",
            &samples,
        );
        print_generated_runtime_diagnostics_table(
            "Combustion 3000: Sparse residual AOT chunking runtime diagnostics",
            &samples,
        );
        print_generated_backend_lifecycle_samples(
            "Combustion 3000: Sparse residual AOT chunking generated backend actions",
            &samples,
        );

        for row in &summary {
            assert!(
                row.ok_runs > 0,
                "{} {} {}: all runs failed",
                row.source,
                row.matrix,
                row.variant
            );
        }
        assert!(
            samples.iter().any(|sample| {
                sample.row.source == "AOT"
                    && sample.row.variant.contains("chunk")
                    && stats_diagnostic_usize(
                        &sample.statistics,
                        "aot.runtime.residual.actual_jobs",
                    ) > 1.0
                    && stats_diagnostic_string(
                        &sample.statistics,
                        "aot.runtime.residual.fallback_reason",
                    ) == "none"
                    && stats_diagnostic_usize(
                        &sample.statistics,
                        "aot.runtime.sparse_jacobian.actual_jobs",
                    ) == 1.0
                    && stats_diagnostic_string(
                        &sample.statistics,
                        "aot.runtime.sparse_jacobian.fallback_reason",
                    ) == "single_chunk"
            }),
            "heavy residual chunking BVP_sci AOT story must prove real residual runtime parallelism"
        );
    }

    /// I.3d: Cross-solver Jacobian memory story for combustion-200.
    ///
    /// Compares the BVP_sci sparse/global-Jacobian footprint against the BVP_Damp
    /// dense-equivalent Jacobian memory telemetry on the same combustion setup.
    /// The intent is to keep the two solver families aligned on the same
    /// production-sized problem while preserving their distinct memory models.
    #[test]
    #[ignore]
    fn combustion_200_bvp_sci_vs_bvp_damp_jacobian_memory_story() {
        let n_steps = 200usize;
        let repetitions = 3usize;

        let mut sci_dense_kib = Vec::new();
        let mut sci_dense_mib = Vec::new();
        let mut sci_sparse_kib = Vec::new();
        let mut sci_sparse_mib = Vec::new();
        let mut sci_nnz = Vec::new();
        let mut damp_dense_kib = Vec::new();
        let mut damp_dense_mib = Vec::new();
        let mut damp_nnz = Vec::new();

        for rep in 0..repetitions {
            println!(
                "[BVP_sci vs BVP_Damp memory] repetition {}/{} for combustion-{}",
                rep + 1,
                repetitions,
                n_steps
            );

            let mut sci_solver =
                make_combustion_solver(n_steps, BvpSciGeneratedBackendConfig::default_lambdify());
            sci_solver
                .try_solve()
                .expect("BVP_sci combustion memory story should solve");
            let sci_stats = sci_solver.get_statistics();
            sci_dense_kib.push(stats_diagnostic_usize(
                &sci_stats,
                "global jacobian dense equivalent kib",
            ));
            sci_dense_mib.push(stats_diagnostic_usize(
                &sci_stats,
                "global jacobian dense equivalent mib",
            ));
            sci_sparse_kib.push(stats_diagnostic_usize(
                &sci_stats,
                "global jacobian sparse storage kib",
            ));
            sci_sparse_mib.push(stats_diagnostic_usize(
                &sci_stats,
                "global jacobian sparse storage mib",
            ));
            sci_nnz.push(stats_diagnostic_usize(&sci_stats, "global jacobian nnz"));

            let mut damp_solver =
                make_damp_combustion_solver(n_steps, damp_combustion_lambdify_config());
            damp_solver
                .try_solve()
                .expect("BVP_Damp combustion memory story should solve");
            let damp_stats = damp_solver.get_statistics();
            damp_dense_kib.push(damp_dense_memory_kib(&damp_stats));
            damp_dense_mib.push(damp_dense_memory_kib(&damp_stats) / 1024.0);
            damp_nnz.push(damp_jacobian_elements(&damp_stats) as f64);
        }

        let rows = vec![
            CrossSolverMemoryRow {
                solver: "BVP_sci",
                dense_kib: mean(&sci_dense_kib),
                dense_mib: mean(&sci_dense_mib),
                sparse_kib: Some(mean(&sci_sparse_kib)),
                sparse_mib: Some(mean(&sci_sparse_mib)),
                nnz: mean(&sci_nnz).round() as usize,
                status: "ok",
            },
            CrossSolverMemoryRow {
                solver: "BVP_Damp",
                dense_kib: mean(&damp_dense_kib),
                dense_mib: mean(&damp_dense_mib),
                sparse_kib: None,
                sparse_mib: None,
                nnz: mean(&damp_nnz).round() as usize,
                status: "ok",
            },
        ];

        print_cross_solver_memory_comparison(
            &format!(
                "Combustion {n_steps}: BVP_sci vs BVP_Damp Jacobian memory ({repetitions} runs)"
            ),
            &rows,
        );

        let sci_sparse_vs_damp_dense = mean(&sci_sparse_kib) / mean(&damp_dense_kib).max(1.0);
        let sci_dense_vs_damp_dense = mean(&sci_dense_kib) / mean(&damp_dense_kib).max(1.0);
        println!(
            "[BVP_sci vs BVP_Damp memory] sparse_vs_damp_dense_ratio={:.3}, sci_dense_vs_damp_dense_ratio={:.3}",
            sci_sparse_vs_damp_dense, sci_dense_vs_damp_dense
        );

        assert!(sci_dense_kib.iter().all(|value| *value > 0.0));
        assert!(damp_dense_kib.iter().all(|value| *value > 0.0));
    }

    /// I.4: Combustion_200 - ExprLegacy stability (5 reps)
    ///
    /// Focused stability check for the ExprLegacy lambdify workflow with
    /// 5 repetitions. Verifies consistent convergence and solution quality.
    /// When AtomView lambdify mode is added to BVP_sci, this test will be
    /// extended to compare ExprLegacy vs AtomView.
    #[test]
    fn combustion_200_exprlegacy_stability_story() {
        let variants: Vec<RaceVariant> = combustion_sparse_variants()
            .into_iter()
            .filter(|v| v.source == "Lambdify")
            .collect();
        let samples = run_race_samples(&variants, 200, 5);
        let summary = summarize_samples(&variants, &samples);
        print_e2e_correctness_table("Combustion 200: ExprLegacy stability (5 reps)", &summary);
        print_e2e_performance_table("Combustion 200: ExprLegacy stability (5 reps)", &summary);
        // Assert: ExprLegacy converges consistently
        for row in &summary {
            assert!(
                row.ok_runs > 0,
                "{} {} {}: all runs failed",
                row.source,
                row.matrix,
                row.variant
            );
        }
    }

    /// I.5: Combustion_200 - AOT lifecycle: BuildIfMissing then RequirePrebuilt
    ///
    /// BVP_sci currently treats RequirePrebuilt as an in-process contract:
    /// a compiled sparse runtime must already be linked in the current process.
    /// This story proves the production handoff sequence that BVP_sci supports
    /// today: BuildIfMissing registers the tcc runtime, then RequirePrebuilt
    /// uses that runtime without falling back to lambdify or rebuilding.
    #[test]
    #[ignore]
    fn combustion_200_tcc_build_then_require_prebuilt_story() {
        let variants = combustion_tcc_lifecycle_variants();
        let samples = run_race_samples(&variants, 200, 3);
        let summary = summarize_samples(&variants, &samples);

        print_e2e_correctness_table(
            "Combustion 200: tcc BuildIfMissing -> RequirePrebuilt correctness",
            &summary,
        );
        print_e2e_performance_table(
            "Combustion 200: tcc BuildIfMissing -> RequirePrebuilt timing",
            &summary,
        );
        print_e2e_lifecycle_table(
            "Combustion 200: tcc BuildIfMissing -> RequirePrebuilt lifecycle",
            &summary,
        );
        print_generated_backend_lifecycle_samples(
            "Combustion 200: tcc BuildIfMissing -> RequirePrebuilt generated backend actions",
            &samples,
        );

        for row in &summary {
            assert!(
                row.ok_runs == row.runs,
                "{} {} {} {}: expected every lifecycle run to pass, got {}",
                row.source,
                row.matrix,
                row.variant,
                row.bootstrap_hint,
                row.status
            );
            if row.source == "AOT" {
                assert!(
                    row.solve_diff.mean < 1e-6,
                    "{} {} {} {}: solve_diff too large: {:e}",
                    row.source,
                    row.matrix,
                    row.variant,
                    row.bootstrap_hint,
                    row.solve_diff.mean
                );
            }
        }

        let aot_samples = samples
            .iter()
            .filter(|sample| sample.row.source == "AOT")
            .collect::<Vec<_>>();
        assert!(
            aot_samples.iter().any(|sample| stats_diagnostic_string(
                &sample.statistics,
                "generated_backend_action"
            ) == "built_and_linked"),
            "BuildIfMissing lifecycle story must include at least one actual build/link event"
        );
        let strict_samples = aot_samples
            .iter()
            .filter(|sample| sample.row.bootstrap_hint == "require_prebuilt_in_process")
            .collect::<Vec<_>>();
        assert!(
            strict_samples.iter().all(|sample| stats_diagnostic_string(
                &sample.statistics,
                "generated_backend_policy"
            ) == "RequirePrebuilt"),
            "RequirePrebuilt lifecycle rows must report the strict policy"
        );
        assert!(
            strict_samples.iter().all(|sample| stats_diagnostic_string(
                &sample.statistics,
                "generated_backend_action"
            ) == "reused_linked"),
            "RequirePrebuilt lifecycle rows must reuse an already linked compiled backend"
        );
    }

    /// I.6: Combustion_200 - Sparse vs promoted AutoBanded vs strict bordered routing
    ///
    /// This is a correctness/diagnostic story, not a performance claim. The
    /// standard endpoint-BC BVP_sci Newton matrix has a compact collocation
    /// body plus boundary rows. AutoBanded must recognize the bordered-banded
    /// candidate and promote this parameter-free endpoint system. The explicit
    /// policy verifies the same native route without fallback.
    #[test]
    fn combustion_200_auto_banded_linear_policy_route_story() {
        let (variants, samples) = run_combustion_linear_policy_samples(200, 1);
        let summary = summarize_samples(&variants, &samples);

        print_e2e_correctness_table(
            "Combustion 200: Sparse vs promoted AutoBanded vs strict bordered correctness",
            &summary,
        );
        print_e2e_performance_table(
            "Combustion 200: Sparse vs promoted AutoBanded vs strict bordered timing",
            &summary,
        );
        print_linear_policy_route_table(
            "Combustion 200: Sparse vs promoted AutoBanded vs strict bordered route counters",
            &samples,
        );

        assert_linear_policy_contract(&summary, &samples);
    }

    /// I.7: Combustion linear-policy release candidate stress.
    ///
    /// This is the performance/correctness companion to the fast combustion-200
    /// route gate. It keeps Sparse as the baseline, verifies AutoBanded
    /// promotion, and measures the strict bordered policy with multi-run
    /// statistics. Defaults are
    /// intentionally moderate; release runs can scale with:
    ///
    /// - `BVP_SCI_LINEAR_POLICY_N_STEPS`
    /// - `BVP_SCI_LINEAR_POLICY_RUNS`
    #[test]
    #[ignore]
    fn combustion_linear_policy_release_candidate_story() {
        let n_steps = env_usize_or("BVP_SCI_LINEAR_POLICY_N_STEPS", 1000);
        let repetitions = env_usize_or("BVP_SCI_LINEAR_POLICY_RUNS", 3);
        println!(
            "[BVP_sci linear policy] release candidate settings: n_steps={n_steps}, runs={repetitions}"
        );

        let (variants, samples) = run_combustion_linear_policy_samples(n_steps, repetitions);
        let summary = summarize_samples(&variants, &samples);

        print_e2e_correctness_table(
            &format!(
                "Combustion {n_steps}: Sparse vs promoted AutoBanded vs strict bordered correctness ({repetitions} runs)"
            ),
            &summary,
        );
        print_e2e_performance_table(
            &format!(
                "Combustion {n_steps}: Sparse vs promoted AutoBanded vs strict bordered timing ({repetitions} runs)"
            ),
            &summary,
        );
        print_linear_policy_route_table(
            &format!(
                "Combustion {n_steps}: Sparse vs promoted AutoBanded vs strict bordered route counters ({repetitions} runs)"
            ),
            &samples,
        );

        assert_linear_policy_contract(&summary, &samples);
    }

    /// I.8: Larger-mesh combustion confirmation for bordered opt-in policy.
    ///
    /// This stress is intentionally separate from the combustion-1000 release
    /// candidate: promotion decisions should not be made from a single mesh size.
    /// Defaults are heavier but still user-configurable:
    ///
    /// - `BVP_SCI_LARGE_LINEAR_POLICY_N_STEPS`
    /// - `BVP_SCI_LARGE_LINEAR_POLICY_RUNS`
    #[test]
    #[ignore]
    fn combustion_large_linear_policy_release_story() {
        let n_steps = env_usize_or("BVP_SCI_LARGE_LINEAR_POLICY_N_STEPS", 3000);
        let repetitions = env_usize_or("BVP_SCI_LARGE_LINEAR_POLICY_RUNS", 3);
        println!(
            "[BVP_sci linear policy] large combustion settings: n_steps={n_steps}, runs={repetitions}"
        );

        let (variants, samples) = run_combustion_linear_policy_samples(n_steps, repetitions);
        let summary = summarize_samples(&variants, &samples);

        print_e2e_correctness_table(
            &format!(
                "Combustion {n_steps}: large Sparse vs promoted AutoBanded vs strict bordered correctness ({repetitions} runs)"
            ),
            &summary,
        );
        print_e2e_performance_table(
            &format!(
                "Combustion {n_steps}: large Sparse vs promoted AutoBanded vs strict bordered timing ({repetitions} runs)"
            ),
            &summary,
        );
        print_jacobian_memory_table(
            &format!(
                "Combustion {n_steps}: large Sparse vs promoted AutoBanded vs strict bordered Jacobian memory ({repetitions} runs)"
            ),
            &samples,
        );
        print_linear_policy_route_table(
            &format!(
                "Combustion {n_steps}: large Sparse vs promoted AutoBanded vs strict bordered route counters ({repetitions} runs)"
            ),
            &samples,
        );

        assert_linear_policy_contract(&summary, &samples);
    }

    /// I.9: Non-combustion endpoint-BC confirmation for bordered opt-in policy.
    ///
    /// Uses the exponential BVP already present in the BVP_sci compare suite.
    /// This guards against accidentally tuning the bordered route only for the
    /// combustion matrix shape.
    ///
    /// - `BVP_SCI_NONCOMB_LINEAR_POLICY_N_STEPS`
    /// - `BVP_SCI_NONCOMB_LINEAR_POLICY_RUNS`
    #[test]
    #[ignore]
    fn exponential_endpoint_linear_policy_release_story() {
        let n_steps = env_usize_or("BVP_SCI_NONCOMB_LINEAR_POLICY_N_STEPS", 1000);
        let repetitions = env_usize_or("BVP_SCI_NONCOMB_LINEAR_POLICY_RUNS", 3);
        println!(
            "[BVP_sci linear policy] exponential endpoint settings: n_steps={n_steps}, runs={repetitions}"
        );

        let (variants, samples) = run_linear_policy_samples_with_options(
            || exponential_endpoint_problem_options(n_steps),
            repetitions,
        );
        let summary = summarize_samples(&variants, &samples);

        print_e2e_correctness_table(
            &format!(
                "Exponential endpoint {n_steps}: Sparse vs promoted AutoBanded vs strict bordered correctness ({repetitions} runs)"
            ),
            &summary,
        );
        print_e2e_performance_table(
            &format!(
                "Exponential endpoint {n_steps}: Sparse vs promoted AutoBanded vs strict bordered timing ({repetitions} runs)"
            ),
            &summary,
        );
        print_linear_policy_route_table(
            &format!(
                "Exponential endpoint {n_steps}: Sparse vs promoted AutoBanded vs strict bordered route counters ({repetitions} runs)"
            ),
            &samples,
        );

        assert_linear_policy_contract(&summary, &samples);
    }

    // ============================================================
    // Section J: Process-isolated cold test infrastructure
    // ============================================================

    const ISOLATED_STRESS_CHILD_INDEX_ENV: &str = "BVP_SCI_ISOLATED_STRESS_CHILD_INDEX";
    const ISOLATED_STRESS_CHILD_REPETITION_ENV: &str = "BVP_SCI_ISOLATED_STRESS_CHILD_REPETITION";
    const ISOLATED_RACE_ROW_MARKER: &str = "[BVP_SCI_ISOLATED_RACE_ROW]";
    const ISOLATED_RACE_SOLUTION_MARKER: &str = "[BVP_SCI_ISOLATED_RACE_SOLUTION]";
    const ISOLATED_RACE_PID_MARKER: &str = "[BVP_SCI_ISOLATED_RACE_PID]";
    const COLD_COOLDOWN_MS_ENV: &str = "BVP_SCI_COLD_COOLDOWN_MS";
    const COLD_CLEAN_ARTIFACTS_ENV: &str = "BVP_SCI_COLD_CLEAN_ARTIFACTS";
    const COMBUSTION_3000_ISOLATED_STRESS_TEST_NAME: &str = "numerical::BVP_sci::legacy::BVP_sci_story_tests::tests::combustion_3000_sparse_isolated_stress_story";

    /// Encode a `RaceRow` as a tab-separated string for child→parent IPC.
    /// Only the 9 hardcoded fields are serialized; solver-specific metrics
    /// are not available in the child process (no `BvpSciStatistics` cross-process).
    fn encode_isolated_race_row(row: &RaceRow) -> String {
        [
            row.total_ms.to_string(),
            row.max_abs_solution.to_string(),
            row.status.clone(),
        ]
        .join("\t")
    }

    /// Parse a single tab-separated field from an iterator.
    fn parse_isolated_field<T: std::str::FromStr>(
        fields: &mut impl Iterator<Item = String>,
        name: &str,
    ) -> T
    where
        T::Err: std::fmt::Debug,
    {
        fields
            .next()
            .unwrap_or_else(|| panic!("isolated race row missing field {name}"))
            .parse::<T>()
            .unwrap_or_else(|err| panic!("isolated race row field {name} could not parse: {err:?}"))
    }

    /// Decode a `RaceRow` from a child process stdout line.
    /// The `variant` provides the static fields (source, matrix, variant, bootstrap_hint).
    fn decode_isolated_race_row(line: &str, variant: &RaceVariant) -> RaceRow {
        let payload = line
            .strip_prefix(ISOLATED_RACE_ROW_MARKER)
            .expect("isolated race row marker should be present")
            .trim_start_matches('\t');
        let mut fields = payload.split('\t').map(str::to_string);
        let row = RaceRow {
            source: variant.source,
            matrix: variant.matrix,
            variant: variant.variant,
            bootstrap_hint: variant.bootstrap_hint,
            total_ms: parse_isolated_field(&mut fields, "total_ms"),
            max_abs_solution: parse_isolated_field(&mut fields, "max_abs_solution"),
            solve_diff: 0.0,
            rel_x_diff: 0.0,
            status: parse_isolated_field(&mut fields, "status"),
        };
        assert!(
            fields.next().is_none(),
            "isolated race row contains unexpected trailing fields"
        );
        row
    }

    /// Encode a `DMatrix<f64>` solution as a tab-separated string for child→parent IPC.
    fn encode_isolated_solution(solution: &DMatrix<f64>) -> String {
        solution
            .iter()
            .map(|value| value.to_string())
            .collect::<Vec<_>>()
            .join("\t")
    }

    /// Decode a `DMatrix<f64>` solution from a child process stdout line.
    /// The combustion problem has 6 variables, so we assert `values.len() % 6 == 0`.
    fn decode_isolated_solution(line: &str) -> DMatrix<f64> {
        let payload = line
            .strip_prefix(ISOLATED_RACE_SOLUTION_MARKER)
            .expect("isolated solution marker should be present")
            .trim_start_matches('\t');
        let values = payload
            .split('\t')
            .filter(|value| !value.is_empty())
            .map(|value| {
                value
                    .parse::<f64>()
                    .unwrap_or_else(|err| panic!("isolated solution value could not parse: {err}"))
            })
            .collect::<Vec<_>>();
        assert_eq!(
            values.len() % 6,
            0,
            "isolated combustion solution should contain six state rows"
        );
        DMatrix::from_column_slice(6, values.len() / 6, &values)
    }

    fn cold_cooldown_ms() -> u64 {
        std::env::var(COLD_COOLDOWN_MS_ENV)
            .ok()
            .and_then(|value| value.parse::<u64>().ok())
            .unwrap_or(0)
    }

    fn clean_cold_artifacts_enabled() -> bool {
        std::env::var(COLD_CLEAN_ARTIFACTS_ENV)
            .map(|value| matches!(value.to_ascii_lowercase().as_str(), "1" | "true" | "yes"))
            .unwrap_or(false)
    }

    fn remove_generated_aot_builds_for_child(child_pid: u32) {
        if !clean_cold_artifacts_enabled() {
            return;
        }
        let root = PathBuf::from(env!("CARGO_MANIFEST_DIR"))
            .join("target")
            .join("generated-aot");
        let prefix = format!("build-{child_pid}-");
        let Ok(problem_dirs) = fs::read_dir(&root) else {
            return;
        };
        for problem_dir in problem_dirs.flatten() {
            let Ok(build_dirs) = fs::read_dir(problem_dir.path()) else {
                continue;
            };
            for build_dir in build_dirs.flatten() {
                let name = build_dir.file_name().to_string_lossy().to_string();
                if name.starts_with(&prefix) {
                    fs::remove_dir_all(build_dir.path()).unwrap_or_else(|err| {
                        panic!(
                            "failed to remove isolated child AOT artifact directory {}: {err}",
                            build_dir.path().display()
                        )
                    });
                }
            }
        }
    }

    /// Run process-isolated cold samples by spawning child processes.
    ///
    /// Each child process runs a single variant (identified by
    /// `ISOLATED_STRESS_CHILD_INDEX_ENV`) and emits the race row + solution
    /// as marker-prefixed lines on stdout. The parent collects these,
    /// decodes them, computes solution diffs, and returns `Vec<RaceRow>`.
    ///
    /// Note: Unlike `run_race_samples`, this returns `Vec<RaceRow>` (not
    /// `Vec<RaceSample>`) because `BvpSciStatistics` cannot be serialized
    /// across process boundaries. The caller can still use `summarize_samples`
    /// by wrapping rows in a dummy `RaceSample` with default statistics.
    fn run_isolated_race_samples(
        test_name: &str,
        variants: &[RaceVariant],
        _n_steps: usize,
        repetitions: usize,
    ) -> Vec<RaceRow> {
        let executable = std::env::current_exe().expect("current test executable should resolve");
        let mut samples = Vec::with_capacity(variants.len() * repetitions);
        for repetition in 0..repetitions {
            let mut rows = Vec::with_capacity(variants.len());
            let mut solutions = Vec::with_capacity(variants.len());
            for (index, variant) in variants.iter().enumerate() {
                println!(
                    "[BVP_sci isolated cold] launching repetition {}/{} source={} variant={}",
                    repetition + 1,
                    repetitions,
                    variant.source,
                    variant.variant
                );
                let output = Command::new(&executable)
                    .arg("--exact")
                    .arg(test_name)
                    .arg("--ignored")
                    .arg("--nocapture")
                    .env(ISOLATED_STRESS_CHILD_INDEX_ENV, index.to_string())
                    .env(ISOLATED_STRESS_CHILD_REPETITION_ENV, repetition.to_string())
                    .output()
                    .expect("isolated cold child process should launch");
                let stdout = String::from_utf8_lossy(&output.stdout);
                let stderr = String::from_utf8_lossy(&output.stderr);
                assert!(
                    output.status.success(),
                    "isolated cold child failed for {}:\nstdout:\n{}\nstderr:\n{}",
                    variant.variant,
                    stdout,
                    stderr
                );
                let row_line = stdout
                    .lines()
                    .find(|line| line.starts_with(ISOLATED_RACE_ROW_MARKER))
                    .unwrap_or_else(|| {
                        panic!(
                            "isolated cold child did not emit metrics for {} using exact test name `{}`:\nstdout:\n{}\nstderr:\n{}",
                            variant.variant, test_name, stdout, stderr
                        )
                    });
                let solution_line = stdout
                    .lines()
                    .find(|line| line.starts_with(ISOLATED_RACE_SOLUTION_MARKER))
                    .unwrap_or_else(|| {
                        panic!(
                            "isolated cold child did not emit a solution for {}:\n{}",
                            variant.variant, stdout
                        )
                    });
                let row = decode_isolated_race_row(row_line, variant);
                let solution = decode_isolated_solution(solution_line);
                let child_pid = stdout
                    .lines()
                    .find(|line| line.starts_with(ISOLATED_RACE_PID_MARKER))
                    .and_then(|line| line.strip_prefix(ISOLATED_RACE_PID_MARKER).map(str::trim))
                    .and_then(|value| value.parse::<u32>().ok())
                    .unwrap_or_else(|| {
                        panic!("isolated cold child emitted no pid for {}", variant.variant)
                    });
                remove_generated_aot_builds_for_child(child_pid);
                let cooldown_ms = cold_cooldown_ms();
                if cooldown_ms > 0 {
                    thread::sleep(Duration::from_millis(cooldown_ms));
                }
                println!(
                    "[BVP_sci isolated cold] finished source={} variant={} total_ms={:.3} status={}",
                    row.source, row.variant, row.total_ms, row.status
                );
                rows.push(row);
                solutions.push(Some(solution));
            }
            fill_solution_diffs_from_rows(&mut rows, &solutions);
            samples.extend(rows);
        }
        samples
    }

    /// Fill solution diffs for `Vec<RaceRow>` (no `BvpSciStatistics` wrapper).
    fn fill_solution_diffs_from_rows(rows: &mut [RaceRow], solutions: &[Option<DMatrix<f64>>]) {
        let baseline = solutions.iter().find_map(|solution| solution.as_ref());
        if let Some(baseline) = baseline {
            for (row, solution) in rows.iter_mut().zip(solutions.iter()) {
                if let Some(solution) = solution {
                    row.solve_diff = solution_linf_diff(solution, baseline);
                    row.rel_x_diff = solution_rel_diff(solution, baseline);
                }
            }
        }
    }

    /// Print a raw isolated cold sample table (one row per observation).
    fn print_isolated_cold_sample_table(title: &str, samples: &[RaceRow], variants_count: usize) {
        println!("\n{}", title);
        println!(
            "source    | matrix | variant        | bootstrap_hint | total_ms     | max_abs_solution | solve_diff    | rel_x_diff    | status"
        );
        println!(
            "----------|--------|----------------|----------------|--------------|------------------|---------------|---------------|--------"
        );
        for (i, row) in samples.iter().enumerate() {
            let rep = i / variants_count;
            println!(
                "{} | {} | {} | {} | {:>12.3} | {:>16.6e} | {:>13.6e} | {:>13.6e} | {}",
                row.source,
                row.matrix,
                row.variant,
                row.bootstrap_hint,
                row.total_ms,
                row.max_abs_solution,
                row.solve_diff,
                row.rel_x_diff,
                row.status
            );
            if (i + 1) % variants_count == 0 && (i + 1) < samples.len() {
                println!("--- repetition {} ---", rep);
            }
        }
    }

    /// I.7: Combustion_1000 - banded production-ready route story
    ///
    /// This is the first canonical banded story for `BVP_sci`.  It keeps the
    /// already-safe `AutoBanded` policy, the explicit experimental bordered
    /// route, and the Sparse baseline in one multi-run matrix so we can track:
    /// - correctness
    /// - total wall-clock timing
    /// - stage timing breakdown
    /// - route counters and memory footprint
    ///
    /// The goal is to make banded a first-class production candidate instead of
    /// a hidden internal branch.
    #[test]
    #[ignore]
    fn combustion_1000_banded_production_story() {
        let (variants, samples) = run_combustion_linear_policy_samples(1_000, 5);
        let summary = summarize_samples(&variants, &samples);

        print_e2e_correctness_table(
            "Combustion 1000: Sparse vs promoted AutoBanded vs strict bordered correctness (5 runs)",
            &summary,
        );
        print_e2e_performance_table(
            "Combustion 1000: Sparse vs promoted AutoBanded vs strict bordered timing (5 runs)",
            &summary,
        );
        print_e2e_stage_breakdown_table(
            "Combustion 1000: Sparse vs promoted AutoBanded vs strict bordered stage breakdown (5 runs)",
            &samples,
        );
        print_linear_policy_route_table(
            "Combustion 1000: Sparse vs promoted AutoBanded vs strict bordered route counters (5 runs)",
            &samples,
        );
        print_jacobian_memory_table(
            "Combustion 1000: Sparse vs promoted AutoBanded vs strict bordered Jacobian memory (5 runs)",
            &samples,
        );

        assert_linear_policy_contract(&summary, &samples);
    }

    /// I.8: Combustion_3000 - larger-mesh banded stress story
    ///
    /// This is the same banded-vs-sparse contract, but on the larger combustion
    /// grid so we can see whether the banded route still behaves sanely at a
    /// heavier scale.
    #[test]
    #[ignore]
    fn combustion_3000_banded_production_story() {
        let (variants, samples) = run_combustion_linear_policy_samples(3_000, 3);
        let summary = summarize_samples(&variants, &samples);

        print_e2e_correctness_table(
            "Combustion 3000: Sparse vs promoted AutoBanded vs strict bordered correctness (3 runs)",
            &summary,
        );
        print_e2e_performance_table(
            "Combustion 3000: Sparse vs promoted AutoBanded vs strict bordered timing (3 runs)",
            &summary,
        );
        print_e2e_stage_breakdown_table(
            "Combustion 3000: Sparse vs promoted AutoBanded vs strict bordered stage breakdown (3 runs)",
            &samples,
        );
        print_linear_policy_route_table(
            "Combustion 3000: Sparse vs promoted AutoBanded vs strict bordered route counters (3 runs)",
            &samples,
        );
        print_jacobian_memory_table(
            "Combustion 3000: Sparse vs promoted AutoBanded vs strict bordered Jacobian memory (3 runs)",
            &samples,
        );

        assert_linear_policy_contract(&summary, &samples);
    }

    // ============================================================
    // Section K: Isolated stress story test
    // ============================================================

    /// J.1: Combustion_3000 — Process-isolated cold stress test
    ///
    /// Runs the combustion problem with 3000 mesh points in separate child
    /// processes to measure cold-start performance (no warm caches, no
    /// reused AOT artifacts from previous in-process runs). Each variant
    /// is run 2 times in isolated processes.
    ///
    /// This test is `#[ignore]` by default because it spawns child processes
    /// and takes significant time. Run with `cargo test -- --ignored`.
    #[test]
    #[ignore]
    fn combustion_3000_sparse_isolated_stress_story() {
        let n_steps = 3_000usize;
        let repetitions = 2usize;
        let variants = combustion_sparse_variants();

        // Child process path: run a single variant and emit metrics
        if let Ok(index) = std::env::var(ISOLATED_STRESS_CHILD_INDEX_ENV) {
            let index = index
                .parse::<usize>()
                .expect("isolated stress child index should be an integer");
            let variant = variants.get(index).unwrap_or_else(|| {
                panic!("isolated stress child variant index {index} is invalid")
            });
            let (sample, solution) = run_race_variant(
                n_steps,
                variant.source,
                variant.matrix,
                variant.variant,
                variant.bootstrap_hint,
                variant.config.clone(),
            );
            assert_eq!(
                sample.row.status, "ok",
                "isolated stress child {} should solve successfully",
                variant.variant
            );
            let solution = solution
                .as_ref()
                .expect("isolated stress child should provide its converged solution");
            println!("{ISOLATED_RACE_PID_MARKER}\t{}", std::process::id());
            println!(
                "{ISOLATED_RACE_ROW_MARKER}\t{}",
                encode_isolated_race_row(&sample.row)
            );
            println!(
                "{ISOLATED_RACE_SOLUTION_MARKER}\t{}",
                encode_isolated_solution(solution)
            );
            return;
        }

        // Parent process path: spawn children and collect results
        println!(
            "[BVP_sci isolated cold] protocol cooldown_ms={}, cleanup_child_artifacts={}",
            cold_cooldown_ms(),
            clean_cold_artifacts_enabled()
        );
        let samples = run_isolated_race_samples(
            COMBUSTION_3000_ISOLATED_STRESS_TEST_NAME,
            &variants,
            n_steps,
            repetitions,
        );

        // Convert RaceRows to RaceSamples with default statistics for summarization
        let default_stats = BvpSciStatistics::default();
        let race_samples: Vec<RaceSample> = samples
            .iter()
            .map(|row| RaceSample {
                row: row.clone(),
                statistics: default_stats.clone(),
            })
            .collect();
        let summary = summarize_samples(&variants, &race_samples);

        println!();
        print_isolated_cold_sample_table(
            "[BVP_sci stress] combustion-3000 raw process-isolated cold observations",
            &samples,
            variants.len(),
        );
        println!();
        print_e2e_correctness_table(
            "[BVP_sci stress] combustion-3000 Lambdify vs AOT correctness (isolated cold)",
            &summary,
        );
        println!();
        print_e2e_performance_table(
            "[BVP_sci stress] combustion-3000 Lambdify vs AOT timing (isolated cold)",
            &summary,
        );

        // Assert: all variants converge in isolated cold runs
        for row in &summary {
            assert!(
                row.ok_runs > 0,
                "{} {} {}: all isolated cold runs failed",
                row.source,
                row.matrix,
                row.variant
            );
        }
    }
}