//! Integration tests for [`cobre_sddp::simulate`] (simulation pipeline).
//!
//! Uses a [`MockSolver`] and [`StubComm`] to exercise the simulation pipeline
//! end-to-end without a real LP solver or MPI communicator. Covers scenario
//! count, error propagation, cost accumulation, event emission, load patching,
//! inflow truncation, baked-template acceptance, and warm-start basis handling.

#![allow(
    clippy::unwrap_used,
    clippy::expect_used,
    clippy::panic,
    clippy::too_many_lines,
    clippy::cast_possible_truncation,
    clippy::cast_possible_wrap,
    clippy::cast_sign_loss,
    clippy::float_cmp,
    clippy::cast_precision_loss
)]

use std::collections::HashMap;
use std::sync::mpsc;

use cobre_comm::{CommData, CommError, Communicator, ReduceOp};
use cobre_core::scenario::SamplingScheme;
use cobre_solver::{
    Basis, LpSolution, RowBatch, SolverError, SolverInterface, SolverStatistics, StageTemplate,
};
use cobre_stochastic::StochasticContext;

use cobre_sddp::{
    CapturedBasis, SimulationError,
    context::{StageContext, TrainingContext},
    cut::FutureCostFunction,
    horizon_mode::HorizonMode,
    indexer::StageIndexer,
    inflow_method::InflowNonNegativityMethod,
    lp_builder::PatchBuffer,
    simulation::{EntityCounts, SimulationConfig, SimulationOutputSpec},
    workspace::{SolverWorkspace, WorkspaceSizing},
};

// ── Stub communicator ────────────────────────────────────────────────────────

/// Single-rank stub communicator for pipeline tests.
struct StubComm {
    rank: usize,
    size: usize,
}

impl Communicator for StubComm {
    fn allgatherv<T: CommData>(
        &self,
        _send: &[T],
        _recv: &mut [T],
        _counts: &[usize],
        _displs: &[usize],
    ) -> Result<(), CommError> {
        unreachable!("StubComm allgatherv not used in simulate tests")
    }

    fn allreduce<T: CommData>(
        &self,
        _send: &[T],
        _recv: &mut [T],
        _op: ReduceOp,
    ) -> Result<(), CommError> {
        unreachable!("StubComm allreduce not used in simulate tests")
    }

    fn broadcast<T: CommData>(&self, _buf: &mut [T], _root: usize) -> Result<(), CommError> {
        unreachable!("StubComm broadcast not used in simulate tests")
    }

    fn barrier(&self) -> Result<(), CommError> {
        Ok(())
    }

    fn rank(&self) -> usize {
        self.rank
    }

    fn size(&self) -> usize {
        self.size
    }

    fn abort(&self, error_code: i32) -> ! {
        std::process::exit(error_code)
    }
}

// ── Mock solver ──────────────────────────────────────────────────────────

/// Mock solver that returns a configurable fixed `LpSolution` on every solve.
///
/// Optionally returns `SolverError::Infeasible` at a specific (0-based) solve
/// call index (counting across both cold-start and warm-start calls).
///
/// `load_count` and `add_rows_count` track how many times `load_model` and
/// `add_rows` were called, used by the baked-path acceptance tests.
///
/// `solve_count` counts calls to the cold-start `solve(None)` path only.
/// `solve_with_basis_count` counts calls to the warm-start `solve(Some(&basis))`
/// path only.  `call_count` tracks the combined total across both paths and
/// is used by the infeasibility injection logic.
///
/// `recorded_basis` captures the last `Basis` passed to `solve(Some(&basis))`,
/// used by the warm-start reconstruction acceptance tests.
///
/// `reconstruction_counter` counts `record_reconstruction_stats` invocations
/// (one per warm-start solve that applied a stored basis via slot reconciliation).
struct MockSolver {
    solution: LpSolution,
    infeasible_at: Option<usize>,
    call_count: usize,
    buf_primal: Vec<f64>,
    buf_dual: Vec<f64>,
    buf_reduced_costs: Vec<f64>,
    /// Number of `load_model` calls since construction.
    load_count: usize,
    /// Number of `add_rows` calls since construction.
    add_rows_count: usize,
    /// Number of cold-start `solve(None)` calls.
    solve_count: usize,
    /// Number of warm-start `solve(Some(&basis))` calls.
    solve_with_basis_count: usize,
    /// Last `Basis` passed to `solve(Some(&basis))`, if any.
    recorded_basis: Option<Basis>,
    /// Cumulative `preserved` argument from `record_reconstruction_stats`.
    reconstruction_counter: u32,
}

impl MockSolver {
    fn always_ok(solution: LpSolution) -> Self {
        let buf_primal = solution.primal.clone();
        let buf_dual = solution.dual.clone();
        let buf_reduced_costs = solution.reduced_costs.clone();
        Self {
            solution,
            infeasible_at: None,
            call_count: 0,
            buf_primal,
            buf_dual,
            buf_reduced_costs,
            load_count: 0,
            add_rows_count: 0,
            solve_count: 0,
            solve_with_basis_count: 0,
            recorded_basis: None,
            reconstruction_counter: 0,
        }
    }

    fn infeasible_on(solution: LpSolution, n: usize) -> Self {
        let buf_primal = solution.primal.clone();
        let buf_dual = solution.dual.clone();
        let buf_reduced_costs = solution.reduced_costs.clone();
        Self {
            solution,
            infeasible_at: Some(n),
            call_count: 0,
            buf_primal,
            buf_dual,
            buf_reduced_costs,
            load_count: 0,
            add_rows_count: 0,
            solve_count: 0,
            solve_with_basis_count: 0,
            recorded_basis: None,
            reconstruction_counter: 0,
        }
    }

    fn do_solve(&mut self) -> Result<cobre_solver::SolutionView<'_>, SolverError> {
        let call = self.call_count;
        self.call_count += 1;
        if self.infeasible_at == Some(call) {
            return Err(SolverError::Infeasible);
        }
        self.buf_primal.clone_from(&self.solution.primal);
        self.buf_dual.clone_from(&self.solution.dual);
        self.buf_reduced_costs
            .clone_from(&self.solution.reduced_costs);
        Ok(cobre_solver::SolutionView {
            objective: self.solution.objective,
            primal: &self.buf_primal,
            dual: &self.buf_dual,
            reduced_costs: &self.buf_reduced_costs,
            iterations: self.solution.iterations,
            solve_time_seconds: self.solution.solve_time_seconds,
        })
    }
}

impl SolverInterface for MockSolver {
    fn solver_name_version(&self) -> String {
        "MockSolver 0.0.0".to_string()
    }
    fn load_model(&mut self, _template: &StageTemplate) {
        self.load_count += 1;
    }
    fn add_rows(&mut self, _cuts: &RowBatch) {
        self.add_rows_count += 1;
    }
    fn set_row_bounds(&mut self, _indices: &[usize], _lower: &[f64], _upper: &[f64]) {}
    fn set_col_bounds(&mut self, _indices: &[usize], _lower: &[f64], _upper: &[f64]) {}
    fn solve(
        &mut self,
        basis: Option<&Basis>,
    ) -> Result<cobre_solver::SolutionView<'_>, SolverError> {
        if let Some(b) = basis {
            self.solve_with_basis_count += 1;
            self.recorded_basis = Some(b.clone());
        } else {
            self.solve_count += 1;
        }
        self.do_solve()
    }
    fn get_basis(&mut self, _out: &mut Basis) {}
    fn record_reconstruction_stats(&mut self) {
        self.reconstruction_counter += 1;
    }
    fn statistics(&self) -> SolverStatistics {
        SolverStatistics::default()
    }
    fn name(&self) -> &'static str {
        "Mock"
    }
}

// ── Helpers ──────────────────────────────────────────────────────────────

/// Minimal valid stage template for N=1 hydro, L=0 PAR order.
///
/// Column layout (N=1, L=0):
/// - col 0: `storage_out` (no NZ in structural rows)
/// - col 1: `z_inflow` (no NZ — `z_inflow` row at row 1)
/// - col 2: `storage_in` (1 NZ: row 0, storage-fixing row)
/// - col 3: `theta` (no NZ)
///
/// Row layout:
/// - row 0: storage-fixing (`storage_out` fixed to incoming state)
/// - row 1: `z_inflow` definition row
fn minimal_template_1_0() -> StageTemplate {
    StageTemplate {
        num_cols: 4,
        num_rows: 2,
        num_nz: 1,
        // CSC col_starts: 4 cols + 1 sentinel = 5 entries.
        // col 0 (storage_out): 0 NZ
        // col 1 (z_inflow):    0 NZ
        // col 2 (storage_in):  1 NZ at row 0
        // col 3 (theta):       0 NZ
        col_starts: vec![0_i32, 0, 0, 1, 1],
        row_indices: vec![0_i32],
        values: vec![1.0],
        col_lower: vec![0.0, f64::NEG_INFINITY, 0.0, 0.0],
        col_upper: vec![f64::INFINITY; 4],
        objective: vec![0.0, 0.0, 0.0, 1.0],
        row_lower: vec![0.0, 0.0],
        row_upper: vec![0.0, 0.0],
        n_state: 1,
        n_transfer: 0,
        n_dual_relevant: 1,
        n_hydro: 1,
        max_par_order: 0,
        col_scale: Vec::new(),
        row_scale: Vec::new(),
    }
}

/// Build a fixed `LpSolution` for the minimal N=1 L=0 template.
///
/// N=1 L=0 column layout: `storage`(0), `z_inflow`(1), `storage_in`(2), `theta`(3).
/// `primal[3] = theta_val`, `objective = objective`.
fn fixed_solution(objective: f64, theta_val: f64) -> LpSolution {
    let num_cols = 4; // storage(0), z_inflow(1), storage_in(2), theta(3)
    let mut primal = vec![0.0_f64; num_cols];
    primal[3] = theta_val; // theta at col 3 (N=1, L=0 → theta = N*(3+L) = 3)
    LpSolution {
        objective,
        primal,
        dual: vec![0.0_f64; 2],
        reduced_costs: vec![0.0_f64; num_cols],
        iterations: 0,
        solve_time_seconds: 0.0,
    }
}

/// Build a minimal `EntityCounts` for 1 hydro, no other entities.
fn entity_counts_1_hydro() -> EntityCounts {
    EntityCounts {
        hydro_ids: vec![1],
        hydro_productivities: vec![1.0],
        thermal_ids: vec![],
        line_ids: vec![],
        bus_ids: vec![],
        pumping_station_ids: vec![],
        contract_ids: vec![],
        non_controllable_ids: vec![],
    }
}

