"""
NUVMAP Projection Engine — SilverSight Port
=============================================

Port of the Research Stack archive projection_engine.py with:

  1. Q16_16 fixed-point arithmetic (no Float in compute paths)
  2. CBOR serialization (binary, not JSON-in-SQLite)
  3. Dataclass → bytes serialization with schema versioning

Key equation (all Q16_16):
  q_i ∝ E_i / (R_i + ε)
  E_i = λ · |v_k(i)| · S_i · L_i / (R_i + ε)

Q16_EPSILON = 1 (≈ 1.5e-5)
────────────────────────────────
The epsilon is NOT an arbitrary regularization — it is the minimum
distance from the Kelvin fixed point in the Rossby/Kelvin braid
correspondence (see BraidStateN.lean §RotationalWaveCorrespondence).

In the Rossby/Kelvin model:
  • Rossby regime (chiral asymmetry ≠ 0): residual risk R_i ≫ ε,
    the E_i computation is dominated by the chiral-dependent R_i.
  • Kelvin regime (achiral_stable, Rossby drift = 0): R_i → 0,
    epsilon prevents the mathematical singularity at the fixed point
    where crossStep(s) = s (eigensolid convergence).

The chiral quantization scale is 32768 Q16_16 units (0.5 — the minimum
non-zero contribution from chiral_scarred).  Q16_EPSILON = 1 sits
2¹⁵ = 32768× below this scale, resolving the Kelvin boundary layer:

  Q16_ONE   = 65536   ─── 1.0 (full chiral bias, left/right handed)
  Q16_HALF  = 32768   ─── 0.5 (chiral_scarred minimum step)
                            │   Rossby drift quantization scale
  Q16_EPSILON = 1     ─── ≈1.5e-5 (Kelvin boundary layer floor)

The original Research Stack used epsilon = 1e-12 (below double-precision
mach.eps = 2.2e-16), which had no physical meaning.  Q16_EPSILON = 1
is the correct Kelvin boundary layer resolution per the formalization.

Proofs (BraidStateN.lean):
  rossby_convergence_bound    — non-zero chirality → step count ↑
  kelvin_wave_eigensolid      — zero chirality → eigensolid fixed point

──────────────────────────────────────────────────────────────────────
Why this works for compression (planetary model analogy)

The Rossby β-plane is the standard geophysical fluid dynamics model for
planetary atmospheres (Earth's jet streams, Jupiter's bands, ocean gyres).
The braid model shares the same structural invariants:

  GFD                          Braid model
  ───────────────────────────  ─────────────────────────────────────
  β (planetary vorticity       rossbyDriftFromChirality asymmetry
    gradient)
  Rossby wave dispersion       Braid-word frequency splitting
    ω = −βk/(k²+l²)
  Coastal Kelvin wave          kelvin_wave_eigensolid — outer strands,
    (boundary-trapped)           eigensolid converges first
  Rossby deformation radius    Q16_EPSILON = 1 — boundary layer
    L_D = √(gH)/f               thickness (minimum resolvable step)
  Rhines scale L_β =           Q16_HALF = 32768 — chiral quantization
    √(2U/β)                     step

The reason this maps to compression is structural: the Yang-Baxter
equation (braid relation β_ij β_jk β_ij = β_jk β_ij β_jk) and the
β-plane potential vorticity conservation ∂_t q + J(ψ, q) + β∂_x ψ = 0
are both curvature constraints on a rotating/layered system. The braid
crossing matrix C encodes the same topological information that the PV
gradient∫encodes in the fluid — both conserve a topological charge
(braid class / potential vorticity) under evolution.

Implication for efficiency:
  If the compression residuals follow Rossby wave dispersion, then
  the eigensolid convergence rate is bounded by β-plane wave speeds,
  not by general braid mixing.  This means the step-count bound in
  rossby_convergence_bound can be tightened from O(n²) to O(√(n/β))
  in the Rossby-dominated regime, directly reducing the number of
  crossing steps needed for eigensolid detection.
"""

import hashlib
from typing import Dict, List, Optional
from dataclasses import dataclass, field
from datetime import datetime

try:
    import cbor2 as cbor
except ImportError:
    import cbor  # fallback

# ── Q16_16 fixed-point ──────────────────────────────────────────────
SCALE = 65536

