#!/usr/bin/env python3
"""Bekenstein Phase-0 spectrometer (panel-greenlit, 2026-07-02).

Extends the Lean-checked traced-marginal computations
(IndisputableMonolith/Holography/SharedCutMarginal.lean, N=2,3 by `decide`)
to arbitrary lattice sizes by exact GF(2) linear algebra, per the
bekenstein_master_plan judge decision:

  "Instrumented Phase 0 ... emit the full marginal(N) and fiber sequence,
   N=1..5+ (Python past Lean's decide limit) ... marginal(N), fiber
   structure, and defect d(N) = 4 - marginal/face ... the Phase-0 defect
   sequence fit vs {1/pi, 2pi, phi^5/pi} plus the gluing law
   capacityDefect(m+n) = capacityDefect m + capacityDefect n + seam."

Mathematical basis (which is why enumeration is unnecessary):
the closed set C is the kernel of the GF(2) constraint matrix, hence a
LINEAR subspace of F_2^V. The marginal of a vertex subset S is the image
of C under the coordinate projection pi_S, itself a linear subspace, so

  marginal_dim(S) = dim C - dim { x in C : x|_S = 0 }.

Fibers of a linear map over its image are kernel cosets, so ALL fibers are
constant automatically (fiber = 2^(dim C - marginal_dim)). The Lean module
verified this by enumeration at N=2,3; here it is structural.

Geometries: 1xN face strips and n x m face patches, open (planar) and
torus (closed-surface, the horizon's topology class). Closure conventions:
GLOBAL (one XOR over all vertices) and LOCAL (one XOR per face).
Subregions: single face, k-face prefixes (strip), a x b sub-patches (2D).

Outputs a plain-text report; every number is exact (integers).
"""

from __future__ import annotations

import math
from fractions import Fraction


# ---------------------------------------------------------------- GF(2) core

def rank_gf2(rows: list[int]) -> int:
    """Rank of a GF(2) matrix given as bitmask rows."""
    r = 0
    pivots: list[int] = []
    for row in rows:
        cur = row
        for p in pivots:
            cur = min(cur, cur ^ p)
        if cur:
            pivots.append(cur)
            # keep pivots reduced (not required for rank, cheap here)
            r += 1
    return r


def kernel_basis_gf2(rows: list[int], nvars: int) -> list[int]:
    """Basis (as bitmasks) of the kernel of the constraint matrix."""
    # Gaussian elimination to row echelon, track pivot columns.
    mat = [r for r in rows if r]
    pivot_col_of_row: list[int] = []
    echelon: list[int] = []
    for row in mat:
        cur = row
        for e, pc in zip(echelon, pivot_col_of_row):
            if (cur >> pc) & 1:
                cur ^= e
        if cur:
            pc = cur.bit_length() - 1
            echelon.append(cur)
            pivot_col_of_row.append(pc)
    pivot_cols = set(pivot_col_of_row)
    free_cols = [c for c in range(nvars) if c not in pivot_cols]
    basis: list[int] = []
    for fc in free_cols:
        vec = 1 << fc
        # back-substitute: for each echelon row, decide pivot value
        for e, pc in sorted(zip(echelon, pivot_col_of_row), key=lambda t: t[1]):
            # value of row e dotted with current vec, excluding pivot column
            dot = bin(e & vec & ~(1 << pc)).count("1") & 1
            if dot:
                vec |= 1 << pc
        basis.append(vec)
    return basis


def marginal_dim(kernel: list[int], subset_mask: int) -> int:
    """dim of the projection of span(kernel) onto the coordinates in subset_mask."""
    projected = [k & subset_mask for k in kernel]
    return rank_gf2(projected)


# ------------------------------------------------------------- lattice builders

def strip_faces(nfaces: int):
    """2 x (nfaces+1) vertex grid; face j has corners (top j, top j+1, bot j+1, bot j).
    Vertex numbering: top row 0..n, bottom row n+1..2n+1  (n = nfaces)."""
    ncols = nfaces + 1
    nverts = 2 * ncols

    def top(j: int) -> int:
        return j

    def bot(j: int) -> int:
        return ncols + j

    faces = []
    for j in range(nfaces):
        faces.append([top(j), top(j + 1), bot(j + 1), bot(j)])
    return nverts, faces


def patch_faces(nrows: int, ncols_faces: int, torus: bool = False):
    """nrows x ncols_faces face grid. Open: (nrows+1) x (ncols+1) vertices.
    Torus: nrows x ncols vertices with wraparound."""
    if torus:
        vr, vc = nrows, ncols_faces
        nverts = vr * vc

        def v(i: int, j: int) -> int:
            return (i % vr) * vc + (j % vc)
    else:
        vr, vc = nrows + 1, ncols_faces + 1
        nverts = vr * vc

        def v(i: int, j: int) -> int:
            return i * vc + j

    faces = []
    for i in range(nrows):
        for j in range(ncols_faces):
            faces.append([v(i, j), v(i, j + 1), v(i + 1, j + 1), v(i + 1, j)])
    return nverts, faces, v


