//! 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 (T1).
//! - `Gf2Backend` (M2): exacto, bit-packed sobre 𝔽₂. Aún no implementado.

/// Escalar del backend activo. Alias que resuelve la impl elegida (§6.1, §7).
pub type Scalar = f64;

/// Matriz densa almacenada por filas (row-major).
#[derive(Debug, Clone, PartialEq)]
pub struct Matrix {
    pub rows: usize,
    pub cols: usize,
    /// Datos row-major: `data[r * cols + c]`.
    pub data: Vec<Scalar>,
}

impl Matrix {
    /// Matriz de ceros de tamaño `rows × cols`.
    pub fn zeros(rows: usize, cols: usize) -> Self {
        Self {
            rows,
            cols,
            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.
///
/// Las dos operaciones que la cohomología necesita (§7): rango del coboundary y
/// una base del cokernel para localizar la obstrucción.
pub trait LinAlg {
    /// `rank(m)` — número de filas/columnas independientes.
    fn rank(&self, m: &Matrix) -> usize;

    /// 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.
#[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 {
        self.rref(m).1.len()
    }

    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:?}");
        }
    }
}