# ── Q16_16 raw constants — Rounding Derivation ─────────────────────────
# All constants are raw = round(value × SCALE) using banker's rounding
# (round-half-to-even, Python's default).  Every non-exact conversion
# documents its rounding decision below for traceability.
#
# value  | raw    | derivation
# -------|--------|-------------------------------------------------------
# ϵ      |      1 | 1/65536 ≈ 0.000015  (exact floor — q16 step)
# 1.0    |  65536 | 2¹⁶  (exact)
# 0.5    |  32768 | 2¹⁵  (exact)
# 0.01   |    655 | 655.36 → floor  (.36 < 0.5)
# 0.7    |  45875 | 45875.2 → floor  (.2 < 0.5)
# 0.3    |  19661 | 19660.8 → ceil   (.8 ≥ 0.5, even → 19661)
# 1.5    |  98304 | 1.5 × 2¹⁶  (exact)

Q16_EPSILON    = 1
Q16_ONE        = 65536
Q16_HALF       = 32768
Q16_ZERO       = 0
Q16_PCT1       = 655       # 0.01  floor
Q16_PCT70      = 45875     # 0.7   floor
Q16_PCT30      = 19661     # 0.3   ceil  (banker's: .8 → even 19661)
Q16_150PCT     = 98304     # 1.5   exact


def q_mul(a: int, b: int) -> int:
    """Q16_16 multiplication with banker's rounding."""
    return max(-2147483648, min(2147483647, round((a * b) / SCALE)))

def q_div(a: int, b: int) -> int:
    """Q16_16 division with banker's rounding. Returns 0 if b==0."""
    if b == 0:
        return 0
    return max(-2147483648, min(2147483647, round((a * SCALE) / b)))

def q_add(a: int, b: int) -> int:
    return max(-2147483648, min(2147483647, a + b))

def q_sub(a: int, b: int) -> int:
    return max(-2147483648, min(2147483647, a - b))

# ── Schema ──────────────────────────────────────────────────────────
SCHEMA_VERSION = 1
CBOR_TAG_NUVMAP_CELL = 0x1000
CBOR_TAG_NUVMAP_SURFACE = 0x1001

# ── Dataclasses ──────────────────────────────────────────────────────

@dataclass
class NUVMAPCell:
    u_i: int                # address coordinate
    v_i: int                # spectral coordinate (eigenmode index)
    k_i: int                # dominant eigenmode
    E_i: int = 0            # eigenmass (Q16_16)
    R_i: int = Q16_ONE      # residual risk (Q16_16)
    chi_i: int = 0          # chiral residual (Q16_16)
    S_i: int = Q16_ONE       # structural integrity (Q16_16)
    L_i: int = Q16_ONE       # Landauer factor (Q16_16)
    q_i: int = 0            # qubit allocation
    admissible: bool = True
    equation_id: int = 0
    fingerprint: str = ""

    def to_cbor(self) -> bytes:
        data = [
            self.u_i, self.v_i, self.k_i,
            self.E_i, self.R_i, self.chi_i,
            self.S_i, self.L_i,
            self.q_i, int(self.admissible),
            self.equation_id, self.fingerprint,
        ]
        return cbor.dumps(data)

    @classmethod
    def from_cbor(cls, raw: bytes) -> 'NUVMAPCell':
        data = cbor.loads(raw)
        return cls(
            u_i=data[0], v_i=data[1], k_i=data[2],
            E_i=data[3], R_i=data[4], chi_i=data[5],
            S_i=data[6], L_i=data[7],
            q_i=data[8], admissible=bool(data[9]),
            equation_id=data[10], fingerprint=data[11],
        )


