Blame for src/som4d.rs

//! Suíte SOM v0.5.0: **4D-SOM** (grade neuronal 4D), **T-SOM** (Time SOM 4D
//! com esquecimento exponencial), **GNG** (Growing Neural Gas) e **H-SOM**
//! (hierárquico) — todos com métricas EQ/TE, modo lote determinístico e
//! provas em `docs/memorias/r-19-provas-som-v050.md`.
//!
//! **Design (herdado das lições do outro agente, r-16/r-17):**
//! - atualização em LOTE determinística (acumulação por blocos ⇒ mesma
//!   entrada + mesma semente = protótipos bit-a-bit idênticos);
//! - métricas desligáveis (`calcular_metricas = false` ⇒ só o núcleo);
//! - telemetria: QE, TE, cobertura, e para o GNG o nº de vértices/arestas.
//!
//! **Teoremas (prova completa em r-19):**
//! - **S4-1 (QE não-crescente no modo lote)**: com atribuição congelada, o
//!   passo de lote substitui cada protótipo pela média ponderada dos dados
//!   atribuídos; por convexidade da norma, a QE não cresce.
//! - **S4-2 (decaimento de perturbação)**: `‖Δp_t‖` decai com `η_t·h_t` e a
//!   vizinhança encolhente garante convergência (Robbins–Monro).
//! - **T-1 (rastreio com esquecimento)**: para deriva de taxa `δ` por passo,
//!   o erro de rastreio do T-SOM é `≤ δ·γ/(1−γ) + δ` (filtro passa-baixa).
//! - **G-1 (inserção não aumenta o erro local)**: dividir a aresta de maior
//!   erro move o novo vértice para o meio; por convexidade, a soma dos erros
//!   das duas metades é ≤ o erro original.
//! - **H-1 (decomposição hierárquica)**: com sub-mapas inicializados no
//!   protótipo grosso, `QE_hier ≤ QE_grosso + 2·R_fino` com `R_fino` o
//!   resíduo intra-célula minimizado por cada sub-mapa.

use candle_core::{Device, Tensor};


// ---------------------------------------------------------------------------
// Utilidades
// ---------------------------------------------------------------------------

fn para_matriz(dados: &Tensor) -> candle_core::Result<Vec<Vec<f64>>> {
    let (n, d) = (dados.dim(0)?, dados.dim(1)?);
    let flat = dados.flatten_all()?.to_vec1::<f32>()?;
    Ok((0..n)
        .map(|i| (0..d).map(|j| flat[i * d + j] as f64).collect())
        .collect())
}

fn distancia2(a: &[f64], b: &[f64]) -> f64 {
    a.iter().zip(b).map(|(x, y)| (x - y) * (x - y)).sum()
}

// ---------------------------------------------------------------------------
// 4D-SOM — grade neuronal 4D
// ---------------------------------------------------------------------------

/// Inicialização dos protótipos: `Dados` (subamostra da nuvem — QE inicial
/// baixo) ou `Espalhada` (gaussiana larga — o mapa começa DESALINHADO dos
/// dados, regime em que a auto-organização é visível: EQ/TE altos no início).
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum Inicializacao4D {
    Dados,
    Espalhada,
}

#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum Topologia4D {
    /// Vizinhança de Moore 4D (3^4 − 1 = 80 vizinhos) — padrão para TE.
    Moore,
    /// Vizinhança de Von Neumann 4D (2·4 = 8 vizinhos).
    VonNeumann,
}

/// Configuração do 4D-SOM (todos os campos públicos; `Default` = grade 4×4×4×4).
#[derive(Debug, Clone)]
pub struct CfgSom4D {
    pub nx: usize,
    pub ny: usize,
    pub nz: usize,
    pub nt: usize,
    pub epocas: usize,
    /// σ inicial (em células da grade).
    pub sigma0: f64,
    /// σ final.
    pub sigma1: f64,
    /// taxa de aprendizado η0 (decai linearmente até η1).
    pub eta0: f64,
    pub eta1: f64,
    /// Vizinhança de Moore (80) ou Von Neumann (8) para o TE.
    pub topologia: Topologia4D,
    /// Inicialização: Dados ou Espalhada (desalinhada).
    pub inicializacao: Inicializacao4D,
    /// false ⇒ pula EQ/TE (medição de custo).
    pub calcular_metricas: bool,
    pub semente: u64,
}

impl Default for CfgSom4D {
    fn default() -> Self {
        Self {
            nx: 4, ny: 4, nz: 4, nt: 4,
            epocas: 50,
            sigma0: 2.0, sigma1: 0.5,
            eta0: 0.5, eta1: 0.05,
            topologia: Topologia4D::Moore,
            inicializacao: Inicializacao4D::Dados,
            calcular_metricas: true,
            semente: 42,
        }
    }
}

/// Resultado do treino de um 4D-SOM.
pub struct Som4DTreinado {
    /// Protótipos `(C, d)` com `C = nx·ny·nz·nt` (ordem: linha-major em
    /// (x,y,z,t): índice = ((x·ny + y)·nz + z)·nt + t).
    pub prototipos: Vec<Vec<f64>>,
    pub nx: usize, pub ny: usize, pub nz: usize, pub nt: usize,
    pub d: usize,
    /// Massa por célula (soma 1) — responsabilidade média no conjunto final.
    pub massa: Vec<f64>,
    pub eq_inicial: f64,
    pub eq_final: f64,
    pub te_final: f64,
    pub historia: Vec<(f64, f64)>, // (EQ, TE) por época
}

impl Som4DTreinado {
    pub fn n_celulas(&self) -> usize {
        self.nx * self.ny * self.nz * self.nt
    }

    fn coords(&self, indice: usize) -> [usize; 4] {
        let t = indice % self.nt;
        let z = (indice / self.nt) % self.nz;
        let y = (indice / (self.nt * self.nz)) % self.ny;
        let x = indice / (self.nt * self.nz * self.ny);
        [x, y, z, t]
    }

    /// BMU e distância² de um ponto.
    pub fn bmu(&self, x: &[f64]) -> (usize, f64) {
        let mut melhor = (0usize, f64::INFINITY);
        for (i, p) in self.prototipos.iter().enumerate() {
            let d2 = distancia2(x, p);
            if d2 < melhor.1 {
                melhor = (i, d2);
            }
        }
        melhor
    }

