//! `Historical` scenario sampling scheme — library type and eta pre-standardization.
//!
//! [`HistoricalScenarioLibrary`] stores pre-standardized eta values for
//! historical inflow windows. During the forward pass, the
//! `ClassSampler::Historical` variant indexes into this library to retrieve
//! noise vectors for a given (window, stage) pair.
//!
//! The [`standardize_historical_windows`] function populates the library by
//! inverting the PAR(p) model: for every valid window and study stage it
//! computes `η = (obs - deterministic_base - Σ ψ[ℓ]·lag[ℓ]) / σ` using a
//! **rolling lag chain seeded from `initial_conditions.past_inflows`** (not
//! the window's own pre-study observations). This guarantees that, starting
//! from `past_inflows` at stage 0, the PAR forward pass exactly reconstructs
//! the raw historical target at every stage of the replay.
//!
//! ## Replay correctness
//!
//! With `x₀ = past_inflows` the lag state at stage 0, the η stored for window
//! `w` at stage `t` satisfies
//!
//! ```text
//! η(w, t) = (target(w, t) - b(t) - Σ ψ(t)[ℓ] · lag(w, t)[ℓ]) / σ(t)
//! ```
//!
//! where `lag(w, t)` advances `x₀` through window `w`'s raw historical targets
//! at stages `0..t` — exactly what the forward accumulator produces from `x₀`.
//! Replay is therefore exact only if the forward pass starts from `past_inflows`
//! verbatim; `ClassSampler::Historical::apply_initial_state` is a no-op, so that
//! is the caller's responsibility, not this module's.
//!
//! The `eta` buffer uses **window-major** layout
//! (`eta[window * n_stages * n_hydros + stage * n_hydros + hydro]`), matching
//! sequential stage iteration within a window (same access pattern as
//! [`PrecomputedPar`]).
//!
//! Initial inflow lags are **NOT** stored on the library nor injected into the
//! solver state vector. Per the hydrological-tendency convention they come from
//! `initial_conditions.past_inflows` for every scenario regardless of the window
//! replayed; η inversion uses that `past_inflows`-seeded rolling chain.
//!
//! [`PrecomputedPar`]: crate::par::precompute::PrecomputedPar

use chrono::{Datelike, NaiveDate};
use cobre_core::{
    EntityId, HydroPastInflows,
    scenario::{HistoricalYears, InflowHistoryRow},
    temporal::{SeasonMap, Stage, StageLagTransition},
};
use siphasher::sip::SipHasher13;
use std::hash::Hasher as _;

use crate::{
    StochasticError,
    par::{evaluate::solve_par_noise, fitting::find_season_for_date, precompute::PrecomputedPar},
};

// ---------------------------------------------------------------------------
// HistoricalScenarioLibrary
// ---------------------------------------------------------------------------

/// Pre-standardized eta store for historical scenario windows.
///
/// A pure data container — no sampling logic (permutation, selection) is
/// included. Population is performed by the eta-standardisation pass after
/// construction; selection is performed by the `ClassSampler::Historical`
/// variant during the forward pass.
///
/// # Construction
///
/// Use [`HistoricalScenarioLibrary::new`], which allocates zero-filled
/// buffers.
///
/// # Examples
///
/// ```
/// use cobre_stochastic::HistoricalScenarioLibrary;
///
/// let mut lib = HistoricalScenarioLibrary::new(3, 12, 5, 2, vec![1990, 1995, 2000]);
/// assert_eq!(lib.n_windows(), 3);
/// assert_eq!(lib.n_stages(), 12);
/// assert_eq!(lib.n_hydros(), 5);
/// assert_eq!(lib.max_order(), 2);
/// assert_eq!(lib.window_year(0), 1990);
/// assert_eq!(lib.window_year(2), 2000);
///
/// // Write and read eta values.
/// lib.eta_slice_mut(1, 3).copy_from_slice(&[0.1, 0.2, 0.3, 0.4, 0.5]);
/// assert_eq!(lib.eta_slice(1, 3), &[0.1, 0.2, 0.3, 0.4, 0.5]);
/// ```
#[derive(Debug, Clone)]
pub struct HistoricalScenarioLibrary {
    eta: Box<[f64]>,
    window_years: Box<[i32]>,
    n_windows: usize,
    n_stages: usize,
    n_hydros: usize,
    max_order: usize,
    /// SipHash-1-3 fingerprint of the seeding `past_inflows`; `0` until
    /// [`standardize_historical_windows`] writes it. A mismatch against a fresh
    /// digest means the library is stale and must be re-standardised.
    past_inflows_digest: u64,
}

impl HistoricalScenarioLibrary {
    /// Construct a new library with zero-filled buffers.
    ///
    /// # Parameters
    ///
    /// - `max_order` — PAR model order (number of pre-window lag stages)
    /// - `window_years` — starting year per window (length must equal `n_windows`)
    ///
    /// # Panics
    ///
    /// Panics in debug builds if `window_years.len() != n_windows`.
    #[must_use]
    pub fn new(
        n_windows: usize,
        n_stages: usize,
        n_hydros: usize,
        max_order: usize,
        window_years: Vec<i32>,
    ) -> Self {
        debug_assert_eq!(
            window_years.len(),
            n_windows,
            "window_years length ({}) must equal n_windows ({})",
            window_years.len(),
            n_windows,
        );
        Self {
            eta: vec![0.0_f64; n_windows * n_stages * n_hydros].into_boxed_slice(),
            window_years: window_years.into_boxed_slice(),
            n_windows,
            n_stages,
            n_hydros,
            max_order,
            past_inflows_digest: 0,
        }
    }

    // -----------------------------------------------------------------------
    // Dimension accessors
    // -----------------------------------------------------------------------

    /// Returns the number of historical windows.
    #[must_use]
    #[inline]
    pub fn n_windows(&self) -> usize {
        self.n_windows
    }

    /// Returns the number of study stages per window.
    #[must_use]
    #[inline]
    pub fn n_stages(&self) -> usize {
        self.n_stages
    }

    /// Returns the number of hydro entities (eta vector width).
    #[must_use]
    #[inline]
    pub fn n_hydros(&self) -> usize {
        self.n_hydros
    }

    /// Returns the PAR model order (number of pre-window lag stages stored).
    #[must_use]
    #[inline]
    pub fn max_order(&self) -> usize {
        self.max_order
    }

    /// Returns the SipHash-1-3 fingerprint of the seeding `past_inflows`, or `0`
    /// before [`standardize_historical_windows`] has run.
    #[must_use]
    #[inline]
    pub fn past_inflows_digest(&self) -> u64 {
        self.past_inflows_digest
    }

    /// Returns the starting year label for window `window`.
    #[must_use]
    #[inline]
    pub fn window_year(&self, window: usize) -> i32 {
        debug_assert!(
            window < self.n_windows,
            "window ({window}) must be < n_windows ({})",
            self.n_windows
        );
        self.window_years[window]
    }

    // -----------------------------------------------------------------------
    // Eta accessors
    // -----------------------------------------------------------------------

    /// Returns the `n_hydros`-length slice of eta values for `(window, stage)`.
    ///
    /// Layout: `eta[window * n_stages * n_hydros + stage * n_hydros + hydro]`.
    ///
    /// # Panics
    ///
    /// Panics if `window >= n_windows` or `stage >= n_stages`.
    #[must_use]
    #[inline]
    pub fn eta_slice(&self, window: usize, stage: usize) -> &[f64] {
        assert!(
            window < self.n_windows,
            "window ({window}) must be < n_windows ({})",
            self.n_windows
        );
        assert!(
            stage < self.n_stages,
            "stage ({stage}) must be < n_stages ({})",
            self.n_stages
        );
        let offset = (window * self.n_stages + stage) * self.n_hydros;
        &self.eta[offset..offset + self.n_hydros]
    }

    /// Returns a mutable `n_hydros`-length slice of eta values for `(window, stage)`.
    ///
    /// # Panics
    ///
    /// Panics if `window >= n_windows` or `stage >= n_stages`.
    #[must_use]
    #[inline]
    pub fn eta_slice_mut(&mut self, window: usize, stage: usize) -> &mut [f64] {
        assert!(
            window < self.n_windows,
            "window ({window}) must be < n_windows ({})",
            self.n_windows
        );
        assert!(
            stage < self.n_stages,
            "stage ({stage}) must be < n_stages ({})",
            self.n_stages
        );
        let offset = (window * self.n_stages + stage) * self.n_hydros;
        &mut self.eta[offset..offset + self.n_hydros]
    }
}

