#!/usr/bin/env python3
# /// script
# requires-python = ">=3.10"
# dependencies = []
# ///
"""
spectral_codebook.py — Gap-aware spectral codebook over the PIST matrix corpus.

Implements the "Next Steps" of docs/SPECTRAL_CODEBOOK_ANALYSIS.md, with
corrections found while building it (see docs/SPECTRAL_CODEBOOK_GENERATOR.md):

  1. Parse the 250 8×8 integer braid adjacency matrices from
     formal/SilverSight/PIST/Matrices250.lean (regex on the auto-generated
     `def <name> : Array (Array Int)` blocks, plus the `findMatrix` cases).
  2. Compute a full spectral profile per matrix:
       - EXACT integer characteristic-polynomial coefficients via the
         Faddeev–LeVerrier algorithm in fractions.Fraction arithmetic
         (the primary, float-free fingerprint),
       - spectral radius λ as the max root modulus of that exact char poly
         (NumPy eigvals fast path; pure-stdlib Durand–Kerner fallback),
       - power-iteration λ retained only for traceability: the analysis
         doc's 9 "mid band" values in (0.5, 1) are power-iteration
         NON-CONVERGENCE artifacts — all 9 matrices have exact ρ = 1.0
         (peripheral spectrum with several eigenvalues of modulus 1).
  3. Gap-aware quantization with a sample-size guard: dedupe byte-identical
     matrices first (the corpus has 238 distinct matrices over 250 ids),
     find gaps > 3× median gap between consecutive distinct λ, then merge
     any segment with fewer than MIN_SUPPORT distinct matrices into its
     neighbor across the smaller gap. Points isolated above a suppressed
     gap with sub-threshold occupancy are flagged as sparse-tail outliers
     rather than treated as clusters.
  4. Round-trip: Codebook.encode(matrix) → (codeword, index) and
     Codebook.decode(codeword, index) → equation_id are mutually inverse
     over the corpus. Neither λ nor the char poly is injective (nilpotent
     x^8 class alone has 31 distinct members), so the decode key is the
     within-cluster rank index; all collision classes are explicit in the
     emitted JSON.

Pure stdlib. NumPy is used only as an optional fast path when importable;
the pure-Python path is authoritative.

Usage:
    python3 python/spectral_codebook.py                   # build + report
    python3 python/spectral_codebook.py --check-manifold  # incl. 278-row corpus
    python3 python/spectral_codebook.py --out data/spectral_codebook.json
"""
from __future__ import annotations

import argparse
import json
import math
import re
import sys
from datetime import date
from fractions import Fraction
from pathlib import Path
from typing import Dict, List, Optional, Sequence, Tuple

try:  # optional fast path only; both paths agree to float rounding
    import numpy as _np
except ImportError:  # pragma: no cover
    _np = None

REPO_ROOT = Path(__file__).resolve().parent.parent
MATRICES_LEAN = REPO_ROOT / "formal" / "SilverSight" / "PIST" / "Matrices250.lean"
MANIFOLD_LEAN = REPO_ROOT / "formal" / "SilverSight" / "RRC" / "Q16_16Manifold.lean"
DEFAULT_OUT = REPO_ROOT / "data" / "spectral_codebook.json"

POWER_ITERS = 300           # matches docs/SPECTRAL_CODEBOOK_ANALYSIS.md
GAP_FACTOR = 3.0            # boundary threshold = GAP_FACTOR × median gap
MIN_SUPPORT = 10            # min distinct matrices per cluster (sample-size guard)
LAMBDA_DP = 4               # rounding used for "unique λ" counting (as in doc)

Matrix = Tuple[Tuple[int, ...], ...]


# ──────────────────────────────────────────────────────────────────────
# Parsing Matrices250.lean
# ──────────────────────────────────────────────────────────────────────

_DEF_RE = re.compile(
    r"^def\s+(\w+)\s*:\s*Array\s*\(Array\s+Int\)\s*:=\s*\n\s*#\[(.*?)\n\s*\]",
    re.MULTILINE | re.DOTALL,
)
_ROW_RE = re.compile(r"#\[([^\]]*)\]")
_FIND_CASE_RE = re.compile(r'\|\s*"([^"]+)"\s*=>\s*some\s+(\w+)')


def parse_matrices_lean(path: Path = MATRICES_LEAN) -> Dict[str, Matrix]:
    """Parse `def <name> : Array (Array Int) := #[#[...], ...]` blocks.

    Returns an equation_id → matrix mapping (row-major tuples of ints),
    resolving ids through the `findMatrix` match cases when present so the
    mapping survives Lean-identifier sanitization of ids.
    """
    text = path.read_text()
    defs: Dict[str, Matrix] = {}
    for name, body in _DEF_RE.findall(text):
        rows = tuple(
            tuple(int(x) for x in row.split(",") if x.strip())
            for row in _ROW_RE.findall(body)
        )
        if rows:
            defs[name] = rows

    cases = _FIND_CASE_RE.findall(text)
    if cases:
        mats = {eid: defs[name] for eid, name in cases if name in defs}
    else:  # fall back: def name == equation id
        mats = dict(defs)

    for eid, mat in mats.items():
        n = len(mat)
        if any(len(r) != n for r in mat):
            raise ValueError(f"non-square matrix for {eid}")
    return mats


