//! Multi-limit vector completion with SafEDMD integration.
//!
//! Provenance: Lean 4 → LambdaIR → MiniC → Rust
//! Source: HeytingLean.Analysis.PMBoundedControl.MultiLimit

use crate::completion_operator::sat_rational;

/// Sum of squares energy proxy
pub fn sq_sum(v: &[f64]) -> f64 {
    v.iter().map(|x| x * x).sum()
}

/// Euclidean norm
///
/// Theorem preserved: euclideanNorm_nonneg
pub fn euclidean_norm(v: &[f64]) -> f64 {
    sq_sum(v).sqrt()
}

/// Componentwise PM completion
pub fn componentwise_completion(q_pm: f64, v: &[f64]) -> Vec<f64> {
    v.iter().map(|&x| sat_rational(q_pm, x)).collect()
}

/// Norm-bounded completion by radial projection
///
/// Theorem preserved: sat_norm_bounded_le
/// Theorem preserved: sat_norm_bounded_identity
pub fn sat_norm_bounded(q_pm: f64, v: &[f64]) -> Vec<f64> {
    let norm = euclidean_norm(v);
    if norm <= q_pm {
        v.to_vec()
    } else {
        let scale = q_pm / norm;
        v.iter().map(|&x| scale * x).collect()
    }
}

/// ReLU nucleus applied componentwise
pub fn relu_nucleus(v: &[f64]) -> Vec<f64> {
    v.iter().map(|&x| x.max(0.0)).collect()
}

/// Componentwise completion implies norm bound: ‖v‖ ≤ √n * Q_PM
///
/// Theorem preserved: componentwise_implies_norm
pub fn componentwise_implies_norm_bound(n: usize, q_pm: f64) -> f64 {
    (n as f64).sqrt() * q_pm
}

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

    #[test]
    fn test_sat_norm_bounded_le() {
        let q = 5.0;
        let v = vec![3.0, 4.0, 5.0, 6.0];
        let out = sat_norm_bounded(q, &v);
        assert!(euclidean_norm(&out) <= q + 1e-12);
    }

    #[test]
    fn test_sat_norm_bounded_identity() {
        let q = 100.0;
        let v = vec![1.0, 2.0, 3.0];
        let out = sat_norm_bounded(q, &v);
        for (a, b) in v.iter().zip(out.iter()) {
            assert!((a - b).abs() < 1e-12);
        }
    }

    #[test]
    fn test_componentwise_sqsum_le() {
        let q = 5.0;
        let v = vec![10.0, -20.0, 30.0, -40.0];
        let comp = componentwise_completion(q, &v);
        let ss = sq_sum(&comp);
        let bound = (v.len() as f64) * q * q;
        assert!(ss <= bound + 1e-12);
    }

    #[test]
    fn test_relu_after_componentwise() {
        let q = 5.0;
        let v = vec![10.0, -20.0, 30.0, -40.0];
        let comp = componentwise_completion(q, &v);
        let relu_comp = relu_nucleus(&comp);
        let ss = sq_sum(&relu_comp);
        let bound = (v.len() as f64) * q * q;
        assert!(ss <= bound + 1e-12);
    }
}