//! `cohomology` — la única cuenta que importa (§6.4, §7).
//!
//! `rank δ`, `dim H⁰`, `dim H¹` y una base del cokernel para localizar el nudo.

use crate::error::Result;
use crate::linalg::{F64Backend, LinAlg, Scalar};
use crate::sheaf::Sheaf;

/// Resultado del cálculo cohomológico sobre un haz.
#[derive(Debug, Clone, PartialEq)]
pub struct Cohomology {
    /// `dim C⁰ = Σ dim F(v)`.
    pub dim_c0: usize,
    /// `dim C¹ = Σ dim F(e)`.
    pub dim_c1: usize,
    /// `rank δ`.
    pub rank_delta: usize,
    /// `dim H⁰ = dim_c0 − rank_delta` (secciones globales coherentes).
    pub dim_h0: usize,
    /// `dim H¹ = dim_c1 − rank_delta` (obstrucciones).
    pub dim_h1: usize,
    /// Si `dim_h1 > 0`, una base del cokernel: cada vector señala aristas
    /// que participan en la obstrucción (para localizar el nudo).
    pub obstruction_basis: Vec<Vec<Scalar>>,
}

/// Calcula la cohomología del haz con el backend `f64` por defecto (§7.1).
pub fn compute(sheaf: &Sheaf) -> Result<Cohomology> {
    compute_with(sheaf, &F64Backend::new())
}

/// Calcula la cohomología del haz con un backend de álgebra lineal explícito.
///
/// Pasos (§7): ensamblar `δ`, `rank(δ)`, `dim H⁰ = dimC⁰ − rank`,
/// `dim H¹ = dimC¹ − rank`, y si `H¹ ≠ 0`, una base del cokernel.
pub fn compute_with<L: LinAlg>(sheaf: &Sheaf, backend: &L) -> Result<Cohomology> {
    let delta = sheaf.coboundary()?;
    let dim_c0 = sheaf.dim_c0();
    let dim_c1 = sheaf.dim_c1();
    let rank_delta = backend.rank(&delta);

    // rank ≤ min(dimC⁰, dimC¹), así que las restas nunca subdesbordan.
    let dim_h0 = dim_c0 - rank_delta;
    let dim_h1 = dim_c1 - rank_delta;

    let obstruction_basis = if dim_h1 > 0 {
        backend.cokernel_basis(&delta)
    } else {
        Vec::new()
    };

    Ok(Cohomology {
        dim_c0,
        dim_c1,
        rank_delta,
        dim_h0,
        dim_h1,
        obstruction_basis,
    })
}

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

    /// Prueba de cordura del haz constante (§7.2): para stalk de dimensión `d`,
    /// `dim H⁰ = (#componentes)·d` y `dim H¹ = (primer nº de Betti)·d`.
    ///
    /// Aquí verificamos con `d = 1`, donde H⁰ = #componentes y H¹ = b₁. Además
    /// cross-validamos que el backend exacto `GF(2)` coincide con `f64` en las
    /// dimensiones (sin epsilon): si difirieran, uno de los dos estaría mal.
    fn coho_constante(nv: usize, edges: &[(usize, usize)]) -> Cohomology {
        let sheaf = Sheaf::constant(nv, edges, 1);
        let f64_coho = compute(&sheaf).expect("haz constante bien formado");
        let gf2_coho = compute_with(&sheaf, &Gf2Backend).expect("haz constante");
        assert_eq!(
            (f64_coho.dim_h0, f64_coho.dim_h1),
            (gf2_coho.dim_h0, gf2_coho.dim_h1),
            "f64 y GF(2) deben coincidir en las dimensiones cohomológicas"
        );
        f64_coho
    }

    #[test]
    fn arbol_no_tiene_obstruccion() {
        // Camino de 4 vértices (árbol conexo): H⁰ = 1, H¹ = 0.
        let coho = coho_constante(4, &[(0, 1), (1, 2), (2, 3)]);
        assert_eq!(coho.dim_h0, 1, "un árbol es conexo");
        assert_eq!(coho.dim_h1, 0, "un árbol no tiene ciclos");
        assert!(coho.obstruction_basis.is_empty());
    }

    #[test]
    fn triangulo_tiene_un_ciclo() {
        // Triángulo: conexo con un ciclo. H⁰ = 1, H¹ = 1.
        let coho = coho_constante(3, &[(0, 1), (1, 2), (2, 0)]);
        assert_eq!(coho.dim_h0, 1);
        assert_eq!(coho.dim_h1, 1, "b₁ del triángulo es 1");
        assert_eq!(coho.obstruction_basis.len(), 1);
        // La obstrucción toca las 3 aristas del ciclo (soporte no nulo en todas).
        let v = &coho.obstruction_basis[0];
        assert_eq!(v.len(), 3);
        assert!(v.iter().all(|x| x.abs() > 1e-9), "el nudo es el ciclo entero");
    }

    #[test]
    fn dos_triangulos_pegados_tienen_dos_ciclos() {
        // Dos triángulos que comparten la arista (1,2): V=4, E=5, C=1.
        // b₁ = E − V + C = 5 − 4 + 1 = 2.
        let edges = [(0, 1), (0, 2), (1, 2), (1, 3), (2, 3)];
        let coho = coho_constante(4, &edges);
        assert_eq!(coho.dim_h0, 1);
        assert_eq!(coho.dim_h1, 2, "dos ciclos independientes");
        assert_eq!(coho.obstruction_basis.len(), 2);
    }

    #[test]
    fn bosque_desconexo_cuenta_componentes() {
        // Dos aristas disjuntas: V=4, E=2, C=2, sin ciclos.
        let coho = coho_constante(4, &[(0, 1), (2, 3)]);
        assert_eq!(coho.dim_h0, 2, "dos componentes conexas");
        assert_eq!(coho.dim_h1, 0);
    }

    #[test]
    fn stalk_de_dimension_mayor_escala_betti() {
        // Mismo triángulo pero con stalk d=2: H⁰ = 1·2, H¹ = 1·2.
        let sheaf = Sheaf::constant(3, &[(0, 1), (1, 2), (2, 0)], 2);
        let coho = compute(&sheaf).unwrap();
        assert_eq!(coho.dim_h0, 2);
        assert_eq!(coho.dim_h1, 2);
    }
}