def parse_manifold_ids(path: Path = MANIFOLD_LEAN) -> List[str]:
    """All equationId occurrences (with multiplicity) in the fixture corpus."""
    text = path.read_text()
    return re.findall(r'equationId\s*:=\s*"([^"]+)"', text)


# ──────────────────────────────────────────────────────────────────────
# Spectral profile
# ──────────────────────────────────────────────────────────────────────

def charpoly_coeffs(mat: Sequence[Sequence[int]]) -> Tuple[int, ...]:
    """Exact characteristic-polynomial coefficients via Faddeev–LeVerrier.

    For an n×n integer matrix M returns (c1, ..., cn) with

        det(λI − M) = λ^n + c1·λ^(n−1) + ... + cn.

    Runs in exact fractions.Fraction arithmetic; every ck is provably an
    integer (char poly of an integer matrix), which is asserted.
    """
    n = len(mat)
    m = [[Fraction(x) for x in row] for row in mat]
    aux = [[Fraction(0)] * n for _ in range(n)]   # M_0 = 0
    coeffs: List[int] = []
    c = Fraction(1)
    for k in range(1, n + 1):
        # M_k = M · (M_{k−1} + c_{k−1}·I)
        shifted = [
            [aux[i][j] + (c if i == j else 0) for j in range(n)]
            for i in range(n)
        ]
        aux = [
            [sum(m[i][t] * shifted[t][j] for t in range(n)) for j in range(n)]
            for i in range(n)
        ]
        c = -sum(aux[i][i] for i in range(n)) / k
        assert c.denominator == 1, "char-poly coefficient must be integral"
        coeffs.append(int(c))
    return tuple(coeffs)


def _durand_kerner_max_modulus(coeffs: Sequence[int]) -> float:
    """Max root modulus of the monic polynomial x^n + c1 x^(n-1) + ... + cn.

    Pure-stdlib root finder (Durand–Kerner). Zero roots are stripped
    exactly first, so nilpotent classes return 0.0 with no iteration.
    """
    poly = [1] + [int(c) for c in coeffs]
    # strip exact zero roots (trailing zero coefficients)
    while len(poly) > 1 and poly[-1] == 0:
        poly.pop()
    n = len(poly) - 1
    if n == 0:
        return 0.0
    # Cauchy bound and standard (0.4 + 0.9i)^k seeds
    bound = 1.0 + max(abs(c) for c in poly[1:])
    roots = [complex(0.4, 0.9) ** k * (bound / 1.1) for k in range(n)]

    def ev(z: complex) -> complex:
        r = complex(0)
        for c in poly:
            r = r * z + c
        return r

    for _ in range(1000):
        delta = 0.0
        new_roots = []
        for i, z in enumerate(roots):
            den = complex(1)
            for j, w in enumerate(roots):
                if j != i:
                    den *= (z - w)
            dz = ev(z) / den if den != 0 else complex(0)
            new_roots.append(z - dz)
            delta = max(delta, abs(dz))
        roots = new_roots
        if delta < 1e-13:
            break
    return max(abs(z) for z in roots)


def spectral_radius(mat: Sequence[Sequence[int]],
                    coeffs: Optional[Sequence[int]] = None) -> float:
    """Exact-spectrum spectral radius ρ(M) = max |eigenvalue|.

    NumPy eigvals fast path when available; otherwise Durand–Kerner on the
    exact integer char poly. Do NOT use power iteration here: it fails to
    converge on peripheral spectra (several eigenvalues of equal modulus),
    which is exactly what the ρ = 1 matrices in this corpus have.
    """
    if _np is not None:
        m = _np.asarray(mat, dtype=float)
        if m.size == 0:
            return 0.0
        return float(max(abs(_np.linalg.eigvals(m))))
    if coeffs is None:
        coeffs = charpoly_coeffs(mat)
    return _durand_kerner_max_modulus(coeffs)


def power_iteration(mat: Sequence[Sequence[int]], iters: int = POWER_ITERS) -> float:
    """Power-iteration λ estimate — UNTRUSTED, kept only for traceability
    against data/spectral_codebook_raw.json. Does not converge for
    peripheral spectra (see spectral_radius)."""
    n = len(mat)
    if n == 0:
        return 0.0
    v = [1.0 / math.sqrt(n)] * n
    lam = 0.0
    for _ in range(iters):
        w = [sum(mat[i][j] * v[j] for j in range(n)) for i in range(n)]
        norm = math.sqrt(sum(x * x for x in w))
        if norm < 1e-12:
            return 0.0
        lam = norm
        v = [x / norm for x in w]
    return lam


