#!/usr/bin/env python3
"""
A MATEMÁTICA DA STACK — X-FACTOR, NANO E NORMALIZE
==================================================

Três peças, três estruturas matemáticas diferentes. Este benchmark mede cada
uma contra a sua própria definição, não contra a intuição de quem a escreveu.

  N1  nano: o isomorfismo ℤ/256ℤ → μ₂₅₆
      Encode: f(b) = e^{i(b+pos+j)·2π/256}
      Decode: b = arg(ψ)·256/2π − pos − j  (mod 256)
      Testado: a homomorfia f(a+b) = f(a)f(b), a injetividade, o round-trip
      sobre os 256 bytes, e quanta informação os 64 modos × 4 componentes
      realmente carregam.

  X1  xfactor: power iteration com deflação contra `eigh` exato
      O rank efetivo é correto? O subespaço é o mesmo? O erro cresce com r?

  X2  xfactor: por que P não é idempotente
      A deflação acumula erro de ortogonalidade em float32. Mede a
      assinatura: o desvio deve crescer com o número de autovetores.

  G1  normalize: o roteador acerta?
      Roteia por cobertura de prova. Testa manifold forte, fraco e ruído
      isotrópico, e compara o recall roteado contra bases fixas.

  G2  normalize: de onde vem o ganho do roteador
      O relatório anuncia basis=pca_corpus, k1=0.05. Isola se o ganho vem
      disso ou de outro parâmetro.

Uso:
    python3 matematica_da_stack.py
    python3 matematica_da_stack.py --rapido

BSL 1.1 | (c) Winnex Brasil Soluções Empresariais LTDA-ME
"""

from __future__ import annotations

import argparse
import json
import os
import sys
import time

os.environ.setdefault("OMP_NUM_THREADS", "4")

import numpy as np

RESULTS: dict[str, object] = {}


def titulo(t: str) -> None:
    print(f"\n{'=' * 74}\n {t}\n{'=' * 74}", flush=True)


def teste(nome: str, afirma: str, falsifica: str) -> None:
    print(f"\n--- {nome} ---", flush=True)
    print(f"    afirmação : {afirma}", flush=True)
    print(f"    falsifica : {falsifica}", flush=True)
    print("    " + "-" * 60, flush=True)


# ===========================================================================
# N1 — o isomorfismo do tokenizador espectral
# ===========================================================================