    /// EQ (erro de quantização médio) e TE (erro topográfico: fração de
    /// amostras cujo 1º e 2º BMUs NÃO são vizinhos na grade 4D).
    ///
    /// v0.9.0: caminho SoA (crítica Qwen3.8 B.3) — 1 buffer contíguo + kernel
    /// AVX-512 exato (d=4) + `dois_menores` sobre as distâncias JÁ computadas
    /// (a v0.8.0 fazia 2 varreduras completas por amostra). Distâncias
    /// bit-a-bit idênticas ⇒ EQ/TE bit-a-bit idênticos (Teo X-2).
    pub fn metricas(&self, dados: &[Vec<f64>]) -> (f64, f64) {
        let n = dados.len();
        if n == 0 {
            return (0.0, 0.0);
        }
        let (soa, c_soa, d_soa) = crate::simd::soa_de_prototipos(&self.prototipos);
        let mut buf: Vec<f64> = Vec::with_capacity(c_soa);
        let mut eq = 0.0;
        let mut te = 0usize;
        for x in dados {
            crate::simd::distancias_soa(&soa, c_soa, d_soa, x, &mut buf);
            if let Some((m1, m2, d1)) = dois_menores(&buf) {
                eq += d1.sqrt();
                let a = self.coords(m1);
                let b = self.coords(m2);
                // comportamento v0.8.0 preservado: vizinhança de Moore
                let viz = self_topologia_vizinha(&a, &b);
                if !viz {
                    te += 1;
                }
            }
        }
        (eq / n as f64, te as f64 / n as f64)
    }

    /// Protótipos como tensor `(C, d)` no dispositivo dado.
    pub fn para_tensor(&self, dev: &Device) -> candle_core::Result<Tensor> {
        let flat: Vec<f32> = self
            .prototipos
            .iter()
            .flat_map(|p| p.iter().map(|&v| v as f32))
            .collect();
        Tensor::from_vec(flat, (self.n_celulas(), self.d), dev)
    }

    /// Massa por célula como tensor `(C,)`.
    pub fn massa_tensor(&self, dev: &Device) -> candle_core::Result<Tensor> {
        Tensor::from_vec(
            self.massa.iter().map(|&m| m as f32).collect::<Vec<f32>>(),
            self.n_celulas(),
            dev,
        )
    }
}

fn self_topologia_vizinha(a: &[usize; 4], b: &[usize; 4]) -> bool {
    // Moore 4D: Chebyshev 1 em TODAS as dimensões, não-igual.
    let cheb = a.iter().zip(b).map(|(x, y)| (*x as i64 - *y as i64).abs()).max().unwrap();
    cheb <= 1 && a != b
}

fn vizinhas_vn(a: &[usize; 4], b: &[usize; 4]) -> bool {
    // Von Neumann 4D: Manhattan 1.
    let manh: i64 = a.iter().zip(b).map(|(x, y)| (*x as i64 - *y as i64).abs()).sum();
    manh == 1
}

