//! Closed-form validation of the generalized eigenvalue solver.
//!
//! `ModalAnalysis` is the intended consumer of this solver, so the solver has
//! to be right before the analysis is wired to it. Every expected value here
//! is derived analytically and stated in the test, not lifted from a previous
//! run of this code.
//!
//! The cases build up deliberately:
//!
//! 1. `M = I` with a matrix whose spectrum has a closed form — isolates the
//!    Lanczos recurrence and the tridiagonal solve.
//! 2. `M != I` — exercises the generalized path `M^-1 K`, which is where a
//!    mass matrix that is merely assumed to be the identity would hide.
//! 3. `n` much larger than the requested mode count — the case a real modal
//!    analysis always hits, and the one that exercises the back-transformation
//!    of eigenvectors out of the Krylov basis.

use approx::assert_relative_eq;
use rtx_fea::assembly::SparseMatrix;
use rtx_fea::solvers::eigenvalue::EigenvalueSolver;

/// Build a `SparseMatrix` from a dense row-major listing.
fn sparse_from_rows(rows: &[&[f64]]) -> SparseMatrix {
    let n = rows.len();
    let mut m = SparseMatrix::new(n, rows[0].len());
    for (i, row) in rows.iter().enumerate() {
        for (j, &v) in row.iter().enumerate() {
            if v != 0.0 {
                m.set_entry(i, j, v).unwrap();
            }
        }
    }
    assert_eq!(m.nrows(), n);
    m
}

/// `K = tridiag(-1, 4, -1)` of order 3, `M = I`.
///
/// The eigenvalues of an `n x n` matrix `tridiag(c, a, c)` are
/// `a + 2c cos(k pi / (n + 1))` for `k = 1..=n`. With `a = 4`, `c = -1`,
/// `n = 3` and `cos(pi/4) = cos(3pi/4).abs() = sqrt(2)/2`, that gives
/// `4 - sqrt(2)`, `4`, `4 + sqrt(2)`.
#[test]
fn tridiagonal_spectrum_matches_closed_form() {
    let k = sparse_from_rows(&[&[4.0, -1.0, 0.0], &[-1.0, 4.0, -1.0], &[0.0, -1.0, 4.0]]);
    let m = sparse_from_rows(&[&[1.0, 0.0, 0.0], &[0.0, 1.0, 0.0], &[0.0, 0.0, 1.0]]);

    let (eigenvalues, _) = EigenvalueSolver::new(3).solve(&k, &m).unwrap();

    let sqrt2 = 2.0_f64.sqrt();
    let expected = [4.0 - sqrt2, 4.0, 4.0 + sqrt2];

    assert_eq!(eigenvalues.len(), 3);
    for (i, &want) in expected.iter().enumerate() {
        assert_relative_eq!(eigenvalues[i], want, epsilon = 1e-9);
    }
}

/// Two-degree-of-freedom spring-mass chain, fixed at one end.
///
/// `K = [[2, -1], [-1, 1]]`, `M = diag(2, 1)`. Then
/// `det(K - lambda M) = 2 lambda^2 - 4 lambda + 1`, whose roots are
/// `1 -+ sqrt(2)/2`. `M` is deliberately not the identity: if the solver
/// ignored the mass matrix it would return the eigenvalues of `K` alone,
/// which are `(3 -+ sqrt(5))/2` and would fail this assertion.
#[test]
fn generalized_problem_accounts_for_the_mass_matrix() {
    let k = sparse_from_rows(&[&[2.0, -1.0], &[-1.0, 1.0]]);
    let m = sparse_from_rows(&[&[2.0, 0.0], &[0.0, 1.0]]);

    let (eigenvalues, _) = EigenvalueSolver::new(2).solve(&k, &m).unwrap();

    let half_sqrt2 = 2.0_f64.sqrt() / 2.0;
    assert_relative_eq!(eigenvalues[0], 1.0 - half_sqrt2, epsilon = 1e-9);
    assert_relative_eq!(eigenvalues[1], 1.0 + half_sqrt2, epsilon = 1e-9);
}