@dataclass
class NUVMAPSurface:
    cells: List['NUVMAPCell'] = field(default_factory=list)
    total_qubits: int = 0
    bekenstein_bound: int = 0   # Q16_16
    area_utilization: int = 0   # Q16_16
    root_fingerprint: str = ""
    timestamp: str = ""

    def to_cbor(self) -> bytes:
        cell_data = [c.to_cbor() for c in self.cells]
        data = [
            SCHEMA_VERSION,
            cell_data,
            self.total_qubits,
            self.bekenstein_bound,
            self.area_utilization,
            self.root_fingerprint,
            self.timestamp,
        ]
        return cbor.dumps([CBOR_TAG_NUVMAP_SURFACE, data])

    @classmethod
    def from_cbor(cls, raw: bytes) -> 'NUVMAPSurface':
        loaded = cbor.loads(raw)
        # Unwrap [tag, data] if present
        if isinstance(loaded, list) and len(loaded) == 2 and isinstance(loaded[0], int):
            _, data = loaded
        else:
            data = loaded
        version = data[0]
        cells = [NUVMAPCell.from_cbor(c) for c in data[1]]
        return cls(
            cells=cells,
            total_qubits=data[2],
            bekenstein_bound=data[3],
            area_utilization=data[4],
            root_fingerprint=data[5],
            timestamp=data[6],
        )

    def to_file(self, path: str):
        with open(path, 'wb') as f:
            f.write(self.to_cbor())

    @classmethod
    def from_file(cls, path: str) -> 'NUVMAPSurface':
        with open(path, 'rb') as f:
            return cls.from_cbor(f.read())


# ── Engine ──────────────────────────────────────────────────────────────