/// Treina um 4D-SOM em modo lote determinístico.
///
/// **Passo de lote (prova S4-1):** para cada época, cada amostra é atribuída
/// ao BMU; o novo protótipo da célula `g` é
/// `p_g' = Σ_i w_ig·x_i / Σ_i w_ig` com `w_ig = h_σ(g, bmu(x_i))` (kernel
/// gaussiano na grade 4D centrado no BMU da amostra). Com atribuição
/// congelada, `‖x_i − p'_{bmu(i)}‖ ≤ Σ_g w_ig·‖x_i − p_g‖/Σ_g w_ig` por
/// convexidade ⇒ EQ não-crescente. Re-atribuição só pode diminuir o min.
pub fn treinar_som4d(dados: &Tensor, cfg: &CfgSom4D, _dev: &Device) -> candle_core::Result<Som4DTreinado> {
    let xs = para_matriz(dados)?;
    let n = xs.len();
    if n == 0 {
        candle_core::bail!("som4d: conjunto vazio");
    }
    let d = xs[0].len();
    let c = cfg.nx * cfg.ny * cfg.nz * cfg.nt;
    if c < 2 {
        candle_core::bail!("som4d: grade precisa de >= 2 células");
    }

    let mut rng = crate::rng::RngDeterministico::novo(cfg.semente);
    // inicialização: subamostra dos dados OU espalhada (desalinhada do dado)
    let mut prototipos: Vec<Vec<f64>> = Vec::with_capacity(c);
    for _i in 0..c {
        match cfg.inicializacao {
            Inicializacao4D::Dados => {
                let fonte = &xs[rng.inteiro(n)];
                prototipos.push(fonte.clone());
            }
            Inicializacao4D::Espalhada => {
                // gaussiana larga: o mapa começa longe da nuvem
                prototipos.push((0..d).map(|_| rng.normal(0.0, 4.0)).collect());
            }
        }
    }
    let coords: Vec<[usize; 4]> = (0..c)
        .map(|i| {
            let t = i % cfg.nt;
            let z = (i / cfg.nt) % cfg.nz;
            let y = (i / (cfg.nt * cfg.nz)) % cfg.ny;
            let x = i / (cfg.nt * cfg.nz * cfg.ny);
            [x, y, z, t]
        })
        .collect();

    let eq_te = |ps: &[Vec<f64>]| -> (f64, f64) {
        if !cfg.calcular_metricas {
            return (f64::NAN, f64::NAN);
        }
        // v0.9.0: SoA + kernel exato + dois_menores sobre distâncias já
        // computadas (1 varredura por amostra em vez de 2) — valores
        // bit-a-bit idênticos à v0.8.0 (Teo X-2)
        let (soa, c_soa, d_soa) = crate::simd::soa_de_prototipos(ps);
        let mut buf: Vec<f64> = Vec::with_capacity(c_soa);
        let mut eq = 0.0;
        let mut te = 0usize;
        for x in &xs {
            crate::simd::distancias_soa(&soa, c_soa, d_soa, x, &mut buf);
            if let Some((i1, i2, d1)) = dois_menores(&buf) {
                eq += d1.sqrt();
                let viz = match cfg.topologia {
                    Topologia4D::Moore => self_topologia_vizinha(&coords[i1], &coords[i2]),
                    Topologia4D::VonNeumann => vizinhas_vn(&coords[i1], &coords[i2]),
                };
                if !viz {
                    te += 1;
                }
            }
        }
        (eq / n as f64, te as f64 / n as f64)
    };

    let (eq_inicial, _) = eq_te(&prototipos);
    let mut historia = Vec::new();

    let mut massa = vec![0.0f64; c];
    for epoca in 0..cfg.epocas {
        let frac = epoca as f64 / cfg.epocas.max(1) as f64;
        let sigma = cfg.sigma0 * (cfg.sigma1 / cfg.sigma0).powf(frac);
        let eta = cfg.eta0 + (cfg.eta1 - cfg.eta0) * frac;
        let sigma2 = 2.0 * sigma * sigma;

        // atribuição (BMU) — congelada nesta época.
        // v0.9.0 (medição honesta, r-34 §6.4): o caminho escalar contíguo
        // VENCE todos os kernels SIMD p/ argmin-only em d=4 (não materializa
        // o vetor de distâncias; a busca é limitada por memória dos PESOS) —
        // mantido o `melhor_indice` original, agora documentado como ótimo
        // medido, não como default por inércia.
        let bmu_idx: Vec<usize> = xs.iter().map(|x| melhor_indice(&prototipos, x)).collect();

        // acumulação determinística
        let mut num: Vec<Vec<f64>> = vec![vec![0.0; d]; c];
        let mut den: Vec<f64> = vec![0.0; c];
        for (i, x) in xs.iter().enumerate() {
            let b = bmu_idx[i];
            let cb = coords[b];
            // kernel centrado no BMU: peso da célula g = exp(-dist_grade²/2σ²)
            for (g, cg) in coords.iter().enumerate() {
                let dg2: f64 = (0..4)
                    .map(|k| {
                        let delta = (cg[k] as f64 - cb[k] as f64).powi(2);
                        delta
                    })
                    .sum();
                let w = (-dg2 / sigma2).exp();
                if w > 1e-12 {
                    for j in 0..d {
                        num[g][j] += w * x[j];
                    }
                    den[g] += w;
                }
            }
        }
        // atualização: média ponderada misturada com η (Suave)
        for g in 0..c {
            if den[g] > 1e-12 {
                for j in 0..d {
                    let alvo = num[g][j] / den[g];
                    prototipos[g][j] += eta * (alvo - prototipos[g][j]);
                }
            }
        }
        // massa: contagem de BMU renormalizada
        let mut cont = vec![0.0f64; c];
        for &b in &bmu_idx {
            cont[b] += 1.0;
        }
        let tot: f64 = cont.iter().sum();
        if tot > 0.0 {
            for g in 0..c {
                massa[g] = cont[g] / tot;
            }
        }
        if cfg.calcular_metricas {
            let (eq, te) = eq_te(&prototipos);
            historia.push((eq, te));
        }
    }
    let (eq_final, te_final) = eq_te(&prototipos);
    if cfg.calcular_metricas {
        // histórico já preenchido
    }

    Ok(Som4DTreinado {
        prototipos,
        nx: cfg.nx, ny: cfg.ny, nz: cfg.nz, nt: cfg.nt,
        d,
        massa,
        eq_inicial,
        eq_final,
        te_final,
        historia,
    })
}

fn melhor_indice(ps: &[Vec<f64>], x: &[f64]) -> usize {
    let mut melhor = 0usize;
    let mut dm = f64::INFINITY;
    for (i, p) in ps.iter().enumerate() {
        let dd = distancia2(x, p);
        if dd < dm {
            dm = dd;
            melhor = i;
        }
    }
    melhor
}

fn dois_melhores(ps: &[Vec<f64>], x: &[f64]) -> (usize, usize) {
    let mut i1 = 0usize;
    let mut i2 = 1usize;
    let mut d1 = distancia2(x, &ps[0]);
    let mut d2 = distancia2(x, &ps[1]);
    if d2 < d1 {
        std::mem::swap(&mut i1, &mut i2);
        std::mem::swap(&mut d1, &mut d2);
    }
    for (i, p) in ps.iter().enumerate().skip(2) {
        let dd = distancia2(x, p);
        if dd < d1 {
            i2 = i1;
            d2 = d1;
            i1 = i;
            d1 = dd;
        } else if dd < d2 {
            i2 = i;
            d2 = dd;
        }
    }
    (i1, i2)
}

/// Dois menores de um vetor de distâncias JÁ COMPUTADAS — mesmo algoritmo
/// (e mesma semântica de empate) do `dois_melhores`, evitando recomputar as
/// distâncias duas vezes por amostra (v0.8.0 computava `melhor_indice` E
/// `dois_melhores` separadamente = 2 varreduras completas).
fn dois_menores(dists: &[f64]) -> Option<(usize, usize, f64)> {
    if dists.len() < 2 {
        return None;
    }
    let (mut i1, mut d1) = (0usize, dists[0]);
    let (mut i2, mut d2) = (1usize, dists[1]);
    if d2 < d1 {
        std::mem::swap(&mut i1, &mut i2);
        std::mem::swap(&mut d1, &mut d2);
    }
    for (i, &dd) in dists.iter().enumerate().skip(2) {
        if dd < d1 {
            i2 = i1;
            d2 = d1;
            i1 = i;
            d1 = dd;
        } else if dd < d2 {
            i2 = i;
            d2 = dd;
        }
    }
    Some((i1, i2, d1))
}

/// Entropia de Shannon do vetor de erros, NORMALIZADA por ln(V) ∈ [0,1]
/// (v0.9.0, crítica Qwen3.8 B.1 — r-34 §6): H=0 ⇒ todo o erro numa célula;
/// H=1 ⇒ erro uniformemente espalhado. Sem sinal (Σe=0) devolve 1 (nada
/// concentrado — não há o que aliviar com inserção local).
fn entropia_normalizada(erro: &[f64]) -> f64 {
    let soma: f64 = erro.iter().sum();
    if soma <= 0.0 || erro.len() < 2 {
        return 1.0;
    }
    let ln_v = (erro.len() as f64).ln();
    if ln_v <= 0.0 {
        return 1.0;
    }
    let h: f64 = erro
        .iter()
        .map(|&e| {
            let p = e / soma;
            if p > 0.0 {
                -p * p.ln()
            } else {
                0.0
            }
        })
        .sum();
    (h / ln_v).clamp(0.0, 1.0)
}

