T1: cohomología f64 + prueba de cordura del haz constante

Implementa el corazón (§7):
- linalg: F64Backend::rank y cokernel_basis vía RREF con pivoteo parcial;
  Matrix con identity/transpose/get/set.
- sheaf: ensamblado del coboundary δ (bloque +R_u / -R_v por arista) con
  validación de dimensiones, y constructor Sheaf::constant.
- cohomology::compute + compute_with(backend): rank δ, dim H⁰, dim H¹ y
  base del cokernel como localización de la obstrucción.

Cordura topológica (§7.2) verde: árbol H¹=0, triángulo H¹ dim1, dos
triángulos pegados H¹ dim2, bosque desconexo H⁰=2, stalk d=2 escala Betti.
9 tests pasan, sin warnings de clippy.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

src/cohomology.rs

@@ -1,10 +1,9 @@
 //! `cohomology` — la única cuenta que importa (§6.4, §7).
 //!
 //! `rank δ`, `dim H⁰`, `dim H¹` y una base del cokernel para localizar el nudo.
-//! M0: forma del resultado y firma de `compute`. La cuenta llega en T1.
 
-use crate::error::{Result, SheafError};
-use crate::linalg::Scalar;
+use crate::error::Result;
+use crate::linalg::{F64Backend, LinAlg, Scalar};
 use crate::sheaf::Sheaf;
 
 /// Resultado del cálculo cohomológico sobre un haz.
@@ -25,9 +24,101 @@ pub struct Cohomology {
     pub obstruction_basis: Vec<Vec<Scalar>>,
 }
 
-/// Calcula la cohomología del haz vía eliminación gaussiana (§7).
+/// 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.
 ///
-/// M0: no implementado; llega en T1.
-pub fn compute(_sheaf: &Sheaf) -> Result<Cohomology> {
-    Err(SheafError::NotImplemented("cohomology::compute llega en T1"))
+/// 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::*;
+
+    /// 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₁.
+    fn coho_constante(nv: usize, edges: &[(usize, usize)]) -> Cohomology {
+        let sheaf = Sheaf::constant(nv, edges, 1);
+        compute(&sheaf).expect("haz constante bien formado")
+    }
+
+    #[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);
+    }
 }

src/linalg.rs

@@ -1,17 +1,13 @@
 //! Aislamiento del álgebra lineal tras un `trait` (§5.1, §7.1).
 //!
 //! Dos aritméticas, misma interfaz:
-//! - `F64Backend`: rápido, con tolerancia epsilon. Para arrancar el MVP.
+//! - `F64Backend`: rápido, con tolerancia epsilon. Para arrancar el MVP (T1).
 //! - `Gf2Backend` (M2): exacto, bit-packed sobre 𝔽₂. Aún no implementado.
-//!
-//! En M0 solo dejamos la *forma*: la impl `f64` existe pero sus operaciones
-//! todavía no calculan nada. La lógica llega en T1.
 
 /// Escalar del backend activo. Alias que resuelve la impl elegida (§6.1, §7).
 pub type Scalar = f64;
 
-/// Matriz densa mínima almacenada por filas. Placeholder de M0; en T1 puede
-/// migrarse a `nalgebra::DMatrix` o mantenerse propia para el backend `GF(2)`.
+/// Matriz densa almacenada por filas (row-major).
 #[derive(Debug, Clone, PartialEq)]
 pub struct Matrix {
     pub rows: usize,
@@ -29,6 +25,54 @@ impl Matrix {
             data: vec![0.0; rows * cols],
         }
     }