def density(mat: Sequence[Sequence[int]]) -> float:
    """Total entry weight / n² (matches SpectralN.lean and the raw JSON)."""
    n = len(mat)
    return sum(sum(row) for row in mat) / (n * n) if n else 0.0


def frobenius(mat: Sequence[Sequence[int]]) -> float:
    return math.sqrt(sum(x * x for row in mat for x in row))


def spectral_profile(mat: Sequence[Sequence[int]]) -> dict:
    coeffs = charpoly_coeffs(mat)
    return {
        "spectral_radius": round(spectral_radius(mat, coeffs), 6),
        "lambda_power_iteration": round(power_iteration(mat), 6),
        "charpoly": list(coeffs),
        "density": round(density(mat), 4),
        "frobenius": round(frobenius(mat), 4),
        "trace": sum(mat[i][i] for i in range(len(mat))),
    }


# ──────────────────────────────────────────────────────────────────────
# Φ-corkscrew: integer spiral-index packing (exact) + f(n) layout
# ──────────────────────────────────────────────────────────────────────
#
# The packing layer is EXACT and provable: with a corpus-wide base
#     B = 2·max|c| + 1   and offset  o = max|c|,
# every offset digit c_i + o lies in [0, B), so
#     n = Σ (c_i + o)·B^i
# is the base-B positional numeral of the digit tuple — injective on
# char-poly fingerprints by uniqueness of positional representation, and
# unpack_spiral_index inverts it exactly. It inherits charpoly's
# NON-injectivity on matrices (cospectral classes share one index).
#
# The layout layer f(n) = (√n·cos(nψ), √n·sin(nψ)), ψ = 2π/φ², is
# DECORATIVE placement only: injective because 1/φ² is irrational
# (distinct n give distinct angles), information-neutral (everything is
# recoverable from n alone), low-discrepancy by the three-distance
# theorem. It carries no information beyond spiral_index.

def packing_params(fingerprints: Sequence[Sequence[int]]) -> Tuple[int, int]:
    """Corpus-wide (base, offset) for spiral-index packing.

    offset = max|c| over every coefficient of every fingerprint;
    base = 2·offset + 1, so offset-encoded digits fill [0, base).
    """
    offset = max((abs(c) for fp in fingerprints for c in fp), default=0)
    return 2 * offset + 1, offset


def pack_spiral_index(coeffs: Sequence[int], base: int, offset: int) -> int:
    """Little-endian base-B packing of offset-encoded coefficients.

    Exact integer arithmetic; injective for digits in range (raises
    otherwise). Replaces the float phinary packing of
    integration_sprint.phi_corkscrew_index, which was not injective.
    """
    n = 0
    for i, c in enumerate(coeffs):
        d = c + offset
        if not 0 <= d < base:
            raise ValueError(
                f"coefficient {c} out of packing range for "
                f"base={base}, offset={offset}")
        n += d * base ** i
    return n


def unpack_spiral_index(n: int, base: int, offset: int,
                        length: int = 8) -> Tuple[int, ...]:
    """Exact inverse of pack_spiral_index."""
    if n < 0:
        raise ValueError("spiral index must be non-negative")
    out: List[int] = []
    for _ in range(length):
        n, d = divmod(n, base)
        out.append(d - offset)
    if n:
        raise ValueError("spiral index exceeds length digits in this base")
    return tuple(out)


# 1/φ² = (3 − √5)/2: the fractional rotation per step of the corkscrew.
def angle_frac(n: int) -> float:
    """frac(n/φ²) — the corkscrew angle of index n as a fraction of 2π.

    Computed in decimal arithmetic at a precision that covers the integer
    part of n·(3−√5)/2, because float64 loses the fractional part
    entirely once n exceeds 2^52 (spiral indices here reach base^8).
    """
    from decimal import Decimal, getcontext
    getcontext().prec = len(str(abs(n))) + 25
    alpha = (Decimal(3) - Decimal(5).sqrt()) / 2
    return float((Decimal(n) * alpha) % 1)


def corkscrew_xy(n: int) -> Tuple[float, float]:
    """f(n) = (√n·cos(nψ), √n·sin(nψ)), ψ = 2π/φ². Placement only."""
    theta = 2.0 * math.pi * angle_frac(n)
    r = math.sqrt(n)
    return r * math.cos(theta), r * math.sin(theta)


def min_angular_separation(indices: Sequence[int]) -> Optional[float]:
    """Minimum circular distance (as a fraction of a full turn) between
    the corkscrew angles of distinct spiral indices."""
    uniq = sorted({angle_frac(n) for n in set(indices)})
    if len(uniq) < 2:
        return None
    gaps = [uniq[i + 1] - uniq[i] for i in range(len(uniq) - 1)]
    gaps.append(1.0 - uniq[-1] + uniq[0])  # wrap-around
    return min(gaps)


