hydro-core 0.1.0

//! Evaluation metrics for hydrological model comparison.
//!
//! Includes Nash-Sutcliffe Efficiency (= DC 确定性系数), Kling-Gupta Efficiency,
//! RMSE, Percent Bias, R², and GB/T 22482-2008 accuracy grading.

use serde::{Serialize, Deserialize};

/// Nash-Sutcliffe Efficiency (= 确定性系数 DC in Chinese standard).
/// 1.0 = perfect; 0.0 = as good as the observed mean; <0 = worse than mean.
pub fn nse(obs: &[f64], sim: &[f64]) -> f64 {
    if obs.is_empty() || obs.len() != sim.len() {
        return f64::NAN;
    }
    let mean_o = obs.iter().sum::<f64>() / obs.len() as f64;
    let ss_res: f64 = obs.iter().zip(sim).map(|(o, s)| (o - s).powi(2)).sum();
    let ss_tot: f64 = obs.iter().map(|o| (o - mean_o).powi(2)).sum();
    if ss_tot < 1e-12 {
        return f64::NAN;
    }
    1.0 - ss_res / ss_tot
}

/// Kling-Gupta Efficiency (Gupta et al. 2009).
/// Decomposes into correlation (r), variability (α), and bias (β).
pub fn kge(obs: &[f64], sim: &[f64]) -> f64 {
    if obs.len() < 2 || obs.len() != sim.len() {
        return f64::NAN;
    }
    let n = obs.len() as f64;
    let mean_o = obs.iter().sum::<f64>() / n;
    let mean_s = sim.iter().sum::<f64>() / n;
    let std_o = (obs.iter().map(|o| (o - mean_o).powi(2)).sum::<f64>() / n).sqrt();
    let std_s = (sim.iter().map(|s| (s - mean_s).powi(2)).sum::<f64>() / n).sqrt();
    if std_o < 1e-12 || std_s < 1e-12 {
        return f64::NAN;
    }
    // Pearson correlation
    let cov: f64 = obs.iter().zip(sim)
        .map(|(o, s)| (o - mean_o) * (s - mean_s))
        .sum::<f64>() / n;
    let r = cov / (std_o * std_s);
    let alpha = std_s / std_o;
    let beta = mean_s / mean_o;
    let ed = (r - 1.0).powi(2) + (alpha - 1.0).powi(2) + (beta - 1.0).powi(2);
    1.0 - ed.sqrt()
}

/// Root Mean Square Error.
pub fn rmse(obs: &[f64], sim: &[f64]) -> f64 {
    if obs.is_empty() || obs.len() != sim.len() {
        return f64::NAN;
    }
    let ssq: f64 = obs.iter().zip(sim).map(|(o, s)| (o - s).powi(2)).sum();
    (ssq / obs.len() as f64).sqrt()
}

/// Percent Bias (positive = overprediction, negative = underprediction).
pub fn pbias(obs: &[f64], sim: &[f64]) -> f64 {
    let sum_o: f64 = obs.iter().sum();
    if sum_o.abs() < 1e-12 {
        return f64::NAN;
    }
    let sum_s: f64 = sim.iter().sum();
    100.0 * (sum_s - sum_o) / sum_o
}

/// Coefficient of Determination (R²).
pub fn r2(obs: &[f64], sim: &[f64]) -> f64 {
    if obs.len() < 2 || obs.len() != sim.len() {
        return f64::NAN;
    }
    let n = obs.len() as f64;
    let mean_o = obs.iter().sum::<f64>() / n;
    let mean_s = sim.iter().sum::<f64>() / n;
    let ss_oo: f64 = obs.iter().map(|o| (o - mean_o).powi(2)).sum();
    let ss_ss: f64 = sim.iter().map(|s| (s - mean_s).powi(2)).sum();
    let ss_os: f64 = obs.iter().zip(sim).map(|(o, s)| (o - mean_o) * (s - mean_s)).sum();
    if ss_oo < 1e-12 || ss_ss < 1e-12 {
        return f64::NAN;
    }
    (ss_os / (ss_oo * ss_ss).sqrt()).powi(2)
}

