//! `pensamento::combinador` — extraído do monólito na v0.10.0
//! (divisão SRP, crítica Qwen3.8 lote 3 nº 1; r-35 §2).

use super::*;


// ---------------------------------------------------------------------------
// Combinador algébrico (T-P3 — projeção ridge, Teo V-1' portado)
// ---------------------------------------------------------------------------

/// Projeção ridge: `r* = (WᵀW + εI)^{−1} Wᵀ r_alvo` — minimizador de
/// `‖W r − r_alvo‖² + ε‖r‖²` (o montador algébrico do "re-encaixe de peças").
///
/// **Prova (T-P3):** a função é quadrática estritamente convexa (ε > 0);
/// o gradiente anula-se exatamente em `r*` acima; o resíduo mínimo é
/// `‖(I − W(WᵀW+εI)^{−1}Wᵀ) r_alvo‖` — análogo direto do Teo V-1'.
pub fn projetar_ridge(w: &Tensor, r_alvo: &Tensor, eps: f64) -> Result<Tensor> {
    // w: (m, n), r_alvo: (m, d) ⇒ r*: (n, d)
    let n = w.dim(1)?;
    let dev = w.device();
    let wt = w.transpose(0, 1)?;
    let wtw = wt.matmul(w)?;
    let ident = Tensor::eye(n, candle_core::DType::F32, &dev)?.affine(eps, 0.0)?;
    let m = (&wtw + &ident)?;
    let rhs = wt.matmul(r_alvo)?;
    resolver_simetrica(&m, &rhs)
}

/// Resolve `M x = B` para M simétrica (n, n) por eliminação de Gauss com
/// pivoteamento parcial (Cholesky não assumido — M pode ser só PSD+εI).
pub fn resolver_simetrica(m: &Tensor, b: &Tensor) -> Result<Tensor> {
    let n = m.dim(0)?;
    let d = b.dim(1)?;
    let dev = m.device();
    let mut a = m.to_vec2::<f32>()?;
    let mut rhs = b.to_vec2::<f32>()?;
    for col in 0..n {
        // pivoteamento parcial
        let (mut pivo, mut linha_p) = (a[col][col].abs(), col);
        for r in col + 1..n {
            if a[r][col].abs() > pivo {
                pivo = a[r][col].abs();
                linha_p = r;
            }
        }
        if pivo < 1e-12 {
            candle_core::bail!("matriz singular (pivo ~ 0 na coluna {col})");
        }
        a.swap(col, linha_p);
        rhs.swap(col, linha_p);
        for r in col + 1..n {
            let f = a[r][col] / a[col][col];
            for c in col..n {
                a[r][c] -= f * a[col][c];
            }
            for c in 0..d {
                rhs[r][c] -= f * rhs[col][c];
            }
        }
    }
    let mut x = vec![vec![0f32; d]; n];
    for r in (0..n).rev() {
        for c in 0..d {
            let mut s = rhs[r][c];
            for k in r + 1..n {
                s -= a[r][k] * x[k][c];
            }
            x[r][c] = s / a[r][r];
        }
    }
    let flat: Vec<f32> = x.into_iter().flatten().collect();
    Ok(Tensor::from_vec(flat, (n, d), &dev)?)
}