//! Connes–Rovelli thermal time.
//!
//! The thermal-time hypothesis: physical time flow is generated by the
//! statistical state itself, not imposed from outside. Given a density matrix
//! `ρ`, the **modular Hamiltonian** is
//!
//! ```text
//!   K = -ln ρ
//! ```
//!
//! and observables flow under modular time `s`:
//!
//! ```text
//!   A(s) = e^{isK} A e^{-isK},      dA/ds = i[K, A].
//! ```
//!
//! For a Gibbs state `ρ = e^{-βH}/Z` we get `K = βH + (ln Z)·I`, so modular flow
//! *is* ordinary Hamiltonian evolution rescaled by `β`: physical time is
//! recovered from the thermodynamic state.

use crate::complex::Complex;
use crate::complex_matrix::{exp_i_symmetric, CMatrix};
use crate::real_matrix::RealMatrix;

const PROB_FLOOR: f64 = 1e-12;

/// Modular Hamiltonian `K = -ln ρ` for a real symmetric density matrix.
/// Eigenvalues below `PROB_FLOOR` are clamped to keep the logarithm finite.
pub fn modular_hamiltonian(rho: &RealMatrix) -> RealMatrix {
    let (probs, vecs) = rho.symmetric_eigen();
    let k_eigs: Vec<f64> = probs.iter().map(|&p| -(p.max(PROB_FLOOR)).ln()).collect();
    RealMatrix::from_spectrum(&k_eigs, &vecs)
}

/// Modular flow `A(s) = e^{isK} A e^{-isK}`.
pub fn modular_flow(k: &RealMatrix, a: &CMatrix, s: f64) -> CMatrix {
    let u = exp_i_symmetric(k, s);
    let u_dag = exp_i_symmetric(k, -s);
    u.matmul(a).matmul(&u_dag)
}

/// Modular flow generator `dA/ds = i[K, A]` at `s = 0`.
pub fn modular_generator(k: &RealMatrix, a: &CMatrix) -> CMatrix {
    let kc = CMatrix::from_real(k);
    CMatrix::commutator(&kc, a).scale(Complex::I)
}

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

    fn sample_h() -> RealMatrix {
        RealMatrix::from_fn(3, |r, c| if r == c { r as f64 } else { 0.2 })
    }

    fn sample_observable() -> CMatrix {
        // A Hermitian observable.
        let mut a = CMatrix::zeros(3);
        a.set(0, 1, Complex::new(0.0, 1.0));
        a.set(1, 0, Complex::new(0.0, -1.0));
        a.set(2, 2, Complex::real(0.5));
        a
    }

    #[test]
    fn modular_hamiltonian_of_gibbs_is_beta_h() {
        let h = sample_h();
        let beta = 0.7;
        let rho = gibbs_density(&h, beta);
        let k = modular_hamiltonian(&rho);
        // K should equal βH + cI for some constant c. Check K - βH is a multiple
        // of the identity by verifying all diagonal-shifted off-diagonals match.
        let diff = RealMatrix::from_fn(3, |r, c| k.get(r, c) - beta * h.get(r, c));
        // off-diagonal entries of diff ≈ 0
        assert!(diff.max_offdiag() < 1e-7);
        // diagonal entries of diff are all equal (the additive constant)
        let c0 = diff.get(0, 0);
        for i in 1..3 {
            assert!((diff.get(i, i) - c0).abs() < 1e-7);
        }
    }

    #[test]
    fn generator_matches_finite_difference() {
        let h = sample_h();
        let rho = gibbs_density(&h, 0.9);
        let k = modular_hamiltonian(&rho);
        let a = sample_observable();

        let gen = modular_generator(&k, &a);
        let ds = 1e-5;
        let fwd = modular_flow(&k, &a, ds);
        let bwd = modular_flow(&k, &a, -ds);
        // central difference of A(s)
        let fd = fwd.sub(&bwd).scale(Complex::real(1.0 / (2.0 * ds)));
        assert!(gen.sub(&fd).frob_norm() < 1e-4);
    }

    #[test]
    fn modular_flow_equals_rescaled_physical_evolution() {
        // For a Gibbs state, modular flow at parameter s == physical evolution
        // e^{iβHs} A e^{-iβHs} (the ln Z phase cancels).
        let h = sample_h();
        let beta = 0.6;
        let rho = gibbs_density(&h, beta);
        let k = modular_hamiltonian(&rho);
        let a = sample_observable();
        let s = 0.8;

        let modular = modular_flow(&k, &a, s);
        let u = exp_i_symmetric(&h, beta * s);
        let u_dag = exp_i_symmetric(&h, -beta * s);
        let physical = u.matmul(&a).matmul(&u_dag);

        assert!(modular.sub(&physical).frob_norm() < 1e-7);
    }
}
