"""
multi_surface_packer.py — ΔΦΓΛ Multi-Surface Packing Compressor

Unified pipeline implementing the multi-surface packing Lagrangian:

    min_{γ ∈ Γ, φ ∈ Φ}  [  ‖Δ(γ, x)‖₁  +  α·|φ|  +  β·K(γ)  ]
    s.t.   ⟨ψ_x | ψ_x̂⟩² ≥ 1 − ε_λ         (WaveProbe coherence gate)
           gccl_swap(|x|, ‖Δ‖₁, risk)    (GCCL admissibility gate)

Each tunable parameter (α, β, ε_λ, risk) is surfaced for the multi-surface
packing optimization — changing any one shifts the Pareto front across all
three surfaces simultaneously.

Metadata extraction is implemented per surface:
  Δ (Data):    PTOS manifest classification + block statistics
  Φ (Spectral): Golden-angle sampling + phase coherence + spectral energy
  Γ (Program):  GCCL program codon distribution + complexity metrics

Lean source of truth references:
  - DeltaGCLCompression.lean   — PTOS dictionary, delta encoding, compression stats
  - GoldenAngleEncoding.lean   — golden-angle phase sampling, WaveProbeSample
  - Waveprobe.lean             — QUBO projector, overlap energy |⟨ψc|ψp⟩|²
  - DeltaPhiGammaKLambda.lean  — DPGState, dpgTransform, DoctrineAdmissible
  - gccl_waveprobe.py          — GCCL gates, WaveProbe sampling in Python

All arithmetic is Q16_16 fixed-point (no Float in compute paths).
"""

from __future__ import annotations

import hashlib
import json
import math
import statistics
import struct
import time
from dataclasses import dataclass, field
from enum import IntEnum
from typing import Dict, List, Optional, Tuple


# ── Q16_16 Fixed-Point ─────────────────────────────────────────────────────

Q16_SCALE = 65536
Q16_MAX = 32767
Q16_MIN = -32768

GOLDEN_ANGLE_STEP = 40503  # 1/φ × 65536 (from GoldenAngleEncoding.lean)


def q16(x: int) -> int:
    raw = x * Q16_SCALE
    return max(Q16_MIN, min(Q16_MAX, raw))


def q16f(x: float) -> int:
    raw = int(x * Q16_SCALE)
    return max(Q16_MIN, min(Q16_MAX, raw))