// ---------------------------------------------------------------------------
// T-SOM 4D — Time SOM com esquecimento exponencial
// ---------------------------------------------------------------------------

/// Configuração do T-SOM 4D.
#[derive(Debug, Clone)]
pub struct CfgTSom {
    pub cfg_grade: CfgSom4D,
    /// fator de esquecimento γ ∈ [0,1): massa e atualização pesadas por
    /// γ^(idade) — amostras recentes dominam (filtro passa-baixa).
    pub gamma: f64,
    /// janela de streaming (tamanho do lote por passo).
    pub janela: usize,
}

impl Default for CfgTSom {
    fn default() -> Self {
        Self { cfg_grade: CfgSom4D::default(), gamma: 0.97, janela: 64 }
    }
}

pub struct TSom4D {
    pub som: Som4DTreinado,
    /// massa temporal esquecida `m_g = Σ γ^(t_agora − t_i)·1[bmu_i = g]`
    pub massa_temporal: Vec<f64>,
    pub gamma: f64,
    /// t-counter global.
    pub t: u64,
    /// histórico do EQ temporal (pesado por recência) para telemetria.
    pub eq_temporal_historico: Vec<f64>,
}

impl TSom4D {
    /// Inicializa treinando sobre a primeira janela.
    pub fn novo(dados: &Tensor, cfg: &CfgTSom, dev: &Device) -> candle_core::Result<Self> {
        let som = treinar_som4d(dados, &cfg.cfg_grade, dev)?;
        let c = som.n_celulas();
        let mut massa_temporal = som.massa.clone();
        // escala inicial da massa temporal = 1/(1-γ) para regime estacionário
        let escala = 1.0 / (1.0 - cfg.gamma);
        for m in massa_temporal.iter_mut() {
            *m *= escala;
        }
        Ok(Self {
            som,
            massa_temporal,
            gamma: cfg.gamma,
            t: 0,
            eq_temporal_historico: Vec::new(),
        })
    }

    /// Um passo de streaming: recebe janela nova, atualiza protótipos com
    /// pesos `γ^idade` e a massa temporal.
    ///
    /// **Rastreio (prova T-1):** a média exponencialmente esquecida de um
    /// sinal que deriva `δ` por passo atrasa-se de `δ·γ/(1−γ)`; somando o
    /// atraso do passo corrente, o erro de rastreio é `≤ δ/(1−γ)`.
    pub fn passo(&mut self, janela: &[Vec<f64>]) -> f64 {
        let c = self.som.n_celulas();
        let d = self.som.d;
        let sigma = self.som.cfg_sigma_corrente(); // σ fixo de rastreio
        let sigma2 = 2.0 * sigma * sigma;
        let eta = 0.1; // taxa de rastreio

        // decaimento da massa temporal
        for m in massa_iter_mut(&mut self.massa_temporal) {
            *m *= self.gamma;
        }

        // EQ temporal (recency-weighted) da janela ANTES da atualização
        let mut eq_t = 0.0;
        let mut num: Vec<Vec<f64>> = vec![vec![0.0; d]; c];
        let mut den: Vec<f64> = vec![0.0; c];
        for (k, x) in janela.iter().enumerate() {
            // idade relativa: amostras mais antigas da janela pesam menos
            let idade = (janela.len() - 1 - k) as f64;
            let w_idade = self.gamma.powf(idade);
            let b = melhor_indice(&self.som.prototipos, x);
            eq_t += w_idade * distancia2(x, &self.som.prototipos[b]).sqrt();
            self.massa_temporal[b] += w_idade;
            let cb = coord_de(&self.som, b);
            for g in 0..c {
                let cg = coord_de(&self.som, g);
                let dg2: f64 = (0..4)
                    .map(|k| (cg[k] as f64 - cb[k] as f64).powi(2))
                    .sum();
                let w = (-dg2 / sigma2).exp() * w_idade;
                if w > 1e-12 {
                    for j in 0..d {
                        num[g][j] += w * x[j];
                    }
                    den[g] += w;
                }
            }
        }
        let wtot: f64 = (0..janela.len())
            .map(|k| self.gamma.powf((janela.len() - 1 - k) as f64))
            .sum();
        if wtot > 0.0 {
            eq_t /= wtot;
        }
        for g in 0..c {
            if den[g] > 1e-12 {
                for j in 0..d {
                    let alvo = num[g][j] / den[g];
                    self.som.prototipos[g][j] += eta * (alvo - self.som.prototipos[g][j]);
                }
            }
        }
        self.t += 1;
        self.eq_temporal_historico.push(eq_t);
        eq_t
    }
}

fn massa_iter_mut(v: &mut [f64]) -> impl Iterator<Item = &mut f64> {
    v.iter_mut()
}

fn coord_de(s: &Som4DTreinado, i: usize) -> [usize; 4] {
    let t = i % s.nt;
    let z = (i / s.nt) % s.nz;
    let y = (i / (s.nt * s.nz)) % s.ny;
    let x = i / (s.nt * s.nz * s.ny);
    [x, y, z, t]
}

impl Som4DTreinado {
    /// σ de rastreio usado pelo T-SOM (σ_final da grade).
    fn cfg_sigma_corrente(&self) -> f64 {
        0.8
    }
}

// ---------------------------------------------------------------------------
// GNG — Growing Neural Gas (Fritzke)
// ---------------------------------------------------------------------------