# ──────────────────────────────────────────────────────────────────────
# Cartan ∆-floor quantization (Fisher metric on Δ₇)
# ──────────────────────────────────────────────────────────────────────
#
# ∆ = 17/1792 is the minimum eigenvalue of the Cartan crossing blocks,
# proven exact in formal/.../CartanConnection.lean:70 (see
# docs/cartan_fingerprint.md). It is adopted here as the resolution
# floor on the Fisher metric: two codewords closer than ∆ are treated
# as indistinguishable at operator resolution. This gives a PRINCIPLED
# alternative to the 3×-median-gap heuristic; both are emitted so the
# rules can be compared, neither replaces the other.

CARTAN_DELTA = Fraction(17, 1792)  # exact; CartanConnection.lean:70


def simplex_project(coeffs: Sequence[int]) -> Tuple[float, ...]:
    """Project a fingerprint onto Δ₇: p_i = |c_i| / Σ|c_j|.

    Convention: the all-zero fingerprint (nilpotent x⁸ class) maps to
    the uniform distribution — flagged in the report, since that choice
    is a convention, not a theorem.
    """
    total = sum(abs(c) for c in coeffs)
    k = len(coeffs)
    if total == 0:
        return tuple(1.0 / k for _ in range(k))
    return tuple(abs(c) / total for c in coeffs)


def fisher_distance(p: Sequence[float], q: Sequence[float]) -> float:
    """d_F(p, q) = 2·arccos(Σ√(p_i·q_i)) — Fisher geodesic distance on
    the probability simplex (Bhattacharyya angle × 2)."""
    bc = sum(math.sqrt(a * b) for a, b in zip(p, q))
    return 2.0 * math.acos(min(1.0, bc))


def cartan_floor_report(profiles: Dict[str, dict]) -> dict:
    """All-pairs Fisher distances between distinct fingerprints; pairs
    below ∆ = 17/1792 are merged into ∆-resolution components (union-
    find), giving the codeword count the Cartan floor supports."""
    fp_rep: Dict[Tuple[int, ...], str] = {}
    for eid in sorted(profiles):
        fp_rep.setdefault(tuple(profiles[eid]["charpoly"]), eid)
    fps = list(fp_rep)
    pts = [simplex_project(fp) for fp in fps]
    delta = float(CARTAN_DELTA)

    parent = list(range(len(fps)))

    def find(i: int) -> int:
        while parent[i] != i:
            parent[i] = parent[parent[i]]
            i = parent[i]
        return i

    sub_delta: List[dict] = []
    min_d, min_pair = None, None
    for i in range(len(fps)):
        for j in range(i + 1, len(fps)):
            d = fisher_distance(pts[i], pts[j])
            if min_d is None or d < min_d:
                min_d, min_pair = d, (fp_rep[fps[i]], fp_rep[fps[j]])
            if d < delta:
                sub_delta.append({
                    "representatives": [fp_rep[fps[i]], fp_rep[fps[j]]],
                    "fisher_distance": round(d, 8),
                })
                ri, rj = find(i), find(j)
                if ri != rj:
                    parent[rj] = ri
    components = len({find(i) for i in range(len(fps))})
    return {
        "delta_exact": "17/1792",
        "delta_source": "CartanConnection.lean:70 (docs/cartan_fingerprint.md)",
        "delta_float": delta,
        "metric": "Fisher on Δ₇: d_F = 2·arccos(Σ√(p_i·q_i)); "
                  "p_i = |c_i|/Σ|c_j| (all-zero fingerprint → uniform, "
                  "by convention)",
        "distinct_fingerprints": len(fps),
        "pairs_below_delta": len(sub_delta),
        "sub_delta_pairs": sub_delta[:50],
        "delta_resolution_codewords": components,
        "min_fisher_distance": round(min_d, 8) if min_d is not None else None,
        "min_fisher_pair": list(min_pair) if min_pair else None,
    }


# ──────────────────────────────────────────────────────────────────────
# Gap-aware quantization (deduped, sample-size guarded)
# ──────────────────────────────────────────────────────────────────────

