"""RUN-014 - Round 13 (v0.13): the five-parameter smooth event-KKT system - the four equal branches, the gradients
and pressures, the 12 x 12 Jacobian, the restarts, the curvature ledger and the full phase audit.

Round 13 frees the layer width eps beside (w, beta, delta, c) and solves twelve equations in twelve unknowns
(w, beta, delta, c, eps, phi1, phi2, s, mu1..mu4): the two smooth low branches and the 120- and 270-degree cusps all
equal to s, the low branches stationary in phase, the pressures mu balancing the four branches' parameter gradients
and summing to one. It reports the common scale 0.9989143432974851, 4.2128524e-9 above Round 08's; a Jacobian of full
rank but "highly ill-conditioned" (condition 3.77e7); 20 restarts back in one basin; the curvature ledger; and 8 local
minima on the phase circle, the four outside the system at 1.0293 and above.

Checks (each returns a list of failures; empty = pass)
  X1 parameters   Round 08's reference is Round 09's decimal curve and its scale the double nearest that curve's
                  s(3pi/2) (200 digits); the candidate's centre length, alpha, degrees and parameter shifts are its own
                  arithmetic
  X2 branches     at 60 digits: each low phase is a smooth minimum of s with the supports L|p0|p2 (to 1e-15); the
                  120- and 270-degree phases; the four branch values (to 2e-16); the spread, the common scale (the
                  least), the improvement over the reference and the distance to one
  X3 gradients    the stated gradients are the branches' partial derivatives at fixed phase (central differences at
                  60 digits, step 1e-20), to 2e-7
  X4 KKT          the residual is the stated pressures times the stated gradients and its norm is its norm; the
                  pressures sum to one within their rounding; this run's own solution of the twelve equations, by chord
                  steps from the stated point, converges (|F| < 1e-35) with positive pressures, within 1e-7 of the
                  stated parameters and phases, 2e-6 of the pressures and 1e-13 of the scale
  X5 Jacobian     at that solution, the Jacobian (second differences at 60 digits, steps 1e-15 and 1e-12, agreeing to
                  1e-9) has rank 12, its smallest singular value above 1e-3; the stated singular values and condition
                  number are its own (fixed-text failure); the Lagrangian's Hessian is negative definite on the plane
                  the four gradients leave free, so the solution is a strict local maximum of the least of the four
  X6 copies       the data files hold the summary's records; the restarts' count, range and largest distance; the
                  manifest's figures and files
  X7 curvature    the total turn is delta; the peak is the exact peak, to 1e-14 relative; the 5-95 % width is the
                  exact width, to 2e-11
  X8 phase        this run's own phase map of the stated curve (2^18 phases, the switches by bisection at 60 digits,
                  the zeros of each interval's derivative): 8 local minima, as stated; the controlling four are the
                  branches; each of the other four is the exact minimum near it, to 5e-6 in phase and 1e-9 in scale,
                  of the stated kind (fixed-text failures)
  X9 documents    the README, the report's figures and tables, the audit note and the handoff state the data's;
                  section 7's projection is the stated gradients' (fixed-text failure); the printed digits of
                  the pressures, the gradients and the width against the quantities (fixed-text failure)

Drill: each plant names the checks that must trip - trip meaning a failure the unplanted package does not already
have; controls must trip nothing.

Usage:
  python code/run014_round13_smooth_kkt.py [--drill [--log PATH] | --blind X#]
"""
import argparse
import copy
import csv
import io
import json
import math
import pathlib
import re

import mpmath as mp
import numpy as np

from corpus import TREE, sha256, source_bytes
import gate
from run001_corpus_identity import restore_escapes
from run004_round03_symmetric_links import load
from run009_round08_smooth_layer import Wing, ledger
import run010_round09_precision_cusp as r10
import run011_round10_phase_map as r11
from run010_round09_precision_cusp import cusp_270, dec, hi, iv, lo

RID = "round-v0.13"
LOG = TREE / "data" / "gate-logs" / "run014-round13-smooth-kkt.json"
KEYS = ("wing_length", "beta_radians", "delta_radians", "transition_center", "epsilon")
DKEYS = ("wing_length", "beta", "delta", "transition_center", "epsilon")
BR = ("B1_low", "B2_low", "B3_120deg", "B4_270deg")
# the supports (min_x, min_y, max_d) of each branch; at the two cusps either side's give the same value for every curve
SIG = ("L|p0|p2", "L|p0|p2", "p2|p1|p3", "p3|p3|p1")


# ---------------------------------------------------------------- the world

def build_world():
    f = load(RID)
    md = lambda n: restore_escapes(f[n].decode("utf-8"))
    js = lambda n: json.loads(f[n])
    rows = lambda n: list(csv.DictReader(io.StringIO(f[n].decode("utf-8"))))
    return {
        "readme": md("README.md"), "report": md("reports/ROUND_13_REPORT.md"), "audit_note": md("reports/ROUND_13_KKT_AUDIT.md"),
        "handoff": md("reports/ROUND_14_CURVATURE_DISTRIBUTION_HANDOFF.md"), "source": f["src/smooth_kkt_replay.py"].decode("utf-8"),
        "manifest": js("MANIFEST.json"), "members": sorted(f),
        "summary": js("data/round13_summary.json"), "candidate": js("data/smooth_kkt_candidate.json"),
        "minima": js("data/full_phase_minima.json"), "basin": rows("data/basin_restarts.csv"),
        "gradients": rows("data/branch_gradient_matrix.csv"),
        "round09": json.loads(load("round-v0.9")["data/round09_summary.json"]),
        "round05": json.loads(load("round-v0.5")["data/round05_summary.json"])["event_solution"]["parameters"],
        "legs": ("1", "1"),
    }


def event(w):
    E = w["round05"]
    return {"l1": repr(float(E["l1"])), "l2": repr(float(E["l2"])), "beta": repr(float(E["beta_radians"])),
            "delta": repr(float(E["delta_radians"]))}


def hp_at(w, pv):
    """the world RUN-010's high-precision tools read, for the curve with parameters pv (decimal strings, KEYS order)"""
    return {"summary": {"definitions": {"candidate": dict(zip(DKEYS, pv)), "event_control": event(w)}},
            "legs": w["legs"], "cusp_formula": cusp_270}


def stated(w):
    P = w["summary"]["round13_candidate"]["parameters"]
    return [repr(float(P[k])) for k in KEYS]


def hp(w):
    """the stated curve: the summary's parameters as their shortest decimals"""
    return hp_at(w, stated(w))


def hp08(w):
    """Round 08's curve as Round 09 fixes it in decimals"""
    return {"summary": {"definitions": {"candidate": w["round09"]["definitions"]["candidate"], "event_control": event(w)}},
            "legs": w["legs"], "cusp_formula": cusp_270}


def value(C, sig, ph):
    """s at a phase with a signature's supports"""
    with mp.workdps(C.dps):
        cs, sn = mp.cos(ph), mp.sin(ph)
        tot = mp.mpf(0)
        for name, lab in zip(r11.NAMES, sig.split("|")):
            a = C.a[name]
            p = r11.point(C, lab, name, ph)
            tot += (a[0] * cs + a[1] * sn) * p[0] + (-a[0] * sn + a[1] * cs) * p[1]
        return tot


def pts(pv):
    with mp.workdps(60):
        return [mp.nstr(x, 55) for x in pv]


def vat(w, pv, r, ph):
    """branch r's value at a phase, for the curve with parameters pv (60-digit numbers)"""
    return value(r10.cand(hp_at(w, pts(pv)), 60), SIG[r], ph)


def low_phase(C, ph0):
    """the smooth minimum of s with the supports L|p0|p2 within 1e-4 of ph0, by bisection on the derivative's sign"""
    with mp.workdps(60):
        a, b = mp.mpf(ph0) - mp.mpf("1e-4"), mp.mpf(ph0) + mp.mpf("1e-4")
        if not r11.deriv(C, SIG[0], a) < 0 < r11.deriv(C, SIG[0], b):
            return None
        for _ in range(140):
            m = (a + b) / 2
            if r11.deriv(C, SIG[0], m) < 0:
                a = m
            else:
                b = m
        return (a + b) / 2


def branches(w):
    """the four branches of the stated curve at 60 digits: the low phases, the values, the low supports"""
    def go():
        C = r10.cand(hp(w), 60)
        R = w["summary"]["round13_candidate"]["branch_phases"]
        with mp.workdps(60):
            ph = [low_phase(C, R["B1_low"]), low_phase(C, R["B2_low"]), 2 * mp.pi / 3, 3 * mp.pi / 2]
            vals = [None if p is None else value(C, SIG[r], p) for r, p in enumerate(ph)]
            sigs = [None if p is None else r11.sig_at(C, p) for p in ph[:2]]
        return {"phases": ph, "values": vals, "signatures": sigs}
    return r10.memo("branches14", hp(w))(go)


