//! PIST Spectral — Rust Port
//!
//! Mirrors `formal/SilverSight/PIST/Spectral.lean`.
//! Minimal fixed-point spectral feature extraction:
//!   - isqrt (integer square root)
//!   - power iteration for dominant eigenvalue
//!   - SpectralProfile
//!   - Fiedler value via Laplacian

use crate::q16::*;

// ── Integer square root ────────────────────────────────────────────

/// Integer square root via Newton's method. Returns floor(√n).
pub fn isqrt(n: i64) -> i64 {
    if n <= 0 { return 0; }
    let mut x = n / 2 + 1;
    for _ in 0..64 {
        let x_new = (x + n / x) / 2;
        if x_new >= x { return x; }
        x = x_new;
    }
    x
}

// ── Matrix helpers ─────────────────────────────────────────────────

type IntMat = Vec<Vec<i64>>;

fn get_entry(mat: &IntMat, i: usize, j: usize) -> i64 {
    mat.get(i).and_then(|row| row.get(j).copied()).unwrap_or(0)
}

fn row_sum(mat: &IntMat, i: usize, n: usize) -> i64 {
    (0..n).map(|j| get_entry(mat, i, j)).sum()
}

fn symmetrize(mat: &IntMat, n: usize) -> IntMat {
    (0..n).map(|i| {
        (0..n).map(|j| {
            (get_entry(mat, i, j) + get_entry(mat, j, i)) / 2
        }).collect()
    }).collect()
}

fn build_laplacian(sym: &IntMat, n: usize) -> IntMat {
    (0..n).map(|i| {
        let deg = row_sum(sym, i, n);
        (0..n).map(|j| {
            if i == j { deg } else { -get_entry(sym, i, j) }
        }).collect()
    }).collect()
}

fn build_ata(mat: &IntMat, n: usize) -> IntMat {
    (0..n).map(|i| {
        (0..n).map(|j| {
            (0..n).map(|k| get_entry(mat, k, i) * get_entry(mat, k, j)).sum()
        }).collect()
    }).collect()
}

// ── SpectralProfile ────────────────────────────────────────────────

#[derive(Debug, Clone)]
pub struct SpectralProfile {
    pub dominant_eigenvalue: f64,
    pub fiedler_value: f64,
    pub spectral_gap: f64,
    pub condition_number: f64,
}

// ── Power iteration (f64 arithmetic) ───────────────────────────────

/// Dominant eigenvalue of a square Int matrix via power iteration.
pub fn power_iteration(mat: &IntMat, max_iter: usize) -> f64 {
    let n = mat.len();
    if n == 0 { return 0.0; }

    let mut v: Vec<f64> = (0..n).map(|i| (i as f64 + 1.0)).collect();

    for _ in 0..max_iter {
        // mat × v (as f64)
        let mv: Vec<f64> = (0..n).map(|i| {
            (0..n).map(|j| get_entry(mat, i, j) as f64 * v[j]).sum()
        }).collect();

        // Rayleigh quotient
        let v_dot_mv: f64 = v.iter().zip(mv.iter()).map(|(vi, mvi)| vi * mvi).sum();
        let v_dot_v: f64 = v.iter().map(|x| x * x).sum();
        if v_dot_v < 1e-15 { break; }

        // Copy for next iteration, normalized
        let norm = (mv.iter().map(|x| x * x).sum::<f64>()).sqrt();
        if norm < 1e-15 { break; }
        for (vi, mvi) in v.iter_mut().zip(mv.iter()) {
            *vi = mvi / norm;
        }

        // Check convergence
        let eig = v_dot_mv / v_dot_v;
        let mx = (0..n).map(|i| {
            (0..n).map(|j| get_entry(mat, i, j) as f64 * v[j]).sum::<f64>()
        }).collect::<Vec<_>>();
        let resid: f64 = mx.iter().zip(v.iter()).map(|(mxi, vi)| (mxi - eig * vi).abs()).sum::<f64>() / n as f64;
        if resid < 1e-8 { return eig; }
    }

    // Final Rayleigh quotient
    let mv: Vec<f64> = (0..n).map(|i| {
        (0..n).map(|j| get_entry(mat, i, j) as f64 * v[j]).sum()
    }).collect();
    let num: f64 = v.iter().zip(mv.iter()).map(|(vi, mvi)| vi * mvi).sum();
    let den: f64 = v.iter().map(|x| x * x).sum();
    if den > 0.0 { num / den } else { 0.0 }
}

