SilverSight/python/nuvmap/projection_engine.py

"""
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 * SCALE) / (R_i + epsilon)   proportional allocation
  E_i = lam * v_abs * S_i * L_i / (R_i + epsilon)
"""

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 ────────────────────────────────────────────
# ── Q16_16 Constant Derivation ──────────────────────────────────────
# All constants are Q16_16 raw integers = round(value × 65536).
# Banker's rounding (round-half-to-even) used throughout.
#
# value  | raw    | formula
# -------|--------|-----------------------------------
# ϵ      |      1 | 1/65536 ≈ 0.000015 (minimum step)
# 1.0    |  65536 | exact (2^16)
# 0.5    |  32768 | exact
# 0.01   |    655 | 655.36 → 655  (floor, .36 < 0.5)
# 0.7    |  45875 | 45875.2 → 45875 (floor, .2 < 0.5)
# 0.3    |  19661 | 19660.8 → 19661 (ceil, .8 ≥ 0.5 → even 19661)
# 1.5    |  98304 | 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 rounding."""
    return max(-2147483648, min(2147483647, round((a * b) / SCALE)))

def q_div(a: int, b: int) -> int:
    """Q16_16 division with 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))

Q16_EPSILON = 1  # ≈ 1.5e-5 in Q16_16
Q16_ONE = SCALE      # 1.0
Q16_HALF = SCALE // 2  # 0.5
Q16_ZERO = 0

# ── 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)