def grads(w, pv, phis, h):
    """the four branches' partial derivatives in the five parameters at fixed phases, by central differences"""
    with mp.workdps(60):
        h = mp.mpf(h)
        out = []
        for r in range(4):
            row = []
            for i in range(5):
                up = [pv[k] + (h if k == i else 0) for k in range(5)]
                dn = [pv[k] - (h if k == i else 0) for k in range(5)]
                row.append((vat(w, up, r, phis[r]) - vat(w, dn, r, phis[r])) / (2 * h))
            out.append(row)
        return out


def exact_gradients(w):
    def go():
        B = branches(w)
        with mp.workdps(60):
            return grads(w, [mp.mpf(v) for v in stated(w)], B["phases"], "1e-20")
    return r10.memo("grads14", hp(w))(go)


# ---------------------------------------------------------------- the event-KKT system

def kkt_F(w, x):
    """the twelve equations: m_r - s, the low branches' phase derivatives, sum mu_r grad m_r, sum mu - 1"""
    with mp.workdps(60):
        p, phis, s, mu = x[:5], [x[5], x[6], 2 * mp.pi / 3, 3 * mp.pi / 2], x[7], x[8:12]
        C = r10.cand(hp_at(w, pts(p)), 60)
        G = grads(w, p, phis, "1e-20")
        out = [value(C, SIG[r], phis[r]) - s for r in range(4)]
        out += [r11.deriv(C, SIG[r], phis[r]) for r in (0, 1)]
        out += [mp.fsum(mu[r] * G[r][i] for r in range(4)) for i in range(5)]
        out += [mp.fsum(mu) - 1]
        return out


def kkt_J(w, x, h):
    """the Jacobian of the twelve equations, by first and second central differences"""
    with mp.workdps(60):
        h = mp.mpf(h)
        p, phis, mu = x[:5], [x[5], x[6], 2 * mp.pi / 3, 3 * mp.pi / 2], x[8:12]
        mv = lambda d: [p[k] + d[k] for k in range(5)]
        e = lambda i, s_: [s_ * h if k == i else 0 for k in range(5)]
        add = lambda a, b: [a[k] + b[k] for k in range(5)]
        v = lambda d, r: vat(w, mv(d), r, phis[r])
        D = lambda d, r, ph: r11.deriv(r10.cand(hp_at(w, pts(mv(d))), 60), SIG[r], ph)
        z = [0] * 5
        G = grads(w, p, phis, h)
        Hs = []
        for r in range(4):
            H = [[None] * 5 for _ in range(5)]
            for i in range(5):
                for j in range(i, 5):
                    if i == j:
                        t = (v(e(i, 1), r) - 2 * v(z, r) + v(e(i, -1), r)) / h ** 2
                    else:
                        t = (v(add(e(i, 1), e(j, 1)), r) - v(add(e(i, 1), e(j, -1)), r)
                             - v(add(e(i, -1), e(j, 1)), r) + v(add(e(i, -1), e(j, -1)), r)) / (4 * h * h)
                    H[i][j] = H[j][i] = t
            Hs.append(H)
        Mx = [[(D(e(i, 1), r, phis[r]) - D(e(i, -1), r, phis[r])) / (2 * h) for i in range(5)] for r in (0, 1)]
        Dpp = [(D(z, r, phis[r] + h) - D(z, r, phis[r] - h)) / (2 * h) for r in (0, 1)]
        J = mp.zeros(12, 12)
        for r in range(4):
            for i in range(5):
                J[r, i] = G[r][i]
            if r < 2:
                J[r, 5 + r] = D(z, r, phis[r])
            J[r, 7] = -1
        for r in (0, 1):
            for i in range(5):
                J[4 + r, i] = Mx[r][i]
            J[4 + r, 5 + r] = Dpp[r]
        for i in range(5):
            for j in range(5):
                J[6 + i, j] = mp.fsum(mu[r] * Hs[r][i][j] for r in range(4))
            for r in (0, 1):
                J[6 + i, 5 + r] = mu[r] * Mx[r][i]
            for r in range(4):
                J[6 + i, 8 + r] = G[r][i]
        for r in range(4):
            J[11, 8 + r] = 1
        return {"J": J, "G": G, "H": Hs, "M": Mx, "D": Dpp}


def singular_values(J):
    return np.linalg.svd(np.array([[float(J[a, b]) for b in range(12)] for a in range(12)]), compute_uv=False)


def kkt(w):
    """this run's solution of the twelve equations, by chord steps from the stated point (the Jacobian there), and the
    Jacobian and the second-order test at it"""
    def go():
        S = w["summary"]
        R = S["round13_candidate"]
        with mp.workdps(60):
            x = [mp.mpf(v) for v in stated(w)] + [mp.mpf(repr(float(R["branch_phases"][n]))) for n in BR[:2]] \
                + [mp.mpf(repr(float(R["common_scale"])))] + [mp.mpf(repr(float(m))) for m in S["smooth_event_kkt"]["branch_pressures"]]
            J0 = kkt_J(w, x, "1e-15")["J"]
            Fx = kkt_F(w, x)
            steps = 0
            while max(abs(f) for f in Fx) > mp.mpf("1e-35") and steps < 12:
                d = mp.lu_solve(J0, mp.matrix(Fx))
                x = [x[k] - d[k] for k in range(12)]
                Fx = kkt_F(w, x)
                steps += 1
            res = max(abs(f) for f in Fx)
            A, B = kkt_J(w, x, "1e-15"), kkt_J(w, x, "1e-12")
            # the low branches are minima over phase, so their values have the Hessian m_pp - m_pphi m_pphi^T / m_phiphi
            mu = x[8:12]
            HL = mp.zeros(5, 5)
            for r in range(4):
                for i in range(5):
                    for j in range(5):
                        c = A["H"][r][i][j] - (A["M"][r][i] * A["M"][r][j] / A["D"][r] if r < 2 else 0)
                        HL[i, j] += mu[r] * c
            Gm = mp.matrix([[A["G"][r][i] for i in range(5)] for r in range(4)])
            _, Sg, V = mp.svd_r(Gm, full_matrices=True)
            T = [V[3, :], V[4, :]]
            Rm = mp.matrix(2, 2)
            for a in range(2):
                for b in range(2):
                    Rm[a, b] = mp.fsum(T[a][0, i] * HL[i, j] * T[b][0, j] for i in range(5) for j in range(5))
            ev = sorted(mp.eigsy(Rm)[0])
            null_ok = max(abs(mp.fsum(Gm[r, i] * t[0, i] for i in range(5))) for r in range(4) for t in T)
        return {"x": x, "steps": steps, "residual": res, "sv": singular_values(A["J"]), "sv_coarse": singular_values(B["J"]),
                "gradient_sv": list(Sg), "null_residual": null_ok, "reduced_hessian": ev, "phase_curvature": A["D"]}
    return r10.memo("kkt14", hp(w))(go)


# ---------------------------------------------------------------- the phase map