/// GB/T 22482-2008 accuracy grade from DC (= NSE).
/// 甲 (excellent) > 0.90; 乙 (good) 0.70–0.90; 丙 (acceptable) 0.50–0.70.
pub fn dc_grade(dc: f64) -> &'static str {
    if dc >= 0.90 { "甲" }
    else if dc >= 0.70 { "乙" }
    else if dc >= 0.50 { "丙" }
    else { "不合格" }
}

// ── GB/T 22482-2008 事件误差 + 合格率 ────────────────────────────────

/// 单场洪水的 GB/T 22482 事件误差(相对误差 %,峰现误差以时段步数计)。
#[derive(Clone, Debug, Serialize, Deserialize, Default)]
pub struct GbtEventErrors {
    /// 洪峰相对误差 % = 100·(Qpk_sim − Qpk_obs)/Qpk_obs
    pub peak_rel_err_pct: f64,
    /// 峰现时间误差(时段步数)= argmax(sim) − argmax(obs);正=滞后,负=超前
    pub time_to_peak_err_steps: i64,
    /// 径流深相对误差 % = 100·(D_sim − D_obs)/D_obs;D 由水量/面积折算(mm)
    pub runoff_depth_rel_err_pct: f64,
}

/// GB/T 22482 许可误差(默认:洪峰 ±20%、峰现 ±1 时段、径流深 ±20%)。
#[derive(Clone, Debug, Serialize, Deserialize)]
pub struct QualificationTolerance {
    pub peak_rel_pct: f64,
    pub time_to_peak_steps: i64,
    pub runoff_depth_rel_pct: f64,
}

impl Default for QualificationTolerance {
    fn default() -> Self {
        Self { peak_rel_pct: 20.0, time_to_peak_steps: 1, runoff_depth_rel_pct: 20.0 }
    }
}

/// 计算单场洪水的 GB/T 22482 事件误差。`area_km2` 用于水量→径流深折算。
pub fn event_errors(obs: &[f64], sim: &[f64], area_km2: f64, dt_h: f64) -> GbtEventErrors {
    let mut e = GbtEventErrors::default();
    if obs.is_empty() || sim.is_empty() || area_km2 <= 0.0 || dt_h <= 0.0 {
        return e;
    }
    let (obs_pk, obs_idx) = argmax(obs);
    let (sim_pk, sim_idx) = argmax(sim);
    if obs_pk.abs() > 1e-9 {
        e.peak_rel_err_pct = 100.0 * (sim_pk - obs_pk) / obs_pk;
    }
    e.time_to_peak_err_steps = sim_idx as i64 - obs_idx as i64;
    let d_obs = runoff_depth_mm(obs, area_km2, dt_h);
    let d_sim = runoff_depth_mm(sim, area_km2, dt_h);
    if d_obs.abs() > 1e-9 {
        e.runoff_depth_rel_err_pct = 100.0 * (d_sim - d_obs) / d_obs;
    }
    e
}

/// 单场洪水是否合格(三项误差均在许可范围内)。
pub fn event_qualified(e: &GbtEventErrors, tol: &QualificationTolerance) -> bool {
    e.peak_rel_err_pct.abs() <= tol.peak_rel_pct
        && e.time_to_peak_err_steps.abs() <= tol.time_to_peak_steps
        && e.runoff_depth_rel_err_pct.abs() <= tol.runoff_depth_rel_pct
}

/// 合格率汇总报告(多场洪水)。
#[derive(Clone, Debug, Serialize, Deserialize, Default)]
pub struct QualificationReport {
    pub total: usize,
    pub qualified: usize,
    /// 合格率 % = 100·qualified/total
    pub rate_pct: f64,
    /// 甲/乙/丙/不合格(甲≥85、乙 70–85、丙 60–70)
    pub grade: String,
}