// ---------------------------------------------------------------------------
// standardize_historical_windows
// ---------------------------------------------------------------------------

/// Populate a [`HistoricalScenarioLibrary`] with pre-standardized η values whose
/// lag chain is seeded from `past_inflows` and advanced by the same
/// accumulate/finalize pattern as `standardize_external_inflow`.
///
/// Replay is exact only if the forward pass starts from the same `past_inflows`;
/// see the module docs for the inductive argument.
/// `ClassSampler::Historical::apply_initial_state` is a no-op, so that is the
/// caller's responsibility. The η values depend on `past_inflows`; a change
/// requires re-standardising, which `past_inflows_digest` lets callers detect.
///
/// Only uniform finalising transitions (`accumulate_weight == 1.0`,
/// `spillover_weight == 0.0`, `finalize_period == true` — all monthly studies)
/// are supported; a `debug_assert!` rejects any other pattern, and release builds
/// treat it as best-effort. See `TODO(historical-replay-non-monthly)` below.
///
/// # Inputs
///
/// - `hydro_ids` — canonical-order hydro entity IDs (must match `par`)
/// - `season_map` — three-tier mapping of observation dates to season IDs (same as
///   [`discover_historical_windows`](super::window::discover_historical_windows)):
///   1. dates within a stage's `[start_date, end_date)` → that stage's season;
///   2. otherwise `season_map.season_for_date(date)` when `Some`;
///   3. when `None`, `month0()`. Unmappable observations are skipped.
/// - `past_inflows` — pre-study history sorted by `hydro_id`, seeding the stage-0
///   chain; absent hydros default to 0.0 for all lag slots.
/// - `stage_lag_transitions` — one per stage, same length as `stages`.
///
/// # Panics
///
/// Panics in debug builds if dimension mismatches between `library`, `par`,
/// `stages`, or `stage_lag_transitions` are detected.
// Rationale: the lag-state and accumulator buffers, the observation table, and the
// digest computation thread through the (window × stage × entity) pass in
// sequence; extracting helpers would pass several mutable buffers across boundaries.
#[allow(clippy::too_many_lines)]
pub fn standardize_historical_windows(
    library: &mut HistoricalScenarioLibrary,
    inflow_history: &[InflowHistoryRow],
    hydro_ids: &[EntityId],
    stages: &[Stage],
    par: &PrecomputedPar,
    window_years: &[i32],
    season_map: Option<&SeasonMap>,
    past_inflows: &[HydroPastInflows],
    stage_lag_transitions: &[StageLagTransition],
) {
    // Fallback for an empty `stage_lag_transitions` or a per-stage noop
    // (accumulate_weight≈0, finalize_period=false) — what an empty SeasonMap
    // emits for monthly stages.
    const UNIFORM_MONTHLY: StageLagTransition = StageLagTransition {
        accumulate_weight: 1.0,
        spillover_weight: 0.0,
        finalize_period: true,
        accumulate_downstream: false,
        downstream_accumulate_weight: 0.0,
        downstream_spillover_weight: 0.0,
        downstream_finalize: false,
        rebuild_from_downstream: false,
    };

    debug_assert_eq!(
        library.n_windows(),
        window_years.len(),
        "library.n_windows() ({}) must equal window_years.len() ({})",
        library.n_windows(),
        window_years.len(),
    );
    debug_assert_eq!(
        library.n_stages(),
        stages.len(),
        "library.n_stages() ({}) must equal stages.len() ({})",
        library.n_stages(),
        stages.len(),
    );
    debug_assert_eq!(
        library.n_hydros(),
        hydro_ids.len(),
        "library.n_hydros() ({}) must equal hydro_ids.len() ({})",
        library.n_hydros(),
        hydro_ids.len(),
    );
    debug_assert_eq!(
        library.max_order(),
        par.max_order(),
        "library.max_order() ({}) must equal par.max_order() ({})",
        library.max_order(),
        par.max_order(),
    );
    debug_assert!(
        stage_lag_transitions.is_empty() || stage_lag_transitions.len() == stages.len(),
        "stage_lag_transitions.len() ({}) must be 0 or equal to stages.len() ({})",
        stage_lag_transitions.len(),
        stages.len(),
    );

    // Only uniform-finalising or noop (empty-SeasonMap) transitions are handled;
    // any other pattern is a non-monthly grid this scheme does not yet support.
    // TODO(historical-replay-non-monthly): relax when historical windows are
    // needed for weekly/sub-monthly study grids.
    debug_assert!(
        stage_lag_transitions.iter().all(|t| {
            let is_uniform_monthly = (t.accumulate_weight - 1.0).abs() < 1e-9
                && t.spillover_weight.abs() < 1e-9
                && t.finalize_period;
            let is_noop = t.accumulate_weight.abs() < 1e-9
                && t.spillover_weight.abs() < 1e-9
                && !t.finalize_period;
            is_uniform_monthly || is_noop
        }),
        "standardize_historical_windows: non-trivial stage_lag_transitions detected \
         (non-monthly study). Historical replay under non-uniform lag periods is \
         not fully validated — see TODO(historical-replay-non-monthly).",
    );

    let n_hydros = library.n_hydros();
    let n_stages = library.n_stages();
    let max_order = library.max_order();

    if n_hydros == 0 || n_stages == 0 || window_years.is_empty() {
        return;
    }

    let n_seasons = stages
        .iter()
        .filter_map(|s| s.season_id)
        .max()
        .map_or(1, |m| m + 1);

    let hydro_id_to_idx: std::collections::HashMap<EntityId, usize> = hydro_ids
        .iter()
        .enumerate()
        .map(|(i, &id)| (id, i))
        .collect();

    let (Some(min_year), Some(max_year)) = (
        inflow_history.iter().map(|r| r.date.year()).min(),
        inflow_history.iter().map(|r| r.date.year()).max(),
    ) else {
        return;
    };
    #[allow(clippy::cast_sign_loss)]
    let n_years = (max_year - min_year + 1) as usize;

    let table_size = n_hydros * n_years * n_seasons;
    let mut obs_table = vec![f64::NAN; table_size];

    let table_idx = |h: usize, year: i32, s: usize| -> Option<usize> {
        if year < min_year || year > max_year || s >= n_seasons {
            return None;
        }
        #[allow(clippy::cast_sign_loss)]
        let y = (year - min_year) as usize;
        Some(h * n_years * n_seasons + y * n_seasons + s)
    };

    let mut stage_index: Vec<(NaiveDate, NaiveDate, i32, usize)> = stages
        .iter()
        .filter_map(|s| s.season_id.map(|sid| (s.start_date, s.end_date, s.id, sid)))
        .collect();
    stage_index.sort_unstable_by_key(|(start, _, _, _)| *start);

    for r in inflow_history {
        let season_id = find_season_for_date(&stage_index, r.date)
            .or_else(|| season_map.and_then(|sm| sm.season_for_date(r.date)))
            .or_else(|| {
                if season_map.is_none() {
                    Some(r.date.month0() as usize)
                } else {
                    None
                }
            });
        if let Some(sid) = season_id
            && let Some(&h) = hydro_id_to_idx.get(&r.hydro_id)
            && let Some(idx) = table_idx(h, r.date.year(), sid)
        {
            obs_table[idx] = r.value_m3s;
        }
    }

    let lookup = |h: usize, year: i32, season_id: usize| -> f64 {
        table_idx(h, year, season_id).map_or(0.0, |idx| {
            let v = obs_table[idx];
            if v.is_nan() { 0.0 } else { v }
        })
    };

    let full_sequence: Vec<(i32, usize)> =
        super::build_observation_sequence(stages, max_order, n_seasons);

    // past_lag_buf[h_idx * safe_max_order + lag]
    let safe_max_order = max_order.max(1);
    let mut past_lag_buf = vec![0.0_f64; n_hydros * safe_max_order];
    for pi in past_inflows {
        if let Some(&h_idx) = hydro_id_to_idx.get(&pi.hydro_id) {
            // Fill up to `max_order` slots so PAR(p)-A annual contributions
            // (widened across the `psi` slice) see real lag values.
            for (lag, &value) in pi.values_m3s.iter().enumerate().take(max_order) {
                past_lag_buf[h_idx * safe_max_order + lag] = value;
            }
        }
    }

    // Digest over little-endian f64 bytes so it is reproducible across runs.
    {
        let mut hasher = SipHasher13::new();
        for pi in past_inflows {
            for &v in &pi.values_m3s {
                hasher.write(&v.to_le_bytes());
            }
        }
        library.past_inflows_digest = hasher.finish();
    }

    // lag_state[h * safe_max_order + l]
    let mut lag_state = vec![0.0_f64; n_hydros * safe_max_order];
    let mut lag_buf = vec![0.0_f64; safe_max_order];
    let mut lag_accum = vec![0.0_f64; n_hydros];

    for (w, &window_year) in window_years.iter().enumerate() {
        // Each window starts from the same past_inflows-seeded lag state.
        lag_state.copy_from_slice(&past_lag_buf);
        lag_accum.fill(0.0);
        let mut lag_weight_accum = 0.0_f64;

        for t in 0..n_stages {
            let (year_offset, season_id) = full_sequence[max_order + t];
            let obs_year = window_year + year_offset;
            debug_assert!(
                table_idx(0, obs_year, season_id).is_some_and(|i| {
                    _ = i;
                    true
                }),
                "missing study observation for year={obs_year}, season={season_id}; \
                 window discovery should have excluded this window",
            );

            let eta_slice = library.eta_slice_mut(w, t);
            for h in 0..n_hydros {
                debug_assert!(
                    max_order == 0
                        || table_idx(h, obs_year, season_id)
                            .is_some_and(|i| !obs_table[i].is_nan()),
                    "missing study observation for hydro={}, year={obs_year}, \
                     season={season_id}; window discovery should have excluded this window",
                    hydro_ids[h].0,
                );
                let target = lookup(h, obs_year, season_id);

                // lag_state is h-major: lag_state[h * safe_max_order + l]. Fill
                // every slot so PAR(p)-A annual contributions (widened across the
                // `psi` slice) participate in the η inversion.
                for (l, slot) in lag_buf.iter_mut().enumerate() {
                    *slot = lag_state[h * safe_max_order + l];
                }

                let det_base = par.deterministic_base(t, h);
                let psi = par.psi_slice(t, h);
                let sigma = par.sigma(t, h);

                let eta = solve_par_noise(det_base, psi, &lag_buf, sigma, target);
                eta_slice[h] = eta;
            }

            // Advance the lag chain, mirroring standardize_external_inflow.
            let stage_lag = if stage_lag_transitions.is_empty() {
                &UNIFORM_MONTHLY
            } else {
                let t_lag = &stage_lag_transitions[t];
                // A noop (empty-SeasonMap) entry is treated as uniform-monthly.
                if t_lag.accumulate_weight.abs() < 1e-9 && !t_lag.finalize_period {
                    &UNIFORM_MONTHLY
                } else {
                    t_lag
                }
            };

            let aw = stage_lag.accumulate_weight;
            for (accum, h) in lag_accum.iter_mut().zip(0..n_hydros) {
                let val = lookup(h, obs_year, season_id);
                *accum += val * aw;
            }
            lag_weight_accum += aw;

            if stage_lag.finalize_period && lag_weight_accum > 0.0 {
                let inv = 1.0 / lag_weight_accum;
                for (h, accum) in (0..n_hydros).zip(lag_accum.iter_mut()) {
                    let avg = *accum * inv;
                    // Shift highest-lag-first so each slot reads its predecessor
                    // before being overwritten.
                    for l in (1..safe_max_order).rev() {
                        lag_state[h * safe_max_order + l] = lag_state[h * safe_max_order + l - 1];
                    }
                    lag_state[h * safe_max_order] = avg;
                }

                // Spillover seeds from the RAW value, not the period average.
                let sw = stage_lag.spillover_weight;
                if sw > 0.0 {
                    for (accum, h) in lag_accum.iter_mut().zip(0..n_hydros) {
                        let raw_val = lookup(h, obs_year, season_id);
                        *accum = raw_val * sw;
                    }
                    lag_weight_accum = sw;
                } else {
                    lag_accum.fill(0.0);
                    lag_weight_accum = 0.0;
                }
            }
            // Non-finalising stages keep lag_state frozen.
        }
    }
}