fn rayleigh_quotient(mat: &IntMat, v: &[f64]) -> f64 {
    let n = mat.len();
    let mv: Vec<f64> = (0..n).map(|i| {
        (0..n).map(|j| get_entry(mat, i, j) as f64 * v[j]).sum()
    }).collect();
    let num: f64 = v.iter().zip(mv.iter()).map(|(vi, mvi)| vi * mvi).sum();
    let den: f64 = v.iter().map(|x| x * x).sum();
    if den > 0.0 { num / den } else { 0.0 }
}

/// Fiedler value (smallest non-zero eigenvalue of Laplacian) via shifted inverse iteration.
pub fn fiedler_value(mat: &IntMat) -> f64 {
    let n = mat.len();
    if n < 2 { return 0.0; }

    let sym = symmetrize(mat, n);
    let lap = build_laplacian(&sym, n);

    // Power iteration for dominant eigenvalue
    let lambda_max = power_iteration(&lap, 100);

    // Shift-invert: solve (L - μI)⁻¹ where μ = lambda_max * 0.9
    // to find the smallest eigenvalue near the upper end
    let mu = lambda_max * 0.9;

    let mut v: Vec<f64> = (0..n).map(|i| if i % 2 == 0 { 1.0 } else { -1.0 }).collect();

    for _ in 0..50 {
        // (L - μI) × v
        let mv: Vec<f64> = (0..n).map(|i| {
            let row: f64 = (0..n).map(|j| {
                let l_ij = if i == j {
                    row_sum(&lap, i, n) as f64
                } else {
                    -get_entry(&lap, i, j) as f64
                };
                l_ij * v[j]
            }).sum();
            row - mu * v[i]
        }).collect();

        // Normalize
        let norm = mv.iter().map(|x| x * x).sum::<f64>().sqrt();
        if norm < 1e-10 { break; }
        for vi in v.iter_mut() { *vi /= norm; }
    }

    // Compute Rayleigh quotient with the converged vector
    rayleigh_quotient(&lap, &v.iter().copied().collect::<Vec<_>>())
}

/// Full spectral profile for an n×n Int matrix.
pub fn compute_profile(mat: &IntMat) -> SpectralProfile {
    let n = mat.len();
    if n == 0 {
        return SpectralProfile {
            dominant_eigenvalue: 0.0, fiedler_value: 0.0,
            spectral_gap: 0.0, condition_number: 0.0,
        };
    }

    let sym = symmetrize(mat, n);
    let lap = build_laplacian(&sym, n);

    let lambda_max = power_iteration(&lap, 100);
    let fv = fiedler_value(mat);

    SpectralProfile {
        dominant_eigenvalue: lambda_max,
        fiedler_value: fv,
        spectral_gap: lambda_max - fv,
        condition_number: if fv.abs() > 1e-10 { lambda_max / fv } else { f64::INFINITY },
    }
}

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

    #[test]
    fn test_isqrt() {
        assert_eq!(isqrt(0), 0);
        assert_eq!(isqrt(1), 1);
        assert_eq!(isqrt(4), 2);
        assert_eq!(isqrt(9), 3);
        assert_eq!(isqrt(16), 4);
        assert_eq!(isqrt(2), 1);   // floor(√2)
        assert_eq!(isqrt(10), 3);  // floor(√10)
    }

    #[test]
    fn test_symmetrize() {
        let mat: IntMat = vec![
            vec![1, 2],
            vec![3, 4],
        ];
        let sym = symmetrize(&mat, 2);
        assert_eq!(sym[0][1], sym[1][0]);  // symmetric
        assert_eq!(sym[0][1], (2 + 3) / 2);
    }

    #[test]
    fn test_power_iteration_small() {
        // 2×2 identity → eigenvalue should be 1
        let mat: IntMat = vec![
            vec![1, 0],
            vec![0, 1],
        ];
        let eig = power_iteration(&mat, 100);
        assert!((eig - 1.0).abs() < 0.1);
    }

    #[test]
    fn test_spectral_profile() {
        // 3×3 matrix with known structure
        let mat: IntMat = vec![
            vec![2, 1, 0],
            vec![1, 2, 1],
            vec![0, 1, 2],
        ];
        let profile = compute_profile(&mat);
        assert!(profile.dominant_eigenvalue > 0.0);
        assert!(profile.spectral_gap >= 0.0);
    }
}