/// 由各场洪水的合格标志聚合合格率。
pub fn qualified_rate(qualified_flags: &[bool]) -> QualificationReport {
    let total = qualified_flags.len();
    let qualified = qualified_flags.iter().filter(|&&q| q).count();
    let rate_pct = if total == 0 {
        0.0
    } else {
        100.0 * qualified as f64 / total as f64
    };
    QualificationReport {
        total,
        qualified,
        rate_pct,
        grade: qualification_grade(rate_pct).to_string(),
    }
}

/// GB/T 22482 合格率评级:甲≥85%、乙 70–85%、丙 60–70%、否则不合格。
pub fn qualification_grade(rate_pct: f64) -> &'static str {
    if rate_pct >= 85.0 { "甲" }
    else if rate_pct >= 70.0 { "乙" }
    else if rate_pct >= 60.0 { "丙" }
    else { "不合格" }
}

/// 方案等级 = DC 等级 与 合格率等级 取低(丙级为洪水预报业务最低要求)。
pub fn scheme_grade(dc_grade: &str, qual_grade: &str) -> &'static str {
    let rank = |g: &str| -> i32 {
        match g { "甲" => 3, "乙" => 2, "丙" => 1, _ => 0 }
    };
    match rank(dc_grade).min(rank(qual_grade)) {
        3 => "甲",
        2 => "乙",
        1 => "丙",
        _ => "不合格",
    }
}

/// Full metrics report for one model's simulation vs observation.
#[derive(Clone, Debug, Serialize, Deserialize, Default)]
pub struct MetricsReport {
    pub nse: f64,
    pub kge: f64,
    pub rmse: f64,
    pub pbias: f64,
    pub r2: f64,
    pub grade: String,
    /// GB/T 22482 事件误差(需提供 area/dt;否则 None)。
    pub gbt: Option<GbtEventErrors>,
    /// 单场洪水是否合格(GB/T 22482 许可误差内)。
    pub qualified: Option<bool>,
}

/// Compute all metrics at once(不含 GB/T 事件误差;gbt/qualified = None)。
pub fn compute_metrics(obs: &[f64], sim: &[f64]) -> MetricsReport {
    let n = nse(obs, sim);
    MetricsReport {
        nse: n,
        kge: kge(obs, sim),
        rmse: rmse(obs, sim),
        pbias: pbias(obs, sim),
        r2: r2(obs, sim),
        grade: dc_grade(n).to_string(),
        gbt: None,
        qualified: None,
    }
}

/// 全量指标 + GB/T 22482 事件误差 + 合格判定。
pub fn compute_metrics_gbt(
    obs: &[f64],
    sim: &[f64],
    area_km2: f64,
    dt_h: f64,
    tol: &QualificationTolerance,
) -> MetricsReport {
    let mut m = compute_metrics(obs, sim);
    let e = event_errors(obs, sim, area_km2, dt_h);
    m.qualified = Some(event_qualified(&e, tol));
    m.gbt = Some(e);
    m
}

// ── 概率/集合指标(CRPS / POD-FAR-CSI / Brier)────────────────────────

/// 集合概率预报验证报告。
#[derive(Clone, Debug, Serialize, Deserialize, Default)]
pub struct ProbabilisticReport {
    /// 连续排名概率分数(fair CRPS,越低越好;完美集合=0)。
    pub crps: f64,
    /// Brier 分数(阈值超限概率的均方误差,0=完美,1=最差)。
    pub brier: f64,
    /// 命中率 POD = hits/(hits+misses)。
    pub pod: f64,
    /// 误警率 FAR = false_alarms/(hits+false_alarms)。
    pub far: f64,
    /// 临界成功指数 CSI = hits/(hits+misses+false_alarms)。
    pub csi: f64,
}

