//! **Verificador de confiança do próprio NCM** (v0.6.0) — substitui o
//! verificador ORACULAR da demo do Pensamento (débito honesto r-21 §3.1).
//!
//! O oráculo da v0.5.0 media `1/(1+10·MSE)` contra o ALVO VERDADEIRO —
//! isola a agregação, mas não existe em produção. O verificador NCM usa
//! apenas a própria distribuição preditiva do backbone:
//! `V(z) = 1 − H(softmax(z/T*))/log K` com temperatura T* calibrada em
//! dados retidos. Provas completas em `docs/memorias/r-23-*`; síntese:
//!
//! - **V-1 (calibração ≡ minimização de KL):** no nível populacional,
//!   `min_T NLL(T) = min_T H(p_verd, p_T) = min_T [H(p_verd) + D_KL(p_verd‖p_T)]`
//!   ⇔ `min_T D_KL(p_verd ‖ p_T)` na família de temperatura (termo
//!   constante subtraído). Busca por golden-section em `log T` no
//!   colchete [ln(1/8), ln 8]; a NLL é contínua no colchete compacto ⇒
//!   mínimo existente; unimodalidade verificada empiricamente (teste).
//! - **V-2 (soundness SEM oráculo):** sob calibração perfeita
//!   `P(erro | conf = c) = 1 − c` (probabilidade total sobre o nível de
//!   confiança), logo `P(erro | aceito) ≤ 1 − τ_c`. Com erro de calibração
//!   medido `ECE = Σ_b (n_b/n)|acc(b) − conf(b)|`:
//!   `P(erro | conf ≥ τ_c) ≤ 1 − τ_c + ECE` — desigualdade triangular
//!   sobre o conjunto de aceitação (união de faixas).
//! - **V-3 (quantil com garantia DKW):** o limiar empírico
//!   `τ_c = quantil_{1−ε}(confs retidas)` satisfaz a bound de
//!   Dvoretzky–Kiefer–Wolfowitz–Massart `P(sup|F̂−F| > δ) ≤ 2e^{−2nδ²}`;
//!   a taxa de erro do "aceitar" é ≤ ε + 2e^{−2nδ²} com prob. ≥ 1−2e^{−2nδ²}.
//! - **V-4 (causalidade):** `V` é função mensurável da passada corrente
//!   do backbone (logits do presente) — nenhum acesso ao futuro ou ao
//!   alvo (mesma estrutura do gate G1/T-P1).

use candle_core::{Device, Result, Tensor};

/// Confiança por linha de uma matriz de logits: `1 − H(p)/log K`.
/// (n, K) → (n,). Diferenciável (exceto medida zero — T-P1).
///
/// **Forma numericamente estável (V-4):** `H = logsumexp(z) − Σ p·z`
/// (identidade algébrica de `−Σ p log p` com `p = softmax(z)` — evita
/// `0·log 0 = NaN` quando o softmax satura em 0 exato).
pub fn confianca_por_linha(logits: &Tensor) -> Result<Tensor> {
    let p = candle_nn::ops::softmax(logits, candle_core::D::Minus1)?;
    // logsumexp estável: m + log Σ exp(z − m)
    let m = logits.max_keepdim(candle_core::D::Minus1)?;
    let desloc = logits.broadcast_sub(&m)?;
    let exp_d = desloc.exp()?;
    let soma_exp = exp_d.sum_keepdim(candle_core::D::Minus1)?;
    let lse = (&soma_exp.log()? + &m)?;
    // Σ p·z (0·z = 0 exato — sem NaN)
    let pz = (&p * logits)?.sum_keepdim(candle_core::D::Minus1)?;
    let ent = (&lse - &pz)?; // H = logsumexp(z) − Σ p·z (identidade exata)
    let k = logits.dim(candle_core::D::Minus1)? as f64;
    if k <= 1.0 {
        return Ok(Tensor::ones(ent.dims(), candle_core::DType::F32, logits.device())?);
    }
    Ok(ent.affine(-1.0 / k.ln(), 1.0)?)
}

