"""RUN-008 - Round 05: the saturation gate.

Round 05 splits "saturation" into two layers. Finite hard-witness families are
complete in the limit: a_Leb = sup Lambda(F) = lim lambda_m. Whether some finite
level already equals a_Leb is a separate question, left OPEN. On the upper side
the round adds cell-local area lower bounds and prune rules for the
candidate-cell branch-and-bound. The theorems are read in the report; here each
formula a computation could get wrong is derived again from its proof's own
steps, measured on this tree's own cells and cell members, and the package's
sanity numbers are rebuilt.

Checks
  A0 anchors       the formulas below occur in Round 05 as written
  A1 citations     its references (all verified before) still match
  X1 external      sec. 10/31: Mishra's test sets and bound, Zeng's limit, as their abstracts state them
  W1 ladder base   sec. 10: lambda_1 = pi/4 and the public gap, as displayed
  W2 two witnesses this tree's own example: Lambda({D, RT}) = pi/8 + sqrt(3)/4, three ways
  W3 achiral       Mishra's three bodies have constant width 1 and a mirror axis, so E(2) = SE(2) for them
  D1 eps_b         sec. 19: eps_b = q + 2 R_b d_M from the proof's two one-sided steps
  D2 LB_outer      sec. 20: from Steiner's formula with Per(U) <= Per(P_b)
  D3 LB_env        sec. 21: the factor r/2 from Minkowski's formula with h >= r; Cauchy; Lipschitz
  D4 Phi_r         sec. 22: kite minus sector; increasing in s
  D5 continuity    sec. 9: Steiner's formula with Per <= 2 pi R, R = 2 A-bar
  N1 cells         sec. 18-19 on this tree's cells: each cover lies in its own cell, members within eps_b of P_b
  N2 Phi measured  sec. 22: Phi_r against the hull area of a sampled disc and a point
  N3 lower bounds  sec. 20-23: every cell bound at most the area of every member found
  N4 perimeter     sec. 9 and 21 on bodies of C_A-bar: Per <= 2 pi R, Area >= (r/2) Per, the bound's increment
  P1 sanity        the package's VALIDATION.json rebuilt from this tree's own cells
  P2 package       the package's self-description of its documents matches the archive
  R1 registry      Round 05's judgments reviewed, first-used statuses defined

Usage: python code/run008_round05_saturation_gate.py [--refresh] [--drill] [--blind CID]
"""
import argparse
import datetime
import json
import math
import re
from fractions import Fraction

import numpy as np
import sympy as sp
from scipy.optimize import minimize_scalar
from scipy.spatial import ConvexHull

import gate
import registry
from corpus import SOURCES, TREE, source_bytes, walk_sources
from run007_round04_cover_compiler import clip, hexagon, perimeter, shoelace, support_of, trig_cover, up

PKG = "main.zip!/AMRAL_LUC_FC_Round_05_Package_v0.1.zip!/"
DOC = PKG + "AMRAL_Lebesgue_Finite_Closure_Round_05_Saturation_Gate_v0.1.md"
CRYSTAL = PKG + "AMRAL_LUC_FC_State_Crystal_Round_05_v0.1.md"
VALIDATION = PKG + "AMRAL_LUC_FC_Round_05_VALIDATION.json"
SNAPSHOT = TREE / "data" / "external" / "round05-citations.json"
LOG = TREE / "data" / "gate-logs" / "run008-round05-saturation-gate.json"
MANIFEST = TREE / "data" / "SOURCE_MANIFEST.v1.json"
ABAR = 0.8440935944
R0 = 0.5                 # every cover holds B_{1/2}(0)
NH = 2 ** 16             # angles for sampled support comparisons

ANCHORS = {
    "W1": [r"\lambda_1=\frac{\pi}{4}\approx0.7853981634.", r"\lambda_3\ge0.8344.", r"\boxed{0.0096935944.}"],
    "X1": [r"-disk；-Reuleauxtriangle；-Reuleauxpentagon；", r"a_{\mathrm{Leb}}=\lim_M\Lambda_M."],
    "D1": [r"\underlineb_j=\max\left\{\frac12,b_j-q\right\}.", r"R_{\mathbfb}=\max_{x\inP_{\mathbfb}}\|x\|.",
           r"\varepsilon_{\mathbfb}:=q+2R_{\mathbfb}d_M."],
    "D2": [r"\operatorname{LB}_{\mathrm{outer}}(\mathbfb)=\operatorname{Area}(P_{\mathbfb})-\operatorname{Per}"
           r"(P_{\mathbfb})\varepsilon_{\mathbfb}-\pi\varepsilon_{\mathbfb}^2."],
    "D3": [r"\ell_{\mathbfb}(\theta)=\max\left\{r,\max_j[\underlineb_j-R\|u_\theta-u_j\|]\right\}.",
           r"\operatorname{LB}_{\mathrm{env}}(\mathbfb)=\frac14\int_0^{2\pi}\ell_{\mathbfb}(\theta)\,d\theta."],
    "D4": [r"\Phi_r(s)=\pir^2+r\sqrt{s^2-r^2}-r^2\arccos\fracrs.",
           r"\operatorname{LB}_{\mathrm{tip}}(\mathbfb)=\max_j\Phi_r(\underlineb_j).",
           r"\operatorname{LB}_{\mathrm{cell}}(\mathbfb)=\max\left\{\pir^2,\operatorname{LB}_{\mathrm{outer}},"
           r"\operatorname{LB}_{\mathrm{env}},\operatorname{LB}_{\mathrm{tip}}\right\}."],
    "D5": [r"\operatorname{Per}(U_K)\le2\piR=4\pi\barA.",
           r"|\Lambda(K)-\Lambda(L)|\le4\pi\barA\varepsilon+\pi\varepsilon^2."],
}
REFS = {1: {"kind": "arxiv", "id": "1502.01251"}, 2: {"kind": "arxiv", "id": "1810.10089"},
        3: {"kind": "arxiv", "id": "2608.30538"}, 4: {"kind": "arxiv", "id": "2609.01284"},
        5: {"kind": "crossref", "doi": "10.1017/CBO9781139003858"}}
CITE_OPTS = {1: {}, 2: {}, 3: {}, 4: {}, 5: {}}
ABSTRACTS = (3, 4)       # sec. 10 and 31 say what Mishra and Zeng proved


# ---------------------------------------------------------------- geometry, this tree's own

def envelope(owner, n, pieces):
    """The upper envelope on the circle of finitely many pieces, each holding one arc of
    it; owner(t) is the maximizing piece at the angles t. Returns (cuts, owners): each
    owner holds from its cut to the next, the last wrapping round. Hand-overs are found by
    bisection on the owner, and every hand-over inside a grid step is followed, so neither
    a piece whose arc is shorter than a step nor a crossing lying exactly on a grid angle
    is lost. (The first version compared two pieces by a sign change of their difference
    and put the cut a step early when the difference vanished on the grid: RUN-009.)
    Angles are taken mod 2 pi before the owner is asked, so the end of the last step is
    the grid's first angle and not a float a hair short of it; a chain of hand-overs that
    does not arrive at the next grid owner is an error, never a silent result."""
    own_at = lambda t: owner(np.mod(t, 2 * np.pi))
    th = 2 * np.pi * np.arange(n) / n
    own = own_at(th)
    idx = np.nonzero(own != np.roll(own, -1))[0]
    if len(idx) == 0:
        return [0.0], [int(own[0])]
    lo, hi, first = th[idx].copy(), th[idx] + 2 * np.pi / n, own[idx]
    for _ in range(50):
        mid = (lo + hi) / 2
        same = own_at(mid) == first
        lo, hi = np.where(same, mid, lo), np.where(same, hi, mid)
    after = own_at(hi)
    cuts, owners = [], []
    for j, i in enumerate(idx):
        cuts.append(float(hi[j]))
        owners.append(int(after[j]))
        end, target = th[i] + 2 * np.pi / n, int(own[(i + 1) % n])
        for _ in range(pieces):
            if owners[-1] == target:
                break
            a, b, cur = cuts[-1], end, owners[-1]
            for _ in range(50):
                m = (a + b) / 2
                if own_at(np.array([m]))[0] == cur:
                    a = m
                else:
                    b = m
            cuts.append(b)
            owners.append(int(own_at(np.array([b]))[0]))
        if owners[-1] != target:
            raise RuntimeError(f"envelope: hand-overs after {th[i]:.6f} do not reach piece {target}")
    return cuts, owners