// ---------------------------------------------------------------------------
// validate_historical_library
// ---------------------------------------------------------------------------

/// Validate a [`HistoricalScenarioLibrary`] against construction inputs.
///
/// The Tier 2 gate, run after window discovery and eta standardization. Fail-fast:
/// the first failed check returns `Err`; the V2.6 warning is emitted via
/// `tracing::warn!` without aborting.
///
/// ## Checks performed
///
/// | ID  | Kind    | Description                                                   |
/// |-----|---------|---------------------------------------------------------------|
/// | V2.1 | Error  | Every study stage must have `season_id: Some(_)`.             |
/// | V2.9 | Error  | `hydro_ids.len()` must equal `library.n_hydros()`.            |
/// | V2.5 | Error  | At least one window must be discovered when `user_pool` is `None`. |
/// | V2.3 | Error  | No eta value in the library may be `f64::NEG_INFINITY`.       |
/// | V2.6 | Warning| `library.n_windows() < forward_passes` — log a warning.      |
/// | V2.2 | Assert | Window contiguity — `debug_assert` only (construction invariant). |
/// | V2.4 | Assert | User pool validity — `debug_assert` only (construction invariant). |
/// | V2.7 | Assert | Lag warmup sufficiency — `debug_assert` only (construction invariant). |
///
/// # Errors
///
/// Returns [`StochasticError::InsufficientData`] with a message prefixed by
/// the check ID (e.g., `"V2.1: ..."`) for the first failed error check.
///
/// # Examples
///
/// ```
/// use cobre_core::{EntityId, scenario::InflowHistoryRow};
/// use cobre_core::temporal::{
///     Block, BlockMode, NoiseMethod, ScenarioSourceConfig, Stage, StageRiskConfig,
///     StageStateConfig,
/// };
/// use chrono::NaiveDate;
/// use cobre_stochastic::{HistoricalScenarioLibrary, sampling::historical::validate_historical_library};
///
/// let lib = HistoricalScenarioLibrary::new(3, 1, 2, 1, vec![1990, 1995, 2000]);
/// let stage = Stage {
///     index: 0,
///     id: 0,
///     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: "B".to_string(), duration_hours: 720.0 }],
///     block_mode: BlockMode::Parallel,
///     state_config: StageStateConfig { storage: true, inflow_lags: false },
///     risk_config: StageRiskConfig::Expectation,
///     scenario_config: ScenarioSourceConfig { branching_factor: 1, noise_method: NoiseMethod::Saa },
/// };
/// let hydro_ids = [EntityId(1), EntityId(2)];
/// let result = validate_historical_library(&lib, &[], &hydro_ids, &[stage], 1, None, 5);
/// assert!(result.is_ok());
/// ```
pub fn validate_historical_library(
    library: &HistoricalScenarioLibrary,
    _inflow_history: &[InflowHistoryRow],
    hydro_ids: &[EntityId],
    stages: &[Stage],
    max_par_order: usize,
    user_pool: Option<&HistoricalYears>,
    forward_passes: u32,
) -> Result<(), StochasticError> {
    // V2.1 — every study stage must have season_id: Some(_).
    for stage in stages {
        if stage.season_id.is_none() {
            return Err(StochasticError::InsufficientData {
                context: format!(
                    "V2.1: stage {} (index {}) has season_id: None; \
                     all study stages must have a season_id assigned",
                    stage.id, stage.index,
                ),
            });
        }
    }

    // V2.9 — hydro_ids length must equal library.n_hydros().
    if hydro_ids.len() != library.n_hydros() {
        return Err(StochasticError::InsufficientData {
            context: format!(
                "V2.9: hydro_ids slice length ({}) does not match \
                 library.n_hydros() ({})",
                hydro_ids.len(),
                library.n_hydros(),
            ),
        });
    }

    // V2.5 — at least 1 window when user_pool is None.
    if user_pool.is_none() && library.n_windows() == 0 {
        return Err(StochasticError::InsufficientData {
            context: "V2.5: historical library has 0 windows after auto-discovery; \
                      at least 1 complete historical window is required"
                .to_string(),
        });
    }

    // V2.3 / V2.8 — no eta may be NEG_INFINITY (the sigma=0 non-matching-observation
    // sentinel from standardize_historical_windows) or NaN.
    for w in 0..library.n_windows() {
        for t in 0..library.n_stages() {
            let eta = library.eta_slice(w, t);
            for (h, &value) in eta.iter().enumerate() {
                if value == f64::NEG_INFINITY || value.is_nan() {
                    return Err(StochasticError::InsufficientData {
                        context: format!(
                            "V2.3: historical library contains non-finite eta (NEG_INFINITY or NaN) \
                             at window {w}, stage {t}, hydro {h} — sigma=0 with \
                             non-matching historical observation or numerical failure",
                        ),
                    });
                }
            }
        }
    }

    // V2.2/V2.4/V2.7 are construction invariants; only V2.7 is re-asserted here.
    debug_assert!(
        library.max_order() == max_par_order,
        "V2.7: library.max_order() ({}) must equal max_par_order ({max_par_order})",
        library.max_order(),
    );

    // V2.6 — warn (do not fail) when windows < forward passes.
    if library.n_windows() < forward_passes as usize {
        tracing::warn!(
            n_windows = library.n_windows(),
            forward_passes = forward_passes,
            "historical library has fewer windows ({}) than forward passes ({}); \
             windows will be reused across forward passes",
            library.n_windows(),
            forward_passes,
        );
    }

    Ok(())
}