def n1_tokenizador(wn) -> dict:
    titulo("N1 — nano: o isomorfismo ℤ/256ℤ → μ₂₅₆")
    out: dict[str, object] = {}
    tk = wn.SpectralTokenizer(embed_dim=64)

    teste("N1.1 Homomorfia e injetividade",
          "f(a+b) = f(a)·f(b) sobre fases, e a → f(a) é injetora em ℤ/256ℤ",
          "dois bytes distintos produzirem a mesma fase")
    # O tokenizador é orientado a BYTES UTF-8 (não a caracteres). A fase
    # codifica (byte + pos) mod 256. Testar com bytes crus, não com str
    # latin-1 (que seria re-codificado em UTF-8 e mudaria os bytes).
    import ctypes

    def fases_de(bs: bytes):
        """Fase de cada byte de uma sequência, via encode() do binding."""
        txt = bs.decode("utf-8", errors="surrogateescape")
        q = tk.encode(txt)
        n = len(q) // 64
        out = []
        for i in range(n):
            s = q[i * 64]
            out.append((np.arctan2(s.x, s.w) % (2 * np.pi)) * 256 / (2 * np.pi))
        return out

    amostra = bytes(range(1, 128))          # ASCII imprimível + controle
    f = fases_de(amostra)
    distintos = {int(round(v)) % 256 for v in f}
    out["bytes_unicos_ascii"] = len(distintos)
    print(f"    {len(amostra)} bytes ASCII → {len(distintos)} fases distintas")

    # A fase observada é (byte + pos) mod 256; o isomorfismo é sobre o BYTE,
    # então a posição precisa ser descontada antes de testar a aditividade.
    byte_rec = [(f[i] - i) % 256 for i in range(len(f))]
    a, b_ = 37, 60
    fa, fb = byte_rec[a - 1], byte_rec[b_ - 1]
    fab = byte_rec[(a + b_) - 1]
    err_hom = abs((fa + fb - fab + 128) % 256 - 128)
    out["erro_homomorfia"] = float(err_hom)
    print(f"    byte recuperado = fase − pos; amostra: {byte_rec[:5]}")
    print(f"    homomorfia byte(a)+byte(b) vs byte(a+b): erro {err_hom:.4f} "
          f"em unidades de 1/256 de volta")
    ok = len(distintos) == len(amostra) and err_hom < 0.5
    if ok:
        print("    ✓ injetora em [1,127] e aditiva em fase — isomorfismo de "
              "grupo sobre ℤ/256ℤ (o domínio completo exige bytes ≥128, que "
              "não existem isolados em UTF-8 válido)")
    else:
        print(f"    ✗ bijeção={len(distintos) == len(amostra)} "
              f"homomorfia={err_hom < 0.5}")

    teste("N1.2 Round-trip sobre todo o alfabeto",
          "decode(encode(c)) == c para cada um dos 256 valores de byte",
          "um único byte falhar")
    ok = 0
    falhas = []
    for b in range(256):
        if b == 0:
            continue          # NUL não sobrevive a uma str Python
        try:
            ch = bytes([b]).decode("latin-1")
            if tk.decode(tk.encode(ch)) == ch:
                ok += 1
            else:
                falhas.append(b)
        except Exception:  # noqa: BLE001
            falhas.append(b)
    out["roundtrip_bytes"] = f"{ok}/255"
    print(f"    {ok}/255 bytes fazem round-trip exato")
    if falhas:
        print(f"    falhas em: {falhas[:8]}{'...' if len(falhas) > 8 else ''}")
    (print("    ✓ round-trip completo" if ok == 255 else "    ✗ incompleto"))

    teste("N1.3 Quanta informação os 64 modos carregam?",
          "cada um dos 64 modos × 4 componentes contribui informação",
          "os componentes forem rotações determinísticas uns dos outros — "
          "nesse caso a representação tem 8 bits por caractere guardados em "
          "1024× o necessário")
    txt = "A superposição é uma característica"
    q = np.array([[x.w, x.x, x.y, x.z] for x in tk.encode(txt)])
    nb = len(q) // 64
    q = q[: nb * 64].reshape(nb, 64, 4)
    w0, x0 = q[:, 0, 0], q[:, 0, 1]

    # (y,z) deveria ser (w,x) rotacionado em π/4
    c45, s45 = np.cos(np.pi / 4), np.sin(np.pi / 4)
    e_y = float(np.abs(q[:, :, 2] - (q[:, :, 0] * c45 - q[:, :, 1] * s45)).max())
    e_z = float(np.abs(q[:, :, 3] - (q[:, :, 0] * s45 + q[:, :, 1] * c45)).max())
    # modos j>0 deveriam ser o modo 0 rotacionado em j·2π/256
    e_j = []
    for j in (1, 2, 7, 31, 63):
        a = j * 2 * np.pi / 256
        e_j.append(float(np.abs(
            q[:, j, 0] - (w0 * np.cos(a) - x0 * np.sin(a))).max()))

    out.update(erro_yz=max(e_y, e_z), erro_modos=max(e_j))
    print(f"    ‖(y,z) − rot_{'{π/4}'}(w,x)‖max = {max(e_y, e_z):.2e}")
    print(f"    ‖modo_j − rot_{{j·2π/256}}(modo_0)‖max = {max(e_j):.2e}")
    print(f"    → 1 byte (8 bits) guardado em 64×4 floats ({64*4*32} bits): "
          f"expansão de {64*4*32 // 8}×")
    if max(e_y, e_z) < 1e-5 and max(e_j) < 1e-5:
        print("    ✓ confirmado: (y,z) e os modos j>0 são rotações "
              "determinísticas do modo 0 — a informação é 1 fase por byte")
    else:
        print("    ~ há conteúdo independente nos modos/componentes extras")

    RESULTS["nano"] = out
    return out