/// Build a minimal stochastic context for 1 hydro, `n_stages` stages.
fn make_stochastic_context(n_stages: usize) -> StochasticContext {
    use std::collections::BTreeMap;

    use chrono::NaiveDate;
    use cobre_core::entities::hydro::{Hydro, HydroGenerationModel, HydroPenalties};
    use cobre_core::scenario::{
        CorrelationEntity, CorrelationGroup, CorrelationModel, CorrelationProfile, InflowModel,
    };
    use cobre_core::temporal::{
        Block, BlockMode, NoiseMethod, ScenarioSourceConfig, Stage, StageRiskConfig,
        StageStateConfig,
    };
    use cobre_core::{Bus, DeficitSegment, EntityId, SystemBuilder};
    use cobre_stochastic::context::{ClassSchemes, OpeningTreeInputs, build_stochastic_context};

    let bus = Bus {
        id: EntityId(0),
        name: "B0".to_string(),
        deficit_segments: vec![DeficitSegment {
            depth_mw: None,
            cost_per_mwh: 1000.0,
        }],
        excess_cost: 0.0,
    };
    let hydro = Hydro {
        id: EntityId(1),
        name: "H1".to_string(),
        bus_id: EntityId(0),
        downstream_id: None,
        entry_stage_id: None,
        exit_stage_id: None,
        min_storage_hm3: 0.0,
        max_storage_hm3: 100.0,
        min_outflow_m3s: 0.0,
        max_outflow_m3s: None,
        generation_model: HydroGenerationModel::ConstantProductivity,
        min_turbined_m3s: 0.0,
        max_turbined_m3s: 100.0,
        specific_productivity_mw_per_m3s_per_m: None,
        min_generation_mw: 0.0,
        max_generation_mw: 100.0,
        tailrace: None,
        hydraulic_losses: None,
        efficiency: None,
        evaporation_coefficients_mm: None,
        evaporation_reference_volumes_hm3: None,
        diversion: None,
        filling: None,
        penalties: HydroPenalties {
            spillage_cost: 0.0,
            diversion_cost: 0.0,
            fpha_turbined_cost: 0.0,
            storage_violation_below_cost: 0.0,
            filling_target_violation_cost: 0.0,
            turbined_violation_below_cost: 0.0,
            outflow_violation_below_cost: 0.0,
            outflow_violation_above_cost: 0.0,
            generation_violation_below_cost: 0.0,
            evaporation_violation_cost: 0.0,
            water_withdrawal_violation_cost: 0.0,
            water_withdrawal_violation_pos_cost: 0.0,
            water_withdrawal_violation_neg_cost: 0.0,
            evaporation_violation_pos_cost: 0.0,
            evaporation_violation_neg_cost: 0.0,
            inflow_nonnegativity_cost: 1000.0,
        },
    };
    let make_stage = |idx: usize, id: i32| Stage {
        index: idx,
        id,
        start_date: NaiveDate::from_ymd_opt(2024, 1, 1).unwrap(),
        end_date: NaiveDate::from_ymd_opt(2024, 2, 1).unwrap(),
        season_id: Some(0),
        blocks: vec![Block {
            index: 0,
            name: "S".to_string(),
            duration_hours: 744.0,
        }],
        block_mode: BlockMode::Parallel,
        state_config: StageStateConfig {
            storage: true,
            inflow_lags: false,
        },
        risk_config: StageRiskConfig::Expectation,
        scenario_config: ScenarioSourceConfig {
            branching_factor: 3,
            noise_method: NoiseMethod::Saa,
        },
    };
    let stages: Vec<Stage> = (0..n_stages)
        .map(|i| make_stage(i, i32::try_from(i).unwrap()))
        .collect();
    let inflow = |stage_id: i32| InflowModel {
        hydro_id: EntityId(1),
        stage_id,
        mean_m3s: 100.0,
        std_m3s: 30.0,
        ar_coefficients: vec![],
        residual_std_ratio: 1.0,
        annual: None,
    };
    let inflow_models: Vec<InflowModel> = (0..n_stages)
        .map(|i| inflow(i32::try_from(i).unwrap()))
        .collect();
    let mut profiles = BTreeMap::new();
    profiles.insert(
        "default".to_string(),
        CorrelationProfile {
            groups: vec![CorrelationGroup {
                name: "g1".to_string(),
                entities: vec![CorrelationEntity {
                    entity_type: "inflow".to_string(),
                    id: EntityId(1),
                }],
                matrix: vec![vec![1.0]],
            }],
        },
    );
    let correlation = CorrelationModel {
        method: "spectral".to_string(),
        profiles,
        schedule: vec![],
    };
    let system = SystemBuilder::new()
        .buses(vec![bus])
        .hydros(vec![hydro])
        .stages(stages)
        .inflow_models(inflow_models)
        .correlation(correlation)
        .build()
        .unwrap();
    build_stochastic_context(
        &system,
        42,
        None,
        &[],
        &[],
        OpeningTreeInputs::default(),
        ClassSchemes {
            inflow: Some(SamplingScheme::InSample),
            load: Some(SamplingScheme::InSample),
            ncs: Some(SamplingScheme::InSample),
        },
    )
    .unwrap()
}

// ── Helpers ──────────────────────────────────────────────────────────────

/// Build per-stage hydro productivities for 1 hydro with productivity 1.0.
///
/// Returns `n_stages` inner vecs each containing a single `1.0` entry,
/// matching `entity_counts_1_hydro().hydro_productivities`.
fn hydro_productivities_1hydro(n_stages: usize) -> Vec<Vec<f64>> {
    vec![vec![1.0]; n_stages]
}

/// Build a zero-valued [`cobre_sddp::energy_conversion::EnergyConversionSet`] for tests
/// that do not assert on energy fields.
fn zero_energy_conversion(
    n_hydros: usize,
    n_stages: usize,
) -> cobre_sddp::energy_conversion::EnergyConversionSet {
    use cobre_sddp::energy_conversion::{EnergyConversion, EnergyConversionSet};
    let zero_ec = EnergyConversion {
        equivalent_productivity_mw_per_m3s: 0.0,
        reference_volume_hm3: 0.0,
        reference_outflow_m3s: 0.0,
    };
    EnergyConversionSet::new(
        vec![vec![zero_ec; n_stages]; n_hydros],
        vec![vec![0.0_f64; n_stages]; n_hydros],
        n_hydros,
        n_stages,
    )
}

/// Wrap a `MockSolver` in a single-workspace slice for `simulate()` calls.
///
/// All tests use a single workspace (serial execution) so that existing
/// assertions about scenario ordering and call counts remain valid.
fn single_workspace(solver: MockSolver) -> Vec<SolverWorkspace<MockSolver>> {
    vec![SolverWorkspace::new(
        0,
        0,
        solver,
        PatchBuffer::new(1, 0, 0, 0),
        1,
        WorkspaceSizing {
            hydro_count: 1,
            ..WorkspaceSizing::default()
        },
    )]
}

// ── Tests ────────────────────────────────────────────────────────────────

/// Acceptance criterion: `n_scenarios=4`, single rank → exactly 4 results in
/// channel and cost buffer has length 4.
#[test]
fn simulate_single_rank_4_scenarios_produces_4_results() {
    let n_stages = 2;
    let templates: Vec<StageTemplate> = (0..n_stages).map(|_| minimal_template_1_0()).collect();
    let base_rows: Vec<usize> = vec![0; n_stages];

    let indexer = StageIndexer::new(1, 0); // N=1, L=0; theta=2
    let fcf = FutureCostFunction::new(n_stages, 1, 1, 10, &vec![0; n_stages]);
    let stochastic = make_stochastic_context(n_stages);
    let config = SimulationConfig {
        n_scenarios: 4,
        io_channel_capacity: 16,
        basis_activity_window: cobre_sddp::DEFAULT_BASIS_ACTIVITY_WINDOW,
    };
    let horizon = HorizonMode::Finite {
        num_stages: n_stages,
    };
    let initial_state = vec![50.0_f64]; // n_state=1

    let solution = fixed_solution(100.0, 30.0); // objective=100, theta=30
    let solver = MockSolver::always_ok(solution);
    let comm = StubComm { rank: 0, size: 1 };
    let entity_counts = entity_counts_1_hydro();

    let (tx, rx) = mpsc::sync_channel(16);

    let hprod = hydro_productivities_1hydro(n_stages);
    let ec = zero_energy_conversion(1, n_stages);
    let mut workspaces = single_workspace(solver);
    let result = cobre_sddp::simulate(
        &mut workspaces,
        &StageContext {
            templates: &templates,
            base_rows: &base_rows,
            noise_scale: &[],
            n_hydros: 0,
            n_load_buses: 0,
            load_balance_row_starts: &[],
            load_bus_indices: &[],
            block_counts_per_stage: &[],
            ncs_max_gen: &[],
            discount_factors: &[],
            cumulative_discount_factors: &[],
            stage_lag_transitions: &[],
            noise_group_ids: &[],
            downstream_par_order: 0,
        },
        &fcf,
        &TrainingContext {
            horizon: &horizon,
            indexer: &indexer,
            inflow_method: &InflowNonNegativityMethod::None,
            stochastic: &stochastic,
            initial_state: &initial_state,
            inflow_scheme: SamplingScheme::InSample,
            load_scheme: SamplingScheme::InSample,
            ncs_scheme: SamplingScheme::InSample,
            stages: &[],
            historical_library: None,
            external_inflow_library: None,
            external_load_library: None,
            external_ncs_library: None,
            recent_accum_seed: &[],
            recent_weight_seed: 0.0,
        },
        &config,
        SimulationOutputSpec {
            result_tx: &tx,
            zeta_per_stage: &[],
            block_hours_per_stage: &[],
            entity_counts: &entity_counts,
            generic_constraint_row_entries: &[],
            ncs_col_starts: &[],
            n_ncs_per_stage: &[],
            ncs_entity_ids_per_stage: &[],
            diversion_upstream: &HashMap::new(),
            hydro_productivities_per_stage: &hprod,
            energy_conversion: &ec,
            hydro_min_storage_hm3: &[0.0],
            event_sender: None,
        },
        None,
        &[],
        &comm,
    );

    assert!(result.is_ok(), "simulate returned error: {result:?}");
    let run_result = result.unwrap();
    assert_eq!(
        run_result.costs.len(),
        4,
        "cost buffer should have 4 entries"
    );

    // Drain the channel and count results.
    let mut received = 0;
    while rx.try_recv().is_ok() {
        received += 1;
    }
    assert_eq!(received, 4, "channel should have received 4 results");
}

/// Acceptance criterion: solver infeasible at scenario 2, stage 1 (0-based)
/// → `SimulationError::LpInfeasible` with correct `scenario_id` and `stage_id`.
///
/// With 4 scenarios and 2 stages, the solve calls are numbered 0..7 in
/// scenario-outer, stage-inner order:
///   scenario 0: solves 0, 1
///   scenario 1: solves 2, 3
///   scenario 2: solves 4 (stage 0), 5 (stage 1)  ← infeasible at call 5
///   scenario 3: solves 6, 7
///
/// Infeasible at call 5 = `scenario_id=2`, `stage_id=1`.
#[test]
fn simulate_infeasible_returns_lp_infeasible_error() {
    let n_stages = 2;
    let templates: Vec<StageTemplate> = (0..n_stages).map(|_| minimal_template_1_0()).collect();
    let base_rows: Vec<usize> = vec![0; n_stages];

    let indexer = StageIndexer::new(1, 0);
    let fcf = FutureCostFunction::new(n_stages, 1, 1, 10, &vec![0; n_stages]);
    let stochastic = make_stochastic_context(n_stages);
    let config = SimulationConfig {
        n_scenarios: 4,
        io_channel_capacity: 16,
        basis_activity_window: cobre_sddp::DEFAULT_BASIS_ACTIVITY_WINDOW,
    };
    let horizon = HorizonMode::Finite {
        num_stages: n_stages,
    };
    let initial_state = vec![50.0_f64];

    let solution = fixed_solution(100.0, 30.0);
    // Call 5 = scenario_id=2 (0-indexed), stage=1 (0-indexed)
    let solver = MockSolver::infeasible_on(solution, 5);
    let comm = StubComm { rank: 0, size: 1 };
    let entity_counts = entity_counts_1_hydro();

    let (tx, _rx) = mpsc::sync_channel(16);

    let hprod = hydro_productivities_1hydro(n_stages);
    let ec = zero_energy_conversion(1, n_stages);
    let mut workspaces = single_workspace(solver);
    let result = cobre_sddp::simulate(
        &mut workspaces,
        &StageContext {
            templates: &templates,
            base_rows: &base_rows,
            noise_scale: &[],
            n_hydros: 0,
            n_load_buses: 0,
            load_balance_row_starts: &[],
            load_bus_indices: &[],
            block_counts_per_stage: &[],
            ncs_max_gen: &[],
            discount_factors: &[],
            cumulative_discount_factors: &[],
            stage_lag_transitions: &[],
            noise_group_ids: &[],
            downstream_par_order: 0,
        },
        &fcf,
        &TrainingContext {
            horizon: &horizon,
            indexer: &indexer,
            inflow_method: &InflowNonNegativityMethod::None,
            stochastic: &stochastic,
            initial_state: &initial_state,
            inflow_scheme: SamplingScheme::InSample,
            load_scheme: SamplingScheme::InSample,
            ncs_scheme: SamplingScheme::InSample,
            stages: &[],
            historical_library: None,
            external_inflow_library: None,
            external_load_library: None,
            external_ncs_library: None,
            recent_accum_seed: &[],
            recent_weight_seed: 0.0,
        },
        &config,
        SimulationOutputSpec {
            result_tx: &tx,
            zeta_per_stage: &[],
            block_hours_per_stage: &[],
            entity_counts: &entity_counts,
            generic_constraint_row_entries: &[],
            ncs_col_starts: &[],
            n_ncs_per_stage: &[],
            ncs_entity_ids_per_stage: &[],
            diversion_upstream: &HashMap::new(),
            hydro_productivities_per_stage: &hprod,
            energy_conversion: &ec,
            hydro_min_storage_hm3: &[0.0],
            event_sender: None,
        },
        None,
        &[],
        &comm,
    );

    match result {
        Err(SimulationError::LpInfeasible {
            scenario_id,
            stage_id,
            ..
        }) => {
            assert_eq!(scenario_id, 2, "expected scenario_id=2, got {scenario_id}");
            assert_eq!(stage_id, 1, "expected stage_id=1, got {stage_id}");
        }
        other => panic!("expected LpInfeasible, got {other:?}"),
    }
}