def gap_quantize(lambda_counts: Dict[float, int],
                 factor: float = GAP_FACTOR,
                 min_support: int = MIN_SUPPORT) -> dict:
    """Gap-aware boundaries over distinct-λ values with occupancy counts.

    `lambda_counts` maps λ → number of DISTINCT matrices at that λ (dedupe
    byte-identical matrices before calling; 12 of the 250 corpus matrices
    are exact duplicates of another).

    Candidate boundaries sit at midpoints of gaps > factor × median gap.
    Sample-size guard: segments holding < min_support distinct matrices are
    merged into the neighbor across the smaller bounding gap, so no cluster
    is asserted on a handful of points. The suppressed high-end boundaries
    define a sparse-tail cutoff: everything above the lowest suppressed
    boundary with < min_support mass above it is an outlier region, not a
    cluster.
    """
    uniq = sorted(lambda_counts)
    if len(uniq) < 2:
        return {"boundaries": [], "median_gap": 0.0, "threshold": 0.0,
                "candidates": [], "suppressed": [], "sparse_tail_cutoff": None}
    gaps = [uniq[i + 1] - uniq[i] for i in range(len(uniq) - 1)]
    sizes = sorted(gaps)
    mid = len(sizes) // 2
    median = sizes[mid] if len(sizes) % 2 else 0.5 * (sizes[mid - 1] + sizes[mid])
    threshold = factor * median
    # candidate boundary index i ↔ gap between uniq[i] and uniq[i+1]
    candidates = [i for i, g in enumerate(gaps) if g > threshold]

    def seg_counts(bnds: List[int]) -> List[int]:
        counts, prev = [], 0
        for b in bnds + [len(uniq) - 1]:
            counts.append(sum(lambda_counts[uniq[k]] for k in range(prev, b + 1)))
            prev = b + 1
        return counts

    kept = list(candidates)
    suppressed: List[int] = []
    while kept:
        counts = seg_counts(kept)
        weak = [s for s, c in enumerate(counts) if c < min_support]
        if not weak:
            break
        s = min(weak, key=lambda i: counts[i])
        # bounding boundaries of segment s in `kept`
        options = []
        if s > 0:
            options.append(kept[s - 1])
        if s < len(kept):
            options.append(kept[s])
        drop = min(options, key=lambda b: gaps[b])  # merge across smaller gap
        kept.remove(drop)
        suppressed.append(drop)

    def midpoint(i: int) -> float:
        return 0.5 * (uniq[i] + uniq[i + 1])

    total = sum(lambda_counts.values())
    sparse_cutoff: Optional[float] = None
    for b in sorted(suppressed):
        above = sum(lambda_counts[v] for v in uniq if v > midpoint(b))
        if above < min_support:
            sparse_cutoff = midpoint(b)
            break

    return {
        "boundaries": [round(midpoint(b), 6) for b in sorted(kept)],
        "median_gap": median,
        "threshold": threshold,
        "candidates": [round(midpoint(b), 6) for b in sorted(candidates)],
        "suppressed": [round(midpoint(b), 6) for b in sorted(suppressed)],
        "sparse_tail_cutoff": sparse_cutoff,
        "total_support": total,
    }