# ===========================================================================
# X1 / X2 — power iteration e a não-idempotência
# ===========================================================================

def x1_power_iteration(wn) -> dict:
    titulo("X1 — xfactor: power iteration vs eigendecomposição exata")
    out: dict[str, object] = {}

    teste("X1.1 O rank efetivo é correto?",
          "r(τ) do power iteration coincide com o de eigh exato",
          "o rank divergir — indicaria que o critério de parada erra a "
          "variância acumulada")
    rng = np.random.default_rng(0)
    linhas = []
    for V, D, tau in ((2000, 64, 0.95), (5000, 128, 0.95), (3000, 32, 0.90)):
        E = rng.standard_normal((V, D)).astype(np.float32)
        P, rank, var = wn.compute_xfactor(E, V, D, tau)
        Ec = E - E.mean(0)
        C = (Ec.T @ Ec) / V
        w = np.clip(np.linalg.eigvalsh(C)[::-1], 0, None)
        r_ex = int(np.searchsorted(np.cumsum(w) / w.sum(), tau) + 1)
        U = np.linalg.eigh(C)[1][:, ::-1][:, :r_ex]
        err = float(np.abs(np.asarray(P) - U @ U.T).max())
        linhas.append({"V": V, "D": D, "rank": int(rank), "rank_exato": r_ex,
                       "var": round(float(var), 4), "err_subespaco": err})
        print(f"    V={V:<5} D={D:<4} τ={tau}: rank {rank} (exato {r_ex}) | "
              f"var {var:.4f} | ‖P−P_exato‖={err:.2e}")
    out["casos"] = linhas
    ok = all(l["rank"] == l["rank_exato"] for l in linhas)
    RESULTS["xfactor_rank"] = out
    (print("    ✓ rank efetivo exato em todos os casos" if ok
           else "    ✗ rank divergiu"))

    teste("X2 A não-idempotência tem a assinatura da deflação?",
          "o desvio de ortogonalidade cresce com o número de autovetores "
          "deflacionados",
          "o desvio for constante no rank — indicaria erro de transcrição, "
          "não acúmulo numérico")
    D, V = 96, 4000
    E = rng.standard_normal((V, D)).astype(np.float32)
    tab = []
    print(f"    {'τ':>6}{'rank':>6}{'var_motor':>12}{'‖P²−P‖':>12}")
    for tau in (0.3, 0.5, 0.7, 0.85, 0.95, 0.99):
        P, rank, var = wn.compute_xfactor(E, V, D, tau)
        P = np.asarray(P)
        idem = float(np.abs(P @ P - P).max())
        tab.append({"tau": tau, "rank": int(rank), "var": round(float(var), 4),
                    "idem": idem})
        print(f"    {tau:>6}{rank:>6}{var:>12.4f}{idem:>12.2e}")
    out["idempotencia"] = tab
    ranks = [t["rank"] for t in tab]
    idems = [t["idem"] for t in tab]
    cresc = idems[-1] > idems[0]
    print(f"\n    rank {min(ranks)}→{max(ranks)}: ‖P²−P‖ {idems[0]:.2e}→{idems[-1]:.2e}")
    if cresc:
        print("    ✓ o desvio cresce com o número de deflações — assinatura de "
              "acúmulo numérico em float32, não de erro de fórmula")
    else:
        print("    ~ o desvio não cresce monotonicamente; a causa pode ser outra")
    return out


# ===========================================================================
# G1 / G2 — o roteador do normalize
# ===========================================================================

def _manifold(rng, n, d, r, ruido=0.05):
    B = rng.standard_normal((r, d))
    Cm = rng.standard_normal((n, r))
    X = Cm @ B
    X += rng.standard_normal((n, d)) * ruido * np.linalg.norm(X) / np.sqrt(n * d)
    return (X / np.linalg.norm(X, axis=1, keepdims=True)).astype(np.float32)