def phase_map(w):
    """this run's phase map of the stated curve: the switches of a 2^18-phase double-precision scan, each fixed by
    bisection on the changing support's label at 60 digits; the zeros of each interval's derivative (65 equal and 30
    geometric steps toward the ends, then bisection); the local minima - cusps (D- < 0 < D+) and smooth minima"""
    def go():
        P = w["summary"]["round13_candidate"]["parameters"]
        ws = {"summary": {"curve_parameters": {d: P[k] for k, d in zip(KEYS, DKEYS)}}, "legs": w["legs"]}
        ph, _, sig = r11.scan(ws, 2 ** 18)
        ch = np.flatnonzero(sig != np.roll(sig, 1))
        C = r10.cand(hp(w), 60)
        sw = []
        with mp.workdps(60):
            for i in ch:
                L, R = str(sig[i - 1]).split("|"), str(sig[i]).split("|")
                k = [j for j in range(3) if L[j] != R[j]][0]
                if i == 0:
                    q = mp.mpf(0)          # min_y passes from p3 to p0 at phase 0 for every curve: the bottom side is level
                else:
                    a, b = mp.mpf(ph[i - 1]), mp.mpf(ph[i])
                    for _ in range(90):
                        m = (a + b) / 2
                        if C.terms(m)[r11.NAMES[k]][1] == [L[k]]:
                            a = m
                        else:
                            b = m
                    q = (a + b) / 2
                corner = L[k].startswith("p") and R[k].startswith("p")
                sw.append({"phase": q, "left": "|".join(L), "right": "|".join(R), "corner": corner,
                           "D-": r11.deriv(C, "|".join(L), q), "D+": r11.deriv(C, "|".join(R), q), "scale": C.s(q)})
            stat = []
            for j, s in enumerate(sw):
                a, b = s["phase"], sw[(j + 1) % len(sw)]["phase"] + (2 * mp.pi if j + 1 == len(sw) else 0)
                fr = sorted({mp.mpf(k) / 64 for k in range(65)} | {mp.mpf(10) ** -e for e in range(1, 16)}
                            | {1 - mp.mpf(10) ** -e for e in range(1, 16)})
                vals = [(a + f * (b - a), r11.deriv(C, s["right"], a + f * (b - a))) for f in fr]
                for (p0, d0), (p1, d1) in zip(vals, vals[1:]):
                    if (d0 > 0) != (d1 > 0):
                        x_, y_ = p0, p1
                        for _ in range(110):
                            m = (x_ + y_) / 2
                            if (r11.deriv(C, s["right"], m) > 0) == (d0 > 0):
                                x_ = m
                            else:
                                y_ = m
                        q = (x_ + y_) / 2
                        stat.append({"phase": q, "scale": C.s(q), "kind": "minimum" if d0 < 0 else "maximum", "signature": s["right"],
                                     "to_next_switch": b - q})
            mins = [{"source": "smooth", "phase": e["phase"], "scale": e["scale"], "to_next_switch": e["to_next_switch"]}
                    for e in stat if e["kind"] == "minimum"]
            mins += [{"source": "contact switch", "phase": s["phase"], "scale": s["scale"]} for s in sw if s["D-"] < 0 < s["D+"]]
        return {"switches": sw, "stationary": stat, "minima": sorted(mins, key=lambda m: m["scale"])}
    return r10.memo("phase14", hp(w))(go)


# ---------------------------------------------------------------- the checks

def x1_parameters(w):
    failures = []
    S = w["summary"]
    ref = S["round8_reference"]
    D = w["round09"]["definitions"]["candidate"]
    for k, d in zip(KEYS, DKEYS):
        if ref["parameters"][k] != float(D[d]):
            failures.append(f"the reference's {k} is not Round 09's decimal")
    if float(r10.s270(hp08(w))) != ref["scale"]:
        failures.append("the reference scale is not the double nearest the reference curve's s(3pi/2)")
    P = S["round13_candidate"]["parameters"]
    if abs(P["center_length"] - (1 - 2 * P["wing_length"])) > 2.3e-16 or abs(P["alpha_radians"] - (P["beta_radians"] + P["delta_radians"])) > 4.5e-16:
        failures.append("the centre length or alpha is not the candidate's")
    for k in ("beta", "delta", "alpha"):
        if abs(P[f"{k}_degrees"] / math.degrees(P[f"{k}_radians"]) - 1) > 1e-15:
            failures.append(f"{k} in degrees is not {k} in radians")
    for k in KEYS:
        if S["round13_candidate"]["parameter_shift"][k] != P[k] - ref["parameters"][k]:
            failures.append(f"the shift in {k} is not the candidate's less the reference's")
    return failures


def x2_branches(w):
    failures = []
    B = branches(w)
    R = w["summary"]["round13_candidate"]
    with mp.workdps(60):
        for j, n in enumerate(BR[:2]):
            if B["phases"][j] is None or abs(mp.mpf(repr(R["branch_phases"][n])) - B["phases"][j]) > mp.mpf("1e-15") \
                    or B["signatures"][j] != SIG[j]:
                failures.append(f"the {n} phase is not a smooth minimum of s with the supports {SIG[j]}")
        # the doubles of 2pi/3 and 3pi/2 as double arithmetic gives them: within one unit in the last place
        for n, e in (("B3_120deg", 2 * mp.pi / 3), ("B4_270deg", 3 * mp.pi / 2)):
            if abs(mp.mpf(R["branch_phases"][n]) - e) > math.ulp(R["branch_phases"][n]):
                failures.append(f"the {n} phase is not its cusp's")
        if None in B["values"]:
            return failures + ["the branch values cannot be computed"]
        for n, v in zip(BR, B["values"]):
            if abs(mp.mpf(repr(R["branch_values"][n])) - v) > mp.mpf("2e-16"):
                failures.append(f"the {n} value is not the branch's")
        least = min(B["values"])
        if abs(mp.mpf(repr(R["branch_spread"])) - (max(B["values"]) - least)) > mp.mpf("2e-16"):
            failures.append("the spread is not the branches'")
        if R["common_scale"] != min(R["branch_values"].values()) or abs(mp.mpf(repr(R["common_scale"])) - least) > mp.mpf("2e-16"):
            failures.append("the common scale is not the least branch")
        ref = w["summary"]["round8_reference"]["scale"]
        if R["improvement_over_round8"] != R["common_scale"] - ref \
                or abs(mp.mpf(repr(R["improvement_over_round8"])) - (least - r10.s270(hp08(w)))) > mp.mpf("3e-16"):
            failures.append("the improvement is not the least branch less the reference's s(3pi/2)")
        if R["distance_to_one"] != 1 - R["common_scale"]:
            failures.append("the distance to one is not one less the common scale")
    return failures


def x3_gradients(w):
    failures = []
    G = exact_gradients(w)
    Gs = w["summary"]["smooth_event_kkt"]["gradient_matrix"]
    for r, n in enumerate(BR):
        for i, k in enumerate(DKEYS):
            if abs(Gs[r][i] - float(G[r][i])) > 2e-7:
                failures.append(f"the {n} gradient in {k} is not the branch's partial derivative")
    return failures


def x4_kkt(w):
    failures = []
    K = w["summary"]["smooth_event_kkt"]
    mu, Gs, res = K["branch_pressures"], K["gradient_matrix"], K["weighted_stationarity_residual"]
    if any(abs(sum(mu[r] * Gs[r][i] for r in range(4)) - res[i]) > 1e-16 for i in range(5)):
        failures.append("the residual is not the pressures times the gradients")
    if abs(K["weighted_stationarity_norm"] / math.sqrt(sum(x * x for x in res)) - 1) > 1e-12:
        failures.append("the residual's norm is not its norm")
    if abs(sum(mu) - 1) > 2e-8 or (K["unknown_count"], K["equation_count"]) != (12, 12):
        failures.append("the pressures do not sum to one within their rounding, or the system is not 12 by 12")
    X = kkt(w)
    x = X["x"]
    R = w["summary"]["round13_candidate"]
    with mp.workdps(60):
        if X["residual"] > mp.mpf("1e-35") or min(x[8:12]) <= 0:
            failures.append("this run's chord steps from the stated point do not reach a solution with positive pressures")
            return failures
        if max(abs(mp.mpf(v) - x[k]) for k, v in enumerate(stated(w))) > mp.mpf("1e-7") \
                or max(abs(mp.mpf(repr(R["branch_phases"][n])) - x[5 + j]) for j, n in enumerate(BR[:2])) > mp.mpf("1e-7"):
            failures.append("the stated parameters or low phases are not within 1e-7 of the solution")
        if max(abs(mp.mpf(repr(m)) - x[8 + r]) for r, m in enumerate(mu)) > mp.mpf("2e-6"):
            failures.append("the stated pressures are not within 2e-6 of the solution's")
        if abs(mp.mpf(repr(R["common_scale"])) - x[7]) > mp.mpf("1e-13"):
            failures.append("the common scale is not within 1e-13 of the solution's")
    return failures


def x5_jacobian(w):
    failures = []
    X = kkt(w)
    K = w["summary"]["smooth_event_kkt"]
    sv, sv2 = X["sv"], X["sv_coarse"]
    if X["residual"] > mp.mpf("1e-35"):
        return ["there is no solution to test"]
    if max(abs(sv - sv2) / sv) > 1e-9 or sv[-1] < 1e-3 or K["jacobian_rank"] != 12 or len(K["jacobian_singular_values"]) != 12:
        failures.append("the Jacobian at the solution is not of rank 12 by steps that agree")
    if len(K["jacobian_singular_values"]) != 12 or max(abs(np.array(K["jacobian_singular_values"]) / sv - 1)) > 1e-3 \
            or abs(K["jacobian_condition_number"] / (sv[0] / sv[-1]) - 1) > 1e-3:
        failures.append("the stated singular values and condition number are not the system's Jacobian's")
    if not all(e < 0 for e in X["reduced_hessian"]) or X["null_residual"] > mp.mpf("1e-25"):
        failures.append("the Lagrangian's Hessian is not negative definite on the plane the gradients leave free")
    return failures