// ---------------------------------------------------------------------------
// Tests
// ---------------------------------------------------------------------------

#[cfg(test)]
#[allow(
    clippy::unwrap_used,
    clippy::expect_used,
    clippy::panic,
    clippy::float_cmp
)]
mod tests {
    use super::HistoricalScenarioLibrary;

    #[test]
    fn test_new_allocates_correct_sizes() {
        let lib = HistoricalScenarioLibrary::new(3, 12, 5, 2, vec![1990, 1995, 2000]);

        assert_eq!(lib.n_windows(), 3);
        assert_eq!(lib.n_stages(), 12);
        assert_eq!(lib.n_hydros(), 5);
        assert_eq!(lib.max_order(), 2);

        assert_eq!(
            lib.eta_slice(0, 0).len(),
            5,
            "eta slice length must equal n_hydros"
        );
        assert_eq!(
            lib.eta_slice(2, 11).len(),
            5,
            "eta slice length must equal n_hydros"
        );
    }

    #[test]
    fn test_eta_roundtrip() {
        let mut lib = HistoricalScenarioLibrary::new(2, 3, 4, 1, vec![2000, 2001]);

        let values = [1.0_f64, 2.0, 3.0, 4.0];
        lib.eta_slice_mut(1, 2).copy_from_slice(&values);

        assert_eq!(
            lib.eta_slice(1, 2),
            &values,
            "eta_slice must return the values written via eta_slice_mut"
        );

        // Confirm that other (window, stage) cells were not disturbed.
        assert_eq!(
            lib.eta_slice(0, 0),
            &[0.0, 0.0, 0.0, 0.0],
            "untouched eta cells must remain zero"
        );
        assert_eq!(
            lib.eta_slice(1, 0),
            &[0.0, 0.0, 0.0, 0.0],
            "untouched eta cells must remain zero"
        );
    }

    #[test]
    fn test_window_years() {
        let lib = HistoricalScenarioLibrary::new(3, 1, 1, 1, vec![1990, 1995, 2000]);

        assert_eq!(lib.window_year(0), 1990);
        assert_eq!(lib.window_year(1), 1995);
        assert_eq!(lib.window_year(2), 2000);
    }

    #[test]
    fn test_send_sync() {
        fn assert_send_sync<T: Send + Sync>() {}
        assert_send_sync::<HistoricalScenarioLibrary>();
    }

    // -----------------------------------------------------------------------
    // Helpers for standardize_historical_windows tests
    // -----------------------------------------------------------------------

    use chrono::{Datelike, NaiveDate};
    use cobre_core::{
        EntityId, HydroPastInflows,
        scenario::{InflowHistoryRow, InflowModel},
        temporal::{
            Block, BlockMode, NoiseMethod, ScenarioSourceConfig, SeasonCycleType, SeasonDefinition,
            SeasonMap, StageRiskConfig, StageStateConfig,
        },
    };

    use super::{Stage, standardize_historical_windows};
    use crate::par::precompute::PrecomputedPar;

    /// Build a monthly stage with the given array index and 0-based `season_id` (0=Jan..11=Dec).
    ///
    /// Uses a 12-season cycle so that wrap-around (Dec→Jan) causes a year-offset
    /// increment in `build_observation_sequence`, matching the window.rs tests.
    #[allow(clippy::cast_possible_truncation, clippy::cast_possible_wrap)]
    fn make_monthly_stage(index: usize, season_id: usize) -> Stage {
        let month = (season_id as u32) + 1;
        Stage {
            index,
            id: index as i32,
            start_date: NaiveDate::from_ymd_opt(2024, month, 1).unwrap(),
            end_date: NaiveDate::from_ymd_opt(2024, month, 28).unwrap(),
            season_id: Some(season_id),
            blocks: vec![Block {
                index: 0,
                name: "SINGLE".to_string(),
                duration_hours: 720.0,
            }],
            block_mode: BlockMode::Parallel,
            state_config: StageStateConfig {
                storage: true,
                inflow_lags: false,
            },
            risk_config: StageRiskConfig::Expectation,
            scenario_config: ScenarioSourceConfig {
                branching_factor: 1,
                noise_method: NoiseMethod::Saa,
            },
        }
    }

    /// Build a row keyed by `month0` (0-based month: 0=Jan, 11=Dec).
    ///
    /// The lookup in `standardize_historical_windows` uses `date.month0()` as
    /// the `season_id`. This helper makes that mapping explicit.
    fn make_row(hydro_id: EntityId, year: i32, month0: u32, value: f64) -> InflowHistoryRow {
        InflowHistoryRow {
            hydro_id,
            date: NaiveDate::from_ymd_opt(year, month0 + 1, 1).unwrap(),
            value_m3s: value,
        }
    }

    /// Build 12 monthly study stages covering a full calendar year (seasons 0-11).
    ///
    /// With `max_order=0` and these 12 stages, the year-offset increments once
    /// at the Jan→Feb transition when the season wraps from 11→0, placing study
    /// stages at year `window_year + 1`.
    fn twelve_monthly_stages() -> Vec<Stage> {
        (0..12).map(|i| make_monthly_stage(i, i)).collect()
    }

    // -----------------------------------------------------------------------
    // Test 1: AR(0) single hydro standardization
    // -----------------------------------------------------------------------

    /// Given 1 hydro with AR(0), mean=100, std=30, 12 study stages (full year),
    /// and window year 1990 with observations [120.0, 90.0, ...], the eta values
    /// at stage 0 and 1 must be (120-100)/30 and (90-100)/30.
    ///
    /// With 12 monthly stages (`season_ids` 0-11), `n_seasons`=12. `max_order`=0.
    /// Full sequence: [(0,0),(0,1),...,(0,11)] — all `year_offset`=0 because the
    /// sequence 0→1→...→11 never wraps backwards. So all study observations
    /// are in year `window_year` itself.
    #[test]
    fn test_ar0_standardization() {
        let hydro = EntityId(1);
        // Use 2 study stages with season_ids 0 and 1 (Jan and Feb).
        // n_seasons = max(0,1)+1 = 2. Full sequence with max_order=0: [(0,0),(0,1)].
        // Year offsets: 0→1 is increasing (no wrap), both stay at offset=0.
        // Observations are at year=window_year=1990.
        let stages = vec![make_monthly_stage(0, 0), make_monthly_stage(1, 1)];
        let models = vec![
            InflowModel {
                hydro_id: hydro,
                stage_id: 0,
                mean_m3s: 100.0,
                std_m3s: 30.0,
                ar_coefficients: vec![],
                residual_std_ratio: 1.0,
                annual: None,
            },
            InflowModel {
                hydro_id: hydro,
                stage_id: 1,
                mean_m3s: 100.0,
                std_m3s: 30.0,
                ar_coefficients: vec![],
                residual_std_ratio: 1.0,
                annual: None,
            },
        ];
        let par = PrecomputedPar::build(&models, &stages, &[hydro], None).unwrap();

        // Observations at window_year=1990, season 0 (Jan) and season 1 (Feb).
        let history = vec![
            make_row(hydro, 1990, 0, 120.0),
            make_row(hydro, 1990, 1, 90.0),
        ];

        let mut lib = HistoricalScenarioLibrary::new(1, 2, 1, 0, vec![1990]);
        standardize_historical_windows(
            &mut lib,
            &history,
            &[hydro],
            &stages,
            &par,
            &[1990],
            None,
            &[],
            &[],
        );

        let expected_0 = (120.0 - 100.0) / 30.0;
        let expected_1 = (90.0 - 100.0) / 30.0;

        assert!(
            (lib.eta_slice(0, 0)[0] - expected_0).abs() < 1e-10,
            "AR(0) eta stage 0: expected {expected_0}, got {}",
            lib.eta_slice(0, 0)[0]
        );
        assert!(
            (lib.eta_slice(0, 1)[0] - expected_1).abs() < 1e-10,
            "AR(0) eta stage 1: expected {expected_1}, got {}",
            lib.eta_slice(0, 1)[0]
        );
    }

