"""RUN-053 - LESR Paper 03: support-function variational calculus.

LESR Paper 03 (2026-09-20) builds the series' calculus layer. Area is the mixed-area form B(h, h), exactly
quadratic along affine support paths h + t phi; its first variation is int phi d mu_h and its second
int (phi^2 - phi'^2), independent of h; the constant direction gives the perimeter (Steiner); translations and
rotations give 0. Convexity along a path is mu_h + t L phi >= 0, with a step cap bounded below by 1 / ||q|| when the
negative part of L phi has a bounded density q against mu_h, and 0 when it does not. A witness W held at a placement
(psi, t) is the constraint h(theta) - h_W(theta - psi) - t.u >= 0; with container direction phi and placement
velocity (omega, v) its first variation at an active contact is phi(theta_j) + omega h_W'(theta_j - psi) - v.u_j.
First-order feasibility and a descent test over a finite mode library are then finite LPs (sec. 29-39), and sec. 61
turns an LP optimum >= 0 into a "finite-mode first-order closure certificate". It ships a schema, a sanity-check
record, a compiler scaffold, a state record and VALIDATION. It is evidence, not a canonical document. Here the
formulas are checked in sympy and on bodies, the step caps and the contact variation are computed, and the LP's
first-order model is tested on the Lebesgue kernel itself, where the container is the hull of the placed witnesses.

Checks
  A0 anchors      the formulas below occur in the document as written
  W1 formulas     sec. 5-13, 15, 55-56 in sympy; the sanity-check record's numbers
  W2 bodies       sec. 9-11 on bodies: the Steiner formula and the perimeter on the kernel's hull and two Reuleaux
                  polygons; the exact quadratic along an admissible path of the Reuleaux triangle, against its drawn
                  boundary
  W3 caps         sec. 17-21: the circle's cap; the triangle's one admissible mode, by sec. 19's density; the square,
                  where the negative part of L phi falls where mu_h is empty
  W4 contacts     sec. 26-29: the first variation of a fixed-placement contact against finite differences
  W5 kernel       sec. 27-39, 61 on the kernel: the first-order area change of witness moves - measured, by edge
                  averages, and by the affine model the LP uses - at a generic placement, at the mirror-symmetric
                  polish, and at q-dagger, where the LP is solved
  W6 records      the schema and the state against the sections that carry them
  W7 scaffold     the compiler scaffold, read after this rebuild and never run
  P1 package      the VALIDATION record against the bytes and the archive

Usage: python code/run053_lesr_paper03.py [--drill] [--blind CID]
"""
import argparse
import copy
import hashlib
import inspect
import json
import math
import re

import numpy as np
import sympy as sp
from scipy.optimize import linprog

import gate
import run044_uesfcm_cycles_01_03 as r44
import run046_uesfcm_cycles_07_09 as r46
from corpus import SOURCES, TREE, source_bytes, walk_sources
from run009_round06_witness_exchange import hull_pieces, placed_reuleaux
from run052_lesr_paper02 import QDAG, P5, integral_2pi, segs_of

FOLDER = "main.zip!/LESR_Skew_Field_Variational/"
PKG = FOLDER + "LESR_Paper03_Package_v0.1.zip!/"
DOC = PKG + "LESR_Paper03_Support_Function_Variational_Calculus_v0.1.md"
SCHEMA = PKG + "LESR_Paper03_Variational_Calculus_Schema_v0.1.json"
CHECKS_JSON = PKG + "LESR_Paper03_Sanity_Checks.json"
SCAFFOLD = PKG + "LESR_Paper03_Variational_Compiler_Scaffold.py"
STATE = PKG + "LESR_State_Paper03_v0.1.json"
VALIDATION = PKG + "LESR_Paper03_VALIDATION.json"
PAPER04 = FOLDER + "LESR_Paper04_Package_v0.1.zip!/LESR_Paper04_Skew_Constrained_Shape_Flow_and_Closure_Graph_v0.1.md"
LOG = TREE / "data" / "gate-logs" / "run053-lesr-paper03.json"
MANIFEST = TREE / "data" / "SOURCE_MANIFEST.v1.json"
CACHE = TREE / "data" / "cache" / "run053-searches.json"     # data/cache is not committed
# RUN-052's mirror-symmetric polish of q-dagger (y3 = y5 = 0, phi5 = 2 pi / 5, x3 and x5 by golden sections)
ZPOL = [-0.012158832628610313, 0.0, P5, 0.022996606917054876, 0.0]
J = np.array([[0.0, -1.0], [1.0, 0.0]])


def load_inputs():
    blobs = {n: source_bytes(n) for n in SOURCES}
    want = (DOC, SCHEMA, CHECKS_JSON, SCAFFOLD, STATE, VALIDATION)
    got = {leaf.id: d for leaf, d in walk_sources(blobs, keep=lambda l: l.id in want + (PAPER04,)) if d is not None}
    rows = {row[0]: row for row in json.loads(MANIFEST.read_text(encoding="utf-8"))["leaves"]}
    return {"text": got[DOC].decode("utf-8"), "schema": json.loads(got[SCHEMA]), "sanity": json.loads(got[CHECKS_JSON]),
            "scaffold": got[SCAFFOLD].decode("utf-8"), "state": json.loads(got[STATE]), "validation": json.loads(got[VALIDATION]),
            "bytes": {k: got[k] for k in want}, "manifest_rows": {k: rows[k] for k in want}, "paper04": got[PAPER04].decode("utf-8")}