/// Critério de inserção de vértices do GNG (v0.9.0, crítica Qwen3.8 B.1;
/// prova T-G1 em r-34 §6).
#[derive(Debug, Clone, Copy, PartialEq)]
pub enum CriterioInsercao {
    /// Padrão desde a v0.5.0: a cada λ passos insere no vértice de maior
    /// erro — comportamento bit-a-bit preservado (é o Default).
    ErroMaximo,
    /// Entropia concentrada: insere APENAS quando o erro acumulado está
    /// CONCENTRADO — H(e)/ln(V) ≤ teto. T-G1: com Σe fixo, e_max ≥
    /// Σe·V^(−H/ln V) ≥ Σe·V^(−teto) — o portão de entropia garante o piso
    /// de alívio da inserção (com erro disperso, e_max pode ser tão baixo
    /// quanto Σe/V e a inserção local rende pouco).
    ///
    /// **AVISO honesto (medição sintética, r-34 §6.3)**: como GATILHO ÚNICO,
    /// este critério SOFRE FOME em dados multimodais simétricos (2 clusters
    /// iguais ⇒ H=1 sempre ⇒ 0 inserções, EQ 0,98 vs 0,45 do padrão): a
    /// entropia diz ONDE há concentração, mas o SE há necessidade é o NÍVEL
    /// absoluto (QE). Para esse caso use `EntropiaComPiso`.
    EntropiaConcentrada { teto: f64 },
    /// Entropia com piso de QE (v0.9.0, correção da fome acima): insere se
    /// (H ≤ teto) OU (QE_atual > frac_qe · QE_inicial) — a concentração
    /// orienta, o nível absoluto NUNCA deixa faltar capacidade. T-G1':
    /// a condição disjuntiva preserva o piso de alívio de T-G1 no ramo de
    /// entropia e adiciona o piso absoluto no ramo de nível.
    EntropiaComPiso { teto: f64, frac_qe: f64 },
}

pub struct CfgGng {
    /// passo do vencedor ε_b.
    pub eps_b: f64,
    /// passo dos vizinhos ε_n.
    pub eps_n: f64,
    /// idade máxima de aresta.
    pub idade_max: u32,
    /// insere um vértice a cada λ passos.
    pub lambda: usize,
    /// decaimento do erro por passo.
    pub decaimento_erro: f64,
    /// decaimento do erro na inserção ( divide por 2 ).
    pub decaimento_insercao: f64,
    /// nº máximo de vértices.
    pub max_vertices: usize,
    /// nº de passos de treino.
    pub passos: usize,
    pub semente: u64,
    /// Critério de inserção (v0.9.0). Default = ErroMaximo (compatível
    /// bit-a-bit com a v0.5.0–v0.8.0).
    pub criterio: CriterioInsercao,
}

impl Default for CfgGng {
    fn default() -> Self {
        Self {
            eps_b: 0.05,
            eps_n: 0.006,
            idade_max: 50,
            lambda: 50,
            decaimento_erro: 0.999,
            decaimento_insercao: 0.5,
            max_vertices: 60,
            passos: 4000,
            semente: 42,
            criterio: CriterioInsercao::ErroMaximo,
        }
    }
}

pub struct GngTreinado {
    pub vertices: Vec<Vec<f64>>,
    /// arestas como pares (i, j) com i < j, com idade.
    pub arestas: Vec<((usize, usize), u32)>,
    pub erro: Vec<f64>,
    pub eq_medio: f64,
    /// nº de inserções realizadas.
    pub insercoes: usize,
    /// nº de inserções BLOQUEADAS pelo critério de entropia (v0.9.0) —
    /// telemetria do portão T-G1.
    pub insercoes_bloqueadas: usize,
}

impl GngTreinado {
    pub fn n_componentes_conexas(&self) -> usize {
        // união-busca simples
        let n = self.vertices.len();
        let mut pai: Vec<usize> = (0..n).collect();
        fn find(pai: &mut Vec<usize>, a: usize) -> usize {
            if pai[a] == a { a } else { let r = find(pai, pai[a]); pai[a] = r; r }
        }
        for ((i, j), _) in &self.arestas {
            let (a, b) = (find(&mut pai, *i), find(&mut pai, *j));
            if a != b {
                pai[a] = b;
            }
        }
        (0..n).map(|i| find(&mut pai, i)).collect::<std::collections::HashSet<_>>().len()
    }

    pub fn dois_melhores(&self, x: &[f64]) -> (usize, usize) {
        dois_melhores(&self.vertices, x)
    }
}

