"""RUN-054 - LESR Paper 04: skew-constrained shape flow and closure graph.

LESR Paper 04 (2026-09-20) builds the series' dynamics layer. A metric turns the first variation into a gradient -
L^2 gives h + h'' (the radius of curvature), the Sobolev H^1_l metric damps mode m by 1 / (1 + l^2 m^2), a general
preconditioner M gives M^-1 (h + h''); the L^2 gradient is orthogonal to translations. The flow is the metric
projection of a raw drive onto the linearized feasible cone K_h; a projection theorem makes any nonzero projected
velocity a first-order descent (delta A[h; V*] <= -||V*||^2), and V* = 0 is "local first-order projected
stationarity". Steps are capped four ways and checked by the exact quadratic along the step. The research is kept
as a typed closure graph (seven edge types, five proof statuses, dominator spine, SCC quotient). It ships a flow
schema, a graph 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, the projection theorem on random cones and
metrics, and the flow on the Lebesgue kernel with RUN-053's first-order models.

Checks
  A0 anchors      the formulas below occur in the document as written
  W1 formulas     sec. 3, 5-6, 8, 10, 33-35 in sympy; the sanity record's numbers
  W2 projection   sec. 15-17 on 200 random polyhedral cones with random metrics (exact projections by NNLS), cones
                  whose polar holds -G, and polytopes that hold 0 but are not cones
  W3 kernel       sec. 14, 17 on the kernel: V* = 0 exactly when the affine LP of Paper 03 finds no descent; that LP
                  at three placements, against the area's measured gradient
  W4 records      the flow and graph schemas and the state against the sections that carry them
  W5 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/run054_lesr_paper04.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 minimize, nnls

import gate
import run046_uesfcm_cycles_07_09 as r46
from corpus import SOURCES, TREE, source_bytes, walk_sources
from run052_lesr_paper02 import QDAG, integral_2pi
from run053_lesr_paper03 import ZPOL, first_order, lp_descent, structure

FOLDER = "main.zip!/LESR_Skew_Field_Variational/"
PKG = FOLDER + "LESR_Paper04_Package_v0.1.zip!/"
DOC = PKG + "LESR_Paper04_Skew_Constrained_Shape_Flow_and_Closure_Graph_v0.1.md"
GRAPH = PKG + "LESR_Paper04_Closure_Graph_Schema_v0.1.json"
FLOW = PKG + "LESR_Paper04_Flow_Schema_v0.1.json"
CHECKS_JSON = PKG + "LESR_Paper04_Sanity_Checks.json"
SCAFFOLD = PKG + "LESR_Paper04_Shape_Flow_Compiler_Scaffold.py"
STATE = PKG + "LESR_State_Paper04_v0.1.json"
VALIDATION = PKG + "LESR_Paper04_VALIDATION.json"
PAPER05 = FOLDER + "LESR_Paper05_Package_v0.1.zip!/LESR_Paper05_Contact_Pressure_Dual_Multipliers_and_No_Slack_Boundary_v0.1.md"
LOG = TREE / "data" / "gate-logs" / "run054-lesr-paper04.json"
MANIFEST = TREE / "data" / "SOURCE_MANIFEST.v1.json"
CACHE = TREE / "data" / "cache" / "run054-searches.json"     # data/cache is not committed
GENERIC = [QDAG[0] + 0.002, QDAG[1] + 0.001, QDAG[2] + 0.05, QDAG[3] - 0.001, QDAG[4] + 0.002]   # RUN-053's generic placement


def load_inputs():
    blobs = {n: source_bytes(n) for n in SOURCES}
    want = (DOC, GRAPH, FLOW, CHECKS_JSON, SCAFFOLD, STATE, VALIDATION)
    got = {leaf.id: d for leaf, d in walk_sources(blobs, keep=lambda l: l.id in want + (PAPER05,)) 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"), "graph": json.loads(got[GRAPH]), "flow": json.loads(got[FLOW]), "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}, "paper05": got[PAPER05].decode("utf-8")}


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

def formulas():
    """sec. 3: the boundary x = h u + h' u' has x' = (h + h'') u', so h + h'' is ds/d theta, the radius of curvature;
    sec. 5-6: the H^1_l gradient mode by mode on e^{i m theta}; sec. 8: a Fourier-diagonal preconditioner; sec. 10:
    int (h + h'') cos = int (h + h'') sin = 0; sec. 34-35: the decrease interval and the line minimizer; the sanity
    record recomputed"""
    t, l, e, g, q = sp.symbols("theta ell epsilon g q", real=True)
    m = sp.symbols("m", integer=True, nonzero=True)
    c = sp.symbols("c0:5", real=True)
    d = sp.symbols("d1:5", real=True)
    h = c[0] + sum(c[k] * sp.cos(k * t) + d[k - 1] * sp.sin(k * t) for k in range(1, 5))
    D = lambda x, n=1: sp.diff(x, t, n)
    u, up = sp.Matrix([sp.cos(t), sp.sin(t)]), sp.Matrix([-sp.sin(t), sp.cos(t)])
    x = h * u + D(h) * up
    wave = sp.exp(sp.I * m * t)
    out = {"radius": sp.simplify(D(x) - (h + D(h, 2)) * up) == sp.zeros(2, 1),
           "h1_multiplier": sp.simplify((wave - l ** 2 * D(wave, 2)) / wave - (1 + l ** 2 * m ** 2)) == 0,
           "translation": [sp.simplify(integral_2pi((h + D(h, 2)) * sp.cos(t), t)) == 0, sp.simplify(integral_2pi((h + D(h, 2)) * sp.sin(t), t)) == 0]}
    # a Fourier-diagonal preconditioner M = 1 + l^2 m^2 + l^4 m^4 (on each mode): G_M has <G_M, phi>_M = delta A
    phi = sum(sp.Symbol(f"p{k}") * sp.cos(k * t) for k in range(0, 4))
    Mk = lambda k: 1 + l ** 2 * k ** 2 + l ** 4 * k ** 4
    rho = h + D(h, 2)
    GM = sum(integral_2pi(rho * sp.cos(k * t), t) / (sp.pi * (2 if k == 0 else 1)) / Mk(k) * sp.cos(k * t) for k in range(0, 5)) + \
        sum(integral_2pi(rho * sp.sin(k * t), t) / sp.pi / Mk(k) * sp.sin(k * t) for k in range(1, 5))
    MGM = GM + l ** 2 * (-D(GM, 2)) + l ** 4 * D(GM, 4)
    out["preconditioner"] = sp.simplify(sp.expand(integral_2pi(MGM * phi, t) - integral_2pi(rho * phi, t))) == 0
    Aeps = g * e + q * e ** 2 / 2
    out["decrease"] = [sp.simplify(Aeps.subs(e, -2 * g / q)) == 0, sp.simplify(sp.diff(Aeps, e).subs(e, -g / q)) == 0]
    out = {k: ([bool(v) for v in val] if isinstance(val, list) else bool(val)) for k, val in out.items()}
    out["sanity"] = {"h1": [1 / (1 + 0.25 * mm * mm) for mm in (0, 1, 2, 4)], "line": [0.2 / 0.5, 1.0 - 0.2 * 0.4 + 0.5 * 0.5 * 0.4 ** 2]}
    return out


# ---------------------------------------------------------------- the projection theorem

def cone_projection(B, R, G):
    """V* = argmin over {V : B V >= 0} of 1/2 ||V + G||_M^2, M = R^T R: in w = R V the Euclidean projection of
    -R G onto {w : C w >= 0}, C = B R^-1, by Moreau - the polar is {-C^T lam : lam >= 0}, and its projection comes
    from nonnegative least squares"""
    C = B @ np.linalg.inv(R)
    y = -R @ G
    lam, _ = nnls(C.T, -y, maxiter=10000)
    w = y + C.T @ lam
    return np.linalg.solve(R, w), w, C


def projection_scan(trials=200, seed=54):
    """sec. 15-17: random cones {B V >= 0} (dimension 3-12, 2-20 rows), random metrics, random G: the identity
    <G, V*>_M = -||V*||_M^2 (Moreau: for a cone the inequality is an equality), feasibility; cones whose polar holds
    -G give V* = 0; polytopes {A V <= b}, b > 0, give the inequality, strictly in some cases"""
    rng = np.random.default_rng(seed)
    worst, infeas, polar_zero, nonzero = 0.0, 0.0, 0.0, 0
    for _ in range(trials):
        n, k = int(rng.integers(3, 13)), int(rng.integers(2, 21))
        B = rng.normal(size=(k, n))
        R = np.eye(n) + 0.3 * np.triu(rng.normal(size=(n, n)))
        R[np.diag_indices(n)] = np.abs(R[np.diag_indices(n)]) + 0.5
        G = rng.normal(size=n)
        V, w, C = cone_projection(B, R, G)
        g_dot = float((R @ G) @ w)
        worst = max(worst, abs(g_dot + float(w @ w)) / max(1.0, float(w @ w)))
        infeas = max(infeas, float(max(0.0, -np.min(C @ w))))
        nonzero += bool(np.linalg.norm(w) > 1e-9)
        # -G in the polar: G = R^-1 C^T lam with lam >= 0
        lam = np.abs(rng.normal(size=k))
        _, w0, _ = cone_projection(B, R, np.linalg.solve(R, C.T @ lam))
        polar_zero = max(polar_zero, float(np.linalg.norm(w0)))
    # polytopes that hold 0 in their interior: the inequality, and how often it is strict
    slack_min, strict = math.inf, 0
    for _ in range(60):
        n, k = 4, 9
        A, b = rng.normal(size=(k, n)), rng.uniform(0.2, 1.0, size=k)
        R = np.eye(n) + 0.3 * np.triu(rng.normal(size=(n, n)))
        R[np.diag_indices(n)] = np.abs(R[np.diag_indices(n)]) + 0.5
        M = R.T @ R
        G = 3 * rng.normal(size=n)
        res = minimize(lambda x: 0.5 * (x + G) @ M @ (x + G), np.zeros(n), jac=lambda x: M @ (x + G), method="SLSQP",
                       constraints=[{"type": "ineq", "fun": lambda x: b - A @ x, "jac": lambda x: -A}], options={"ftol": 1e-15, "maxiter": 500})
        V = res.x
        slack = -float(V @ M @ V) - float(G @ M @ V)       # >= 0 is the theorem
        slack_min = min(slack_min, slack)
        strict += slack > 1e-6
    return {"identity": worst, "infeasible": infeas, "nonzero": nonzero, "trials": trials, "polar_zero": polar_zero,
            "polytope_slack_min": slack_min, "polytope_strict": strict}


# ---------------------------------------------------------------- the kernel

def kernel_flow():
    """sec. 14, 17 on the kernel. For a cone, V* = P_K(-G) is 0 exactly when <G, V> >= 0 on K (the polar), i.e. when
    the LP min of the affine first-order change over K is 0 - RUN-053's LP (container free, witnesses re-placed).
    Solved here at the generic placement, at q-dagger and at the polish Z; the area's gradient at the generic one
    by central differences, its descent along the steepest direction, and what the affine model charges there"""
    A = r46.area
    out = {}
    for name, v in (("generic", GENERIC), ("qdag", QDAG), ("z", ZPOL)):
        st = structure(v)
        val, arg = lp_descent(st, v)
        out[name] = {"lp": val, "zero_argmin": all(x == 0.0 for x in arg)}
    E = np.eye(5)
    grad = np.array([(A([GENERIC[j] + 1e-7 * E[i][j] for j in range(5)]) - A([GENERIC[j] - 1e-7 * E[i][j] for j in range(5)])) / 2e-7 for i in range(5)])
    dn = -grad / np.linalg.norm(grad)
    st = structure(GENERIC)
    out["gradient"] = float(np.linalg.norm(grad))
    out["steepest"] = {"measured": (A([GENERIC[j] + 1e-7 * dn[j] for j in range(5)]) - A(GENERIC)) / 1e-7, "affine": first_order(st, GENERIC, dn)[1]}
    return out


# ---------------------------------------------------------------- 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
    import run053_lesr_paper03 as r53
    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__, r53.__file__):
        digest.update(open(path, "rb").read())
    digest.update(extra.encode())
    return digest.hexdigest()


HELPERS = {"formulas": (), "projection": (cone_projection,), "kernel": (r46.area, r46.hull)}


def searches():
    cached = json.loads(CACHE.read_text(encoding="utf-8")) if CACHE.exists() else {}
    out = {}
    for name, fn in (("formulas", formulas), ("projection", projection_scan), ("kernel", kernel_flow)):
        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"G_{L^2}=h+h''.", r"h+h''=\rho=\frac1\kappa.", r"(1-\ell^2\partial_\theta^2)G_{H^1}=h+h''.", r"(G_{H^1})_m=\frac{\rho_m}{1+\ell^2m^2}.",
           r"G_M=M_h^{-1}(h+h'').", r"\int_0^{2\pi}\rho\cos\theta\,d\theta=0,", r"V^\star=\arg\min_{V\inK_h}\frac12\|V+G_{\mathcalM}\|_{\mathcalM}^2.",
           r"\deltaA[h;V^\star]\le-\|V^\star\|_{\mathcalM}^2.", r"\deltaA[h;V^\star]<0.", r"0<\varepsilon<\min\left(\varepsilon_{\max},-\frac{2g}{q}\right).",
           r"\varepsilon=\min(\varepsilon^\star,\varepsilon_{\max})", r"\varepsilon_k=\varepsilon_02^{-k}.", r"\min_c\frac12c^TMc+g^Tc",
           r"\mathcalG_{\mathrm{closure}}=(V,E,\tau).", r"\DeltaB_n=\widetildeB_n^+-B_n^-.", r"\text{directedacyclicclosuregraph}."],
}
STATE_MAP = {"L2 shape metric gradient": (3, "$L^2$Metric"), "H1 Sobolev preconditioned gradient": (5, "$H^1$/SobolevMetric"),
             "general metric/preconditioner form": (8, "一般PreconditionedMetric"), "metric projection onto local feasible cone": (14, "MetricProjection"),
             "projected descent inequality": (15, "Theorem—ProjectedDescentInequality"), "finite-mode QP form": (19, "QPForm"),
             "certified finite-step caps": (28, "四個StepCap"), "exact quadratic area line-step rule": (33, "ExactAreaStepTest"),
             "backtracking certificate policy": (36, "BacktrackingCertificate"), "hybrid continuous/discrete shape dynamics": (41, "HybridFlow"),
             "expansion vs closure operating modes": (23, "OptimizationMode與ExplorationMode分開"), "Shape Closure Graph": (44, "ShapeClosureGraph"),
             "typed graph edges": (46, "EdgeTypes"), "reachable possibility graph": (48, "ReachablePossibilitySubgraph"),
             "mandatory/dominator spine": (49, "MandatorySpine"), "computation-deferred state": (62, "HeavyCompute不阻塞理論主線"),
             "local recomputation handoff schema": (61, "ComputationHandoff"), "visualization three-view plan": (65, "三層最終視圖")}
CLAIMS = {
    "formulas_claim": {"radius": True, "h1_multiplier": True, "translation": [True, True], "preconditioner": True, "decrease": [True, True],
                       "sanity": [[True, True, True, True], True, True]},
    "projection_claim": {"cone_identity": True, "feasible": True, "nonzero": 126, "trials": 200, "polar_zero": True, "polytope_inequality": True,
                         "polytope_strict": 60},
    "kernel_claim": {"lp": {"generic": [True, True], "qdag": [True, True], "z": [True, True]}, "gradient": 0.39, "steepest": [-0.39, 0.107]},
    "records_claim": {"edge_types": True, "proof_statuses": True, "views": True, "node_required": True, "events_in_text": ["CERT_PRUNE", "COMPUTE_DEFERRED"],
                      "compute_in_text": ["COMPUTE_DEFERRED"], "step_caps": True, "unmapped": [], "completed": 18, "not_claimed": False, "next": [True, False]},
    "scaffold_facts": {"lines": 76, "defs": ["h1_precondition", "area_step", "guaranteed_decrease_cap", "qp_payload"],
                       "classes": ["FourierComponent", "GraphNode", "GraphEdge", "ClosureGraph"], "solver": False, "graph_algorithms": False,
                       "edge_types_checked": False, "decrease_rule": True},
    "package_claim": {"digests": [True] * 6, "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):
    s, rec = w["search"]["formulas"], w["sanity"]
    same = lambda x, y: abs(x - y) <= 1e-12
    ex = rec["exact_area_decrease_example"]
    got = {k: s[k] for k in ("radius", "h1_multiplier", "translation", "preconditioner", "decrease")}
    got["sanity"] = [[same(r["multiplier"], v) for r, v in zip(rec["H1_multiplier_examples"], s["sanity"]["h1"])],
                     same(ex["unconstrained_line_minimizer"], s["sanity"]["line"][0]), same(ex["A_at_minimizer"], s["sanity"]["line"][1])]
    w["w1"] = got
    return [] if got == w["formulas_claim"] else [f"the formulas: {got}"]


def w2_projection(w):
    """sec. 15: <G, V*> <= -||V*||^2 for 0 in a closed convex set; for a cone it is an equality"""
    s = w["search"]["projection"]
    got = {"cone_identity": s["identity"] < 1e-9, "feasible": s["infeasible"] < 1e-9, "nonzero": s["nonzero"], "trials": s["trials"],
           "polar_zero": s["polar_zero"] < 1e-9, "polytope_inequality": s["polytope_slack_min"] > -1e-7, "polytope_strict": s["polytope_strict"]}
    w["w2"] = got
    return [] if got == w["projection_claim"] else [f"the projection: {got}"]


def w3_kernel(w):
    s = w["search"]["kernel"]
    got = {"lp": {k: [s[k]["lp"] == 0.0, s[k]["zero_argmin"]] for k in ("generic", "qdag", "z")}, "gradient": sig(s["gradient"], 3),
           "steepest": [sig(s["steepest"]["measured"], 3), sig(s["steepest"]["affine"], 3)]}
    w["w3"] = got
    return [] if got == w["kernel_claim"] else [f"the kernel: {got}"]


def w4_records(w):
    """the schemas' lists against the sections that name them; the state's items against their sections"""
    fl, gr, st, raw, text = w["flow"], w["graph"], w["state"], w["text"], gate.squash(w["text"])
    edge_types = re.findall(r"^##\s*`([A-Z_]+)`", raw, re.M)
    proof = re.search(r"proof_status:[\s]*(.*?)```", raw, re.S)
    got = {"edge_types": gr["edge_types"] == edge_types, "proof_statuses": proof is not None and gr["proof_statuses"] == proof.group(1).split(),
           "views": [v.replace("_", "") for v in gr["views"]] == ["SHAPEVIEW", "STATEGRAPH", "CLOSURESPINE"]
           and all(x in text for x in ("ViewA—ShapeView", "ViewB—StateGraph", "ViewC—ClosureSpine")),
           "node_required": all(f in (re.search(r"node:[\s]*(.*?)edge:", raw, re.S).group(1).split()) for f in gr["node"]["required"]),
           "events_in_text": [e for e in fl["events"] if e in raw], "compute_in_text": [c for c in gr["compute_statuses"] if c in raw],
           "step_caps": fl["step_caps"] == ["convexity", "topology", "corridor", "witness"],
           "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"]), "not_claimed": "not_claimed" in st,
           # the state names Paper 05 as sec. 74 does; Paper 05's own title reads "... and the No-Slack Boundary Principle"
           "next": [st["next"] == "Paper 05 — Contact Pressure, Dual Multipliers, and No-Slack Boundary" and "ContactPressure,DualMultipliers,andNo-SlackBoundary}" in text,
                    "Contact Pressure, Dual Multipliers, and No-Slack Boundary" in w["paper05"][:300]]}
    w["w4"] = got
    return [] if got == w["records_claim"] else [f"the records: {got}"]


def w5_scaffold(w):
    """the scaffold, read: whether it solves the QP; whether it computes dominators or the SCC quotient (sec. 49-50,
    70); whether it checks an edge's type against sec. 46's seven; sec. 34's decrease cap"""
    code = w["scaffold"]
    seven = ("ORBIT", "SKEW_EXPAND", "FLOW", "TOPOLOGY_JUMP", "CERT_PRUNE", "MERGE", "REOPEN")
    got = {"lines": code.count("\n"), "defs": re.findall(r"^def (\w+)", code, re.M), "classes": re.findall(r"^class (\w+)", code, re.M),
           "solver": bool(re.search(r"scipy|cvxpy|quadprog|osqp|minimize|linprog", code)),
           "graph_algorithms": bool(re.search(r"domin|scc|strongly|topolog", code, re.I)), "edge_types_checked": any(f'"{e}"' in code for e in seven),
           "decrease_rule": bool(re.search(r"return min\(eps_max, -2\.0\*dA/d2A\)", code))}
    w["w5"] = 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["graph_schema_sha256"] == sha(B[GRAPH]), v["flow_schema_sha256"] == sha(B[FLOW]),
                       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_projection, "W3": w3_kernel, "W4": w4_records, "W5": w5_scaffold, "P1": p1_package}