/// solver infeasible at scenario 2, stage 3
/// with 4 scenarios and 4 stages → `SimulationError::LpInfeasible { scenario_id: 2, stage_id: 3 }`.
///
/// Solve call index for (scenario=2, stage=3) = 2*4 + 3 = 11 (0-based).
#[test]
fn simulate_infeasible_at_scenario2_stage3() {
    let n_stages = 4;
    let templates: Vec<StageTemplate> = (0..n_stages).map(|_| minimal_template_1_0()).collect();
    let base_rows: Vec<usize> = vec![0; n_stages];

    let indexer = StageIndexer::new(1, 0);
    let fcf = FutureCostFunction::new(n_stages, 1, 1, 10, &vec![0; n_stages]);
    let stochastic = make_stochastic_context(n_stages);
    let config = SimulationConfig {
        n_scenarios: 4,
        io_channel_capacity: 16,
        basis_activity_window: cobre_sddp::DEFAULT_BASIS_ACTIVITY_WINDOW,
    };
    let horizon = HorizonMode::Finite {
        num_stages: n_stages,
    };
    let initial_state = vec![50.0_f64];

    let solution = fixed_solution(100.0, 30.0);
    // Call 11 = scenario 2 (0-based), stage 3 (0-based): 2*4 + 3 = 11.
    let solver = MockSolver::infeasible_on(solution, 11);
    let comm = StubComm { rank: 0, size: 1 };
    let entity_counts = entity_counts_1_hydro();

    let (tx, _rx) = mpsc::sync_channel(16);

    let hprod = hydro_productivities_1hydro(n_stages);
    let ec = zero_energy_conversion(1, n_stages);
    let mut workspaces = single_workspace(solver);
    let result = cobre_sddp::simulate(
        &mut workspaces,
        &StageContext {
            templates: &templates,
            base_rows: &base_rows,
            noise_scale: &[],
            n_hydros: 0,
            n_load_buses: 0,
            load_balance_row_starts: &[],
            load_bus_indices: &[],
            block_counts_per_stage: &[],
            ncs_max_gen: &[],
            discount_factors: &[],
            cumulative_discount_factors: &[],
            stage_lag_transitions: &[],
            noise_group_ids: &[],
            downstream_par_order: 0,
        },
        &fcf,
        &TrainingContext {
            horizon: &horizon,
            indexer: &indexer,
            inflow_method: &InflowNonNegativityMethod::None,
            stochastic: &stochastic,
            initial_state: &initial_state,
            inflow_scheme: SamplingScheme::InSample,
            load_scheme: SamplingScheme::InSample,
            ncs_scheme: SamplingScheme::InSample,
            stages: &[],
            historical_library: None,
            external_inflow_library: None,
            external_load_library: None,
            external_ncs_library: None,
            recent_accum_seed: &[],
            recent_weight_seed: 0.0,
        },
        &config,
        SimulationOutputSpec {
            result_tx: &tx,
            zeta_per_stage: &[],
            block_hours_per_stage: &[],
            entity_counts: &entity_counts,
            generic_constraint_row_entries: &[],
            ncs_col_starts: &[],
            n_ncs_per_stage: &[],
            ncs_entity_ids_per_stage: &[],
            diversion_upstream: &HashMap::new(),
            hydro_productivities_per_stage: &hprod,
            energy_conversion: &ec,
            hydro_min_storage_hm3: &[0.0],
            event_sender: None,
        },
        None,
        &[],
        &comm,
    );

    match result {
        Err(SimulationError::LpInfeasible {
            scenario_id,
            stage_id,
            ..
        }) => {
            assert_eq!(scenario_id, 2, "expected scenario_id=2, got {scenario_id}");
            assert_eq!(stage_id, 3, "expected stage_id=3, got {stage_id}");
        }
        other => panic!("expected LpInfeasible, got {other:?}"),
    }
}

/// Acceptance criterion: drop receiver before calling simulate → `ChannelClosed`.
#[test]
fn simulate_channel_closed_returns_error() {
    let n_stages = 2;
    let templates: Vec<StageTemplate> = (0..n_stages).map(|_| minimal_template_1_0()).collect();
    let base_rows: Vec<usize> = vec![0; n_stages];

    let indexer = StageIndexer::new(1, 0);
    let fcf = FutureCostFunction::new(n_stages, 1, 1, 10, &vec![0; n_stages]);
    let stochastic = make_stochastic_context(n_stages);
    let config = SimulationConfig {
        n_scenarios: 2,
        io_channel_capacity: 1,
        basis_activity_window: cobre_sddp::DEFAULT_BASIS_ACTIVITY_WINDOW,
    };
    let horizon = HorizonMode::Finite {
        num_stages: n_stages,
    };
    let initial_state = vec![50.0_f64];

    let solution = fixed_solution(100.0, 30.0);
    let solver = MockSolver::always_ok(solution);
    let comm = StubComm { rank: 0, size: 1 };
    let entity_counts = entity_counts_1_hydro();

    let (tx, rx) = mpsc::sync_channel(1);
    // Drop the receiver immediately so send() will fail.
    drop(rx);

    let hprod = hydro_productivities_1hydro(n_stages);
    let ec = zero_energy_conversion(1, n_stages);
    let mut workspaces = single_workspace(solver);
    let result = cobre_sddp::simulate(
        &mut workspaces,
        &StageContext {
            templates: &templates,
            base_rows: &base_rows,
            noise_scale: &[],
            n_hydros: 0,
            n_load_buses: 0,
            load_balance_row_starts: &[],
            load_bus_indices: &[],
            block_counts_per_stage: &[],
            ncs_max_gen: &[],
            discount_factors: &[],
            cumulative_discount_factors: &[],
            stage_lag_transitions: &[],
            noise_group_ids: &[],
            downstream_par_order: 0,
        },
        &fcf,
        &TrainingContext {
            horizon: &horizon,
            indexer: &indexer,
            inflow_method: &InflowNonNegativityMethod::None,
            stochastic: &stochastic,
            initial_state: &initial_state,
            inflow_scheme: SamplingScheme::InSample,
            load_scheme: SamplingScheme::InSample,
            ncs_scheme: SamplingScheme::InSample,
            stages: &[],
            historical_library: None,
            external_inflow_library: None,
            external_load_library: None,
            external_ncs_library: None,
            recent_accum_seed: &[],
            recent_weight_seed: 0.0,
        },
        &config,
        SimulationOutputSpec {
            result_tx: &tx,
            zeta_per_stage: &[],
            block_hours_per_stage: &[],
            entity_counts: &entity_counts,
            generic_constraint_row_entries: &[],
            ncs_col_starts: &[],
            n_ncs_per_stage: &[],
            ncs_entity_ids_per_stage: &[],
            diversion_upstream: &HashMap::new(),
            hydro_productivities_per_stage: &hprod,
            energy_conversion: &ec,
            hydro_min_storage_hm3: &[0.0],
            event_sender: None,
        },
        None,
        &[],
        &comm,
    );

    assert!(
        matches!(result, Err(SimulationError::ChannelClosed)),
        "expected ChannelClosed, got {result:?}"
    );
}

/// Acceptance criterion: `total_cost` in cost buffer equals sum of
/// `(objective - primal[theta])` across all stages for each scenario.
///
/// With objective=100.0 and theta=30.0: `stage_cost` = (100 - 30) * `COST_SCALE_FACTOR` = `70_000` per stage.
/// For 3 stages: `total_cost` = 3 \* `70_000` = `210_000`.
#[test]
fn simulate_total_cost_equals_sum_of_stage_costs() {
    let n_stages = 3;
    let templates: Vec<StageTemplate> = (0..n_stages).map(|_| minimal_template_1_0()).collect();
    let base_rows: Vec<usize> = vec![0; n_stages];

    let indexer = StageIndexer::new(1, 0); // theta=2
    let fcf = FutureCostFunction::new(n_stages, 1, 1, 10, &vec![0; n_stages]);
    let stochastic = make_stochastic_context(n_stages);
    let config = SimulationConfig {
        n_scenarios: 2,
        io_channel_capacity: 16,
        basis_activity_window: cobre_sddp::DEFAULT_BASIS_ACTIVITY_WINDOW,
    };
    let horizon = HorizonMode::Finite {
        num_stages: n_stages,
    };
    let initial_state = vec![50.0_f64];

    let objective = 100.0_f64;
    let theta_val = 30.0_f64;
    // stage_cost = (objective - theta) * COST_SCALE_FACTOR = 70 * 1000 = 70_000
    let expected_stage_cost = (objective - theta_val) * 1000.0; // 70_000.0
    #[allow(clippy::cast_precision_loss)]
    let expected_total_cost = expected_stage_cost * n_stages as f64; // 210_000.0

    let solution = fixed_solution(objective, theta_val);
    let solver = MockSolver::always_ok(solution);
    let comm = StubComm { rank: 0, size: 1 };
    let entity_counts = entity_counts_1_hydro();

    let (tx, _rx) = mpsc::sync_channel(16);

    let hprod = hydro_productivities_1hydro(n_stages);
    let ec = zero_energy_conversion(1, n_stages);
    let mut workspaces = single_workspace(solver);
    let run_result = cobre_sddp::simulate(
        &mut workspaces,
        &StageContext {
            templates: &templates,
            base_rows: &base_rows,
            noise_scale: &[],
            n_hydros: 0,
            n_load_buses: 0,
            load_balance_row_starts: &[],
            load_bus_indices: &[],
            block_counts_per_stage: &[],
            ncs_max_gen: &[],
            discount_factors: &[],
            cumulative_discount_factors: &[],
            stage_lag_transitions: &[],
            noise_group_ids: &[],
            downstream_par_order: 0,
        },
        &fcf,
        &TrainingContext {
            horizon: &horizon,
            indexer: &indexer,
            inflow_method: &InflowNonNegativityMethod::None,
            stochastic: &stochastic,
            initial_state: &initial_state,
            inflow_scheme: SamplingScheme::InSample,
            load_scheme: SamplingScheme::InSample,
            ncs_scheme: SamplingScheme::InSample,
            stages: &[],
            historical_library: None,
            external_inflow_library: None,
            external_load_library: None,
            external_ncs_library: None,
            recent_accum_seed: &[],
            recent_weight_seed: 0.0,
        },
        &config,
        SimulationOutputSpec {
            result_tx: &tx,
            zeta_per_stage: &[],
            block_hours_per_stage: &[],
            entity_counts: &entity_counts,
            generic_constraint_row_entries: &[],
            ncs_col_starts: &[],
            n_ncs_per_stage: &[],
            ncs_entity_ids_per_stage: &[],
            diversion_upstream: &HashMap::new(),
            hydro_productivities_per_stage: &hprod,
            energy_conversion: &ec,
            hydro_min_storage_hm3: &[0.0],
            event_sender: None,
        },
        None,
        &[],
        &comm,
    )
    .unwrap();

    assert_eq!(run_result.costs.len(), 2);
    for (scenario_id, total_cost, _) in &run_result.costs {
        assert!(
            (total_cost - expected_total_cost).abs() < 1e-9,
            "scenario {scenario_id}: expected total_cost={expected_total_cost}, got {total_cost}"
        );
    }
}

