#!/usr/bin/env python3
"""
NUVMAP Equivalence Test — SilverSight vs Research Stack Archive

Compares the Q16_16 SilverSight port against the original float-based
Research Stack implementation.

Key difference accounted for:
  Original: epsilon = 1e-12 (unphysically small)
  Ported:   Q16_EPSILON = 1 ≈ 1.5e-5  (smallest Q16_16 step)

This creates a systematic ~±1 Q16_16 step divergence in E_i.
All other divergence is pure Q16_16 quantization (±0.5 LSB).

Run: python -m tests.test_nuvmap_equivalence
"""

import sys
import importlib.util

# Import original (float) from Research Stack
spec_orig = importlib.util.spec_from_file_location(
    "original_nuvmap",
    "/home/allaun/Research Stack/5-Applications/cff/nuvmap/projection_engine.py"
)
orig_mod = importlib.util.module_from_spec(spec_orig)
spec_orig.loader.exec_module(orig_mod)

OriginalEngine = orig_mod.NUVMAPProjectionEngine
OriginalSurface = orig_mod.NUVMAPSurface

# Import ported (Q16_16) from SilverSight
spec_ptd = importlib.util.spec_from_file_location(
    "ported_nuvmap",
    "/home/allaun/SilverSight/python/nuvmap/projection_engine.py"
)
ptd_mod = importlib.util.module_from_spec(spec_ptd)
spec_ptd.loader.exec_module(ptd_mod)

PortedEngine = ptd_mod.NUVMAPProjectionEngine
PortedSurface = ptd_mod.NUVMAPSurface

Q16_ONE   = ptd_mod.Q16_ONE
Q16_HALF  = ptd_mod.Q16_HALF
Q16_PCT1  = ptd_mod.Q16_PCT1
Q16_EPSILON = ptd_mod.Q16_EPSILON
SCALE     = ptd_mod.SCALE

# ── helpers ──────────────────────────────────────────────────────────

def f2q(f: float) -> int:
    import math
    if math.isnan(f) or math.isinf(f):
        return 0
    return max(-2147483648, min(2147483647, round(f * SCALE)))

def q2f(q: int) -> float:
    return q / SCALE


TEST_DATA = [
    {"equation_id": 1, "amvr": 0.8,  "avmr": 0.75, "chiral_residual": 0.05, "chiral_state": "achiral_stable"},
    {"equation_id": 2, "amvr": 0.3,  "avmr": 0.4,  "chiral_residual": 0.2,  "chiral_state": "left_handed_mass_bias"},
    {"equation_id": 3, "amvr": 0.1,  "avmr": 0.15, "chiral_residual": 0.6,  "chiral_state": "chiral_scarred"},
    {"equation_id": 4, "amvr": 0.5,  "avmr": 0.5,  "chiral_residual": 0.1,  "chiral_state": "right_handed_vector_bias"},
    {"equation_id": 5, "amvr": 0.9,  "avmr": 0.85, "chiral_residual": 0.02, "chiral_state": "achiral_stable"},
]


def to_q16(data: list) -> list:
    return [{
        "equation_id": d["equation_id"],
        "amvr_q16": f2q(d["amvr"]),
        "avmr_q16": f2q(d["avmr"]),
        "chiral_residual_q16": f2q(d["chiral_residual"]),
        "chiral_state": d["chiral_state"],
    } for d in data]


# ── Comparators ─────────────────────────────────────────────────────

def cmp_float_vs_q16(oc, pc, label="", rel_tol=0.02) -> bool:
    """Compare float cell with Q16_16 cell. Tolerates ε quantization."""
    fields = [
        ("E_i",   oc.E_i,   pc.E_i),
        ("R_i",   oc.R_i,   pc.R_i),
        ("chi_i", oc.chi_i, pc.chi_i),
        ("S_i",   oc.S_i,   pc.S_i),
        ("L_i",   oc.L_i,   pc.L_i),
    ]
    ok = True
    for name, fval, qval in fields:
        f2 = q2f(qval)
        if fval == 0.0:
            diff = abs(f2)
        else:
            diff = abs(f2 - fval)
        reld = diff / max(abs(fval), 1e-12)
        if reld > rel_tol and diff > 2 / SCALE:
            print(f"  {label}{name}: float={fval:.6f}  q16={f2:.6f}  diff={diff:.6f}  rel={reld:.4f}")
            ok = False

    q_diff = abs(pc.q_i - oc.q_i)
    q_reld = q_diff / max(oc.q_i, 1)
    if q_diff > 1 and q_reld > 0.01:
        print(f"  {label}q_i: float={oc.q_i}  q16={pc.q_i}  rel={q_reld:.4f}")
        ok = False
    if oc.admissible != pc.admissible:
        print(f"  {label}admissible: float={oc.admissible}  q16={pc.admissible}")
        ok = False
    return ok