/// Eigenvalues must come back in ascending order.
///
/// Modal analysis names modes by index — mode 1 is the fundamental — so the
/// ordering is part of the contract, not a convenience. `nalgebra`'s
/// `symmetric_eigen` gives no ordering guarantee, so this has to be imposed.
#[test]
fn eigenvalues_are_returned_in_ascending_order() {
    let k = sparse_from_rows(&[&[4.0, -1.0, 0.0], &[-1.0, 4.0, -1.0], &[0.0, -1.0, 4.0]]);
    let m = sparse_from_rows(&[&[1.0, 0.0, 0.0], &[0.0, 1.0, 0.0], &[0.0, 0.0, 1.0]]);

    let (eigenvalues, _) = EigenvalueSolver::new(3).solve(&k, &m).unwrap();

    for i in 1..eigenvalues.len() {
        assert!(
            eigenvalues[i] >= eigenvalues[i - 1],
            "eigenvalue {} ({}) is smaller than eigenvalue {} ({})",
            i,
            eigenvalues[i],
            i - 1,
            eigenvalues[i - 1]
        );
    }
}

/// Ask for fewer modes than the problem has degrees of freedom.
///
/// This is the normal case — a modal analysis wants the lowest handful of an
/// otherwise large system — and it is the case that exercises the
/// back-transformation from the Krylov basis to the original space, where the
/// Lanczos vector block and the tridiagonal eigenvector block must agree in
/// dimension.
///
/// `K = tridiag(-1, 2, -1)` of order 10 has eigenvalues
/// `2 - 2 cos(k pi / 11)`, so the two smallest are `2 - 2 cos(pi/11)` and
/// `2 - 2 cos(2 pi / 11)`.
#[test]
fn requesting_a_subset_of_modes_returns_the_lowest_ones() {
    const N: usize = 10;

    let mut k = SparseMatrix::new(N, N);
    let mut m = SparseMatrix::new(N, N);
    for i in 0..N {
        k.set_entry(i, i, 2.0).unwrap();
        if i > 0 {
            k.set_entry(i, i - 1, -1.0).unwrap();
            k.set_entry(i - 1, i, -1.0).unwrap();
        }
        m.set_entry(i, i, 1.0).unwrap();
    }

    let num_modes = 2;
    let (eigenvalues, eigenvectors) = EigenvalueSolver::new(num_modes).solve(&k, &m).unwrap();

    assert_eq!(eigenvalues.len(), num_modes);
    assert_eq!(
        eigenvectors.nrows(),
        N,
        "eigenvectors must be returned in the original space, not the Krylov basis"
    );
    assert_eq!(eigenvectors.ncols(), num_modes);

    let pi = std::f64::consts::PI;
    for (i, want) in [
        2.0 - 2.0 * (pi / 11.0).cos(),
        2.0 - 2.0 * (2.0 * pi / 11.0).cos(),
    ]
    .into_iter()
    .enumerate()
    {
        assert_relative_eq!(eigenvalues[i], want, epsilon = 1e-8);
    }

    for i in 0..num_modes {
        assert_relative_eq!(eigenvectors.column(i).norm(), 1.0, epsilon = 1e-10);
    }
}

/// `f = sqrt(lambda) / (2 pi)`, not `sqrt(lambda / (2 pi))`.
///
/// `lambda` is `omega^2` in rad^2/s^2, so `omega = sqrt(lambda)` and
/// `f = omega / (2 pi)`. The two forms happen to agree only at
/// `lambda = 2 pi`, which is why a smoke test would not catch the difference.
#[test]
fn frequencies_convert_from_eigenvalues_correctly() {
    use nalgebra::{DMatrix, DVector};
    use rtx_fea::solvers::eigenvalue::ModalResults;

    let eigenvalues = DVector::from_vec(vec![1.0, 4.0, 9.0]);
    let modal = ModalResults::from_eigenvalues(eigenvalues, DMatrix::identity(3, 3));

    let two_pi = 2.0 * std::f64::consts::PI;
    for (i, omega) in [1.0_f64, 2.0, 3.0].into_iter().enumerate() {
        assert_relative_eq!(modal.frequencies[i], omega / two_pi, epsilon = 1e-12);
        assert_relative_eq!(modal.periods[i], two_pi / omega, epsilon = 1e-12);
    }
}