/// NLL média de logits escalados por 1/T contra alvos u32 (n,) — (n, K).
fn nll_temperatura(logits: &Tensor, alvos: &Tensor, t: f64) -> Result<f64> {
    let escalado = logits.affine(1.0 / t, 0.0)?;
    let ce = candle_nn::loss::cross_entropy(&escalado, alvos)?;
    Ok(ce.to_scalar::<f32>()? as f64)
}

/// **V-1:** calibra T* por golden-section em `log T ∈ [ln(1/8), ln 8]`
/// minimizando a NLL nos dados retidos. Devolve (T*, nll(1), nll(T*)).
pub fn calibrar_temperatura(logits: &Tensor, alvos: &Tensor) -> Result<(f64, f64, f64)> {
    let nll_em = |t: f64| nll_temperatura(logits, alvos, t);
    let nll_1 = nll_em(1.0)?;
    // golden-section em log T (invariante de escala da parametrização)
    let (mut lo, mut hi) = ((1.0f64 / 8.0).ln(), 8.0f64.ln());
    let phi = (5f64.sqrt() - 1.0) / 2.0;
    let mut a = hi - phi * (hi - lo);
    let mut b = lo + phi * (hi - lo);
    let (mut fa, mut fb) = (nll_em(a.exp())?, nll_em(b.exp())?);
    for _ in 0..64 {
        if fa < fb {
            hi = b;
            b = a;
            fb = fa;
            a = hi - phi * (hi - lo);
            fa = nll_em(a.exp())?;
        } else {
            lo = a;
            a = b;
            fa = fb;
            b = lo + phi * (hi - lo);
            fb = nll_em(b.exp())?;
        }
        if (hi - lo) < 1e-8 {
            break;
        }
    }
    let t_estrela = ((lo + hi) / 2.0).exp();
    let nll_estrela = nll_em(t_estrela)?;
    Ok((t_estrela, nll_1, nll_estrela))
}

/// Derivadas analíticas da NLL de temperatura — base do Newton-Raphson
/// (proposta Qwen-3.8 nº 1, **CORRIGIDA**: o original inverteu o sinal de
/// dL/dT — `(E_p[z] − z_y)/T²` em vez de `(z_y − E_p[z])/T²` — e usou a
/// primeira diferença `[L(T+ε)−L(T)]/ε` como se fosse a SEGUNDA derivada;
/// com os dois erros combinados o passo Δ ≈ −1 empurra T para cima a cada
/// iteração até o clamp — errata completa em r-33 §2). Com p = softmax(z/T):
/// `L(T)  = média_i [−z_{i,y_i}/T + log Σ_j e^{z_{i,j}/T}]`
/// `L'(T) = média_i ( z_{i,y_i} − E_{p_i}[z_i] ) / T²`
/// `L''(T)= média_i ( Var_{p_i}(z_i) − 2T·(z_{i,y_i} − E_{p_i}[z_i]) ) / T⁴`
/// (Teo N-1 — derivação direta; no ótimo L'(T*)=0 ⇒ L''(T*) = Var/T⁴ > 0.)
fn grad_hess_temperatura(logits: &Tensor, alvos: &Tensor, t: f64) -> Result<(f64, f64, f64)> {
    let n = logits.dim(0)?;
    let escalado = logits.affine(1.0 / t, 0.0)?;
    let ce = candle_nn::loss::cross_entropy(&escalado, alvos)?;
    let l = ce.to_scalar::<f32>()? as f64;
    let p = candle_nn::ops::softmax(&escalado, candle_core::D::Minus1)?;
    // E_p[z] e E_p[z²] por linha — Var_p(z) = E[z²] − (E[z])² em UMA passada
    let ez = (&p * logits)?.sum_keepdim(candle_core::D::Minus1)?;
    let ez2 = (&p * &logits.sqr()?)?.sum_keepdim(candle_core::D::Minus1)?;
    let var = (&ez2 - &(&ez * &ez)?)?;
    // z_y por linha via gather
    let idx = alvos.reshape((n, 1))?;
    let zt = logits.gather(&idx, candle_core::D::Minus1)?;
    let ez_m = ez.mean_all()?.to_scalar::<f32>()? as f64;
    let var_m = var.mean_all()?.to_scalar::<f32>()? as f64;
    let zt_m = zt.mean_all()?.to_scalar::<f32>()? as f64;
    let dif = zt_m - ez_m;
    let t2 = t * t;
    let dl = dif / t2;
    let d2l = (var_m - 2.0 * t * dif) / (t2 * t2);
    Ok((l, dl, d2l))
}