/// Verify that the `scenario_ids` in the cost buffer match the assigned range.
///
/// With `n_scenarios=6`, `world_size=2`, rank=0: `assign_scenarios(6, 0, 2) = 0..3`.
/// The cost buffer must contain `scenario_ids` 0, 1, 2 in that order.
#[test]
fn simulate_cost_buffer_scenario_ids_match_assigned_range() {
    let n_stages = 1;
    let templates: Vec<StageTemplate> = (0..n_stages).map(|_| minimal_template_1_0()).collect();
    let base_rows: Vec<usize> = vec![0; n_stages];

    let indexer = StageIndexer::new(1, 0);
    let fcf = FutureCostFunction::new(n_stages, 1, 1, 10, &vec![0; n_stages]);
    let stochastic = make_stochastic_context(n_stages);
    let config = SimulationConfig {
        n_scenarios: 6,
        io_channel_capacity: 16,
        basis_activity_window: cobre_sddp::DEFAULT_BASIS_ACTIVITY_WINDOW,
    };
    let horizon = HorizonMode::Finite {
        num_stages: n_stages,
    };
    let initial_state = vec![50.0_f64];

    let solution = fixed_solution(50.0, 10.0);
    let solver = MockSolver::always_ok(solution);
    // rank=0 of 2: assign_scenarios(6, 0, 2) = 0..3
    let comm = StubComm { rank: 0, size: 2 };
    let entity_counts = entity_counts_1_hydro();

    let (tx, _rx) = mpsc::sync_channel(16);

    let hprod = hydro_productivities_1hydro(n_stages);
    let ec = zero_energy_conversion(1, n_stages);
    let mut workspaces = single_workspace(solver);
    let run_result = cobre_sddp::simulate(
        &mut workspaces,
        &StageContext {
            templates: &templates,
            base_rows: &base_rows,
            noise_scale: &[],
            n_hydros: 0,
            n_load_buses: 0,
            load_balance_row_starts: &[],
            load_bus_indices: &[],
            block_counts_per_stage: &[],
            ncs_max_gen: &[],
            discount_factors: &[],
            cumulative_discount_factors: &[],
            stage_lag_transitions: &[],
            noise_group_ids: &[],
            downstream_par_order: 0,
        },
        &fcf,
        &TrainingContext {
            horizon: &horizon,
            indexer: &indexer,
            inflow_method: &InflowNonNegativityMethod::None,
            stochastic: &stochastic,
            initial_state: &initial_state,
            inflow_scheme: SamplingScheme::InSample,
            load_scheme: SamplingScheme::InSample,
            ncs_scheme: SamplingScheme::InSample,
            stages: &[],
            historical_library: None,
            external_inflow_library: None,
            external_load_library: None,
            external_ncs_library: None,
            recent_accum_seed: &[],
            recent_weight_seed: 0.0,
        },
        &config,
        SimulationOutputSpec {
            result_tx: &tx,
            zeta_per_stage: &[],
            block_hours_per_stage: &[],
            entity_counts: &entity_counts,
            generic_constraint_row_entries: &[],
            ncs_col_starts: &[],
            n_ncs_per_stage: &[],
            ncs_entity_ids_per_stage: &[],
            diversion_upstream: &HashMap::new(),
            hydro_productivities_per_stage: &hprod,
            energy_conversion: &ec,
            hydro_min_storage_hm3: &[0.0],
            event_sender: None,
        },
        None,
        &[],
        &comm,
    )
    .unwrap();

    assert_eq!(
        run_result.costs.len(),
        3,
        "rank 0 should process 3 scenarios"
    );
    let ids: Vec<u32> = run_result.costs.iter().map(|(id, _, _)| *id).collect();
    assert_eq!(
        ids,
        vec![0, 1, 2],
        "scenario IDs must match assigned range 0..3"
    );
}

/// Verify channel receives results in scenario order for single rank.
#[test]
fn simulate_channel_receives_results_in_scenario_order() {
    let n_stages = 1;
    let templates: Vec<StageTemplate> = (0..n_stages).map(|_| minimal_template_1_0()).collect();
    let base_rows: Vec<usize> = vec![0; n_stages];

    let indexer = StageIndexer::new(1, 0);
    let fcf = FutureCostFunction::new(n_stages, 1, 1, 10, &vec![0; n_stages]);
    let stochastic = make_stochastic_context(n_stages);
    let config = SimulationConfig {
        n_scenarios: 3,
        io_channel_capacity: 16,
        basis_activity_window: cobre_sddp::DEFAULT_BASIS_ACTIVITY_WINDOW,
    };
    let horizon = HorizonMode::Finite {
        num_stages: n_stages,
    };
    let initial_state = vec![50.0_f64];

    let solution = fixed_solution(100.0, 20.0);
    let solver = MockSolver::always_ok(solution);
    let comm = StubComm { rank: 0, size: 1 };
    let entity_counts = entity_counts_1_hydro();

    let (tx, rx) = mpsc::sync_channel(16);

    let hprod = hydro_productivities_1hydro(n_stages);
    let ec = zero_energy_conversion(1, n_stages);
    let mut workspaces = single_workspace(solver);
    cobre_sddp::simulate(
        &mut workspaces,
        &StageContext {
            templates: &templates,
            base_rows: &base_rows,
            noise_scale: &[],
            n_hydros: 0,
            n_load_buses: 0,
            load_balance_row_starts: &[],
            load_bus_indices: &[],
            block_counts_per_stage: &[],
            ncs_max_gen: &[],
            discount_factors: &[],
            cumulative_discount_factors: &[],
            stage_lag_transitions: &[],
            noise_group_ids: &[],
            downstream_par_order: 0,
        },
        &fcf,
        &TrainingContext {
            horizon: &horizon,
            indexer: &indexer,
            inflow_method: &InflowNonNegativityMethod::None,
            stochastic: &stochastic,
            initial_state: &initial_state,
            inflow_scheme: SamplingScheme::InSample,
            load_scheme: SamplingScheme::InSample,
            ncs_scheme: SamplingScheme::InSample,
            stages: &[],
            historical_library: None,
            external_inflow_library: None,
            external_load_library: None,
            external_ncs_library: None,
            recent_accum_seed: &[],
            recent_weight_seed: 0.0,
        },
        &config,
        SimulationOutputSpec {
            result_tx: &tx,
            zeta_per_stage: &[],
            block_hours_per_stage: &[],
            entity_counts: &entity_counts,
            generic_constraint_row_entries: &[],
            ncs_col_starts: &[],
            n_ncs_per_stage: &[],
            ncs_entity_ids_per_stage: &[],
            diversion_upstream: &HashMap::new(),
            hydro_productivities_per_stage: &hprod,
            energy_conversion: &ec,
            hydro_min_storage_hm3: &[0.0],
            event_sender: None,
        },
        None,
        &[],
        &comm,
    )
    .unwrap();

    let received: Vec<u32> = (0..3).map(|_| rx.recv().unwrap().scenario_id).collect();
    assert_eq!(received, vec![0, 1, 2]);
}

/// New acceptance criterion: cost buffer from 1 workspace equals cost buffer from 4 workspaces.
///
/// Both runs must produce identical `(scenario_id, total_cost, category_costs)` tuples for all
/// 20 scenarios. The cost buffer must be sorted by `scenario_id` in both cases.
#[test]
fn test_simulation_parallel_cost_determinism() {
    let n_stages = 2;
    let n_scenarios = 20u32;
    let templates: Vec<StageTemplate> = (0..n_stages).map(|_| minimal_template_1_0()).collect();
    let base_rows: Vec<usize> = vec![0; n_stages];

    let indexer = StageIndexer::new(1, 0);
    let fcf = FutureCostFunction::new(n_stages, 1, 1, 10, &vec![0; n_stages]);
    let stochastic = make_stochastic_context(n_stages);
    let config = SimulationConfig {
        n_scenarios,
        io_channel_capacity: 64,
        basis_activity_window: cobre_sddp::DEFAULT_BASIS_ACTIVITY_WINDOW,
    };
    let horizon = HorizonMode::Finite {
        num_stages: n_stages,
    };
    let initial_state = vec![50.0_f64];

    let objective = 100.0_f64;
    let theta_val = 30.0_f64;
    let solution = fixed_solution(objective, theta_val);
    let comm = StubComm { rank: 0, size: 1 };
    let entity_counts = entity_counts_1_hydro();

    let hprod = hydro_productivities_1hydro(n_stages);
    let ec = zero_energy_conversion(1, n_stages);

    // Run with 1 workspace.
    let (tx1, _rx1) = mpsc::sync_channel(64);
    let mut workspaces_1 = single_workspace(MockSolver::always_ok(solution.clone()));
    let result_1 = cobre_sddp::simulate(
        &mut workspaces_1,
        &StageContext {
            templates: &templates,
            base_rows: &base_rows,
            noise_scale: &[],
            n_hydros: 0,
            n_load_buses: 0,
            load_balance_row_starts: &[],
            load_bus_indices: &[],
            block_counts_per_stage: &[],
            ncs_max_gen: &[],
            discount_factors: &[],
            cumulative_discount_factors: &[],
            stage_lag_transitions: &[],
            noise_group_ids: &[],
            downstream_par_order: 0,
        },
        &fcf,
        &TrainingContext {
            horizon: &horizon,
            indexer: &indexer,
            inflow_method: &InflowNonNegativityMethod::None,
            stochastic: &stochastic,
            initial_state: &initial_state,
            inflow_scheme: SamplingScheme::InSample,
            load_scheme: SamplingScheme::InSample,
            ncs_scheme: SamplingScheme::InSample,
            stages: &[],
            historical_library: None,
            external_inflow_library: None,
            external_load_library: None,
            external_ncs_library: None,
            recent_accum_seed: &[],
            recent_weight_seed: 0.0,
        },
        &config,
        SimulationOutputSpec {
            result_tx: &tx1,
            zeta_per_stage: &[],
            block_hours_per_stage: &[],
            entity_counts: &entity_counts,
            generic_constraint_row_entries: &[],
            ncs_col_starts: &[],
            n_ncs_per_stage: &[],
            ncs_entity_ids_per_stage: &[],
            diversion_upstream: &HashMap::new(),
            hydro_productivities_per_stage: &hprod,
            energy_conversion: &ec,
            hydro_min_storage_hm3: &[0.0],
            event_sender: None,
        },
        None,
        &[],
        &comm,
    )
    .unwrap();

    // Run with 4 workspaces.
    let (tx4, _rx4) = mpsc::sync_channel(64);
    let mut workspaces_4: Vec<SolverWorkspace<MockSolver>> = (0..4_i32)
        .map(|idx| {
            SolverWorkspace::new(
                0,
                idx,
                MockSolver::always_ok(solution.clone()),
                PatchBuffer::new(1, 0, 0, 0),
                1,
                WorkspaceSizing {
                    hydro_count: 1,
                    ..WorkspaceSizing::default()
                },
            )
        })
        .collect();
    let result_4 = cobre_sddp::simulate(
        &mut workspaces_4,
        &StageContext {
            templates: &templates,
            base_rows: &base_rows,
            noise_scale: &[],
            n_hydros: 0,
            n_load_buses: 0,
            load_balance_row_starts: &[],
            load_bus_indices: &[],
            block_counts_per_stage: &[],
            ncs_max_gen: &[],
            discount_factors: &[],
            cumulative_discount_factors: &[],
            stage_lag_transitions: &[],
            noise_group_ids: &[],
            downstream_par_order: 0,
        },
        &fcf,
        &TrainingContext {
            horizon: &horizon,
            indexer: &indexer,
            inflow_method: &InflowNonNegativityMethod::None,
            stochastic: &stochastic,
            initial_state: &initial_state,
            inflow_scheme: SamplingScheme::InSample,
            load_scheme: SamplingScheme::InSample,
            ncs_scheme: SamplingScheme::InSample,
            stages: &[],
            historical_library: None,
            external_inflow_library: None,
            external_load_library: None,
            external_ncs_library: None,
            recent_accum_seed: &[],
            recent_weight_seed: 0.0,
        },
        &config,
        SimulationOutputSpec {
            result_tx: &tx4,
            zeta_per_stage: &[],
            block_hours_per_stage: &[],
            entity_counts: &entity_counts,
            generic_constraint_row_entries: &[],
            ncs_col_starts: &[],
            n_ncs_per_stage: &[],
            ncs_entity_ids_per_stage: &[],
            diversion_upstream: &HashMap::new(),
            hydro_productivities_per_stage: &hprod,
            energy_conversion: &ec,
            hydro_min_storage_hm3: &[0.0],
            event_sender: None,
        },
        None,
        &[],
        &comm,
    )
    .unwrap();

    let costs_1 = &result_1.costs;
    let costs_4 = &result_4.costs;

    // Both cost buffers must have exactly 20 entries sorted by scenario_id.
    assert_eq!(
        costs_1.len(),
        n_scenarios as usize,
        "1-workspace: 20 entries"
    );
    assert_eq!(
        costs_4.len(),
        n_scenarios as usize,
        "4-workspace: 20 entries"
    );

    let ids_1: Vec<u32> = costs_1.iter().map(|(id, _, _)| *id).collect();
    let ids_4: Vec<u32> = costs_4.iter().map(|(id, _, _)| *id).collect();
    let expected_ids: Vec<u32> = (0..n_scenarios).collect();
    assert_eq!(ids_1, expected_ids, "1-workspace: sorted scenario IDs");
    assert_eq!(ids_4, expected_ids, "4-workspace: sorted scenario IDs");

    // Cost values must be identical.
    for i in 0..n_scenarios as usize {
        let (id1, cost1, _) = &costs_1[i];
        let (id4, cost4, _) = &costs_4[i];
        assert_eq!(id1, id4, "scenario_id mismatch at index {i}");
        assert!(
            (cost1 - cost4).abs() < 1e-9,
            "cost mismatch for scenario {id1}: 1-ws={cost1}, 4-ws={cost4}"
        );
    }
}