def segments(cuts, owners):
    """(start, end, owner) for each arc of the envelope"""
    ends = cuts[1:] + [cuts[0] + 2 * np.pi]
    return zip(cuts, ends, owners)


def verify_segments(C, rho, ranges, segs, tol=1e-14):
    """The envelope's segments checked in closed form (added 2026-09-26, RUN-031). Hand-overs are found
    between grid angles, so a piece that wins only inside one grid step, both ends of the step keeping
    one owner, was lost: on one leaf of Round 24's certificate cell a corner of the B7 core stands
    1.7e-6 outside the disc over 5e-3 rad, and the hull came out 1.47e-9 low. Pieces have support
    c.u + rho, and ranges[k] is None or (mid, half), the normals an arc carries. h_j - h_k =
    |c_j - c_k| cos(t - a) + (rho_j - rho_k) peaks at t = a, so its maximum on a segment, within j's
    range, is at an end or at a. A piece that beats the owner by more than tol takes the sub-arc where
    it does, and the new segments are checked in turn, until no piece beats any owner."""
    C, rho = np.asarray(C, dtype=float), np.asarray(rho, dtype=float)

    def windows(j, a, b):
        if ranges[j] is None:
            return [(a, b)]
        mid, half = ranges[j]
        m0 = math.floor((a - mid - half) / (2 * math.pi))
        out = []
        for m in range(m0, m0 + 3):
            lo, hi = mid - half + 2 * math.pi * m, mid + half + 2 * math.pi * m
            if max(lo, a) < min(hi, b):
                out.append((max(lo, a), min(hi, b)))
        return out

    segs, changed = list(segs), True
    while changed:
        changed, nxt = False, []
        for a, b, k in segs:
            best = None
            for j in range(len(C)):
                if j == k:
                    continue
                d = C[j] - C[k]
                A, al, B = math.hypot(d[0], d[1]), math.atan2(d[1], d[0]), rho[j] - rho[k]
                for lo, hi in windows(j, a, b):
                    peaks = [al + 2 * math.pi * m for m in range(math.floor((lo - al) / (2 * math.pi)), math.floor((hi - al) / (2 * math.pi)) + 2)]
                    for t in [lo, hi] + [p for p in peaks if lo < p < hi]:
                        v = A * math.cos(t - al) + B
                        if v > tol and (best is None or v > best[0]):
                            best = (v, j, t, lo, hi, A, al, B)
            if best is None:
                nxt.append((a, b, k))
                continue
            _, j, t, lo, hi, A, al, B = best
            if A == 0:
                s, e = lo, hi
            else:
                beta = math.acos(max(-1.0, min(1.0, -B / A)))
                centre = al + 2 * math.pi * round((t - al) / (2 * math.pi))
                s, e = max(lo, centre - beta), min(hi, centre + beta)
            nxt += [x for x in ((a, s, k), (s, e, j), (e, b, k)) if x[1] > x[0]]
            changed = True
        segs = nxt
    return segs


def hull_disc(pts, r=R0, per_piece=64):
    """Area and perimeter of conv(B_r(0) and the points), exactly from the support
    function h = max(r, max_i p_i.u): Area = 1/2 int (h^2 - h'^2), Per = int h
    (Cauchy), in closed form between the angles where the maximum changes hands."""
    P = np.asarray(pts, dtype=float).reshape(-1, 2)
    rho = np.hypot(P[:, 0], P[:, 1])
    P, rho = P[rho > r], rho[rho > r]
    if len(rho) == 0:
        return math.pi * r * r, 2 * math.pi * r
    al = np.arctan2(P[:, 1], P[:, 0])
    # piece 0 is the disc, piece i + 1 the point i
    owner = lambda t: np.argmax(np.vstack([np.full((1, len(t)), r), rho[:, None] * np.cos(t[None, :] - al[:, None])]),
                                axis=0)
    cuts, owners = envelope(owner, max(4096, per_piece * len(rho)), len(rho) + 1)
    segs = verify_segments(np.vstack([np.zeros((1, 2)), P]), np.concatenate([[r], np.zeros(len(P))]), [None] * (len(P) + 1),
                           segments(cuts, owners))
    twice, per = 0.0, 0.0
    for a, b, k in segs:
        if k == 0:
            twice += r * r * (b - a)
            per += r * (b - a)
        else:
            q, c = rho[k - 1], al[k - 1]
            twice += q * q / 2 * (math.sin(2 * (b - c)) - math.sin(2 * (a - c)))
            per += q * (math.sin(b - c) - math.sin(a - c))
    return twice / 2, per


def support_disc_pts(P, r, th):
    out = np.full(np.shape(th), r, dtype=float)
    for p in np.asarray(P).reshape(-1, 2):
        out = np.maximum(out, p[0] * np.cos(th) + p[1] * np.sin(th))
    return out


def wrap(x):
    return (x + np.pi) % (2 * np.pi) - np.pi


def _F(x):
    """an antiderivative of |sin(x/2)| on [-pi, pi]"""
    return math.copysign(2 - 2 * math.cos(x / 2), x)


def ell_integral(bl, R, M, r=R0):
    """int_0^{2 pi} max{r, max_j [bl_j - R |u_theta - u_j|]} dtheta, piece by piece;
    |u_theta - u_j| = 2 |sin((theta - theta_j)/2)|."""
    tj = 2 * np.pi * np.arange(M) / M

    def owner(t):                         # piece 0 is the floor r, piece j + 1 the tent about u_j
        best, arg = np.full(np.shape(t), r, dtype=float), np.zeros(np.shape(t), dtype=int)
        for j in range(M):
            v = bl[j] - 2 * R * np.abs(np.sin(wrap(t - tj[j]) / 2))
            m = v > best
            best[m], arg[m] = v[m], j + 1
        return arg

    cuts, owners = envelope(owner, max(8192, 64 * M), M + 1)
    total = 0.0
    for a, b, k in segments(cuts, owners):
        if k == 0:
            total += r * (b - a)
        else:
            xa = float(wrap(a - tj[k - 1]))     # a tent holds an arc away from its own antipode
            total += bl[k - 1] * (b - a) - 2 * R * (_F(xa + (b - a)) - _F(xa))
    return total


def ell_rectangle(bl, R, M, n=200000, r=R0):
    """the same integral as an n-point rectangle rule"""
    td = np.linspace(0.0, 2 * np.pi, n, endpoint=False)
    tj = 2 * np.pi * np.arange(M) / M
    ell = np.full(n, r)
    for j in range(M):
        ell = np.maximum(ell, bl[j] - 2 * R * np.abs(np.sin((td - tj[j]) / 2)))
    return float(np.mean(ell) * 2 * np.pi)


def outer(b, M):
    th = 2 * np.pi * np.arange(M) / M
    P = np.array([[-100.0, -100.0], [100.0, -100.0], [100.0, 100.0], [-100.0, 100.0]])
    for t, bb in zip(th, b):
        P = clip(P, np.array([math.cos(t), math.sin(t)]), bb)
    return P