/// 逐时步集合超限概率(成员 > 阈值 的比例)。成员长度不齐时按各自是否可达该时步计。
pub fn ensemble_exceedance(ensemble: &[Vec<f64>], threshold: f64) -> Vec<f64> {
    if ensemble.is_empty() {
        return Vec::new();
    }
    let n_t = ensemble.iter().map(|m| m.len()).max().unwrap_or(0);
    let m = ensemble.len() as f64;
    (0..n_t)
        .map(|t| {
            let cnt = ensemble.iter()
                .filter(|mem| mem.get(t).map_or(false, |v| *v > threshold))
                .count() as f64;
            cnt / m
        })
        .collect()
}

/// fair CRPS(Gneiting & Raftery 2007;Ferro 2017):逐时步平均。
/// 每步 `CRPS_t = (1/M)Σ|x_m−y| − (1/(2M²))ΣΣ|x_m−x_n|`;完美集合(成员全=y)→ 0。
/// 单成员集合退化为 MAE。成员该时步缺失则跳过该步。
pub fn crps(ensemble: &[Vec<f64>], observed: &[f64]) -> f64 {
    if ensemble.is_empty() || observed.is_empty() {
        return f64::NAN;
    }
    let m = ensemble.len();
    let mf = m as f64;
    let mut sum = 0.0;
    let mut count = 0usize;
    for (t, &y) in observed.iter().enumerate() {
        let vals: Vec<f64> = ensemble.iter().filter_map(|mem| mem.get(t).copied()).collect();
        if vals.len() != m {
            continue; // 该时步成员不全,跳过
        }
        let mae: f64 = vals.iter().map(|x| (x - y).abs()).sum::<f64>() / mf;
        let mut pair = 0.0;
        for i in 0..m {
            for j in 0..m {
                pair += (vals[i] - vals[j]).abs();
            }
        }
        let spread = pair / (2.0 * mf * mf);
        sum += (mae - spread).max(0.0);
        count += 1;
    }
    if count == 0 { f64::NAN } else { sum / count as f64 }
}

/// Brier 分数:超限概率 vs 超限指示的均方误差。0=完美,1=最差。
pub fn brier(prob_exceed: &[f64], observed_exceed: &[bool]) -> f64 {
    if prob_exceed.is_empty() || prob_exceed.len() != observed_exceed.len() {
        return f64::NAN;
    }
    let n = prob_exceed.len() as f64;
    prob_exceed.iter().zip(observed_exceed.iter())
        .map(|(p, o)| (p - if *o { 1.0 } else { 0.0 }).powi(2))
        .sum::<f64>() / n
}

/// 阈值超限列联表 → (POD, FAR, CSI)。
/// POD = a/(a+b)、FAR = c/(a+c)、CSI = a/(a+b+c);对应分母为 0 时该量返回 NAN。
pub fn pod_far_csi(forecast_exceed: &[bool], observed_exceed: &[bool]) -> (f64, f64, f64) {
    let mut hits = 0.0;         // a: 预报超 ∧ 观测超
    let mut misses = 0.0;       // b: 预报未超 ∧ 观测超
    let mut false_alarms = 0.0; // c: 预报超 ∧ 观测未超
    for (f, o) in forecast_exceed.iter().zip(observed_exceed.iter()) {
        match (*f, *o) {
            (true, true) => hits += 1.0,
            (false, true) => misses += 1.0,
            (true, false) => false_alarms += 1.0,
            (false, false) => {} // 正确拒绝
        }
    }
    let pod = if hits + misses > 0.0 { hits / (hits + misses) } else { f64::NAN };
    let far = if hits + false_alarms > 0.0 { false_alarms / (hits + false_alarms) } else { f64::NAN };
    let csi = if hits + misses + false_alarms > 0.0 {
        hits / (hits + misses + false_alarms)
    } else { f64::NAN };
    (pod, far, csi)
}