    // -----------------------------------------------------------------------
    // Test 2: AR(1) uses RAW historical lags (not reconstructed)
    // -----------------------------------------------------------------------

    /// Single hydro, AR(1), `psi_orig`=0.5 in original units, base=80, sigma=25.
    ///
    /// Uses 12 monthly stages so that the lag season (one step before Jan) is Dec.
    /// With `n_seasons`=12 and `max_order`=1, under the new convention:
    ///   - Lag season: (0 - 1 + 12) % 12 = 11 (Dec), `year_offset` = -1
    ///   - Study seasons: 0..11, `year_offset` = 0
    ///   - Full sequence: [(-1,11),(0,0),(0,1),...,(0,11)]
    ///   - Lag observation: (`window_year` - 1, season 11) = (1989, Dec) → month0=11
    ///   - Stage 0 observation: (`window_year` + 0, season 0) = (1990, Jan) → month0=0
    ///   - Stage 1 observation: (1990, Feb) → month0=1
    ///
    /// For stage 0: eta = (130 - 80 - 0.5*110) / 25 = -5/25 = -0.2
    /// For stage 1: lags are RAW → lag[0] = 130.0 (Jan 1990 raw observation)
    ///              eta = (95 - 80 - 0.5*130) / 25 = -50/25 = -2.0
    ///
    /// PAR parametrisation: mean=160, std=25, `psi_star`=0.5 (when stds equal,
    /// `psi_orig` = `psi_star`). base = mean - `psi_orig`*`mean_lag` = 160 - 0.5*160 = 80.
    #[test]
    fn test_ar1_standardization_uses_raw_lags() {
        let hydro = EntityId(1);
        // 12 monthly stages. n_seasons=12.
        let stages = twelve_monthly_stages();

        // Build PAR models for study stages (stage_id 0-11) plus one pre-study
        // stage (stage_id=-1, season 11) needed for coefficient unit conversion.
        let all_stage_ids: Vec<i32> = std::iter::once(-1_i32).chain(0..12_i32).collect();
        let models: Vec<InflowModel> = all_stage_ids
            .iter()
            .map(|&sid| InflowModel {
                hydro_id: hydro,
                stage_id: sid,
                mean_m3s: 160.0,
                std_m3s: 25.0,
                ar_coefficients: vec![0.5],
                residual_std_ratio: 1.0,
                annual: None,
            })
            .collect();
        let par = PrecomputedPar::build(&models, &stages, &[hydro], None).unwrap();

        // Verify precomputed values: psi_orig = 0.5 * 25/25 = 0.5;
        // base = 160 - 0.5*160 = 80; sigma = 25.
        assert!(
            (par.deterministic_base(0, 0) - 80.0).abs() < 1e-10,
            "expected base=80, got {}",
            par.deterministic_base(0, 0)
        );
        assert!(
            (par.sigma(0, 0) - 25.0).abs() < 1e-10,
            "expected sigma=25, got {}",
            par.sigma(0, 0)
        );

        // Window year 1990, max_order=1 (new convention: study at window_year):
        //   lag: (1989, season 11 = Dec) → 110.0
        //   stage 0: (1990, season 0 = Jan) → 130.0
        //   stage 1: (1990, season 1 = Feb) → 95.0
        //   (remaining study stages: use 100.0, not used in assertions)
        let mut history = vec![
            make_row(hydro, 1989, 11, 110.0), // Dec 1989 = pre-study lag
            make_row(hydro, 1990, 0, 130.0),  // Jan 1990 = stage 0
            make_row(hydro, 1990, 1, 95.0),   // Feb 1990 = stage 1
        ];
        for m in 2..12_u32 {
            history.push(make_row(hydro, 1990, m, 100.0));
        }

        // past_inflows: provide lag-1 = Dec 1989 = 110.0 so the rolling chain
        // starts from the same value as the window's own pre-study observation.
        // This keeps all existing eta assertions numerically unchanged while
        // exercising the new past_inflows-seeded chain.
        let past = vec![HydroPastInflows {
            hydro_id: hydro,
            values_m3s: vec![110.0], // lag-1 = Dec 1989
            season_ids: None,
        }];

        let mut lib = HistoricalScenarioLibrary::new(1, 12, 1, 1, vec![1990]);
        standardize_historical_windows(
            &mut lib,
            &history,
            &[hydro],
            &stages,
            &par,
            &[1990],
            None,
            &past,
            &[],
        );

        // Stage 0: lag_state seeded from past_inflows → lag = 110.0.
        // eta = (130 - 80 - 0.5*110) / 25 = -5/25 = -0.2
        let eta_0 = lib.eta_slice(0, 0)[0];
        let expected_0 = (130.0 - 80.0 - 0.5 * 110.0) / 25.0;
        assert!(
            (eta_0 - expected_0).abs() < 1e-10,
            "eta stage 0: expected {expected_0}, got {eta_0}"
        );

        // Stage 1: lag_state advanced by finalize after stage 0 → lag = 130.0 (Jan 1990).
        // eta = (95 - 80 - 0.5*130) / 25 = (95 - 145) / 25 = -2.0
        let eta_1 = lib.eta_slice(0, 1)[0];
        let expected_1 = (95.0 - 80.0 - 0.5 * 130.0) / 25.0;
        assert!(
            (eta_1 - expected_1).abs() < 1e-10,
            "eta stage 1: expected {expected_1}, got {eta_1} (rolling chain: lag=130.0)"
        );
    }

    // -----------------------------------------------------------------------
    // Test 3: 2 hydros, 2 windows — all slices populated independently
    // -----------------------------------------------------------------------

    /// Two hydros, two window years, AR(0). All four (window,stage) slices must
    /// be independently populated with the correct eta values.
    ///
    /// Uses `season_ids` 0 and 1 (`n_seasons`=2). With `max_order`=0 and seasons 0→1
    /// (increasing, no wrap), all `year_offset`s are 0. Observations are at
    /// year = `window_year` itself.
    #[test]
    fn test_multi_hydro_multi_window() {
        let h1 = EntityId(1);
        let h2 = EntityId(2);
        // 2 stages, season_ids 0 and 1 (n_seasons=2, no wrap → year_offset=0).
        let stages = vec![make_monthly_stage(0, 0), make_monthly_stage(1, 1)];

        let models = vec![
            InflowModel {
                hydro_id: h1,
                stage_id: 0,
                mean_m3s: 100.0,
                std_m3s: 10.0,
                ar_coefficients: vec![],
                residual_std_ratio: 1.0,
                annual: None,
            },
            InflowModel {
                hydro_id: h1,
                stage_id: 1,
                mean_m3s: 100.0,
                std_m3s: 10.0,
                ar_coefficients: vec![],
                residual_std_ratio: 1.0,
                annual: None,
            },
            InflowModel {
                hydro_id: h2,
                stage_id: 0,
                mean_m3s: 200.0,
                std_m3s: 20.0,
                ar_coefficients: vec![],
                residual_std_ratio: 1.0,
                annual: None,
            },
            InflowModel {
                hydro_id: h2,
                stage_id: 1,
                mean_m3s: 200.0,
                std_m3s: 20.0,
                ar_coefficients: vec![],
                residual_std_ratio: 1.0,
                annual: None,
            },
        ];
        let par = PrecomputedPar::build(&models, &stages, &[h1, h2], None).unwrap();

        // Observations at year=window_year (year_offset=0 for both stages).
        let history = vec![
            // Window 1990: year_offset=0 → obs at 1990
            make_row(h1, 1990, 0, 110.0),
            make_row(h1, 1990, 1, 90.0),
            make_row(h2, 1990, 0, 220.0),
            make_row(h2, 1990, 1, 180.0),
            // Window 1991: obs at 1991
            make_row(h1, 1991, 0, 105.0),
            make_row(h1, 1991, 1, 95.0),
            make_row(h2, 1991, 0, 210.0),
            make_row(h2, 1991, 1, 190.0),
        ];

        let mut lib = HistoricalScenarioLibrary::new(2, 2, 2, 0, vec![1990, 1991]);
        standardize_historical_windows(
            &mut lib,
            &history,
            &[h1, h2],
            &stages,
            &par,
            &[1990, 1991],
            None,
            &[],
            &[],
        );

        // All 4 (window, stage) slices have length 2 (n_hydros).
        for w in 0..2 {
            for t in 0..2 {
                assert_eq!(
                    lib.eta_slice(w, t).len(),
                    2,
                    "eta slice (w={w}, t={t}) must have length n_hydros=2"
                );
            }
        }

        // Window 0, stage 0: h1=(110-100)/10=1.0, h2=(220-200)/20=1.0
        let e00 = lib.eta_slice(0, 0);
        assert!((e00[0] - 1.0).abs() < 1e-10, "w=0,t=0,h=0: {}", e00[0]);
        assert!((e00[1] - 1.0).abs() < 1e-10, "w=0,t=0,h=1: {}", e00[1]);

        // Window 1, stage 1: h1=(95-100)/10=-0.5, h2=(190-200)/20=-0.5
        let e11 = lib.eta_slice(1, 1);
        assert!((e11[0] - (-0.5)).abs() < 1e-10, "w=1,t=1,h=0: {}", e11[0]);
        assert!((e11[1] - (-0.5)).abs() < 1e-10, "w=1,t=1,h=1: {}", e11[1]);
    }