# ---------------------------------------------------------------- the formulas

def formulas():
    """sec. 5-13, 15, 55-56 as identities over trigonometric polynomials with free coefficients; the sanity record's
    numbers recomputed"""
    t, s, R, a, l = sp.symbols("theta t R a ell", real=True)
    c = sp.symbols("c0:5", real=True)
    d = sp.symbols("d1:5", real=True)
    f = sp.symbols("f0:5", real=True)
    g = sp.symbols("g1:5", real=True)
    h = c[0] + sum(c[m] * sp.cos(m * t) + d[m - 1] * sp.sin(m * t) for m in range(1, 5))
    p = f[0] + sum(f[m] * sp.cos(m * t) + g[m - 1] * sp.sin(m * t) for m in range(1, 5))
    I = lambda e: integral_2pi(e, t)
    D = lambda e, n=1: sp.diff(e, t, n)
    B = lambda x, y: I(x * y - D(x) * D(y)) / 2
    A = lambda x: B(x, x)
    z = lambda e: sp.simplify(sp.expand(e)) == 0
    out = {"symmetric": z(B(h, p) - B(p, h)), "by_parts": z(B(h, p) - I(h * (p + D(p, 2))) / 2),
           "quadratic": z(A(h + s * p) - A(h) - 2 * s * B(h, p) - s ** 2 * B(p, p)),
           "first": z(2 * B(h, p) - I(p * (h + D(h, 2)))), "second": z(2 * B(p, p) - I(p ** 2 - D(p) ** 2)),
           "circle": z(A(R + s) - (sp.pi * R ** 2 + 2 * sp.pi * R * s + sp.pi * s ** 2)), "perimeter": z(I(h + D(h, 2)) - I(h)),
           "translation": [z(2 * B(h, sp.cos(t))), z(2 * B(sp.cos(t) + sp.sin(t), sp.cos(t) + sp.sin(t)))],
           "rotation": z(2 * B(h, -D(h))), "mixed": z(I(p * (h + D(h, 2))) - I(h * (p + D(p, 2))))}
    # sec. 13: the affine path along the rotation tangent is not the rotation - it lowers the area at second order
    h3 = 1 + a * sp.cos(3 * t)
    out["rotation_affine"] = str(sp.factor(sp.expand(A(h3 - s * D(h3)) - A(h3))))
    out["rotation_exact"] = z(A(h3.subs(t, t - s)) - A(h3))
    # sec. 55-56: the L2 gradient is h + h''; the H1 gradient g solves (1 - l^2 d^2) g = h + h'', mode by mode
    grad = c[0] + sum((1 - m * m) * (c[m] * sp.cos(m * t) + d[m - 1] * sp.sin(m * t)) / (1 + l ** 2 * m * m) for m in range(1, 5))
    out["l2_gradient"] = z(I((h + D(h, 2)) * p) - 2 * B(h, p))
    out["h1_gradient"] = [z(grad - l ** 2 * D(grad, 2) - (h + D(h, 2))), z(I(grad * p + l ** 2 * D(grad) * D(p)) - 2 * B(h, p))]
    out = {k: ([bool(x) for x in v] if isinstance(v, list) else (v if isinstance(v, str) else bool(v))) for k, v in out.items()}
    # the sanity record, recomputed: a circle of radius 2 with phi = 1; the translation mode; cos 2 theta on the unit circle
    Rv = 2
    out["sanity"] = {"circle": [math.pi * Rv ** 2, 2 * math.pi * Rv, 2 * math.pi], "m1": [1 - 1, 0.0],
                     "m2": [1 - 4, math.pi * (1 - 4), 1 / 3]}
    return out


# ---------------------------------------------------------------- bodies

def hull_ap(pieces):
    """the verified hull's area and perimeter, with RUN-046's grid retries"""
    for n in (16384, 16381, 12289, 8191, 4096):
        try:
            return hull_pieces(pieces, n)
        except RuntimeError:
            continue
    raise RuntimeError("no grid")