/// 集合概率预报全量验证。forecast 超限 = 集合超限概率 > 0.5;observed 超限 = obs > 阈值。
/// 各时步按 observed 与集合的最短长度对齐。
pub fn compute_probabilistic(
    ensemble: &[Vec<f64>],
    observed: &[f64],
    threshold: f64,
) -> ProbabilisticReport {
    let prob = ensemble_exceedance(ensemble, threshold);
    let n = observed.len().min(prob.len());
    let prob_aligned: Vec<f64> = prob.iter().take(n).copied().collect();
    let obs_exceed: Vec<bool> = observed[..n].iter().map(|&q| q > threshold).collect();
    let fcst_exceed: Vec<bool> = prob_aligned.iter().map(|&p| p > 0.5).collect();
    let (pod, far, csi) = pod_far_csi(&fcst_exceed, &obs_exceed);
    ProbabilisticReport {
        crps: crps(ensemble, observed),
        brier: brier(&prob_aligned, &obs_exceed),
        pod, far, csi,
    }
}

/// 集合秩直方图(Talagrand):逐时步把观测在 M 个成员中的秩(0..=M)累计成 M+1 个桶。
/// 桶越均匀 → 集合越可靠(可信);U 形 → 欠散布(过自信);钟形 → 过散布。
pub fn rank_histogram(ensemble: &[Vec<f64>], observed: &[f64]) -> Vec<usize> {
    if ensemble.is_empty() || observed.is_empty() {
        return Vec::new();
    }
    let m = ensemble.len();
    let mut hist = vec![0usize; m + 1];
    for (t, &y) in observed.iter().enumerate() {
        let vals: Vec<f64> = ensemble.iter().filter_map(|mem| mem.get(t).copied()).collect();
        if vals.len() != m {
            continue;
        }
        // rank = 严格小于观测的成员数(观测落在该秩区间)
        let rank = vals.iter().filter(|&&x| x < y).count();
        hist[rank] += 1;
    }
    hist
}

// ── helpers ──

fn argmax(xs: &[f64]) -> (f64, usize) {
    let mut best = f64::NEG_INFINITY;
    let mut idx = 0;
    for (i, &v) in xs.iter().enumerate() {
        if v > best {
            best = v;
            idx = i;
        }
    }
    (best, idx)
}