def g1_g2_normalize(rapido: bool) -> dict:
    titulo("G1/G2 — winnex-ai-normalize: o roteador e de onde vem o ganho")
    try:
        from winnex_ai_normalize.core.quality import (
            audit_corpus, build_quality_engine)
        from winnex_madhava import build_engine
    except ImportError as exc:
        print(f"    (dependência ausente: {exc})", flush=True)
        return {"disponivel": False}

    rng = np.random.default_rng(0)
    n, d = (1500, 96) if rapido else (3000, 128)
    casos = {
        "manifold r=32": _manifold(rng, n, d, 32),
        "manifold r=96": _manifold(rng, n, d, max(2, d - 32)),
        "isotrópico": (lambda A: A / np.linalg.norm(A, axis=1, keepdims=True))(
            rng.standard_normal((n, d))).astype(np.float32),
    }

    teste("G1 O roteador escolhe a base certa?",
          "o gate distingue corpus com manifold de ruído isotrópico e roteia "
          "cada um para a configuração que maximiza o recall",
          "rotas idênticas para os três, ou recall roteado abaixo da melhor "
          "base fixa")

    out: dict[str, object] = {}
    for nome, X in casos.items():
        q = X[:40]
        gt = np.argsort(-(q @ X.T), axis=1)[:, :10]
        rep = audit_corpus(X, dim=d, n_seed_queries=6)

        def rec(eng):
            p = []
            for x in q:
                r_ = eng.search(x)
                p.append(np.asarray(getattr(r_, "indices", r_))[:10])
            return float(np.mean([len(set(a[:10]) & set(b[:10])) / 10
                                  for a, b in zip(p, gt)]))

        eng_rot, _ = build_quality_engine(X, dim=d, return_report=True)
        r_rot = rec(eng_rot)
        r_rnd = rec(build_engine(X, dim=d, metric="cosine", k=10, basis="random"))
        r_pca = rec(build_engine(X, dim=d, metric="cosine", k=10,
                                 basis="pca_corpus"))
        melhor = max(r_rnd, r_pca)
        out[nome] = {"basis": str(rep.basis), "k1": float(rep.k1_fraction),
                     "flags": [(f.code, f.severity) for f in rep.flags],
                     "recall_roteado": round(r_rot, 4),
                     "recall_random": round(r_rnd, 4),
                     "recall_pca": round(r_pca, 4)}
        print(f"\n    {nome}")
        print(f"      rotas: basis={rep.basis} k1={rep.k1_fraction} "
              f"flags={[(f.code, f.severity) for f in rep.flags]}")
        print(f"      recall: roteado {r_rot:.4f} | random {r_rnd:.4f} | "
              f"pca_corpus {r_pca:.4f}")
        if r_rot >= melhor - 1e-9:
            print("      ✓ o roteado iguala ou supera a melhor base fixa")
        else:
            print(f"      ✗ o roteado fica {melhor - r_rot:.4f} abaixo da "
                  f"melhor base fixa ({'pca' if r_pca > r_rnd else 'random'})")

    # ------------------------------------------------------------------
    teste("G2 De onde vem o ganho do roteador?",
          "o relatório atribui o ganho a basis=pca_corpus e k1=0.05, e esses "
          "parâmetros de fato explicam o recall obtido",
          "reproduzir o build direto com os parâmetros anunciados der recall "
          "diferente do roteado — indicaria que outro parâmetro é o agente")

    X = casos["isotrópico"]
    q = X[:40]
    gt = np.argsort(-(q @ X.T), axis=1)[:, :10]

    def rec2(eng):
        p = []
        for x in q:
            r_ = eng.search(x)
            p.append(np.asarray(getattr(r_, "indices", r_))[:10])
        return float(np.mean([len(set(a[:10]) & set(b[:10])) / 10
                              for a, b in zip(p, gt)]))

    eng_rot, rep = build_quality_engine(X, dim=d, return_report=True)
    r_rot = rec2(eng_rot)
    print(f"    relatório anuncia : basis={rep.basis} k1={rep.k1_fraction} "
          f"stage1_dim={rep.stage1_dim}")
    print(f"    recall roteado    : {r_rot:.4f}")

    variantes = {
        "basis+k1 anunciados": dict(basis=str(rep.basis),
                                    k1_fraction=float(rep.k1_fraction)),
        "+ stage1_dim do gate": dict(basis=str(rep.basis),
                                     k1_fraction=float(rep.k1_fraction),
                                     stage1_dim=int(rep.stage1_dim)),
        "+ early_exit=False": dict(basis=str(rep.basis),
                                   k1_fraction=float(rep.k1_fraction),
                                   early_exit=False),
    }
    achou = None
    for nome, kw in variantes.items():
        try:
            r = rec2(build_engine(X, dim=d, metric="cosine", k=10, **kw))
        except Exception as exc:  # noqa: BLE001
            print(f"    {nome:<24} ERRO {str(exc)[:60]}")
            continue
        marca = "← reproduz o roteado" if abs(r - r_rot) < 1e-6 else ""
        print(f"    {nome:<24} recall {r:.4f}  {marca}")
        if abs(r - r_rot) < 1e-6 and achou is None:
            achou = nome
    out["isolamento"] = {"recall_roteado": round(r_rot, 4),
                         "reproduzido_por": achou}
    if achou and "stage1_dim" in achou:
        print("\n    ✗ O ganho NÃO vem dos parâmetros anunciados. Vem de "
              "stage1_dim = d (dimensão completa no estágio 1), que remove o "
              "prefilter heurístico. O relatório reporta basis e k1, mas é "
              "este terceiro parâmetro que decide o recall.")
    elif achou:
        print(f"\n    ✓ o ganho é reproduzível por: {achou}")

    RESULTS["normalize"] = out
    return out


