SilverSight/python/nr_bracket_validation.py

#!/usr/bin/env python3
"""
nr_bracket_validation.py — Real-data validation of d_CE μ = 0 (Gate C)

Tests the Nijenhuis–Richardson bracket breakglass on actual braid state data
from the RRC EntropyCandidates pipeline.  Two test suites:

  Suite A — Basis vectors v_i = e_i − e_7 (Fin 7 basis of V).
            Expected: ‖d_CE μ‖_∞ = 0 (exact) — mirrors Lean Jacobiator_basis_all.

  Suite B — Real strand phase vectors from RRC entropy candidates.
            Extracts phaseAcc.x (and .y) from each of the 8 strands,
            projects into V = ker(Σ), and evaluates d_CE μ.
            Should be near machine-zero if the formal proof matches reality.

Reference: SilverSight/formal/SilverSight/PIST/CartanConnection.lean
"""

from __future__ import annotations
import json
import math
import re
import sys

# ─── Crossing matrix C (exact rationals) ─────────────────────────────────────

SIGMA = 39 / 256  # diagonal weight
TAU = 1 / 7       # same-block off-diagonal weight


def crossing_matrix() -> list[list[float]]:
    C = [[0.0] * 8 for _ in range(8)]
    for i in range(8):
        for j in range(8):
            if i == j:
                C[i][j] = SIGMA
            elif i // 2 == j // 2:
                C[i][j] = TAU
    return C


C = crossing_matrix()


def mat_vec_mul(M: list[list[float]], v: list[float]) -> list[float]:
    return [sum(M[i][j] * v[j] for j in range(8)) for i in range(8)]


def vec_add(v: list[float], w: list[float]) -> list[float]:
    return [a + b for a, b in zip(v, w)]


def vec_sub(v: list[float], w: list[float]) -> list[float]:
    return [a - b for a, b in zip(v, w)]


def inf_norm(v: list[float]) -> float:
    return max(abs(x) for x in v)


def project_into_V(v: list[float]) -> list[float]:
    mean = sum(v) / 8.0
    return [x - mean for x in v]


def mu(X: list[float], Y: list[float]) -> list[float]:
    CX = mat_vec_mul(C, X)
    CY = mat_vec_mul(C, Y)
    return [CX[k] * Y[k] - X[k] * CY[k] for k in range(8)]


def jacobiator(X: list[float], Y: list[float], Z: list[float]) -> list[float]:
    return vec_add(
        mu(mu(X, Y), Z),
        vec_add(mu(mu(Y, Z), X), mu(mu(Z, X), Y))
    )


# ─── Suite A: Basis vectors ──────────────────────────────────────────────────

def basis_vec(k: int) -> list[float]:
    v = [0.0] * 8
    v[k] = 1.0
    v[7] = -1.0
    return v


def run_suite_a() -> dict:
    non_zero: list[tuple[int, int, int, float]] = []
    worst = 0.0
    for i in range(7):
        vi = basis_vec(i)
        for j in range(7):
            vj = basis_vec(j)
            for k in range(7):
                vk = basis_vec(k)
                J = jacobiator(vi, vj, vk)
                nrm = inf_norm(J)
                if nrm > worst:
                    worst = nrm
                if nrm > 1e-12:
                    non_zero.append((i, j, k, nrm))
    return {
        "total_triples": 7 ** 3,
        "non_zero_count": len(non_zero),
        "worst_inf_norm": worst,
        "all_zero": worst < 1e-12,
    }


# ─── Suite B: Parse Candidates.lean (stateful line-based) ────────────────────

def parse_candidates(lines: list[str]) -> list[list[tuple[int, int, int, int]]]:
    """
    Extract (x_raw, y_raw, slot, kappa_raw) per strand per candidate.

    Each candidate has 8 strand blocks like:
      | ⟨0, _⟩ => { phaseAcc := { x := Q16_16.ofRawInt 37813, y := Q16_16.ofRawInt 45787 }
        , parity := false
        , slot := 1
        , residue := Q16_16.ofRawInt 0
        , jitter := Q16_16.ofRawInt 0
        , bracket := { lower := ...
                        , kappa := Q16_16.ofRawInt 44921
    """
    x_pat = re.compile(r"x\s*:=\s*Q16_16\.ofRawInt\s+(-?\d+)")
    y_pat = re.compile(r"y\s*:=\s*Q16_16\.ofRawInt\s+(-?\d+)")
    slot_pat = re.compile(r"slot\s*:=\s*(\d+)")
    kappa_pat = re.compile(r"kappa\s*:=\s*Q16_16\.ofRawInt\s+(-?\d+)")

    candidates: list[list[tuple[int, int, int, int]]] = []
    current: list[tuple[int, int, int, int]] = []
    buf: dict[str, int | None] = {"x": None, "y": None, "slot": None, "kappa": None}

    def flush():
        nonlocal buf, current
        if all(v is not None for v in buf.values()):
            current.append((
                buf["x"], buf["y"], buf["slot"], buf["kappa"]
            ))
        buf = {"x": None, "y": None, "slot": None, "kappa": None}

    for line in lines:
        if "phaseAcc" in line:
            flush()  # commit previous strand if pending

        x_m = x_pat.search(line)
        y_m = y_pat.search(line)
        slot_m = slot_pat.search(line)
        kappa_m = kappa_pat.search(line)

        if x_m:
            buf["x"] = int(x_m.group(1))
        if y_m:
            buf["y"] = int(y_m.group(1))
        if slot_m:
            buf["slot"] = int(slot_m.group(1))
        if kappa_m:
            buf["kappa"] = int(kappa_m.group(1))

        # When we have all 4 fields, make a strand entry
        if all(v is not None for v in buf.values()):
            flush()
            if len(current) == 8:
                candidates.append(current)
                current = []

    return candidates


