lapl 0.2.0

//! Spectral analysis of graph Laplacians using Random Matrix Theory.
//!
//! Computes the Laplacian of random and structured graphs, extracts eigenvalues,
//! and uses RMT-style spacing and density summaries alongside graph-specific
//! low-eigenvalue diagnostics.
//!
//! Key idea: for normalized graph Laplacians, community structure appears near
//! the bottom of the spectrum. A graph with k weakly connected communities has
//! k small eigenvalues followed by an eigengap.
//!
//! Run: cargo run --example spectral_rmt

use lapl::{
    adjacency_to_laplacian, gaussian_similarity, normalized_laplacian, symmetric_eigenvalues,
};
use ndarray::Array2;
use rmt::{empirical_spectral_density, level_spacing_ratios, mean_spacing_ratio};

fn main() -> lapl::Result<()> {
    println!("=== Graph Laplacian Spectral Analysis via RMT ===\n");

    // 1. Random graph (Erdos-Renyi style via Gaussian similarity)
    let n = 60;
    let random_points = random_points_2d(n, 42);
    let adj_random = gaussian_similarity(&random_points, 1.0);
    let lap_random = normalized_laplacian(&adj_random);
    let eigs_random = symmetric_eigenvalues(&lap_random, 1e-12, 100 * n * n)?;

    // 2. Structured graph (3 well-separated clusters)
    let structured_points = clustered_points_2d(20, 3, 5.0, 42);
    let adj_struct = gaussian_similarity(&structured_points, 1.0);
    let lap_struct = normalized_laplacian(&adj_struct);
    let eigs_struct = symmetric_eigenvalues(&lap_struct, 1e-12, 100 * n * n)?;

    // Analyze with RMT
    println!("--- Random graph ({n} nodes, Gaussian similarity) ---\n");
    analyze_spectrum("Random", &eigs_random);

    println!("\n--- Structured graph (3 clusters of 20, separation=5.0) ---\n");
    analyze_spectrum("Structured", &eigs_struct);

    // Compare unnormalized Laplacian
    let lap_unnorm = adjacency_to_laplacian(&adj_struct);
    let eigs_unnorm = symmetric_eigenvalues(&lap_unnorm, 1e-12, 100 * n * n)?;
    println!("\n--- Structured (unnormalized Laplacian) ---\n");
    analyze_spectrum("Unnormalized", &eigs_unnorm);

    Ok(())
}

fn analyze_spectrum(label: &str, eigenvalues: &[f64]) {
    let n = eigenvalues.len();

    // Spacing ratio: GOE ~ 0.5307, Poisson ~ 0.3863
    let msr = mean_spacing_ratio(eigenvalues);
    let ratios = level_spacing_ratios(eigenvalues);
    let regime = if msr > 0.48 {
        "GOE (correlated/repulsive)"
    } else if msr < 0.42 {
        "Poisson (uncorrelated/independent)"
    } else {
        "intermediate"
    };

    println!("  {label}: {n} eigenvalues");
    println!(
        "  Range: [{:.4}, {:.4}]",
        eigenvalues[0],
        eigenvalues[n - 1]
    );
    println!("  Mean spacing ratio: {msr:.4} -> {regime}");

    // Count near-zero eigenvalues (connected components indicator)
    let near_zero = eigenvalues.iter().filter(|&&e| e.abs() < 1e-6).count();
    println!("  Near-zero eigenvalues: {near_zero} (= connected components)");

    let small_threshold = 0.1;
    let small_count = eigenvalues.iter().filter(|&&e| e < small_threshold).count();
    println!("  Eigenvalues < {small_threshold:.1}: {small_count}");
    println!(
        "  First eigenvalues: {}",
        format_first_eigenvalues(eigenvalues, 8)
    );

    if let Some((idx, gap)) = largest_low_end_gap(eigenvalues, 8) {
        println!(
            "  Largest low-end gap: lambda_{idx}->lambda_{} = {gap:.4}",
            idx + 1
        );
    }

    // MP outlier count is a rough high-end diagnostic, not the community count.
    let high_end_outliers = rmt::effective_dimension(eigenvalues, n, n);
    println!("  MP high-end outliers: {high_end_outliers}");

    // Histogram summary
    let (centers, densities) = empirical_spectral_density(eigenvalues, 10);
    println!("  Spectral density (10 bins):");
    for (c, d) in centers.iter().zip(densities.iter()) {
        let bar_len = (d * 20.0) as usize;
        let bar = "#".repeat(bar_len.min(40));
        println!("    {c:6.3} | {d:5.2} {bar}");
    }

    // Spacing ratio distribution summary
    if ratios.len() >= 5 {
        let mut sorted_ratios = ratios.clone();
        sorted_ratios.sort_by(|a, b| a.partial_cmp(b).unwrap());
        let q1 = sorted_ratios[sorted_ratios.len() / 4];
        let median = sorted_ratios[sorted_ratios.len() / 2];
        let q3 = sorted_ratios[3 * sorted_ratios.len() / 4];
        println!("  Spacing ratio quartiles: Q1={q1:.3} med={median:.3} Q3={q3:.3}");
    }
}

fn format_first_eigenvalues(eigenvalues: &[f64], limit: usize) -> String {
    eigenvalues
        .iter()
        .take(limit)
        .map(|e| format!("{e:.4}"))
        .collect::<Vec<_>>()
        .join(", ")
}

fn largest_low_end_gap(eigenvalues: &[f64], limit: usize) -> Option<(usize, f64)> {
    let upper = eigenvalues.len().min(limit);
    if upper < 2 {
        return None;
    }

    (0..upper - 1)
        .map(|idx| (idx, eigenvalues[idx + 1] - eigenvalues[idx]))
        .max_by(|(_, a), (_, b)| a.partial_cmp(b).unwrap())
}

fn random_points_2d(n: usize, seed: u64) -> Array2<f64> {
    // Simple LCG for reproducible random points
    let mut state = seed;
    let mut points = Array2::zeros((n, 2));
    for i in 0..n {
        for j in 0..2 {
            state = state.wrapping_mul(6364136223846793005).wrapping_add(1);
            points[[i, j]] = (state >> 33) as f64 / (1u64 << 31) as f64;
        }
    }
    points
}

fn clustered_points_2d(
    per_cluster: usize,
    n_clusters: usize,
    separation: f64,
    seed: u64,
) -> Array2<f64> {
    let n = per_cluster * n_clusters;
    let mut points = Array2::zeros((n, 2));
    let mut state = seed;
    for c in 0..n_clusters {
        let cx = (c as f64) * separation;
        let cy = 0.0;
        for i in 0..per_cluster {
            let idx = c * per_cluster + i;
            for j in 0..2 {
                state = state.wrapping_mul(6364136223846793005).wrapping_add(1);
                let noise = ((state >> 33) as f64 / (1u64 << 31) as f64 - 0.5) * 0.5;
                points[[idx, j]] = if j == 0 { cx } else { cy } + noise;
            }
        }
    }
    points
}