// ── Integration tests for event emission ─────────────────────────────────

/// Acceptance criterion: with `event_sender: Some(&tx)` and 10 scenarios,
/// at least 1 `SimulationProgress` event is received with `scenarios_complete > 0`
/// and a finite non-NaN `scenario_cost`.
#[test]
fn simulate_emits_progress_events() {
    use cobre_core::TrainingEvent;

    let n_stages = 2;
    let templates: Vec<StageTemplate> = (0..n_stages).map(|_| minimal_template_1_0()).collect();
    let base_rows: Vec<usize> = vec![0; n_stages];

    let indexer = StageIndexer::new(1, 0);
    let fcf = FutureCostFunction::new(n_stages, 1, 1, 10, &vec![0; n_stages]);
    let stochastic = make_stochastic_context(n_stages);
    let config = SimulationConfig {
        n_scenarios: 10,
        io_channel_capacity: 32,
        basis_activity_window: cobre_sddp::DEFAULT_BASIS_ACTIVITY_WINDOW,
    };
    let horizon = HorizonMode::Finite {
        num_stages: n_stages,
    };
    let initial_state = vec![50.0_f64];

    let solution = fixed_solution(100.0, 30.0);
    let solver = MockSolver::always_ok(solution);
    let comm = StubComm { rank: 0, size: 1 };
    let entity_counts = entity_counts_1_hydro();

    let (result_tx, _result_rx) = mpsc::sync_channel(32);
    let (event_tx, event_rx) = mpsc::channel::<TrainingEvent>();

    let hprod = hydro_productivities_1hydro(n_stages);
    let ec = zero_energy_conversion(1, n_stages);
    let mut workspaces = single_workspace(solver);
    let result = cobre_sddp::simulate(
        &mut workspaces,
        &StageContext {
            templates: &templates,
            base_rows: &base_rows,
            noise_scale: &[],
            n_hydros: 0,
            n_load_buses: 0,
            load_balance_row_starts: &[],
            load_bus_indices: &[],
            block_counts_per_stage: &[],
            ncs_max_gen: &[],
            discount_factors: &[],
            cumulative_discount_factors: &[],
            stage_lag_transitions: &[],
            noise_group_ids: &[],
            downstream_par_order: 0,
        },
        &fcf,
        &TrainingContext {
            horizon: &horizon,
            indexer: &indexer,
            inflow_method: &InflowNonNegativityMethod::None,
            stochastic: &stochastic,
            initial_state: &initial_state,
            inflow_scheme: SamplingScheme::InSample,
            load_scheme: SamplingScheme::InSample,
            ncs_scheme: SamplingScheme::InSample,
            stages: &[],
            historical_library: None,
            external_inflow_library: None,
            external_load_library: None,
            external_ncs_library: None,
            recent_accum_seed: &[],
            recent_weight_seed: 0.0,
        },
        &config,
        SimulationOutputSpec {
            result_tx: &result_tx,
            zeta_per_stage: &[],
            block_hours_per_stage: &[],
            entity_counts: &entity_counts,
            generic_constraint_row_entries: &[],
            ncs_col_starts: &[],
            n_ncs_per_stage: &[],
            ncs_entity_ids_per_stage: &[],
            diversion_upstream: &HashMap::new(),
            hydro_productivities_per_stage: &hprod,
            energy_conversion: &ec,
            hydro_min_storage_hm3: &[0.0],
            event_sender: Some(event_tx),
        },
        None,
        &[],
        &comm,
    );
    assert!(result.is_ok(), "simulate returned error: {result:?}");

    let events: Vec<TrainingEvent> = event_rx.iter().collect();

    let progress_events: Vec<_> = events
        .iter()
        .filter(|e| matches!(e, TrainingEvent::SimulationProgress { .. }))
        .collect();

    assert!(
        !progress_events.is_empty(),
        "at least 1 SimulationProgress event expected"
    );

    for event in &progress_events {
        let TrainingEvent::SimulationProgress {
            scenarios_complete,
            scenario_cost,
            ..
        } = event
        else {
            continue;
        };

        assert!(
            *scenarios_complete > 0,
            "scenarios_complete must be > 0, got {scenarios_complete}"
        );
        assert!(
            scenario_cost.is_finite() && !scenario_cost.is_nan(),
            "scenario_cost must be finite and non-NaN, got {scenario_cost}"
        );
    }
}

/// Acceptance criterion: with `event_sender: None`, no events are sent and
/// the function returns the same cost buffer as before.
#[test]
fn simulate_no_events_when_sender_is_none() {
    let n_stages = 2;
    let templates: Vec<StageTemplate> = (0..n_stages).map(|_| minimal_template_1_0()).collect();
    let base_rows: Vec<usize> = vec![0; n_stages];

    let indexer = StageIndexer::new(1, 0);
    let fcf = FutureCostFunction::new(n_stages, 1, 1, 10, &vec![0; n_stages]);
    let stochastic = make_stochastic_context(n_stages);
    let config = SimulationConfig {
        n_scenarios: 4,
        io_channel_capacity: 16,
        basis_activity_window: cobre_sddp::DEFAULT_BASIS_ACTIVITY_WINDOW,
    };
    let horizon = HorizonMode::Finite {
        num_stages: n_stages,
    };
    let initial_state = vec![50.0_f64];

    let solution = fixed_solution(100.0, 30.0);
    let solver = MockSolver::always_ok(solution);
    let comm = StubComm { rank: 0, size: 1 };
    let entity_counts = entity_counts_1_hydro();

    let (result_tx, _result_rx) = mpsc::sync_channel(16);

    let hprod = hydro_productivities_1hydro(n_stages);
    let ec = zero_energy_conversion(1, n_stages);
    let mut workspaces = single_workspace(solver);
    let result = cobre_sddp::simulate(
        &mut workspaces,
        &StageContext {
            templates: &templates,
            base_rows: &base_rows,
            noise_scale: &[],
            n_hydros: 0,
            n_load_buses: 0,
            load_balance_row_starts: &[],
            load_bus_indices: &[],
            block_counts_per_stage: &[],
            ncs_max_gen: &[],
            discount_factors: &[],
            cumulative_discount_factors: &[],
            stage_lag_transitions: &[],
            noise_group_ids: &[],
            downstream_par_order: 0,
        },
        &fcf,
        &TrainingContext {
            horizon: &horizon,
            indexer: &indexer,
            inflow_method: &InflowNonNegativityMethod::None,
            stochastic: &stochastic,
            initial_state: &initial_state,
            inflow_scheme: SamplingScheme::InSample,
            load_scheme: SamplingScheme::InSample,
            ncs_scheme: SamplingScheme::InSample,
            stages: &[],
            historical_library: None,
            external_inflow_library: None,
            external_load_library: None,
            external_ncs_library: None,
            recent_accum_seed: &[],
            recent_weight_seed: 0.0,
        },
        &config,
        SimulationOutputSpec {
            result_tx: &result_tx,
            zeta_per_stage: &[],
            block_hours_per_stage: &[],
            entity_counts: &entity_counts,
            generic_constraint_row_entries: &[],
            ncs_col_starts: &[],
            n_ncs_per_stage: &[],
            ncs_entity_ids_per_stage: &[],
            diversion_upstream: &HashMap::new(),
            hydro_productivities_per_stage: &hprod,
            energy_conversion: &ec,
            hydro_min_storage_hm3: &[0.0],
            event_sender: None,
        },
        None,
        &[],
        &comm,
    );

    assert!(result.is_ok(), "simulate returned error: {result:?}");
    let run_result = result.unwrap();
    assert_eq!(
        run_result.costs.len(),
        4,
        "cost buffer must have 4 entries when event_sender is None"
    );
}

/// `SimulationProgress` events are
/// received in the channel BEFORE `simulate()` returns (during the
/// parallel region).
///
/// With a single workspace (serial rayon execution), the worker emits
/// progress events as each scenario completes. Because events are sent
/// from the closure rather than the post-collect loop, the receiver
/// contains events by the time `simulate()` returns.
#[test]
fn simulate_progress_events_received_before_return() {
    use cobre_core::TrainingEvent;

    let n_stages = 1;
    let n_scenarios = 10;
    let templates: Vec<StageTemplate> = (0..n_stages).map(|_| minimal_template_1_0()).collect();
    let base_rows: Vec<usize> = vec![0; n_stages];

    let indexer = StageIndexer::new(1, 0);
    let fcf = FutureCostFunction::new(n_stages, 1, 1, 10, &vec![0; n_stages]);
    let stochastic = make_stochastic_context(n_stages);
    let config = SimulationConfig {
        n_scenarios,
        io_channel_capacity: 32,
        basis_activity_window: cobre_sddp::DEFAULT_BASIS_ACTIVITY_WINDOW,
    };
    let horizon = HorizonMode::Finite {
        num_stages: n_stages,
    };
    let initial_state = vec![50.0_f64];

    let solution = fixed_solution(100.0, 30.0);
    let solver = MockSolver::always_ok(solution);
    let comm = StubComm { rank: 0, size: 1 };
    let entity_counts = entity_counts_1_hydro();

    let (result_tx, _result_rx) = mpsc::sync_channel(32);
    let (event_tx, event_rx) = mpsc::channel::<TrainingEvent>();

    let hprod = hydro_productivities_1hydro(n_stages);
    let ec = zero_energy_conversion(1, n_stages);
    let mut workspaces = single_workspace(solver);
    cobre_sddp::simulate(
        &mut workspaces,
        &StageContext {
            templates: &templates,
            base_rows: &base_rows,
            noise_scale: &[],
            n_hydros: 0,
            n_load_buses: 0,
            load_balance_row_starts: &[],
            load_bus_indices: &[],
            block_counts_per_stage: &[],
            ncs_max_gen: &[],
            discount_factors: &[],
            cumulative_discount_factors: &[],
            stage_lag_transitions: &[],
            noise_group_ids: &[],
            downstream_par_order: 0,
        },
        &fcf,
        &TrainingContext {
            horizon: &horizon,
            indexer: &indexer,
            inflow_method: &InflowNonNegativityMethod::None,
            stochastic: &stochastic,
            initial_state: &initial_state,
            inflow_scheme: SamplingScheme::InSample,
            load_scheme: SamplingScheme::InSample,
            ncs_scheme: SamplingScheme::InSample,
            stages: &[],
            historical_library: None,
            external_inflow_library: None,
            external_load_library: None,
            external_ncs_library: None,
            recent_accum_seed: &[],
            recent_weight_seed: 0.0,
        },
        &config,
        SimulationOutputSpec {
            result_tx: &result_tx,
            zeta_per_stage: &[],
            block_hours_per_stage: &[],
            entity_counts: &entity_counts,
            generic_constraint_row_entries: &[],
            ncs_col_starts: &[],
            n_ncs_per_stage: &[],
            ncs_entity_ids_per_stage: &[],
            diversion_upstream: &HashMap::new(),
            hydro_productivities_per_stage: &hprod,
            energy_conversion: &ec,
            hydro_min_storage_hm3: &[0.0],
            event_sender: Some(event_tx),
        },
        None,
        &[],
        &comm,
    )
    .unwrap();

    // Because the sender was moved into simulate() and dropped when it
    // returns, the channel is now closed. Collect all events.
    let events: Vec<TrainingEvent> = event_rx.iter().collect();
    let progress_count = events
        .iter()
        .filter(|e| matches!(e, TrainingEvent::SimulationProgress { .. }))
        .count();

    assert!(
        progress_count > 0,
        "expected SimulationProgress events in channel after simulate() returns, got 0"
    );
    assert_eq!(
        progress_count, n_scenarios as usize,
        "expected {n_scenarios} SimulationProgress events (one per scenario), got {progress_count}"
    );
}