def face_constraint_rows(faces: list[list[int]]) -> list[int]:
    return [sum(1 << c for c in set(f)) for f in faces]


def global_constraint_row(nverts: int) -> list[int]:
    return [(1 << nverts) - 1]


# ------------------------------------------------------------------- reporting

PHI = (1 + 5 ** 0.5) / 2
CONSTANTS = {
    "1/pi": 1 / math.pi,
    "2pi": 2 * math.pi,
    "phi^5/pi": PHI ** 5 / math.pi,
    "1": 1.0,
    "2": 2.0,
    "3": 3.0,
    "4": 4.0,
}


def fit_constant(x: float) -> str:
    best = min(CONSTANTS.items(), key=lambda kv: abs(kv[1] - x))
    return f"nearest {best[0]} (|err| = {abs(best[1] - x):.4f})"


def line(s: str = "") -> None:
    print(s)


def analyze_strip(max_n: int = 12) -> None:
    line("=" * 78)
    line("STRIP 1xN (open, 2x(N+1) vertices). Subregion = k-face prefix.")
    line("=" * 78)
    for convention in ("GLOBAL", "LOCAL"):
        line(f"\n-- closure = {convention} --")
        line(f"{'N':>3} {'dimC':>5} {'S/face(tot)':>11} "
             f"{'marg(1)':>8} {'marg(k=N-1)':>11} {'perface(k=N-1)':>14} "
             f"{'defect d(N-1)':>13}")
        for n in range(2, max_n + 1):
            nverts, faces = strip_faces(n)
            rows = (global_constraint_row(nverts) if convention == "GLOBAL"
                    else face_constraint_rows(faces))
            ker = kernel_basis_gf2(rows, nverts)
            dim_c = len(ker)
            # single face 0 marginal
            m1 = marginal_dim(ker, sum(1 << c for c in faces[0]))
            # prefix of k = n-1 faces (proper subregion)
            k = n - 1
            sub_verts = set()
            for f in faces[:k]:
                sub_verts.update(f)
            mk = marginal_dim(ker, sum(1 << c for c in sub_verts))
            perface = Fraction(mk, k)
            defect = 4 - perface
            line(f"{n:>3} {dim_c:>5} {Fraction(dim_c, n)!s:>11} "
                 f"{m1:>8} {mk:>11} {perface!s:>14} {defect!s:>13}")
        # gluing law on capacity defect D(k) = 4k - marginal(k)
        n = max_n
        nverts, faces = strip_faces(n)
        rows = (global_constraint_row(nverts) if convention == "GLOBAL"
                else face_constraint_rows(faces))
        ker = kernel_basis_gf2(rows, nverts)

        def cap_defect(k: int) -> int:
            sv = set()
            for f in faces[:k]:
                sv.update(f)
            return 4 * k - marginal_dim(ker, sum(1 << c for c in sv))

        line(f"\n  gluing law check (N = {n} host):  D(m+n) - D(m) - D(n) = seam")
        for m in range(1, 5):
            for k2 in range(1, 5):
                if m + k2 >= n:
                    continue
                seam = cap_defect(m + k2) - cap_defect(m) - cap_defect(k2)
                line(f"    D({m}+{k2}) - D({m}) - D({k2}) = {seam}")
        dvals = [4 - Fraction(4 * k - cap_defect(k), k) for k in (n - 1,)]
        line(f"  asymptotic defect/face -> {float(dvals[-1]):.4f}  "
             f"[{fit_constant(float(dvals[-1]))}]")


def analyze_patch(max_n: int = 8) -> None:
    line("")
    line("=" * 78)
    line("2D PATCH n x n faces (open). Subregion = a x a sub-patch (corner).")
    line("=" * 78)
    for convention in ("GLOBAL", "LOCAL"):
        line(f"\n-- closure = {convention} --")
        line(f"{'n':>3} {'faces':>6} {'dimC':>5} {'S/face(tot)':>11} "
             f"{'marg(1)':>8} {'marg(axa)':>10} {'perface':>9} {'defect':>8}")
        for n in range(2, max_n + 1):
            nverts, faces, v = patch_faces(n, n, torus=False)
            rows = (global_constraint_row(nverts) if convention == "GLOBAL"
                    else face_constraint_rows(faces))
            ker = kernel_basis_gf2(rows, nverts)
            dim_c = len(ker)
            m1 = marginal_dim(ker, sum(1 << c for c in faces[0]))
            a = n - 1  # proper (a x a) corner sub-patch
            sub_verts = set()
            for i in range(a):
                for j in range(a):
                    sub_verts.update(faces[i * n + j])
            mk = marginal_dim(ker, sum(1 << c for c in sub_verts))
            nf = a * a
            perface = Fraction(mk, nf)
            defect = 4 - perface
            line(f"{n:>3} {n*n:>6} {dim_c:>5} {Fraction(dim_c, n*n)!s:>11} "
                 f"{m1:>8} {mk:>10} {perface!s:>9} {defect!s:>8}")
        line(f"  (perface column is the JOINT marginal per face of the "
             f"(n-1)x(n-1) corner sub-patch)")