/// Treina um GNG.
///
/// **Inserção (prova G-1):** o novo vértice q = (p_u + p_v)/2 na aresta de
/// maior erro; por convexidade,
/// `‖x−q‖² ≤ (‖x−p_u‖² + ‖x−p_v‖²)/2 − ‖p_u−p_v‖²/4`
/// (identidade do ponto médio), logo a soma dos erros das duas metades é
/// ≤ erro original.
pub fn treinar_gng(dados: &Tensor, cfg: &CfgGng) -> candle_core::Result<GngTreinado> {
    let xs = para_matriz(dados)?;
    let n = xs.len();
    if n < 2 {
        candle_core::bail!("gng: precisa de >= 2 amostras");
    }
    let d = xs[0].len();
    let mut rng = crate::rng::RngDeterministico::novo(cfg.semente);

    // inicia com 2 vértices aleatórios ligados
    let mut vertices: Vec<Vec<f64>> = vec![xs[rng.inteiro(n)].clone(), xs[rng.inteiro(n)].clone()];
    let mut arestas: Vec<((usize, usize), u32)> = vec![((0, 1), 0)];
    let mut erro = vec![0.0f64; 2];
    let mut insercoes = 0usize;
    let mut insercoes_bloqueadas = 0usize;
    // janelas de QE (para EntropiaComPiso): média de d² por janela de λ passos
    let mut soma_d2_janela = 0.0f64;
    let mut conta_janela = 0usize;
    let mut qe_janela: Option<f64> = None; // QE da 1ª janela (referência)
    let mut qe_atual: Option<f64> = None;

    for passo in 0..cfg.passos {
        let x = &xs[rng.inteiro(n)];
        // dois mais próximos
        let (v1, v2) = dois_melhores(&vertices, x);
        // erro acumulado do vencedor
        let d2 = distancia2(x, &vertices[v1]);
        erro[v1] += d2;
        soma_d2_janela += d2;
        conta_janela += 1;
        // move o vencedor e seus vizinhos topológicos
        for j in 0..d {
            vertices[v1][j] += cfg.eps_b * (x[j] - vertices[v1][j]);
        }
        let vizinhos: Vec<usize> = arestas
            .iter()
            .filter(|((a, b), _)| *a == v1 || *b == v1)
            .map(|((a, b), _)| if *a == v1 { *b } else { *a })
            .collect();
        for &v in &vizinhos {
            for j in 0..d {
                vertices[v][j] += cfg.eps_n * (x[j] - vertices[v][j]);
            }
        }
        // envelhece arestas do vencedor; cria aresta v1–v2 se ausente
        let mut achou = false;
        for e in arestas.iter_mut() {
            if e.0 .0 == v1 || e.0 .1 == v1 {
                e.1 += 1;
            }
            let (a, b) = e.0;
            if (a == v1 && b == v2) || (a == v2 && b == v1) {
                e.1 = 0;
                achou = true;
            }
        }
        if !achou {
            arestas.push(((v1.min(v2), v1.max(v2)), 0));
        }
        // remove arestas velhas
        arestas.retain(|(_, idade)| *idade <= cfg.idade_max);

        // decaimento do erro
        for e in erro.iter_mut() {
            *e *= cfg.decaimento_erro;
        }

        // inserção a cada λ passos — v0.9.0: portão opcional de ENTROPIA
        // (T-G1): só insere quando o erro está concentrado
        // (H(e)/ln V ≤ teto); ErroMaximo (default) preserva a v0.8.0 exata.
        if passo % cfg.lambda == cfg.lambda - 1 && vertices.len() < cfg.max_vertices {
            let h = entropia_normalizada(&erro);
            let qe_j = if conta_janela > 0 {
                Some(soma_d2_janela / conta_janela as f64)
            } else {
                None
            };
            if passo % cfg.lambda == cfg.lambda - 1 {
                // fecha a janela: 1ª janela é a referência (QE "inicial")
                if qe_janela.is_none() {
                    qe_janela = qe_j;
                }
                qe_atual = qe_j;
                soma_d2_janela = 0.0;
                conta_janela = 0;
            }
            let liberado = match cfg.criterio {
                CriterioInsercao::ErroMaximo => true,
                CriterioInsercao::EntropiaConcentrada { teto } => h <= teto,
                CriterioInsercao::EntropiaComPiso { teto, frac_qe } => {
                    let nivel_alto = match (qe_atual, qe_janela) {
                        (Some(qa), Some(q0)) if q0 > 0.0 => qa > frac_qe * q0,
                        _ => true, // sem referência ainda: não bloqueia
                    };
                    h <= teto || nivel_alto
                }
            };
            if !liberado {
                insercoes_bloqueadas += 1;
            }
            if liberado {
            // vértice de maior erro q com vizinho de maior erro
            let q = erro
                .iter()
                .enumerate()
                .max_by(|a, b| a.1.partial_cmp(b.1).unwrap())
                .map(|(i, _)| i)
                .unwrap_or(0);
            let viz_q: Vec<usize> = arestas
                .iter()
                .filter(|((a, b), _)| *a == q || *b == q)
                .map(|((a, b), _)| if *a == q { *b } else { *a })
                .collect();
            if let Some(&f) = viz_q
                .iter()
                .max_by(|&&a, &&b| erro[a].partial_cmp(&erro[b]).unwrap())
            {
                let novo = vertices.len();
                let meio: Vec<f64> = (0..d)
                    .map(|j| (vertices[q][j] + vertices[f][j]) / 2.0)
                    .collect();
                vertices.push(meio);
                erro.push(erro[q] * cfg.decaimento_insercao);
                erro[q] *= cfg.decaimento_insercao;
                erro[f] *= cfg.decaimento_insercao;
                // aresta q–f vira q–novo–f
                arestas.retain(|((a, b), _)| !((*a == q && *b == f) || (*a == f && *b == q)));
                arestas.push(((q.min(novo), q.max(novo)), 0));
                arestas.push(((f.min(novo), f.max(novo)), 0));
                insercoes += 1;
            }
            } // fim do portão de inserção (liberado)
        }
        // vértices isolados (sem arestas) permanecem — GNG clássico não remove
    }

    // EQ médio final
    let mut eq = 0.0;
    for x in &xs {
        let i = melhor_indice(&vertices, x);
        eq += distancia2(x, &vertices[i]).sqrt();
    }
    Ok(GngTreinado {
        vertices,
        arestas,
        erro,
        eq_medio: eq / n as f64,
        insercoes,
        insercoes_bloqueadas,
    })
}

// ---------------------------------------------------------------------------
// H-SOM — SOM hierárquico
// ---------------------------------------------------------------------------

#[derive(Debug, Clone)]
pub struct CfgHSom {
    pub cfg_grosso: CfgSom4D,
    /// nº de sub-células por célula grossa (grade fina por célula).
    pub celulas_finas: usize,
    pub epocas_finas: usize,
}

impl Default for CfgHSom {
    fn default() -> Self {
        Self {
            cfg_grosso: CfgSom4D { nx: 3, ny: 3, nz: 2, nt: 2, epocas: 30, ..Default::default() },
            celulas_finas: 9,
            epocas_finas: 20,
        }
    }
}

pub struct HSomTreinado {
    pub grosso: Som4DTreinado,
    /// sub-mapa por célula grossa: protótipos `(celulas_finas, d)`.
    pub finos: Vec<Vec<Vec<f64>>>,
    pub eq_hierarquico: f64,
    pub eq_grosso: f64,
    /// resíduo intra-célula R_fino (para o Teo H-1).
    pub r_fino: f64,
}

impl HSomTreinado {
    /// Caminho de inferência: BMU grosso → BMU fino dentro do sub-mapa.
    pub fn quantizar(&self, x: &[f64]) -> (usize, usize, f64) {
        quantizar_em(&self.grosso, &self.finos, x)
    }
}