    // -----------------------------------------------------------------------
    // Test 5: sigma=0 matching deterministic value stores eta=0.0
    // -----------------------------------------------------------------------

    /// Single hydro, single stage, sigma=0, observation matches deterministic value.
    ///
    /// With 1 stage (`season_id`=0, `n_seasons`=1) and `max_order`=0:
    ///   Full sequence: [(0,0)]. `year_offset`=0. Observation at year=`window_year`.
    #[test]
    fn test_sigma_zero_returns_zero_eta() {
        let hydro = EntityId(1);
        // Single stage, season_id=0.
        let stages = vec![make_monthly_stage(0, 0)];
        let models = vec![InflowModel {
            hydro_id: hydro,
            stage_id: 0,
            mean_m3s: 50.0,
            std_m3s: 0.0,
            ar_coefficients: vec![],
            residual_std_ratio: 1.0,
            annual: None,
        }];
        let par = PrecomputedPar::build(&models, &stages, &[hydro], None).unwrap();

        // n_seasons=1, max_order=0: full_sequence = [(0,0)].
        // Observation at (window_year + 0, season 0) = (2000, month 0 = Jan 2000).
        let history = vec![make_row(hydro, 2000, 0, 50.0)];

        let mut lib = HistoricalScenarioLibrary::new(1, 1, 1, 0, vec![2000]);
        standardize_historical_windows(
            &mut lib,
            &history,
            &[hydro],
            &stages,
            &par,
            &[2000],
            None,
            &[],
            &[],
        );

        let eta = lib.eta_slice(0, 0)[0];
        assert!(
            eta == 0.0,
            "sigma=0 with obs matching deterministic value must give eta=0.0, got {eta}"
        );
    }

    // -----------------------------------------------------------------------
    // validate_historical_library tests
    // -----------------------------------------------------------------------

    use super::validate_historical_library;
    use crate::StochasticError;

    /// Build a minimal valid Stage with `season_id`: `Some(season)`.
    #[allow(clippy::cast_possible_truncation, clippy::cast_possible_wrap)]
    fn make_validate_stage(index: usize, season_id: Option<usize>) -> Stage {
        Stage {
            index,
            id: index as i32,
            start_date: NaiveDate::from_ymd_opt(2024, 1, 1).unwrap(),
            end_date: NaiveDate::from_ymd_opt(2024, 2, 1).unwrap(),
            season_id,
            blocks: vec![Block {
                index: 0,
                name: "B".to_string(),
                duration_hours: 720.0,
            }],
            block_mode: BlockMode::Parallel,
            state_config: StageStateConfig {
                storage: true,
                inflow_lags: false,
            },
            risk_config: StageRiskConfig::Expectation,
            scenario_config: ScenarioSourceConfig {
                branching_factor: 1,
                noise_method: NoiseMethod::Saa,
            },
        }
    }

    /// Given a valid library with 5 windows, 12 stages, 3 hydros, all finite
    /// eta values, `validate_historical_library` returns `Ok(())`.
    #[test]
    fn test_valid_library_passes() {
        let n_windows = 5;
        let n_stages = 12;
        let n_hydros = 3;
        let lib = HistoricalScenarioLibrary::new(
            n_windows,
            n_stages,
            n_hydros,
            1,
            (1990..1995).collect(),
        );
        let stages: Vec<Stage> = (0..n_stages)
            .map(|i| make_validate_stage(i, Some(i % 12)))
            .collect();
        let hydro_ids: Vec<EntityId> = (1..=3).map(EntityId).collect();

        let result = validate_historical_library(&lib, &[], &hydro_ids, &stages, 1, None, 5);
        assert!(result.is_ok(), "expected Ok(()), got: {result:?}");
    }

    /// Given a library where `eta_slice(2, 5)[1]` is `f64::NEG_INFINITY`,
    /// `validate_historical_library` returns `Err` with a message containing
    /// "V2.3" and `"NEG_INFINITY"`.
    #[test]
    fn test_neg_infinity_eta_fails_v2_3() {
        let n_windows = 5;
        let n_stages = 12;
        let n_hydros = 3;
        let mut lib = HistoricalScenarioLibrary::new(
            n_windows,
            n_stages,
            n_hydros,
            1,
            (1990..1995).collect(),
        );
        // Inject NEG_INFINITY at window=2, stage=5, hydro=1.
        lib.eta_slice_mut(2, 5)[1] = f64::NEG_INFINITY;

        let stages: Vec<Stage> = (0..n_stages)
            .map(|i| make_validate_stage(i, Some(i % 12)))
            .collect();
        let hydro_ids: Vec<EntityId> = (1..=3).map(EntityId).collect();

        let result = validate_historical_library(&lib, &[], &hydro_ids, &stages, 1, None, 5);
        match result {
            Err(StochasticError::InsufficientData { context }) => {
                assert!(
                    context.contains("V2.3"),
                    "expected message to contain 'V2.3', got: {context}"
                );
                assert!(
                    context.contains("NEG_INFINITY"),
                    "expected message to contain 'NEG_INFINITY', got: {context}"
                );
            }
            other => panic!("expected Err(InsufficientData), got: {other:?}"),
        }
    }

    /// Given a stage with `season_id: None`, `validate_historical_library`
    /// returns `Err` with a message containing "V2.1" and `"season_id"`.
    #[test]
    fn test_missing_season_id_fails_v2_1() {
        let lib = HistoricalScenarioLibrary::new(1, 2, 2, 0, vec![1990]);
        let stages = vec![
            make_validate_stage(0, Some(0)),
            make_validate_stage(1, None), // missing season_id
        ];
        let hydro_ids = vec![EntityId(1), EntityId(2)];

        let result = validate_historical_library(&lib, &[], &hydro_ids, &stages, 0, None, 1);
        match result {
            Err(StochasticError::InsufficientData { context }) => {
                assert!(
                    context.contains("V2.1"),
                    "expected message to contain 'V2.1', got: {context}"
                );
                assert!(
                    context.contains("season_id"),
                    "expected message to contain 'season_id', got: {context}"
                );
            }
            other => panic!("expected Err(InsufficientData), got: {other:?}"),
        }
    }

    /// Given `library.n_hydros() = 3` but `hydro_ids.len() = 4`,
    /// `validate_historical_library` returns `Err` with a message containing
    /// "V2.9".
    #[test]
    fn test_hydro_count_mismatch_fails_v2_9() {
        let lib = HistoricalScenarioLibrary::new(1, 1, 3, 0, vec![1990]);
        let stages = vec![make_validate_stage(0, Some(0))];
        // 4 hydro IDs but library has n_hydros=3.
        let hydro_ids = vec![EntityId(1), EntityId(2), EntityId(3), EntityId(4)];

        let result = validate_historical_library(&lib, &[], &hydro_ids, &stages, 0, None, 1);
        match result {
            Err(StochasticError::InsufficientData { context }) => {
                assert!(
                    context.contains("V2.9"),
                    "expected message to contain 'V2.9', got: {context}"
                );
            }
            other => panic!("expected Err(InsufficientData), got: {other:?}"),
        }
    }