def x6_copies(w):
    failures = []
    S = w["summary"]
    if w["candidate"] != S["round13_candidate"] or w["minima"] != S["full_phase_audit"]:
        failures.append("the candidate or minima file is not the summary's")
    G = [[float(r[k]) for k in ("d_w", "d_beta", "d_delta", "d_center", "d_epsilon")] for r in w["gradients"]]
    if G != S["smooth_event_kkt"]["gradient_matrix"]:
        failures.append("the gradient file is not the summary's matrix")
    B, R = w["basin"], S["basin_restarts"]
    sc = [float(r["scale"]) for r in B]
    if [int(r["trial"]) for r in B] != list(range(R["trials"])) or R["successful"] != len(B) \
            or (min(sc), max(sc)) != (R["minimum_scale"], R["maximum_scale"]) \
            or max(float(r["parameter_distance"]) for r in B) != R["maximum_parameter_distance"]:
        failures.append("the restarts file is not the summary's count, range and largest distance")
    M = w["manifest"]
    if (M["common_scale"], M["improvement_over_round8"], M["jacobian_rank"], M["status"], M["version"]) != (
            S["round13_candidate"]["common_scale"], S["round13_candidate"]["improvement_over_round8"],
            S["smooth_event_kkt"]["jacobian_rank"], S["status"], "0.13") \
            or sorted(M["files"]) != sorted(set(w["members"]) - {"MANIFEST.json"}):
        failures.append("the manifest is not the package's")
    return failures


def wing(w):
    P = w["summary"]["round13_candidate"]["parameters"]
    return Wing(P["wing_length"], P["beta_radians"], P["delta_radians"], P["transition_center"], P["epsilon"])


def x7_curvature(w):
    failures = []
    L = w["summary"]["curvature_ledger"]
    mine = ledger(wing(w))
    if L["total_turn_radians"] != w["summary"]["round13_candidate"]["parameters"]["delta_radians"]:
        failures.append("the total turn is not delta")
    if abs(L["peak_curvature"] / mine["peak_abs_curvature"] - 1) > 1e-14:
        failures.append("the peak curvature is not the layer's")
    if abs(L["transition_width_5_95"] - mine["effective_width_5_95"]) > 2e-11:
        failures.append("the 5-95 % width is not the layer's")
    return failures


def x8_phase(w):
    failures = []
    F = w["summary"]["full_phase_audit"]
    M = phase_map(w)["minima"]
    B = branches(w)
    if F["local_minimum_count"] != len(M) or len(F["controlling_minima"]) + len(F["other_minima"]) != len(M):
        failures.append("the count of local minima is not the phase map's")
    near = lambda ph: min(M, key=lambda m: abs(m["phase"] - mp.mpf(repr(ph))))
    with mp.workdps(60):
        for c, n in zip(F["controlling_minima"], BR):
            r = BR.index(n)
            m = near(c["phase"])
            if c["name"] != n or abs(m["phase"] - B["phases"][r]) > mp.mpf("1e-12") or c["scale"] != w["summary"]["round13_candidate"]["branch_values"][n] \
                    or c["phase"] != w["summary"]["round13_candidate"]["branch_phases"][n]:
                failures.append(f"the controlling minimum {n} is not its branch")
        for o in F["other_minima"]:
            m = near(o["phase"])
            if abs(m["phase"] - mp.mpf(repr(o["phase"]))) > mp.mpf("5e-6") or abs(m["scale"] - mp.mpf(repr(o["scale"]))) > mp.mpf("1e-9"):
                failures.append(f"the other minimum at {o['phase']} is not the exact minimum near it")
            if o["source"] != m["source"]:
                failures.append(f"the other minimum at {o['phase']} is not a {o['source']} minimum")
        others = [m for m in M if all(abs(m["phase"] - B["phases"][r]) > mp.mpf("1e-9") for r in range(4))]
        if not F["no_new_lower_branch_found"] or min(m["scale"] for m in others) <= max(B["values"]):
            failures.append("a minimum outside the four is not above them")
    return failures


def projection(G):
    """the norm of the common gradient's projection onto the plane on which the four branch gradients agree"""
    with mp.workdps(60):
        Gm = [[mp.mpf(repr(float(z))) if not isinstance(z, mp.mpf) else z for z in row] for row in G]
        Dm = mp.matrix([[Gm[r][i] - Gm[0][i] for i in range(5)] for r in (1, 2, 3)])
        _, _, V = mp.svd_r(Dm, full_matrices=True)
        return mp.sqrt(mp.fsum(mp.fsum(Gm[0][i] * V[t, i] for i in range(5)) ** 2 for t in (3, 4)))


def digits_report(w):
    """how many decimals of the report's pressures, gradients and width are the quantities'"""
    rp = w["report"].replace("\r\n", "\n")
    X = kkt(w)
    G = exact_gradients(w)
    out = {}
    pr = re.findall(r"^\|\s*(?:B[12] low|120°|270°)\s*\|\s*(0\.\d{10})\s*\|$", rp, re.M)
    out["pressures"] = [(p, len(p.split(".")[1]), r10.digits_held(p, X["x"][8 + r])) for r, p in enumerate(pr)]
    rows = re.findall(r"^\|\s*(?:B[12] low|120°|270°)\s*\|((?:\s*-?\d\.\d{9}\s*\|){5})$", rp, re.M)
    out["gradients"] = []
    for r, row in enumerate(rows):
        for i, txt in enumerate(re.findall(r"-?\d\.\d{9}", row)):
            out["gradients"].append((txt, 9, r10.digits_held(txt.lstrip("-"), abs(G[r][i]))))
    m = re.search(r"W_\{5\\%-95\\%\}\\approx(0\.\d+)", gate.squash(rp))
    out["width"] = [] if not m else [(m.group(1), len(m.group(1).split(".")[1]),
                                     r10.digits_held(m.group(1), mp.mpf(repr(ledger(wing(w))["effective_width_5_95"]))))]
    return out