/// **N-4:** calibra T* por **Newton-Raphson com derivadas analíticas**
/// (convergência quadrática — Teo N-3) e **salvaguarda de bisseção
/// geométrica** sobre o sinal de L̃'. A iteração é feita em **γ = 1/T**:
/// L̃(γ) = NLL(softmax(γz)) é CONVEXA (N-2: L̃''(γ) = Var_p(z) ≥ 0
/// SEMPRE — não apenas perto do ótimo), com
/// `L̃'(γ) = E_p[z] − z_y` e `L̃''(γ) = Var_p(z)` — o Newton em γ nunca
/// encontra curvatura negativa e o colchão γ ∈ [1/8, 8] (⇔ T ∈ [1/8, 8],
/// o mesmo da V-1) é mantido como invariante com a bisseção de salvaguarda.
/// Pior caso: bisseção pura, ⌈log₂(ln 64/ε)⌉ ≈ 35 passos; típico: 4–12
/// avaliações (vs ~44 do golden-section). Devolve (T*, nll(1), nll(T*), n_aval).
pub fn calibrar_temperatura_newton(
    logits: &Tensor,
    alvos: &Tensor,
    max_iters: usize,
) -> Result<(f64, f64, f64, usize)> {
    let nll_em = |t: f64| nll_temperatura(logits, alvos, t);
    let nll_1 = nll_em(1.0)?;
    // (L, L̃'(γ), L̃''(γ)) com escalado = γ·z (p = softmax(γz)):
    // L̃' = média(E_p[z] − z_y); L̃'' = média(Var_p(z)) ≥ 0 (N-2)
    let gh_gamma = |g: f64| -> Result<(f64, f64, f64)> {
        let n = logits.dim(0)?;
        let escalado = logits.affine(g, 0.0)?;
        let ce = candle_nn::loss::cross_entropy(&escalado, alvos)?;
        let l = ce.to_scalar::<f32>()? as f64;
        let p = candle_nn::ops::softmax(&escalado, candle_core::D::Minus1)?;
        let ez = (&p * logits)?.sum_keepdim(candle_core::D::Minus1)?;
        let ez2 = (&p * &logits.sqr()?)?.sum_keepdim(candle_core::D::Minus1)?;
        let var = (&ez2 - &(&ez * &ez)?)?;
        let idx = alvos.reshape((n, 1))?;
        let zt = logits.gather(&idx, candle_core::D::Minus1)?;
        let ez_m = ez.mean_all()?.to_scalar::<f32>()? as f64;
        let var_m = var.mean_all()?.to_scalar::<f32>()? as f64;
        let zt_m = zt.mean_all()?.to_scalar::<f32>()? as f64;
        Ok((l, ez_m - zt_m, var_m))
    };
    let (lo, hi) = (1.0f64 / 8.0, 8.0f64); // colchão em γ ⇔ T ∈ [1/8, 8]
    // mínimo na borda? (convexa: L̃' não troca sinal dentro do colchão)
    let (_, dlo, _) = gh_gamma(lo)?;
    if dlo >= 0.0 {
        // mínimo em γ=lo ⇒ T* = 8 (esfriar ao máximo)
        return Ok((1.0 / lo, nll_1, nll_em(1.0 / lo)?, 1));
    }
    let (_, dhi, _) = gh_gamma(hi)?;
    if dhi <= 0.0 {
        // mínimo em γ=hi ⇒ T* = 1/8 (aquecer ao máximo)
        return Ok((1.0 / hi, nll_1, nll_em(1.0 / hi)?, 1));
    }
    // invariante do colchão: L̃'(a) < 0 < L̃'(b) (L̃ convexa ⇒ válida)
    let (mut a, mut b) = (lo, hi);
    let mut g = 1.0f64; // chute inicial γ=1 (T=1, incalibrado)
    let mut melhor = (1.0f64, nll_1); // (T, NLL)
    let mut n_aval = 2usize;
    for _ in 0..max_iters {
        let (l, dp, d2p) = gh_gamma(g)?;
        n_aval += 1;
        let t = 1.0 / g;
        if l < melhor.1 {
            melhor = (t, l);
        }
        if dp.abs() < 1e-10 || (b - a) < 1e-12 {
            break;
        }
        // atualiza o colchão pelo sinal de L̃' (bisseção válida — N-2)
        if dp < 0.0 {
            a = a.max(g);
        } else {
            b = b.min(g);
        }
        // Newton em γ: L̃'' = Var ≥ 0 SEMPRE (convexidade global);
        // passo aceito só DENTRO do colchão, senão bisseção geométrica
        let gn = if d2p > 1e-12 { g - dp / d2p } else { f64::NAN };
        g = if gn.is_finite() && gn > a && gn < b { gn } else { (a * b).sqrt() };
    }
    Ok((melhor.0, nll_1, melhor.1, n_aval))
}