    /// Given `library.n_windows() = 5` and `forward_passes = 20`,
    /// `validate_historical_library` returns `Ok(())` (warning is emitted
    /// via tracing but does not abort construction).
    #[test]
    fn test_pool_warning_path_returns_ok() {
        let lib = HistoricalScenarioLibrary::new(5, 1, 2, 0, (1990..1995).collect());
        let stages = vec![make_validate_stage(0, Some(0))];
        let hydro_ids = vec![EntityId(1), EntityId(2)];

        // 5 windows < 20 forward passes triggers warn! but must still return Ok(()).
        let result = validate_historical_library(&lib, &[], &hydro_ids, &stages, 0, None, 20);
        assert!(
            result.is_ok(),
            "warning path must return Ok(()), got: {result:?}"
        );
    }

    // -----------------------------------------------------------------------
    // Helpers for SeasonMap-aware standardize tests
    // -----------------------------------------------------------------------

    /// Build a standard monthly `SeasonMap` (12 seasons, IDs 0–11).
    fn monthly_season_map() -> SeasonMap {
        let seasons = (0_usize..12)
            .map(|i| SeasonDefinition {
                id: i,
                label: format!("Month{i}"),
                #[allow(clippy::cast_possible_truncation)]
                month_start: (i as u32) + 1,
                day_start: None,
                month_end: None,
                day_end: None,
            })
            .collect();
        SeasonMap {
            cycle_type: SeasonCycleType::Monthly,
            seasons,
        }
    }

    // -----------------------------------------------------------------------
    // Test 6: monthly SeasonMap produces bit-for-bit identical eta values
    // -----------------------------------------------------------------------

    /// Given a monthly `SeasonMap` with IDs 0–11, calling
    /// `standardize_historical_windows` with `Some(&sm)` must produce
    /// bit-for-bit identical eta values to calling with `None` (month0 fallback).
    ///
    /// Uses 1 hydro, AR(0), 2 monthly stages (seasons 0 and 1).
    #[test]
    fn test_standardize_monthly_season_map_identical() {
        let hydro = EntityId(1);
        let stages = vec![make_monthly_stage(0, 0), make_monthly_stage(1, 1)];
        let models = vec![
            InflowModel {
                hydro_id: hydro,
                stage_id: 0,
                mean_m3s: 80.0,
                std_m3s: 20.0,
                ar_coefficients: vec![],
                residual_std_ratio: 1.0,
                annual: None,
            },
            InflowModel {
                hydro_id: hydro,
                stage_id: 1,
                mean_m3s: 60.0,
                std_m3s: 15.0,
                ar_coefficients: vec![],
                residual_std_ratio: 1.0,
                annual: None,
            },
        ];
        let par = PrecomputedPar::build(&models, &stages, &[hydro], None).unwrap();

        // Observations at window_year=2000, seasons 0 (Jan) and 1 (Feb).
        let history = vec![
            make_row(hydro, 2000, 0, 100.0),
            make_row(hydro, 2000, 1, 45.0),
        ];

        let sm = monthly_season_map();

        // Run with None (month0 fallback).
        let mut lib_none = HistoricalScenarioLibrary::new(1, 2, 1, 0, vec![2000]);
        standardize_historical_windows(
            &mut lib_none,
            &history,
            &[hydro],
            &stages,
            &par,
            &[2000],
            None,
            &[],
            &[],
        );

        // Run with monthly SeasonMap.
        let mut lib_sm = HistoricalScenarioLibrary::new(1, 2, 1, 0, vec![2000]);
        standardize_historical_windows(
            &mut lib_sm,
            &history,
            &[hydro],
            &stages,
            &par,
            &[2000],
            Some(&sm),
            &[],
            &[],
        );

        // Bit-for-bit identical for all (window, stage) slices.
        for t in 0..2 {
            assert_eq!(
                lib_none.eta_slice(0, t),
                lib_sm.eta_slice(0, t),
                "eta values at stage {t} must be identical between None and monthly SeasonMap"
            );
        }
    }

    // -----------------------------------------------------------------------
    // Test 7: quarterly SeasonMap maps observations to correct season IDs
    // -----------------------------------------------------------------------

    #[allow(clippy::cast_possible_truncation, clippy::cast_possible_wrap)]
    fn make_quarterly_stage(index: usize, season_id: usize) -> Stage {
        let month = (season_id as u32) * 3 + 1; // 1, 4, 7, 10
        Stage {
            index,
            id: index as i32,
            start_date: NaiveDate::from_ymd_opt(2024, month, 1).unwrap(),
            end_date: NaiveDate::from_ymd_opt(
                2024,
                if month + 2 <= 12 { month + 2 } else { 12 },
                28,
            )
            .unwrap(),
            season_id: Some(season_id),
            blocks: vec![Block {
                index: 0,
                name: "SINGLE".to_string(),
                duration_hours: 2160.0,
            }],
            block_mode: BlockMode::Parallel,
            state_config: StageStateConfig {
                storage: true,
                inflow_lags: false,
            },
            risk_config: StageRiskConfig::Expectation,
            scenario_config: ScenarioSourceConfig {
                branching_factor: 1,
                noise_method: NoiseMethod::Saa,
            },
        }
    }

    fn quarterly_season_map() -> SeasonMap {
        SeasonMap {
            cycle_type: SeasonCycleType::Custom,
            seasons: vec![
                SeasonDefinition {
                    id: 0,
                    label: "Q1".to_string(),
                    month_start: 1,
                    day_start: Some(1),
                    month_end: Some(3),
                    day_end: Some(31),
                },
                SeasonDefinition {
                    id: 1,
                    label: "Q2".to_string(),
                    month_start: 4,
                    day_start: Some(1),
                    month_end: Some(6),
                    day_end: Some(30),
                },
                SeasonDefinition {
                    id: 2,
                    label: "Q3".to_string(),
                    month_start: 7,
                    day_start: Some(1),
                    month_end: Some(9),
                    day_end: Some(30),
                },
                SeasonDefinition {
                    id: 3,
                    label: "Q4".to_string(),
                    month_start: 10,
                    day_start: Some(1),
                    month_end: Some(12),
                    day_end: Some(31),
                },
            ],
        }
    }

    #[test]
    fn test_standardize_quarterly_season_map_correct() {
        let hydro = EntityId(1);
        let stages: Vec<Stage> = (0..4).map(|i| make_quarterly_stage(i, i)).collect();
        let models: Vec<InflowModel> = (0_i32..4)
            .map(|i| InflowModel {
                hydro_id: hydro,
                stage_id: i,
                mean_m3s: 100.0 + f64::from(i) * 10.0,
                std_m3s: 10.0,
                ar_coefficients: vec![],
                residual_std_ratio: 1.0,
                annual: None,
            })
            .collect();
        let par = PrecomputedPar::build(&models, &stages, &[hydro], None).unwrap();

        let sm = quarterly_season_map();

        // Quarterly observations: one per quarter for year 2000.
        // Jan 1 -> Q1 (season 0), Apr 1 -> Q2 (season 1),
        // Jul 1 -> Q3 (season 2), Oct 1 -> Q4 (season 3).
        let history = vec![
            InflowHistoryRow {
                hydro_id: hydro,
                date: NaiveDate::from_ymd_opt(2000, 1, 1).unwrap(),
                value_m3s: 120.0,
            },
            InflowHistoryRow {
                hydro_id: hydro,
                date: NaiveDate::from_ymd_opt(2000, 4, 1).unwrap(),
                value_m3s: 130.0,
            },
            InflowHistoryRow {
                hydro_id: hydro,
                date: NaiveDate::from_ymd_opt(2000, 7, 1).unwrap(),
                value_m3s: 140.0,
            },
            InflowHistoryRow {
                hydro_id: hydro,
                date: NaiveDate::from_ymd_opt(2000, 10, 1).unwrap(),
                value_m3s: 150.0,
            },
        ];

        let mut lib = HistoricalScenarioLibrary::new(1, 4, 1, 0, vec![2000]);
        standardize_historical_windows(
            &mut lib,
            &history,
            &[hydro],
            &stages,
            &par,
            &[2000],
            Some(&sm),
            &[],
            &[],
        );

        // eta = (obs - mean) / std for each season (AR(0)):
        // Q1: (120 - 100) / 10 = 2.0
        // Q2: (130 - 110) / 10 = 2.0
        // Q3: (140 - 120) / 10 = 2.0
        // Q4: (150 - 130) / 10 = 2.0
        for t in 0..4 {
            let eta = lib.eta_slice(0, t)[0];
            assert!(
                (eta - 2.0).abs() < 1e-10,
                "stage {t}: expected eta=2.0, got {eta}"
            );
        }
    }

    // -----------------------------------------------------------------------
    // Helper: quarterly_history
    // -----------------------------------------------------------------------

