tawasuyu/hifas

1//! `reach` — invariantes acoplados y ruteo por el **Laplaciano de haces vectorial**
2//! (EXPERIMENTAL, no toca M4).
3//!
4//! El reach que Balegas señala como no resuelto: restricciones sobre **varios
5//! contadores a la vez** ("constraints over multiple counters"). La región válida
6//! `V` es un **politopo** (intersección de semiespacios `C·x ≤ slack`).
7//!
8//! **Separación honesta de qué hace qué (la lección de T5–T11):**
9//! - **La corrección NO necesita haces.** Es **escrow-por-restricción** sobre el
10//!   politopo: presupuesto vectorial por dispositivo `b_v ∈ ℝ^{d_r}`, con
11//!   `C·Σ b_v ≤ slack`. Un gasto `a ≤ b_v` mantiene `C·Σ gastos ≤ slack` para todo
12//!   `C ≥ 0` (restricciones monótonas). Ese es el certificado.
13//! - **El haz SÍ es la herramienta correcta para RUTEAR** el presupuesto vectorial
14//!   sobre el grafo: el Laplaciano de haces (Hansen–Ghrist) generaliza el
15//!   Laplaciano escalar de `router.rs` a stalks vectoriales. Es **eficiencia de
16//!   ruteo, no corrección** — y con restricciones identidad colapsa a `d_r`
17//!   difusiones escalares independientes (ver `REACH.md` para el veredicto).
18//!
19//! Cierre de arco: se reusa `nalgebra` y la idea del coboundary `δ` del MVP. El
20//! Laplaciano `L = δᵀδ` que en el motor lineal era un **detector con falso-positivo**
21//! reaparece aquí como **ruteador legítimo** de derechos acoplados.
22
23use nalgebra::{DMatrix, DVector};
24
25/// Restricciones lineales acopladas `C·x ≤ slack`. `C` con entradas ≥ 0
26/// (restricciones monótonas: consumir más recurso consume más slack).
27#[derive(Debug, Clone)]
28pub struct CoupledInvariants {
29    /// Matriz de restricciones `d_c × d_r` (constraints × resources), `C ≥ 0`.
30    pub c: DMatrix<f64>,
31    /// Slack por restricción (`d_c`).
32    pub slack: DVector<f64>,
33}
34
35impl CoupledInvariants {
36    /// Número de recursos (dimensión del presupuesto por dispositivo).
37    pub fn resources(&self) -> usize {
38        self.c.ncols()
39    }
40
41    /// ¿El consumo total `consumed ∈ ℝ^{d_r}` respeta todas las restricciones?
42    pub fn holds(&self, consumed: &DVector<f64>) -> bool {
43        let used = &self.c * consumed;
44        (0..used.len()).all(|i| used[i] <= self.slack[i] + 1e-9)
45    }
46}
47
48/// Presupuesto vectorial repartido: `per_device[v] ∈ ℝ^{d_r}`.
49#[derive(Debug, Clone)]
50pub struct VectorBudget {
51    /// Un vector de recurso por dispositivo.
52    pub per_device: Vec<DVector<f64>>,
53}
54
55impl VectorBudget {
56    /// Suma de los presupuestos (el recurso total repartido).
57    pub fn total(&self) -> DVector<f64> {
58        let d = self.per_device[0].len();
59        self.per_device.iter().fold(DVector::zeros(d), |acc, b| acc + b)
60    }
61
62    /// El certificado del reparto: `b_v ≥ 0` y `C·Σ b_v ≤ slack`.
63    pub fn is_valid_split(&self, inv: &CoupledInvariants) -> bool {
64        let no_neg = self.per_device.iter().all(|b| b.iter().all(|&x| x >= -1e-9));
65        no_neg && inv.holds(&self.total())
66    }
67}
68
69/// Ensambla el **coboundary de haces** `δ: C⁰ → C¹` con stalks vectoriales de
70/// dimensión `d` y restricciones **identidad** (haz constante). Para arista
71/// `e=(u,v)`, bloque `+I` en `u`, `−I` en `v`. (La generalización a restricciones
72/// no triviales `R_e` sustituye `I` por `R_e`; ver `REACH.md`.)
73pub fn coboundary(num_vertices: usize, edges: &[(usize, usize)], d: usize) -> DMatrix<f64> {
74    let mut delta = DMatrix::zeros(edges.len() * d, num_vertices * d);
75    for (e, &(u, v)) in edges.iter().enumerate() {
76        for k in 0..d {
77            delta[(e * d + k, u * d + k)] += 1.0;
78            delta[(e * d + k, v * d + k)] -= 1.0;
79        }
80    }
81    delta
82}
83
84/// El **Laplaciano de haces** `L = δᵀδ` (bloque `(N·d) × (N·d)`). Para el haz
85/// constante es el Laplaciano de grafo `⊗ Iₐ`; su núcleo son las secciones
86/// globales (los `d` "totales conservados", uno por recurso y componente conexa).
87pub fn sheaf_laplacian(num_vertices: usize, edges: &[(usize, usize)], d: usize) -> DMatrix<f64> {
88    let delta = coboundary(num_vertices, edges, d);
89    delta.transpose() * delta
90}
91
92/// Rutea un presupuesto vectorial sobre el grafo por difusión del Laplaciano de
93/// haces. Un paso `x ← x − α·L·x` con `α ≤ 1/(1+grado_máx)`: **conserva** cada
94/// total de recurso (`1` por componente está en `ker L`) y **suaviza** (reduce la
95/// energía de desacuerdo `xᵀLx`). Escala la reciprocidad del router escalar (T9b)
96/// a stalks vectoriales.
97pub struct SheafRouter {
98    laplacian: DMatrix<f64>,
99    alpha: f64,
100    d: usize,
101    n: usize,
102}
103
104impl SheafRouter {
105    /// Construye el ruteador para `num_vertices` dispositivos con presupuesto de
106    /// dimensión `d` sobre el grafo `edges`.
107    pub fn new(num_vertices: usize, edges: &[(usize, usize)], d: usize) -> Self {
108        let laplacian = sheaf_laplacian(num_vertices, edges, d);
109        let max_diag = (0..laplacian.nrows())
110            .map(|i| laplacian[(i, i)])
111            .fold(0.0_f64, f64::max);
112        Self {
113            laplacian,
114            alpha: 1.0 / (2.0 * (1.0 + max_diag)),
115            d,
116            n: num_vertices,
117        }
118    }
119
120    /// Apila los presupuestos por dispositivo en un vector `(N·d)`.
121    fn stack(&self, budget: &VectorBudget) -> DVector<f64> {
122        let mut x = DVector::zeros(self.n * self.d);
123        for (v, b) in budget.per_device.iter().enumerate() {
124            for k in 0..self.d {
125                x[v * self.d + k] = b[k];
126            }
127        }
128        x
129    }
130
131    /// La energía de desacuerdo `xᵀLx` (0 sii todos los dispositivos coinciden).
132    pub fn energy(&self, budget: &VectorBudget) -> f64 {
133        let x = self.stack(budget);
134        (x.transpose() * &self.laplacian * &x)[(0, 0)]
135    }
136
137    /// Un paso de difusión, in situ. Conserva cada total de recurso.
138    pub fn diffuse(&self, budget: &mut VectorBudget) {
139        let x = self.stack(budget);
140        let next = &x - self.alpha * (&self.laplacian * &x);
141        for v in 0..self.n {
142            for k in 0..self.d {
143                budget.per_device[v][k] = next[v * self.d + k];
144            }
145        }
146    }
147}
148
149#[cfg(test)]
150mod tests {
151    use super::*;
152    use proptest::prelude::*;
153
154    fn vb(rows: &[&[f64]]) -> VectorBudget {
155        VectorBudget {
156            per_device: rows.iter().map(|r| DVector::from_row_slice(r)).collect(),
157        }
158    }
159
160    /// El núcleo del Laplaciano de haces = las secciones globales: para un grafo
161    /// conexo, `dim ker L = d` (los `d` totales conservados). Es el mismo `H⁰` del
162    /// MVP (§7.2), ahora leído como "cantidades que la difusión conserva".
163    #[test]
164    fn el_nucleo_son_los_totales_conservados() {
165        // Triángulo conexo, presupuesto de dimensión 2.
166        let l = sheaf_laplacian(3, &[(0, 1), (1, 2), (2, 0)], 2);
167        // El vector constante por componente (recurso 0 = 1 en todos) está en ker.
168        let mut ones0 = DVector::zeros(6);
169        for v in 0..3 {
170            ones0[v * 2] = 1.0; // recurso 0 uniforme
171        }
172        let lx = &l * &ones0;
173        assert!(lx.iter().all(|&y| y.abs() < 1e-9), "el total del recurso 0 se conserva");
174        // Rango: N·d − dim(ker) = 6 − 2 = 4 (conexo ⟹ ker de dimensión d=2).
175        assert_eq!(l.rank(1e-9), 4);
176    }
177
178    /// La difusión conserva cada total de recurso y no aumenta la energía.
179    #[test]
180    fn la_difusion_conserva_los_totales() {
181        let router = SheafRouter::new(4, &[(0, 1), (1, 2), (2, 3), (3, 0)], 2);
182        let mut budget = vb(&[&[10.0, 0.0], &[0.0, 6.0], &[0.0, 0.0], &[2.0, 4.0]]);
183        let antes = budget.total();
184        let energia_antes = router.energy(&budget);
185
186        for _ in 0..50 {
187            router.diffuse(&mut budget);
188        }
189        let despues = budget.total();
190        for k in 0..2 {
191            assert!((antes[k] - despues[k]).abs() < 1e-6, "recurso {k} conservado");
192        }
193        assert!(router.energy(&budget) <= energia_antes + 1e-9, "energía no aumenta");
194        // Y suaviza: tras difundir, los dispositivos casi coinciden (energía ~0).
195        assert!(router.energy(&budget) < 1e-3, "la difusión reconcilia el reparto");
196    }
197
198    /// Con restricciones identidad, el Laplaciano de haces son `d` Laplacianos
199    /// escalares independientes: difundir el vector = difundir cada recurso por
200    /// separado. (El veredicto honesto de `REACH.md`: el haz no aporta aquí más
201    /// que el router escalar; sólo importaría con restricciones no triviales.)
202    #[test]
203    fn identidad_colapsa_a_d_difusiones_escalares() {
204        let l2 = sheaf_laplacian(3, &[(0, 1), (1, 2)], 2);
205        let l1 = sheaf_laplacian(3, &[(0, 1), (1, 2)], 1);
206        // El bloque del recurso 0 de L2 coincide con el Laplaciano escalar L1.
207        for i in 0..3 {
208            for j in 0..3 {
209                assert_eq!(l2[(i * 2, j * 2)], l1[(i, j)]);
210            }
211        }
212    }
213
214    proptest! {
215        /// **EL CERTIFICADO del politopo** (la corrección, sin haces): para todo
216        /// reparto válido y todo gasto `s_v ≤ b_v`, el consumo total respeta TODAS
217        /// las restricciones acopladas `C·consumed ≤ slack`. Vale para `C ≥ 0`.
218        #[test]
219        fn certificado_politopo(
220            c_raw in prop::collection::vec(0.0f64..5.0, 4), // 2×2 constraint matrix ≥ 0
221            slack0 in 10.0f64..100.0,
222            slack1 in 10.0f64..100.0,
223            b in prop::collection::vec(0.0f64..30.0, 6),     // 3 devices × 2 resources
224            frac in prop::collection::vec(0.0f64..1.0, 6),   // gasto = frac · budget
225        ) {
226            let c = DMatrix::from_row_slice(2, 2, &c_raw);
227            let slack = DVector::from_row_slice(&[slack0, slack1]);
228            let inv = CoupledInvariants { c: c.clone(), slack };
229
230            // Reparto proporcional válido: escalar b para que C·Σb ≤ slack.
231            let raw = VectorBudget {
232                per_device: (0..3).map(|v| DVector::from_row_slice(&b[v*2..v*2+2])).collect(),
233            };
234            let used = &c * raw.total();
235            let escala = [0usize, 1].iter()
236                .map(|&i| if used[i] > 1e-9 { inv.slack[i] / used[i] } else { 1.0 })
237                .fold(1.0_f64, f64::min)
238                .min(1.0);
239            let budget = VectorBudget {
240                per_device: raw.per_device.iter().map(|bv| bv * escala).collect(),
241            };
242            prop_assert!(budget.is_valid_split(&inv), "el reparto escalado es válido");
243
244            // Gastos dentro del presupuesto de cada dispositivo.
245            let d = inv.resources();
246            let mut consumed = DVector::zeros(d);
247            for (v, bv) in budget.per_device.iter().enumerate() {
248                for k in 0..d {
249                    consumed[k] += bv[k] * frac[v * d + k]; // s_v[k] = frac · b_v[k] ≤ b_v[k]
250                }
251            }
252            prop_assert!(inv.holds(&consumed), "el consumo respeta todas las restricciones");
253        }
254    }
255}