def slab(P, u, lo):
    """the part of P with x.u >= lo; when lo equals P's support in u the slab has
    zero width and clipping can lose it to rounding - then it is the face of P
    with normal u"""
    S = clip(P, -u, -lo)
    if len(S):
        return S
    hp = P @ u
    return P[hp >= hp.max() - 1e-12] if lo <= hp.max() + 1e-12 else S


def closest_point(S):
    """the point of the convex polygon S nearest the origin"""
    best, bd = None, np.inf
    for i in range(len(S)):
        a, d = S[i], S[(i + 1) % len(S)] - S[i]
        L = float(d @ d)
        t = 0.0 if L == 0 else min(1.0, max(0.0, -float(a @ d) / L))
        x = a + t * d
        if x @ x < bd:
            best, bd = x, float(x @ x)
    return best


def reuleaux_vertices(n):
    rho = 1 / (2 * math.cos(math.pi / (2 * n)))       # the diagonal of the regular n-gon is the width 1
    return rho * np.array([[math.cos(2 * math.pi * k / n), math.sin(2 * math.pi * k / n)] for k in range(n)])


def reuleaux_support(n):
    """support function of the Reuleaux n-gon of width 1 with a vertex at angle 0:
    the arc of radius 1 about vertex k carries the normals within pi/(2n) of the
    direction opposite that vertex; elsewhere a vertex supports"""
    V, tk = reuleaux_vertices(n), 2 * np.pi * np.arange(n) / n

    def h(th):
        th = np.asarray(th, dtype=float)
        vert = V @ np.stack([np.cos(th), np.sin(th)])
        out = np.max(vert, axis=0)
        for k in range(n):
            out = np.where(np.abs(wrap(th - tk[k] - np.pi)) <= np.pi / (2 * n), vert[k] + 1.0, out)
        return out
    return h