+
+    /// Matriz identidad `n × n`.
+    pub fn identity(n: usize) -> Self {
+        let mut m = Self::zeros(n, n);
+        for i in 0..n {
+            m.set(i, i, 1.0);
+        }
+        m
+    }
+
+    /// Lee la entrada `(r, c)`.
+    #[inline]
+    pub fn get(&self, r: usize, c: usize) -> Scalar {
+        self.data[r * self.cols + c]
+    }
+
+    /// Escribe la entrada `(r, c)`.
+    #[inline]
+    pub fn set(&mut self, r: usize, c: usize, v: Scalar) {
+        self.data[r * self.cols + c] = v;
+    }
+
+    /// Suma `v` a la entrada `(r, c)` (para acumular bloques solapados).
+    #[inline]
+    pub fn add_to(&mut self, r: usize, c: usize, v: Scalar) {
+        self.data[r * self.cols + c] += v;
+    }
+
+    /// Intercambia dos filas in situ.
+    fn swap_rows(&mut self, a: usize, b: usize) {
+        if a == b {
+            return;
+        }
+        for c in 0..self.cols {
+            self.data.swap(a * self.cols + c, b * self.cols + c);
+        }
+    }
+
+    /// Transpuesta.
+    pub fn transpose(&self) -> Matrix {
+        let mut t = Matrix::zeros(self.cols, self.rows);
+        for r in 0..self.rows {
+            for c in 0..self.cols {
+                t.set(c, r, self.get(r, c));
+            }
+        }
+        t
+    }
 }
 
 /// Interfaz de álgebra lineal que el motor usa sin conocer la aritmética concreta.
@@ -41,33 +85,161 @@ pub trait LinAlg {
 
     /// Base del cokernel `C¹ / im δ`: vectores de longitud `m.rows` que no
     /// provienen de ninguna sección global. Vacío si `H¹ = 0`.
+    ///
+    /// Se calcula como el espacio nulo izquierdo de `m` (= `ker(mᵀ)`), que en un
+    /// espacio con producto interno representa `coker(m) ≅ (im m)^⊥`. Cada vector
+    /// vive en el espacio de aristas `C¹`; sus componentes no nulas son el nudo.
     fn cokernel_basis(&self, m: &Matrix) -> Vec<Vec<Scalar>>;
 }
 
 /// Backend `f64` del MVP (§7.1). Rango con tolerancia epsilon.
-///
-/// M0: shell vacío. La eliminación gaussiana llega en T1.
-#[derive(Debug, Default, Clone, Copy)]
+#[derive(Debug, Clone, Copy)]
 pub struct F64Backend {
     /// Tolerancia para considerar un pivote como cero.
     pub epsilon: Scalar,
 }
 