DEPENDS = {"A0": {"text"}, "W1": {"search", "sanity", "formulas_claim"}, "W2": {"search", "projection_claim"}, "W3": {"search", "kernel_claim"},
           "W4": {"flow", "graph", "state", "text", "paper05", "records_claim"}, "W5": {"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-projection-sign", "mutation", {"text": text.replace("\\delta A[h;V^\\star]\n\\le\n", "\\delta A[h;V^\\star]\n\\ge\n", 1)}, {"A0"}),
        ("M2-radius", "mutation", {"search": edit(s, lambda x: x["formulas"].update(radius=False))}, {"W1"}),
        ("M3-sanity-multiplier", "mutation", {"sanity": edit(base["sanity"], lambda x: x["H1_multiplier_examples"][2].update(multiplier=0.4))}, {"W1"}),
        ("M4-cone-identity", "mutation", {"search": edit(s, lambda x: x["projection"].update(identity=1e-3))}, {"W2"}),
        ("M5-polar-nonzero", "mutation", {"search": edit(s, lambda x: x["projection"].update(polar_zero=0.1))}, {"W2"}),
        ("M6-generic-descent", "mutation", {"search": edit(s, lambda x: x["kernel"]["generic"].update(lp=-0.3))}, {"W3"}),
        ("M7-flat-gradient", "mutation", {"search": edit(s, lambda x: x["kernel"].update(gradient=0.0))}, {"W3"}),
        ("M8-edge-type", "mutation", {"graph": edit(base["graph"], lambda x: x["edge_types"].append("SPLIT"))}, {"W4"}),
        ("M9-state-item", "mutation", {"state": edit(base["state"], lambda x: x["completed"].append("global convergence proof"))}, {"W4"}),
        ("M10-scaffold-dominators", "mutation", {"scaffold": base["scaffold"] + "\ndef dominators(graph, root):\n    return {root: {root}}\n"}, {"W5"}),
        ("M11-scaffold-cap", "mutation", {"scaffold": base["scaffold"].replace("-2.0*dA/d2A", "-dA/d2A", 1)}, {"W5"}),
        ("M12-validation-digest", "mutation", {"validation": {**base["validation"], "flow_schema_sha256": "0" * 64}}, {"P1"}),
        ("M13-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", "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-054", "subject": "LESR Paper 04: skew-constrained shape flow and closure graph", "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", "p1")}, "drill": summary}
    LOG.write_text(json.dumps(log, ensure_ascii=False, indent=1, default=str) + "\n", encoding="utf-8")


if __name__ == "__main__":
    main()