//! Statistical micro-library for the tau-Epoch framework.
//!
//! CAB-certified extraction from:
//!   HeytingLean.Bridge.AlMayahi.TauEpoch.Stats
//!
//! Provenance: Lean 4 kernel-verified -> LambdaIR -> MiniC -> C -> Rust

/// Compute the weighted mean of `values` with uncertainties `sigmas`.
///
/// Theorem preserved: `weighted_mean_in_range`
pub fn weighted_mean(values: &[f64], sigmas: &[f64]) -> f64 {
    let (num, den) = values.iter().zip(sigmas.iter()).fold((0.0, 0.0), |(n, d), (&v, &s)| {
        let w = 1.0 / (s * s);
        (n + w * v, d + w)
    });
    if den > 0.0 { num / den } else { 0.0 }
}

/// Chi-squared statistic for a set of measurements.
pub fn chi_squared(values: &[f64], sigmas: &[f64]) -> f64 {
    let wm = weighted_mean(values, sigmas);
    values.iter().zip(sigmas.iter()).map(|(&v, &s)| {
        let d = (v - wm) / s;
        d * d
    }).sum()
}

/// Birge ratio: sqrt(chi2 / (n-1)).
///
/// Theorem preserved: `birge_one_iff_consistent`
pub fn birge_ratio(values: &[f64], sigmas: &[f64]) -> f64 {
    let n = values.len();
    if n <= 1 { return 0.0; }
    (chi_squared(values, sigmas) / (n - 1) as f64).sqrt()
}

/// Mid-rank tie-handling rank assignment.
pub fn ranks(xs: &[f64]) -> Vec<f64> {
    xs.iter().map(|&x| {
        let lt = xs.iter().filter(|&&v| v < x).count();
        let eq = xs.iter().filter(|&&v| (v - x).abs() < f64::EPSILON).count();
        lt as f64 + (eq as f64 + 1.0) / 2.0
    }).collect()
}

fn mean(xs: &[f64]) -> f64 {
    xs.iter().sum::<f64>() / xs.len() as f64
}

fn variance(xs: &[f64]) -> f64 {
    let m = mean(xs);
    xs.iter().map(|&x| (x - m) * (x - m)).sum::<f64>() / xs.len() as f64
}

fn covariance(xs: &[f64], ys: &[f64]) -> f64 {
    let mx = mean(xs);
    let my = mean(ys);
    xs.iter().zip(ys.iter()).map(|(&x, &y)| (x - mx) * (y - my)).sum::<f64>() / xs.len() as f64
}

/// Pearson product-moment correlation coefficient.
pub fn pearson(xs: &[f64], ys: &[f64]) -> f64 {
    let sx = variance(xs).sqrt();
    let sy = variance(ys).sqrt();
    if sx == 0.0 || sy == 0.0 { return 0.0; }
    (covariance(xs, ys) / (sx * sy)).clamp(-1.0, 1.0)
}

/// Spearman rank correlation (Pearson on ranks).
pub fn spearman(xs: &[f64], ys: &[f64]) -> f64 {
    pearson(&ranks(xs), &ranks(ys))
}

/// Linear regression: returns (slope, intercept, r_squared).
///
/// Theorem preserved: `regression_residual_orthogonal`
pub fn linear_regression(xs: &[f64], ys: &[f64]) -> (f64, f64, f64) {
    let vx = variance(xs);
    if vx == 0.0 { return (0.0, 0.0, 0.0); }

    let slope = covariance(xs, ys) / vx;
    let intercept = mean(ys) - slope * mean(xs);

    let my = mean(ys);
    let ss_res: f64 = xs.iter().zip(ys.iter())
        .map(|(&x, &y)| { let r = y - (slope * x + intercept); r * r })
        .sum();
    let ss_tot: f64 = ys.iter().map(|&y| (y - my) * (y - my)).sum();
    let r_sq = if ss_tot > 0.0 { 1.0 - ss_res / ss_tot } else { 0.0 };

    (slope, intercept, r_sq)
}

/// Binomial upper-tail p-value: P(X >= k) under Binomial(n, 0.5).
pub fn binomial_sign_pvalue(positive_count: usize, total: usize) -> f64 {
    let mut num: u64 = 0;
    for i in positive_count..=total {
        num += choose(total, i);
    }
    num as f64 / (1u64 << total) as f64
}

fn choose(n: usize, k: usize) -> u64 {
    if k > n { return 0; }
    let k = k.min(n - k);
    (0..k).fold(1u64, |acc, i| acc * (n - i) as u64 / (i + 1) as u64)
}

/// Two-Clock Projection Operator: Pi(x) = 1 + beta*(x-1) + gamma*(x-1)^2
///
/// Theorem preserved: `projection_eval_one`
pub fn projection_eval(beta: f64, gamma: f64, x: f64) -> f64 {
    let d = x - 1.0;
    1.0 + beta * d + gamma * d * d
}

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

    #[test]
    fn test_weighted_mean() {
        let vals = [67.4, 73.04];
        let sigs = [0.5, 1.04];
        let wm = weighted_mean(&vals, &sigs);
        // w1 = 1/(0.5^2) = 4.0, w2 = 1/(1.04^2) ~ 0.925
        // wm ~ (4.0*67.4 + 0.925*73.04) / (4.0 + 0.925) ~ 68.46
        assert!((wm - 68.46).abs() < 0.1);
    }

    #[test]
    fn test_pearson_perfect() {
        let xs = [1.0, 2.0, 3.0, 4.0, 5.0];
        let ys = [2.0, 4.0, 6.0, 8.0, 10.0];
        assert!((pearson(&xs, &ys) - 1.0).abs() < 1e-10);
    }

    #[test]
    fn test_projection_at_one() {
        assert_eq!(projection_eval(5.0, 3.0, 1.0), 1.0);
    }

    #[test]
    fn test_binomial() {
        // 7/7 positive: p = 1/128 = 0.0078125
        let p = binomial_sign_pvalue(7, 7);
        assert!((p - 0.0078125).abs() < 1e-10);
    }

    #[test]
    fn test_ranks_with_ties() {
        let xs = [3.0, 1.0, 1.0, 2.0];
        let r = ranks(&xs);
        assert_eq!(r[0], 4.0); // 3.0 is largest
        assert_eq!(r[1], 1.5); // tied 1.0s get mid-rank
        assert_eq!(r[2], 1.5);
        assert_eq!(r[3], 3.0); // 2.0 is third
    }
}