/// ECE (Expected Calibration Error) por faixas de confiança.
pub fn ece(confs: &[f64], corretos: &[bool], n_faixas: usize) -> f64 {
    if confs.is_empty() || confs.len() != corretos.len() {
        return f64::NAN;
    }
    let n = confs.len() as f64;
    let mut total = 0f64;
    for b in 0..n_faixas {
        let (lo, hi) = (b as f64 / n_faixas as f64, (b + 1) as f64 / n_faixas as f64);
        let mut soma_conf = 0f64;
        let mut soma_acerto = 0f64;
        let mut m = 0usize;
        for (c, ok) in confs.iter().zip(corretos.iter()) {
            let c1 = c.min(0.999_999);
            if *c >= lo && c1 < hi {
                soma_conf += c;
                soma_acerto += if *ok { 1.0 } else { 0.0 };
                m += 1;
            }
        }
        if m > 0 {
            total += (m as f64 / n) * (soma_conf / m as f64 - soma_acerto / m as f64).abs();
        }
    }
    total
}

/// **V-3:** limiar τ_c = quantil 1−ε das confianças retidas (ordem estatística).
pub fn limiar_por_erro(confs: &[f64], epsilon: f64) -> f64 {
    let mut v = confs.to_vec();
    if v.is_empty() {
        return 1.0;
    }
    v.sort_by(|a, b| a.partial_cmp(b).unwrap());
    let idx = (((1.0 - epsilon.clamp(0.0, 1.0)) * v.len() as f64).ceil() as usize)
        .saturating_sub(1)
        .min(v.len() - 1);
    v[idx]
}

/// Verificador NCM calibrado — o substituto de produção do oráculo.
#[derive(Debug, Clone)]
pub struct VerificadorNcm {
    /// temperatura calibrada T* (1 = incalibrado).
    pub temperatura: f64,
    pub ece_treino: f64,
    pub n_ll_antes: f64,
    pub n_ll_depois: f64,
}

impl VerificadorNcm {
    pub fn novo() -> Self {
        VerificadorNcm { temperatura: 1.0, ece_treino: 0.0, n_ll_antes: 0.0, n_ll_depois: 0.0 }
    }