/// Acceptance criterion: each `SimulationProgress` event carries the raw
/// `scenario_cost` of the completed scenario.
///
/// With `MockSolver` returning a fixed solution, all scenarios have the same
/// `total_cost`. Validates that every `SimulationProgress.scenario_cost`
/// equals the expected per-scenario cost.
#[test]
fn simulate_progress_scenario_cost_equals_total_cost() {
    use cobre_core::TrainingEvent;

    let n_stages = 1;
    let n_scenarios = 5_u32;
    let templates: Vec<StageTemplate> = (0..n_stages).map(|_| minimal_template_1_0()).collect();
    let base_rows: Vec<usize> = vec![0; n_stages];

    let indexer = StageIndexer::new(1, 0);
    let fcf = FutureCostFunction::new(n_stages, 1, 1, 10, &vec![0; n_stages]);
    let stochastic = make_stochastic_context(n_stages);
    let config = SimulationConfig {
        n_scenarios,
        io_channel_capacity: 32,
        basis_activity_window: cobre_sddp::DEFAULT_BASIS_ACTIVITY_WINDOW,
    };
    let horizon = HorizonMode::Finite {
        num_stages: n_stages,
    };
    let initial_state = vec![50.0_f64];

    // objective=100, theta=30 → stage_cost = (100-30)*COST_SCALE_FACTOR = 70_000.0 every scenario.
    let solution = fixed_solution(100.0, 30.0);
    let expected_stage_cost = 70_000.0_f64;

    let solver = MockSolver::always_ok(solution);
    let comm = StubComm { rank: 0, size: 1 };
    let entity_counts = entity_counts_1_hydro();

    let (result_tx, _result_rx) = mpsc::sync_channel(32);
    let (event_tx, event_rx) = mpsc::channel::<TrainingEvent>();

    let hprod = hydro_productivities_1hydro(n_stages);
    let ec = zero_energy_conversion(1, n_stages);
    let mut workspaces = single_workspace(solver);
    cobre_sddp::simulate(
        &mut workspaces,
        &StageContext {
            templates: &templates,
            base_rows: &base_rows,
            noise_scale: &[],
            n_hydros: 0,
            n_load_buses: 0,
            load_balance_row_starts: &[],
            load_bus_indices: &[],
            block_counts_per_stage: &[],
            ncs_max_gen: &[],
            discount_factors: &[],
            cumulative_discount_factors: &[],
            stage_lag_transitions: &[],
            noise_group_ids: &[],
            downstream_par_order: 0,
        },
        &fcf,
        &TrainingContext {
            horizon: &horizon,
            indexer: &indexer,
            inflow_method: &InflowNonNegativityMethod::None,
            stochastic: &stochastic,
            initial_state: &initial_state,
            inflow_scheme: SamplingScheme::InSample,
            load_scheme: SamplingScheme::InSample,
            ncs_scheme: SamplingScheme::InSample,
            stages: &[],
            historical_library: None,
            external_inflow_library: None,
            external_load_library: None,
            external_ncs_library: None,
            recent_accum_seed: &[],
            recent_weight_seed: 0.0,
        },
        &config,
        SimulationOutputSpec {
            result_tx: &result_tx,
            zeta_per_stage: &[],
            block_hours_per_stage: &[],
            entity_counts: &entity_counts,
            generic_constraint_row_entries: &[],
            ncs_col_starts: &[],
            n_ncs_per_stage: &[],
            ncs_entity_ids_per_stage: &[],
            diversion_upstream: &HashMap::new(),
            hydro_productivities_per_stage: &hprod,
            energy_conversion: &ec,
            hydro_min_storage_hm3: &[0.0],
            event_sender: Some(event_tx),
        },
        None,
        &[],
        &comm,
    )
    .unwrap();

    let events: Vec<TrainingEvent> = event_rx.iter().collect();
    let progress_events: Vec<_> = events
        .iter()
        .filter(|e| matches!(e, TrainingEvent::SimulationProgress { .. }))
        .collect();

    assert_eq!(
        progress_events.len(),
        n_scenarios as usize,
        "expected {n_scenarios} progress events"
    );

    // Every progress event must carry the raw scenario cost.
    for event in &progress_events {
        let TrainingEvent::SimulationProgress { scenario_cost, .. } = event else {
            continue;
        };
        assert!(
            (scenario_cost - expected_stage_cost).abs() < 1e-9,
            "scenario_cost must equal expected cost {expected_stage_cost}, got {scenario_cost}"
        );
    }
}

/// Acceptance criterion: `SimulationFinished` event is the last event
/// emitted after all `SimulationProgress` events.
#[test]
fn simulate_emits_simulation_finished_as_last_event() {
    use cobre_core::TrainingEvent;

    let n_stages = 1;
    let n_scenarios = 6_u32;
    let templates: Vec<StageTemplate> = (0..n_stages).map(|_| minimal_template_1_0()).collect();
    let base_rows: Vec<usize> = vec![0; n_stages];

    let indexer = StageIndexer::new(1, 0);
    let fcf = FutureCostFunction::new(n_stages, 1, 1, 10, &vec![0; n_stages]);
    let stochastic = make_stochastic_context(n_stages);
    let config = SimulationConfig {
        n_scenarios,
        io_channel_capacity: 32,
        basis_activity_window: cobre_sddp::DEFAULT_BASIS_ACTIVITY_WINDOW,
    };
    let horizon = HorizonMode::Finite {
        num_stages: n_stages,
    };
    let initial_state = vec![50.0_f64];

    let solution = fixed_solution(100.0, 30.0);
    let solver = MockSolver::always_ok(solution);
    let comm = StubComm { rank: 0, size: 1 };
    let entity_counts = entity_counts_1_hydro();

    let (result_tx, _result_rx) = mpsc::sync_channel(32);
    let (event_tx, event_rx) = mpsc::channel::<TrainingEvent>();

    let hprod = hydro_productivities_1hydro(n_stages);
    let ec = zero_energy_conversion(1, n_stages);
    let mut workspaces = single_workspace(solver);
    cobre_sddp::simulate(
        &mut workspaces,
        &StageContext {
            templates: &templates,
            base_rows: &base_rows,
            noise_scale: &[],
            n_hydros: 0,
            n_load_buses: 0,
            load_balance_row_starts: &[],
            load_bus_indices: &[],
            block_counts_per_stage: &[],
            ncs_max_gen: &[],
            discount_factors: &[],
            cumulative_discount_factors: &[],
            stage_lag_transitions: &[],
            noise_group_ids: &[],
            downstream_par_order: 0,
        },
        &fcf,
        &TrainingContext {
            horizon: &horizon,
            indexer: &indexer,
            inflow_method: &InflowNonNegativityMethod::None,
            stochastic: &stochastic,
            initial_state: &initial_state,
            inflow_scheme: SamplingScheme::InSample,
            load_scheme: SamplingScheme::InSample,
            ncs_scheme: SamplingScheme::InSample,
            stages: &[],
            historical_library: None,
            external_inflow_library: None,
            external_load_library: None,
            external_ncs_library: None,
            recent_accum_seed: &[],
            recent_weight_seed: 0.0,
        },
        &config,
        SimulationOutputSpec {
            result_tx: &result_tx,
            zeta_per_stage: &[],
            block_hours_per_stage: &[],
            entity_counts: &entity_counts,
            generic_constraint_row_entries: &[],
            ncs_col_starts: &[],
            n_ncs_per_stage: &[],
            ncs_entity_ids_per_stage: &[],
            diversion_upstream: &HashMap::new(),
            hydro_productivities_per_stage: &hprod,
            energy_conversion: &ec,
            hydro_min_storage_hm3: &[0.0],
            event_sender: Some(event_tx),
        },
        None,
        &[],
        &comm,
    )
    .unwrap();

    let events: Vec<TrainingEvent> = event_rx.iter().collect();

    // Must have at least n_scenarios progress events + 1 finished event.
    assert!(
        events.len() > n_scenarios as usize,
        "expected at least {} events, got {}",
        n_scenarios + 1,
        events.len()
    );

    // The last event must be SimulationFinished.
    let last = events.last().unwrap();
    assert!(
        matches!(last, TrainingEvent::SimulationFinished { .. }),
        "last event must be SimulationFinished, got {last:?}"
    );

    // SimulationFinished must carry the correct scenario count.
    let TrainingEvent::SimulationFinished { scenarios, .. } = last else {
        panic!("last event is not SimulationFinished");
    };
    assert_eq!(
        *scenarios, n_scenarios,
        "SimulationFinished.scenarios must equal n_scenarios={n_scenarios}, got {scenarios}"
    );

    // All events before the last must be SimulationProgress.
    let progress_count = events
        .iter()
        .filter(|e| matches!(e, TrainingEvent::SimulationProgress { .. }))
        .count();
    assert_eq!(
        progress_count, n_scenarios as usize,
        "expected {n_scenarios} SimulationProgress events before SimulationFinished"
    );
}

/// Acceptance criterion: each `SimulationProgress` event carries a finite,
/// non-NaN `scenario_cost`. Statistics accumulation is deferred to the
/// progress thread; this test verifies the per-scenario cost
/// field is always valid.
#[test]
fn simulate_progress_scenario_cost_is_finite() {
    use cobre_core::TrainingEvent;

    let n_stages = 1;
    let templates: Vec<StageTemplate> = (0..n_stages).map(|_| minimal_template_1_0()).collect();
    let base_rows: Vec<usize> = vec![0; n_stages];

    let indexer = StageIndexer::new(1, 0);
    let fcf = FutureCostFunction::new(n_stages, 1, 1, 10, &vec![0; n_stages]);
    let stochastic = make_stochastic_context(n_stages);
    let config = SimulationConfig {
        n_scenarios: 5,
        io_channel_capacity: 16,
        basis_activity_window: cobre_sddp::DEFAULT_BASIS_ACTIVITY_WINDOW,
    };
    let horizon = HorizonMode::Finite {
        num_stages: n_stages,
    };
    let initial_state = vec![50.0_f64];

    // All scenarios have cost = objective - theta = 100 - 30 = 70.
    let solution = fixed_solution(100.0, 30.0);
    let solver = MockSolver::always_ok(solution);
    let comm = StubComm { rank: 0, size: 1 };
    let entity_counts = entity_counts_1_hydro();

    let (result_tx, _result_rx) = mpsc::sync_channel(16);
    let (event_tx, event_rx) = mpsc::channel::<TrainingEvent>();

    let hprod = hydro_productivities_1hydro(n_stages);
    let ec = zero_energy_conversion(1, n_stages);
    let mut workspaces = single_workspace(solver);
    cobre_sddp::simulate(
        &mut workspaces,
        &StageContext {
            templates: &templates,
            base_rows: &base_rows,
            noise_scale: &[],
            n_hydros: 0,
            n_load_buses: 0,
            load_balance_row_starts: &[],
            load_bus_indices: &[],
            block_counts_per_stage: &[],
            ncs_max_gen: &[],
            discount_factors: &[],
            cumulative_discount_factors: &[],
            stage_lag_transitions: &[],
            noise_group_ids: &[],
            downstream_par_order: 0,
        },
        &fcf,
        &TrainingContext {
            horizon: &horizon,
            indexer: &indexer,
            inflow_method: &InflowNonNegativityMethod::None,
            stochastic: &stochastic,
            initial_state: &initial_state,
            inflow_scheme: SamplingScheme::InSample,
            load_scheme: SamplingScheme::InSample,
            ncs_scheme: SamplingScheme::InSample,
            stages: &[],
            historical_library: None,
            external_inflow_library: None,
            external_load_library: None,
            external_ncs_library: None,
            recent_accum_seed: &[],
            recent_weight_seed: 0.0,
        },
        &config,
        SimulationOutputSpec {
            result_tx: &result_tx,
            zeta_per_stage: &[],
            block_hours_per_stage: &[],
            entity_counts: &entity_counts,
            generic_constraint_row_entries: &[],
            ncs_col_starts: &[],
            n_ncs_per_stage: &[],
            ncs_entity_ids_per_stage: &[],
            diversion_upstream: &HashMap::new(),
            hydro_productivities_per_stage: &hprod,
            energy_conversion: &ec,
            hydro_min_storage_hm3: &[0.0],
            event_sender: Some(event_tx),
        },
        None,
        &[],
        &comm,
    )
    .unwrap();

    let events: Vec<TrainingEvent> = event_rx.iter().collect();
    let progress_events: Vec<_> = events
        .iter()
        .filter(|e| matches!(e, TrainingEvent::SimulationProgress { .. }))
        .collect();

    // Every scenario_cost must be finite and non-NaN.
    for event in &progress_events {
        let TrainingEvent::SimulationProgress { scenario_cost, .. } = event else {
            continue;
        };
        assert!(
            scenario_cost.is_finite() && !scenario_cost.is_nan(),
            "scenario_cost must be finite and non-NaN, got {scenario_cost}"
        );
    }
}