/// 径流深(mm)= ΣQ·dt_h·3.6 / area_km2(Q m³/s, dt_h 小时, area km²)。
fn runoff_depth_mm(q: &[f64], area_km2: f64, dt_h: f64) -> f64 {
    let sum_q: f64 = q.iter().sum();
    sum_q * dt_h * 3.6 / area_km2
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_nse_perfect() {
        let obs = vec![1.0, 2.0, 3.0, 4.0, 5.0];
        assert!((nse(&obs, &obs) - 1.0).abs() < 1e-10);
    }

    #[test]
    fn test_nse_mean_prediction() {
        let obs = vec![1.0, 2.0, 3.0, 4.0, 5.0];
        let sim = vec![3.0; 5]; // mean
        assert!((nse(&obs, &sim) - 0.0).abs() < 1e-10);
    }

    #[test]
    fn test_kge_perfect() {
        let obs = vec![10.0, 20.0, 30.0, 40.0];
        assert!((kge(&obs, &obs) - 1.0).abs() < 1e-10);
    }

    #[test]
    fn test_rmse() {
        let obs = vec![1.0, 2.0, 3.0];
        let sim = vec![1.0, 2.0, 4.0];
        // RMSE = sqrt((0+0+1)/3) = sqrt(1/3) ≈ 0.577
        assert!((rmse(&obs, &sim) - (1.0_f64 / 3.0).sqrt()).abs() < 1e-10);
    }

    #[test]
    fn test_pbias() {
        let obs = vec![100.0, 200.0];
        let sim = vec![110.0, 220.0];
        // PBIAS = 100*(330-300)/300 = 10%
        assert!((pbias(&obs, &sim) - 10.0).abs() < 1e-10);
    }

    #[test]
    fn test_r2_perfect() {
        let obs = vec![1.0, 2.0, 3.0, 4.0, 5.0];
        assert!((r2(&obs, &obs) - 1.0).abs() < 1e-10);
    }

    #[test]
    fn test_dc_grade() {
        assert_eq!(dc_grade(0.95), "甲");
        assert_eq!(dc_grade(0.85), "乙");
        assert_eq!(dc_grade(0.60), "丙");
        assert_eq!(dc_grade(0.30), "不合格");
    }

    #[test]
    fn test_compute_metrics() {
        let obs = vec![10.0, 20.0, 30.0, 40.0, 50.0];
        let m = compute_metrics(&obs, &obs);
        assert!((m.nse - 1.0).abs() < 1e-10);
        assert_eq!(m.grade, "甲");
        assert!(m.gbt.is_none() && m.qualified.is_none(), "compute_metrics 不应填 GB/T");
    }

    #[test]
    fn test_event_errors_perfect() {
        // obs==sim → 三项误差全 0;径流深 = ΣQ·1·3.6/1000 = 150·3.6/1000 = 0.54 mm
        let obs = vec![10.0, 20.0, 30.0, 40.0, 50.0];
        let e = event_errors(&obs, &obs, 1000.0, 1.0);
        assert!(e.peak_rel_err_pct.abs() < 1e-9);
        assert_eq!(e.time_to_peak_err_steps, 0);
        assert!(e.runoff_depth_rel_err_pct.abs() < 1e-9);
    }

    #[test]
    fn test_event_errors_peak_and_timing() {
        // obs 峰在第 4 步=50;sim 峰在第 2 步=45(超前 2 步,峰值偏小 10%)
        let obs = vec![10.0, 20.0, 30.0, 40.0, 50.0];
        let sim = vec![10.0, 20.0, 45.0, 30.0, 10.0];
        let e = event_errors(&obs, &sim, 1000.0, 1.0);
        assert!((e.peak_rel_err_pct - (-10.0)).abs() < 1e-6, "峰偏 -10%, 得 {}", e.peak_rel_err_pct);
        assert_eq!(e.time_to_peak_err_steps, -2, "超前 2 步");
    }

    #[test]
    fn test_event_qualified_boundary() {
        let tol = QualificationTolerance::default();
        // 三项都在限内 → 合格
        let ok = GbtEventErrors { peak_rel_err_pct: 19.0, time_to_peak_err_steps: 1, runoff_depth_rel_err_pct: -15.0 };
        assert!(event_qualified(&ok, &tol));
        // 峰超限 → 不合格
        let bad = GbtEventErrors { peak_rel_err_pct: 21.0, time_to_peak_err_steps: 0, runoff_depth_rel_err_pct: 0.0 };
        assert!(!event_qualified(&bad, &tol));
    }

    #[test]
    fn test_qualified_rate_and_grade() {
        // 9/10 合格 = 90% → 甲
        let flags = vec![true; 9].into_iter().chain(std::iter::once(false)).collect::<Vec<_>>();
        let r = qualified_rate(&flags);
        assert_eq!(r.total, 10);
        assert_eq!(r.qualified, 9);
        assert!((r.rate_pct - 90.0).abs() < 1e-9);
        assert_eq!(r.grade, "甲");
        assert_eq!(qualification_grade(75.0), "乙");
        assert_eq!(qualification_grade(65.0), "丙");
        assert_eq!(qualification_grade(50.0), "不合格");
    }

    #[test]
    fn test_scheme_grade_takes_lower() {
        assert_eq!(scheme_grade("甲", "乙"), "乙"); // 取低
        assert_eq!(scheme_grade("丙", "甲"), "丙");
        assert_eq!(scheme_grade("甲", "甲"), "甲");
        assert_eq!(scheme_grade("乙", "不合格"), "不合格");
    }

    #[test]
    fn test_compute_metrics_gbt_qualified() {
        let obs = vec![10.0, 20.0, 30.0, 40.0, 50.0];
        // sim 完全一致 → 合格
        let m = compute_metrics_gbt(&obs, &obs, 1000.0, 1.0, &Default::default());
        assert_eq!(m.qualified, Some(true));
        assert!(m.gbt.is_some());
    }

    #[test]
    fn test_crps_perfect_ensemble_zero() {
        // 成员全 = obs → CRPS = 0
        let obs = vec![10.0, 20.0, 30.0];
        let ens = vec![obs.clone(), obs.clone(), obs.clone()];
        assert!(crps(&ens, &obs).abs() < 1e-9);
    }

    #[test]
    fn test_crps_single_member_is_mae_and_spread_lowers() {
        let obs = vec![0.0];
        // 单成员 [10] → CRPS = MAE = 10(spread 项=0)
        assert!((crps(&[vec![10.0]], &obs) - 10.0).abs() < 1e-9);
        // 双成员 [0,10]:mae=(0+10)/2=5;spread=(1/8)·20=2.5 → CRPS=2.5
        assert!((crps(&[vec![0.0], vec![10.0]], &obs) - 2.5).abs() < 1e-9);
    }

    #[test]
    fn test_brier_perfect_and_worst() {
        let obs = vec![true, false, true];
        assert!(brier(&[1.0, 0.0, 1.0], &obs).abs() < 1e-9); // 完美
        assert!((brier(&[0.0, 1.0, 0.0], &obs) - 1.0).abs() < 1e-9); // 最差
    }

    #[test]
    fn test_pod_far_csi_contingency() {
        // 3 命中 + 1 漏报 + 2 误报(+ 2 正确拒绝)
        let fcst = vec![true, true, true, false, true, true, false, false];
        let obs = vec![true, true, true, true, false, false, false, false];
        let (pod, far, csi) = pod_far_csi(&fcst, &obs);
        assert!((pod - 0.75).abs() < 1e-9, "POD=3/4, got {}", pod);
        assert!((far - 0.4).abs() < 1e-9, "FAR=2/5, got {}", far);
        assert!((csi - 0.5).abs() < 1e-9, "CSI=3/6, got {}", csi);
    }

    #[test]
    fn test_ensemble_exceedance_fraction() {
        let ens = vec![vec![10.0, 20.0, 5.0], vec![15.0, 5.0, 5.0], vec![20.0, 25.0, 5.0]];
        let prob = ensemble_exceedance(&ens, 12.0);
        // t0:{10,15,20}→2/3; t1:{20,5,25}→2/3; t2:{5,5,5}→0
        assert!((prob[0] - 2.0 / 3.0).abs() < 1e-9);
        assert!((prob[1] - 2.0 / 3.0).abs() < 1e-9);
        assert!(prob[2].abs() < 1e-9);
    }

    #[test]
    fn test_compute_probabilistic_end_to_end() {
        // 2 成员,obs=[10,30];阈值=20。成员1=[12,28],成员2=[8,32]
        let ens = vec![vec![12.0, 28.0], vec![8.0, 32.0]];
        let obs = vec![10.0, 30.0];
        let r = compute_probabilistic(&ens, &obs, 20.0);
        // 超限概率:t0={12,8}→0(都≤20); t1={28,32}→1(都>20)
        // fcst 超限(>0.5):[false,true]; obs 超限(>20):[false,true] → 完美:POD=1,FAR=0,CSI=1
        assert!((r.pod - 1.0).abs() < 1e-9);
        assert!(r.far.abs() < 1e-9);
        assert!((r.csi - 1.0).abs() < 1e-9);
        // Brier:prob=[0,1],obs指示=[0,1] → 0
        assert!(r.brier.abs() < 1e-9);
    }

    #[test]
    fn test_rank_histogram_bins() {
        // 2 成员,t0:成员[1,3] obs=2 → rank=1(只有1<2);t1:成员[2,4] obs=3 → rank=1(只有2<3)
        let ens = vec![vec![1.0, 2.0], vec![3.0, 4.0]];
        let obs = vec![2.0, 3.0];
        let h = rank_histogram(&ens, &obs);
        assert_eq!(h, vec![0, 2, 0], "两步都落 bin1, got {:?}", h);
        // 观测恒低于所有成员 → 全落 bin0
        let h2 = rank_histogram(&vec![vec![10.0], vec![20.0]], &[5.0]);
        assert_eq!(h2, vec![1, 0, 0]);
    }
}