def x9_documents(w):
    failures = []
    S = w["summary"]
    R = S["round13_candidate"]
    P = R["parameters"]
    K = S["smooth_event_kkt"]
    rp = w["report"].replace("\r\n", "\n")
    rs = gate.squash(rp)
    one = lambda pat: (re.search(pat, rs) or [None, None])[1]
    rt = lambda txt, x: txt is not None and r10.rounds_to(txt, x)
    # section 2: the parameters, the angles in degrees, the centre length
    for pat, x in ((r"\$\$w=([\d.]+),\$\$", P["wing_length"]), (r"\$\$\\beta=(1\.\d+),\$\$", P["beta_radians"]),
                   (r"\$\$\\delta=(0\.\d+),\$\$", P["delta_radians"]), (r"\$\$c=([\d.]+),\$\$", P["transition_center"]),
                   (r"\$\$\\varepsilon=([\d.]+)\.\$\$", P["epsilon"]), (r"\\beta=([\d.]+)\^\\circ", P["beta_degrees"]),
                   (r"\\delta=([\d.]+)\^\\circ", P["delta_degrees"]), (r"\\alpha=([\d.]+)\^\\circ", P["alpha_degrees"]),
                   (r"l_0=([\d.]+)\.", P["center_length"])):
        if not rt(one(pat), x):
            failures.append(f"section 2's {pat[2:14]} is not the data's")
    # section 3's table, the common scale, the spread, the gain, the distance
    rows = re.findall(r"^\|\s*(B1 low|B2 low|120°|270°)\s*\|\s*([\d.]+)\s*\|\s*([\d.]+)\s*\|$", rp, re.M)
    if len(rows) != 4 or any(not rt(ph, R["branch_phases"][n]) or not rt(v, R["branch_values"][n]) for (_, ph, v), n in zip(rows, BR)):
        failures.append("section 3's table is not the data's")
    for pat, x in ((r"\\boxed\{s_\{13\}=([\d.]+)\}", R["common_scale"]), (r"分支差距：\$\$([\d.e-]+)\.\$\$", R["branch_spread"]),
                   (r"s_\{13\}-s_8=([\d.e-]+)\}", R["improvement_over_round8"]), (r"1-s_\{13\}=([\d.e-]+)\.", R["distance_to_one"])):
        if not rt(one(pat), x):
            failures.append(f"section 3's {pat[:16]} is not the data's")
    if len(re.findall(r"\\boxed\{s_\{13\}=0\.998914343297485\}", rs)) != 2:
        failures.append("the common scale is not stated alike in sections 3 and 15")
    # section 4's shifts and the epsilon line
    for k, sym in zip(KEYS, ("w", "\\beta", "\\delta", "c", "\\varepsilon")):
        if not rt(one(r"\\Delta" + re.escape(sym) + r"=(-?\d\.\d+e[+-]\d+)[,.]\$\$"), R["parameter_shift"][k]):
            failures.append(f"section 4's shift in {k} is not the data's")
    if not rt(one(r"0\.037\\longrightarrow([\d.]+)\."), P["epsilon"]):
        failures.append("section 4's epsilon is not the data's")
    # section 5's pressures, section 6's gradients, residual and norm
    pr = re.findall(r"^\|\s*(?:B[12] low|120°|270°)\s*\|\s*(0\.\d+)\s*\|$", rp, re.M)
    if len(pr) != 4 or any(float(p) != m for p, m in zip(pr, K["branch_pressures"])):
        failures.append("section 5's pressures are not the data's")
    mu = K["branch_pressures"]
    if f"承擔約 ${100 * mu[2]:.1f}\\%$ 的壓力" not in rp or f"分支約 ${100 * mu[3]:.1f}\\%$" not in rp:
        failures.append("section 5's shares are not the data's")
    rows = re.findall(r"^\|\s*(?:B[12] low|120°|270°)\s*\|((?:\s*-?\d\.\d{9}\s*\|){5})$", rp, re.M)
    if len(rows) != 4 or any(not rt(t, x) for row, g in zip(rows, K["gradient_matrix"]) for t, x in zip(re.findall(r"-?\d\.\d{9}", row), g)):
        failures.append("section 6's gradients are not the data's")
    if str(K["weighted_stationarity_residual"]).replace(" ", "") not in rs or not rt(one(r"範數：\$\$([\d.e-]+)\.\$\$"), K["weighted_stationarity_norm"]):
        failures.append("section 6's residual or norm is not the data's")
    # section 8, 9, 10, 11, 12
    if rs.count(f"\\operatorname{{rank}}J={K['jacobian_rank']}.") != 2 \
            or not rt(one(r"最小奇異值：\$\$([\d.e+-]+)\.\$\$"), K["jacobian_singular_values"][-1]) \
            or not rt(one(r"\\kappa\(J\)\\approx([\d.e+-]+)\.\$\$"), K["jacobian_condition_number"]):
        failures.append("section 8 is not the data's")
    B = S["basin_restarts"]
    m = re.search(r"\[([\d.]+),\\([\d.]+)\]", rs)
    if f"\\boxed{{{B['successful']}/{B['trials']}}}" not in rs or not m or not rt(m.group(1), B["minimum_scale"]) or not rt(m.group(2), B["maximum_scale"]) \
            or not rt(one(r"最大參數距離：\$\$([\d.e-]+)\.\$\$"), B["maximum_parameter_distance"]):
        failures.append("section 9 is not the data's")
    L = S["curvature_ledger"]
    if not rt(one(r"\\int\|\\kappa\|\\,ds\\approx([\d.]+)\."), L["total_turn_radians"]) or not rt(one(r"\\infty\\approx([\d.]+)\."), L["peak_curvature"]) \
            or not rt(one(r"W_\{5\\%-95\\%\}\\approx([\d.]+)\."), L["transition_width_5_95"]):
        failures.append("section 10 is not the data's")
    others = S["full_phase_audit"]["other_minima"]
    if f"完整相位景觀仍有 ${S['full_phase_audit']['local_minimum_count']}$ 個局部極小" not in rp \
            or not rt(one(r"最低的其他局部極小約為：?\$\$([\d.]+),\$\$"), min(o["scale"] for o in others)):
        failures.append("section 11 is not the data's")
    if not rt(one(r"其尺度約為：\$\$\\boxed\{([\d.]+)\.\}"), R["common_scale"]):
        failures.append("section 12's scale is not the data's")
    # section 7: the projection of the common gradient onto the plane the branches leave
    if not rt(one(r"投影約為：\$\$([\d.]+)\\times10\^\{-9\}"), float(projection(K["gradient_matrix"]) * 10 ** 9)):
        failures.append("section 7's projection is not the stated gradients'")
    # the README, the audit note, the handoff
    rd = w["readme"]
    for label, x in (("第 8 輪候選：", S["round8_reference"]["scale"]), ("第 13 輪候選：", R["common_scale"]),
                     ("提升：", R["improvement_over_round8"]), ("四分支差距：", R["branch_spread"])):
        m = re.search(re.escape(label) + r"([\d.e-]+)", rd)
        if not m or not rt(m.group(1), x):
            failures.append(f"the README's {label} is not the data's")
    if f"KKT Jacobian：{K['jacobian_rank']}/12 滿秩" not in rd or f"局部重啟：{B['successful']}/{B['trials']} 返回同一盆地" not in rd:
        failures.append("the README's rank or restarts are not the data's")
    for k, sym in zip(KEYS, ("w", "\\beta", "\\delta", "c", "\\varepsilon")):
        m = re.search(r"\$" + re.escape(sym) + r"=([\d.]+)\$", rd)
        if not m or not rt(m.group(1), P[k]):
            failures.append(f"the README's {k} is not the data's")
    an = w["audit_note"]
    m = re.search(r"```json\n(.*?)\n```", an, re.S)
    if not m or json.loads(m.group(1)) != R:
        failures.append("the audit note's candidate is not the data's")
    if f"[{', '.join(repr(x) for x in K['branch_pressures'])}]" not in an or str(K["weighted_stationarity_residual"]) not in an \
            or f"rank = {K['jacobian_rank']}" not in an or str(K["jacobian_singular_values"]) not in an \
            or not re.search(r"condition = ([\d.e+]+)", an) or float(re.search(r"condition = ([\d.e+]+)", an).group(1)) != K["jacobian_condition_number"] \
            or float(re.search(r"norm = ([\d.e-]+)", an).group(1)) != K["weighted_stationarity_norm"]:
        failures.append("the audit note's pressures, residual or Jacobian are not the data's")
    m = re.search(r"scale range = \[([\d.]+), ([\d.]+)\]", an)
    if f"successful = {B['successful']} / {B['trials']}" not in an or not m or (float(m.group(1)), float(m.group(2))) != (B["minimum_scale"], B["maximum_scale"]):
        failures.append("the audit note's restarts are not the data's")
    if not rt((re.search(r"s_\{\\mathrm\{base\}\}=([\d.]+)\.\$\$", gate.squash(w["handoff"])) or [None, None])[1], R["common_scale"]) \
            or not rt(one(r"s_\{\\mathrm\{base\}\}=([\d.]+)\$\$"), R["common_scale"]):
        failures.append("the handoff's base scale is not the data's")
    # the printed digits against the quantities
    D = digits_report(w)
    if any(held < n for _, n, held in D["pressures"] + D["gradients"] + D["width"]):
        failures.append("the report prints the pressures, gradients and width past the decimals the quantities hold")
    return failures


CHECKS = {"X1": x1_parameters, "X2": x2_branches, "X3": x3_gradients, "X4": x4_kkt, "X5": x5_jacobian, "X6": x6_copies,
          "X7": x7_curvature, "X8": x8_phase, "X9": x9_documents}
CURVE = {"summary", "round05", "legs"}
DEPENDS = {
    "X1": {"summary", "round09", "round05", "legs"},
    "X2": CURVE | {"round09"},
    "X3": CURVE,
    "X4": CURVE,
    "X5": CURVE,
    "X6": {"summary", "candidate", "minima", "gradients", "basin", "manifest", "members"},
    "X7": {"summary"},
    "X8": CURVE,
    "X9": CURVE | {"report", "readme", "audit_note", "handoff"},
}


# ---------------------------------------------------------------- the certificate: the least is at least T

