tawasuyu/hifas

Exploración del reach que Balegas marca no resuelto (constraints over multiple
counters), en módulo aislado, sin tocar la API que Tawasuyu consume.

- reach::CoupledInvariants (politopo C·x ≤ slack, C ≥ 0) + VectorBudget:
  escrow-por-restricción. certificado_politopo (proptest): reparto válido +
  gastos ≤ presupuesto ⟹ todas las restricciones aguantan. La corrección.
- reach::{coboundary, sheaf_laplacian, SheafRouter}: el Laplaciano de haces
  vectorial (L=δᵀδ) reusando nalgebra y la idea del coboundary del MVP.
  diffuse conserva cada recurso y reduce la energía; ker L = totales
  conservados = el H⁰ del MVP.

VEREDICTO HONESTO (REACH.md): la corrección es escrow sobre el politopo, NO
necesita haces. Con restricciones identidad (lo que un presupuesto real tiene),
el Laplaciano de haces colapsa a d_r Laplacianos escalares = el router de T9b
d_r veces (test identidad_colapsa_...). El haz sólo importaría con coupling en
bases incompatibles, que no es requisito real. Consistente con todo el arco: el
haz nunca fue lo que sostuvo la garantía. Recomendación: si Tawasuyu necesita
multi-invariante, usar el certificado del politopo + routers escalares; NO
promover reach.rs a estable ni venderlo como 'potenciado por haces'. Frontera
abierta dicha: restricciones no monótonas (C con signos mixtos) no cubiertas.

85 tests, clippy limpio.

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