class Codebook:
    """Gap-aware spectral codebook: matrix → (codeword, index) → equation_id."""

    def __init__(self, matrices: Dict[str, Matrix],
                 gap_factor: float = GAP_FACTOR,
                 min_support: int = MIN_SUPPORT) -> None:
        self.matrices = matrices
        self.gap_factor = gap_factor
        self.min_support = min_support
        self.profiles: Dict[str, dict] = {
            eid: spectral_profile(mat) for eid, mat in matrices.items()
        }

        # Φ-corkscrew integer packing (corpus-wide base/offset, exact)
        self.pack_base, self.pack_offset = packing_params(
            [p["charpoly"] for p in self.profiles.values()])
        self.spiral_indices: Dict[str, int] = {
            eid: pack_spiral_index(p["charpoly"], self.pack_base, self.pack_offset)
            for eid, p in self.profiles.items()
        }

        # dedupe byte-identical matrices BEFORE gap analysis
        self.matrix_groups: Dict[Matrix, List[str]] = {}
        for eid in sorted(self.profiles):
            self.matrix_groups.setdefault(matrices[eid], []).append(eid)
        lambda_counts: Dict[float, int] = {}
        for mat, ids in self.matrix_groups.items():
            lam = self.profiles[ids[0]]["spectral_radius"]
            lambda_counts[lam] = lambda_counts.get(lam, 0) + 1

        self.quant = gap_quantize(lambda_counts, gap_factor, min_support)
        self.boundaries: List[float] = self.quant["boundaries"]

        # cluster assignment + deterministic within-cluster rank
        by_cluster: Dict[int, List[str]] = {}
        for eid, p in self.profiles.items():
            by_cluster.setdefault(self.cluster_of(p["spectral_radius"]), []).append(eid)
        self._assign: Dict[str, Tuple[str, int]] = {}
        self._decode: Dict[Tuple[str, int], str] = {}
        self._by_matrix: Dict[Matrix, Tuple[str, int]] = {}
        for ci in sorted(by_cluster):
            members = sorted(
                by_cluster[ci],
                key=lambda e: (
                    self.profiles[e]["spectral_radius"],
                    self.profiles[e]["charpoly"],
                    self.matrices[e],
                    e,
                ),
            )
            for idx, eid in enumerate(members):
                cw = f"C{ci}"
                self._assign[eid] = (cw, idx)
                self._decode[(cw, idx)] = eid
                # first (canonical) entry wins for byte-identical duplicates
                self._by_matrix.setdefault(self.matrices[eid], (cw, idx))

    def cluster_of(self, lam: float) -> int:
        return sum(1 for b in self.boundaries if lam > b)

    def is_sparse_tail(self, lam: float) -> bool:
        cut = self.quant["sparse_tail_cutoff"]
        return cut is not None and lam > cut

    def codeword_of(self, equation_id: str) -> Tuple[str, int]:
        return self._assign[equation_id]

    def encode(self, mat: Sequence[Sequence[int]]) -> Tuple[str, Optional[int]]:
        """(codeword, index) for a matrix. Known matrices get their exact
        index (canonical entry for byte-identical duplicates); novel
        matrices get a cluster codeword with index None."""
        key = tuple(tuple(int(x) for x in row) for row in mat)
        if key in self._by_matrix:
            return self._by_matrix[key]
        return f"C{self.cluster_of(spectral_radius(key))}", None

    def decode(self, codeword: str, index: int) -> str:
        return self._decode[(codeword, index)]

    # ── reporting ────────────────────────────────────────────────────

    def cluster_table(self) -> List[dict]:
        clusters: Dict[str, List[str]] = {}
        for eid, (cw, _) in self._assign.items():
            clusters.setdefault(cw, []).append(eid)
        rows = []
        for cw in sorted(clusters, key=lambda c: int(c[1:])):
            ids = clusters[cw]
            lams = [self.profiles[e]["spectral_radius"] for e in ids]
            distinct = len({self.matrices[e] for e in ids})
            rows.append({
                "codeword": cw,
                "count": len(ids),
                "distinct_matrices": distinct,
                "lambda_min": round(min(lams), 6),
                "lambda_max": round(max(lams), 6),
                "sparse_tail_members": sum(1 for l in lams if self.is_sparse_tail(l)),
            })
        return rows

    def collision_report(self) -> dict:
        lam_groups: Dict[float, List[str]] = {}
        cp_groups: Dict[Tuple[int, ...], List[str]] = {}
        for eid, p in self.profiles.items():
            lam_groups.setdefault(round(p["spectral_radius"], LAMBDA_DP), []).append(eid)
            cp_groups.setdefault(tuple(p["charpoly"]), []).append(eid)

        def collide(groups):  # groups with >1 member
            return {k: sorted(v) for k, v in groups.items() if len(v) > 1}

        lam_c, cp_c = collide(lam_groups), collide(cp_groups)
        mat_c = {m: ids for m, ids in self.matrix_groups.items() if len(ids) > 1}
        # cospectral = same char poly but NOT byte-identical matrices
        dup_sets = [frozenset(v) for v in mat_c.values()]
        cospectral = {
            ",".join(map(str, k)): v
            for k, v in cp_c.items() if frozenset(v) not in dup_sets
        }
        pi_mismatch = [
            eid for eid, p in self.profiles.items()
            if abs(p["spectral_radius"] - p["lambda_power_iteration"]) > 1e-3
        ]
        n = len(self.profiles)
        return {
            "corpus_size": n,
            "distinct_matrices": len(self.matrix_groups),
            "duplicate_matrix_groups": len(mat_c),
            "duplicate_matrix_members": sum(len(v) for v in mat_c.values()),
            "duplicate_matrix_classes": [sorted(v) for v in mat_c.values()],
            "unique_lambda_4dp": len(lam_groups),
            "lambda_collision_groups": len(lam_c),
            "lambda_collision_members": sum(len(v) for v in lam_c.values()),
            "unique_charpoly": len(cp_groups),
            "charpoly_collision_groups": len(cp_c),
            "charpoly_collision_members": sum(len(v) for v in cp_c.values()),
            "charpoly_collision_classes": {
                ",".join(map(str, k)): v for k, v in cp_c.items()
            },
            "cospectral_nonidentical_groups": len(cospectral),
            "cospectral_nonidentical": cospectral,
            "identified_by_lambda_4dp": sum(1 for v in lam_groups.values() if len(v) == 1),
            "identified_by_charpoly": sum(1 for v in cp_groups.values() if len(v) == 1),
            "power_iteration_artifacts": sorted(pi_mismatch),
            "power_iteration_artifact_count": len(pi_mismatch),
        }

    def to_json(self) -> dict:
        entries = []
        for eid in sorted(self.profiles):
            cw, idx = self._assign[eid]
            p = self.profiles[eid]
            n = self.spiral_indices[eid]
            entry = {
                "equation_id": eid,
                "codeword": cw,
                "index": idx,
                "spectral_radius": p["spectral_radius"],
                "lambda_power_iteration": p["lambda_power_iteration"],
                "charpoly": p["charpoly"],
                "density": p["density"],
                "frobenius": p["frobenius"],
                "trace": p["trace"],
                "spiral_index": n,
                "layout": {
                    "radius_sq": n,
                    "angle_frac": round(angle_frac(n), 12),
                },
            }
            if self.is_sparse_tail(p["spectral_radius"]):
                entry["sparse_tail"] = True
            twins = [t for t in self.matrix_groups[self.matrices[eid]] if t != eid]
            if twins:
                entry["identical_matrix_ids"] = twins
            entries.append(entry)
        return {
            "schema": "spectral_codebook_v2",
            "generated": date.today().isoformat(),
            "source": "formal/SilverSight/PIST/Matrices250.lean",
            "matrix_count": len(entries),
            "lambda_semantics": (
                "spectral_radius is the exact max root modulus of the integer "
                "characteristic polynomial (charpoly); lambda_power_iteration is "
                "the legacy estimate and is unreliable on peripheral spectra"
            ),
            "phi_corkscrew": {
                "packing": {
                    "base": self.pack_base,
                    "offset": self.pack_offset,
                    "digits": 8,
                    "order": "little-endian: charpoly c1 is digit 0",
                    "semantics": (
                        "spiral_index = Σ (c_i + offset)·base^i — exact "
                        "integer, injective on charpoly fingerprints "
                        "(unpack_spiral_index inverts it exactly); inherits "
                        "charpoly's non-injectivity on matrices, so identity "
                        "still needs the within-cluster index"
                    ),
                },
                "layout": {
                    "psi": "2π/φ² (golden angle)",
                    "angle_frac": "frac(n·(3−√5)/2), decimal-exact precision",
                    "min_angular_separation_frac": (
                        min_angular_separation(list(self.spiral_indices.values()))
                    ),
                    "semantics": (
                        "decorative placement only: f(n) is injective because "
                        "1/φ² is irrational, and carries no information beyond "
                        "spiral_index"
                    ),
                },
            },
            "quantization": {
                "method": "gap_aware_deduped_min_support",
                "gap_factor": self.gap_factor,
                "min_support": self.min_support,
                "median_gap": round(self.quant["median_gap"], 6),
                "threshold": round(self.quant["threshold"], 6),
                "boundaries": self.boundaries,
                "candidate_boundaries": self.quant["candidates"],
                "suppressed_boundaries": self.quant["suppressed"],
                "sparse_tail_cutoff": self.quant["sparse_tail_cutoff"],
            },
            "clusters": self.cluster_table(),
            "collisions": self.collision_report(),
            "entries": entries,
        }