def cmp_q16_vs_q16(a, b, label="", tol=0) -> bool:
    """Compare two Q16_16 cells (exact match for round-trip)."""
    for name in ("u_i", "v_i", "k_i", "E_i", "R_i", "chi_i", "S_i", "L_i", "q_i", "equation_id"):
        va, vb = getattr(a, name), getattr(b, name)
        if abs(va - vb) > tol:
            print(f"  {label}{name}: {va} vs {vb}")
            return False
    if a.admissible != b.admissible:
        print(f"  {label}admissible: {a.admissible} vs {b.admissible}")
        return False
    if a.fingerprint != b.fingerprint:
        print(f"  {label}fingerprint: {a.fingerprint[:16]} vs {b.fingerprint[:16]}")
        return False
    return True


def cmp_surface_float_vs_q16(o, p, rel_tol=0.02) -> bool:
    if len(o.cells) != len(p.cells):
        print(f"Cell count: float={len(o.cells)}  q16={len(p.cells)}")
        return False
    ok = True
    for i, (oc, pc) in enumerate(zip(o.cells, p.cells)):
        if not cmp_float_vs_q16(oc, pc, f"cell[{i}] ", rel_tol):
            ok = False
    if o.total_qubits != p.total_qubits:
        reld = abs(p.total_qubits - o.total_qubits) / max(o.total_qubits, 1)
        if reld > 0.01:
            print(f"total_qubits: float={o.total_qubits}  q16={p.total_qubits}  rel={reld:.4f}")
            ok = False
    for name in ("bekenstein_bound", "area_utilization"):
        fv = getattr(o, name)
        qv = q2f(getattr(p, name))
        reld = abs(qv - fv) / max(abs(fv), 1e-12)
        if reld > rel_tol:
            print(f"{name}: float={fv:.6f}  q16={qv:.6f}  rel={reld:.4f}")
            ok = False
    return ok


def cmp_surface_q16_vs_q16(a, b, tol=0) -> bool:
    if len(a.cells) != len(b.cells):
        print(f"Cell count: {len(a.cells)} vs {len(b.cells)}")
        return False
    ok = True
    for i, (ac, bc) in enumerate(zip(a.cells, b.cells)):
        if not cmp_q16_vs_q16(ac, bc, f"cell[{i}] ", tol):
            ok = False
    for name in ("total_qubits", "bekenstein_bound", "area_utilization", "root_fingerprint"):
        if getattr(a, name) != getattr(b, name):
            print(f"{name}: {getattr(a, name)} vs {getattr(b, name)}")
            ok = False
    return ok


# ── Tests ───────────────────────────────────────────────────────────