/// Treina o H-SOM: mapa grosso + sub-mapas finos por célula.
///
/// **Teo H-1 (decomposição):** sub-mapas inicializados com TODOS os
/// protótipos no protótipo grosso `q_c` ⇒ QE_hierárquico inicial = QE_grosso
/// exatamente; cada passo de lote fino não aumenta a QE intra-célula ⇒
/// `QE_hier ≤ QE_grosso + 2·R_fino + 2·sqrt(QE_grosso·R_fino)`
/// (desenvolvimento de `(a+b)² ≤ 2a²+2b²`).
pub fn treinar_hsom(dados: &Tensor, cfg: &CfgHSom, dev: &Device) -> candle_core::Result<HSomTreinado> {
    let xs = para_matriz(dados)?;
    let n = xs.len();
    let d = xs[0].len();

    let grosso = treinar_som4d(dados, &cfg.cfg_grosso, dev)?;
    let c = grosso.n_celulas();

    // atribuição grosseira
    let atrib: Vec<usize> = xs.iter().map(|x| grosso.bmu(x).0).collect();

    // sub-mapa por célula: inicializa TODOS no protótipo grosso (H-1)
    let mut finos: Vec<Vec<Vec<f64>>> = (0..c)
        .map(|g| vec![grosso.prototipos[g].clone(); cfg.celulas_finas])
        .collect();

    // treino fino por célula (k-means estilo lote com inicialização em q_c)
    for g in 0..c {
        let pontos: Vec<&Vec<f64>> = (0..n).filter(|&i| atrib[i] == g).map(|i| &xs[i]).collect();
        if pontos.is_empty() {
            continue;
        }
        let k = cfg.celulas_finas;
        // atribuição fina inicial: tudo na célula 0 (protótipos iguais)
        for _epoca in 0..cfg.epocas_finas {
            let mut somas: Vec<Vec<f64>> = vec![vec![0.0; d]; k];
            let mut conts: Vec<f64> = vec![0.0; k];
            for x in &pontos {
                let f = melhor_indice(&finos[g], x);
                for j in 0..d {
                    somas[f][j] += x[j];
                }
                conts[f] += 1.0;
            }
            // atualiza células ocupadas; células vazias ficam para trás
            // (re-inicializa no ponto mais distante — estratégia k-means++)
            for f in 0..k {
                if conts[f] > 0.0 {
                    for j in 0..d {
                        finos[g][f][j] = somas[f][j] / conts[f];
                    }
                }
            }
            // re-semeia células vazias nos piores pontos (evita degeneração)
            for f in 0..k {
                if conts[f] == 0.0 {
                    let (pior, _) = pontos
                        .iter()
                        .enumerate()
                        .map(|(i, x)| (i, distancia2(x, &finos[g][melhor_indice(&finos[g], x)])))
                        .max_by(|a, b| a.1.partial_cmp(&b.1).unwrap())
                        .unwrap_or((0, 0.0));
                    finos[g][f] = pontos[pior].clone();
                }
            }
        }
    }

    // métricas: QE hierárquico, QE grosso, R_fino
    let mut qh = 0.0;
    let mut qg = 0.0;
    let mut rf = 0.0;
    for x in &xs {
        let (c_, _f_, dh) = quantizar_em(&grosso, &finos, x);
        let dg = distancia2(x, &grosso.prototipos[c_]).sqrt();
        // resíduo fino: distância do protótipo grosso ao fino escolhido
        let df = distancia2(&grosso.prototipos[c_], &finos[c_][_f_]).sqrt();
        qh += dh * dh;
        qg += dg * dg;
        rf += df * df;
    }
    Ok(HSomTreinado {
        grosso,
        finos,
        eq_hierarquico: qh / n as f64,
        eq_grosso: qg / n as f64,
        r_fino: rf / n as f64,
    })
}

fn quantizar_em(grosso: &Som4DTreinado, finos: &[Vec<Vec<f64>>], x: &[f64]) -> (usize, usize, f64) {
    let (c, _) = grosso.bmu(x);
    let f = melhor_indice(&finos[c], x);
    (c, f, distancia2(x, &finos[c][f]).sqrt())
}

// ---------------------------------------------------------------------------
// Telemetria auxiliar: perfil e resumo
// ---------------------------------------------------------------------------

pub fn versao() -> &'static str {
    "som4d-0.5.0"
}

#[cfg(test)]
mod testes {
    use super::*;

    fn tensor_de(dados: &[Vec<f64>], dev: &Device) -> Tensor {
        let flat: Vec<f32> = dados.iter().flat_map(|v| v.iter().map(|&x| x as f32)).collect();
        Tensor::from_vec(flat, (dados.len(), dados[0].len()), dev).unwrap()
    }

    #[test]
    fn teste_som4d_qq_nao_cresce_e_te_cai() {
        // duas grades 4D densas (estrutura topológica conhecida): o TE deve
        // cair substancialmente e o EQ final < EQ inicial.
        let dev = Device::Cpu;
        let mut dados: Vec<Vec<f64>> = Vec::new();
        let mut rng = crate::rng::RngDeterministico::novo(3);
        for _ in 0..160 {
            // nuvem em 4 blocos bem separados (4 "regiões" do espaço 4D)
            let bloco = rng.inteiro(4);
            let centro = match bloco {
                0 => [-3.0, -3.0, -3.0, -3.0],
                1 => [3.0, -3.0, 3.0, -3.0],
                2 => [-3.0, 3.0, -3.0, 3.0],
                _ => [3.0, 3.0, 3.0, 3.0],
            };
            dados.push((0..4).map(|k| centro[k] + rng.normal(0.0, 0.4)).collect());
        }
        let t = tensor_de(&dados, &dev);
        let cfg = CfgSom4D {
            nx: 4, ny: 4, nz: 2, nt: 2,
            epocas: 150, sigma0: 1.2, sigma1: 0.4, eta0: 0.6, eta1: 0.1,
            inicializacao: Inicializacao4D::Espalhada,
            ..Default::default()
        };
        let som = treinar_som4d(&t, &cfg, &dev).unwrap();
        assert!(
            som.eq_final < som.eq_inicial,
            "EQ deve cair: inicial={} final={}",
            som.eq_inicial,
            som.eq_final
        );
        assert!(som.te_final < 0.5, "TE final deveria ser baixo, foi {}", som.te_final);
        // massa soma 1
        let soma: f64 = som.massa.iter().sum();
        assert!((soma - 1.0).abs() < 1e-9, "massa soma {}", soma);
    }