// ── baked-template acceptance tests ────────────────────────────

/// When `baked_templates` is `Some`,
/// `add_rows` is never called (zero `add_rows_count`) and `load_model` is
/// called exactly `n_scenarios * n_stages` times.
#[test]
fn simulate_baked_path_issues_zero_add_rows() {
    let n_stages = 2;
    let n_scenarios = 3u32;
    let templates: Vec<StageTemplate> = (0..n_stages).map(|_| minimal_template_1_0()).collect();
    // Baked templates: for MockSolver the content does not matter; use the
    // same minimal template as a stand-in for a baked (pre-merged) template.
    let baked: Vec<StageTemplate> = (0..n_stages).map(|_| minimal_template_1_0()).collect();
    let base_rows: Vec<usize> = vec![0; n_stages];

    let indexer = StageIndexer::new(1, 0);
    let fcf = FutureCostFunction::new(n_stages, 1, 1, 10, &vec![0; n_stages]);
    let stochastic = make_stochastic_context(n_stages);
    let config = SimulationConfig {
        n_scenarios,
        io_channel_capacity: 16,
        basis_activity_window: cobre_sddp::DEFAULT_BASIS_ACTIVITY_WINDOW,
    };
    let horizon = HorizonMode::Finite {
        num_stages: n_stages,
    };
    let initial_state = vec![50.0_f64];

    let solution = fixed_solution(100.0, 30.0);
    let solver = MockSolver::always_ok(solution);
    let comm = StubComm { rank: 0, size: 1 };
    let entity_counts = entity_counts_1_hydro();
    let (tx, _rx) = mpsc::sync_channel(32);
    let hprod = hydro_productivities_1hydro(n_stages);
    let ec = zero_energy_conversion(1, n_stages);

    let mut workspaces = single_workspace(solver);
    let result = cobre_sddp::simulate(
        &mut workspaces,
        &StageContext {
            templates: &templates,
            base_rows: &base_rows,
            noise_scale: &[],
            n_hydros: 0,
            n_load_buses: 0,
            load_balance_row_starts: &[],
            load_bus_indices: &[],
            block_counts_per_stage: &[],
            ncs_max_gen: &[],
            discount_factors: &[],
            cumulative_discount_factors: &[],
            stage_lag_transitions: &[],
            noise_group_ids: &[],
            downstream_par_order: 0,
        },
        &fcf,
        &TrainingContext {
            horizon: &horizon,
            indexer: &indexer,
            inflow_method: &InflowNonNegativityMethod::None,
            stochastic: &stochastic,
            initial_state: &initial_state,
            inflow_scheme: SamplingScheme::InSample,
            load_scheme: SamplingScheme::InSample,
            ncs_scheme: SamplingScheme::InSample,
            stages: &[],
            historical_library: None,
            external_inflow_library: None,
            external_load_library: None,
            external_ncs_library: None,
            recent_accum_seed: &[],
            recent_weight_seed: 0.0,
        },
        &config,
        SimulationOutputSpec {
            result_tx: &tx,
            zeta_per_stage: &[],
            block_hours_per_stage: &[],
            entity_counts: &entity_counts,
            generic_constraint_row_entries: &[],
            ncs_col_starts: &[],
            n_ncs_per_stage: &[],
            ncs_entity_ids_per_stage: &[],
            diversion_upstream: &HashMap::new(),
            hydro_productivities_per_stage: &hprod,
            energy_conversion: &ec,
            hydro_min_storage_hm3: &[0.0],
            event_sender: None,
        },
        Some(baked.as_slice()),
        &[],
        &comm,
    );

    assert!(result.is_ok(), "baked path must succeed: {result:?}");
    let expected_load_count = n_scenarios as usize * n_stages;
    let solver = &workspaces[0].solver;
    assert_eq!(
        solver.add_rows_count, 0,
        "baked path must call add_rows 0 times; got {}",
        solver.add_rows_count
    );
    assert_eq!(
        solver.load_count, expected_load_count,
        "baked path must call load_model {} times; got {}",
        expected_load_count, solver.load_count
    );
}

/// When `baked_templates` is `None`
/// (fallback path), `add_rows` is called only when cuts exist.
///
/// The FCF has 0 active cuts, so `cut_batch.num_rows == 0` for every stage,
/// meaning `add_rows` is gated by `if cut_batch.num_rows > 0`. This test
/// verifies the structural contract: with no cuts `add_rows_count == 0` on
/// the fallback path. When `baked_templates` is provided, `load_model` is
/// called against the baked template; when `None`, it is called against
/// `ctx.templates[t]`. The absence of `add_rows` when `num_rows == 0` holds
/// in both cases.
///
/// To test the fallback counter path, we verify that `add_rows_count`
/// remains 0 on a zero-cut fallback run and that `load_count` equals
/// `n_scenarios * n_stages`, confirming the fallback path calls `load_model`
/// the same number of times as the baked path.
#[test]
fn simulate_fallback_path_issues_expected_add_rows() {
    let n_stages = 2;
    let n_scenarios = 3u32;
    let templates: Vec<StageTemplate> = (0..n_stages).map(|_| minimal_template_1_0()).collect();
    let base_rows: Vec<usize> = vec![0; n_stages];

    let indexer = StageIndexer::new(1, 0);
    // FCF with 0 cuts — cut_batch.num_rows will be 0 for every stage.
    let fcf = FutureCostFunction::new(n_stages, 1, 1, 10, &vec![0; n_stages]);
    let stochastic = make_stochastic_context(n_stages);
    let config = SimulationConfig {
        n_scenarios,
        io_channel_capacity: 16,
        basis_activity_window: cobre_sddp::DEFAULT_BASIS_ACTIVITY_WINDOW,
    };
    let horizon = HorizonMode::Finite {
        num_stages: n_stages,
    };
    let initial_state = vec![50.0_f64];

    let solution = fixed_solution(100.0, 30.0);
    let solver = MockSolver::always_ok(solution);
    let comm = StubComm { rank: 0, size: 1 };
    let entity_counts = entity_counts_1_hydro();
    let (tx, _rx) = mpsc::sync_channel(32);
    let hprod = hydro_productivities_1hydro(n_stages);
    let ec = zero_energy_conversion(1, n_stages);

    let mut workspaces = single_workspace(solver);
    let result = cobre_sddp::simulate(
        &mut workspaces,
        &StageContext {
            templates: &templates,
            base_rows: &base_rows,
            noise_scale: &[],
            n_hydros: 0,
            n_load_buses: 0,
            load_balance_row_starts: &[],
            load_bus_indices: &[],
            block_counts_per_stage: &[],
            ncs_max_gen: &[],
            discount_factors: &[],
            cumulative_discount_factors: &[],
            stage_lag_transitions: &[],
            noise_group_ids: &[],
            downstream_par_order: 0,
        },
        &fcf,
        &TrainingContext {
            horizon: &horizon,
            indexer: &indexer,
            inflow_method: &InflowNonNegativityMethod::None,
            stochastic: &stochastic,
            initial_state: &initial_state,
            inflow_scheme: SamplingScheme::InSample,
            load_scheme: SamplingScheme::InSample,
            ncs_scheme: SamplingScheme::InSample,
            stages: &[],
            historical_library: None,
            external_inflow_library: None,
            external_load_library: None,
            external_ncs_library: None,
            recent_accum_seed: &[],
            recent_weight_seed: 0.0,
        },
        &config,
        SimulationOutputSpec {
            result_tx: &tx,
            zeta_per_stage: &[],
            block_hours_per_stage: &[],
            entity_counts: &entity_counts,
            generic_constraint_row_entries: &[],
            ncs_col_starts: &[],
            n_ncs_per_stage: &[],
            ncs_entity_ids_per_stage: &[],
            diversion_upstream: &HashMap::new(),
            hydro_productivities_per_stage: &hprod,
            energy_conversion: &ec,
            hydro_min_storage_hm3: &[0.0],
            event_sender: None,
        },
        // fallback path
        None,
        &[],
        &comm,
    );

    assert!(result.is_ok(), "fallback path must succeed: {result:?}");
    let expected_load_count = n_scenarios as usize * n_stages;
    let solver = &workspaces[0].solver;
    // With zero cuts, `cut_batch.num_rows == 0` so the guard
    // `if cut_batch.num_rows > 0` prevents any `add_rows` call.
    assert_eq!(
        solver.add_rows_count, 0,
        "fallback path with zero cuts must call add_rows 0 times; got {}",
        solver.add_rows_count
    );
    assert_eq!(
        solver.load_count, expected_load_count,
        "fallback path must call load_model {} times; got {}",
        expected_load_count, solver.load_count
    );
}

/// When `baked_templates` is `Some`
/// but the slice length differs from `num_stages`, `simulate` returns
/// `SimulationError::InvalidConfiguration` whose message contains both lengths.
#[test]
fn simulate_baked_length_mismatch_returns_error() {
    let n_stages = 3;
    let templates: Vec<StageTemplate> = (0..n_stages).map(|_| minimal_template_1_0()).collect();
    let base_rows: Vec<usize> = vec![0; n_stages];

    let indexer = StageIndexer::new(1, 0);
    let fcf = FutureCostFunction::new(n_stages, 1, 1, 10, &vec![0; n_stages]);
    let stochastic = make_stochastic_context(n_stages);
    let config = SimulationConfig {
        n_scenarios: 2,
        io_channel_capacity: 8,
        basis_activity_window: cobre_sddp::DEFAULT_BASIS_ACTIVITY_WINDOW,
    };
    let horizon = HorizonMode::Finite {
        num_stages: n_stages,
    };
    let initial_state = vec![50.0_f64];

    let solution = fixed_solution(100.0, 30.0);
    let solver = MockSolver::always_ok(solution);
    let comm = StubComm { rank: 0, size: 1 };
    let entity_counts = entity_counts_1_hydro();
    let (tx, _rx) = mpsc::sync_channel(8);
    let hprod = hydro_productivities_1hydro(n_stages);
    let ec = zero_energy_conversion(1, n_stages);

    // Baked templates with wrong length (2 instead of 3).
    let wrong_baked: Vec<StageTemplate> =
        (0..n_stages - 1).map(|_| minimal_template_1_0()).collect();

    let mut workspaces = single_workspace(solver);
    let result = cobre_sddp::simulate(
        &mut workspaces,
        &StageContext {
            templates: &templates,
            base_rows: &base_rows,
            noise_scale: &[],
            n_hydros: 0,
            n_load_buses: 0,
            load_balance_row_starts: &[],
            load_bus_indices: &[],
            block_counts_per_stage: &[],
            ncs_max_gen: &[],
            discount_factors: &[],
            cumulative_discount_factors: &[],
            stage_lag_transitions: &[],
            noise_group_ids: &[],
            downstream_par_order: 0,
        },
        &fcf,
        &TrainingContext {
            horizon: &horizon,
            indexer: &indexer,
            inflow_method: &InflowNonNegativityMethod::None,
            stochastic: &stochastic,
            initial_state: &initial_state,
            inflow_scheme: SamplingScheme::InSample,
            load_scheme: SamplingScheme::InSample,
            ncs_scheme: SamplingScheme::InSample,
            stages: &[],
            historical_library: None,
            external_inflow_library: None,
            external_load_library: None,
            external_ncs_library: None,
            recent_accum_seed: &[],
            recent_weight_seed: 0.0,
        },
        &config,
        SimulationOutputSpec {
            result_tx: &tx,
            zeta_per_stage: &[],
            block_hours_per_stage: &[],
            entity_counts: &entity_counts,
            generic_constraint_row_entries: &[],
            ncs_col_starts: &[],
            n_ncs_per_stage: &[],
            ncs_entity_ids_per_stage: &[],
            diversion_upstream: &HashMap::new(),
            hydro_productivities_per_stage: &hprod,
            energy_conversion: &ec,
            hydro_min_storage_hm3: &[0.0],
            event_sender: None,
        },
        Some(wrong_baked.as_slice()),
        &[],
        &comm,
    );

    match &result {
        Err(SimulationError::InvalidConfiguration(msg)) => {
            assert!(
                msg.contains('2') && msg.contains('3'),
                "error message must contain both lengths (2 and 3), got: {msg}"
            );
        }
        other => panic!("expected InvalidConfiguration error, got: {other:?}"),
    }
}