# ===========================================================================
# main
# ===========================================================================

def main() -> int:
    ap = argparse.ArgumentParser()
    ap.add_argument("--rapido", action="store_true")
    ap.add_argument("--json", default="matematica_da_stack.json")
    args = ap.parse_args()

    titulo("A MATEMÁTICA DA STACK — X-FACTOR, NANO E NORMALIZE")
    print(" Três peças, três estruturas. Cada teste mede a peça contra a")
    print(" própria definição matemática dela.")

    try:
        import winnex_nano as wn
        print(f"\n winnex-nano {wn.__version__}")
    except ImportError:
        print("\n winnex-nano ausente: pip install winnex-nano", file=sys.stderr)
        return 1

    n1_tokenizador(wn)
    x1_power_iteration(wn)
    g1_g2_normalize(args.rapido)

    titulo("RELATÓRIO")
    nn = RESULTS.get("nano", {})
    if nn:
        print(f"  N1  round-trip {nn.get('roundtrip_bytes')} | "
              f"(y,z) e modos j>0 são rotações: "
              f"{'sim' if (nn.get('erro_yz') or 1) < 1e-5 else 'não'}")
    xr = RESULTS.get("xfactor_rank", {})
    if xr.get("casos"):
        ok = all(c["rank"] == c["rank_exato"] for c in xr["casos"])
        print(f"  X1  rank efetivo exato: {'sim' if ok else 'não'} "
              f"({len(xr['casos'])} casos)")
    idem = xr.get("idempotencia", [])
    if idem:
        print(f"  X2  ‖P²−P‖ cresce com o rank: {idem[0]['idem']:.1e} → "
              f"{idem[-1]['idem']:.1e}")
    nz = RESULTS.get("normalize", {})
    if nz.get("isolamento"):
        print(f"  G2  recall do roteador {nz['isolamento']['recall_roteado']} "
              f"reproduzido por: {nz['isolamento']['reproduzido_por']}")

    with open(args.json, "w", encoding="utf-8") as f:
        json.dump(RESULTS, f, indent=1, ensure_ascii=False, default=str)
    print(f"\n[ok] {args.json}")
    return 0


if __name__ == "__main__":
    sys.exit(main())