def reuleaux_boundary(n, per_arc):
    V = reuleaux_vertices(n)
    arcs = []
    for k in range(n):
        a, b = V[(k + (n - 1) // 2) % n] - V[k], V[(k + (n + 1) // 2) % n] - V[k]
        t0 = math.atan2(a[1], a[0])
        ts = t0 + wrap(math.atan2(b[1], b[0]) - t0) * np.linspace(0, 1, per_arc)
        arcs.append(V[k] + np.stack([np.cos(ts), np.sin(ts)], axis=1))
    return np.vstack(arcs)


def disc_samples(n, r=R0):
    t = 2 * np.pi * np.arange(n) / n
    return r * np.stack([np.cos(t), np.sin(t)], axis=1)


# test covers for the cells: all hold B_{1/2}(0). The heptagon circumscribed about the disc
# (area 7 tan(pi/7)/4 = 0.842756 <= A-bar) has long straight edges off every grid; Pal's
# hexagon (area sqrt(3)/2) does too but lies outside C_A-bar, which no bound's proof uses
def _disc_rt():
    hrt = reuleaux_support(3)
    return (lambda t: np.maximum(R0, hrt(t))), math.pi / 8 + math.sqrt(3) / 4


def _heptagon():
    V = (R0 / math.cos(math.pi / 7)) * np.array([[math.cos((2 * k + 1) * math.pi / 7), math.sin((2 * k + 1) * math.pi / 7)]
                                                  for k in range(7)])
    return (lambda t: support_disc_pts(V, R0, np.asarray(t, dtype=float))), 7 * R0 * R0 * math.tan(math.pi / 7)


def test_covers():
    covers = {"sec33": trig_cover(0.515, {2: (0.01, 0.0)}),
              "asym": trig_cover(0.517, {2: (0.008 * math.cos(0.3), -0.008 * math.sin(0.3)),
                                         3: (0.004 * math.cos(1.1), -0.004 * math.sin(1.1))}),
              "disc+RT": _disc_rt(), "heptagon": _heptagon(), "hexagon": hexagon()}
    return covers


GRIDS = ((8, 1 / 64), (16, 1 / 256), (32, 1 / 1024), (64, 1 / 4096))


def make_cell(h, area_U, M, q, quantize, b_lower, rng):
    th = 2 * np.pi * np.arange(M) / M
    U = np.stack([np.cos(th), np.sin(th)], axis=1)
    hj = h(th)
    b = quantize(hj, q)
    bl = b_lower(b, q)
    P = outer(b, M)
    Rb = float(np.max(np.linalg.norm(P, axis=1)))
    tt = 2 * np.pi * np.arange(NH) / NH
    hP, hU = support_of(P, tt), h(tt)
    cell = {"M": M, "q": q, "b": b, "b_lower": bl, "area_P": shoelace(P), "per_P": perimeter(P), "R_b": Rb,
            "cover_in_cell": bool(np.all(bl - 1e-12 <= hj) and np.all(hj <= b + 1e-12)),
            "cover_pokes_out": float(np.max(hU - hP)), "members": []}
    if cell["cover_in_cell"]:
        cell["members"].append({"name": "cover", "area": area_U, "dH": float(np.max(hP - hU)), "in_cell": True})
    cell["members"].append({"name": "P_b", "area": cell["area_P"], "dH": 0.0, "in_cell": True})

    def member(pts, label):
        hm_j = support_disc_pts(pts, R0, th)
        return {"name": label, "area": hull_disc(pts)[0], "dH": float(np.max(hP - support_disc_pts(pts, R0, tt))),
                "in_cell": bool(np.all(bl - 1e-12 <= hm_j) and np.all(hm_j <= b + 1e-12))}

    # members conv(B_{1/2} and x_j) with x_j in the slab P_b and x.u_j >= b_lower_j
    slabs = [slab(P, U[j], bl[j]) for j in range(M)]
    if all(len(s) for s in slabs):
        inner = np.array([closest_point(s) for s in slabs])
        cell["members"].append(member(inner, "inner"))
        for k in range(8):
            cell["members"].append(member(np.array([rng.dirichlet(np.ones(len(s))) @ s for s in slabs]), f"random-{k}"))
        if M <= 16:              # a small descent for the smallest member
            pts, a = inner.copy(), hull_disc(inner)[0]
            for _ in range(300):
                j = int(rng.integers(M))
                cand = pts.copy()
                cand[j] = (rng.dirichlet(np.ones(len(slabs[j]))) @ slabs[j] if rng.random() < 0.5
                           else 0.5 * (pts[j] + closest_point(slabs[j])))
                a2 = hull_disc(cand)[0]
                if a2 < a:
                    pts, a = cand, a2
            cell["members"].append(member(pts, "descent"))
    cell["int_ell_Rb"] = ell_integral(bl, Rb, M)
    cell["int_ell_Rglobal"] = ell_integral(bl, 2 * ABAR, M)
    return cell


_CELLS, _PACKAGE = {}, {}


def cells(w):
    key = (id(w["quantize"]), id(w["b_lower"]))
    if key not in _CELLS:
        out = []
        for ci, (name, (h, aU)) in enumerate(test_covers().items()):
            for M, q in GRIDS:
                c = make_cell(h, aU, M, q, w["quantize"], w["b_lower"], np.random.default_rng(20260924 + 100 * ci + M))
                c["cover"] = name
                out.append(c)
        _CELLS[key] = out
    return _CELLS[key]


def package_cells(w):
    """the package's five cells: the sec. 33 candidate of Round 04 at M = 64..1024, q = 1/M^2"""
    key = (id(w["quantize"]), id(w["b_lower"]))
    if key not in _PACKAGE:
        h, _ = test_covers()["sec33"]
        out = []
        for M in (64, 128, 256, 512, 1024):
            q = 1 / (M * M)
            th = 2 * np.pi * np.arange(M) / M
            b = w["quantize"](h(th), q)
            bl = w["b_lower"](b, q)
            P = outer(b, M)
            Rb = float(np.max(np.linalg.norm(P, axis=1)))
            # distance of (h - 1/2)/q from the nearest integer: how close a height sits to a level
            lev = (h(th) - 0.5) / q
            out.append({"M": M, "q": q, "vertices": len(P), "area": shoelace(P), "per": perimeter(P), "R_b": Rb,
                        "b_lower": bl, "int_rect": ell_rectangle(bl, Rb, M), "int_exact": ell_integral(bl, Rb, M),
                        "nearest_level_gap": float(np.min(np.abs(lev - np.round(lev))))})
        _PACKAGE[key] = out
    return _PACKAGE[key]


# ---------------------------------------------------------------- the round's formulas, as written

def _sym(*xs):
    return any(isinstance(x, sp.Basic) for x in xs)


def eps_b_text(q, R, M):
    """sec. 19: eps_b = q + 2 R_b d_M, with Round 04's d_M = 2 sin(pi/(2M))"""
    sin, pi = (sp.sin, sp.pi) if _sym(q, R, M) else (math.sin, math.pi)
    return q + 2 * R * (2 * sin(pi / (2 * M)))


def lb_outer_text(area, per, e):
    pi = sp.pi if _sym(area, per, e) else math.pi
    return area - per * e - pi * e ** 2


def phi_text(r, s):
    sqrt, acos, pi = (sp.sqrt, sp.acos, sp.pi) if _sym(r, s) else (math.sqrt, math.acos, math.pi)
    return pi * r ** 2 + r * sqrt(s ** 2 - r ** 2) - r ** 2 * acos(r / s)


def cont_bound_text(Abar, e):
    pi = sp.pi if _sym(Abar, e) else math.pi
    return 4 * pi * Abar * e + pi * e ** 2


def b_lower_text(b, q):
    return np.maximum(0.5, b - q)


CLAIMS = {
    "lambda1_display": "0.7853981634", "gap_display": "0.0096935944",
    "lower_anchor": "0.8344", "upper_anchor": "0.8440935944",
    "mishra_family": ("disk", "Reuleaux triangle", "Reuleaux pentagon"),
    "witnesses": {"disc": lambda t: np.full(np.shape(t), R0), "Reuleaux triangle": reuleaux_support(3),
                  "Reuleaux pentagon": reuleaux_support(5)},
    "phi": phi_text, "eps_b": eps_b_text, "lb_outer": lb_outer_text, "env_factor": sp.Rational(1, 4),
    "cont_bound": cont_bound_text, "quantize": up, "b_lower": b_lower_text,
}


# ---------------------------------------------------------------- checks

def x1_external(w):
    out, ab = [], w["snapshot"]["abstracts"]
    mishra = ab["3"]["abstract"]
    sentence = re.search(r"[^.]*test sets[^.]*\.", mishra)
    named = sentence.group(0).lower().replace("disk", "disc") if sentence else ""
    for body in w["mishra_family"]:
        if body.lower().replace("disk", "disc") not in named:
            out.append(f"Mishra's abstract does not list {body!r} among its test sets")
    if not re.search(r"area at least 0\.8344\b", mishra):
        out.append("Mishra's abstract does not state the bound 0.8344")
    if r"a_{\mathrm{Leb}}=\lim_{M\to\infty}\Lambda_M" not in gate.squash(ab["4"]["abstract"]):
        out.append("Zeng's abstract does not state a_Leb = lim Lambda_M")
    return out


def w1_ladder_base(w):
    out = []
    if f"{math.pi / 4:.10f}" != w["lambda1_display"]:
        out.append(f"pi/4 = {math.pi / 4:.12f}, displayed {w['lambda1_display']}")
    if abs(hull_disc(np.zeros((0, 2)))[0] - math.pi / 4) > 1e-15:
        out.append("the disc alone does not have area pi/4")
    lo, hi = Fraction(w["lower_anchor"]), Fraction(w["upper_anchor"])
    if hi - lo != Fraction(w["gap_display"]):
        out.append(f"gap is {float(hi - lo)!r}, displayed {w['gap_display']}")
    if not (math.pi / 4 < lo < hi):
        out.append("lambda_1 <= 0.8344 <= 0.8440935944 fails")
    return out


def rt_hull_area(t, n_disc, per_arc):
    pts = np.vstack([disc_samples(n_disc), reuleaux_boundary(3, per_arc) + np.asarray(t)])
    return ConvexHull(pts).volume


def w2_two_witnesses(w):
    out = []
    exact = math.pi / 8 + math.sqrt(3) / 4
    s = 1 / math.sqrt(3)           # the Reuleaux triangle's vertices sit at 1/sqrt(3) from its centre
    via_phi = math.pi * R0 ** 2 + 3 * (float(w["phi"](R0, s)) - math.pi * R0 ** 2)
    if abs(via_phi - exact) > 1e-12:
        out.append(f"disc + three tips via Phi = {via_phi!r}, pi/8 + sqrt(3)/4 = {exact!r}")
    dense = rt_hull_area((0.0, 0.0), 200000, 100000)
    if abs(dense - exact) > 1e-8:
        out.append(f"hull of disc and Reuleaux triangle = {dense!r}, not pi/8 + sqrt(3)/4")
    # the concentric hull is the disc and three tips: the arcs leave each vertex along lines
    # tangent to the disc (the tangent at A of the arc about B is at distance A.(A - B) = 1/2)
    V = reuleaux_vertices(3)
    tangent = [abs(float(V[0] @ (V[0] - V[k]))) for k in (1, 2)]
    if max(abs(x - R0) for x in tangent) > 1e-15:
        out.append(f"arc tangents at a vertex lie at {tangent} from the centre, not 1/2")
    # the lower bound's two numerical facts: a vertex 0.7 out already costs more than the value,
    # and inside 0.7 two vertices a width apart are more than two tip half-widths apart
    far = float(w["phi"](R0, 0.7))
    if not far > exact:
        out.append(f"Phi(0.7) = {far!r} does not exceed the value")
    if not 2 * math.acos(R0 / 0.7) < math.pi / 2:
        out.append("tips of vertices within 0.7 can overlap")
    base = rt_hull_area((0.0, 0.0), 20000, 5000)
    worst = min(rt_hull_area((d * math.cos(a), d * math.sin(a)), 20000, 5000) - base
                for d in (0.002, 0.01, 0.03, 0.1) for a in np.linspace(0, 2 * math.pi, 24, endpoint=False))
    if worst < -1e-9:
        out.append(f"an offset placement is smaller by {-worst:.2e}")
    w["w2"] = {"exact": exact, "via_phi": via_phi, "hull_dense": dense, "phi_at_0.7": far,
               "two_tip_half_widths_at_0.7_deg": math.degrees(2 * math.acos(R0 / 0.7)),
               "least_offset_increase": worst}
    return out


def mirror_residual(h, n=4096):
    th = 2 * np.pi * np.arange(n) / n
    res = lambda a: float(np.max(np.abs(h(2 * a - th) - h(th))))
    grid = np.linspace(0, np.pi, 721, endpoint=False)
    a0 = grid[int(np.argmin([res(a) for a in grid]))]
    fine = minimize_scalar(res, bounds=(a0 - np.pi / 720, a0 + np.pi / 720), method="bounded",
                           options={"xatol": 1e-12})
    return min(res(a0), fine.fun)


def w3_achiral(w):
    out, rows = [], {}
    th = 2 * np.pi * np.arange(8192) / 8192
    for name, h in w["witnesses"].items():
        width = float(np.max(np.abs(h(th) + h(th + np.pi) - 1)))
        mirror = mirror_residual(h)
        rows[name] = {"width_defect": width, "mirror_residual": mirror}
        if width > 1e-12:
            out.append(f"{name}: not of constant width 1 ({width:.2e})")
        if mirror > 1e-9:
            out.append(f"{name}: no mirror axis (residual {mirror:.2e}) - reflections would matter")
    w["w3"] = rows
    return out


q_, R_, M_, bj_ = sp.symbols("q R M b", positive=True)


def d1_eps(w):
    out = []
    d = 2 * sp.sin(sp.pi / (2 * M_))
    upper = bj_ + R_ * d               # x.v <= x.u_j + |x||v - u_j| <= b_j + R_b d_M
    lower = bj_ - q_ - R_ * d          # h_U(v) >= h_U(u_j) - R_b d_M >= b_j - q - R_b d_M
    if sp.simplify(w["eps_b"](q_, R_, M_) - (upper - lower)) != 0:
        out.append("eps_b is not (b_j + R_b d_M) - (b_j - q - R_b d_M)")
    # d_M is the farthest a direction sits from the grid
    for M in (8, 64, 1024):
        t = np.linspace(0, 2 * np.pi / M, 100001)
        far = float(np.max(np.minimum(2 * np.abs(np.sin(t / 2)), 2 * np.abs(np.sin((t - 2 * np.pi / M) / 2)))))
        if abs(far - 2 * math.sin(math.pi / (2 * M))) > 1e-12:
            out.append(f"M={M}: farthest direction {far!r} vs d_M")
    return out


AP_, PerP_, e_ = sp.symbols("A_P P e", positive=True)


def d2_outer(w):
    steiner = lambda per, eps: per * eps + sp.pi * eps ** 2        # Area(U + eps B) - Area(U)
    derived = AP_ - steiner(PerP_, e_)                             # with Per(U) <= Per(P_b)
    return [] if sp.simplify(w["lb_outer"](AP_, PerP_, e_) - derived) == 0 else ["LB_outer is not Area(P) - (Per(P) eps + pi eps^2)"]


def d3_env(w):
    out = []
    r = sp.Rational(1, 2)
    # Minkowski: Area(U) = 1/2 int h dS_U >= (r/2) S_U(S^1) = (r/2) Per(U) when h >= r
    if sp.nsimplify(w["env_factor"]) != r / 2:
        out.append(f"LB_env factor {w['env_factor']} is not r/2 = 1/4")
    rng = np.random.default_rng(58)
    for _ in range(200):
        P = rng.normal(size=(int(rng.integers(3, 9)), 2))
        R = float(np.max(np.linalg.norm(P, axis=1)))
        a, b = rng.uniform(0, 2 * np.pi, 2)
        ha, hb = support_of(P, np.array([a])), support_of(P, np.array([b]))
        if abs(ha - hb)[0] > R * 2 * abs(math.sin((a - b) / 2)) + 1e-12:
            out.append("a support function is not R-Lipschitz on the circle")
            break
    # Cauchy: Per = int h, on random polygons
    for _ in range(20):
        P = rng.normal(size=(int(rng.integers(3, 9)), 2))
        P = P[ConvexHull(P).vertices]
        t = 2 * np.pi * np.arange(400000) / 400000
        if abs(float(np.mean(support_of(P, t)) * 2 * np.pi) - perimeter(P)) > 1e-8:
            out.append("Cauchy's formula fails on a polygon")
            break
    return out


r_, s_ = sp.symbols("r s", positive=True)


def d4_phi(w):
    out = []
    kite = 2 * (r_ * sp.sqrt(s_ ** 2 - r_ ** 2) / 2)       # two right triangles, legs r and sqrt(s^2 - r^2)
    sector = r_ ** 2 * sp.acos(r_ / s_)                     # (1/2) r^2 times the angle 2 arccos(r/s)
    derived = sp.pi * r_ ** 2 + kite - sector
    stated = w["phi"](r_, s_)
    slope = sp.diff(stated, s_) - r_ * sp.sqrt(s_ ** 2 - r_ ** 2) / s_
    for rv, sv in ((0.5, 0.5001), (0.5, 0.57735), (0.5, 1.2), (0.3, 0.9), (1.0, 3.0)):
        at = {r_: sp.Float(rv, 50), s_: sp.Float(sv, 50)}
        if abs(sp.N((stated - derived).subs(at), 50)) > 1e-40:
            out.append(f"Phi({rv}, {sv}) differs from disc + kite - sector")
        if abs(sp.N(slope.subs(at), 50)) > 1e-40:
            out.append(f"dPhi/ds at ({rv}, {sv}) is not r sqrt(s^2 - r^2)/s > 0")
    return out


Ab_, ep_ = sp.symbols("Abar e", positive=True)


def d5_continuity(w):
    R = Ab_ / sp.Rational(1, 2)                  # Round 04: U inside B_{A/r}, r = 1/2
    derived = ep_ * (2 * sp.pi * R) + sp.pi * ep_ ** 2   # Steiner increment with Per(U_K) <= Per(B_R)
    return [] if sp.simplify(w["cont_bound"](Ab_, ep_) - derived) == 0 else ["sec. 9 bound is not eps 2 pi R + pi eps^2"]


def n1_cells(w):
    out, rows = [], []
    for c in cells(w):
        tag = f"{c['cover']} M={c['M']}"
        if not c["cover_in_cell"]:
            out.append(f"{tag}: the cover is not in the cell of its own quantized supports")
        if c["cover_pokes_out"] > 1e-12:
            out.append(f"{tag}: the cover pokes out of P_b by {c['cover_pokes_out']:.2e}")
        e = float(w["eps_b"](c["q"], c["R_b"], c["M"]))
        margin = 2 * c["R_b"] * math.pi / NH                 # sampled support gap, Lipschitz half-step
        for m in c["members"]:
            if not m["in_cell"]:
                out.append(f"{tag}: member {m['name']} left the cell")
            elif m["dH"] + margin > e:
                out.append(f"{tag}: member {m['name']} is {m['dH']:.4f} from P_b > eps_b {e:.4f}")
        dmax = max(m["dH"] for m in c["members"])
        rows.append({"cover": c["cover"], "M": c["M"], "q": c["q"], "R_b": c["R_b"], "eps_b": e,
                     "eps_round04": float(eps_round04(c["q"], c["M"])),
                     "max_member_dH": dmax, "members": len(c["members"])})
    w["n1_rows"] = rows
    return out


def eps_round04(q, M):
    D, R = math.pi / M, 2 * ABAR
    return q + (R + (R + q) / math.cos(D)) * 2 * math.sin(D / 2)


def n2_phi_measured(w):
    out, rows, prev = [], {}, -1.0
    disc = disc_samples(400000)
    for s in (0.5, 0.505, 0.525, 1 / math.sqrt(3), 0.7, 1.0, 2 * ABAR):
        hull = ConvexHull(np.vstack([disc, [[s, 0.0]]])).volume
        stated = float(w["phi"](R0, s))
        rows[f"{s:.6f}"] = {"hull": hull, "phi": stated}
        if abs(hull - stated) > 1e-9:
            out.append(f"s={s:.4f}: hull area {hull:.12f}, Phi {stated:.12f}")
        if stated <= prev:
            out.append(f"Phi not increasing at s={s:.4f}")
        prev = stated
    w["n2_rows"] = rows
    return out


def lower_bounds(w, c):
    e = float(w["eps_b"](c["q"], c["R_b"], c["M"]))
    k = float(w["env_factor"])
    lb = {"outer": float(w["lb_outer"](c["area_P"], c["per_P"], e)), "env": k * c["int_ell_Rb"],
          "env_global_R": k * c["int_ell_Rglobal"], "tip": max(float(w["phi"](R0, x)) for x in c["b_lower"])}
    lb["cell"] = max(math.pi * R0 ** 2, lb["outer"], lb["env"], lb["tip"])
    return lb


def n3_lower_bounds(w):
    out, rows = [], []
    for c in cells(w):
        lb = lower_bounds(w, c)
        live = [m for m in c["members"] if m["in_cell"]]
        amin = min(m["area"] for m in live)
        for name, v in lb.items():
            if v > amin + 1e-12:
                out.append(f"{c['cover']} M={c['M']}: LB_{name} = {v:.6f} > smallest member area {amin:.6f}")
        e = float(w["eps_b"](c["q"], c["R_b"], c["M"]))
        rows.append({"cover": c["cover"], "M": c["M"], **{f"LB_{k}": v for k, v in lb.items()},
                     "smallest_member": amin, "smallest_is": min(live, key=lambda m: m["area"])["name"],
                     "P_b_over_smallest": c["area_P"] - amin, "pi_eps_b_squared": math.pi * e * e})
    w["n3_rows"] = rows
    return out


def n4_perimeter(w):
    out, rng = [], np.random.default_rng(20260925)
    bodies = []
    for n in (7, 8, 12, 24):            # regular polygons circumscribed about B_{1/2} with area <= A-bar
        V = (R0 / math.cos(math.pi / n)) * np.array([[math.cos((2 * k + 1) * math.pi / n), math.sin((2 * k + 1) * math.pi / n)]
                                                       for k in range(n)])
        bodies.append((f"tangential-{n}", V))
    bodies.append(("disc", np.zeros((0, 2))))
    for i in range(300):
        k = int(rng.integers(1, 9))
        ang = rng.uniform(0, 2 * np.pi, k)
        bodies.append((f"random-{i}", rng.uniform(0.5, 0.8, k)[:, None] * np.stack([np.cos(ang), np.sin(ang)], axis=1)))
    ratio, used = 0.0, 0
    for name, pts in bodies:
        a, p = hull_disc(pts)
        if a > ABAR:
            continue
        used += 1
        R = max([R0] + [float(np.linalg.norm(x)) for x in pts])
        if p > 2 * math.pi * R + 1e-12 or R > 2 * a + 1e-12:
            out.append(f"{name}: Per <= 2 pi R <= 4 pi A fails")
        if a < R0 / 2 * p - 1e-12:
            out.append(f"{name}: Area {a:.6f} < (r/2) Per {R0 / 2 * p:.6f}")
        for e in (1e-3, 1e-2, 1e-1):
            if e * p + math.pi * e * e > float(w["cont_bound"](ABAR, e)) + 1e-12:
                out.append(f"{name}: Steiner increment at eps={e} exceeds the sec. 9 bound")
                break
        ratio = max(ratio, p / (4 * ABAR))
    # this tree's exact hull areas against sampled hulls
    for i in range(5):
        pts = bodies[5 + i][1]
        a = hull_disc(pts)[0]
        s = ConvexHull(np.vstack([disc_samples(400000), pts])).volume
        if abs(a - s) > 1e-8:
            out.append(f"exact hull area {a!r} vs sampled {s!r}")
    tang7 = hull_disc(bodies[0][1])
    w["n4"] = {"bodies_in_C_Abar": used, "max_Per_over_4Abar": ratio,
               "heptagon_area": tang7[0], "heptagon_Per_over_Area": tang7[1] / tang7[0]}
    return out


_TOOL = {}


def integrator_check():
    """this tree's exact envelope integral against a closed form (one tent, b = 0.6,
    R = 1/2) and against a 2^22-point rectangle rule on the package's M = 64 cell"""
    if not _TOOL:
        a = 2 * math.asin(0.1)
        closed = math.pi + 0.2 * a - 4 * (1 - math.cos(a / 2))
        h, _ = test_covers()["sec33"]
        M, q = 64, 1 / 4096
        b = up(h(2 * np.pi * np.arange(M) / M), q)
        bl = b_lower_text(b, q)
        Rb = float(np.max(np.linalg.norm(outer(b, M), axis=1)))
        exact = ell_integral(bl, Rb, M)
        _TOOL.update({"one_tent_minus_closed_form": ell_integral(np.array([0.6]), 0.5, 1) - closed,
                      "rectangle_2^22_minus_exact": ell_rectangle(bl, Rb, M, n=2 ** 22) - exact})
    return _TOOL


def p1_sanity(w):
    out, v = [], w["validation"]["sanity_check"]
    tool = integrator_check()
    if abs(tool["one_tent_minus_closed_form"]) > 1e-14 or abs(tool["rectangle_2^22_minus_exact"]) > 1e-10:
        out.append(f"this tree's exact integral fails its own checks: {tool}")
    w["p1_tool"] = dict(tool)
    h, A = test_covers()["sec33"]
    fixed = {"true_test_candidate_area": A, "disk_only_lambda1": math.pi / 4, "current_public_lower_anchor": 0.8344,
             "current_public_upper_anchor": 0.8440935944, "public_gap": 0.8440935944 - 0.8344}
    for k, mine in fixed.items():
        if abs(v[k] - mine) > 1e-15:
            out.append(f"{k}: recorded {v[k]!r}, recomputed {mine!r}")
    rows = []
    for rec, pc in zip(v["cell_lower_bound_convergence"], package_cells(w)):
        M, q = pc["M"], pc["q"]
        if rec["M"] != M or rec["q"] != q or rec["vertices"] != pc["vertices"]:
            out.append(f"M={M}: grid, step or vertex count differ")
        e = float(w["eps_b"](q, pc["R_b"], M))
        k = float(w["env_factor"])
        mine = {"area_polygon": pc["area"], "perimeter_polygon": pc["per"], "R_cell": pc["R_b"], "epsilon_cell": e,
                "LB_outer": float(w["lb_outer"](pc["area"], pc["per"], e)), "LB_env": k * pc["int_rect"]}
        mine["LB_combined"] = max(math.pi * R0 ** 2, mine["LB_outer"], mine["LB_env"])
        for key, val in mine.items():
            if abs(rec[key] - val) > 1e-12:
                out.append(f"M={M} {key}: recorded {rec[key]!r}, rebuilt {val!r}")
        tip = max(float(w["phi"](R0, x)) for x in pc["b_lower"])
        rows.append({"M": M, "LB_env_rectangle": k * pc["int_rect"], "LB_env_exact": k * pc["int_exact"],
                     "rectangle_minus_exact": k * (pc["int_rect"] - pc["int_exact"]), "LB_tip": tip,
                     "LB_cell_text": max(math.pi * R0 ** 2, mine["LB_outer"], k * pc["int_exact"], tip),
                     "eps_round04": eps_round04(q, M), "nearest_level_gap": pc["nearest_level_gap"]})
    safe = all(r["LB_combined"] <= A + 1e-8 for r in v["cell_lower_bound_convergence"])
    if v["all_cell_lower_bounds_safe_for_test_candidate"] != safe:
        out.append("the recorded safety flag differs")
    w["p1_rows"] = rows
    return out


def p2_package(w):
    out, leaves = [], {row[0]: row for row in w["manifest"]["leaves"]}
    for rec, leaf, text in ((w["validation"]["round05"], DOC, w["text"]),
                            (w["validation"]["state_crystal"], CRYSTAL, w["crystal"])):
        _, size, digest = leaves[leaf]
        if rec["sha256"] != digest or rec["bytes"] != size:
            out.append(f"{rec['file']}: recorded sha256/bytes differ from the archive")
        n = text.count("$$")
        if rec["double_dollar_count"] != n or rec["double_dollar_even"] != (n % 2 == 0):
            out.append(f"{rec['file']}: recorded {rec['double_dollar_count']} $$, document has {n}")
        if rec["contains_unicode_escape_sequence"] != bool(re.search(r"\\u[0-9a-fA-F]{4}", text)):
            out.append(f"{rec['file']}: unicode-escape flag differs")
        if rec["contains_bracket_math_delimiter"] != bool(re.search(r"\\[\[\]()]", text)):
            out.append(f"{rec['file']}: bracket-delimiter flag differs")
    return out


# ---------------------------------------------------------------- registry (Round 05)

R05_DEFINITIONS = {
    "FORMALIZED": ("contextual", "sec. 0 and sec. 38: the saturation theory is written out as definitions and proved "
                                 "statements; nothing is computed, and finite attainment itself stays OPEN"),
    "FALSE-POSITIVE": ("not-a-status", "sec. 35: one of five witness-ancestry labels (REUSED, REDUCED, EXTENSION, "
                                       "NOVEL, FALSE-POSITIVE) for Round 06's witness records; the census meets it "
                                       "as a status only in Round 06"),
}
PBCT = "PASS by cited theorem"
R05_REVIEWS = [
    ("PROVED", "Finite witness lower theorem", "analytic", "PASS", "Thm 3.1: a universal cover holds every placed hull"),
    ("PROVED", "Finite witness lower value attainment", "analytic", PBCT,
     "sec. 4: C_A-bar compact (Round 04 sec. 5, i.e. Blaschke selection), fit sets closed (Thm 5.1), area continuous"),
    ("PROVED", "Finite-witness completeness", "analytic", PBCT,
     "Thm 6.1: finite intersection property in the compact C_A, then Round 01's constant-width reduction; uses sec. 4"),
    ("PROVED", "Cardinality hierarchy convergence", "analytic", "PASS", "sec. 7: immediate from Thm 6.1"),
    ("PROVED", "Level attainment", "analytic", PBCT,
     "sec. 8-9: Steiner's formula and Per <= 2 pi R (D5); the bound holds but is loose by pi, since "
     "Per <= 4 Area <= 4 A-bar and the constant 4 is nearly attained in C_A-bar (N4)"),
    ("PROVED", "No false permanent plateau", "analytic", "PASS",
     "Thm 11.1 and sec. 12; the statement already follows from sec. 7, the proof adds the finite separating batch"),
    ("PROVED", "Exact finite saturation criterion", "analytic", "PASS", "Thm 14.1"),
    ("PROVED", "Candidate-cell area lower bound", "analytic", PBCT,
     "sec. 18-23: Steiner, Cauchy and Minkowski formulas (D1-D4); safe on every member of this tree's cells (N1-N3)"),
    ("PROVED", "Outer non-universal prune", "analytic", "PASS", "sec. 24 P1: U inside P_b (Round 04 Cor 8.2)"),
    ("PROVED", "Strict non-universal eventual prune", "analytic", "PASS", "sec. 26 with Round 04 Thm 17.1"),
    ("PROVED", "Strict area-suboptimal eventual prune", "analytic", "PASS", "sec. 27: LB_outer tends to Area(U)"),
    ("CLOSED", "Convergence vs attainment separation", "analytic", "PASS", "sec. 1: floor(n alpha)/n, alpha irrational"),
    ("COMPUTE-DEFERRED", "Finite Witness Exchange Compiler and Minimizer-Set Separation", "n/a", "n/a",
     "sec. 35: a conditional instruction for Round 06, not a judgment on this round"),
    ("OPEN", "FINITE ATTAINMENT", "n/a", "n/a", "sec. 0; = sec. 32's conjecture"),
    ("OPEN", "CROSS-MAP", "n/a", "n/a", "sec. 31: Zeng's Lambda_M against lambda_m"),
    ("OPEN", "Exact finite witness attainment", "n/a", "n/a", ""),
    ("FORMALIZED", "SATURATION GATE", "reading", "PASS", "every theorem of sec. 2-28 read and checked"),
    ("NONE", "New numerical bound", "reading", "PASS", "none claimed (sec. 38)"),
]


def r1_registry(w):
    return registry.completeness(w["registry"], "R05", w["definitions"], w["reviews"])


# ---------------------------------------------------------------- world, checks, drill

CHECKS = {"A0": gate.anchors_check(ANCHORS), "A1": gate.citations_check(CITE_OPTS), "X1": x1_external,
          "W1": w1_ladder_base, "W2": w2_two_witnesses, "W3": w3_achiral,
          "D1": d1_eps, "D2": d2_outer, "D3": d3_env, "D4": d4_phi, "D5": d5_continuity,
          "N1": n1_cells, "N2": n2_phi_measured, "N3": n3_lower_bounds, "N4": n4_perimeter,
          "P1": p1_sanity, "P2": p2_package, "R1": r1_registry}
DEPENDS = {"A0": {"text"}, "A1": {"text", "snapshot"}, "X1": {"snapshot", "mishra_family"},
           "W1": {"lambda1_display", "gap_display", "lower_anchor", "upper_anchor"}, "W2": {"phi"},
           "W3": {"witnesses"}, "D1": {"eps_b"}, "D2": {"lb_outer"}, "D3": {"env_factor"}, "D4": {"phi"},
           "D5": {"cont_bound"}, "N1": {"eps_b", "quantize", "b_lower"}, "N2": {"phi"},
           "N3": {"eps_b", "lb_outer", "env_factor", "phi", "quantize", "b_lower"}, "N4": {"cont_bound"},
           "P1": {"validation", "eps_b", "lb_outer", "env_factor", "quantize", "b_lower"},
           "P2": {"validation", "text", "crystal", "manifest"}, "R1": {"registry", "definitions", "reviews"}}


def load_inputs():
    blobs = {n: source_bytes(n) for n in SOURCES}
    want = {DOC, CRYSTAL, VALIDATION}
    got = {leaf.id: d for leaf, d in walk_sources(blobs, keep=lambda l: l.id in want) if d is not None}
    return (got[DOC].decode("utf-8"), got[CRYSTAL].decode("utf-8"), json.loads(got[VALIDATION]),
            json.loads(MANIFEST.read_text(encoding="utf-8")))


def package_assertion_power(w):
    """which wrong formulas the package's own assertion (every LB_combined <= the
    candidate's area) would notice on its five cells"""
    A = test_covers()["sec33"][1]
    tests = {
        "LB_outer without the Per*eps term": lambda pc, e: pc["area"] - math.pi * e * e,
        "eps_b with one R_b d_M": lambda pc, e: lb_outer_text(pc["area"], pc["per"], e - pc["R_b"] * 2 * math.sin(math.pi / (2 * pc["M"]))),
        "LB_env with factor 1/2": lambda pc, e: 0.5 * pc["int_rect"],
        "LB = Area(P_b)": lambda pc, e: pc["area"],
    }
    res = {}
    for name, lb in tests.items():
        vals = [lb(pc, eps_b_text(pc["q"], pc["R_b"], pc["M"])) for pc in package_cells(w)]
        res[name] = {"noticed": any(x > A + 1e-8 for x in vals), "largest_LB": max(vals), "candidate_area": A}
    return res


def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("--refresh", action="store_true")
    ap.add_argument("--drill", action="store_true")
    ap.add_argument("--blind", metavar="CID")
    args = ap.parse_args()
    if args.refresh:
        recs = {str(n): (gate.fetch_arxiv(r["id"]) if r["kind"] == "arxiv" else gate.fetch_crossref(doi=r["doi"]))
                for n, r in REFS.items()}
        abstracts = {str(n): gate.fetch_arxiv_abstract(REFS[n]["id"]) for n in ABSTRACTS}
        SNAPSHOT.write_text(json.dumps({"retrieved": datetime.date.today().isoformat(), "records": recs,
                                        "abstracts": abstracts}, ensure_ascii=False, indent=1) + "\n", encoding="utf-8")
    snapshot = json.loads(SNAPSHOT.read_text(encoding="utf-8"))
    text, crystal, validation, manifest = load_inputs()

    def build():
        return {"text": text, "crystal": crystal, "validation": json.loads(json.dumps(validation)),
                "manifest": manifest, "snapshot": json.loads(json.dumps(snapshot)), "registry": registry.load(),
                "definitions": dict(R05_DEFINITIONS), "reviews": list(R05_REVIEWS), **CLAIMS}

    if args.blind:
        CHECKS[args.blind] = lambda w: []
        args.drill = True
    w = build()
    results = {c: fn(w) for c, fn in CHECKS.items()}
    for c, f in results.items():
        print(f"{c} {'PASS' if not f else 'FAIL'}", *(f"\n     {x}" for x in f[:6]))
    power = package_assertion_power(w)
    for key in ("w2", "w3", "n2_rows", "n4", "p1_tool"):
        print(key, w.get(key))
    for key in ("n1_rows", "n3_rows", "p1_rows"):
        for row in w.get(key) or []:
            print(key, row)
    print("package assertion power", power)
    if not args.drill:
        return
    val = validation["sanity_check"]
    row0 = val["cell_lower_bound_convergence"][0]
    chiral = lambda t: 0.5 - 0.0375 * np.cos(3 * np.asarray(t)) - 0.005 * np.sin(5 * np.asarray(t))
    cases = [
        ("K0-unmodified", "control", {}, set()),
        ("K1-eps-rewritten", "control",
         {"eps_b": lambda q, R, M: q + 4 * R * (sp.sin(sp.pi / (2 * M)) if _sym(q, R, M) else math.sin(math.pi / (2 * M)))},
         set()),
        ("K2-phi-by-arctan", "control",
         {"phi": lambda r, s: (sp.pi * r ** 2 + r * sp.sqrt(s ** 2 - r ** 2) - r ** 2 * sp.atan(sp.sqrt(s ** 2 - r ** 2) / r))
          if _sym(r, s) else (math.pi * r ** 2 + r * math.sqrt(s * s - r * r) - r * r * math.atan(math.sqrt(s * s - r * r) / r))},
         set()),
        ("K3-doc-crlf", "control", {"text": text.replace("\n", "\r\n")}, set()),
        ("K4-bound-factored", "control",
         {"cont_bound": lambda A, e: (sp.pi if _sym(A, e) else math.pi) * e * (4 * A + e)}, set()),
        ("K5-outer-factored", "control",
         {"lb_outer": lambda a, p, e: a - e * (p + (sp.pi if _sym(a, p, e) else math.pi) * e)}, set()),
        ("M1-doc-eps-edited", "mutation", {"text": text.replace(r"2R_{\mathbf b}d_M", r"R_{\mathbf b}d_M")}, {"A0"}),
        ("M2-ref-year", "mutation", {"snapshot": {**snapshot, "records": {
            **snapshot["records"], "2": {**snapshot["records"]["2"], "years": [2019]}}}}, {"A1"}),
        ("M3-mishra-bound", "mutation", {"snapshot": {**snapshot, "abstracts": {**snapshot["abstracts"], "3": {
            **snapshot["abstracts"]["3"], "abstract": snapshot["abstracts"]["3"]["abstract"].replace("0.8344", "0.8343")}}}},
         {"X1"}),
        ("M4-family-misnamed", "mutation", {"mishra_family": ("disk", "Reuleaux triangle", "regular pentagon")}, {"X1"}),
        ("M5-lambda1-display", "mutation", {"lambda1_display": "0.7853981633"}, {"W1"}),
        ("M6-gap-display", "mutation", {"gap_display": "0.0096935945"}, {"W1"}),
        ("M7-chiral-witness", "mutation", {"witnesses": {**CLAIMS["witnesses"], "Reuleaux pentagon": chiral}}, {"W3"}),
        ("M8-phi-sector-sign", "mutation",
         {"phi": lambda r, s: phi_text(r, s) + 2 * r ** 2 * (sp.acos(r / s) if _sym(r, s) else math.acos(r / s))},
         {"D4", "N2", "N3", "W2"}),
        ("M9-phi-no-protrusion", "mutation", {"phi": lambda r, s: (sp.pi if _sym(r, s) else math.pi) * r ** 2 + 0 * s},
         {"D4", "N2", "W2"}),
        ("M10-eps-one-sided", "mutation",
         {"eps_b": lambda q, R, M: q + R * 2 * (sp.sin(sp.pi / (2 * M)) if _sym(q, R, M) else math.sin(math.pi / (2 * M)))},
         {"D1", "P1"}),
        ("M11-eps-tiny", "mutation", {"eps_b": lambda q, R, M: 1e-9 + 0 * q}, {"D1", "N1", "N3", "P1"}),
        # visible to measurement only on cells of covers with long straight edges off the grid
        ("M12-outer-without-perimeter", "mutation",
         {"lb_outer": lambda a, p, e: a - (sp.pi if _sym(a, p, e) else math.pi) * e ** 2}, {"D2", "N3", "P1"}),
        ("M13-outer-is-area", "mutation", {"lb_outer": lambda a, p, e: a + 0 * e}, {"D2", "N3", "P1"}),
        ("M14-env-factor-half", "mutation", {"env_factor": sp.Rational(1, 2)}, {"D3", "N3", "P1"}),
        ("M15-bound-perimeter-pi-Abar", "mutation",
         {"cont_bound": lambda A, e: (sp.pi if _sym(A, e) else math.pi) * (A * e + e ** 2)}, {"D5", "N4"}),
        ("M16-bound-without-square", "mutation",
         {"cont_bound": lambda A, e: 4 * (sp.pi if _sym(A, e) else math.pi) * A * e}, {"D5"}),
        ("M17-cell-lower-is-upper", "mutation", {"b_lower": lambda b, q: np.maximum(0.5, b + 0 * q)}, {"N1", "P1"}),
        ("M18-quantize-downward", "mutation", {"quantize": lambda hj, q: 0.5 + q * np.floor((hj - 0.5) / q)}, {"N1", "P1"}),
        ("M19-recorded-LB-env", "mutation", {"validation": {**validation, "sanity_check": {
            **val, "cell_lower_bound_convergence": [{**row0, "LB_env": row0["LB_env"] + 1e-9}]
            + val["cell_lower_bound_convergence"][1:]}}}, {"P1"}),
        ("M20-package-sha", "mutation",
         {"validation": {**validation, "round05": {**validation["round05"], "sha256": "0" * 64}}}, {"P2"}),
        ("M21-review-dropped", "mutation", {"reviews": list(R05_REVIEWS[1:])}, {"R1"}),
    ]
    if args.blind:
        cases = [c for c in cases if args.blind in c[3]]
    d = gate.run_drill(build, CHECKS, DEPENDS, cases)
    summary = gate.summarize(d, CHECKS)
    print(gate.report_drill(summary))
    for x in d:
        if not x["ok"]:
            print("   DRILL FAIL", x)
    if args.blind:
        n = sum(1 for x in d if args.blind in x["expected"] and args.blind not in x["tripped"])
        print(f"blind {args.blind}: drill reported it blind in {n} cases (nothing written)")
        return
    log = {"run": "RUN-008", "subject": "Round 05: the saturation gate", "document": DOC,
           "snapshot_retrieved": snapshot["retrieved"],
           "checks": {c: {"pass": not f, "failures": f} for c, f in results.items()},
           "w2": w.get("w2"), "w3": w.get("w3"), "n1_rows": w.get("n1_rows"), "n2_rows": w.get("n2_rows"),
           "n3_rows": w.get("n3_rows"), "n4": w.get("n4"), "p1_rows": w.get("p1_rows"), "p1_tool": w.get("p1_tool"),
           "package_assertion_power": power, "drill": summary}
    LOG.write_text(json.dumps(log, ensure_ascii=False, indent=1, default=float) + "\n", encoding="utf-8")
    if not results["R1"]:
        registry.record("R05", "RUN-008", R05_DEFINITIONS, R05_REVIEWS)
        print(f"registry: {len(R05_DEFINITIONS)} definitions, {len(R05_REVIEWS)} reviews recorded for R05")


if __name__ == "__main__":
    main()