+impl Default for F64Backend {
+    fn default() -> Self {
+        Self::new()
+    }
+}
+
 impl F64Backend {
     /// Backend con la tolerancia por defecto.
     pub fn new() -> Self {
         Self { epsilon: 1e-9 }
     }
+
+    /// Forma escalonada reducida por filas (RREF) con pivoteo parcial.
+    /// Devuelve la matriz reducida y la lista de columnas pivote (en orden).
+    fn rref(&self, m: &Matrix) -> (Matrix, Vec<usize>) {
+        let mut a = m.clone();
+        let mut pivot_cols = Vec::new();
+        let mut row = 0;
+
+        for col in 0..a.cols {
+            if row >= a.rows {
+                break;
+            }
+            // Pivoteo parcial: mayor magnitud en la columna, filas >= row.
+            let mut sel = row;
+            let mut best = a.get(row, col).abs();
+            for r in (row + 1)..a.rows {
+                let v = a.get(r, col).abs();
+                if v > best {
+                    best = v;
+                    sel = r;
+                }
+            }
+            if best <= self.epsilon {
+                continue; // columna sin pivote → variable libre
+            }
+            a.swap_rows(row, sel);
+
+            // Normaliza la fila pivote.
+            let piv = a.get(row, col);
+            for c in 0..a.cols {
+                let v = a.get(row, c) / piv;
+                a.set(row, c, v);
+            }
+
+            // Elimina la columna en el resto de filas.
+            for r in 0..a.rows {
+                if r == row {
+                    continue;
+                }
+                let factor = a.get(r, col);
+                if factor.abs() > 0.0 {
+                    for c in 0..a.cols {
+                        let v = a.get(r, c) - factor * a.get(row, c);
+                        a.set(r, c, v);
+                    }
+                }
+            }
+
+            pivot_cols.push(col);
+            row += 1;
+        }
+
+        (a, pivot_cols)
+    }
+
+    /// Base del espacio nulo `{x : m·x = 0}`, una por columna libre.
+    fn null_space(&self, m: &Matrix) -> Vec<Vec<Scalar>> {
+        let (rref, pivots) = self.rref(m);
+        let ncols = m.cols;
+        let is_pivot = |c: usize| pivots.contains(&c);
+
+        let mut basis = Vec::new();
+        for free in 0..ncols {
+            if is_pivot(free) {
+                continue;
+            }
+            // Variable libre = 1; las pivote se despejan desde la RREF.
+            let mut vec = vec![0.0; ncols];
+            vec[free] = 1.0;
+            for (i, &pc) in pivots.iter().enumerate() {
+                vec[pc] = -rref.get(i, free);
+            }
+            basis.push(vec);
+        }
+        basis
+    }
 }
 
 impl LinAlg for F64Backend {
-    fn rank(&self, _m: &Matrix) -> usize {
-        // TODO(T1): eliminación gaussiana con pivoteo, contar pivotes > epsilon.
-        unimplemented!("F64Backend::rank llega en T1")
+    fn rank(&self, m: &Matrix) -> usize {
+        self.rref(m).1.len()
     }
 
-    fn cokernel_basis(&self, _m: &Matrix) -> Vec<Vec<Scalar>> {
-        // TODO(T1): extraer base del cokernel a partir de la forma escalonada.
-        unimplemented!("F64Backend::cokernel_basis llega en T1")
+    fn cokernel_basis(&self, m: &Matrix) -> Vec<Vec<Scalar>> {
+        // coker(m) ≅ ker(mᵀ): vectores del espacio de aristas ortogonales a im(m).
+        self.null_space(&m.transpose())
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    /// Matriz a partir de filas, para tests legibles.
+    fn mat(rows: &[&[Scalar]]) -> Matrix {
+        let r = rows.len();
+        let c = rows[0].len();
+        let mut m = Matrix::zeros(r, c);
+        for (i, row) in rows.iter().enumerate() {
+            for (j, &v) in row.iter().enumerate() {
+                m.set(i, j, v);
+            }
+        }
+        m
+    }
+
+    #[test]
+    fn rank_de_matriz_conocida() {
+        // Rango 2: la tercera fila es suma de las dos primeras.
+        let m = mat(&[&[1.0, 0.0, 1.0], &[0.0, 1.0, 1.0], &[1.0, 1.0, 2.0]]);
+        assert_eq!(F64Backend::new().rank(&m), 2);
+    }
+
+    #[test]
+    fn rank_identidad_y_cero() {
+        let be = F64Backend::new();
+        assert_eq!(be.rank(&Matrix::identity(4)), 4);
+        assert_eq!(be.rank(&Matrix::zeros(3, 5)), 0);
+    }
+
+    #[test]
+    fn null_space_dimension_correcta() {
+        // 1×3 de rango 1 ⇒ espacio nulo de dimensión 2.
+        let m = mat(&[&[1.0, 2.0, 3.0]]);
+        let ns = F64Backend::new().null_space(&m);
+        assert_eq!(ns.len(), 2);
+        // Cada vector debe anular la fila: 1·x0 + 2·x1 + 3·x2 = 0.
+        for v in &ns {
+            let dot = v[0] + 2.0 * v[1] + 3.0 * v[2];
+            assert!(dot.abs() < 1e-9, "vector no está en el núcleo: {v:?}");
+        }
     }
 }

src/sheaf.rs

@@ -1,7 +1,4 @@
 //! `sheaf` — el haz celular y el ensamblado del coboundary `δ` (§6.3, §7).
-//!
-//! M0: forma de los tipos y firma de `coboundary`. El ensamblado real
-//! (bloque por arista, `+R_{u⊴e}` / `−R_{v⊴e}`) llega en T1.
 
 use crate::error::{Result, SheafError};
 use crate::linalg::Matrix;
@@ -11,7 +8,7 @@ use crate::linalg::Matrix;
 pub type RestrictionMap = Matrix;
 
 /// Un haz celular sobre el nervio.
-#[derive(Debug, Default)]
+#[derive(Debug, Default, Clone)]
 pub struct Sheaf {
     /// Para cada réplica: la dimensión de su stalk `F(v)`.
     pub vertex_stalks: Vec<usize>,
@@ -34,13 +31,106 @@ impl Sheaf {
         self.edge_stalks.iter().sum()
     }
 
+    /// Desplazamiento de la columna donde empieza el bloque del vértice `v` en `C⁰`.
+    fn vertex_offset(&self, v: usize) -> usize {
+        self.vertex_stalks[..v].iter().sum()
+    }
+
+    /// Desplazamiento de la fila donde empieza el bloque de la arista `e` en `C¹`.
+    fn edge_offset(&self, e: usize) -> usize {
+        self.edge_stalks[..e].iter().sum()
+    }
+
+    /// Comprueba que los tamaños de las matrices de restricción cuadran con las
+    /// dimensiones declaradas de stalks (fila = arista, columna = vértice).
+    fn validate(&self) -> Result<()> {
+        let ne = self.edge_stalks.len();
+        if self.restrictions.len() != ne || self.edge_endpoints.len() != ne {
+            return Err(SheafError::DimensionMismatch(format!(
+                "aristas: {ne} stalks, {} restricciones, {} extremos",
+                self.restrictions.len(),
+                self.edge_endpoints.len()
+            )));
+        }
+        let nv = self.vertex_stalks.len();
+        for (e, ((ru, rv), &(u, v))) in self
+            .restrictions
+            .iter()
+            .zip(self.edge_endpoints.iter())
+            .enumerate()
+        {
+            if u >= nv || v >= nv {
+                return Err(SheafError::MissingEntity(format!(
+                    "arista {e} referencia vértice inexistente ({u},{v}) de {nv}"
+                )));
+            }
+            let ed = self.edge_stalks[e];
+            if ru.rows != ed || ru.cols != self.vertex_stalks[u] {
+                return Err(SheafError::DimensionMismatch(format!(
+                    "arista {e}: R_u es {}×{}, se esperaba {ed}×{}",
+                    ru.rows, ru.cols, self.vertex_stalks[u]
+                )));
+            }
+            if rv.rows != ed || rv.cols != self.vertex_stalks[v] {
+                return Err(SheafError::DimensionMismatch(format!(
+                    "arista {e}: R_v es {}×{}, se esperaba {ed}×{}",
+                    rv.rows, rv.cols, self.vertex_stalks[v]
+                )));
+            }
+        }
+        Ok(())
+    }
+
     /// Ensambla la matriz global de coboundary `δ: C⁰ → C¹`.
     ///
     /// Convención (§6.3): para arista `e=(u,v)`,
     /// `(δx)_e = R_{u⊴e}·x_u − R_{v⊴e}·x_v`.
     ///
-    /// M0: no implementado; llega en T1.
+    /// Bloque de la arista `e`: `+R_u` en las columnas de `u`, `−R_v` en las de
+    /// `v`. Se acumula (`+=`) para tolerar lazos o multiaristas.
     pub fn coboundary(&self) -> Result<Matrix> {
-        Err(SheafError::NotImplemented("Sheaf::coboundary llega en T1"))
+        self.validate()?;
+        let mut delta = Matrix::zeros(self.dim_c1(), self.dim_c0());
+
+        for (e, ((ru, rv), &(u, v))) in self
+            .restrictions
+            .iter()
+            .zip(self.edge_endpoints.iter())
+            .enumerate()
+        {
+            let row0 = self.edge_offset(e);
+            let ed = self.edge_stalks[e];
+
+            let ucol0 = self.vertex_offset(u);
+            for i in 0..ed {
+                for j in 0..ru.cols {
+                    delta.add_to(row0 + i, ucol0 + j, ru.get(i, j));
+                }
+            }
+
+            let vcol0 = self.vertex_offset(v);
+            for i in 0..ed {
+                for j in 0..rv.cols {
+                    delta.add_to(row0 + i, vcol0 + j, -rv.get(i, j));
+                }
+            }
+        }
+
+        Ok(delta)
+    }
+
+    /// Construye el **haz constante** de dimensión de stalk `d` sobre un grafo:
+    /// todas las stalks (vértice y arista) valen `d`, y toda restricción es la
+    /// identidad. Es el andamio de la prueba de cordura topológica (§7.2).
+    ///
+    /// `num_vertices` réplicas y `edges` como pares `(u, v)` de índices.
+    pub fn constant(num_vertices: usize, edges: &[(usize, usize)], d: usize) -> Self {
+        let id = Matrix::identity(d);
+        Self {
+            vertex_stalks: vec![d; num_vertices],
+            edge_stalks: vec![d; edges.len()],
+            restrictions: edges.iter().map(|_| (id.clone(), id.clone())).collect(),
+            edge_endpoints: edges.to_vec(),
+        }
     }
 }