def test_equivalence():
    print("═══ Equivalence: float vs Q16_16 ═══")
    data_q16 = to_q16(TEST_DATA)

    oe = OriginalEngine(total_qubit_budget=0, chi_max=0.5, R_max=0.5, landauer_threshold=0.1)
    os = oe.project(TEST_DATA)

    pe = PortedEngine(total_qubit_budget=0,
                      chi_max_q16=Q16_HALF,
                      R_max_q16=Q16_HALF,
                      landauer_threshold_q16=Q16_ONE // 10)
    ps = pe.project(data_q16)

    print(f"  Float:  cells={len(os.cells)}  qubits={os.total_qubits}  B={os.bekenstein_bound:.4f}  U={os.area_utilization:.4f}")
    print(f"  Q16_16: cells={len(ps.cells)}  qubits={ps.total_qubits}  B={q2f(ps.bekenstein_bound):.4f}  U={q2f(ps.area_utilization):.4f}")

    # rel_tol=0.02 = 2% — eps is 1e-12 vs 1.5e-5 which propagates through E_i formula
    ok = cmp_surface_float_vs_q16(os, ps, rel_tol=0.02)
    print(f"  → {'PASS' if ok else 'FAIL'}")
    return ok


def test_cbor_roundtrip():
    print("═══ CBOR round-trip ═══")
    data_q16 = to_q16(TEST_DATA)
    engine = PortedEngine(total_qubit_budget=100)
    s1 = engine.project(data_q16)
    blob = s1.to_cbor()
    print(f"  Size: {len(blob)} bytes")
    s2 = PortedSurface.from_cbor(blob)
    ok = cmp_surface_q16_vs_q16(s1, s2, tol=0)
    print(f"  → {'PASS' if ok else 'FAIL'}")
    return ok


def test_file_roundtrip():
    import tempfile, os
    print("═══ File round-trip ═══")
    data_q16 = to_q16(TEST_DATA)
    engine = PortedEngine(total_qubit_budget=100)
    s1 = engine.project(data_q16)
    with tempfile.NamedTemporaryFile(suffix='.cbor', delete=False) as f:
        path = f.name
    try:
        s1.to_file(path)
        s2 = PortedSurface.from_file(path)
        ok = cmp_surface_q16_vs_q16(s1, s2, tol=0)
        print(f"  → {'PASS' if ok else 'FAIL'}")
        return ok
    finally:
        os.unlink(path)


def test_admissibility_gate():
    print("═══ Admissibility gate ═══")
    data_q16 = to_q16(TEST_DATA)
    oe = OriginalEngine(chi_max=0.5, R_max=0.5, landauer_threshold=0.1)
    oe.project(TEST_DATA)
    pe = PortedEngine(chi_max_q16=Q16_HALF, R_max_q16=Q16_HALF,
                      landauer_threshold_q16=Q16_ONE // 10)
    pe.project(data_q16)
    ok = True
    for i in range(len(TEST_DATA)):
        r1 = oe.quantum_storage_admissible(i, 1.0)
        r2 = pe.quantum_storage_admissible(i, Q16_ONE)
        if r1 != r2:
            print(f"  Gate mismatch cell[{i}]: float={r1}  q16={r2}")
            ok = False
    print(f"  → {'PASS' if ok else 'FAIL'}")
    return ok


def test_numerical_analysis():
    """Detailed breakdown of what diverges and why."""
    print("═══ Numerical analysis ═══")
    data_q16 = to_q16(TEST_DATA)

    oe = OriginalEngine(total_qubit_budget=0, chi_max=0.5, R_max=0.5, landauer_threshold=0.1)
    os = oe.project(TEST_DATA)
    pe = PortedEngine(total_qubit_budget=0,
                      chi_max_q16=Q16_HALF,
                      R_max_q16=Q16_HALF,
                      landauer_threshold_q16=Q16_ONE // 10)
    ps = pe.project(data_q16)

    print(f"  epsilon float: {oe.epsilon}     (unphysically small)")
    print(f"  epsilon Q16_16: {Q16_EPSILON} = {q2f(Q16_EPSILON):.8f}")
    print(f"  E_i formula: lam * v_abs * S * L / (R + ε)")
    print()

    chiral_states = [d["chiral_state"] for d in TEST_DATA]
    for i, (oc, pc) in enumerate(zip(os.cells, ps.cells)):
        cs = chiral_states[i]
        print(f"  cell[{i}] eq={oc.equation_id} {cs}")
        print(f"    float: E_i={oc.E_i:.6f}  R_i={oc.R_i:.6f}  L_i={oc.L_i:.6f}  q_i={oc.q_i}")
        print(f"    q16:   E_i={pc.E_i} ({q2f(pc.E_i):.6f})  R_i={pc.R_i} ({q2f(pc.R_i):.6f})  L_i={pc.L_i} ({q2f(pc.L_i):.6f})  q_i={pc.q_i}")
        r_q16 = q2f(pc.R_i) + q2f(Q16_EPSILON)
        r_float = oc.R_i + oe.epsilon
        print(f"    R+ε: float={r_float:.8f}  q16={r_q16:.8f}")
        print()
    return True


if __name__ == "__main__":
    results = {}
    results["equivalence"] = test_equivalence()
    results["cbor"] = test_cbor_roundtrip()
    results["file"] = test_file_roundtrip()
    results["gate"] = test_admissibility_gate()
    results["analysis"] = test_numerical_analysis()

    print()
    n_pass = sum(1 for v in results.values() if v)
    n_total = len(results)
    print(f"{n_pass}/{n_total} passed")
    if all(results.values()):
        sys.exit(0)
    else:
        sys.exit(1)