# ──────────────────────────────────────────────────────────────────────
# 278-row manifold cross-check
# ──────────────────────────────────────────────────────────────────────

def manifold_check(codebook: Codebook,
                   manifold_path: Path = MANIFOLD_LEAN) -> dict:
    """Run the pipeline over the full 278-row fixture corpus.

    The manifold repeats equation ids (278 rows over 250 unique ids), so the
    28 extra rows reuse matrices already in the codebook. After dedupe the
    λ multiset is unchanged, so boundaries must match the 250-matrix
    codebook exactly; this is reported rather than assumed.
    """
    ids = parse_manifold_ids(manifold_path)
    missing = sorted({e for e in ids if e not in codebook.profiles})
    seen: Dict[Matrix, str] = {}
    lambda_counts: Dict[float, int] = {}
    for eid in ids:
        if eid in codebook.profiles:
            mat = codebook.matrices[eid]
            if mat not in seen:  # dedupe, matching the codebook build
                seen[mat] = eid
                lam = codebook.profiles[eid]["spectral_radius"]
                lambda_counts[lam] = lambda_counts.get(lam, 0) + 1
    quant = gap_quantize(lambda_counts, codebook.gap_factor, codebook.min_support)
    return {
        "rows": len(ids),
        "unique_ids": len(set(ids)),
        "extra_rows": len(ids) - len(set(ids)),
        "ids_without_matrix": missing,
        "distinct_matrices": len(seen),
        "boundaries": quant["boundaries"],
        "median_gap": round(quant["median_gap"], 6),
        "boundaries_match_250": quant["boundaries"] == codebook.boundaries,
        "clusters": len(quant["boundaries"]) + 1,
    }


# ──────────────────────────────────────────────────────────────────────
# CLI
# ──────────────────────────────────────────────────────────────────────

def build_codebook(matrices_path: Path = MATRICES_LEAN,
                   gap_factor: float = GAP_FACTOR,
                   min_support: int = MIN_SUPPORT) -> Codebook:
    return Codebook(parse_matrices_lean(matrices_path), gap_factor, min_support)