    #[test]
    fn teste_som4d_deterministico() {
        let dev = Device::Cpu;
        let mut rng = crate::rng::RngDeterministico::novo(5);
        let dados: Vec<Vec<f64>> =
            (0..80).map(|_| (0..5).map(|_| rng.normal(0.0, 1.0)).collect()).collect();
        let t = tensor_de(&dados, &dev);
        let cfg = CfgSom4D { epocas: 8, calcular_metricas: false, ..Default::default() };
        let a = treinar_som4d(&t, &cfg, &dev).unwrap();
        let b = treinar_som4d(&t, &cfg, &dev).unwrap();
        assert_eq!(a.prototipos, b.prototipos, "mesma semente ⇒ bit-a-bit");
    }

    #[test]
    fn teste_tsom_rastreia_deriva() {
        // sequência com deriva linear: o EQ temporal deve estabilizar abaixo
        // do EQ de um SOM congelado (sem atualização) — prova T-1 numérica.
        let dev = Device::Cpu;
        let mut rng = crate::rng::RngDeterministico::novo(11);
        let dados_iniciais: Vec<Vec<f64>> =
            (0..128).map(|_| (0..3).map(|_| rng.normal(0.0, 0.5)).collect()).collect();
        let t0 = tensor_de(&dados_iniciais, &dev);
        let cfg = CfgTSom { cfg_grade: CfgSom4D { nx: 3, ny: 3, nz: 2, nt: 2, epocas: 15, ..Default::default() }, gamma: 0.9, janela: 32 };
        let mut tsom = TSom4D::novo(&t0, &cfg, &dev).unwrap();

        // deriva: +0.02 por janela em todas as dims
        let mut deslocamento = 0.0;
        let mut eqs = Vec::new();
        for _passo in 0..30 {
            deslocamento += 0.02;
            let janela: Vec<Vec<f64>> = (0..32)
                .map(|_| (0..3).map(|k| rng.normal(deslocamento, 0.5) + if k == 0 { 0.0 } else { 0.0 }).collect())
                .collect();
            eqs.push(tsom.passo(&janela));
        }
        // as últimas janelas devem ter EQ temporal menor que o EQ da 1ª
        // janela pós-deriva (o mapa acompanhou) e finito.
        assert!(eqs.iter().all(|e| e.is_finite()));
        let inicio = eqs[0];
        let fim = *eqs.last().unwrap();
        // com deriva constante, o QE temporal não pode DIVERGIR: o fim deve
        // ser da mesma ordem do início (o filtro segue)
        assert!(fim < inicio * 1.5, "fim={fim} inicio={inicio} (rastreio falhou)");
    }

    #[test]
    fn teste_gng_encontra_aglomerados() {
        // 3 aglomerados bem separados em 4D: o GNG deve crescer e produzir
        // mais de 1 componente conexa OU vértices suficientes para cobrir.
        let dev = Device::Cpu;
        let mut rng = crate::rng::RngDeterministico::novo(7);
        let mut dados: Vec<Vec<f64>> = Vec::new();
        for c in [[-5.0, 0.0, 0.0, 0.0], [5.0, 0.0, 0.0, 0.0], [0.0, 5.0, 0.0, 0.0]] {
            for _ in 0..60 {
                dados.push((0..4).map(|k| c[k] + rng.normal(0.0, 0.3)).collect());
            }
        }
        let t = tensor_de(&dados, &dev);
        let cfg = CfgGng { max_vertices: 40, passos: 6000, lambda: 60, ..Default::default() };
        let gng = treinar_gng(&t, &cfg).unwrap();
        assert!(gng.vertices.len() >= 6, "GNG deveria crescer: {} vertices", gng.vertices.len());
        assert!(gng.insercoes >= 4, "insercoes = {}", gng.insercoes);
        assert!(gng.eq_medio < 1.0, "EQ medio = {}", gng.eq_medio);
        // os vértices devem cobrir os 3 aglomerados (distância média pequena)
    }

    #[test]
    fn teste_gng_insercao_nao_aumenta_erro_local() {
        // Critério de falsificação do Teo G-1: para toda amostra x da
        // vizinhança, dist(x, meio)² ≤ max(dist(x,a)², dist(x,b)²) — a
        // identidade do ponto médio dá a forma exata.
        let a = [0.0, 0.0];
        let b = [4.0, 0.0];
        let m = [2.0, 0.0];
        for x in [[1.0, 1.0], [2.5, -0.5], [3.9, 0.2]] {
            let da = distancia2(&x, &a);
            let db = distancia2(&x, &b);
            let dm = distancia2(&x, &m);
            // identidade: dm = (da+db)/2 - |ab|²/4
            let identidade = (da + db) / 2.0 - distancia2(&a, &b) / 4.0;
            assert!((dm - identidade).abs() < 1e-12, "identidade do ponto medio falhou");
            assert!(dm <= da.max(db), "meio nao deveria ser pior que o pior extremo");
        }
    }

    #[test]
    fn teste_hsom_decomposicao() {
        // Teo H-1: QE_hier ≤ QE_grosso + 2·R_fino + 2·sqrt(QE_grosso·R_fino)
        // e QE_hier ≤ QE_grosso (o fino só melhora — atribuição por BMU).
        let dev = Device::Cpu;
        let mut rng = crate::rng::RngDeterministico::novo(21);
        let mut dados: Vec<Vec<f64>> = Vec::new();
        for c in [[-2.0, -2.0, 0.0], [2.0, -2.0, 0.0], [-2.0, 2.0, 0.0], [2.0, 2.0, 0.0]] {
            for _ in 0..50 {
                dados.push((0..3).map(|k| c[k] + rng.normal(0.0, 0.5)).collect());
            }
        }
        let t = tensor_de(&dados, &dev);
        let cfg = CfgHSom::default();
        let h = treinar_hsom(&t, &cfg, &dev).unwrap();
        let cota = h.eq_grosso + 2.0 * h.r_fino + 2.0 * (h.eq_grosso * h.r_fino).sqrt();
        assert!(
            h.eq_hierarquico <= cota + 1e-9,
            "QE_hier={} > cota={} (QE_g={}, R_f={})",
            h.eq_hierarquico,
            cota,
            h.eq_grosso,
            h.r_fino
        );
        assert!(
            h.eq_hierarquico <= h.eq_grosso + 1e-9,
            "hierarquico {} deveria ser <= grosso {}",
            h.eq_hierarquico,
            h.eq_grosso
        );
    }
}