    /// Ajusta T* em (logits, alvos) retidos e mede o ECE antes/depois.
    /// v0.8.0: Newton-Raphson analítico (N-4) com salvaguarda — o
    /// golden-section (V-1) permanece como verificação/rede de segurança.
    pub fn ajustar(&mut self, logits: &Tensor, alvos: &Tensor) -> Result<(f64, f64)> {
        let (mut t, mut nll1, mut nllt) = (1.0, 0.0, 0.0);
        let mut ok = false;
        // Newton com salvaguarda (4× o orçamento típico de avaliações)
        if let Ok((tn, n1, nt, _)) = calibrar_temperatura_newton(logits, alvos, 40) {
            if nt.is_finite() && nt <= n1 + 1e-6 {
                (t, nll1, nllt) = (tn, n1, nt);
                ok = true;
            }
        }
        if !ok {
            // rede de segurança: golden-section (V-1) — nunca removida
            let (tg, n1g, ntg) = calibrar_temperatura(logits, alvos)?;
            (t, nll1, nllt) = (tg, n1g, ntg);
        }
        self.temperatura = t;
        self.n_ll_antes = nll1;
        self.n_ll_depois = nllt;
        let preds = logits.argmax(candle_core::D::Minus1)?;
        let corretos_v = preds
            .eq(alvos)?
            .to_dtype(candle_core::DType::U32)?
            .flatten_all()?
            .to_vec1::<u32>()?;
        let confs_antes: Vec<f64> = confianca_por_linha(logits)?
            .flatten_all()?.to_vec1::<f32>()?.iter().map(|x| *x as f64).collect();
        let corretos: Vec<bool> = corretos_v.iter().map(|x| *x == 1).collect();
        let e_antes = ece(&confs_antes, &corretos, 10);
        let escalado = logits.affine(1.0 / t, 0.0)?;
        let confs_depois: Vec<f64> = confianca_por_linha(&escalado)?
            .flatten_all()?.to_vec1::<f32>()?.iter().map(|x| *x as f64).collect();
        let e_depois = ece(&confs_depois, &corretos, 10);
        self.ece_treino = e_depois;
        Ok((e_antes, e_depois))
    }

    /// Confiança calibrada (n, K) → (n,) — a decisão do ciclo d usa isto.
    pub fn confianca(&self, logits: &Tensor) -> Result<Tensor> {
        let escalado = logits.affine(1.0 / self.temperatura, 0.0)?;
        confianca_por_linha(&escalado)
    }

    /// Confiança calibrada de UMA linha (K,) → f64.
    pub fn confianca_1d(&self, logits: &Tensor) -> Result<f64> {
        let c = self.confianca(&logits.reshape((1, logits.elem_count()))?)?;
        Ok(c.flatten_all()?.to_vec1::<f32>()?[0] as f64)
    }
}

impl Default for VerificadorNcm {
    fn default() -> Self {
        Self::novo()
    }
}

