tawasuyu/hifas

1//! Aislamiento del álgebra lineal tras un `trait` (§5.1, §7.1).
2//!
3//! Dos aritméticas, misma interfaz:
4//! - `F64Backend`: rápido, con tolerancia epsilon. Para arrancar el MVP (T1).
5//! - `Gf2Backend` (M2): exacto, bit-packed sobre 𝔽₂. Aún no implementado.
6
7/// Escalar del backend activo. Alias que resuelve la impl elegida (§6.1, §7).
8pub type Scalar = f64;
9
10/// Matriz densa almacenada por filas (row-major).
11#[derive(Debug, Clone, PartialEq)]
12pub struct Matrix {
13    pub rows: usize,
14    pub cols: usize,
15    /// Datos row-major: `data[r * cols + c]`.
16    pub data: Vec<Scalar>,
17}
18
19impl Matrix {
20    /// Matriz de ceros de tamaño `rows × cols`.
21    pub fn zeros(rows: usize, cols: usize) -> Self {
22        Self {
23            rows,
24            cols,
25            data: vec![0.0; rows * cols],
26        }
27    }
28
29    /// Matriz identidad `n × n`.
30    pub fn identity(n: usize) -> Self {
31        let mut m = Self::zeros(n, n);
32        for i in 0..n {
33            m.set(i, i, 1.0);
34        }
35        m
36    }
37
38    /// Lee la entrada `(r, c)`.
39    #[inline]
40    pub fn get(&self, r: usize, c: usize) -> Scalar {
41        self.data[r * self.cols + c]
42    }
43
44    /// Escribe la entrada `(r, c)`.
45    #[inline]
46    pub fn set(&mut self, r: usize, c: usize, v: Scalar) {
47        self.data[r * self.cols + c] = v;
48    }
49
50    /// Suma `v` a la entrada `(r, c)` (para acumular bloques solapados).
51    #[inline]
52    pub fn add_to(&mut self, r: usize, c: usize, v: Scalar) {
53        self.data[r * self.cols + c] += v;
54    }
55
56    /// Intercambia dos filas in situ.
57    fn swap_rows(&mut self, a: usize, b: usize) {
58        if a == b {
59            return;
60        }
61        for c in 0..self.cols {
62            self.data.swap(a * self.cols + c, b * self.cols + c);
63        }
64    }
65
66    /// Transpuesta.
67    pub fn transpose(&self) -> Matrix {
68        let mut t = Matrix::zeros(self.cols, self.rows);
69        for r in 0..self.rows {
70            for c in 0..self.cols {
71                t.set(c, r, self.get(r, c));
72            }
73        }
74        t
75    }
76}
77
78/// Interfaz de álgebra lineal que el motor usa sin conocer la aritmética concreta.
79///
80/// Las dos operaciones que la cohomología necesita (§7): rango del coboundary y
81/// una base del cokernel para localizar la obstrucción.
82pub trait LinAlg {
83    /// `rank(m)` — número de filas/columnas independientes.
84    fn rank(&self, m: &Matrix) -> usize;
85
86    /// Base del cokernel `C¹ / im δ`: vectores de longitud `m.rows` que no
87    /// provienen de ninguna sección global. Vacío si `H¹ = 0`.
88    ///
89    /// Se calcula como el espacio nulo izquierdo de `m` (= `ker(mᵀ)`), que en un
90    /// espacio con producto interno representa `coker(m) ≅ (im m)^⊥`. Cada vector
91    /// vive en el espacio de aristas `C¹`; sus componentes no nulas son el nudo.
92    fn cokernel_basis(&self, m: &Matrix) -> Vec<Vec<Scalar>>;
93}
94
95/// Backend `f64` del MVP (§7.1). Rango con tolerancia epsilon.
96#[derive(Debug, Clone, Copy)]
97pub struct F64Backend {
98    /// Tolerancia para considerar un pivote como cero.
99    pub epsilon: Scalar,
100}
101
102impl Default for F64Backend {
103    fn default() -> Self {
104        Self::new()
105    }
106}
107
108impl F64Backend {
109    /// Backend con la tolerancia por defecto.
110    pub fn new() -> Self {
111        Self { epsilon: 1e-9 }
112    }
113
114    /// Forma escalonada reducida por filas (RREF) con pivoteo parcial.
115    /// Devuelve la matriz reducida y la lista de columnas pivote (en orden).
116    fn rref(&self, m: &Matrix) -> (Matrix, Vec<usize>) {
117        let mut a = m.clone();
118        let mut pivot_cols = Vec::new();
119        let mut row = 0;
120
121        for col in 0..a.cols {
122            if row >= a.rows {
123                break;
124            }
125            // Pivoteo parcial: mayor magnitud en la columna, filas >= row.
126            let mut sel = row;
127            let mut best = a.get(row, col).abs();
128            for r in (row + 1)..a.rows {
129                let v = a.get(r, col).abs();
130                if v > best {
131                    best = v;
132                    sel = r;
133                }
134            }
135            if best <= self.epsilon {
136                continue; // columna sin pivote → variable libre
137            }
138            a.swap_rows(row, sel);
139
140            // Normaliza la fila pivote.
141            let piv = a.get(row, col);
142            for c in 0..a.cols {
143                let v = a.get(row, c) / piv;
144                a.set(row, c, v);
145            }
146
147            // Elimina la columna en el resto de filas.
148            for r in 0..a.rows {
149                if r == row {
150                    continue;
151                }
152                let factor = a.get(r, col);
153                if factor.abs() > 0.0 {
154                    for c in 0..a.cols {
155                        let v = a.get(r, c) - factor * a.get(row, c);
156                        a.set(r, c, v);
157                    }
158                }
159            }
160
161            pivot_cols.push(col);
162            row += 1;
163        }
164
165        (a, pivot_cols)
166    }
167
168    /// Base del espacio nulo `{x : m·x = 0}`, una por columna libre.
169    fn null_space(&self, m: &Matrix) -> Vec<Vec<Scalar>> {
170        let (rref, pivots) = self.rref(m);
171        let ncols = m.cols;
172        let is_pivot = |c: usize| pivots.contains(&c);
173
174        let mut basis = Vec::new();
175        for free in 0..ncols {
176            if is_pivot(free) {
177                continue;
178            }
179            // Variable libre = 1; las pivote se despejan desde la RREF.
180            let mut vec = vec![0.0; ncols];
181            vec[free] = 1.0;
182            for (i, &pc) in pivots.iter().enumerate() {
183                vec[pc] = -rref.get(i, free);
184            }
185            basis.push(vec);
186        }
187        basis
188    }
189}
190
191impl LinAlg for F64Backend {
192    fn rank(&self, m: &Matrix) -> usize {
193        self.rref(m).1.len()
194    }
195
196    fn cokernel_basis(&self, m: &Matrix) -> Vec<Vec<Scalar>> {
197        // coker(m) ≅ ker(mᵀ): vectores del espacio de aristas ortogonales a im(m).
198        self.null_space(&m.transpose())
199    }
200}
201
202#[cfg(test)]
203mod tests {
204    use super::*;
205
206    /// Matriz a partir de filas, para tests legibles.
207    fn mat(rows: &[&[Scalar]]) -> Matrix {
208        let r = rows.len();
209        let c = rows[0].len();
210        let mut m = Matrix::zeros(r, c);
211        for (i, row) in rows.iter().enumerate() {
212            for (j, &v) in row.iter().enumerate() {
213                m.set(i, j, v);
214            }
215        }
216        m
217    }
218
219    #[test]
220    fn rank_de_matriz_conocida() {
221        // Rango 2: la tercera fila es suma de las dos primeras.
222        let m = mat(&[&[1.0, 0.0, 1.0], &[0.0, 1.0, 1.0], &[1.0, 1.0, 2.0]]);
223        assert_eq!(F64Backend::new().rank(&m), 2);
224    }
225
226    #[test]
227    fn rank_identidad_y_cero() {
228        let be = F64Backend::new();
229        assert_eq!(be.rank(&Matrix::identity(4)), 4);
230        assert_eq!(be.rank(&Matrix::zeros(3, 5)), 0);
231    }
232
233    #[test]
234    fn null_space_dimension_correcta() {
235        // 1×3 de rango 1 ⇒ espacio nulo de dimensión 2.
236        let m = mat(&[&[1.0, 2.0, 3.0]]);
237        let ns = F64Backend::new().null_space(&m);
238        assert_eq!(ns.len(), 2);
239        // Cada vector debe anular la fila: 1·x0 + 2·x1 + 3·x2 = 0.
240        for v in &ns {
241            let dot = v[0] + 2.0 * v[1] + 3.0 * v[2];
242            assert!(dot.abs() < 1e-9, "vector no está en el núcleo: {v:?}");
243        }
244    }
245}