def analyze_torus(max_n: int = 10) -> None:
    line("")
    line("=" * 78)
    line("TORUS n x n faces (closed surface, the horizon topology class).")
    line("Subregion = a x a patch on the torus (a = n-1: proper, all but a frame).")
    line("=" * 78)
    for convention in ("GLOBAL", "LOCAL"):
        line(f"\n-- closure = {convention} --")
        line(f"{'n':>3} {'faces':>6} {'dimC':>5} {'S/face(tot)':>11} "
             f"{'marg(1)':>8} {'marg(axa)':>10} {'perface':>9} {'defect':>8}")
        for n in range(3, max_n + 1):
            nverts, faces, v = patch_faces(n, n, torus=True)
            rows = (global_constraint_row(nverts) if convention == "GLOBAL"
                    else face_constraint_rows(faces))
            ker = kernel_basis_gf2(rows, nverts)
            dim_c = len(ker)
            m1 = marginal_dim(ker, sum(1 << c for c in faces[0]))
            a = n - 1
            sub_verts = set()
            for i in range(a):
                for j in range(a):
                    sub_verts.update(faces[i * n + j])
            mk = marginal_dim(ker, sum(1 << c for c in sub_verts))
            nf = a * a
            perface = Fraction(mk, nf)
            defect = 4 - perface
            line(f"{n:>3} {n*n:>6} {dim_c:>5} {Fraction(dim_c, n*n)!s:>11} "
                 f"{m1:>8} {mk:>10} {perface!s:>9} {defect!s:>8}")


def analyze_per_pixel_sum(max_n: int = 8) -> None:
    line("")
    line("=" * 78)
    line("PER-PIXEL SUM vs JOINT (the PerPixelRecordAdditivity fork), 2D open patch.")
    line("sum_i marg(face_i) vs marg(all faces jointly), whole n x n patch as the")
    line("region inside a larger (n+2) x (n+2) host so the region is proper.")
    line("=" * 78)
    for convention in ("GLOBAL", "LOCAL"):
        line(f"\n-- closure = {convention} --")
        line(f"{'n':>3} {'sum marg':>9} {'perface':>9} {'joint':>7} "
             f"{'perface':>9} {'kappa(sum)':>10} {'kappa(joint)':>12}")
        for n in range(2, max_n + 1):
            host = n + 2
            nverts, faces, v = patch_faces(host, host, torus=False)
            rows = (global_constraint_row(nverts) if convention == "GLOBAL"
                    else face_constraint_rows(faces))
            ker = kernel_basis_gf2(rows, nverts)
            # region = central-ish n x n faces (offset 1,1)
            region_faces = [faces[(i + 1) * host + (j + 1)]
                            for i in range(n) for j in range(n)]
            summ = 0
            for f in region_faces:
                summ += marginal_dim(ker, sum(1 << c for c in set(f)))
            sub_verts = set()
            for f in region_faces:
                sub_verts.update(f)
            joint = marginal_dim(ker, sum(1 << c for c in sub_verts))
            nf = n * n
            line(f"{n:>3} {summ:>9} {Fraction(summ, nf)!s:>9} {joint:>7} "
                 f"{Fraction(joint, nf)!s:>9} {float(summ/nf):>10.3f} "
                 f"{float(joint/nf):>12.3f}")


def main() -> None:
    line("BEKENSTEIN PHASE-0 SPECTROMETER  (exact GF(2) linear algebra)")
    line(f"reference constants: 1/pi = {1/math.pi:.5f}   2pi = {2*math.pi:.5f}   "
         f"phi^5/pi = {PHI**5/math.pi:.5f}")
    line("fibers are constant BY LINEARITY (kernel cosets); Lean verified this by")
    line("enumeration at N=2,3 (SharedCutMarginal.lean). All values exact.\n")
    analyze_strip(max_n=12)
    analyze_patch(max_n=8)
    analyze_torus(max_n=10)
    analyze_per_pixel_sum(max_n=8)
    line("")
    line("=" * 78)
    line("READING GUIDE")
    line("=" * 78)
    line("""
kappa candidates on the table (bits/pixel):
  4    = per-pixel traced marginal, GLOBAL closure  (LEG-A kappa=4)
  3    = per-pixel traced marginal, LOCAL closure   (the nullity reading)
  ~1-2 = joint marginal per face                    (kappa -> 1 fork)
  ->0  = joint marginal per face, LOCAL 2D          (perimeter law)
defect d = 4 - kappa. The judge's Live Bet 1 (defect calibrates 1/pi in G)
requires the defect sequence to approach a transcendental; a clean integer
asymptote with an integer seam gluing law is the kill condition ("structureless
rational") UNLESS the gluing seam itself carries the structure.
""")


if __name__ == "__main__":
    main()