pub fn versao() -> &'static str {
    "verificador-0.6.0"
}

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

    fn logits_sinteticos(n: usize, k: usize, escala: f32, semente: u64) -> (Tensor, Tensor) {
        let mut rng = crate::rng::RngDeterministico::novo(semente);
        let mut l = Vec::with_capacity(n * k);
        let mut a = Vec::with_capacity(n);
        for _ in 0..n {
            // classe verdadeira com margem; logits = one-hot(alvo)·margem + ruído
            let alvo = rng.inteiro(k) as u32;
            a.push(alvo);
            for j in 0..k {
                let base = if j as u32 == alvo { 3.0 } else { 0.0 };
                l.push(base + rng.normal(0.0, 1.0) as f32 * escala);
            }
        }
        let dev = Device::Cpu;
        let lt = Tensor::from_vec(l, (n, k), &dev).unwrap();
        let at = Tensor::from_vec(a, (n,), &dev).unwrap();
        (lt, at)
    }

    #[test]
    fn v1_calibracao_reduz_ece_e_nll() {
        // superconfiante (escala 3×): T* > 1 deve reduzir ECE
        let (lt, at) = logits_sinteticos(512, 16, 3.0, 7);
        let mut v = VerificadorNcm::novo();
        let (e_antes, e_depois) = v.ajustar(&lt, &at).unwrap();
        assert!(v.temperatura > 1.0, "superconfiante exige T*>1 (T*={})", v.temperatura);
        assert!(e_depois < e_antes, "ECE deve cair: {e_antes:.4} -> {e_depois:.4}");
        assert!(v.n_ll_depois <= v.n_ll_antes + 1e-6, "NLL calibrada ≤ NLL bruta");
    }

    #[test]
    fn v2_soundness_sem_oraculo() {
        // P(erro | aceito) ≤ 1 − τ_c + ECE + folga estatística
        let (lt, at) = logits_sinteticos(4096, 8, 2.0, 11);
        let mut v = VerificadorNcm::novo();
        let (_, e_depois) = v.ajustar(&lt, &at).unwrap();
        let confs: Vec<f64> = v.confianca(&lt).unwrap().flatten_all().unwrap()
            .to_vec1::<f32>().unwrap().iter().map(|x| *x as f64).collect();
        let preds = lt.argmax(candle_core::D::Minus1).unwrap();
        let corretos: Vec<bool> = preds.eq(&at).unwrap().to_dtype(candle_core::DType::U32).unwrap()
            .flatten_all().unwrap().to_vec1::<u32>().unwrap().iter().map(|x| *x == 1).collect();
        let tau = limiar_por_erro(&confs, 0.20); // aceitar ~80%
        let mut erro_aceito = 0usize;
        let mut n_aceito = 0usize;
        for (c, ok) in confs.iter().zip(corretos.iter()) {
            if *c >= tau {
                n_aceito += 1;
                if !ok {
                    erro_aceito += 1;
                }
            }
        }
        let p_erro = erro_aceito as f64 / n_aceito.max(1) as f64;
        let cota = (1.0 - tau) + e_depois + 0.05; // folga de amostragem finita
        assert!(p_erro <= cota, "P(erro|aceito)={p_erro:.4} ≤ 1−τ_c+ECE+folga={cota:.4} (τ_c={tau:.3}, ECE={e_depois:.4})");
    }

    #[test]
    fn v3_unimodalidade_da_nll_em_t() {
        // verificação empírica da unimodalidade (condição do golden-section)
        let (lt, at) = logits_sinteticos(1024, 32, 2.5, 13);
        let mut prev = f64::INFINITY;
        let mut descendo = true;
        let mut subindo = false;
        let ts: Vec<f64> = (0..40).map(|i| 0.25 * (i as f64 / 39.0 * 6.0).exp()).collect();
        for t in &ts {
            let n = nll_temperatura(&lt, &at, *t).unwrap();
            if n > prev {
                descendo = false;
                subindo = true;
            } else if subindo && n < prev {
                panic!("NLL não unimodal em T (re-queda após subida)");
            }
            prev = n;
        }
        assert!(!descendo || true);
    }

    #[test]
    fn v4_causalidade_somente_presente() {
        // V é função apenas dos logits correntes: mesmos logits ⇒ mesma conf.
        let (lt, _) = logits_sinteticos(64, 8, 1.0, 17);
        let v = VerificadorNcm { temperatura: 1.37, ..Default::default() };
        let c1 = v.confianca(&lt).unwrap();
        let c2 = v.confianca(&lt).unwrap();
        let d = (&c1 - &c2).unwrap().sqr().unwrap().sum_all().unwrap()
            .to_scalar::<f32>().unwrap();
        assert!(d == 0.0, "determinístico nos logits do presente");
    }

    /// **N-1:** gradiente/curvatura analíticos batem com diferenças finitas
    /// centralmente — validação numérica da derivação do Teo N-1 (e errata
    /// da proposta Qwen-3.8 nº 1: sinal de L' e ordem da diferença finita).
    #[test]
    fn n1_grad_hess_analiticos_vs_finitas() {
        let (lt, at) = logits_sinteticos(512, 16, 2.0, 23);
        let t = 1.3f64;
        let (_, dl, d2l) = grad_hess_temperatura(&lt, &at, t).unwrap();
        let h = 1e-2; // f32: equilíbrio ruído/truncamento (r-33 §2.4)
        let l_mais = nll_temperatura(&lt, &at, t + h).unwrap();
        let l_zero = nll_temperatura(&lt, &at, t).unwrap();
        let l_menos = nll_temperatura(&lt, &at, t - h).unwrap();
        let dl_fd = (l_mais - l_menos) / (2.0 * h);
        // SEGUNDA diferença (a proposta original usava a PRIMEIRA)
        let d2l_fd = (l_mais - 2.0 * l_zero + l_menos) / (h * h);
        assert!(
            (dl - dl_fd).abs() < 5e-3,
            "N-1: L' analítico {dl:.6} vs finita {dl_fd:.6}"
        );
        assert!(
            (d2l - d2l_fd).abs() < 5e-2,
            "N-1: L'' analítico {d2l:.6} vs finita {d2l_fd:.6}"
        );
        // REGRESSÃO DA ERRATA: superconfiante (escala 2) ⇒ T* > 1 ⇒ L'(T=1)
        // deve ser NEGATIVO (L desce em direção a T*). O sinal invertido da
        // proposta daria positivo — e o Newton divergiria para o clamp 10.
        let (_, dl_em_1, _) = grad_hess_temperatura(&lt, &at, 1.0).unwrap();
        assert!(dl_em_1 < 0.0, "N-1: L'(1)={dl_em_1} deve ser < 0 (superconfiante)");
    }

    /// **N-2/N-4:** Newton com salvaguarda encontra o MESMO T* do
    /// golden-section (V-1) com MENOS avaliações de NLL.
    #[test]
    fn n2_newton_igual_golden_mais_rapido() {
        for (escala, semente) in [(1.0f32, 31u64), (3.0, 37), (0.3, 41)] {
            let (lt, at) = logits_sinteticos(1024, 32, escala, semente);
            let (t_g, _, nll_g) = calibrar_temperatura(&lt, &at).unwrap();
            let (t_n, nll_1, nll_n, n_aval) =
                calibrar_temperatura_newton(&lt, &at, 40).unwrap();
            assert!(
                (t_g - t_n).abs() < 1e-3,
                "N-2: T* golden {t_g:.6} vs newton {t_n:.6} (escala {escala})"
            );
            assert!(
                nll_n <= nll_g + 1e-4,
                "N-4: NLL newton {nll_n:.6} ≤ golden {nll_g:.6} + 1e-4"
            );
            assert!(
                nll_n <= nll_1 + 1e-6,
                "N-4: NLL calibrada {nll_n:.6} ≤ bruta {nll_1:.6}"
            );
            // golden ~44 avaliações típicas (66 no teto); Newton ≤ 36
            assert!(
                n_aval < 36,
                "N-4: newton usou {n_aval} avaliações (esperado < 36)"
            );
        }
    }

    /// **N-2 (verificação do teorema):** L' troca o sinal EXATAMENTE uma
    /// vez no colchão (unimodalidade — NLL convexa em γ = 1/T).
    #[test]
    fn n3_sinal_de_l_troca_uma_vez() {
        let (lt, at) = logits_sinteticos(2048, 16, 2.5, 43);
        let mut trocas = 0usize;
        let mut sinal_ant: Option<f64> = None;
        for i in 0..64 {
            let t = ((i as f64 / 63.0 * (8.0f64.ln() - (1.0f64 / 8.0).ln()))
                + (1.0f64 / 8.0).ln())
                .exp();
            let (_, dl, _) = grad_hess_temperatura(&lt, &at, t).unwrap();
            let s: f64 = if dl < 0.0 { -1.0 } else { 1.0 };
            if let Some(sa) = sinal_ant {
                if s != sa {
                    trocas += 1;
                }
            }
            sinal_ant = Some(s);
        }
        assert_eq!(trocas, 1, "N-2: L' deve trocar o sinal exatamente 1 vez (trocas={trocas})");
    }

    /// **N-4 (casos extremos):** direção correta da calibração nos dois
    /// regimes — superconfiante (ERROS confiantes: margem alta com ruído
    /// moderado ⇒ L̃'(1) = média(E_p[z] − z_y) > 0 ⇒ γ* < 1 ⇒ T* > 1,
    /// esfria) e subconfiante (margem modesta com ruído mínimo: corretos
    /// dominam com p_y < 1 ⇒ T* < 1, aquece). Nota: a contribuição de cada
    /// exemplo para L̃' é Σ_{j≠y} p_j(z_j − z_y) ≤ 0 para acertos e ≫ 0
    /// para erros confiantes — é por isso que a NLL "escuta" os erros.
    #[test]
    fn n4b_casos_extremos() {
        let dev = Device::Cpu;
        // (a) superconfiante (padrão de overfitting): 95% acertos com
        // margem 5,4 (p_y ≈ 0,97) e 5% ERROS confiantes com margem 8,4
        // (p_errado ≈ 1) ⇒ conf ≈ 0,97 > acc = 0,95 ⇒ T* > 1. A
        // contribuição dos erros para L̃'(1) é ≈ +8,4·0,05 = +0,42 vs
        // −0,15 dos acertos ⇒ L̃'(1) > 0 ⇒ γ* < 1 ⇒ T* > 1 (esfria).
        let mut rng = crate::rng::RngDeterministico::novo(47);
        let (n, k) = (2048usize, 8usize);
        let mut l = Vec::with_capacity(n * k);
        let mut a = Vec::with_capacity(n);
        for _ in 0..n {
            let alvo = rng.inteiro(k) as u32;
            a.push(alvo);
            if rng.uniforme() < 0.05 {
                let mut w = rng.inteiro(k) as u32;
                if w == alvo {
                    w = (w + 1) % k as u32;
                }
                for j in 0..k {
                    l.push(if j as u32 == w { 8.4 } else { 0.0 });
                }
            } else {
                for j in 0..k {
                    let base = if j as u32 == alvo { 5.4 } else { 0.0 };
                    l.push(base + rng.normal(0.0, 0.25) as f32);
                }
            }
        }
        let lt = Tensor::from_vec(l, (n, k), &dev).unwrap();
        let at = Tensor::from_vec(a, (n,), &dev).unwrap();
        let (t_n, nll_1, nll_n, _) = calibrar_temperatura_newton(&lt, &at, 40).unwrap();
        assert!(t_n > 1.0, "N-4b: superconfiante exige T*>1 (T*={t_n})");
        assert!(nll_n <= nll_1 + 1e-6);
        assert!((1.0 / 8.0..=8.0).contains(&t_n), "N-4b: T* dentro do colchão");
        // (b) subconfiante: margem 3 com ruído mínimo ⇒ p_y ≈ 0,74 < acc
        // ≈ 1 ⇒ T* < 1 (aquece; medido na borda quente 1/8)
        let (lt2, at2) = logits_sinteticos(512, 8, 0.1, 53);
        let (t_n2, nll_12, nll_n2, _) = calibrar_temperatura_newton(&lt2, &at2, 40).unwrap();
        assert!(t_n2 < 1.0, "N-4b: subconfiante exige T*<1 (T*={t_n2})");
        assert!(nll_n2 <= nll_12 + 1e-6, "N-4b: NLL não piora");
        assert!((1.0 / 8.0..=8.0).contains(&t_n2));
    }
}