def least_bounds(H, T, K=262144, eta=0.05):
    """s(phi) >= T for every phi of the curve H, in interval arithmetic - RUN-011's certificate with a target T in
    place of s(3pi/2), and the 2pi/3 corners' value compared with T:
    - within eta of 3pi/2 and of 2pi/3: s >= g, the sum of the three functionals on fixed corners, and g' has one sign
      on each side, so g moves away from its value at the cusp: at 3pi/2 that value is s(3pi/2) (the support identity,
      with theta in (pi/6, pi/2)), at 2pi/3 a linear form in x1, y1 with each box entering once; both must exceed T
    - elsewhere: on an interval of half-width h about m, s >= G(m) - |G'(m)| h - Lip h^2/2 for the points chosen at m
      (a float estimate, lowered by the boxes' reach, picks the intervals; interval arithmetic decides each, and an
      interval it refuses is split)
    and s(3pi/2) <= the upper end of its enclosure, an upper bound for the least"""
    iv.dps = 30
    A, B = r10.legs_iv(H)
    funs = [(1 / A, 1 / B), (-1 / A, iv.mpf(0)), (iv.mpf(0), -1 / B)]           # max_d, min_x, min_y
    bx = r10.certificate_boxes(H, K)
    PX, PY = bx["PX"], bx["PY"]
    n_pts = len(PX)
    x1, y1, l0 = bx["x1"], bx["y1"], bx["l0"]
    D = H["summary"]["definitions"]["candidate"]
    beta, delta = iv.mpf(D["beta"]), iv.mpf(D["delta"])
    premise = bool(lo(beta) > hi(iv.pi / 6) and hi(beta + delta) < lo(iv.pi / 2))
    S270 = cusp_270(x1, y1, l0, A, B)
    with mp.workdps(40):
        Tm = mp.mpf(T)
    K1 = (len(PX) - 2) // 2
    p = {"p0": (PX[0], PY[0]), "p1": (PX[K1], PY[K1]), "p2": (PX[K1 + 1], PY[K1 + 1]), "p3": (PX[-1], PY[-1])}

    def G(ph, pts_, d=0):
        c, s = iv.cos(ph), iv.sin(ph)
        tot = iv.mpf(0)
        for (a0, a1), (px, py) in zip(funs, pts_):
            if d == 0:
                tot += (a0 * c + a1 * s) * px + (-a0 * s + a1 * c) * py
            else:
                tot += (-a0 * s + a1 * c) * px + (-a0 * c - a1 * s) * py
        return tot

    def one_sign(pts_, a_, b_, sign, pieces=64):
        xs = [a_] + [a_ + (b_ - a_) * j / pieces for j in range(1, pieces)] + [b_]
        for j in range(pieces):
            g = G(iv.mpf([xs[j], xs[j + 1]]), pts_, 1)
            if not ((sign < 0 and hi(g) < 0) or (sign > 0 and lo(g) > 0)):
                return False
        return True
    with mp.workdps(40):
        e270, e120 = mp.mpf(3) * mp.pi / 2, 2 * mp.pi / 3
        E270, E120 = 3 * iv.pi / 2, 2 * iv.pi / 3
        wide = eta + mp.mpf("1e-12")
        signs = {"left_of_3pi_2": one_sign((p["p1"], p["p3"], p["p3"]), lo(E270) - wide, hi(E270), -1),
                 "right_of_3pi_2": one_sign((p["p1"], p["p0"], p["p3"]), lo(E270), hi(E270) + wide, +1),
                 "left_of_2pi_3": one_sign((p["p3"], p["p2"], p["p1"]), lo(E120) - wide, hi(E120), -1),
                 "right_of_2pi_3": one_sign((p["p0"], p["p2"], p["p1"]), lo(E120), hi(E120) + wide, +1)}
    ph = 2 * iv.pi / 3
    c, s = iv.cos(ph), iv.sin(ph)

    def g120(md):
        """g at 2pi/3 with max_d on md (p3 or p0), min_x on p2 and min_y on p1, uncentred (the functionals sum to zero)"""
        nd = (c / A + s / B, -s / A + c / B)
        nx = (-c / A, s / A)
        ny = (-s / B, -c / B)
        return (nx[0] + ny[0] + (2 * nd[0] if md == "p3" else 0)) * x1 + (nx[1] + ny[1]) * y1 + nx[0] * l0 + (nd[0] * l0 if md == "p3" else 0)
    G120 = [g120("p3"), g120("p0")]
    cx = np.array([float(v.mid) for v in PX])
    cy = np.array([float(v.mid) for v in PY])
    rx = np.array([float(v.delta) / 2 for v in PX])
    ry = np.array([float(v.delta) / 2 for v in PY])
    with mp.workdps(40):
        mxe = max(max(abs(mp.mpf(v._mpi_[0])), abs(mp.mpf(v._mpi_[1]))) for v in PX)
        mye = max(max(abs(mp.mpf(v._mpi_[0])), abs(mp.mpf(v._mpi_[1]))) for v in PY)
    R_up = hi(iv.sqrt(iv.mpf(mxe) ** 2 + iv.mpf(mye) ** 2))
    Lip_up = hi(iv.mpf(R_up) * sum((iv.mpf(hi(iv.sqrt(a0 ** 2 + a1 ** 2))) for a0, a1 in funs), iv.mpf(0)))
    Lip = float(Lip_up) * (1 + 1e-12)
    af = [(float(a0.mid), float(a1.mid)) for a0, a1 in funs]
    neg_tang = -np.unwrap(np.arctan2(np.diff(cy), np.diff(cx)))
    T_f = float(Tm)

    def best_idx(ph_):
        c_, s_ = math.cos(ph_), math.sin(ph_)
        out = []
        for a0, a1 in af:
            n0, n1 = a0 * c_ + a1 * s_, -a0 * s_ + a1 * c_
            psi = math.atan2(n1, n0)
            ci = [0, n_pts - 1]
            for bnd in (psi - math.pi / 2, psi + math.pi / 2):
                k = int(np.searchsorted(neg_tang, -math.remainder(bnd, 2 * math.pi)))
                ci += [i for i in range(k - 2, k + 3) if 0 <= i < n_pts]
            ci = np.array(ci)
            vals = cx[ci] * n0 + cy[ci] * n1
            i = int(ci[int(np.argmax(vals))])
            out.append((i, float(vals.max()) - abs(n0) * rx[i] - abs(n1) * ry[i],
                        float(cx[i] * (-a0 * s_ + a1 * c_) + cy[i] * (-a0 * c_ - a1 * s_))))
        return out

    def iv_bound(u, v, m, idx):
        M = iv.mpf(m)
        h = iv.mpf(max(hi(M - iv.mpf(u)), hi(iv.mpf(v) - M)))
        g0, g1 = G(M, [(PX[i], PY[i]) for i in idx]), G(M, [(PX[i], PY[i]) for i in idx], 1)
        return lo(g0 - iv.mpf(max(abs(lo(g1)), abs(hi(g1)))) * h - iv.mpf(Lip_up) * h * h / 2)
    pieces = [(0.0, float(e120) - eta), (float(e120) + eta, float(e270) - eta), (float(e270) + eta, 2 * math.pi + 1e-9)]
    leaves, refused, worst = [], 0, None
    for a_, b_ in pieces:
        todo = [(a_ + (b_ - a_) * k / 512, a_ + (b_ - a_) * (k + 1) / 512) for k in range(512)][::-1]
        while todo:
            u, v = todo.pop()
            m = (u + v) / 2
            h = (v - u) / 2
            bi = best_idx(m)
            if sum(x[1] for x in bi) - abs(sum(x[2] for x in bi)) * h - Lip * h * h / 2 - 1e-12 > T_f:
                lb = iv_bound(u, v, m, [x[0] for x in bi])
                if lb > Tm:
                    leaves.append((u, v))
                    if worst is None or lb < worst[0]:
                        worst = (lb, m)
                    continue
                refused += 1
            if h < 1e-13:
                raise ValueError(f"cannot cover near {m}")
            todo += [(m, v), (u, m)]
    sl = sorted(leaves)
    tiled = all(x[1] == y[0] or (x[1], y[0]) in ((pieces[0][1], pieces[1][0]), (pieces[1][1], pieces[2][0]))
                for x, y in zip(sl, sl[1:])) and sl[0][0] == 0.0 and sl[-1][1] >= 2 * math.pi
    with mp.workdps(40):
        nb_ok = abs(mp.mpf(pieces[0][1]) - (e120 - eta)) < 1e-12 and abs(mp.mpf(pieces[1][0]) - (e120 + eta)) < 1e-12 \
            and abs(mp.mpf(pieces[1][1]) - (e270 - eta)) < 1e-12 and abs(mp.mpf(pieces[2][0]) - (e270 + eta)) < 1e-12
    out = {"target": dec(Tm, 20), "wing_steps": K, "interval_digits": 30, "premise_beta_above_pi_6_alpha_below_pi_2": premise,
           "largest_box": [max(float(v.delta) for v in PX), max(float(v.delta) for v in PY)],
           "s270_at_least": dec(lo(S270), 18), "s270_at_most": dec(hi(S270), 18, up=True),
           "neighbourhood_half_width": eta, "derivative_signs": signs,
           "s_2pi_3_corners_at_least": [dec(lo(g), 18) for g in G120], "lipschitz_at_most": dec(Lip_up, 8, up=True),
           "phase_intervals": len(leaves), "intervals_refused_and_split": refused, "pieces_tiled": bool(tiled and nb_ok),
           "elsewhere_at_least": dec(worst[0], 18), "elsewhere_worst_midpoint": worst[1]}
    out["certified"] = bool(premise and all(signs.values()) and tiled and nb_ok and lo(S270) > Tm
                            and all(lo(g) > Tm for g in G120) and worst[0] > Tm)
    out["_lower"], out["_upper"] = Tm, hi(S270)
    return out


def s270_enclosure(H, K=262144):
    """s(3pi/2) of a curve, enclosed from the same boxes (RUN-011's enclosure of Round 08's least)"""
    iv.dps = 30
    A, B = r10.legs_iv(H)
    bx = r10.certificate_boxes(H, K)
    S = cusp_270(bx["x1"], bx["y1"], bx["l0"], A, B)
    return lo(S), hi(S)