diff --git a/REACH.md b/REACH.md
new file mode 100644
index 0000000..4c582ee
--- /dev/null
+++ b/REACH.md
@@ -0,0 +1,78 @@
+# REACH — invariantes acoplados y el Laplaciano de haces vectorial (EXPERIMENTAL)
+
+Estado: **experimental** (`src/reach.rs`). **No toca M4** ni la API que Tawasuyu
+consume. Este archivo es el veredicto honesto de si el haz aporta algo en el reach
+multi-invariante — el sitio que Balegas (SRDS 2015) marca como no resuelto
+("constraints over multiple counters... store them in a single object").
+
+Fuimos "directo al Laplaciano de haces vectorial" a propósito, aceptando el riesgo
+de repetir el patrón "el haz era la trampa". Aquí está lo que salió.
+
+## Qué se construyó (y qué prueba cada test)
+
+| Pieza | Qué es | Test |
+|---|---|---|
+| Certificado del politopo | escrow-por-restricción: `C·Σ b_v ≤ slack`, gasto `s_v ≤ b_v` ⟹ todas las restricciones aguantan (`C ≥ 0`) | `certificado_politopo` |
+| Laplaciano de haces `L = δᵀδ` | stalks vectoriales `ℝ^{d_r}`; ker = totales conservados | `el_nucleo_son_los_totales_conservados` |
+| Difusión vectorial | `x ← x − αLx`: conserva cada recurso, reduce la energía | `la_difusion_conserva_los_totales` |
+| Colapso identidad | con restricciones `I`, `L` = `d_r` Laplacianos escalares | `identidad_colapsa_a_d_difusiones_escalares` |
+
+## El veredicto (sin adornos)
+
+**1. La corrección NO necesita haces.** El multi-invariante coordination-free se
+resuelve con **escrow-por-restricción sobre el politopo**: presupuesto vectorial
+por dispositivo, un gasto es local-seguro sii cabe en el presupuesto de *cada*
+restricción que toca. El certificado (`C·consumed ≤ slack`) sale de que `C ≥ 0` y
+`Σ gastos ≤ Σ presupuestos`. Es la extensión natural de T9–T11 a varios
+contadores, y es **exactamente** el "razonar sobre todos juntos / un solo objeto"
+que Balegas dice que hace falta. Real, usable, sobrio — y **prior-art conceptual**,
+no novedad.
+
+**2. El haz, con las restricciones que un presupuesto de verdad tiene (identidad),
+NO aporta corrección ni nada nuevo.** El Laplaciano de haces del haz constante es
+`Laplaciano_de_grafo ⊗ I`: difundir el presupuesto vectorial es, literalmente,
+correr `d_r` veces el **router escalar de T9b** (`identidad_colapsa_...` lo prueba
+entrada por entrada). El haz sólo se separaría del caso escalar con **restricciones
+no triviales `R_e ≠ I`** — dispositivos que miden recursos acoplados en **bases
+incompatibles**. Eso no es un requisito de ningún sistema de presupuestos real que
+conozcamos; es una solución buscándose un problema.
+
+**3. Lo único genuinamente bonito es una conexión, no una capacidad.** `ker L` =
+secciones globales = los `d_r` totales conservados = el `H⁰` del MVP (§7.2). La
+misma álgebra lineal que en el motor lineal era un **detector con falso-positivo**
+aquí lee las **cantidades que la difusión conserva**. Es un cierre de arco elegante
+y honesto — pero es interpretación, no un teorema nuevo.
+
+## Consistente con todo el proyecto
+
+El patrón se repite por última vez, y conviene nombrarlo: **cada vez que el haz fue
+la protagonista, no fue lo que sostuvo la garantía.** La corrección vino del
+politopo/escrow (T9–T11, prior-art bien hecho); el haz fue, en el mejor caso, una
+lente o una comodidad de ruteo que colapsa al caso escalar. El reach multi-
+invariante **no** cambia eso.
+
+## Recomendación (para no ser contraproducente)
+
+- **Si Tawasuyu necesita multi-invariante:** dale **escrow-por-restricción** (el
+  certificado del politopo) y rutea con `d_r` routers escalares de T9b. **No** lo
+  vendas ni lo cablees como "multi-invariante potenciado por haces": no lo es.
+- **No promuevas `reach.rs` a la API estable** hasta que aparezca un caso real con
+  coupling-en-bases-incompatibles que justifique el Laplaciano de haces sobre `d_r`
+  difusiones escalares. Hasta entonces es más máquina de la que el problema pide.
+- **Lo que sí es sólido y adoptable:** el `CoupledInvariants` + certificado del
+  politopo. Eso extiende M4 a varios contadores de forma correcta y probada.
+
+## Frontera abierta (donde NO llegamos, dicho)
+
+- **Restricciones no monótonas (`C` con signos mixtos):** el certificado de escrow
+  **no aplica** (consumir un recurso puede *liberar* slack de otra restricción). Es
+  otro problema, probablemente necesita reservas de otro tipo; no está resuelto aquí.
+- **Reparto óptimo del presupuesto vectorial:** el certificado aguanta cualquier
+  reparto válido; cuál minimiza rebalanceos futuros es optimización sobre el
+  politopo — ahí la difusión ayuda, pero es rendimiento, no corrección (R-T9a
+  elevado a vectores).
+
+**Conclusión de una línea:** el reach multi-invariante es real y lo resolvimos
+—con escrow sobre el politopo, no con haces—; el Laplaciano de haces vectorial es
+correcto pero, con restricciones identidad, no es más que el router escalar `d_r`
+veces. Registrado como tal.
diff --git a/src/lib.rs b/src/lib.rs
index 0039688..ba7262b 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -25,6 +25,7 @@ pub mod linalg;
 pub mod native;
 pub mod nerve;
 pub mod oracle;
+pub mod reach;
 pub mod router;
 pub mod sheaf;
 pub mod tarski;
