//! `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");
        }
    }
}