def body_checks():
    """sec. 10-11 on bodies: A(K + sB) = A + P s + pi s^2 (every piece's radius raised by s) and delta A[h; 1] =
    mu_h(S^1) = P, from the verified segments; sec. 9 along RT - t cos 3 theta, t in [0, 1/8] (admissible: RUN-052),
    its area from the boundary x = h u + h' u' drawn at 2^16 normals against A0 + 2t - 4 pi t^2, and its width"""
    bodies = {"kernel": r44.base_pieces(QDAG), "reuleaux3": placed_reuleaux(3, 0.0, (0.0, 0.0)), "reuleaux5": placed_reuleaux(5, 0.0, (0.0, 0.0))}
    out = {"steiner": 0.0, "perimeter": 0.0}
    for pieces in bodies.values():
        A0, P0 = hull_ap(pieces)
        for sv in (0.01, 0.1, 0.5):
            As, _ = hull_ap([(cc, rr + sv, rg) for cc, rr, rg in pieces])
            out["steiner"] = max(out["steiner"], abs(As - A0 - P0 * sv - math.pi * sv * sv))
        C, rho, segs = segs_of(pieces)
        dh = lambda k, x: -C[k, 0] * math.sin(x) + C[k, 1] * math.cos(x)
        mass = sum(rho[k] * (b - a) for a, b, k in segs) + sum(dh(segs[(i + 1) % len(segs)][2], b) - dh(k, b) for i, (a, b, k) in enumerate(segs))
        out["perimeter"] = max(out["perimeter"], abs(mass - P0))
    # the Reuleaux triangle and the one pure mode it admits
    N = 1 << 16
    th = 2 * np.pi * np.arange(N) / N
    R3 = 1 / math.sqrt(3)
    V = [np.array([R3 * math.cos(2 * math.pi * k / 3), R3 * math.sin(2 * math.pi * k / 3)]) for k in range(3)]
    U, Up = np.stack([np.cos(th), np.sin(th)], 1), np.stack([-np.sin(th), np.cos(th)], 1)
    vals = np.stack([U @ v for v in V] + [U @ V[k] + 1.0 for k in range(3)])
    ders = np.stack([Up @ v for v in V] + [Up @ V[k] for k in range(3)])
    for k in range(3):     # an arc counts only on the normals within pi/6 of the direction opposite its corner
        off = np.abs((th - (2 * math.pi * k / 3 + math.pi) + math.pi) % (2 * math.pi) - math.pi) > math.pi / 6
        vals[3 + k, off] = -np.inf
    own = np.argmax(vals, 0)
    h0, h0p = vals[own, np.arange(N)], ders[own, np.arange(N)]
    A0 = (math.pi - math.sqrt(3)) / 2
    path = []
    for tv in (1 / 16, 1 / 8):
        hh, hp = h0 - tv * np.cos(3 * th), h0p + 3 * tv * np.sin(3 * th)
        X = hh[:, None] * U + hp[:, None] * Up
        poly = 0.5 * float(np.sum(X[:, 0] * np.roll(X[:, 1], -1) - np.roll(X[:, 0], -1) * X[:, 1]))
        path.append(abs(poly - (A0 + 2 * tv - 4 * math.pi * tv * tv)))
        width = float(np.max(np.abs(hh + np.roll(hh, -N // 2) - 1.0)))
    out["path"] = max(path)
    out["gain"] = [A0 + 2 * tv - 4 * math.pi * tv * tv - A0 for tv in (1 / 16, 1 / 8)]
    out["width"] = width
    return out


# ---------------------------------------------------------------- step caps

def cap_checks():
    """sec. 17-21. The unit circle and cos 2 theta: the density cap of sec. 20. The Reuleaux triangle and
    phi = -cos 3 theta: L phi = 8 cos 3 theta, whose negative part lives on the arcs (density 1), q = 8 |cos 3 theta|
    there, so sec. 19 gives 1/8 - which is also the largest admissible step found directly (RUN-052). The unit square
    (atoms at its four edge normals) and phi = +-cos 2 theta: L phi = -+3 cos 2 theta is negative on open arcs where
    mu_h is empty, so the cap is 0 for both signs (sec. 18)"""
    N = 1 << 16
    th = 2 * np.pi * np.arange(N) / N
    circle = float(np.min(1.0 / np.maximum(3 * np.cos(2 * th), 1e-300)[3 * np.cos(2 * th) > 0]))
    arcs = np.zeros(N, dtype=bool)
    for k in range(3):
        arcs |= np.abs((th - (2 * math.pi * k / 3 + math.pi) + math.pi) % (2 * math.pi) - math.pi) <= math.pi / 6
    Lphi = 8 * np.cos(3 * th)
    q = np.where(arcs, np.maximum(-Lphi, 0.0), 0.0)
    corner_negative = bool(np.any((Lphi < -1e-12) & ~arcs))
    square = []
    for sgn in (1, -1):
        Lp = -3 * sgn * np.cos(2 * th)
        atoms = np.isclose(np.mod(th, math.pi / 2), 0.0) | np.isclose(np.mod(th, math.pi / 2), math.pi / 2)
        square.append(bool(np.any((Lp < -1e-9) & ~atoms)))
    return {"circle": circle, "reuleaux": 1.0 / float(np.max(q)), "corners_clear": not corner_negative, "square_blocked": square}


# ---------------------------------------------------------------- a fixed-placement contact

def contact_check(eps=1e-7):
    """sec. 26-28: the unit disc holds a triangle placed at psi = 0, t = (0.2, 0.1), two of its corners on the circle
    (normals 0.7 and 3.0) and one inside. For several container directions phi and placement velocities (omega, v),
    the least of C_eps near each contact, over eps, against phi(theta_j) + omega h_W'(theta_j - psi) - v.u_j"""
    tv = np.array([0.2, 0.1])
    world = [np.array([math.cos(0.7), math.sin(0.7)]), np.array([math.cos(3.0), math.sin(3.0)]), 0.6 * np.array([math.cos(4.4), math.sin(4.4)])]
    q = [w - tv for w in world]                      # the triangle in its own frame
    hW = lambda x: max(float(p @ np.array([math.cos(x), math.sin(x)])) for p in q)
    hWp = lambda x: float(q[int(np.argmax([p @ np.array([math.cos(x), math.sin(x)]) for p in q]))] @ np.array([-math.sin(x), math.cos(x)]))
    rng = np.random.default_rng(53)
    worst, formula_uses_omega = 0.0, 0.0
    for _ in range(6):
        a2, b1, om = rng.normal(size=3)
        vv = rng.normal(size=2)
        phi = lambda x: a2 * math.cos(2 * x) + b1 * math.sin(x)
        for tj in (0.7, 3.0):
            pred = phi(tj) + om * hWp(tj) - float(vv @ np.array([math.cos(tj), math.sin(tj)]))
            grid = tj + np.linspace(-1e-3, 1e-3, 4001)
            Ce = [1.0 + eps * phi(x) - hW(x - eps * om) - float((tv + eps * vv) @ np.array([math.cos(x), math.sin(x)])) for x in grid]
            worst = max(worst, abs(min(Ce) / eps - pred))
            formula_uses_omega = max(formula_uses_omega, abs(om * hWp(tj)))
    return {"worst": worst, "omega_term": formula_uses_omega}


# ---------------------------------------------------------------- the kernel: witness moves at first order

def structure(v, tie=1e-9, act=1e-10):
    """the hull at v: segments (tie-level ones dropped), and at each hand-over its normal, its atom (the edge's
    length), the two owners either side and every piece that touches the hull there"""
    P = r44.base_pieces(v)
    C, rho, segs = segs_of(P)
    ranges = [p[2] for p in P]
    S = [(a, b, k) for a, b, k in segs if b - a > tie]
    h = lambda k, x: float(C[k] @ np.array([math.cos(x), math.sin(x)]) + rho[k])
    dh = lambda k, x: float(C[k] @ np.array([-math.sin(x), math.cos(x)]))
    inr = lambda k, x: ranges[k] is None or abs((x - ranges[k][0] + math.pi) % (2 * math.pi) - math.pi) <= ranges[k][1] + 1e-12
    atoms = []
    for i, (a, b, k) in enumerate(S):
        k2 = S[(i + 1) % len(S)][2]
        H = max(h(j, b) for j in range(len(P)) if inr(j, b))
        # an edge's length is >= 0; a tangent hand-over gives 0, which rounding can make -1e-17
        atoms.append((b, max(dh(k2, b) - dh(k, b), 0.0), k, k2, [j for j in range(len(P)) if inr(j, b) and h(j, b) >= H - act]))
    return {"C": C, "rho": rho, "S": S, "atoms": atoms}


def velocities(C, v, xi):
    """each piece's centre velocity for xi = (v3x, v3y, w5, v5x, v5y): the disc is fixed (the translation gauge),
    B3 translates (its turn is the rotation gauge), B5 turns about its centre and translates"""
    c5 = np.array([v[3], v[4]])
    return [np.zeros(2)] + [np.array(xi[:2], dtype=float)] * 6 + [np.array(xi[3:5], dtype=float) + xi[2] * (J @ (C[k] - c5)) for k in range(7, 17)]


def first_order(st, v, xi):
    """the first-order area change of the witness move xi, two ways. 'average': each arc moves with its piece
    (int rho c'.u), each edge by its length times the mean of its two ends' normal speeds - the trapezoid, since the
    container is the hull and its edges turn. 'affine': the least int phi d mu_h over container directions phi
    (continuous, with L phi a measure) that keep every witness inside - on each arc the piece's speed, at each edge
    normal the largest speed of the pieces touching there: the edge moves parallel to itself"""
    cd = velocities(st["C"], v, xi)
    sp_ = lambda k, x: float(cd[k] @ np.array([math.cos(x), math.sin(x)]))
    arcs = sum(st["rho"][k] * (cd[k][0] * (math.sin(b) - math.sin(a)) - cd[k][1] * (math.cos(b) - math.cos(a))) for a, b, k in st["S"] if st["rho"][k] > 0)
    avg = arcs + sum(m * (sp_(k, x) + sp_(k2, x)) / 2 for x, m, k, k2, _ in st["atoms"])
    aff = arcs + sum(m * max(sp_(j, x) for j in act) for x, m, _, _, act in st["atoms"])
    return avg, aff


def lp_descent(st, v):
    """sec. 39 with the container free and the witnesses re-placed: min over |xi| <= 1 of the affine first-order
    change - variables xi and one epigraph value per edge normal"""
    n = len(st["atoms"])
    unit = [first_order_parts(st, v, e) for e in np.eye(5)]
    cobj = np.concatenate([[u[0] for u in unit], [m for _, m, _, _, _ in st["atoms"]]])
    rows = []
    for i, (x, m, _, _, act) in enumerate(st["atoms"]):
        for j in act:
            r = np.zeros(5 + n)
            r[:5] = [float(velocities(st["C"], v, e)[j] @ np.array([math.cos(x), math.sin(x)])) for e in np.eye(5)]
            r[5 + i] = -1.0
            rows.append(r)
    res = linprog(cobj, A_ub=np.array(rows), b_ub=np.zeros(len(rows)), bounds=[(-1, 1)] * 5 + [(None, None)] * n, method="highs")
    return float(res.fun), [float(x) for x in res.x[:5]]


def first_order_parts(st, v, xi):
    """the arcs' share of the first-order change (linear in xi), for the LP's objective"""
    cd = velocities(st["C"], v, xi)
    return (sum(st["rho"][k] * (cd[k][0] * (math.sin(b) - math.sin(a)) - cd[k][1] * (math.cos(b) - math.cos(a))) for a, b, k in st["S"] if st["rho"][k] > 0),)


def kernel_scan(steps=(1e-7, 1e-8)):
    """the family area's one-sided rates (A(v + t xi) - A(v)) / t against the two first-order models: at a generic
    placement near q-dagger (four random directions, the gradient by central differences against the edge-average
    one, the steepest descent direction, and the LP); at the polish Z along y3, phi5, y5 (both ways); at q-dagger
    toward Z; and the LP at q-dagger"""
    A = r46.area
    rate = lambda v, xi, t: (A([v[i] + t * xi[i] for i in range(5)]) - A(v)) / t
    rng = np.random.default_rng(53)
    gen = [QDAG[0] + 0.002, QDAG[1] + 0.001, QDAG[2] + 0.05, QDAG[3] - 0.001, QDAG[4] + 0.002]
    st = structure(gen)
    generic = []
    for _ in range(4):
        xi = rng.normal(size=5)
        avg, aff = first_order(st, gen, xi)
        generic.append({"measured": [rate(gen, xi, t) for t in steps], "average": avg, "affine": aff})
    E = np.eye(5)
    grad = [(A([gen[j] + 1e-7 * E[i][j] for j in range(5)]) - A([gen[j] - 1e-7 * E[i][j] for j in range(5)])) / 2e-7 for i in range(5)]
    grad_avg = [first_order(st, gen, E[i])[0] for i in range(5)]
    dn = -np.array(grad) / np.linalg.norm(grad)
    steepest = {"measured": rate(gen, dn, 1e-7), "average": first_order(st, gen, dn)[0], "affine": first_order(st, gen, dn)[1]}
    lp_gen = lp_descent(st, gen)
    stz = structure(ZPOL)
    atz = {}
    for name, i in (("y3", 1), ("phi5", 2), ("y5", 4)):
        atz[name] = []
        for sgn in (1, -1):
            xi = [0.0] * 5
            xi[i] = float(sgn)
            atz[name].append({"measured": rate(ZPOL, xi, 1e-7), "affine": first_order(stz, ZPOL, xi)[1]})
    stq = structure(QDAG)
    dvec = np.array(ZPOL) - np.array(QDAG)
    dvec /= np.linalg.norm(dvec)
    toward = {"measured": [rate(QDAG, dvec, t) for t in steps], "affine": first_order(stq, QDAG, dvec)[1],
              "ties": sum(1 for a in stq["atoms"] if len(a[4]) > 2)}
    lp, arg = lp_descent(stq, QDAG)
    return {"generic": generic, "ties_generic": sum(1 for a in st["atoms"] if len(a[4]) > 2), "gradient": grad, "gradient_average": grad_avg,
            "steepest": steepest, "lp_generic": lp_gen[0], "lp_generic_argmin": lp_gen[1], "z": atz, "toward": toward, "lp": lp, "lp_argmin": arg}


# ---------------------------------------------------------------- the cache

def _key(name, fn, extra=""):
    import run008_round05_saturation_gate as r8
    import run009_round06_witness_exchange as r9
    import run052_lesr_paper02 as r52
    digest = hashlib.sha256(inspect.getsource(fn).encode())
    for helper in HELPERS[name]:
        digest.update(inspect.getsource(helper).encode())
    for path in (r8.__file__, r9.__file__, r52.__file__):
        digest.update(open(path, "rb").read())
    digest.update(extra.encode())
    return digest.hexdigest()


HELPERS = {"formulas": (), "bodies": (hull_ap, r44.base_pieces), "caps": (), "contact": (),
           "kernel": (structure, velocities, first_order, lp_descent, first_order_parts, r46.area, r46.hull, r44.base_pieces)}


def searches():
    cached = json.loads(CACHE.read_text(encoding="utf-8")) if CACHE.exists() else {}
    out = {}
    for name, fn in (("formulas", formulas), ("bodies", body_checks), ("caps", cap_checks), ("contact", contact_check), ("kernel", kernel_scan)):
        key = _key(name, fn)
        if cached.get(name, {}).get("key") != key:
            cached[name] = {"key": key, "value": json.loads(json.dumps(fn()))}
            CACHE.parent.mkdir(parents=True, exist_ok=True)
            CACHE.write_text(json.dumps(cached, indent=1) + "\n", encoding="utf-8")
        out[name] = cached[name]["value"]
    return out


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

ANCHORS = {
    "A0": [r"B(f,g)=\frac12\int_0^{2\pi}\left(fg-f'g'\right)d\theta.", r"A[h+t\phi]=A[h]+2tB(h,\phi)+t^2B(\phi,\phi).", r"\deltaA[h;\phi]=2B(h,\phi).",
           r"\delta^2A[\phi]=2B(\phi,\phi).", r"A(K+tB)=A(K)+P(K)t+\pit^2", r"\deltaA[h;1]=\mu_h(S^1)=P(K).", r"\deltaA[h;-h']=0.",
           r"\mu_{h+t\phi}=\mu_h+t\nu_\phi.", r"\deltaA[h;\phi]=\inth\,d\nu_\phi", r"\varepsilon_{\mathrm{conv}}^+\ge\frac1{\|q\|_\infty}.",
           r"\varepsilon_{\mathrm{conv}}^+=0.", r"\deltaC_j=\phi(\theta_j)+\omegah_W'(\theta_j-\psi)-v\cdotu_{\theta_j}.",
           r"\min\sum_kc_k\deltaA[\phi_k]", r"\nabla_{L^2}A=h+h''.", r"(1-\ell^2\partial_\theta^2)\nabla_{H^1}A=h+h''.",
           r"\text{finite-modestationarity}\neq\text{continuumstationarity}."],
}
STATE_MAP = {"mixed-area bilinear form": (5, "BilinearMixed-AreaForm"), "exact quadratic support path area": (6, "ExactQuadraticPathFormula"),
             "first area variation": (7, "FirstVariation"), "second area variation": (8, "SecondVariation"),
             "smooth/non-smooth split": (14, "Non-SmoothMeasureLane"), "measure-valued convexity tangent cone": (16, "ConvexityTangentCone"),
             "convexity step cap": (17, "ConvexityStepCap"), "translation/rotation gauge sanity checks": (12, "TranslationGaugeSanityCheck"),
             "fixed-placement witness constraint variation": (28, "FirstVariationofFixed-ContactConstraint"),
             "infinitesimal placement compensation cone": (29, "InfinitesimalPlacement-CompensationCone"),
             "finite witness first-order feasible cone": (32, "First-OrderFeasibleShapeCone"), "finite mode LP descent test": (39, "FiniteLPDescentTest"),
             "shape velocity/integration": (42, "ShapeVelocity"), "hybrid continuous/discrete shape dynamics": (47, "HybridShapeDynamics"),
             "metric-free vs metric-dependent calculus boundary": (53, "Metric-FreeCalculus到此為止")}
DISCLAIMERS = {"unique shape gradient without choosing a metric": "不同metric給不同gradient",
               "smooth calculus across contact topology transitions": "停止smoothflow並重新分類",
               "full derivative formula for placement-optimized continuum universal-cover constraint": "Paper03只建立它的finitestable-contactlocalversion",
               "finite-mode stationarity implies continuum stationarity": r"\text{finite-modestationarity}\neq\text{continuumstationarity}",
               "negative unconstrained second area variation implies feasible descent": "就宣稱feasibledescent"}
CLAIMS = {
    "formulas_claim": {"symmetric": True, "by_parts": True, "quadratic": True, "first": True, "second": True, "circle": True, "perimeter": True,
                       "translation": [True, True], "rotation": True, "mixed": True, "rotation_affine": "-36*pi*a**2*t**2", "rotation_exact": True,
                       "l2_gradient": True, "h1_gradient": [True, True], "sanity": [True] * 8},
    "bodies_claim": {"steiner": True, "perimeter": True, "path": True, "width": True, "gain": [0.075912615, 0.053650459]},
    "caps_claim": {"circle": 0.333333333333, "reuleaux": 0.125, "corners_clear": True, "square_blocked": [True, True]},
    "contact_claim": {"agrees": True, "omega_term": True},
    "kernel_claim": {"average_matches": True, "affine_gaps": [0.66, 0.94, 0.7, 0.85], "ties_generic": 0, "gradient": [0.39, True],
                     "steepest": [-0.39, 0.107], "lp_generic": True,
                     "z": {"y3": [[0.1327, 0.732], [0.1327, 0.732]], "phi5": [[0.0442, 0.1274], [0.0442, 0.1274]], "y5": [[0.2733, 0.795], [0.2733, 0.795]]},
                     "toward": [-6.3e-06, 0.02], "lp": True},
    "records_claim": {"formulas": [True, True, True, True], "cones": True, "unmapped": [], "completed": 15, "disclaimed": [True] * 5, "next": [True, True]},
    "scaffold_facts": {"lines": 61, "defs": ["contact_linear_row", "exact_area_along_direction", "exact_line_minimizer"],
                       "classes": ["FourierMode", "ActiveContact"], "lp_solver": False, "row_signs": True, "m0_variation": True},
    "package_claim": {"digests": [True, True, True, True, True], "utf8": True, "delimiters": ["PASS", []], "manifest": True},
}


def sig(x, n=2):
    return float(f"{x:.{n - 1}e}")


def a0_anchors(w):
    text = gate.squash(w["text"])
    return [f"anchor not found: {a[:60]}" for a in ANCHORS["A0"] if a not in text]


def w1_formulas(w):
    """sec. 5-13, 15, 55-56; and the sanity record's values against the recomputation"""
    s = w["search"]["formulas"]
    rec = w["sanity"]
    same = lambda x, y: abs(x - y) <= 1e-12 * max(1.0, abs(y))
    sanity = [same(rec["circle_R2_constant_variation"]["expected_area"], s["sanity"]["circle"][0]),
              same(rec["circle_R2_constant_variation"]["expected_first_variation_phi1"], s["sanity"]["circle"][1]),
              same(rec["circle_R2_constant_variation"]["expected_second_variation_phi1"], s["sanity"]["circle"][2]),
              rec["translation_mode_m1"]["convexity_eigenvalue"] == s["sanity"]["m1"][0], same(rec["translation_mode_m1"]["second_area_variation_for_cos"], s["sanity"]["m1"][1]),
              rec["mode_m2_unit_cos"]["convexity_eigenvalue"] == s["sanity"]["m2"][0], same(rec["mode_m2_unit_cos"]["second_area_variation"], s["sanity"]["m2"][1]),
              same(rec["mode_m2_unit_cos"]["circle_R1_convexity_abs_epsilon_cap"], s["sanity"]["m2"][2])]
    got = {k: s[k] for k in ("symmetric", "by_parts", "quadratic", "first", "second", "circle", "perimeter", "translation", "rotation", "mixed",
                             "rotation_affine", "rotation_exact", "l2_gradient", "h1_gradient")}
    got["sanity"] = sanity
    w["w1"] = got
    return [] if got == w["formulas_claim"] else [f"the formulas: {got}"]


def w2_bodies(w):
    s = w["search"]["bodies"]
    got = {"steiner": s["steiner"] < 1e-13, "perimeter": s["perimeter"] < 1e-13, "path": s["path"] < 1e-8, "width": s["width"] < 1e-12,
           "gain": [round(x, 9) for x in s["gain"]]}
    w["w2"] = got
    return [] if got == w["bodies_claim"] else [f"the bodies: {got}"]


def w3_caps(w):
    s = w["search"]["caps"]
    got = {"circle": round(s["circle"], 12), "reuleaux": round(s["reuleaux"], 12), "corners_clear": s["corners_clear"], "square_blocked": s["square_blocked"]}
    w["w3"] = got
    return [] if got == w["caps_claim"] else [f"the caps: {got}"]


def w4_contacts(w):
    s = w["search"]["contact"]
    got = {"agrees": s["worst"] < 1e-5, "omega_term": s["omega_term"] > 0.05}
    w["w4"] = got
    return [] if got == w["contact_claim"] else [f"the contacts: {got}"]


def w5_kernel(w):
    """the kernel: the edge-average rule is the measured first-order change; the affine model is larger; at a
    generic placement the area has a nonzero gradient and at q-dagger it falls toward Z, while the LP over affine
    moves finds nothing at either"""
    s = w["search"]["kernel"]
    g = s["generic"]
    zero = lambda v, arg: v == 0.0 and all(x == 0.0 for x in arg)
    got = {"average_matches": all(abs(x["measured"][0] - x["average"]) < 1e-5 for x in g),
           "affine_gaps": [sig(x["affine"] - x["measured"][0]) for x in g], "ties_generic": s["ties_generic"],
           "gradient": [sig(float(np.linalg.norm(s["gradient"])), 3), max(abs(a - b) for a, b in zip(s["gradient"], s["gradient_average"])) < 1e-8],
           "steepest": [sig(s["steepest"]["measured"], 3), sig(s["steepest"]["affine"], 3)], "lp_generic": zero(s["lp_generic"], s["lp_generic_argmin"]),
           "z": {k: [[round(r["measured"], 4), round(r["affine"], 4)] for r in v] for k, v in s["z"].items()},
           "toward": [sig(s["toward"]["measured"][0]), sig(s["toward"]["affine"])], "lp": zero(s["lp"], s["lp_argmin"])}
    w["w5"] = got
    return [] if got == w["kernel_claim"] else [f"the kernel: {got}"]


def w6_records(w):
    sc, st, text = w["schema"], w["state"], gate.squash(w["text"])
    sq = gate.squash
    got = {"formulas": [sq(sc["first_variation"]) == sq("delta A[h;phi] = integral phi dmu_h"), sq(sc["second_variation"]) == sq("delta2 A[phi] = integral (phi^2 - phi'^2) dtheta"),
                        sq(sc["exact_linear_path"]) == sq("A[h+t phi] = A[h] + t deltaA + t^2/2 delta2A"),
                        sq(sc["fixed_witness_contact_variation"]) == sq("phi(theta_j) + omega h_W'(theta_j-psi) - v dot u_theta_j >= 0")],
           "cones": sorted(sc["cones"]) == ["convexity", "corridor", "finite_witness", "full_local_feasible", "topology"],
           "unmapped": [x for x in st["completed"] if x not in STATE_MAP or f"#{STATE_MAP[x][0]}.{STATE_MAP[x][1]}" not in text],
           "completed": len(st["completed"]), "disclaimed": [x in DISCLAIMERS and sq(DISCLAIMERS[x]) in text for x in st["not_claimed"]],
           "next": [sc["next_paper"] in w["paper04"][:300], st["next"] == "Paper 04 — " + sc["next_paper"]]}
    w["w6"] = got
    return [] if got == w["records_claim"] else [f"the records: {got}"]


def w7_scaffold(w):
    """the scaffold, read: whether it solves any LP (sec. 39); whether its contact row has sec. 29's signs; its
    constant mode's second variation"""
    code = w["scaffold"]
    got = {"lines": code.count("\n"), "defs": re.findall(r"^def (\w+)", code, re.M), "classes": re.findall(r"^class (\w+)", code, re.M),
           "lp_solver": bool(re.search(r"linprog|simplex|highs|cvxpy|pulp", code, re.I)),
           "row_signs": bool(re.search(r"return \(\s*contact\.witness_support_derivative,\s*-math\.cos\(contact\.theta\),\s*-math\.sin\(contact\.theta\),\s*\)", code)),
           "m0_variation": bool(re.search(r"if self\.m == 0:.*?return 2\.0 \* math\.pi \* self\.a \* self\.a", code, re.S))}
    w["w7"] = got
    return [] if got == w["scaffold_facts"] else [f"the scaffold: {got}"]


def p1_package(w):
    v, B = w["validation"], w["bytes"]
    sha = lambda b: hashlib.sha256(b).hexdigest()
    body = B[DOC].decode("utf-8")
    got = {"digests": [v["paper_sha256"] == sha(B[DOC]), v["schema_sha256"] == sha(B[SCHEMA]), v["scaffold_sha256"] == sha(B[SCAFFOLD]),
                       v["checks_sha256"] == sha(B[CHECKS_JSON]), v["state_sha256"] == sha(B[STATE])],
           "utf8": v["utf8"] is True, "delimiters": [v["math_delimiters"], sorted(set(re.findall(r"\\\(|\\\)|\\\[|\\\]", body)))],
           "manifest": all(w["manifest_rows"][k][2] == sha(b) for k, b in B.items())}
    w["p1"] = got
    return [] if got == w["package_claim"] else [f"the package: {got}"]


CHECKS = {"A0": a0_anchors, "W1": w1_formulas, "W2": w2_bodies, "W3": w3_caps, "W4": w4_contacts, "W5": w5_kernel, "W6": w6_records,
          "W7": w7_scaffold, "P1": p1_package}
DEPENDS = {"A0": {"text"}, "W1": {"search", "sanity", "formulas_claim"}, "W2": {"search", "bodies_claim"}, "W3": {"search", "caps_claim"},
           "W4": {"search", "contact_claim"}, "W5": {"search", "kernel_claim"}, "W6": {"schema", "state", "text", "paper04", "records_claim"},
           "W7": {"scaffold", "scaffold_facts"}, "P1": {"validation", "bytes", "manifest_rows", "package_claim"}}


def drill_cases(base):
    def edit(obj, fn):
        x = copy.deepcopy(obj)
        fn(x)
        return x

    text, s = base["text"], base["search"]
    return [
        ("K0-unmodified", "control", {}, set()),
        ("K1-crlf", "control", {"text": text.replace("\n", "\r\n"), "scaffold": base["scaffold"].replace("\n", "\r\n")}, set()),
        ("K2-validation-keys-sorted", "control", {"validation": json.loads(json.dumps(base["validation"], sort_keys=True))}, set()),
        ("M1-contact-sign", "mutation", {"text": text.replace("\\omega h_W'(\\theta_j-\\psi)\n-\nv\\cdot u_{\\theta_j}.", "\\omega h_W'(\\theta_j-\\psi)\n+\nv\\cdot u_{\\theta_j}.", 1)}, {"A0"}),
        ("M2-no-fake-descent", "mutation", {"search": edit(s, lambda x: x["formulas"].update(rotation_affine="0"))}, {"W1"}),
        ("M3-sanity-area", "mutation", {"sanity": edit(base["sanity"], lambda x: x["circle_R2_constant_variation"].update(expected_area=12.0))}, {"W1"}),
        ("M4-steiner", "mutation", {"search": edit(s, lambda x: x["bodies"].update(steiner=1e-6))}, {"W2"}),
        ("M5-path", "mutation", {"search": edit(s, lambda x: x["bodies"].update(path=1e-6))}, {"W2"}),
        ("M6-reuleaux-cap", "mutation", {"search": edit(s, lambda x: x["caps"].update(reuleaux=0.1))}, {"W3"}),
        ("M7-square-open", "mutation", {"search": edit(s, lambda x: x["caps"]["square_blocked"].__setitem__(0, False))}, {"W3"}),
        ("M8-contact-off", "mutation", {"search": edit(s, lambda x: x["contact"].update(worst=1e-3))}, {"W4"}),
        ("M9-omega-idle", "mutation", {"search": edit(s, lambda x: x["contact"].update(omega_term=0.0))}, {"W4"}),
        ("M10-average-off", "mutation", {"search": edit(s, lambda x: x["kernel"]["generic"][0].update(average=x["kernel"]["generic"][0]["average"] + 0.01))}, {"W5"}),
        ("M11-lp-descent", "mutation", {"search": edit(s, lambda x: x["kernel"].update(lp_generic=-0.3))}, {"W5"}),
        ("M12-schema-variation", "mutation", {"schema": edit(base["schema"], lambda x: x.update(first_variation="delta A[h;phi] = integral h dmu_phi"))}, {"W6"}),
        ("M13-state-not-claimed", "mutation", {"state": edit(base["state"], lambda x: x["not_claimed"].append("global convergence of the shape flow"))}, {"W6"}),
        ("M14-scaffold-lp", "mutation", {"scaffold": base["scaffold"] + "\ndef descent_lp(rows):\n    from scipy.optimize import linprog\n    return linprog\n"}, {"W7"}),
        ("M15-scaffold-row-sign", "mutation", {"scaffold": base["scaffold"].replace("-math.cos(contact.theta)", "math.cos(contact.theta)", 1)}, {"W7"}),
        ("M16-validation-digest", "mutation", {"validation": {**base["validation"], "checks_sha256": "0" * 64}}, {"P1"}),
        ("M17-sanity-bytes", "mutation", {"bytes": {**base["bytes"], CHECKS_JSON: base["bytes"][CHECKS_JSON] + b" "}}, {"P1"}),
    ]


def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("--drill", action="store_true")
    ap.add_argument("--blind", metavar="CID")
    args = ap.parse_args()
    base = {**load_inputs(), "search": searches()}

    def build():
        return {**copy.deepcopy(base), **copy.deepcopy(CLAIMS)}

    if args.blind:
        CHECKS[args.blind] = lambda w: []
        args.drill = True
    w = build()
    results = {}
    for c, fn in CHECKS.items():
        try:
            results[c] = fn(w)
        except KeyError as e:
            results[c] = [f"claim missing: {e}"]
    for c, f in results.items():
        print(f"{c} {'PASS' if not f else 'FAIL'}", *(f"\n     {x}" for x in f[:6]))
    for key in ("w1", "w2", "w3", "w4", "w5", "w6", "w7", "p1"):
        print(key, "=", json.dumps(w.get(key), default=str, ensure_ascii=False))
    if not args.drill:
        return
    cases = drill_cases(base)
    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-053", "subject": "LESR Paper 03: support-function variational calculus", "document": DOC,
           "checks": {c: {"pass": not f, "failures": f} for c, f in results.items()},
           **{k: w.get(k) for k in ("w1", "w2", "w3", "w4", "w5", "w6", "w7", "p1")}, "drill": summary}
    LOG.write_text(json.dumps(log, ensure_ascii=False, indent=1, default=str) + "\n", encoding="utf-8")


if __name__ == "__main__":
    main()