def run_suite_b(candidates: list[list[tuple[int, int, int, int]]]) -> dict:
    Q16 = 65536.0
    results = []

    for idx, cand in enumerate(candidates):
        X = [s[0] / Q16 for s in cand]
        Y = [s[1] / Q16 for s in cand]
        Z = [s[3] / Q16 for s in cand]

        Xv = project_into_V(X)
        Yv = project_into_V(Y)
        Zv = project_into_V(Z)

        tests: list[tuple[str, list[float], list[float], list[float]]] = [
            ("X,X,X", Xv, Xv, Xv),
            ("Y,Y,Y", Yv, Yv, Yv),
            ("Z,Z,Z", Zv, Zv, Zv),
            ("X,Y,Z", Xv, Yv, Zv),
            ("X,X,Y", Xv, Xv, Yv),
            ("Y,Y,X", Yv, Yv, Xv),
            ("X,Y,Y", Xv, Yv, Yv),
            ("X,Z,Y", Xv, Zv, Yv),
        ]

        worst = 0.0
        details = []
        for label, A, B, D in tests:
            J = jacobiator(A, B, D)
            nrm = inf_norm(J)
            if nrm > worst:
                worst = nrm
            details.append({"triple": label, "inf_norm": round(nrm, 12)})

        results.append({
            "candidate_idx": idx,
            "worst_inf_norm": round(worst, 12),
            "all_near_zero": worst < 1e-6,
            "details": details,
        })

    overall_worst = max((r["worst_inf_norm"] for r in results), default=0.0)
    return {
        "candidates_tested": len(candidates),
        "results": results,
        "overall_worst": overall_worst,
        "overall_pass": all(r["all_near_zero"] for r in results),
    }


# ─── Main ─────────────────────────────────────────────────────────────────────

def main():
    print("=" * 60)
    print("NR Bracket Validation — d_CE μ = 0 on real data")
    print("=" * 60)

    # Suite A
    print("\n─── Suite A: Basis vectors (343 triples) ───")
    sa = run_suite_a()
    print(f"  Total triples:     {sa['total_triples']}")
    print(f"  Non-zero found:    {sa['non_zero_count']}")
    print(f"  Worst ‖J‖_∞:       {sa['worst_inf_norm']:.2e}")
    print(f"  All zero (exact):  {sa['all_zero']}")

    candidates_file = (
        "/home/allaun/research-stack/lean/Semantics/"
        "Semantics/RRC/EntropyCandidates/Candidates.lean"
    )
    print(f"\n─── Suite B: Real candidates ({candidates_file}) ───")

    sb = None
    try:
        with open(candidates_file) as f:
            lines = f.readlines()
        candidates = parse_candidates(lines)
        print(f"  Candidates parsed: {len(candidates)}")
        if len(candidates) == 0:
            print("  ✗ WARNING: No candidates parsed — check regex patterns")
            # Debug: show first 30 lines containing relevant keywords
            for i, ln in enumerate(lines[:200]):
                if any(kw in ln for kw in ["phaseAcc", "slot", "kappa", "ofRawInt"]):
                    print(f"    L{i+1}: {ln.rstrip()}")
        else:
            sb = run_suite_b(candidates)
            print(f"  Tested:            {sb['candidates_tested']} candidates")
            print(f"  Overall worst ‖J‖_∞: {sb['overall_worst']:.6e}")
            print(f"  Overall pass:      {sb['overall_pass']}")
            for cr in sb["results"]:
                status = "✓" if cr["all_near_zero"] else "✗"
                print(f"    {status} Candidate {cr['candidate_idx']}: "
                      f"worst ‖J‖_∞ = {cr['worst_inf_norm']:.6e}")
                for d in cr["details"]:
                    if d["inf_norm"] > 1e-9:
                        print(f"       ⚠ {d['triple']}: {d['inf_norm']:.6e}")
    except FileNotFoundError:
        print(f"  ✗ Candidates.lean not found at {candidates_file}")

    print("\n" + "=" * 60)
    verdict_a = "✓ PASS" if sa["all_zero"] else "✗ FAIL"
    print(f"  Suite A (basis):  {verdict_a}")
    if len(candidates) > 0:
        verdict_b = "✓ PASS" if sb["overall_pass"] else "✗ FAIL"
        print(f"  Suite B (real):   {verdict_b}")
    print("=" * 60)

    # Save JSON receipt
    receipt = {
        "schema": "nr_bracket_validation_v1",
        "suite_a": sa,
        "suite_b": sb if len(candidates) > 0 else None,
        "verdict": {
            "suite_a": "PASS" if sa["all_zero"] else "FAIL",
            "suite_b": "PASS" if (len(candidates) > 0 and sb["overall_pass"]) else "FAIL" if len(candidates) > 0 else "SKIP",
        },
    }
    receipt_path = "/home/allaun/SilverSight/python/nr_bracket_validation_receipt.json"
    with open(receipt_path, "w") as f:
        json.dump(receipt, f, indent=2)
    print(f"\n  Receipt saved: {receipt_path}")


if __name__ == "__main__":
    main()