# ---------------------------------------------------------------- findings

def findings(w):
    out = {}
    S = w["summary"]
    R = S["round13_candidate"]
    B = branches(w)
    s08 = r10.s270(hp08(w))
    with mp.workdps(60):
        out["branch_values_60_digits"] = {n: mp.nstr(v, 20) for n, v in zip(BR, B["values"])}
        out["branch_values_stated_minus_exact"] = {n: mp.nstr(mp.mpf(repr(R["branch_values"][n])) - v, 3) for n, v in zip(BR, B["values"])}
        out["low_phases_stated_minus_exact"] = [mp.nstr(mp.mpf(repr(R["branch_phases"][n])) - B["phases"][j], 3) for j, n in enumerate(BR[:2])]
        out["spread_exact"] = mp.nstr(max(B["values"]) - min(B["values"]), 6)
        out["improvement_exact"] = mp.nstr(min(B["values"]) - s08, 12)
        G = exact_gradients(w)
        Gs = S["smooth_event_kkt"]["gradient_matrix"]
        out["gradients_stated_minus_exact_largest"] = mp.nstr(max(abs(mp.mpf(repr(Gs[r][i])) - G[r][i]) for r in range(4) for i in range(5)), 3)
        mu = S["smooth_event_kkt"]["branch_pressures"]
        out["residual_of_stated_pressures_with_exact_gradients"] = [mp.nstr(mp.fsum(mp.mpf(repr(mu[r])) * G[r][i] for r in range(4)), 3) for i in range(5)]
        X = kkt(w)
        x = X["x"]
        names = list(DKEYS) + ["phi1", "phi2", "s", "mu1", "mu2", "mu3", "mu4"]
        st = [mp.mpf(v) for v in stated(w)] + [mp.mpf(repr(R["branch_phases"][n])) for n in BR[:2]] + [mp.mpf(repr(R["common_scale"]))] \
            + [mp.mpf(repr(m)) for m in mu]
        out["kkt_solution"] = {n: mp.nstr(v, 25) for n, v in zip(names, x)}
        out["kkt_stated_minus_solution"] = {n: mp.nstr(a - v, 3) for n, a, v in zip(names, st, x)}
        out["kkt_chord_steps"], out["kkt_residual"] = X["steps"], mp.nstr(X["residual"], 3)
        out["jacobian_singular_values"] = [float(f"{v:.6g}") for v in X["sv"]]
        out["jacobian_condition"] = float(f"{X['sv'][0] / X['sv'][-1]:.6g}")
        out["jacobian_steps_relative_difference"] = float(f"{max(abs(X['sv'] - X['sv_coarse']) / X['sv']):.3g}")
        out["gradient_singular_values_at_solution"] = [mp.nstr(v, 5) for v in X["gradient_sv"]]
        out["reduced_hessian_eigenvalues"] = [mp.nstr(v, 6) for v in X["reduced_hessian"]]
        out["low_branch_phase_curvature"] = [mp.nstr(v, 6) for v in X["phase_curvature"]]
        out["projection_stated_gradients"] = mp.nstr(projection(Gs), 4)
        out["projection_exact_gradients_at_stated_point"] = mp.nstr(projection(G), 4)
    out["printed_digits"] = digits_report(w)
    sc = [float(r["scale"]) for r in w["basin"]]
    out["restart_scales_below_candidate"] = [float(f"{R['common_scale'] - max(sc):.4g}"), float(f"{R['common_scale'] - min(sc):.4g}")]
    W = wing(w)
    L = ledger(W)
    out["curvature_exact"] = {"peak": L["peak_abs_curvature"], "width": L["effective_width_5_95"],
                              "stated_width_minus_exact": S["curvature_ledger"]["transition_width_5_95"] - L["effective_width_5_95"]}
    M = phase_map(w)
    C = r10.cand(hp(w), 60)
    with mp.workdps(60):
        out["switches"] = len(M["switches"])
        out["corner_switches"] = [mp.nstr(s["phase"], 12) for s in M["switches"] if s["corner"]]
        out["stationary_points"] = len(M["stationary"])
        out["local_minima"] = [{"source": m["source"], "phase": mp.nstr(m["phase"], 12), "scale": mp.nstr(m["scale"], 13),
                                **({"to_next_switch": mp.nstr(m["to_next_switch"], 3)} if "to_next_switch" in m else {})} for m in M["minima"]]
        out["stated_other_minima_against_s"] = [{"phase": o["phase"], "stated": o["scale"], "s_at_phase": mp.nstr(C.s(mp.mpf(repr(o["phase"]))), 13)}
                                               for o in S["full_phase_audit"]["other_minima"]]
    # the certificate: the stated curve's least is at least the common scale less 1.5e-10; Round 08's least is s(3pi/2)
    # (RUN-011), enclosed from the same boxes
    with mp.workdps(40):
        T = mp.mpf(repr(R["common_scale"])) - mp.mpf("1.5e-10")
    cert = least_bounds(hp(w), T)
    lo8, hi8 = s270_enclosure(hp08(w))
    with mp.workdps(40):
        cert["round08_least_within"] = [dec(lo8, 18), dec(hi8, 18, up=True)]
        cert["gain_over_round08_at_least"] = dec(cert["_lower"] - hi8, 18)
        cert["gain_over_round08_at_most"] = dec(cert["_upper"] - lo8, 18, up=True)
    del cert["_lower"], cert["_upper"]
    out["certificate"] = cert
    return out


# ---------------------------------------------------------------- the drill