// ── Warm-start CapturedBasis acceptance tests ─────────────────

/// When `stage_bases` contains a
/// `CapturedBasis` with known slots, the basis passed to `solve(Some(&basis))`
/// has `row_status.len() == base_row_count + active_cuts_count` and the tail
/// entries match the stored cut statuses verbatim (preservation path).
///
/// Setup:
/// - `CapturedBasis` at stage 0: `base_row_count=2`, `cut_row_slots=[10,11,12]`,
///   `state_at_capture=[1.0]`, `row_status` length 5 with distinct sentinel values.
/// - FCF pool at stage 0 has exactly those 3 slots active (slots 10, 11, 12).
/// - `MockSolver::recorded_basis` captures the last basis received by
///   `solve(Some(&basis))`.
///
/// The reconstruction maps all 3 stored cut rows (slots 10, 11, 12) into
/// the output basis — the preservation path.  The tail of the output
/// `row_status` must equal the tail of the stored `row_status`.
#[test]
fn simulate_with_captured_basis_preserves_row_statuses() {
    // Sentinel status values for the stored cut rows so we can verify exact
    // preservation.  These are arbitrary non-zero i32 values; the test only
    // checks that they are passed through unchanged.
    const CUT_STATUS_0: i32 = 7;
    const CUT_STATUS_1: i32 = 11;
    const CUT_STATUS_2: i32 = 13;
    // Base rows use HIGHS_BASIS_STATUS_BASIC = 1.
    const BASE_STATUS: i32 = 1;

    let n_stages = 1;
    let n_scenarios = 1u32;
    let templates: Vec<StageTemplate> = vec![minimal_template_1_0()];
    let base_rows: Vec<usize> = vec![2]; // 2 structural rows in the template

    let indexer = StageIndexer::new(1, 0);

    // Build an FCF with 3 active cuts at slots 10, 11, 12 for stage 0.
    // warm_start_count=10, forward_passes=1 →
    //   add_cut(iter=0, fwd=0) → slot 10
    //   add_cut(iter=1, fwd=0) → slot 11
    //   add_cut(iter=2, fwd=0) → slot 12
    let mut fcf = FutureCostFunction::new(n_stages, 1, 1, 5, &[10]);
    fcf.pools[0].add_cut(0, 0, 50.0, &[1.0]);
    fcf.pools[0].add_cut(1, 0, 60.0, &[1.0]);
    fcf.pools[0].add_cut(2, 0, 70.0, &[1.0]);
    assert_eq!(
        fcf.pools[0].active_count(),
        3,
        "pool must have exactly 3 active cuts at slots 10, 11, 12"
    );
    assert_eq!(
        fcf.pools[0].populated_count, 13,
        "populated_count must be 13 (slot 12 + 1)"
    );

    // Build the CapturedBasis with matching slot metadata.
    let mut cb = CapturedBasis::new(4, 5, 2, 3, 1);
    // row_status: 2 base rows + 3 cut rows with sentinel values.
    cb.basis.row_status = vec![
        BASE_STATUS,
        BASE_STATUS,
        CUT_STATUS_0,
        CUT_STATUS_1,
        CUT_STATUS_2,
    ];
    cb.basis.col_status = vec![1_i32; 4];
    cb.cut_row_slots.extend_from_slice(&[10u32, 11, 12]);
    cb.state_at_capture.push(1.0);

    let stage_bases: Vec<Option<CapturedBasis>> = vec![Some(cb)];

    let stochastic = make_stochastic_context(n_stages);
    let config = SimulationConfig {
        n_scenarios,
        io_channel_capacity: 8,
        basis_activity_window: cobre_sddp::DEFAULT_BASIS_ACTIVITY_WINDOW,
    };
    let horizon = HorizonMode::Finite {
        num_stages: n_stages,
    };
    let initial_state = vec![50.0_f64];

    let solution = fixed_solution(100.0, 30.0);
    let solver = MockSolver::always_ok(solution);
    let comm = StubComm { rank: 0, size: 1 };
    let entity_counts = entity_counts_1_hydro();
    let (tx, _rx) = mpsc::sync_channel(16);
    let hprod = hydro_productivities_1hydro(n_stages);
    let ec = zero_energy_conversion(1, n_stages);

    let mut workspaces = single_workspace(solver);
    let result = cobre_sddp::simulate(
        &mut workspaces,
        &StageContext {
            templates: &templates,
            base_rows: &base_rows,
            noise_scale: &[],
            n_hydros: 0,
            n_load_buses: 0,
            load_balance_row_starts: &[],
            load_bus_indices: &[],
            block_counts_per_stage: &[],
            ncs_max_gen: &[],
            discount_factors: &[],
            cumulative_discount_factors: &[],
            stage_lag_transitions: &[],
            noise_group_ids: &[],
            downstream_par_order: 0,
        },
        &fcf,
        &TrainingContext {
            horizon: &horizon,
            indexer: &indexer,
            inflow_method: &InflowNonNegativityMethod::None,
            stochastic: &stochastic,
            initial_state: &initial_state,
            inflow_scheme: SamplingScheme::InSample,
            load_scheme: SamplingScheme::InSample,
            ncs_scheme: SamplingScheme::InSample,
            stages: &[],
            historical_library: None,
            external_inflow_library: None,
            external_load_library: None,
            external_ncs_library: None,
            recent_accum_seed: &[],
            recent_weight_seed: 0.0,
        },
        &config,
        SimulationOutputSpec {
            result_tx: &tx,
            zeta_per_stage: &[],
            block_hours_per_stage: &[],
            entity_counts: &entity_counts,
            generic_constraint_row_entries: &[],
            ncs_col_starts: &[],
            n_ncs_per_stage: &[],
            ncs_entity_ids_per_stage: &[],
            diversion_upstream: &HashMap::new(),
            hydro_productivities_per_stage: &hprod,
            energy_conversion: &ec,
            hydro_min_storage_hm3: &[0.0],
            event_sender: None,
        },
        // fallback path (no baked templates); reconstruction uses pool.active_cuts()
        None,
        &stage_bases,
        &comm,
    );

    assert!(
        result.is_ok(),
        "simulate must succeed with CapturedBasis warm-start: {result:?}"
    );

    // Verify that solve(Some(&basis)) was called (warm-start path taken).
    let solver = &workspaces[0].solver;
    assert_eq!(
        solver.solve_with_basis_count, 1,
        "warm-start solve must be called exactly once (1 scenario × 1 stage)"
    );
    assert_eq!(
        solver.solve_count, 0,
        "cold-start solve must not be called when a CapturedBasis is provided"
    );

    // Verify the basis passed to solve(Some(&basis)) has the correct structure.
    let recorded = solver
        .recorded_basis
        .as_ref()
        .expect("recorded_basis must be Some after a warm-start solve");

    // Row count: base_row_count=2 + 3 active cut slots = 5.
    assert_eq!(
        recorded.row_status.len(),
        2 + 3,
        "reconstructed basis row_status must have length base_row_count(2) + cuts(3) = 5, \
         got {}",
        recorded.row_status.len()
    );

    // The last 3 entries must match the stored cut statuses verbatim
    // (all 3 slots were in the stored basis → preservation path).
    let tail = &recorded.row_status[2..];
    assert_eq!(
        tail,
        &[CUT_STATUS_0, CUT_STATUS_1, CUT_STATUS_2],
        "reconstructed basis tail must match stored cut statuses verbatim \
         (preservation path): expected [{CUT_STATUS_0}, {CUT_STATUS_1}, {CUT_STATUS_2}], \
         got {tail:?}"
    );
}

/// When `stage_bases` is `&[]`
/// (cold-start), every LP solve must go through `solver.solve(None)` and
/// `solve(Some(&basis))` must never be called.
///
/// Uses `solve_count` and `solve_with_basis_count` split on `MockSolver`.
#[test]
fn simulate_with_empty_stage_bases_cold_starts() {
    let n_stages = 2;
    let n_scenarios = 3u32;
    let templates: Vec<StageTemplate> = (0..n_stages).map(|_| minimal_template_1_0()).collect();
    let base_rows: Vec<usize> = vec![0; n_stages];

    let indexer = StageIndexer::new(1, 0);
    let fcf = FutureCostFunction::new(n_stages, 1, 1, 10, &vec![0; n_stages]);
    let stochastic = make_stochastic_context(n_stages);
    let config = SimulationConfig {
        n_scenarios,
        io_channel_capacity: 16,
        basis_activity_window: cobre_sddp::DEFAULT_BASIS_ACTIVITY_WINDOW,
    };
    let horizon = HorizonMode::Finite {
        num_stages: n_stages,
    };
    let initial_state = vec![50.0_f64];

    let solution = fixed_solution(100.0, 30.0);
    let solver = MockSolver::always_ok(solution);
    let comm = StubComm { rank: 0, size: 1 };
    let entity_counts = entity_counts_1_hydro();
    let (tx, _rx) = mpsc::sync_channel(32);
    let hprod = hydro_productivities_1hydro(n_stages);
    let ec = zero_energy_conversion(1, n_stages);

    let mut workspaces = single_workspace(solver);
    let result = cobre_sddp::simulate(
        &mut workspaces,
        &StageContext {
            templates: &templates,
            base_rows: &base_rows,
            noise_scale: &[],
            n_hydros: 0,
            n_load_buses: 0,
            load_balance_row_starts: &[],
            load_bus_indices: &[],
            block_counts_per_stage: &[],
            ncs_max_gen: &[],
            discount_factors: &[],
            cumulative_discount_factors: &[],
            stage_lag_transitions: &[],
            noise_group_ids: &[],
            downstream_par_order: 0,
        },
        &fcf,
        &TrainingContext {
            horizon: &horizon,
            indexer: &indexer,
            inflow_method: &InflowNonNegativityMethod::None,
            stochastic: &stochastic,
            initial_state: &initial_state,
            inflow_scheme: SamplingScheme::InSample,
            load_scheme: SamplingScheme::InSample,
            ncs_scheme: SamplingScheme::InSample,
            stages: &[],
            historical_library: None,
            external_inflow_library: None,
            external_load_library: None,
            external_ncs_library: None,
            recent_accum_seed: &[],
            recent_weight_seed: 0.0,
        },
        &config,
        SimulationOutputSpec {
            result_tx: &tx,
            zeta_per_stage: &[],
            block_hours_per_stage: &[],
            entity_counts: &entity_counts,
            generic_constraint_row_entries: &[],
            ncs_col_starts: &[],
            n_ncs_per_stage: &[],
            ncs_entity_ids_per_stage: &[],
            diversion_upstream: &HashMap::new(),
            hydro_productivities_per_stage: &hprod,
            energy_conversion: &ec,
            hydro_min_storage_hm3: &[0.0],
            event_sender: None,
        },
        None,
        // empty stage_bases → cold-start for every stage
        &[],
        &comm,
    );

    assert!(
        result.is_ok(),
        "cold-start simulate must succeed: {result:?}"
    );

    let solver = &workspaces[0].solver;
    let expected_solves = n_scenarios as usize * n_stages;

    assert_eq!(
        solver.solve_with_basis_count, 0,
        "warm-start solve must not be called when stage_bases is empty; \
         got solve_with_basis_count={}",
        solver.solve_with_basis_count
    );
    assert_eq!(
        solver.solve_count, expected_solves,
        "cold-start solve must be called exactly n_scenarios({n_scenarios}) × \
         n_stages({n_stages}) = {expected_solves} times; got {}",
        solver.solve_count
    );
}
