mirror of
https://github.com/allaunthefox/SilverSight.git
synced 2026-07-31 01:25:21 +00:00
allaun
4ea68bcecc
Implements the Next Steps of docs/SPECTRAL_CODEBOOK_ANALYSIS.md as python/spectral_codebook.py (stdlib-only; NumPy optional fast path): - Parses the 250 8x8 braid adjacency matrices from PIST/Matrices250.lean. - Primary fingerprint: exact integer characteristic-polynomial coefficients via Faddeev-LeVerrier in Fraction arithmetic (196 unique over 238 distinct matrices, vs 179 unique lambda at 4dp; 6 cospectral non-identical groups; 11 exact-duplicate matrix groups / 23 ids). - Corrects the analysis doc: the 9 'mid band' lambda in (0.5,1) are power-iteration non-convergence artifacts - exact rho = 1.0 for all 9 (peripheral spectra). spectral_radius is now the exact max root modulus (numpy eigvals or Durand-Kerner on the exact char poly); power-iteration lambda kept only for traceability. - Gap-aware quantization: dedupe to 238 distinct matrices, boundaries at gaps > 3x median gap, min-support guard (>=10 distinct per cluster), sparse-tail outlier flagging above lambda ~= 7.66. Result: 9 clusters. - Round trip encode(matrix) -> (codeword, index) -> decode -> equation_id verified bijective over all 250 in tests/test_spectral_codebook.py. - 278-row RRC/Q16_16Manifold corpus: 28 extra rows are repeated ids; boundaries reproduce exactly, no new gaps or clusters. - Emits data/spectral_codebook.json (schema spectral_codebook_v2) with explicit collision classes; docs/SPECTRAL_CODEBOOK_GENERATOR.md notes the hashMatrix base-5 injectivity gap and the ClassifyN threshold (1.5/4.0 Q16.16) vs analysis-doc (0.5/1.0) mismatch. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
644 lines
28 KiB
Python
644 lines
28 KiB
Python
#!/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))),
|
||
}
|
||
|
||
|
||
# ──────────────────────────────────────────────────────────────────────
|
||
# 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()
|
||
}
|
||
|
||
# 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]
|
||
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"],
|
||
}
|
||
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"
|
||
),
|
||
"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("--no-write", action="store_true",
|
||
help="Report only; do not write the JSON")
|
||
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 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())
|