    /// Build a quarterly history for `hydro_id` spanning years [`from_year`, `to_year`].
    ///
    /// One row per quarter per year, dated on the 1st of Jan, Apr, Jul, Oct.
    /// All `value_m3s` values default to 100.0; callers overwrite specific years
    /// when they need distinctive per-quarter observations.
    fn quarterly_history(
        hydro_id: EntityId,
        from_year: i32,
        to_year: i32,
    ) -> Vec<InflowHistoryRow> {
        let quarter_months = [1u32, 4, 7, 10];
        (from_year..=to_year)
            .flat_map(|y| {
                quarter_months.iter().map(move |&m| InflowHistoryRow {
                    hydro_id,
                    date: NaiveDate::from_ymd_opt(y, m, 1).unwrap(),
                    value_m3s: 100.0,
                })
            })
            .collect()
    }

    // -----------------------------------------------------------------------
    // Test 8: None season_map backward compatibility
    // -----------------------------------------------------------------------

    /// Given `season_map = None`, calling `standardize_historical_windows`
    /// falls back to `month0()` and produces the expected eta values.
    ///
    /// Uses 1 hydro, AR(0), 1 monthly stage (season 0, Jan).
    /// Expected: eta = (obs - mean) / std = (110 - 90) / 10 = 2.0.
    #[test]
    fn test_standardize_none_season_map_backward_compat() {
        let hydro = EntityId(1);
        let stages = vec![make_monthly_stage(0, 0)];
        let models = vec![InflowModel {
            hydro_id: hydro,
            stage_id: 0,
            mean_m3s: 90.0,
            std_m3s: 10.0,
            ar_coefficients: vec![],
            residual_std_ratio: 1.0,
            annual: None,
        }];
        let par = PrecomputedPar::build(&models, &stages, &[hydro], None).unwrap();

        // Observation at (window_year=1995, season 0 = Jan) via month0() = 0.
        let history = vec![make_row(hydro, 1995, 0, 110.0)];

        let mut lib = HistoricalScenarioLibrary::new(1, 1, 1, 0, vec![1995]);
        standardize_historical_windows(
            &mut lib,
            &history,
            &[hydro],
            &stages,
            &par,
            &[1995],
            None,
            &[],
            &[],
        );

        let eta = lib.eta_slice(0, 0)[0];
        let expected = (110.0 - 90.0) / 10.0; // = 2.0
        assert!(
            (eta - expected).abs() < 1e-10,
            "None season_map backward compat: expected eta={expected}, got {eta}"
        );
    }

    // -----------------------------------------------------------------------
    // Test 9: quarterly SeasonMap produces correctly season-aligned eta values
    //         for 2 hydros using `standardize_historical_windows`
    // -----------------------------------------------------------------------

    /// Given a quarterly (4-season Custom) `SeasonMap` and quarterly history for
    /// 2 hydros from 1990–2010, calling `standardize_historical_windows` with
    /// `season_map = Some(&quarterly_sm)` must produce eta values correctly
    /// aligned to quarterly season boundaries.
    ///
    /// This verifies that the primary `season_map.season_for_date()` path (tier
    /// 2) is taken for lag observations outside the study stage date ranges, and
    /// that the `month0()` fallback (tier 3) is NOT reached.
    ///
    /// Concretely: quarterly observations are dated on Jan 1 (Q1), Apr 1 (Q2),
    /// Jul 1 (Q3), and Oct 1 (Q4). With the quarterly `SeasonMap`, these map to
    /// season IDs 0, 1, 2, 3 respectively. If `month0()` were erroneously used,
    /// Apr 1 → month0=3, Jul 1 → month0=6 (out of range for `n_seasons`=4) and
    /// Oct 1 → month0=9 (out of range), so Q2/Q3/Q4 observations would be
    /// silently dropped and eta would equal 0.0 (missing-observation fallback)
    /// rather than the correct season-aligned value.
    ///
    /// Setup:
    ///   - 2 hydros (h1 and h2), quarterly history 1990–2010 (`max_order`=0).
    ///   - 4 quarterly stages (`season_ids` 0–3, AR(0)).
    ///   - Distinct mean per season so each quarter's eta is independently
    ///     verifiable.
    ///   - Window year 2000: h1 obs (90, 110, 115, 120), h2 obs (85, 95, 105, 125).
    ///
    /// Expected eta (AR(0), eta = (obs - mean) / std):
    ///   h1: Q1=(90-80)/10=1.0, Q2=(110-90)/10=2.0, Q3=(115-100)/10=1.5, Q4=(120-110)/10=1.0
    ///   h2: Q1=(85-80)/10=0.5, Q2=(95-90)/10=0.5, Q3=(105-100)/10=0.5, Q4=(125-110)/10=1.5
    #[allow(clippy::cast_possible_truncation, clippy::cast_possible_wrap)]
    #[test]
    fn test_quarterly_standardize_historical_windows() {
        let h1 = EntityId(1);
        let h2 = EntityId(2);

        // 4 quarterly stages: Q1–Q4, season_ids 0–3, AR(0).
        let stages: Vec<Stage> = (0..4).map(|i| make_quarterly_stage(i, i)).collect();

        // Distinct mean per season (identical std=10) so quarter-alignment can
        // be confirmed independently for each season.
        let mean_per_season = [80.0_f64, 90.0, 100.0, 110.0];
        let std_per_season = [10.0_f64; 4];

        // Build AR(0) PAR models for both hydros (same parameters).
        let mut models: Vec<InflowModel> = Vec::new();
        for &hydro in &[h1, h2] {
            for (i, &mean) in mean_per_season.iter().enumerate() {
                models.push(InflowModel {
                    hydro_id: hydro,
                    stage_id: i as i32,
                    mean_m3s: mean,
                    std_m3s: std_per_season[i],
                    ar_coefficients: vec![],
                    residual_std_ratio: 1.0,
                    annual: None,
                });
            }
        }
        let par = PrecomputedPar::build(&models, &stages, &[h1, h2], None).unwrap();

        let sm = quarterly_season_map();

        // Build quarterly history 1990–2010 for both hydros.
        // Default value 100.0 for all observations except window year 2000,
        // where we overwrite with distinctive per-quarter values.
        let quarter_months = [1u32, 4, 7, 10];
        let h1_obs_2000 = [90.0_f64, 110.0, 115.0, 120.0]; // Q1–Q4
        let h2_obs_2000 = [85.0_f64, 95.0, 105.0, 125.0]; // Q1–Q4

        let mut history: Vec<InflowHistoryRow> = quarterly_history(h1, 1990, 2010);
        history.extend(quarterly_history(h2, 1990, 2010));

        // Overwrite the 2000 observations for h1 and h2 with the test values.
        for row in &mut history {
            if row.date.year() == 2000 {
                let q = quarter_months
                    .iter()
                    .position(|&m| m == row.date.month())
                    .expect("date month must be a quarter month");
                if row.hydro_id == h1 {
                    row.value_m3s = h1_obs_2000[q];
                } else if row.hydro_id == h2 {
                    row.value_m3s = h2_obs_2000[q];
                }
            }
        }

        // Use a single window year (2000) with max_order=0 for clarity.
        let window_years = vec![2000_i32];
        let n_hydros = 2;
        let n_stages = 4;
        let mut lib =
            HistoricalScenarioLibrary::new(1, n_stages, n_hydros, 0, window_years.clone());

        standardize_historical_windows(
            &mut lib,
            &history,
            &[h1, h2],
            &stages,
            &par,
            &window_years,
            Some(&sm),
            &[],
            &[],
        );

        // Verify eta values for the single window (w=0).
        // h1 expected: Q1=1.0, Q2=2.0, Q3=1.5, Q4=1.0
        let h1_expected = [1.0_f64, 2.0, 1.5, 1.0];
        // h2 expected: Q1=0.5, Q2=0.5, Q3=0.5, Q4=1.5
        let h2_expected = [0.5_f64, 0.5, 0.5, 1.5];

        for t in 0..n_stages {
            let eta = lib.eta_slice(0, t);
            assert!(
                (eta[0] - h1_expected[t]).abs() < 1e-10,
                "h1 stage {t} (Q{}): expected eta={}, got {} — \
                 seasonal alignment via quarterly SeasonMap must be correct",
                t + 1,
                h1_expected[t],
                eta[0]
            );
            assert!(
                (eta[1] - h2_expected[t]).abs() < 1e-10,
                "h2 stage {t} (Q{}): expected eta={}, got {} — \
                 seasonal alignment via quarterly SeasonMap must be correct",
                t + 1,
                h2_expected[t],
                eta[1]
            );
        }
    }
}