def q16_mul(a: int, b: int) -> int:
    return max(Q16_MIN, min(Q16_MAX, (a * b) // Q16_SCALE))


def q16_add(a: int, b: int) -> int:
    return max(Q16_MIN, min(Q16_MAX, a + b))


def q16_sub(a: int, b: int) -> int:
    return max(Q16_MIN, min(Q16_MAX, a - b))


def q16_div(a: int, b: int) -> int:
    if b == 0:
        return 0
    return max(Q16_MIN, min(Q16_MAX, (a * Q16_SCALE) // b))


# ── §1  Metadata: PTOS Manifest Classification (Δ surface) ─────────────────
# Maps raw bytes to PTOS dictionary fields per DeltaGCLCompression.lean §1

class PTOSLayer(IntEnum):
    CORE = 0
    CARRY = 1
    RULE = 2
    STORE = 3
    EXTERNAL = 4


class PTOSDomain(IntEnum):
    COMPUTE = 0
    TOKEN = 1
    RULE = 2
    STORE = 3
    POWER = 4
    COMMS = 5
    MATERIAL = 6
    DATA = 7
    CLOCK = 8
    TEST = 9


class PTOSTier(IntEnum):
    SINGULARITY = 0
    PLASMA = 1
    CRYSTALLINE = 2
    FOAM = 3
    GOVERNANCE = 4
    RESEARCH = 5


class PTOSCondition(IntEnum):
    STABLE = 0
    EXPERIMENTAL = 1
    EXTREME = 2
    DRAFT = 3
    ARCHIVED = 4
    STERILE = 5


PTOS_LAYER_INDEX = {v: k for k, v in PTOSLayer.__members__.items()}
PTOS_DOMAIN_INDEX = {v: k for k, v in PTOSDomain.__members__.items()}
PTOS_TIER_INDEX = {v: k for k, v in PTOSTier.__members__.items()}
PTOS_CONDITION_INDEX = {v: k for k, v in PTOSCondition.__members__.items()}


@dataclass
class PTOSManifest:
    """PTOS manifest classification — metadata extracted from raw data.

    Matches DeltaGCLCompression.lean: structure PTOSManifest
    """
    layer: PTOSLayer
    domain: PTOSDomain
    tier: PTOSTier
    condition: PTOSCondition

    def to_bytes(self) -> bytes:
        """PTOS dictionary encoding: 4 bytes (matches applyPTOSDictionary)."""
        return bytes([
            self.layer.value,
            self.domain.value,
            self.tier.value,
            self.condition.value,
        ])

    def to_dict(self) -> Dict:
        return {
            'layer': self.layer.name.lower(),
            'domain': self.domain.name.lower(),
            'tier': self.tier.name.lower(),
            'condition': self.condition.name.lower(),
        }


def extract_ptos_manifest(data: bytes) -> PTOSManifest:
    """Classify raw bytes into PTOS manifest metadata.

    Heuristic extraction based on statistical features of the data.
    Each PTOS field is determined by a different feature axis:
      Layer   ← structural depth (entropy range)
      Domain  ← byte distribution moments
      Tier    ← complexity (repetition ratio)
      Condition ← stability (delta from uniform)

    This is the Δ-surface metadata extraction.
    """
    n = len(data)
    if n == 0:
        return PTOSManifest(
            layer=PTOSLayer.CORE,
            domain=PTOSDomain.DATA,
            tier=PTOSTier.FOAM,
            condition=PTOSCondition.STABLE,
        )

    # Byte histogram
    hist = [0] * 256
    for b in data:
        hist[b] += 1

    # Entropy
    entropy = 0.0
    for c in hist:
        if c > 0:
            p = c / n
            entropy -= p * math.log2(p)

    # Mean and variance of byte values
    mean_b = sum(data) / n
    var_b = sum((b - mean_b) ** 2 for b in data) / n
    std_b = math.sqrt(var_b)

    # Repetition: fraction of bytes that repeat the previous byte
    reps = sum(1 for i in range(1, n) if data[i] == data[i - 1])
    rep_ratio = reps / max(n - 1, 1)

    # Unique byte ratio
    unique = sum(1 for c in hist if c > 0)
    unique_ratio = unique / 256

    # ── Layer: structural depth ← entropy ──
    if entropy < 1.0:
        layer = PTOSLayer.CORE
    elif entropy < 3.0:
        layer = PTOSLayer.CARRY
    elif entropy < 5.0:
        layer = PTOSLayer.RULE
    elif entropy < 7.0:
        layer = PTOSLayer.STORE
    else:
        layer = PTOSLayer.EXTERNAL

    # ── Domain: byte distribution moments ──
    skewness = sum((b - mean_b) ** 3 for b in data) / (n * std_b ** 3) if std_b > 0 else 0
    if abs(skewness) < 0.1 and entropy > 6.0:
        domain = PTOSDomain.COMPUTE
    elif unique_ratio > 0.8:
        domain = PTOSDomain.TOKEN
    elif rep_ratio > 0.5:
        domain = PTOSDomain.STORE
    elif abs(skewness) > 2.0:
        domain = PTOSDomain.POWER
    elif entropy > 5.0:
        domain = PTOSDomain.DATA
    else:
        domain = PTOSDomain.RULE

    # ── Tier: complexity ← repetition ──
    if rep_ratio > 0.9:
        tier = PTOSTier.SINGULARITY
    elif rep_ratio > 0.7:
        tier = PTOSTier.PLASMA
    elif rep_ratio > 0.4:
        tier = PTOSTier.CRYSTALLINE
    elif rep_ratio > 0.2:
        tier = PTOSTier.FOAM
    elif entropy > 6.0:
        tier = PTOSTier.RESEARCH
    else:
        tier = PTOSTier.GOVERNANCE

    # ── Condition: stability ← deviation from uniform distribution ──
    uniform_expected = n / 256
    chi_sq = sum((c - uniform_expected) ** 2 / max(uniform_expected, 1) for c in hist)
    normalized_chi = chi_sq / max(n, 1)

    if normalized_chi < 0.5:
        condition = PTOSCondition.STABLE
    elif normalized_chi < 2.0:
        condition = PTOSCondition.EXPERIMENTAL
    elif normalized_chi < 5.0:
        condition = PTOSCondition.EXTREME
    elif n < 64:
        condition = PTOSCondition.DRAFT
    elif unique_ratio < 0.05:
        condition = PTOSCondition.STERILE
    else:
        condition = PTOSCondition.ARCHIVED

    return PTOSManifest(
        layer=layer,
        domain=domain,
        tier=tier,
        condition=condition,
    )


# ── §2  Metadata: Block Statistics (Δ surface) ─────────────────────────────

@dataclass
class BlockMetadata:
    """Per-block statistical metadata for compression planning."""
    offset: int
    size: int
    entropy: float
    mean: float
    std: float
    repetition_ratio: float
    unique_ratio: float
    compressibility_estimate: float  # 0 = incompressible, 1 = highly compressible


def extract_block_metadata(data: bytes, block_size: int = 4096) -> List[BlockMetadata]:
    """Extract per-block metadata for compression surface planning.

    Each block's statistical profile determines how it packs on the Δ surface.
    """
    blocks: List[BlockMetadata] = []
    for offset in range(0, len(data), block_size):
        chunk = data[offset:offset + block_size]
        n = len(chunk)
        if n == 0:
            continue

        hist = [0] * 256
        for b in chunk:
            hist[b] += 1

        entropy = 0.0
        for c in hist:
            if c > 0:
                p = c / n
                entropy -= p * math.log2(p)

        mean_b = sum(chunk) / n
        var_b = sum((b - mean_b) ** 2 for b in chunk) / n
        reps = sum(1 for i in range(1, n) if chunk[i] == chunk[i - 1])
        unique = sum(1 for c in hist if c > 0)

        # Compressibility estimate based on entropy and repetition
        # Lower entropy + higher repetition = more compressible
        max_entropy = 8.0
        entropy_factor = 1.0 - (entropy / max_entropy)
        rep_factor = reps / max(n - 1, 1)
        compressibility = 0.5 * entropy_factor + 0.5 * rep_factor

        blocks.append(BlockMetadata(
            offset=offset,
            size=n,
            entropy=entropy,
            mean=mean_b,
            std=math.sqrt(var_b),
            repetition_ratio=reps / max(n - 1, 1),
            unique_ratio=unique / 256,
            compressibility_estimate=compressibility,
        ))

    return blocks


# ── §3  Metadata: Spectral Features (Φ surface) ────────────────────────────
# Golden-angle sampling per GoldenAngleEncoding.lean

@dataclass
class WaveProbeSample:
    """A single golden-angle sample from the Φ surface.

    Matches GoldenAngleEncoding.lean: structure WaveProbeSample
    """
    index: int
    phase_q16: int        # Q16_16 golden-angle phase
    theta_q16: int        # Q16_16 spherical theta
    phi_q16: int           # Q16_16 spherical phi
    amplitude: int         # Raw amplitude at sample point
    timestamp: int = 0


@dataclass
class SpectralFeatures:
    """Extracted spectral metadata from WaveProbe sampling.

    These features describe the Φ-surface geometry.
    """
    samples: List[WaveProbeSample]
    mean_phase: float
    phase_variance: float
    mean_amplitude: float
    amplitude_variance: float
    spectral_energy: float       # Σ amplitude²
    spectral_centroid: float     # Energy-weighted mean phase
    spectral_spread: float       # Energy-weighted variance
    coherence_q16: int           # Q16_16 phase coherence [0, 65536]


class WaveProbe:
    """Φ-surface probe: golden-angle sampler of spectral features.

    Matches WebRTCWaveformSync.lean: structure WaveProbe
    and gccl_waveprobe.py: class WaveProbe
    """
    def __init__(self, probe_id: int = 1, sample_rate: int = 48000, buffer_size: int = 256):
        self.id = probe_id
        self.sample_rate = sample_rate
        self.buffer_size = buffer_size

    def sample(self, data: List[int]) -> List[int]:
        """Golden-angle stepping over data indices.

        From GoldenAngleEncoding.lean: goldenAngleStep = 40503
        φ⁻¹ ≈ 0.618 in Q16_16 ensures maximum coverage with minimum overlap.
        """
        n = len(data)
        if n == 0:
            return []
        step = max(1, (GOLDEN_ANGLE_STEP * n) // Q16_SCALE)
        pos = 0
        result = []
        for _ in range(min(self.buffer_size, n)):
            result.append(data[pos % n])
            pos = (pos + step) % n
        return result

    def signal_overlap(self, data_a: List[int], data_b: List[int]) -> int:
        """Q16_16 cosine similarity: 0 = no overlap, 65536 = identical.

        Matches Waveprobe.lean: overlap = |⟨ψc|ψp⟩|²
        """
        if not data_a or not data_b:
            return 0

        samp_a = self.sample(data_a)
        samp_b = self.sample(data_b)
        m = min(len(samp_a), len(samp_b))
        if m == 0:
            return 0

        mean_a = sum(samp_a[:m]) // m
        mean_b = sum(samp_b[:m]) // m

        dot = 0
        nrm_a = 0
        nrm_b = 0
        for i in range(m):
            da = samp_a[i] - mean_a
            db = samp_b[i] - mean_b
            dot += da * db
            nrm_a += da * da
            nrm_b += db * db

        denom = nrm_a * nrm_b
        if denom == 0:
            return Q16_SCALE if nrm_a == 0 and nrm_b == 0 else 0

        return max(0, min(Q16_SCALE, (dot * dot * Q16_SCALE) // denom))


def extract_spectral_features(data: bytes, probe: Optional[WaveProbe] = None) -> SpectralFeatures:
    """Extract Φ-surface metadata via WaveProbe golden-angle sampling.

    Returns spectral features including phase coherence and spectral energy.
    """
    if probe is None:
        probe = WaveProbe()

    data_ints = list(data)
    n = len(data_ints)
    if n == 0:
        return SpectralFeatures(
            samples=[], mean_phase=0, phase_variance=0,
            mean_amplitude=0, amplitude_variance=0,
            spectral_energy=0, spectral_centroid=0, spectral_spread=0,
            coherence_q16=Q16_SCALE,
        )

    # Golden-angle samples
    step = max(1, (GOLDEN_ANGLE_STEP * n) // Q16_SCALE)
    pos = 0
    samples: List[WaveProbeSample] = []
    amplitudes: List[float] = []
    phases: List[float] = []

    for idx in range(min(probe.buffer_size, n)):
        amp = float(data_ints[pos % n])
        phase_turns = (pos * GOLDEN_ANGLE_STEP % Q16_SCALE) / Q16_SCALE
        theta = (pos * Q16_SCALE // max(n, 1)) / Q16_SCALE

        amp_q16 = q16f(amp)
        phase_q16 = q16f(phase_turns)
        theta_q16 = q16f(theta)
        phi_q16 = q16f((pos * GOLDEN_ANGLE_STEP * Q16_SCALE // Q16_SCALE) % Q16_SCALE / Q16_SCALE)

        samples.append(WaveProbeSample(
            index=idx,
            phase_q16=phase_q16,
            theta_q16=theta_q16,
            phi_q16=phi_q16,
            amplitude=amp_q16,
        ))
        amplitudes.append(amp)
        phases.append(phase_turns * 2 * math.pi)
        pos = (pos + step) % n

    # Spectral statistics
    mean_amp = statistics.mean(amplitudes) if amplitudes else 0
    var_amp = statistics.variance(amplitudes) if len(amplitudes) > 1 else 0
    mean_phase_v = statistics.mean(phases) if phases else 0
    var_phase = statistics.variance(phases) if len(phases) > 1 else 0
    spectral_energy = sum(a * a for a in amplitudes)

    # Spectral centroid: energy-weighted mean phase
    total_energy = spectral_energy if spectral_energy > 0 else 1
    spectral_centroid = sum(amplitudes[i] * amplitudes[i] * phases[i]
                           for i in range(len(amplitudes))) / total_energy if amplitudes else 0
    spectral_spread = sum(amplitudes[i] * amplitudes[i] *
                          (phases[i] - spectral_centroid) ** 2
                          for i in range(len(amplitudes))) / total_energy if amplitudes else 0

    # Phase coherence: Rayleigh test for phase uniformity
    # R = |Σ exp(i·φ)| / N  — 1.0 = perfectly coherent, 0.0 = uniform noise
    complex_sum = complex(0, 0)
    for p in phases:
        complex_sum += cmath.exp(complex(0, p))
    r_value = abs(complex_sum) / max(len(phases), 1)
    coherence_q16 = max(0, min(Q16_SCALE, q16f(r_value)))

    return SpectralFeatures(
        samples=samples,
        mean_phase=mean_phase_v,
        phase_variance=var_phase,
        mean_amplitude=mean_amp,
        amplitude_variance=var_amp,
        spectral_energy=spectral_energy,
        spectral_centroid=spectral_centroid,
        spectral_spread=spectral_spread,
        coherence_q16=coherence_q16,
    )


# ── §4  Metadata: GCL Program Features (Γ surface) ─────────────────────────

GCL_CODONS = {
    'ATG': 'start', 'TAA': 'stop', 'TAG': 'stop', 'TGA': 'stop',
    'CTU': 'store', 'GCU': 'foam',
    'AAA': 'read', 'TTT': 'write', 'GGG': 'evolve', 'CCC': 'prune',
    'AGA': 'route', 'TCT': 'gate', 'GAG': 'bind', 'CTC': 'split',
    'ACA': 'merge', 'TGT': 'probe', 'GTA': 'emit', 'CAC': 'receipt',
}


@dataclass
class GCLProgramFeatures:
    """Γ-surface metadata: GCL program complexity and codon profile."""
    program_length: int
    codon_count: int
    unique_codons: int
    codon_frequencies: Dict[str, int]
    known_codon_ratio: float
    structural_depth: int
    complexity_estimate_q16: int   # Q16_16: 0 = trivial, 65536 = maximal
    mutation_count: int = 0
    generation: int = 0


def extract_gcl_program_features(program: str) -> GCLProgramFeatures:
    """Extract Γ-surface metadata from a GCL program string.

    Analyzes codon distribution, structural depth, and complexity.
    """
    # Count codons (3-char tokens, whitespace-separated)
    tokens = program.split()
    codons = [t.upper() for t in tokens if len(t) == 3 and all(c in 'ACGTU' for c in t.upper())]

    codon_freq: Dict[str, int] = {}
    for c in codons:
        codon_freq[c] = codon_freq.get(c, 0) + 1

    known_count = sum(v for k, v in codon_freq.items() if k in GCL_CODONS)
    total_codons = sum(codon_freq.values())
    known_ratio = known_count / max(total_codons, 1)

    # Structural depth: max nesting of braces/parens
    depth = 0
    max_depth = 0
    for ch in program:
        if ch in '({[':
            depth += 1
            max_depth = max(max_depth, depth)
        elif ch in ')}]':
            depth = max(0, depth - 1)

    # Complexity: blend of structural depth, codon variety, and known ratio
    unique_count = len(codon_freq)
    complexity_raw = (
        0.4 * (max_depth / 10.0) +
        0.3 * (unique_count / max(len(GCL_CODONS), 1)) +
        0.3 * (1.0 - known_ratio)
    )
    complexity_q16 = q16f(min(1.0, complexity_raw))

    return GCLProgramFeatures(
        program_length=len(program),
        codon_count=total_codons,
        unique_codons=unique_count,
        codon_frequencies=codon_freq,
        known_codon_ratio=known_ratio,
        structural_depth=max_depth,
        complexity_estimate_q16=complexity_q16,
    )


# ── §5  Surface Costs: The Multi-Surface Lagrangian ─────────────────────────

@dataclass
class SurfaceCosts:
    """The three surface costs in the multi-surface packing Lagrangian.

    Total:  L_total = ‖Δ‖₁ + α·|φ| + β·K(γ)
    """
    delta_cost: int           # ‖Δ(γ,x)‖₁ — delta residual byte length
    spectral_cost: int        # |φ| — WaveProbe sample count (α-weighted)
    program_cost: int         # K(γ) — program description length (β-weighted)
    alpha: int                # α scalar (Q16_16)
    beta: int                 # β scalar (Q16_16)
    total_lagrangian: int     # L_total (Q16_16)

    def to_dict(self) -> Dict:
        return {
            'delta_cost': self.delta_cost,
            'spectral_cost': self.spectral_cost,
            'program_cost': self.program_cost,
            'alpha': self.alpha,
            'beta': self.beta,
            'total_lagrangian': self.total_lagrangian,
            'total_lagrangian_float': self.total_lagrangian / Q16_SCALE,
        }


def compute_delta_surface_cost(data: bytes, reference: Optional[bytes],
                               program_features: GCLProgramFeatures) -> int:
    """‖Δ(γ,x)‖₁: delta residual byte length under program γ.

    Uses copy-if pattern (zero deltas are skipped). The GCL program's
    structural depth modulates the delta granularity.
    """
    if reference is None or len(reference) < len(data):
        return q16(len(data))  # No reference = full cost

    block_size = max(64, min(4096, 256 * (program_features.structural_depth + 1)))

    delta_count = 0
    for i in range(0, len(data), block_size):
        chunk = data[i:min(i + block_size, len(data))]
        ref_chunk = reference[i:min(i + block_size, len(reference))]

        for j in range(min(len(chunk), len(ref_chunk))):
            if chunk[j] != ref_chunk[j]:
                delta_count += 1

    return q16(delta_count)


def compute_spectral_surface_cost(probe_config: WaveProbe) -> int:
    """|φ|: measurement overhead cost.

    Proportional to buffer_size × sample_rate. Higher resolution = higher cost.
    """
    raw = probe_config.buffer_size * max(probe_config.sample_rate // 48000, 1)
    return q16(raw)


def compute_program_surface_cost(features: GCLProgramFeatures) -> int:
    """K(γ): program description length proxy.

    Based on codon count and structural complexity.
    """
    raw = features.codon_count * (features.structural_depth + 1) + features.program_length
    return q16(raw)


def compute_lagrangian(data: bytes, reference: Optional[bytes],
                       probe_config: WaveProbe, program_features: GCLProgramFeatures,
                       alpha_q16: int = q16f(0.1), beta_q16: int = q16f(0.05)) -> SurfaceCosts:
    """Compute the full multi-surface Lagrangian.

    L_total = ‖Δ‖₁ + α·|φ| + β·K(γ)
    """
    delta_cost = compute_delta_surface_cost(data, reference, program_features)
    spectral_cost = q16_mul(alpha_q16, compute_spectral_surface_cost(probe_config))
    program_cost = q16_mul(beta_q16, compute_program_surface_cost(program_features))

    total = q16_add(delta_cost, q16_add(spectral_cost, program_cost))

    return SurfaceCosts(
        delta_cost=delta_cost,
        spectral_cost=spectral_cost,
        program_cost=program_cost,
        alpha=alpha_q16,
        beta=beta_q16,
        total_lagrangian=total,
    )


# ── §6  GCCL Gates (admissibility on the product manifold) ─────────────────

class GCCLDecision(IntEnum):
    ACCEPT = 0
    REJECT = 1
    HOLD = 2
    QUARANTINE = 3


def gccl_swap_gate(old_cost_q16: int, new_cost_q16: int, recon_risk_q16: int) -> bool:
    """GCCL swap gate: accept if new cost < old cost and risk is bounded.

    Matches GCCL.lean: def gcclSwapGate
    """
    if old_cost_q16 > new_cost_q16:
        admissible = q16_sub(old_cost_q16, new_cost_q16)
    else:
        admissible = 0
    return admissible <= recon_risk_q16


def coherence_gate(overlap_q16: int, epsilon_q16: int) -> bool:
    """WaveProbe coherence gate: ⟨ψ_x|ψ_x̂⟩² ≥ 1 - ε_λ.

    The structure-preservation constraint on the Φ surface.
    """
    threshold = q16_sub(Q16_SCALE, epsilon_q16)
    return overlap_q16 >= threshold


# ── §7  Full Pipeline ──────────────────────────────────────────────────────

@dataclass
class PackingReceipt:
    """Complete multi-surface packing receipt.

    Records the Lagrangian costs, metadata, and gate decisions
    for every surface in the product manifold.
    """
    schema: str = "multi_surface_packing_v1"

    # Input provenance
    input_size: int = 0
    input_hash: str = ""

    # Metadata per surface
    ptos_manifest: Optional[Dict] = None
    block_metadata: List[Dict] = field(default_factory=list)
    spectral_features: Optional[Dict] = None
    program_features: Optional[Dict] = None

    # Lagrangian costs
    surface_costs: Optional[Dict] = None

    # Gate decisions
    coherence_passed: bool = False
    gccl_passed: bool = False
    overall_decision: str = "HOLD"

    # Compression result
    compressed_size: int = 0
    compression_ratio: float = 0.0
    compression_percentage: float = 0.0

    # Timing
    elapsed_ms: float = 0.0

    # Chain
    parent_hash: str = ""
    receipt_hash: str = ""

    def finalize(self) -> str:
        """Compute receipt hash and return JSON."""
        d = self.to_dict()
        raw = json.dumps(d, sort_keys=True, separators=(',', ':'))
        self.receipt_hash = hashlib.sha256(raw.encode()).hexdigest()[:16]
        d['receipt_hash'] = self.receipt_hash
        return json.dumps(d, indent=2)

    def to_dict(self) -> Dict:
        return {
            'schema': self.schema,
            'input_size': self.input_size,
            'input_hash': self.input_hash,
            'ptos_manifest': self.ptos_manifest or {},
            'block_metadata_count': len(self.block_metadata),
            'spectral_features': self.spectral_features or {},
            'program_features': self.program_features or {},
            'surface_costs': self.surface_costs or {},
            'coherence_passed': self.coherence_passed,
            'gccl_passed': self.gccl_passed,
            'overall_decision': self.overall_decision,
            'compressed_size': self.compressed_size,
            'compression_ratio': self.compression_ratio,
            'compression_percentage': self.compression_percentage,
            'elapsed_ms': self.elapsed_ms,
            'parent_hash': self.parent_hash,
            'receipt_hash': self.receipt_hash,
            'claim_boundary': 'multi-surface-packing-assembly-only',
        }


def multi_surface_pack(
    data: bytes,
    reference: Optional[bytes] = None,
    gcl_program: str = "",
    probe_id: int = 1,
    sample_rate: int = 48000,
    buffer_size: int = 256,
    alpha_q16: int = q16f(0.1),
    beta_q16: int = q16f(0.05),
    epsilon_q16: int = q16f(0.05),
    recon_risk_q16: int = q16f(0.5),
    parent_hash: str = "",
    block_size: int = 4096,
) -> Tuple[bytes, PackingReceipt]:
    """Run the full multi-surface packing pipeline.

    Pipeline:
      1. Extract metadata from all three surfaces (Δ, Φ, Γ)
      2. Compute Lagrangian costs
      3. Evaluate coherence and admissibility gates
      4. Apply delta+RLE encoding
      5. Generate receipt

    Args:
        data: Input bytes to compress.
        reference: Optional reference for delta encoding.
        gcl_program: GCCL program string for Γ-surface analysis.
        probe_id/sample_rate/buffer_size: WaveProbe Φ-surface configuration.
        alpha_q16: Φ-surface cost weight in the Lagrangian.
        beta_q16: Γ-surface cost weight in the Lagrangian.
        epsilon_q16: Coherence tolerance (1 - ε_λ threshold).
        recon_risk_q16: GCCL swap gate risk bound.
        parent_hash: Previous receipt hash for chain continuity.
        block_size: Δ-surface block size for metadata extraction.

    Returns:
        Tuple of (compressed_bytes, receipt).
    """
    t0 = time.time()

    # ── Step 1: Δ-surface metadata ──
    ptos = extract_ptos_manifest(data)
    block_meta = extract_block_metadata(data, block_size)

    # ── Step 2: Φ-surface metadata ──
    probe = WaveProbe(probe_id=probe_id, sample_rate=sample_rate, buffer_size=buffer_size)
    spectral = extract_spectral_features(data, probe=probe)

    # Signal overlap for coherence gate
    overlap_q16 = Q16_SCALE
    if reference:
        overlap_q16 = probe.signal_overlap(list(data), list(reference))

    # ── Step 3: Γ-surface metadata ──
    program_features = extract_gcl_program_features(gcl_program)

    # ── Step 4: Lagrangian costs ──
    costs = compute_lagrangian(
        data, reference, probe, program_features,
        alpha_q16=alpha_q16, beta_q16=beta_q16,
    )

    # ── Step 5: Gates ──
    coh_pass = coherence_gate(overlap_q16, epsilon_q16)
    gccl_pass = gccl_swap_gate(
        old_cost_q16=q16(len(data)),
        new_cost_q16=costs.delta_cost,
        recon_risk_q16=recon_risk_q16,
    )

    if coh_pass and gccl_pass:
        decision = "ACCEPT"
    elif not coh_pass:
        decision = "QUARANTINE"
    else:
        decision = "HOLD"

    # ── Step 6: Apply compression (delta + RLE) ──
    if reference and len(reference) >= len(data):
        compressed = bytes(
            (data[i] - reference[i]) & 0xFF
            for i in range(len(data))
        )
    else:
        compressed = data

    # Simple RLE pass
    if len(compressed) > 0:
        rle: List[int] = []
        run_val = compressed[0]
        run_len = 1
        for b in compressed[1:]:
            if b == run_val and run_len < 255:
                run_len += 1
            else:
                rle.append(run_len)
                rle.append(run_val)
                run_val = b
                run_len = 1
        rle.append(run_len)
        rle.append(run_val)
        compressed = bytes(rle)

    elapsed = (time.time() - t0) * 1000

    # ── Step 7: Build receipt ──
    receipt = PackingReceipt(
        input_size=len(data),
        input_hash=hashlib.sha256(data[:1024]).hexdigest()[:16],
        ptos_manifest=ptos.to_dict(),
        block_metadata=[{
            'offset': b.offset, 'size': b.size,
            'entropy': round(b.entropy, 4),
            'compressibility': round(b.compressibility_estimate, 4),
        } for b in block_meta[:16]],  # summary only
        spectral_features={
            'sample_count': len(spectral.samples),
            'coherence_q16': spectral.coherence_q16,
            'coherence': spectral.coherence_q16 / Q16_SCALE,
            'spectral_energy': round(spectral.spectral_energy, 2),
            'mean_amplitude': round(spectral.mean_amplitude, 4),
        },
        program_features={
            'codon_count': program_features.codon_count,
            'unique_codons': program_features.unique_codons,
            'structural_depth': program_features.structural_depth,
            'complexity_q16': program_features.complexity_estimate_q16,
            'known_codon_ratio': round(program_features.known_codon_ratio, 4),
        },
        surface_costs=costs.to_dict(),
        coherence_passed=coh_pass,
        gccl_passed=gccl_pass,
        overall_decision=decision,
        compressed_size=len(compressed),
        compression_ratio=len(data) / max(len(compressed), 1),
        compression_percentage=(len(data) - len(compressed)) / max(len(data), 1) * 100,
        elapsed_ms=round(elapsed, 2),
        parent_hash=parent_hash,
    )

    return compressed, receipt


# ── §8  CLI ─────────────────────────────────────────────────────────────────

if __name__ == '__main__':
    import sys
    import cmath  # needed for spectral phase coherence

    if len(sys.argv) < 2:
        print(__doc__)
        print("Usage: python multi_surface_packer.py <file> [--reference <ref>] [--alpha 0.1] [--beta 0.05]")
        sys.exit(1)

    filepath = sys.argv[1]
    ref_path = None
    alpha = 0.1
    beta = 0.05
    buffer_size = 256
    gcl_program = ""

    i = 2
    while i < len(sys.argv):
        if sys.argv[i] == '--reference' and i + 1 < len(sys.argv):
            ref_path = sys.argv[i + 1]
            i += 2
        elif sys.argv[i] == '--alpha' and i + 1 < len(sys.argv):
            alpha = float(sys.argv[i + 1])
            i += 2
        elif sys.argv[i] == '--beta' and i + 1 < len(sys.argv):
            beta = float(sys.argv[i + 1])
            i += 2
        elif sys.argv[i] == '--buffer' and i + 1 < len(sys.argv):
            buffer_size = int(sys.argv[i + 1])
            i += 2
        elif sys.argv[i] == '--program' and i + 1 < len(sys.argv):
            gcl_program = sys.argv[i + 1]
            i += 2
        else:
            i += 1

    with open(filepath, 'rb') as f:
        data = f.read()

    ref_data = None
    if ref_path:
        with open(ref_path, 'rb') as f:
            ref_data = f.read()

    print(f"Input:  {filepath} ({len(data)} bytes)")
    if ref_data:
        print(f"Ref:    {ref_path} ({len(ref_data)} bytes)")
    print(f"α={alpha}, β={beta}, buffer={buffer_size}")

    compressed, receipt = multi_surface_pack(
        data,
        reference=ref_data,
        gcl_program=gcl_program,
        buffer_size=buffer_size,
        alpha_q16=q16f(alpha),
        beta_q16=q16f(beta),
    )

    receipt_json = receipt.finalize()
    print(f"\nReceipt:")
    print(receipt_json)
    print(f"\nCompressed: {len(compressed)} bytes (ratio: {receipt.compression_ratio:.2f}x)")
    print(f"Decision:   {receipt.overall_decision}")
    print(f"Lagrangian: {receipt.surface_costs['total_lagrangian_float']:.4f}")