diff --git a/python/pist_braid_bridge.py b/python/pist_braid_bridge.py new file mode 100644 index 00000000..508913f4 --- /dev/null +++ b/python/pist_braid_bridge.py @@ -0,0 +1,630 @@ +#!/usr/bin/env python3 +""" +PIST + BRAID BRIDGE — Connecting Sprint Results to Lean Formalism +=================================================================== + +This module bridges the integration_sprint.py results with the Lean 4 +formal infrastructure: + - PIST.Spectral (Q16_16 fixed-point spectral analysis) + - BraidEigensolid (Torus carrier, golden centering, eigensolid) + - BraidField (PIST operator B/G/A/P, RG flow) + - BraidSpherionBridge (Cross-mode ↔ eigensolid correspondence) + +Mappings implemented: + 1. Spectral coefficients → Q16_16 → PIST SpectralProfile + 2. Spiral index → TorusWinding (a, b) on T² + 3. 4-mode receipt → PIST field (B, G, A, P) + 4. Cross-mode agreement → Eigensolid convergence check + 5. Φ-corkscrew ↔ Golden centering (φ⁻¹ = 40560/65536) + +Author: SilverSight Integration Layer +""" + +import json +import math +import struct +from typing import Dict, List, Tuple, Optional +from dataclasses import dataclass, asdict + +import numpy as np + +# ======================================================================== +# §1. Q16_16 FIXED-POINT ARITHMETIC (matches Lean SilverSight.FixedPoint) +# ======================================================================== + +Q16_SCALE = 65536 # 2^16 +Q16_MAX_RAW = 2147483647 # 2^31 - 1 +Q16_MIN_RAW = -2147483648 # -2^31 + + +def q16_clamp(x: int) -> int: + """Saturating clamp to Q16_16 range.""" + if x > Q16_MAX_RAW: + return Q16_MAX_RAW + if x < Q16_MIN_RAW: + return Q16_MIN_RAW + return x + + +def float_to_q16(x: float) -> int: + """Convert float to Q16_16 raw value.""" + return q16_clamp(int(round(x * Q16_SCALE))) + + +def q16_to_float(raw: int) -> float: + """Convert Q16_16 raw value to float.""" + return raw / Q16_SCALE + + +def q16_add(a: int, b: int) -> int: + return q16_clamp(a + b) + + +def q16_sub(a: int, b: int) -> int: + return q16_clamp(a - b) + + +def q16_mul(a: int, b: int) -> int: + """Multiply two Q16_16 values, keeping Q16_16 precision.""" + return q16_clamp((a * b) // Q16_SCALE) + + +def q16_div(a: int, b: int) -> int: + """Divide two Q16_16 values.""" + if b == 0: + return Q16_MAX_RAW + return q16_clamp((a * Q16_SCALE) // b) + + +def q16_sqrt(raw: int) -> int: + """Integer square root of Q16_16 value. Returns floor(√x) in Q16_16.""" + if raw <= 0: + return 0 + # √x in Q16_16 = √(x / 65536) * 65536 = √(x * 65536) + val = raw if raw >= 0 else 0 + return int(math.isqrt(val * Q16_SCALE)) + + +# Golden centering constant: φ⁻¹ ≈ 0.618033988... +# In Lean: Q16_16.ofRawInt 40560 +# 40560 / 65536 = 0.618896484375 +PHI_INV_Q16 = 40560 +PHI = (1.0 + 5.0 ** 0.5) / 2.0 +PSI = 2.0 * math.pi / (PHI ** 2) + + +# ======================================================================== +# §2. PIST.SPECTRAL BRIDGE +# ======================================================================== + +@dataclass +class PISTSpectralProfile: + """Python mirror of SilverSight.PIST.Spectral.SpectralProfile.""" + matrix_size: int # Q16_16 raw + rank: int # Q16_16 raw + spectral_gap: int # Q16_16 raw + density: int # Q16_16 raw + trace_val: int # Q16_16 raw + frobenius_norm: int # Q16_16 raw + laplacian_zero_count: int # Q16_16 raw + adjacency_eigenvalue_max: int # Q16_16 raw + laplacian_eigenvalue_max: int # Q16_16 raw + singular_value_max: int # Q16_16 raw + + def to_float_dict(self) -> Dict[str, float]: + return { + "matrix_size": q16_to_float(self.matrix_size), + "rank": q16_to_float(self.rank), + "spectral_gap": q16_to_float(self.spectral_gap), + "density": q16_to_float(self.density), + "trace_val": q16_to_float(self.trace_val), + "frobenius_norm": q16_to_float(self.frobenius_norm), + "laplacian_zero_count": q16_to_float(self.laplacian_zero_count), + "adjacency_eigenvalue_max": q16_to_float(self.adjacency_eigenvalue_max), + "laplacian_eigenvalue_max": q16_to_float(self.laplacian_eigenvalue_max), + "singular_value_max": q16_to_float(self.singular_value_max), + } + + +def power_iteration_q16(mat: List[List[int]], max_iter: int = 100) -> int: + """Dominant eigenvalue via power iteration in Q16_16. + + Matches: SilverSight.PIST.Spectral.powerIteration + """ + n = len(mat) + if n == 0: + return 0 + + init_val = Q16_SCALE // n + v = [init_val] * n + + for _ in range(max_iter): + # Matrix-vector multiply + mv = [0] * n + for i in range(n): + s = 0 + for j in range(n): + s += mat[i][j] * (v[j] if j < len(v) else 0) + mv[i] = s # Keep in raw integer space + + # Norm squared + nm2 = sum(x * x for x in mv) + if nm2 <= 0: + break + + nm = int(math.isqrt(nm2)) + if nm == 0: + break + + # Normalize + v = [q16_clamp((x * Q16_SCALE) // nm) for x in mv] + + # Rayleigh quotient: v^T A v / v^T v + mv = [0] * n + for i in range(n): + s = 0 + for j in range(n): + s += mat[i][j] * (v[j] if j < len(v) else 0) + mv[i] = s + + num = sum(v[i] * mv[i] for i in range(n)) + den = sum(x * x for x in v) + if den <= 0: + return 0 + return q16_clamp((num * Q16_SCALE) // den) + + +def compute_pist_spectral(eigenvalues: np.ndarray) -> PISTSpectralProfile: + """Compute PIST spectral profile from eigenvalues (float → Q16_16). + + Mirrors: SilverSight.PIST.Spectral.computeSpectral + """ + if len(eigenvalues) == 0: + return PISTSpectralProfile(*([0] * 10)) + + k = min(8, len(eigenvalues)) + evals = np.sort(eigenvalues)[::-1][:k] + + # Build 8x8 symmetric matrix from eigenvalues (diagonalization approximation) + n = 8 + mat = [[0] * n for _ in range(n)] + for i in range(min(len(evals), n)): + val = float_to_q16(float(evals[i])) + mat[i][i] = val + # Add some off-diagonal coupling for realism + if i + 1 < n: + coupling = float_to_q16(float(evals[i]) * 0.1) + mat[i][i + 1] = coupling + mat[i + 1][i] = coupling + + # Symmetrize + sym = [[(mat[i][j] + mat[j][i]) // 2 for j in range(n)] for i in range(n)] + + # Laplacian + lap = [[0] * n for _ in range(n)] + for i in range(n): + deg = sum(sym[i]) + for j in range(n): + lap[i][j] = deg if i == j else -sym[i][j] + + ev_max = power_iteration_q16(sym) + + # Shift-deflation for second eigenvalue + shift = (ev_max * 58982) // Q16_SCALE # 0.9 * ev_max + shifted = [[sym[i][j] for j in range(n)] for i in range(n)] + for i in range(n): + shifted[i][i] = q16_sub(shifted[i][i], shift) + ev_shift = power_iteration_q16(shifted) + ev_second = q16_sub(ev_max, ev_shift) if ev_shift < ev_max else ev_max + spectral_gap = q16_sub(ev_max, ev_second) + + lap_max = power_iteration_q16(lap) + + # Density + total = sum(sum(row) for row in mat) + density = q16_clamp((total * Q16_SCALE) // (n * n)) if n > 0 else 0 + + # Trace + trace_val = sum(mat[i][i] for i in range(n)) + + # Frobenius norm + frob_sq = sum(sum(x * x for x in row) for row in mat) + frob_norm = int(math.isqrt(frob_sq)) + + # Rank (count non-zero rows) + rank = sum(1 for row in mat if any(x != 0 for x in row)) + + # Laplacian zero count + lap_zero = sum(1 for i in range(n) if sum(mat[i]) - mat[i][i] == 0) + + # ATA max singular value + ata = [[0] * n for _ in range(n)] + for i in range(n): + for j in range(n): + s = 0 + for k_idx in range(n): + s += mat[k_idx][i] * mat[k_idx][j] + ata[i][j] = s + ata_ev = power_iteration_q16(ata) + sv_max = q16_sqrt(ata_ev) + + return PISTSpectralProfile( + matrix_size=float_to_q16(n), + rank=float_to_q16(rank), + spectral_gap=spectral_gap, + density=density, + trace_val=float_to_q16(trace_val), + frobenius_norm=float_to_q16(frob_norm), + laplacian_zero_count=float_to_q16(lap_zero), + adjacency_eigenvalue_max=ev_max, + laplacian_eigenvalue_max=lap_max, + singular_value_max=sv_max, + ) + + +# ======================================================================== +# §3. TORUS WINDING ↔ Φ-CORKSCREW +# ======================================================================== + +@dataclass +class TorusWinding: + """Torus winding counts: H₁(T²; Z) = Z⟨a⟩ ⊕ Z⟨b⟩. + + Matches: SilverSight.BraidEigensolid.TorusWinding + a = spatial winding (latitude cycle) + b = phase/torsion winding (longitude cycle) + """ + a: int # Q16_16 raw + b: int # Q16_16 raw + + def to_float(self) -> Tuple[float, float]: + return q16_to_float(self.a), q16_to_float(self.b) + + +def spiral_to_torus_winding(spiral_index: int) -> TorusWinding: + """Map spiral index to torus winding counts. + + The golden spiral f(n) = (√n·cos(nψ), √n·sin(nψ)) with ψ = 2π/φ² + lives on a torus because ψ/2π = 1/φ² is irrational. + + Winding counts: + a = spiral_index mod φ_step (spatial winding) + b = spiral_index // φ_step (phase winding) + + where φ_step = round(φ) = 2 (since φ ≈ 1.618) + + This mirrors: TorusWinding in BraidEigensolid.lean + """ + phi_step = int(round(PHI)) # 2 + a_raw = float_to_q16(float(spiral_index % phi_step)) + b_raw = float_to_q16(float(spiral_index // phi_step)) + return TorusWinding(a=a_raw, b=b_raw) + + +def torus_to_spiral_index(w: TorusWinding) -> int: + """Inverse: recover spiral index from torus winding. + + n = φ_step · b + a (in integer approximation) + """ + phi_step = int(round(PHI)) + b_int = int(round(q16_to_float(w.b))) + a_int = int(round(q16_to_float(w.a))) + return phi_step * b_int + a_int + + +# ======================================================================== +# §4. PIST FIELD ↔ 4-MODE SPRINT RECEIPT +# ======================================================================== + +@dataclass +class PISTField: + """PIST operator: q_{t+1} = PIST(q_t; B, G, A, P). + + Matches: SilverSight.BraidField.PISTField + B = Burden (load, cost, attention) + G = Geometry (basins, manifolds, gradients) + A = Adaptation (sorting rate, pacing, convergence) + P = Protection (compression, thresholding, overload) + """ + burden: int # Q16_16 raw + geometry: int # Q16_16 raw + adaptation: int # Q16_16 raw + protection: int # Q16_16 raw + + def to_float_dict(self) -> Dict[str, float]: + return { + "burden": q16_to_float(self.burden), + "geometry": q16_to_float(self.geometry), + "adaptation": q16_to_float(self.adaptation), + "protection": q16_to_float(self.protection), + } + + +def sprint_receipt_to_pist(result: Dict) -> PISTField: + """Convert a single sprint mode result to PIST field. + + B = execution_time_ms (how expensive was this mode?) + G = dominant_eigenvalue + spectral_gap (what did we learn about the graph?) + A = checks_passed / checks_total (how well did we adapt?) + P = 1.0 if checks_passed == checks_total else 0.0 (did we protect against errors?) + """ + B = float_to_q16(result.get("executionTimeMs", 0.0) / 100.0) + + geom = result.get("dominantEigenvalue", 0.0) + result.get("spectralGap", 0.0) + G = float_to_q16(geom / 10.0) + + checks_total = result.get("checks_total", 7) + checks_passed = result.get("checks_passed", 7) + A = float_to_q16(checks_passed / checks_total if checks_total > 0 else 0.0) + + P = Q16_SCALE if checks_passed == checks_total else Q16_SCALE // 2 + + return PISTField(burden=B, geometry=G, adaptation=A, protection=P) + + +def compare_pist_fields(p1: PISTField, p2: PISTField, tolerance: float = 0.05) -> bool: + """Check if two PIST fields agree within Q16_16 tolerance.""" + for field in ["burden", "geometry", "adaptation", "protection"]: + v1 = q16_to_float(getattr(p1, field)) + v2 = q16_to_float(getattr(p2, field)) + if abs(v1 - v2) > tolerance: + return False + return True + + +# ======================================================================== +# §5. EIGENSOLID CONVERGENCE CHECK +# ======================================================================== + +def check_eigensolid_convergence(results: List[Dict]) -> Dict: + """Verify that all modes converge to the same eigensolid. + + Theorem analog: SilverSight.BraidEigensolid.eigensolid_convergence + + All modes should produce: + - Same dominant eigenvalue (within numerical tolerance) + - Same spiral index (exact agreement) + - Same DNA sequence (exact agreement) + - Same PIST field (within Q16_16 tolerance) + """ + if len(results) < 2: + return {"converged": False, "reason": "need_at_least_2_modes"} + + # Check exact agreement on discrete invariants + spirals = [r.get("spiral_index") for r in results] + dnas = [r.get("dnaSequence") for r in results] + + spiral_agree = len(set(spirals)) == 1 + dna_agree = len(set(dnas)) == 1 + + # Check approximate agreement on continuous invariants + lambdas = [r.get("dominantEigenvalue", 0.0) for r in results] + lambda_cv = np.std(lambdas) / np.mean(lambdas) if np.mean(lambdas) > 0 else 0 + + gaps = [r.get("spectralGap", 0.0) for r in results] + gap_cv = np.std(gaps) / np.mean(gaps) if np.mean(gaps) > 0 else 0 + + # PIST field agreement + pist_fields = [sprint_receipt_to_pist(r) for r in results] + pist_agree = all( + compare_pist_fields(pist_fields[0], pf) for pf in pist_fields[1:] + ) + + converged = spiral_agree and dna_agree and lambda_cv < 0.01 and gap_cv < 0.01 and pist_agree + + return { + "converged": converged, + "spiral_agree": spiral_agree, + "dna_agree": dna_agree, + "lambda_cv": float(lambda_cv), + "gap_cv": float(gap_cv), + "pist_agree": pist_agree, + "n_modes": len(results), + "pist_fields": [pf.to_float_dict() for pf in pist_fields], + } + + +# ======================================================================== +# §6. GOLDEN CENTERING ↔ Φ-CORKSCREW +# ======================================================================== + +def golden_centering_compress(q16_val: int) -> int: + """Apply golden centering contraction: x → φ⁻¹ · x. + + Matches: goldenCentering in BraidEigensolid.lean + Uses the fixed-point constant 40560 = φ⁻¹ in Q16_16. + + This is the compressor's contraction envelope — it keeps all crossing + weights in the Q0_2 range [0, 16384] after finite iteration. + """ + return q16_mul(q16_val, PHI_INV_Q16) + + +def phi_corkscrew_to_q16(spiral_index: int, n_coeffs: int = 9) -> List[int]: + """Convert spiral index to Q16_16 spectral coefficients. + + The Φ-corkscrew packs spectral info into a single integer. + Unpack: recover coefficients via base-φ decomposition. + """ + coeffs = [] + n = spiral_index + for i in range(n_coeffs): + digit = n % 8 # Hachimoji base-8 + coeffs.append(float_to_q16(digit / 8.0)) + n //= 8 + return coeffs + + +def q16_to_phi_corkscrew(coeffs: List[int]) -> int: + """Pack Q16_16 spectral coefficients into spiral index (inverse).""" + spiral = 0 + for i, c in enumerate(coeffs): + digit = int(round(q16_to_float(c) * 8)) % 8 + spiral += digit * (8 ** i) + return spiral + + +# ======================================================================== +# §7. OCTAGONAL NORM ↔ FISHER-RAO BOUND +# ======================================================================== + +def octagonal_norm(x: float, y: float) -> float: + """Octagonal norm: κ ≈ max(|x|,|y|) + 3/8·min(|x|,|y|). + + Matches: PhaseVec.normApprox in BraidBracket.lean + """ + ax, ay = abs(x), abs(y) + hi, lo = max(ax, ay), min(ax, ay) + return hi + (3.0 / 8.0) * lo + + +def fisher_rao_distance(p: np.ndarray, q: np.ndarray) -> float: + """Fisher-Rao distance on probability simplex: arccos(Σ√(pᵢqᵢ)).""" + bc = np.sum(np.sqrt(p * q)) + bc = np.clip(bc, -1.0, 1.0) + return np.arccos(bc) + + +def octagonal_fisher_bound(braid_state_8d: np.ndarray) -> float: + """Show that octagonal norm upper-bounds Fisher-Rao distance. + + For an 8-strand braid state interpreted as a point on Δ₇, + the octagonal norm of the phase vector upper-bounds the + Fisher-Rao distance from the uniform distribution. + + Returns: ratio octagonal_norm / fisher_rao_distance (should be ≥ 1) + """ + # Normalize to probability simplex + s = braid_state_8d.sum() + if s <= 0: + return float('inf') + p = braid_state_8d / s + uniform = np.ones(8) / 8.0 + + fisher_dist = fisher_rao_distance(p, uniform) + + # Octagonal norm on first two coordinates (phase vector) + x, y = p[0] - 1.0 / 8.0, p[1] - 1.0 / 8.0 + oct_norm = octagonal_norm(x, y) + + ratio = oct_norm / fisher_dist if fisher_dist > 1e-12 else float('inf') + return float(ratio) + + +# ======================================================================== +# §8. MAIN — Load sprint receipt and run full bridge +# ======================================================================== + +def run_bridge(sprint_receipt_path: str = "sprint_receipt.json") -> Dict: + """Load sprint receipt and run all bridge analyses.""" + try: + with open(sprint_receipt_path, 'r') as f: + receipt = json.load(f) + except FileNotFoundError: + return {"error": f"Receipt not found: {sprint_receipt_path}"} + + results = receipt.get("results", []) + if not results: + return {"error": "No results in receipt"} + + output = { + "quimb_backend": receipt.get("quimb_backend", "unknown"), + "n_modes": len(results), + "bridge_analyses": {}, + } + + # 1. PIST Spectral for each mode + spectral_profiles = [] + for r in results: + # Create synthetic 8x8 from the eigenvalues + ev_max = r.get("eigenvalue_max", 1.0) + ev_min = r.get("eigenvalue_min", 0.0) + gap = r.get("spectralGap", 0.0) + # Create a spectrum: dominant, dominant-gap, and the rest + spectrum = np.array([ev_max, ev_max - gap] + [ev_min] * 6) + profile = compute_pist_spectral(spectrum) + spectral_profiles.append(profile) + + output["bridge_analyses"]["pist_spectral"] = { + "profiles": [p.to_float_dict() for p in spectral_profiles], + "mean_gap_q16": q16_to_float( + sum(p.spectral_gap for p in spectral_profiles) // len(spectral_profiles) + ), + } + + # 2. Torus winding for each mode + torus_windings = [] + for r in results: + tw = spiral_to_torus_winding(r.get("spiral_index", 0)) + torus_windings.append(tw) + + output["bridge_analyses"]["torus_winding"] = { + "windings": [{"a": tw.to_float()[0], "b": tw.to_float()[1]} + for tw in torus_windings], + "inverse_check": [ + torus_to_spiral_index(tw) == r.get("spiral_index", 0) + for tw, r in zip(torus_windings, results) + ], + } + + # 3. PIST field for each mode + pist_fields = [sprint_receipt_to_pist(r) for r in results] + output["bridge_analyses"]["pist_field"] = { + "fields": [pf.to_float_dict() for pf in pist_fields], + "all_agree": all(compare_pist_fields(pist_fields[0], pf) + for pf in pist_fields[1:]), + } + + # 4. Eigensolid convergence check + output["bridge_analyses"]["eigensolid_convergence"] = \ + check_eigensolid_convergence(results) + + # 5. Golden centering compression check + max_eigenvalue_q16 = float_to_q16( + max(r.get("dominantEigenvalue", 0.0) for r in results) + ) + compressed = golden_centering_compress(max_eigenvalue_q16) + output["bridge_analyses"]["golden_centering"] = { + "original_q16": max_eigenvalue_q16, + "compressed_q16": compressed, + "phi_inv": PHI_INV_Q16, + "phi_inv_float": q16_to_float(PHI_INV_Q16), + } + + # 6. Φ-corkscrew ↔ Q16_16 roundtrip + if results: + spiral = results[0].get("spiral_index", 0) + coeffs = phi_corkscrew_to_q16(spiral) + recovered = q16_to_phi_corkscrew(coeffs) + output["bridge_analyses"]["phi_corkscrew_roundtrip"] = { + "original": spiral, + "recovered": recovered, + "match": spiral == recovered, + "n_coeffs": len(coeffs), + } + + # Overall verdict + conv = output["bridge_analyses"]["eigensolid_convergence"] + pist_agree = output["bridge_analyses"]["pist_field"]["all_agree"] + corkscrew_ok = output["bridge_analyses"]["phi_corkscrew_roundtrip"]["match"] + output["bridge_verdict"] = { + "all_checks_pass": conv.get("converged", False) and pist_agree and corkscrew_ok, + "eigensolid_converged": conv.get("converged", False), + "pist_fields_agree": pist_agree, + "corkscrew_roundtrip_ok": corkscrew_ok, + "formal_backing": [ + "BraidEigensolid.eigensolid_convergence", + "BraidEigensolid.receipt_invertible", + "BraidSpherionBridge.receipt_correspondence", + "BraidSpherionBridge.braidCross_merge_correspondence", + ], + } + + return output + + +if __name__ == "__main__": + import os + receipt_path = os.path.join(os.path.dirname(__file__), "sprint_receipt.json") + result = run_bridge(receipt_path) + print(json.dumps(result, indent=2, default=str))