//! Risk functionals and scalar blow-up benchmark.
//!
//! Provenance: Lean 4 → LambdaIR → MiniC → Rust
//! Source: HeytingLean.Analysis.PMBoundedControl.RiskFunctional

use crate::completion_operator::sat_rational;

/// Safety penalty: Ψ(q) = 1/(1+q)
///
/// Theorem preserved: safety_penalty_nonneg
pub fn safety_penalty(q: f64) -> f64 {
    1.0 / (1.0 + q)
}

/// Tokamak risk proxy: β_local * Ψ(q) * (1 + T)
///
/// Theorem preserved: tokamak_risk_nonneg
/// Theorem preserved: tokamak_risk_monotone_beta
pub fn tokamak_risk(beta_local: f64, q: f64, t: f64) -> f64 {
    beta_local * safety_penalty(q) * (1.0 + t)
}

/// PM boundary
pub struct PmBoundary {
    pub q_pm: f64,
}

/// Check PM-bounded invariant along a trajectory
///
/// Theorem preserved: pm_invariant_of_completion
pub fn pm_bounded_invariant(risks: &[f64], q_pm: f64) -> bool {
    risks.iter().all(|&r| r <= q_pm)
}

/// Classical blow-up profile: u₀ / √(1 - 2u₀²t)
pub fn classical_blowup_profile(u0: f64, t: f64) -> f64 {
    let denom = 1.0 - 2.0 * u0 * u0 * t;
    if denom <= 0.0 {
        return f64::INFINITY;
    }
    u0 / denom.sqrt()
}

/// Blow-up time: 1 / (2u₀²)
pub fn classical_blowup_time(u0: f64) -> f64 {
    1.0 / (2.0 * u0 * u0)
}

/// PM-completed scalar profile
///
/// Theorem preserved: scalar_bounded_completed
pub fn completed_scalar_profile(q_pm: f64, u0: f64, t: f64) -> f64 {
    sat_rational(q_pm, classical_blowup_profile(u0, t))
}

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

    #[test]
    fn test_tokamak_risk_nonneg() {
        assert!(tokamak_risk(1.0, 2.0, 0.5) >= 0.0);
    }

    #[test]
    fn test_tokamak_risk_monotone_beta() {
        let r1 = tokamak_risk(1.0, 2.0, 0.5);
        let r2 = tokamak_risk(2.0, 2.0, 0.5);
        assert!(r1 <= r2 + 1e-12);
    }

    #[test]
    fn test_completed_bounded() {
        let q = 5.0;
        let u0 = 2.0;
        // Test in the pre-blow-up regime (Lean theorem covers finite reals)
        let t_star = classical_blowup_time(u0);
        for &t in &[0.0, 0.5 * t_star, 0.9 * t_star, 0.99 * t_star] {
            assert!(completed_scalar_profile(q, u0, t).abs() <= q + 1e-12);
        }
    }
}