def drill_cases(base):
    def edit(key, fn):
        v = copy.deepcopy(base[key])
        fn(v)
        return v
    rp = base["report"]

    def rrep(text, old, new, n=1):
        assert text.count(old) == n, old
        return text.replace(old, new)

    def kkt_plant(fn):
        """a plant on the gradients or pressures, with the residual and its norm restated from them - in the summary,
        the report and the audit note - so that only the comparison with the system can see it"""
        S = copy.deepcopy(base["summary"])
        K = S["smooth_event_kkt"]
        old_res, old_norm = K["weighted_stationarity_residual"], K["weighted_stationarity_norm"]
        fn(K)
        mu, G = K["branch_pressures"], K["gradient_matrix"]
        K["weighted_stationarity_residual"] = [sum(mu[r] * G[r][i] for r in range(4)) for i in range(5)]
        K["weighted_stationarity_norm"] = math.sqrt(sum(x * x for x in K["weighted_stationarity_residual"]))
        new_res, new_norm = K["weighted_stationarity_residual"], K["weighted_stationarity_norm"]
        report = rrep(rrep(rp, str(old_res), str(new_res)), f"{old_norm:.12e}", f"{new_norm:.12e}")
        audit = rrep(rrep(base["audit_note"], str(old_res), str(new_res)), f"norm = {old_norm:.16e}", f"norm = {new_norm:.16e}")
        return S, report, audit

    def grad_plant(r, i, old, new):
        S, report, audit = kkt_plant(lambda K: K["gradient_matrix"][r].__setitem__(i, float(new)))
        assert base["summary"]["smooth_event_kkt"]["gradient_matrix"][r][i] == float(old)
        rows = copy.deepcopy(base["gradients"])
        key = ("d_w", "d_beta", "d_delta", "d_center", "d_epsilon")[i]
        rows[r][key] = f"{float(new):.18e}"
        report = rrep(report, f" {float(old):.9f} |", f" {float(new):.9f} |")
        return {"summary": S, "gradients": rows, "report": report, "audit_note": audit}

    # the branch value B1, everywhere it is stated
    b1 = edit("summary", lambda S: (S["round13_candidate"]["branch_values"].__setitem__("B1_low", 0.9989143432974897),
                                    S["full_phase_audit"]["controlling_minima"][0].__setitem__("scale", 0.9989143432974897)))
    b1c = edit("candidate", lambda C: C["branch_values"].__setitem__("B1_low", 0.9989143432974897))
    b1m = edit("minima", lambda M: M["controlling_minima"][0].__setitem__("scale", 0.9989143432974897))
    b1r = rrep(rp, "| B1 low | 0.154865921088826 | 0.998914343297489 |", "| B1 low | 0.154865921088826 | 0.998914343297490 |")
    b1a = rrep(base["audit_note"], '"B1_low": 0.9989143432974892', '"B1_low": 0.9989143432974897')
    # the distance to one
    d1 = edit("summary", lambda S: S["round13_candidate"].__setitem__("distance_to_one", 0.0010856567025149))
    d1c = edit("candidate", lambda C: C.__setitem__("distance_to_one", 0.0010856567025149))
    d1a = rrep(base["audit_note"], '"distance_to_one": 0.0010856567025149255', '"distance_to_one": 0.0010856567025149')
    # the centre length
    cl = edit("summary", lambda S: S["round13_candidate"]["parameters"].__setitem__("center_length", 0.3277914308833731))
    clc = edit("candidate", lambda C: C["parameters"].__setitem__("center_length", 0.3277914308833731))
    cla = rrep(base["audit_note"], '"center_length": 0.3277914308833727', '"center_length": 0.3277914308833731')
    # the reference's epsilon
    ref = edit("summary", lambda S: S["round8_reference"]["parameters"].__setitem__("epsilon", 0.0371))
    # two gradients and a pressure, each with the residual restated
    g1 = grad_plant(3, 2, "-0.1629602289", "-0.1629612289")
    g2 = grad_plant(1, 4, "-0.0077729156", "-0.0077739156")
    pS, pR, pA = kkt_plant(lambda K: K["branch_pressures"].__setitem__(0, 0.02826233))
    pR = rrep(pR, "| B1 low | 0.0282523300 |", "| B1 low | 0.0282623300 |")
    pA = rrep(pA, "[0.02825233, 0.02837576", "[0.02826233, 0.02837576")
    unk = edit("summary", lambda S: S["smooth_event_kkt"].__setitem__("unknown_count", 13))
    # the rank, everywhere it is stated
    rk = edit("summary", lambda S: S["smooth_event_kkt"].__setitem__("jacobian_rank", 11))
    rkm = edit("manifest", lambda M: M.__setitem__("jacobian_rank", 11))
    rkr = rrep(rp, "\\operatorname{rank}J=12.", "\\operatorname{rank}J=11.", 2)
    rkd = rrep(base["readme"], "KKT Jacobian：12/12 滿秩", "KKT Jacobian：11/12 滿秩")
    rka = rrep(base["audit_note"], "rank = 12", "rank = 11")
    # eleven singular values
    sv = edit("summary", lambda S: S["smooth_event_kkt"]["jacobian_singular_values"].pop(0))
    sva = rrep(base["audit_note"], "[429267986.0, 135912243.0,", "[135912243.0,")
    basin = edit("basin", lambda B: B[5].__setitem__("scale", "0.99891434329721174"))
    files = edit("manifest", lambda M: M["files"].remove("figures/05_basin_restarts.png"))
    peak = edit("summary", lambda S: S["curvature_ledger"].__setitem__("peak_curvature", 2.1011016211131104))
    width = edit("summary", lambda S: S["curvature_ledger"].__setitem__("transition_width_5_95", 0.03644180145175618))
    width_r = rrep(rp, "0.036441801352", "0.036441801452")
    count = edit("summary", lambda S: S["full_phase_audit"].__setitem__("local_minimum_count", 9))
    count_m = edit("minima", lambda M: M.__setitem__("local_minimum_count", 9))
    count_r = rrep(rp, "完整相位景觀仍有 $8$ 個局部極小", "完整相位景觀仍有 $9$ 個局部極小")
    src = edit("summary", lambda S: S["full_phase_audit"]["other_minima"][3].__setitem__("source", "smooth"))
    src_m = edit("minima", lambda M: M["other_minima"][3].__setitem__("source", "smooth"))
    share = rrep(rp, "承擔約 $73.6\\%$ 的壓力", "承擔約 $73.7\\%$ 的壓力")
    readme = rrep(base["readme"], "提升：4.212852400265e-09", "提升：4.212852400275e-09")
    hand = rrep(base["handoff"], "s_{\\mathrm{base}}=0.99891434329749.", "s_{\\mathrm{base}}=0.99891434329759.")
    audit = rrep(base["audit_note"], "successful = 20 / 20", "successful = 19 / 20")
    return [
        ("K0-unmodified", "control", {}, set()),
        ("K1-report-crlf", "control", {"report": rp.replace("\n", "\r\n")}, set()),
        ("K2-readme-trailing-spaces", "control", {"readme": base["readme"] + "  "}, set()),
        ("K3-restarts-copied", "control", {"basin": copy.deepcopy(base["basin"])}, set()),
        ("M1-triangle-leg", "mutation", {"legs": ("1.000000005", "1")}, {"X1", "X2", "X4", "X8"}),
        ("M2-reference-epsilon", "mutation", {"summary": ref}, {"X1"}),
        ("M3-centre-length", "mutation", {"summary": cl, "candidate": clc, "audit_note": cla}, {"X1"}),
        ("M4-branch-value", "mutation", {"summary": b1, "candidate": b1c, "minima": b1m, "report": b1r, "audit_note": b1a}, {"X2"}),
        ("M5-distance-to-one", "mutation", {"summary": d1, "candidate": d1c, "audit_note": d1a}, {"X2"}),
        ("M6-gradient-270-delta", "mutation", g1, {"X3"}),
        ("M7-gradient-B2-epsilon", "mutation", g2, {"X3"}),
        ("M8-pressure", "mutation", {"summary": pS, "report": pR, "audit_note": pA}, {"X4"}),
        ("M9-unknown-count", "mutation", {"summary": unk}, {"X4"}),
        ("M10-rank", "mutation", {"summary": rk, "manifest": rkm, "report": rkr, "readme": rkd, "audit_note": rka}, {"X5"}),
        ("M11-singular-values", "mutation", {"summary": sv, "audit_note": sva}, {"X5"}),
        ("M12-gradient-file", "mutation", {"gradients": g1["gradients"]}, {"X6"}),
        ("M13-restart-scale", "mutation", {"basin": basin}, {"X6"}),
        ("M14-manifest-files", "mutation", {"manifest": files}, {"X6"}),
        ("M15-peak", "mutation", {"summary": peak}, {"X7"}),
        ("M16-width", "mutation", {"summary": width, "report": width_r}, {"X7"}),
        ("M17-minimum-count", "mutation", {"summary": count, "minima": count_m, "report": count_r}, {"X8"}),
        ("M18-minimum-source", "mutation", {"summary": src, "minima": src_m}, {"X8"}),
        ("M19-report-share", "mutation", {"report": share}, {"X9"}),
        ("M20-readme-gain", "mutation", {"readme": readme}, {"X9"}),
        ("M21-handoff-base", "mutation", {"handoff": hand}, {"X9"}),
        ("M22-audit-restarts", "mutation", {"audit_note": audit}, {"X9"}),
    ]


def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("--drill", action="store_true")
    ap.add_argument("--blind", metavar="CID")
    ap.add_argument("--log", default=None)
    ap.add_argument("--no-findings", action="store_true")
    ap.add_argument("--dry", action="store_true")
    args = ap.parse_args()
    if args.blind:
        CHECKS[args.blind] = lambda w: []
        args.drill = True
    w = build_world()
    results = {cid: fn(w) for cid, fn in CHECKS.items()}
    for cid, f in results.items():
        print(f"{cid} {'PASS' if not f else 'FAIL'}  ({len(f)} failures)")
        for x in f[:12]:
            print("    ", x)
    log = {"run": "RUN-014",
           "subject": "Round 13 (v0.13): the five-parameter smooth event-KKT system - branches, gradients, pressures, Jacobian, restarts, ledger, phase audit",
           "package": {"id": RID, "sha256": sha256(source_bytes(RID))},
           "checks": {cid: {"pass": not f, "failures": f} for cid, f in results.items()},
           "findings": None if (args.blind or args.no_findings) else findings(w)}
    if log["findings"]:
        print(json.dumps(log["findings"], indent=1, ensure_ascii=False, default=str)[:12000])
    if args.dry:
        # the dry run: every expected set empty, to read what each plant trips against the sets reasoned beforehand
        drill = gate.run_drill(build_world, CHECKS, DEPENDS, [(c, k, ch, set()) for c, k, ch, _ in drill_cases(build_world())])
        for d in drill:
            print(d["case"], "tripped", d["tripped"])
        return
    if args.drill:
        drill = gate.run_drill(build_world, CHECKS, DEPENDS, drill_cases(build_world()))
        log["drill"] = gate.summarize(drill, CHECKS)
        print(gate.report_drill(log["drill"]))
        for d in drill:
            if not d["ok"]:
                print("   DRILL FAIL", d)
    if args.blind:
        seen = sum(1 for d in log["drill"]["cases"] if args.blind in d["expected"] and args.blind not in d["tripped"])
        print(f"blind {args.blind}: reported blind in {seen} cases (log not written)")
        return
    if args.drill:
        path = pathlib.Path(args.log) if args.log else LOG
        path.parent.mkdir(parents=True, exist_ok=True)
        path.write_text(json.dumps(log, ensure_ascii=False, indent=1, default=str) + "\n", encoding="utf-8", newline="\n")


if __name__ == "__main__":
    main()