def main(argv: Optional[List[str]] = None) -> int:
    ap = argparse.ArgumentParser(description="Gap-aware spectral codebook generator")
    ap.add_argument("--matrices", type=Path, default=MATRICES_LEAN,
                    help="Matrices250.lean path")
    ap.add_argument("--out", type=Path, default=DEFAULT_OUT,
                    help="Output codebook JSON (default data/spectral_codebook.json)")
    ap.add_argument("--gap-factor", type=float, default=GAP_FACTOR,
                    help="Boundary threshold = factor × median gap (default 3)")
    ap.add_argument("--min-support", type=int, default=MIN_SUPPORT,
                    help="Minimum distinct matrices per cluster (default 10)")
    ap.add_argument("--check-manifold", action="store_true",
                    help="Also run over the 278-row RRC/Q16_16Manifold corpus")
    ap.add_argument("--cartan-floor", action="store_true",
                    help="Also emit the Cartan ∆=17/1792 Fisher-floor report "
                         "(alongside, not replacing, the gap rule)")
    ap.add_argument("--no-write", action="store_true",
                    help="Report only; do not write the JSON")
    ap.add_argument("--sync-db", action="store_true",
                    help="DRY-RUN sync preview to ene.rrc_predictions on the "
                         "neon-64gb Postgres ($NEON_PG). Always a dry run from "
                         "this entry point; use python/spectral_codebook_db.py "
                         "--apply to actually insert.")
    args = ap.parse_args(argv)

    cb = build_codebook(args.matrices, args.gap_factor, args.min_support)
    doc = cb.to_json()

    col = doc["collisions"]
    q = doc["quantization"]
    print(f"Spectral codebook over {doc['matrix_count']} matrices "
          f"({col['distinct_matrices']} distinct)")
    print(f"  median gap {q['median_gap']:.6f} × {args.gap_factor:g} → "
          f"threshold {q['threshold']:.6f}; min support {args.min_support} "
          f"(candidates {len(q['candidate_boundaries'])}, "
          f"kept {len(q['boundaries'])}, "
          f"suppressed {len(q['suppressed_boundaries'])})")
    print(f"  boundaries: {q['boundaries']}")
    print(f"  sparse-tail cutoff: {q['sparse_tail_cutoff']}")
    print("\n  codeword   ids  distinct   lambda range          sparse-tail")
    for row in doc["clusters"]:
        print(f"  {row['codeword']:<8} {row['count']:>4}  {row['distinct_matrices']:>8}"
              f"   [{row['lambda_min']:.4f}, {row['lambda_max']:.4f}]"
              f"   {row['sparse_tail_members'] or ''}")
    print(f"\n  fingerprint uniqueness (of {col['corpus_size']} ids / "
          f"{col['distinct_matrices']} distinct matrices):")
    print(f"    λ @4dp      : {col['unique_lambda_4dp']} unique values, "
          f"{col['identified_by_lambda_4dp']} ids uniquely identified, "
          f"{col['lambda_collision_groups']} collision groups "
          f"({col['lambda_collision_members']} ids)")
    print(f"    char-poly ℤ : {col['unique_charpoly']} unique fingerprints, "
          f"{col['identified_by_charpoly']} ids uniquely identified, "
          f"{col['charpoly_collision_groups']} collision groups "
          f"({col['charpoly_collision_members']} ids)")
    print(f"    exact duplicate matrices: {col['duplicate_matrix_groups']} groups "
          f"({col['duplicate_matrix_members']} ids)")
    print(f"    cospectral non-identical groups: "
          f"{col['cospectral_nonidentical_groups']}")
    print(f"    power-iteration artifacts (|λ_pi − ρ| > 1e-3): "
          f"{col['power_iteration_artifact_count']}")

    if args.check_manifold:
        mc = manifold_check(cb)
        doc["manifold_278"] = mc
        print(f"\n  278-row manifold: {mc['rows']} rows / {mc['unique_ids']} unique ids "
              f"({mc['extra_rows']} duplicate rows, "
              f"{mc['distinct_matrices']} distinct matrices)")
        print(f"    boundaries match 250-matrix codebook: {mc['boundaries_match_250']} "
              f"({mc['clusters']} clusters)")

    if args.cartan_floor:
        cf = cartan_floor_report(cb.profiles)
        doc["cartan_floor"] = cf
        gap_clusters = len(q["boundaries"]) + 1
        print(f"\n  Cartan ∆-floor (∆ = {cf['delta_exact']} ≈ "
              f"{cf['delta_float']:.6f}, {cf['delta_source']}):")
        print(f"    distinct fingerprints on Δ₇ : {cf['distinct_fingerprints']}")
        print(f"    pairs below ∆ (Fisher)      : {cf['pairs_below_delta']}")
        print(f"    min Fisher distance         : {cf['min_fisher_distance']} "
              f"{cf['min_fisher_pair']}")
        print(f"    rule comparison — codewords : "
              f"∆-floor {cf['delta_resolution_codewords']} vs "
              f"3×-median-gap {gap_clusters} (λ clusters); both emitted, "
              f"neither replaces the other")

    if args.sync_db:
        import spectral_codebook_db
        print("\n  DB sync (dry run):")
        spectral_codebook_db.sync_db(codebook=cb)  # never applies from here

    if not args.no_write:
        args.out.parent.mkdir(parents=True, exist_ok=True)
        args.out.write_text(json.dumps(doc, indent=2) + "\n")
        print(f"\nWrote {args.out}")
    return 0


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