diff --git a/src/reach.rs b/src/reach.rs
new file mode 100644
index 0000000..22f1672
--- /dev/null
+++ b/src/reach.rs
@@ -0,0 +1,255 @@
+//! `reach` — invariantes acoplados y ruteo por el **Laplaciano de haces vectorial**
+//! (EXPERIMENTAL, no toca M4).
+//!
+//! El reach que Balegas señala como no resuelto: restricciones sobre **varios
+//! contadores a la vez** ("constraints over multiple counters"). La región válida
+//! `V` es un **politopo** (intersección de semiespacios `C·x ≤ slack`).
+//!
+//! **Separación honesta de qué hace qué (la lección de T5–T11):**
+//! - **La corrección NO necesita haces.** Es **escrow-por-restricción** sobre el
+//!   politopo: presupuesto vectorial por dispositivo `b_v ∈ ℝ^{d_r}`, con
+//!   `C·Σ b_v ≤ slack`. Un gasto `a ≤ b_v` mantiene `C·Σ gastos ≤ slack` para todo
+//!   `C ≥ 0` (restricciones monótonas). Ese es el certificado.
+//! - **El haz SÍ es la herramienta correcta para RUTEAR** el presupuesto vectorial
+//!   sobre el grafo: el Laplaciano de haces (Hansen–Ghrist) generaliza el
+//!   Laplaciano escalar de `router.rs` a stalks vectoriales. Es **eficiencia de
+//!   ruteo, no corrección** — y con restricciones identidad colapsa a `d_r`
+//!   difusiones escalares independientes (ver `REACH.md` para el veredicto).
+//!
+//! Cierre de arco: se reusa `nalgebra` y la idea del coboundary `δ` del MVP. El
+//! Laplaciano `L = δᵀδ` que en el motor lineal era un **detector con falso-positivo**
+//! reaparece aquí como **ruteador legítimo** de derechos acoplados.
+
+use nalgebra::{DMatrix, DVector};
+
+/// Restricciones lineales acopladas `C·x ≤ slack`. `C` con entradas ≥ 0
+/// (restricciones monótonas: consumir más recurso consume más slack).
+#[derive(Debug, Clone)]
+pub struct CoupledInvariants {
+    /// Matriz de restricciones `d_c × d_r` (constraints × resources), `C ≥ 0`.
+    pub c: DMatrix<f64>,
+    /// Slack por restricción (`d_c`).
+    pub slack: DVector<f64>,
+}
+
+impl CoupledInvariants {
+    /// Número de recursos (dimensión del presupuesto por dispositivo).
+    pub fn resources(&self) -> usize {
+        self.c.ncols()
+    }
+
+    /// ¿El consumo total `consumed ∈ ℝ^{d_r}` respeta todas las restricciones?
+    pub fn holds(&self, consumed: &DVector<f64>) -> bool {
+        let used = &self.c * consumed;
+        (0..used.len()).all(|i| used[i] <= self.slack[i] + 1e-9)
+    }
+}
+
+/// Presupuesto vectorial repartido: `per_device[v] ∈ ℝ^{d_r}`.
+#[derive(Debug, Clone)]
+pub struct VectorBudget {
+    /// Un vector de recurso por dispositivo.
+    pub per_device: Vec<DVector<f64>>,
+}
+
+impl VectorBudget {
+    /// Suma de los presupuestos (el recurso total repartido).
+    pub fn total(&self) -> DVector<f64> {
+        let d = self.per_device[0].len();
+        self.per_device.iter().fold(DVector::zeros(d), |acc, b| acc + b)
+    }
+
+    /// El certificado del reparto: `b_v ≥ 0` y `C·Σ b_v ≤ slack`.
+    pub fn is_valid_split(&self, inv: &CoupledInvariants) -> bool {
+        let no_neg = self.per_device.iter().all(|b| b.iter().all(|&x| x >= -1e-9));
+        no_neg && inv.holds(&self.total())
+    }
+}
+
+/// Ensambla el **coboundary de haces** `δ: C⁰ → C¹` con stalks vectoriales de
+/// dimensión `d` y restricciones **identidad** (haz constante). Para arista
+/// `e=(u,v)`, bloque `+I` en `u`, `−I` en `v`. (La generalización a restricciones
+/// no triviales `R_e` sustituye `I` por `R_e`; ver `REACH.md`.)
+pub fn coboundary(num_vertices: usize, edges: &[(usize, usize)], d: usize) -> DMatrix<f64> {
+    let mut delta = DMatrix::zeros(edges.len() * d, num_vertices * d);
+    for (e, &(u, v)) in edges.iter().enumerate() {
+        for k in 0..d {
+            delta[(e * d + k, u * d + k)] += 1.0;
+            delta[(e * d + k, v * d + k)] -= 1.0;
+        }
+    }
+    delta
+}
+
+/// El **Laplaciano de haces** `L = δᵀδ` (bloque `(N·d) × (N·d)`). Para el haz
+/// constante es el Laplaciano de grafo `⊗ Iₐ`; su núcleo son las secciones
+/// globales (los `d` "totales conservados", uno por recurso y componente conexa).
+pub fn sheaf_laplacian(num_vertices: usize, edges: &[(usize, usize)], d: usize) -> DMatrix<f64> {
+    let delta = coboundary(num_vertices, edges, d);
+    delta.transpose() * delta
+}
+
+/// Rutea un presupuesto vectorial sobre el grafo por difusión del Laplaciano de
+/// haces. Un paso `x ← x − α·L·x` con `α ≤ 1/(1+grado_máx)`: **conserva** cada
+/// total de recurso (`1` por componente está en `ker L`) y **suaviza** (reduce la
+/// energía de desacuerdo `xᵀLx`). Escala la reciprocidad del router escalar (T9b)
+/// a stalks vectoriales.
+pub struct SheafRouter {
+    laplacian: DMatrix<f64>,
+    alpha: f64,
+    d: usize,
+    n: usize,
+}
+
+impl SheafRouter {
+    /// Construye el ruteador para `num_vertices` dispositivos con presupuesto de
+    /// dimensión `d` sobre el grafo `edges`.
+    pub fn new(num_vertices: usize, edges: &[(usize, usize)], d: usize) -> Self {
+        let laplacian = sheaf_laplacian(num_vertices, edges, d);
+        let max_diag = (0..laplacian.nrows())
+            .map(|i| laplacian[(i, i)])
+            .fold(0.0_f64, f64::max);
+        Self {
+            laplacian,
+            alpha: 1.0 / (2.0 * (1.0 + max_diag)),
+            d,
+            n: num_vertices,
+        }
+    }
+
+    /// Apila los presupuestos por dispositivo en un vector `(N·d)`.
+    fn stack(&self, budget: &VectorBudget) -> DVector<f64> {
+        let mut x = DVector::zeros(self.n * self.d);
+        for (v, b) in budget.per_device.iter().enumerate() {
+            for k in 0..self.d {
+                x[v * self.d + k] = b[k];
+            }
+        }
+        x
+    }
+
+    /// La energía de desacuerdo `xᵀLx` (0 sii todos los dispositivos coinciden).
+    pub fn energy(&self, budget: &VectorBudget) -> f64 {
+        let x = self.stack(budget);
+        (x.transpose() * &self.laplacian * &x)[(0, 0)]
+    }
+
+    /// Un paso de difusión, in situ. Conserva cada total de recurso.
+    pub fn diffuse(&self, budget: &mut VectorBudget) {
+        let x = self.stack(budget);
+        let next = &x - self.alpha * (&self.laplacian * &x);
+        for v in 0..self.n {
+            for k in 0..self.d {
+                budget.per_device[v][k] = next[v * self.d + k];
+            }
+        }
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use proptest::prelude::*;
+
+    fn vb(rows: &[&[f64]]) -> VectorBudget {
+        VectorBudget {
+            per_device: rows.iter().map(|r| DVector::from_row_slice(r)).collect(),
+        }
+    }
+
+    /// El núcleo del Laplaciano de haces = las secciones globales: para un grafo
+    /// conexo, `dim ker L = d` (los `d` totales conservados). Es el mismo `H⁰` del
+    /// MVP (§7.2), ahora leído como "cantidades que la difusión conserva".
+    #[test]
+    fn el_nucleo_son_los_totales_conservados() {
+        // Triángulo conexo, presupuesto de dimensión 2.
+        let l = sheaf_laplacian(3, &[(0, 1), (1, 2), (2, 0)], 2);
+        // El vector constante por componente (recurso 0 = 1 en todos) está en ker.
+        let mut ones0 = DVector::zeros(6);
+        for v in 0..3 {
+            ones0[v * 2] = 1.0; // recurso 0 uniforme
+        }
+        let lx = &l * &ones0;
+        assert!(lx.iter().all(|&y| y.abs() < 1e-9), "el total del recurso 0 se conserva");
+        // Rango: N·d − dim(ker) = 6 − 2 = 4 (conexo ⟹ ker de dimensión d=2).
+        assert_eq!(l.rank(1e-9), 4);
+    }
+
+    /// La difusión conserva cada total de recurso y no aumenta la energía.
+    #[test]
+    fn la_difusion_conserva_los_totales() {
+        let router = SheafRouter::new(4, &[(0, 1), (1, 2), (2, 3), (3, 0)], 2);
+        let mut budget = vb(&[&[10.0, 0.0], &[0.0, 6.0], &[0.0, 0.0], &[2.0, 4.0]]);
+        let antes = budget.total();
+        let energia_antes = router.energy(&budget);
+
+        for _ in 0..50 {
+            router.diffuse(&mut budget);
+        }
+        let despues = budget.total();
+        for k in 0..2 {
+            assert!((antes[k] - despues[k]).abs() < 1e-6, "recurso {k} conservado");
+        }
+        assert!(router.energy(&budget) <= energia_antes + 1e-9, "energía no aumenta");
+        // Y suaviza: tras difundir, los dispositivos casi coinciden (energía ~0).
+        assert!(router.energy(&budget) < 1e-3, "la difusión reconcilia el reparto");
+    }
+
+    /// Con restricciones identidad, el Laplaciano de haces son `d` Laplacianos
+    /// escalares independientes: difundir el vector = difundir cada recurso por
+    /// separado. (El veredicto honesto de `REACH.md`: el haz no aporta aquí más
+    /// que el router escalar; sólo importaría con restricciones no triviales.)
+    #[test]
+    fn identidad_colapsa_a_d_difusiones_escalares() {
+        let l2 = sheaf_laplacian(3, &[(0, 1), (1, 2)], 2);
+        let l1 = sheaf_laplacian(3, &[(0, 1), (1, 2)], 1);
+        // El bloque del recurso 0 de L2 coincide con el Laplaciano escalar L1.
+        for i in 0..3 {
+            for j in 0..3 {
+                assert_eq!(l2[(i * 2, j * 2)], l1[(i, j)]);
+            }
+        }
+    }
+
+    proptest! {
+        /// **EL CERTIFICADO del politopo** (la corrección, sin haces): para todo
+        /// reparto válido y todo gasto `s_v ≤ b_v`, el consumo total respeta TODAS
+        /// las restricciones acopladas `C·consumed ≤ slack`. Vale para `C ≥ 0`.
+        #[test]
+        fn certificado_politopo(
+            c_raw in prop::collection::vec(0.0f64..5.0, 4), // 2×2 constraint matrix ≥ 0
+            slack0 in 10.0f64..100.0,
+            slack1 in 10.0f64..100.0,
+            b in prop::collection::vec(0.0f64..30.0, 6),     // 3 devices × 2 resources
+            frac in prop::collection::vec(0.0f64..1.0, 6),   // gasto = frac · budget
+        ) {
+            let c = DMatrix::from_row_slice(2, 2, &c_raw);
+            let slack = DVector::from_row_slice(&[slack0, slack1]);
+            let inv = CoupledInvariants { c: c.clone(), slack };
+
+            // Reparto proporcional válido: escalar b para que C·Σb ≤ slack.
+            let raw = VectorBudget {
+                per_device: (0..3).map(|v| DVector::from_row_slice(&b[v*2..v*2+2])).collect(),
+            };
+            let used = &c * raw.total();
+            let escala = [0usize, 1].iter()
+                .map(|&i| if used[i] > 1e-9 { inv.slack[i] / used[i] } else { 1.0 })
+                .fold(1.0_f64, f64::min)
+                .min(1.0);
+            let budget = VectorBudget {
+                per_device: raw.per_device.iter().map(|bv| bv * escala).collect(),
+            };
+            prop_assert!(budget.is_valid_split(&inv), "el reparto escalado es válido");
+
+            // Gastos dentro del presupuesto de cada dispositivo.
+            let d = inv.resources();
+            let mut consumed = DVector::zeros(d);
+            for (v, bv) in budget.per_device.iter().enumerate() {
+                for k in 0..d {
+                    consumed[k] += bv[k] * frac[v * d + k]; // s_v[k] = frac · b_v[k] ≤ b_v[k]
+                }
+            }
+            prop_assert!(inv.holds(&consumed), "el consumo respeta todas las restricciones");
+        }
+    }
+}