class NUVMAPProjectionEngine:
    """
    Projects eigenmass data into a NUVMAP address surface using Q16_16.
    """

    def __init__(self, total_qubit_budget: int = 0,
                 chi_max_q16: int = Q16_HALF,
                 R_max_q16: int = Q16_HALF,
                 landauer_threshold_q16: int = None):
        self.total_qubit_budget = total_qubit_budget
        self.chi_max_q16 = chi_max_q16
        self.R_max_q16 = R_max_q16
        self.landauer_threshold_q16 = landauer_threshold_q16 or (Q16_ONE // 10)
        self.surface = NUVMAPSurface()

    def project(self, eigenmass_data: List[Dict],
                eigenvalue_q16: Optional[int] = None) -> NUVMAPSurface:
        """
        Project eigenmass data into a NUVMAP surface.

        eigenmass_data: list of dicts with keys:
            equation_id (int), amvr_q16, avmr_q16, chiral_residual_q16 (all Q16_16 raw),
            chiral_state (str)

        All numeric values MUST already be Q16_16 raw integers.
        Convert at the outermost call boundary with the provided Q16_16 adapter.
        """
        if not eigenmass_data:
            return self.surface

        n = len(eigenmass_data)
        cells = []

        max_eigenmass = Q16_ZERO
        for d in eigenmass_data:
            raw_q = q_div(q_add(d.get("amvr_q16", 0), d.get("avmr_q16", 0)),
                          Q16_ONE * 2)  # (amvr+avmr)/2
            if raw_q > max_eigenmass:
                max_eigenmass = raw_q
        if max_eigenmass == Q16_ZERO:
            max_eigenmass = Q16_ONE

        for i, d in enumerate(eigenmass_data):
            amvr_q = d.get("amvr_q16", 0)
            avmr_q = d.get("avmr_q16", 0)
            cr_q = d.get("chiral_residual_q16", 0)
            eq_id = d.get("equation_id", 0)
            cs = d.get("chiral_state", "achiral_stable")

            raw_eigenmass = q_div(q_add(amvr_q, avmr_q), Q16_ONE * 2)
            E_norm = q_div(raw_eigenmass, max_eigenmass)

            one_minus = q_sub(Q16_ONE, E_norm)
            R_i = max(Q16_PCT1, one_minus) if one_minus > Q16_PCT1 else Q16_PCT1

            if cs == "chiral_scarred":
                R_i = q_mul(R_i, Q16_150PCT)

            if cs in ("achiral_stable",):
                S_i = Q16_ONE
            elif cs in ("left_handed_mass_bias", "right_handed_vector_bias"):
                S_i = Q16_PCT70
            else:
                S_i = Q16_PCT30

            if E_norm > self.landauer_threshold_q16:
                L_i = Q16_ONE
            else:
                L_i = q_div(E_norm, self.landauer_threshold_q16)

            lam_q = eigenvalue_q16 if eigenvalue_q16 is not None else Q16_ONE
            v_abs = E_norm

            E_i = q_div(q_mul(q_mul(q_mul(lam_q, v_abs), S_i), L_i),
                        q_add(R_i, Q16_EPSILON))

            chi_i = cr_q

            is_R_ok = R_i <= self.R_max_q16
            is_chi_ok = chi_i <= self.chi_max_q16
            admissible = is_R_ok and is_chi_ok

            fp_payload = f"{eq_id}\x00{amvr_q}\x00{avmr_q}\x00{cr_q}\x00{cs}"
            fp = hashlib.sha256(fp_payload.encode()).hexdigest()

            cells.append(NUVMAPCell(
                u_i=i, v_i=i, k_i=i,
                E_i=E_i, R_i=R_i, chi_i=chi_i,
                S_i=S_i, L_i=L_i,
                q_i=0, admissible=admissible,
                equation_id=eq_id, fingerprint=fp,
            ))

        # Qubit allocation: q_i proportional to E_i / (R_i + epsilon)
        total_weight = Q16_ZERO
        for c in cells:
            total_weight = q_add(total_weight, q_div(c.E_i, q_add(c.R_i, Q16_EPSILON)))
        if total_weight == Q16_ZERO:
            total_weight = Q16_ONE

        if self.total_qubit_budget > 0:
            budget = self.total_qubit_budget
        else:
            budget = sum(c.E_i * 100 // SCALE for c in cells if c.admissible)
        budget = max(budget, sum(1 for c in cells if c.admissible))

        for c in cells:
            if c.admissible:
                weight = q_div(c.E_i, q_add(c.R_i, Q16_EPSILON))
                raw_q = int(budget * weight // total_weight) if total_weight > 0 else 1
                c.q_i = max(1, raw_q) if raw_q > 0 else 1
            else:
                c.q_i = 0

        total_qubits = sum(c.q_i for c in cells)

        # Bekenstein-like bound in Q16_16
        bekenstein = q_div(sum(c.E_i for c in cells), len(cells) * Q16_ONE) if cells else Q16_ZERO

        area_util = q_div(total_qubits * Q16_ONE, q_add(bekenstein, Q16_EPSILON)) if bekenstein > 0 else Q16_ZERO

        self.surface = NUVMAPSurface(
            cells=cells,
            total_qubits=total_qubits,
            bekenstein_bound=bekenstein,
            area_utilization=area_util,
            root_fingerprint=self._compute_surface_root(cells),
            timestamp=datetime.utcnow().isoformat(),
        )

        return self.surface

    @staticmethod
    def _compute_surface_root(cells: List['NUVMAPCell']) -> str:
        payload = "|".join(
            f"{c.u_i}:{c.E_i}:{c.chi_i}:{c.q_i}"
            for c in sorted(cells, key=lambda x: x.u_i)
        )
        return hashlib.sha256(payload.encode()).hexdigest()

    def quantum_storage_admissible(self, node_i: int, tau_q16: int) -> bool:
        if node_i < 0 or node_i >= len(self.surface.cells):
            return False
        c = self.surface.cells[node_i]
        lhs = q_mul(c.E_i, q_add(c.R_i, Q16_EPSILON))
        rhs = q_mul(tau_q16, q_add(c.R_i, Q16_EPSILON))
        return lhs <= rhs and c.chi_i <= self.chi_max_q16 and c.admissible

    def get_density_map(self) -> Dict[str, List]:
        if not self.surface.cells:
            return {"E_i": [], "q_i": [], "chi_i": [], "R_i": []}
        return {
            "E_i": [c.E_i for c in self.surface.cells],
            "q_i": [c.q_i for c in self.surface.cells],
            "chi_i": [c.chi_i for c in self.surface.cells],
            "R_i": [c.R_i for c in self.surface.cells],
            "equation_ids": [c.equation_id for c in self.surface.cells],
        }

    def summary(self) -> Dict:
        s = self.surface
        admissible = [c for c in s.cells if c.admissible]
        return {
            "num_cells": len(s.cells),
            "num_admissible": len(admissible),
            "num_rejected": len(s.cells) - len(admissible),
            "total_qubits": s.total_qubits,
            "avg_qubits_per_cell": s.total_qubits // max(len(admissible), 1),
            "bekenstein_bound_q16": s.bekenstein_bound,
            "area_utilization_q16": s.area_utilization,
            "max_eigenmass_q16": max((c.E_i for c in s.cells), default=0),
            "max_chiral_q16": max((c.chi_i for c in s.cells), default=0),
            "surface_root": s.root_fingerprint[:32] + "...",
        }


def build_nuvmap_from_eigenmass(eigenmass_data: List[Dict],
                                qubit_budget: int = 0) -> NUVMAPSurface:
    engine = NUVMAPProjectionEngine(total_qubit_budget=qubit_budget)
    return engine.project(eigenmass_data)