"""RUN-006 - Round 05 (v0.5): the contact-event equations and the nonsmooth KKT system of the five-link.

Round 05 writes Round 04's four branch minima (RUN-005) at exact event phases - phi1 = pi/2 - alpha, phi2 = pi/2 -
beta, phi3 = 2pi/3, phi4 = 3pi/2 - gives each branch value as a closed form m_r(p) in p = (l1, l2, beta, delta), and
solves nine equations in nine unknowns: m_r(p) = s (r = 1..4), sum mu_r = 1, sum mu_r grad m_r(p) = 0. It reports
the common event scale 0.998903757132509, the four branches level to 2.2e-16, a 9x9 Jacobian of rank 9 with
smallest singular value 3.607219048175e-02 and condition number 162.78, 40 of 40 perturbed starts returning to the
same root, one-sided derivatives D- < 0 < D+ at each event, and a 262,144-phase audit whose grid minimum is the
event scale.

Checks (each returns a list of failures; empty = pass)
  V1 setup          the triangle's A and B are the certified legs divided by 1.0048 and B/A = sqrt 3; the logged
                    curve is the event parameters (from the origin, alpha = beta + delta, l3 = 1 - 2(l1 + l2),
                    degrees the radians'), of length 1, its own mirror image traversed backwards
  V2 phases         the event phases are pi/2 - alpha, pi/2 - beta, 2pi/3 and 3pi/2, in radians and degrees
  V3 values         the report's four closed forms are the exact scale at the event phases, and both are the stated
                    branch values; the common scale and the spread are theirs
  V4 global         the least exact scale over all phases is the common scale, reached at all four events, and the
                    stated 262,144-phase audit is that grid's minimum, at its stated phase
  V5 gradients      the logged gradients are the closed forms' (complex-step), to 1e-12
  V6 kkt            mu is on the simplex and the nine equations hold at the stated solution, with the stated
                    residual's order; the improvement over Round 04 and the distance to 1 are the differences they
                    name
  V7 jacobian       the logged Jacobian is this run's, to 1e-7; its rank, singular values and condition number are
                    the stated ones
  V8 roots          the stated perturbation trials all succeed, with the stated largest distance; this run's Newton
                    iteration from 40 starts, each unknown perturbed by a Gaussian of standard deviation 1e-2, returns
                    to the stated root, to 1e-12
  V9 cusps          at each event the contact ledger is the recorded one, and the exact one-sided derivatives (the
                    active pieces' derivatives) are the stated ones, to 1e-6, with D- < 0 < D+
  V10 consistency   the README and the report state the data's numbers and tables, and event_system.py states the
                    report's nine equations and unknowns

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/run006_round05_contact_events.py [--drill [--log PATH] | --blind V#]
"""
import argparse
import copy
import csv
import io
import json
import math
import pathlib
import re

import numpy as np

from corpus import TREE, sha256, source_bytes
import gate
from run001_corpus_identity import restore_escapes
from run003_round02_phase_scan_chirality import exact_min, exact_scale, legs
from run004_round03_symmetric_links import branch_minima, contacts, grid_values, load, same_B
from run005_round04_branch_equalization import eta_of, from_eta, num, special_values, table

RID = "round-v0.5"
LOG = TREE / "data" / "gate-logs" / "run006-round05-contact-events.json"
RUN005_LOG = TREE / "data" / "gate-logs" / "run005-round04-branch-equalization.json"
LABELS = ("B1", "B2", "B3", "B4")
ROUND04_GLOBAL = 0.9989037504763583     # Round 04's audited global scale (RUN-005 E2)


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

def build_world():
    f = load(RID)
    md = lambda n: restore_escapes(f[n].decode("utf-8"))
    rows = lambda n: list(csv.reader(io.StringIO(f[n].decode("utf-8"))))
    return {
        "readme": md("README.md"), "report": md("reports/ROUND_05_REPORT.md"),
        "source": f["src/event_system.py"].decode("utf-8"),
        "summary": json.loads(f["data/round05_summary.json"]),
        "gradients": [[float(v) for v in r] for r in rows("data/event_branch_gradients.csv")[1:]],
        "jacobian": [[float(v) for v in r] for r in rows("data/event_kkt_jacobian.csv")],
        "trials": json.loads(f["data/perturbation_convergence.json"]),
        "curve": np.array([[float(v) for v in r] for r in rows("data/event_candidate_curve.csv")[1:]]),
        "legs": ((3 + 4 * math.sqrt(3)) / 18, (4 + math.sqrt(3)) / 6),
        "cert_scale_divisor": 1.0048,
        "rot_sign": 1,
        "formula_B": same_B,
    }


def p_of(S):
    P = S["event_solution"]["parameters"]
    return np.array([P["l1"], P["l2"], P["beta_radians"], P["delta_radians"]])


def x_of(S):
    E = S["event_solution"]
    return np.concatenate([p_of(S), [E["common_scale"]], E["branch_pressures"]])


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

def closed(w, p):
    """the report's four closed forms (section 3), complex-safe for complex-step derivatives"""
    A, B = legs(w)
    l1, l2, b, d = p
    a = b + d
    l3 = 1 - 2 * (l1 + l2)
    c, s = np.cos, np.sin
    q1 = (l1 * c(a), l1 * s(a))
    q2 = (q1[0] + l2 * c(b), q1[1] + l2 * s(b))
    q3 = (q2[0] + l3, q2[1])
    D = 2 * l1 * c(a) + 2 * l2 * c(b) + l3
    X = lambda ph, q: q[0] * c(ph) - q[1] * s(ph)
    Y = lambda ph, q: q[0] * s(ph) + q[1] * c(ph)
    H = lambda ph, q: X(ph, q) / A + Y(ph, q) / B
    f1, f2, f3 = math.pi / 2 - a, math.pi / 2 - b, 2 * math.pi / 3
    return np.array([H(f1, q3), H(f2, q3) - X(f2, q1) / A, -X(f3, q3) / A - Y(f3, q2) / w["formula_B"](B),
                     q2[1] / A + (D - q2[0]) / B])


def grad(w, p, h=1e-30):
    """complex-step gradients of the closed forms: exact to rounding"""
    G = np.zeros((4, 4))
    for j in range(4):
        z = np.array(p, dtype=complex)
        z[j] += 1j * h
        G[:, j] = closed(w, z).imag / h
    return G


def F(w, x):
    p, s, mu = x[:4], x[4], x[5:]
    return np.concatenate([closed(w, p) - s, [mu.sum() - 1], mu @ grad(w, p)])


def jacobian(w, x, h=1e-7):
    out = np.zeros((9, 9))
    for j in range(9):
        a, b = x.copy(), x.copy()
        a[j] += h
        b[j] -= h
        out[:, j] = (F(w, a) - F(w, b)) / (2 * h)
    return out


def newton(w, x, iters=60):
    for _ in range(iters):
        dx = np.linalg.solve(jacobian(w, x), -F(w, x))
        x = x + dx
        if np.linalg.norm(dx) < 1e-15:
            break
    return x


def exact_one_sided(w, p, phi, eps=1e-9):
    """the derivatives in phi of the pieces active just left and just right of phi: for a fixed vertex,
    d x/d phi = -y and d y/d phi = x after the rotation"""
    A, B = legs(w)
    Bm = w["formula_B"](B)          # the divisor of the min-y term, as exact_scale uses it
    sg = w["rot_sign"]              # exact_scale turns by rot_sign * phi
    out = []
    for side in (-1, 1):
        t = sg * (phi + side * eps)
        c, s = math.cos(t), math.sin(t)
        q = p @ np.array([[c, s], [-s, c]])
        x, y = q[:, 0], q[:, 1]
        i, j, k = int(np.argmax(x / A + y / B)), int(np.argmin(x)), int(np.argmin(y))
        c, s = math.cos(sg * phi), math.sin(sg * phi)
        q = p @ np.array([[c, s], [-s, c]])
        x, y = q[:, 0], q[:, 1]
        out.append(sg * ((-y[i] / A + x[i] / B) + y[j] / A - x[k] / Bm))
    return out


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

def v1_setup(w):
    failures = []
    S = w["summary"]
    A, B = legs(w)
    T = S["triangle"]
    if abs(T["A"] - A) > 1e-15 or abs(T["B"] - B) > 1e-15:
        failures.append(f"triangle ({T['A']!r}, {T['B']!r}), the certified legs over 1.0048 give ({A!r}, {B!r})")
    if abs(T["B_over_A"] - math.sqrt(3)) > 1e-15 or abs(B / A - math.sqrt(3)) > 1e-15:
        failures.append(f"B/A stated {T['B_over_A']!r}, computed {B / A!r}, sqrt 3 {math.sqrt(3)!r}")
    P = S["event_solution"]["parameters"]
    if abs(P["l3"] - (1 - 2 * (P["l1"] + P["l2"]))) > 1e-15:
        failures.append("l3 is not 1 - 2(l1 + l2)")
    if abs(P["alpha_radians"] - (P["beta_radians"] + P["delta_radians"])) > 1e-15:
        failures.append("alpha is not beta + delta")
    for k in ("alpha", "beta", "delta"):
        if abs(math.degrees(P[f"{k}_radians"]) - P[f"{k}_degrees"]) > 1e-12:
            failures.append(f"{k}: degrees are not the radians'")
    p = w["curve"]
    q = from_eta(p_of(S))
    if p.shape != q.shape or np.abs(p - q).max() > 1e-15:
        failures.append("the logged curve is not the event parameters")
    total = np.linalg.norm(np.diff(p, axis=0), axis=1).sum()
    if abs(total - 1) > 1e-12:
        failures.append(f"the logged curve has length {total!r}")
    r = (p * np.array([-1.0, 1.0]))[::-1]
    if np.abs((r - r[0]) - (p - p[0])).max() > 1e-12:
        failures.append("the logged curve is not its own mirror image traversed backwards")
    return failures


def event_phases(p):
    return [math.pi / 2 - (p[2] + p[3]), math.pi / 2 - p[2], 2 * math.pi / 3, 3 * math.pi / 2]


def v2_phases(w):
    failures = []
    E = w["summary"]["event_solution"]
    for lab, got, deg, want in zip(LABELS, E["event_phases_radians"], E["event_phases_degrees"], event_phases(p_of(w["summary"]))):
        if abs(got - want) > 1e-15:
            failures.append(f"{lab}: event phase {got!r}, the event's {want!r}")
        if abs(math.degrees(got) - deg) > 1e-12:
            failures.append(f"{lab}: degrees {deg!r} are not the radians'")
    for led, ph in zip(w["summary"]["contact_ledgers"], E["event_phases_radians"]):
        if led["phase"] != ph:
            failures.append(f"{led['label']}: the ledger's phase is not the event's")
    return failures


def v3_values(w):
    failures = []
    E = w["summary"]["event_solution"]
    p = p_of(w["summary"])
    cf, ex = closed(w, p), special_values(w, p)
    for lab, a, b, v in zip(LABELS, cf, ex, E["branch_values"]):
        if abs(a - b) > 1e-14:
            failures.append(f"{lab}: the closed form gives {a!r}, the exact scale at the event {b!r}")
        if abs(a - v) > 1e-14:
            failures.append(f"{lab}: the closed form gives {a!r}, stated {v!r}")
    bv = E["branch_values"]
    if abs(E["common_scale"] - np.mean(bv)) > 1e-14:
        failures.append(f"the common scale {E['common_scale']!r} is not the branch values'")
    if abs(E["branch_spread"] - (max(bv) - min(bv))) > 1e-18:
        failures.append(f"the spread {E['branch_spread']!r} is not the branch values' {max(bv) - min(bv)!r}")
    return failures


def v4_global(w):
    failures = []
    S = w["summary"]
    E = S["event_solution"]
    p = from_eta(p_of(S))
    s, at = exact_min(w, p)
    if abs(s - E["common_scale"]) > 1e-13:
        failures.append(f"the least exact scale is {s!r}, the common scale {E['common_scale']!r}")
    mins = sorted(branch_minima(w, p, k=5))
    for lab, ph in zip(LABELS, E["event_phases_radians"]):
        v, x = min(mins, key=lambda m: abs(math.remainder(m[1] * math.pi - ph, 2 * math.pi)))
        if abs(math.remainder(x * math.pi - ph, 2 * math.pi)) > 1e-12 or abs(v - E["common_scale"]) > 1e-13:
            failures.append(f"{lab}: the nearest local minimum is ({x * math.pi:.13f}, {v!r})")
    if mins[4][0] - E["common_scale"] < 1e-6:
        failures.append(f"a fifth local minimum {mins[4]} lies within 1e-6 of the common scale")
    audit = S["full_phase_audit"]
    vals = grid_values(w, p, audit["phase_grid"])
    k = int(vals.argmin())
    if abs(vals.min() - audit["grid_minimum"]) > 1e-15 or abs(2 * math.pi * k / audit["phase_grid"] - audit["grid_phase"]) > 1e-12:
        failures.append(f"the {audit['phase_grid']}-phase grid's minimum is {vals.min()!r} at {2 * math.pi * k / audit['phase_grid']!r}")
    if audit["event_scale_minus_grid_minimum"] != E["common_scale"] - audit["grid_minimum"]:
        failures.append("the audit's difference is not the difference it names")
    return failures


def v5_gradients(w):
    G = grad(w, p_of(w["summary"]))
    d = np.abs(G - np.array(w["gradients"])).max()
    return [] if d <= 1e-12 else [f"the logged gradients differ from the closed forms' by {d:.2e}"]


def v6_kkt(w):
    failures = []
    S = w["summary"]
    E = S["event_solution"]
    mu = np.array(E["branch_pressures"])
    if (mu < 0).any() or abs(mu.sum() - 1) > 1e-14:
        failures.append(f"mu {mu.tolist()} is not on the simplex")
    r = float(np.linalg.norm(F(w, x_of(S))))
    if r > 1e-14:
        failures.append(f"the nine equations leave {r:.2e} at the stated solution")
    if E["equation_residual_norm"] > 1e-14:
        failures.append(f"the stated residual {E['equation_residual_norm']!r} is not at rounding level")
    for label, got, want in (("improvement over Round 04", E["improvement_over_round4"], E["common_scale"] - ROUND04_GLOBAL),
                             ("distance to one", E["distance_to_one"], 1 - E["common_scale"])):
        if abs(got - want) > 1e-18:
            failures.append(f"{label}: stated {got!r}, the difference is {want!r}")
    return failures


def v7_jacobian(w):
    failures = []
    J = np.array(w["jacobian"])
    D = w["summary"]["isolation_diagnostics"]
    mine = jacobian(w, x_of(w["summary"]))
    d = np.abs(mine - J).max() if J.shape == (9, 9) else float("inf")
    if d > 1e-7:
        failures.append(f"the logged Jacobian differs from this run's by {d:.2e}")
    sv = np.linalg.svd(J, compute_uv=False)
    if list(D["jacobian_shape"]) != [9, 9] or D["jacobian_rank"] != np.linalg.matrix_rank(J):
        failures.append(f"shape {D['jacobian_shape']} and rank {D['jacobian_rank']}, the logged Jacobian's {J.shape} and {np.linalg.matrix_rank(J)}")
    if len(D["singular_values"]) != 9 or np.abs(sv - np.array(D["singular_values"])).max() > 1e-12:
        failures.append("the stated singular values are not the logged Jacobian's")
    if abs(D["condition_number"] - D["singular_values"][0] / D["singular_values"][-1]) > 1e-9:
        failures.append("the stated condition number is not the singular values' ratio")
    svm = np.linalg.svd(mine, compute_uv=False)
    if np.abs(svm - sv).max() > 1e-7:
        failures.append(f"this run's singular values differ from the stated by {np.abs(svm - sv).max():.2e}")
    return failures


def v8_roots(w):
    failures = []
    D = w["summary"]["isolation_diagnostics"]
    T = w["trials"]
    ok = [t for t in T if t["success"]]
    if len(T) != D["perturbation_trials"] or len(ok) != D["successful_trials"]:
        failures.append(f"{len(ok)} of {len(T)} trials succeed, stated {D['successful_trials']} of {D['perturbation_trials']}")
    if max(t["distance_from_reference"] for t in T) != D["maximum_root_distance"]:
        failures.append("the stated largest distance is not the trials'")
    if any(t["residual_norm"] > 1e-14 for t in ok):
        failures.append("a successful trial leaves a residual above rounding level")
    x0 = x_of(w["summary"])
    rng = np.random.default_rng(6)
    far = 0
    for _ in range(40):
        try:
            far += np.abs(newton(w, x0 + 1e-2 * rng.normal(size=9)) - x0).max() > 1e-12
        except np.linalg.LinAlgError:           # a singular step: this start does not return
            far += 1
    if far:
        failures.append(f"{far} of 40 Newton runs from starts perturbed by 1e-2 do not return to the stated root")
    return failures


def v9_cusps(w):
    failures = []
    S = w["summary"]
    p = from_eta(p_of(S))
    for led in S["contact_ledgers"]:
        got = list(contacts(w, p, led["phase"]))
        want = [led["left_contact_vertices"], led["bottom_contact_vertices"], led["hypotenuse_contact_vertices"]]
        if got != want:
            failures.append(f"{led['label']}: contacts {got}, recorded {want}")
    for c in S["cusp_slopes"]:
        lo, hi = exact_one_sided(w, p, c["phase"])
        if abs(lo - c["left_derivative"]) > 1e-6 or abs(hi - c["right_derivative"]) > 1e-6:
            failures.append(f"{c['label']}: one-sided derivatives ({lo:.9f}, {hi:.9f}), stated "
                            f"({c['left_derivative']:.9f}, {c['right_derivative']:.9f})")
        if not lo < 0 < hi:
            failures.append(f"{c['label']}: not a kink minimum ({lo}, {hi})")
    return failures


def v10_consistency(w):
    failures = []
    S = w["summary"]
    E, D = S["event_solution"], S["isolation_diagnostics"]
    P = E["parameters"]
    rd, rp = gate.squash(w["readme"]), gate.squash(w["report"])

    def same(label, got, want, tol):
        if got is None or abs(got - want) > tol:
            failures.append(f"{label}: stated {got}, the data {want!r}")
    same("README scale", num(r"事件共同尺度：(\d+\.\d+)", rd), E["common_scale"], 6e-16)
    same("README distance", num(r"距離Wetzel認證尺度1：(\d+\.\d+e-?\d+)", rd), E["distance_to_one"], 6e-16)
    same("README spread", num(r"四事件尺度差距：(\d+\.\d+e-?\d+)", rd), E["branch_spread"], 6e-28)
    same("README rank", num(r"完整9×9Jacobian秩：(\d+)", rd), D["jacobian_rank"], 0)
    same("README sigma min", num(r"最小奇異值：(\d+\.\d+e-?\d+)", rd), D["singular_values"][-1], 6e-15)
    m = re.search(r"附近擾動初值回到同一根：(\d+)/(\d+)", rd)
    if not m or (int(m.group(1)), int(m.group(2))) != (D["successful_trials"], D["perturbation_trials"]):
        failures.append(f"README trials {m.groups() if m else None}")
    m = re.search(r"\(\\phi_1,\\phi_2,\\phi_3,\\phi_4\)=\(([\d.]+),([\d.]+),([\d.]+),([\d.]+)\)", rp)
    if not m or max(abs(float(a) - b) for a, b in zip(m.groups(), E["event_phases_radians"])) > 6e-13:
        failures.append("the report's event phases are not the data's")
    m = re.search(r"\(([\d.]+)\^\\circ,([\d.]+)\^\\circ,120\^\\circ,270\^\\circ\)", rp)
    if not m or max(abs(float(a) - b) for a, b in zip(m.groups(), E["event_phases_degrees"])) > 6e-10:
        failures.append("the report's event degrees are not the data's")
    for k, tol in (("l1", 6e-16), ("l2", 6e-16), ("l3", 6e-16)):
        same(f"report {k}", num(rf"{k[0]}_{k[1]}=(\d+\.\d+)", rp), P[k], tol)
    for k in ("beta", "delta", "alpha"):
        same(f"report {k}", num(rf"\\{k}=(\d+\.\d+)\^\\circ", rp), P[f"{k}_degrees"], 6e-13)
    same("report s_event", num(r"s_\{\\mathrm\{event\}\}=(\d+\.\d+)", rp), E["common_scale"], 6e-16)
    same("report spread", num(r"\\max_rm_r-\\min_rm_r=(\d+\.\d+e-?\d+)", rp), E["branch_spread"], 6e-28)
    same("report improvement", num(r"s_\{\\mathrm\{event\}\}-s_5\^\{\(4\)\}=(\d+\.\d+e-?\d+)", rp), E["improvement_over_round4"], 6e-22)
    same("report distance", num(r"1-s_\{\\mathrm\{event\}\}=(\d+\.\d+e-?\d+)", rp), E["distance_to_one"], 6e-16)
    same("report residual", num(r"\\\|F\\\|=(\d+\.\d+e-?\d+)", rp), E["equation_residual_norm"], 6e-29)
    same("report sigma min", num(r"\\sigma_\{\\min\}=(\d+\.\d+e-?\d+)", rp), D["singular_values"][-1], 6e-15)
    same("report condition", num(r"\\kappa\(J\)=(\d+\.\d+)", rp), D["condition_number"], 6e-10)
    same("report rank", num(r"\\operatorname\{rank\}J=(\d+)", rp), D["jacobian_rank"], 0)
    same("report max distance", num(r"最大解距離：(\d+\.\d+e-?\d+)", rp), D["maximum_root_distance"], 6e-28)
    same("report grid", num(r"使用\$(\d+)\$個相位點", rp), S["full_phase_audit"]["phase_grid"], 0)
    for row, mu in zip(table(w["report"], "| 分支 | 壓力 |"), E["branch_pressures"]):
        same(f"report pressure {row[0]}", float(row[1]), mu, 6e-13)
    for row, g in zip(table(w["report"], "\\partial_{l_1}"), w["gradients"]):
        if max(abs(float(c) - x) for c, x in zip(row[1:5], g)) > 6e-10:
            failures.append(f"report gradient row {row[0]} is not the data's")
    for row, led in zip(table(w["report"], "| 分支 | 相位 | 角度 |"), S["contact_ledgers"]):
        got = [json.loads(c) for c in row[3:6]]
        want = [led["left_contact_vertices"], led["bottom_contact_vertices"], led["hypotenuse_contact_vertices"]]
        if row[0] != led["label"] or got != want or abs(float(row[1]) - led["phase"]) > 6e-10 or \
                abs(float(row[2]) - led["phase_degrees"]) > 6e-7:
            failures.append(f"report contact row {row[0]} is not the data's")
    for row, c in zip(table(w["report"], "左方向導數"), S["cusp_slopes"]):
        if abs(float(row[1]) - c["left_derivative"]) > 6e-10 or abs(float(row[2]) - c["right_derivative"]) > 6e-10:
            failures.append(f"report slope row {row[0]} is not the data's")
    src = gate.squash(w["source"])
    for phrase in ("p=(l1,l2,beta,delta)", "s=commonscale", "mu_1,...,mu_4=branchpressures", "m_r(p)=s,r=1,...,4",
                   "sum_rmu_r=1", "sum_rmu_rgradm_r(p)=0", "nineequationsinnineunknowns"):
        if phrase not in src:
            failures.append(f"event_system.py does not state {phrase}")
    for phrase in ("m_r(p)-s=0", "\\sum_r\\mu_r-1=0", "\\sum_r\\mu_r\\nabla_pm_r=0", "(l_1,l_2,\\beta,\\delta,s,\\mu_1,\\ldots,\\mu_4)"):
        if phrase not in rp:
            failures.append(f"the report does not state {phrase}")
    return failures


CHECKS = {"V1": v1_setup, "V2": v2_phases, "V3": v3_values, "V4": v4_global, "V5": v5_gradients, "V6": v6_kkt,
          "V7": v7_jacobian, "V8": v8_roots, "V9": v9_cusps, "V10": v10_consistency}
GEOM = {"legs", "cert_scale_divisor", "rot_sign", "formula_B"}
DEPENDS = {
    "V1": {"summary", "curve", "legs", "cert_scale_divisor"},
    "V2": {"summary"},
    "V3": {"summary"} | GEOM,
    "V4": {"summary"} | GEOM,
    "V5": {"summary", "gradients"} | GEOM,
    "V6": {"summary"} | GEOM,
    "V7": {"summary", "jacobian"} | GEOM,
    "V8": {"summary", "trials"} | GEOM,
    "V9": {"summary"} | GEOM,
    "V10": {"readme", "report", "source", "summary", "gradients"},
}


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

def findings(w):
    S = w["summary"]
    E = S["event_solution"]
    x0 = x_of(S)
    p = p_of(S)
    G = grad(w, p)
    U, sv, Vt = np.linalg.svd(G)
    v = Vt[-1]
    mu = np.array(E["branch_pressures"])
    # second-order sufficient condition for max z s.t. m_r(p) >= z: with every mu_r > 0 and G of rank 3, the critical
    # cone is the line through G's right null vector v (since mu . 1 = 1, G d = dz 1 forces dz = 0); the Lagrangian's
    # curvature along it is v' (sum mu_r Hess m_r) v, from central differences of the complex-step gradients
    h = 1e-5
    curv = float(v @ ((mu @ grad(w, p + h * v)) - (mu @ grad(w, p - h * v))) / (2 * h))
    curv_closed = float((mu @ closed(w, p + 1e-4 * v) + mu @ closed(w, p - 1e-4 * v) - 2 * mu @ closed(w, p)) / 1e-8)
    rng = np.random.default_rng(0)
    radius = {}
    for scale in (1e-3, 1e-2, 3e-2):
        back = 0
        for _ in range(40):
            try:
                back += np.abs(newton(w, x0 + scale * rng.normal(size=9)) - x0).max() <= 1e-12
            except np.linalg.LinAlgError:
                pass
        radius[scale] = int(back)
    # near the root, is the least scale over all phases the least of the four closed forms (the identities stable),
    # and never above the root's?
    stab = {}
    for size in (1e-5, 1e-4):
        gaps, rises = [], []
        dists = []
        for _ in range(60):
            step = size * rng.normal(size=4)        # Gaussian, standard deviation `size` in each coordinate
            q = p + step
            full = exact_min(w, from_eta(q))[0]
            gaps.append(abs(full - closed(w, q).min()))
            rises.append(full - E["common_scale"])
            dists.append(float(np.linalg.norm(step)))
        stab[str(size)] = {"max_full_minus_closed": float(max(gaps)), "max_rise": float(max(rises)),
                           "distance_min_median_max": [min(dists), float(np.median(dists)), max(dists)]}
    r5 = json.loads(RUN005_LOG.read_text(encoding="utf-8"))["findings"]["exact_phase_equalization"]
    r4 = json.loads(load("round-v0.4.1")["data/round04_summary.json"])["after"]
    curve = from_eta(p)
    start = lambda eta, s: np.concatenate([eta, [s], np.ones(4) / 4])
    A, B = legs(w)
    T = S["triangle"]
    return {
        "precision": {
            "triangle_A_B": max(abs(T["A"] - A), abs(T["B"] - B)),
            "B_over_A_minus_sqrt3": max(abs(B / A - math.sqrt(3)), abs(T["B_over_A"] - math.sqrt(3))),
            "curve_minus_parameters": float(np.abs(w["curve"] - curve).max()),
            "closed_minus_exact_at_events": float(np.abs(closed(w, p) - special_values(w, p)).max()),
            "closed_minus_stated": float(np.abs(closed(w, p) - np.array(E["branch_values"])).max()),
            "least_scale": list(exact_min(w, curve)),
            "gradients_minus_logged": float(np.abs(G - np.array(w["gradients"])).max()),
            "equations_at_stated": float(np.linalg.norm(F(w, x0))),
            "jacobian_minus_logged": float(np.abs(jacobian(w, x0) - np.array(w["jacobian"])).max())},
        "newton_from_stated": float(np.abs(newton(w, x0) - x0).max()),
        "newton_from_run005": float(np.abs(newton(w, start(np.array(r5["eta"]), r5["least_scale_3600"][0])) - x0).max()),
        "newton_from_round04": float(np.abs(newton(w, start(eta_of(r4["parameters"]), r4["global_scale"])) - x0).max()),
        "gradient_rank": {"singular_values": sv.tolist(), "left_null_vs_mu": float(np.abs(U[:, -1] / U[:, -1].sum() - mu).max())},
        "second_order": {"curvature_along_null": curv, "curvature_from_values": curv_closed,
                         "null_direction": v.tolist()},
        "newton_returns_of_40": {str(k): v_ for k, v_ in radius.items()},
        "identity_stability": stab,
        "run005": {"value_minus_event": r5["least_scale_3600"][0] - E["common_scale"],
                   "eta_minus_event": (np.array(r5["eta"]) - p).tolist()},
        "one_sided_exact_minus_stated": [[a - c["left_derivative"], b - c["right_derivative"]] for (a, b), c in
                                         zip((exact_one_sided(w, curve, c["phase"]) for c in S["cusp_slopes"]), S["cusp_slopes"])],
        "grid_holds_event_phase": {lab: (S["full_phase_audit"]["phase_grid"] * ph / (2 * math.pi)) % 1 == 0
                                   for lab, ph in zip(LABELS, E["event_phases_radians"])},
        "fifth_minimum_above": sorted(branch_minima(w, curve, k=5))[4][0] - E["common_scale"],
        "exact_jacobian": exact_jacobian_report(w, x0),
        "one_sided_differences_step_1e-8": one_sided_differences(w, curve, S["cusp_slopes"]),
        "critical_line_probes": critical_line_probes(w, p, v, E["common_scale"]),
        "identity_margins": identity_margins(w, curve, E["event_phases_radians"]),
    }


def hessians(w, p, h=1e-4):
    """the four Hessians: Richardson-extrapolated central differences of the complex-step gradients"""
    Hs = np.zeros((4, 4, 4))
    for k in range(4):
        e = np.zeros(4)
        e[k] = 1
        d = lambda hh: (grad(w, p + hh * e) - grad(w, p - hh * e)) / (2 * hh)
        Hs[:, :, k] = (4 * d(h / 2) - d(h)) / 3
    return Hs


def exact_jacobian(w, x, h=1e-4):
    """the nine equations' Jacobian with exact first-derivative blocks (complex step, and the exact 0, 1 and -1
    entries) and the second-derivative block from hessians()"""
    p, mu = x[:4], x[5:]
    G = grad(w, p)
    J = np.zeros((9, 9))
    J[:4, :4] = G
    J[:4, 4] = -1
    J[4, 5:] = 1
    J[5:, :4] = np.einsum("r,rjk->jk", mu, hessians(w, p, h))
    J[5:, 5:] = G.T
    return J


def exact_jacobian_report(w, x):
    Je = exact_jacobian(w, x)
    sv = np.linalg.svd(Je, compute_uv=False)
    D = w["summary"]["isolation_diagnostics"]
    return {"hessian_self_consistency": float(np.abs(hessians(w, x[:4], 1e-4) - hessians(w, x[:4], 2e-4)).max()),
            "singular_values": sv.tolist(), "condition_number": float(sv[0] / sv[-1]),
            "stated_sigma_min_minus_exact": D["singular_values"][-1] - float(sv[-1]),
            "stated_singular_values_max_minus_exact": float(np.abs(np.array(D["singular_values"]) - sv).max()),
            "stated_condition_minus_exact": D["condition_number"] - float(sv[0] / sv[-1]),
            "logged_minus_exact": float(np.abs(np.array(w["jacobian"]) - Je).max()),
            "this_run_difference_minus_exact": float(np.abs(jacobian(w, x) - Je).max())}


def one_sided_differences(w, curve, slopes, h=1e-8):
    """one-sided differences of the exact scale at step h, against the stated derivatives"""
    out = []
    for c in slopes:
        s0 = exact_scale(w, curve, c["phase"])
        lo = (s0 - exact_scale(w, curve, c["phase"] - h)) / h
        hi = (exact_scale(w, curve, c["phase"] + h) - s0) / h
        out.append({"label": c["label"], "left_minus_stated": lo - c["left_derivative"],
                    "right_minus_stated": hi - c["right_derivative"]})
    return out


def critical_line_probes(w, p, v, s0):
    """along the straight line p + t v, and along the level curve m1 = m2 = m3 = m4 through the root (each point
    corrected back onto it within the plane orthogonal to v by Newton's method)"""
    basis = np.linalg.svd(v[None, :])[2][1:]          # an orthonormal basis of the complement of v
    line, level = {}, {}
    for t in (1e-4, 1e-3, 1e-2):
        vals = []
        for sgn in (1, -1):
            q = p + sgn * t * v
            vals.append([float(closed(w, q).min() - s0), float(abs(exact_min(w, from_eta(q))[0] - closed(w, q).min()))])
        line[str(t)] = vals
    for t in (1e-3, 3e-3, 1e-2):
        vals = []
        for sgn in (1, -1):
            c = np.zeros(3)
            for _ in range(30):
                q = p + sgn * t * v + c @ basis
                r = closed(w, q)
                g = grad(w, q)
                dc = np.linalg.solve((g[:3] - g[3]) @ basis.T, -(r[:3] - r[3]))
                c = c + dc
                if np.abs(dc).max() < 1e-16:
                    break
            q = p + sgn * t * v + c @ basis
            r = closed(w, q)
            vals.append([float(r.min() - s0), float(r.max() - r.min()), float(abs(exact_min(w, from_eta(q))[0] - r.min()))])
        level[str(t)] = vals
    return {"straight_line_drop_and_full_minus_closed": line, "level_curve_drop_spread_and_full_minus_closed": level}


def identity_margins(w, curve, phases):
    """at each event, how far each term's active vertex (or tied pair) lies from the next vertex: positive margins
    keep the identities, and with them the event phases, as the parameters move a little"""
    A, B = legs(w)
    out = []
    for ph in phases:
        c, s = math.cos(ph), math.sin(ph)
        q = curve @ np.array([[c, s], [-s, c]])
        terms = {"hypotenuse": -(q[:, 0] / A + q[:, 1] / B), "left": q[:, 0], "bottom": q[:, 1]}
        m = {}
        for k, vals in terms.items():
            srt = np.sort(vals)
            tied = int(np.sum(srt - srt[0] < 1e-9))
            m[k] = {"active": tied, "margin": float(srt[tied] - srt[0])}
        out.append(m)
    return out


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

def drill_cases(base):
    def edit(key, fn):
        v = copy.deepcopy(base[key])
        fn(v)
        return v
    E = lambda s: s["event_solution"]
    s_scale = edit("summary", lambda s: E(s).__setitem__("common_scale", E(s)["common_scale"] + 1e-11))
    s_phase = edit("summary", lambda s: E(s)["event_phases_radians"].__setitem__(0, E(s)["event_phases_radians"][0] + 1e-9))
    s_l1 = edit("summary", lambda s: E(s)["parameters"].__setitem__("l1", E(s)["parameters"]["l1"] + 1e-9))
    s_mu = edit("summary", lambda s: E(s)["branch_pressures"].__setitem__(1, E(s)["branch_pressures"][1] + 1e-3))
    s_sv = edit("summary", lambda s: s["isolation_diagnostics"]["singular_values"].__setitem__(8, 3.6072190e-02))
    s_led = edit("summary", lambda s: s["contact_ledgers"][2].__setitem__("bottom_contact_vertices", [3]))
    s_slope = edit("summary", lambda s: s["cusp_slopes"][1].__setitem__("right_derivative", s["cusp_slopes"][1]["right_derivative"] + 1e-4))
    s_imp = edit("summary", lambda s: E(s).__setitem__("improvement_over_round4", 6.6561505e-09))
    s_deg = edit("summary", lambda s: E(s)["event_phases_degrees"].__setitem__(1, E(s)["event_phases_degrees"][1] + 1e-6))
    s_grid = edit("summary", lambda s: s["full_phase_audit"].__setitem__("grid_minimum", s["full_phase_audit"]["grid_minimum"] + 1e-14))
    grads = copy.deepcopy(base["gradients"])
    grads[2][1] += 1e-6
    jac = copy.deepcopy(base["jacobian"])
    jac[6][3] += 1e-4
    trials = copy.deepcopy(base["trials"])
    trials[27]["distance_from_reference"] = 8.312345999720775e-14
    cv = base["curve"].copy()
    cv[3, 0] += 1e-9
    rep = base["report"].replace("s_{\\mathrm{event}}\n=\n0.998903757132509\n}", "s_{\\mathrm{event}}\n=\n0.998903757132519\n}")
    readme = base["readme"].replace("最小奇異值：3.607219048175e-02", "最小奇異值：3.607219048185e-02")
    source = base["source"].replace("This yields nine equations", "This yields eight equations")
    for a, b in ((rep, base["report"]), (readme, base["readme"]), (source, base["source"])):
        assert a != b, "a plant must hit the raw text"
    a, b = base["legs"]
    return [
        ("K0-unmodified", "control", {}, set()),
        ("K1-report-crlf", "control", {"report": base["report"].replace("\n", "\r\n")}, set()),
        ("K2-curve-copied", "control", {"curve": base["curve"].copy()}, set()),
        ("K3-readme-trailing-spaces", "control", {"readme": "\n".join(l + " " for l in base["readme"].split("\n"))}, set()),
        ("M1-triangle-leg", "mutation", {"legs": (a * 1.0005, b)}, {"V1", "V3", "V4", "V5", "V6", "V7", "V8", "V9"}),
        ("M2-common-scale", "mutation", {"summary": s_scale}, {"V3", "V4", "V6", "V8", "V10"}),
        ("M3-event-phase", "mutation", {"summary": s_phase}, {"V2", "V4", "V10"}),
        ("M4-length-l1", "mutation", {"summary": s_l1}, {"V1", "V3", "V4", "V5", "V6", "V8", "V10"}),
        ("M5-pressure", "mutation", {"summary": s_mu}, {"V6", "V7", "V8", "V10"}),
        ("M6-singular-value", "mutation", {"summary": s_sv}, {"V7", "V10"}),
        ("M7-contact-ledger", "mutation", {"summary": s_led}, {"V9", "V10"}),
        ("M8-cusp-slope", "mutation", {"summary": s_slope}, {"V9", "V10"}),
        ("M9-improvement", "mutation", {"summary": s_imp}, {"V6", "V10"}),
        ("M10-grid-minimum", "mutation", {"summary": s_grid}, {"V4"}),
        ("M11-gradient-entry", "mutation", {"gradients": grads}, {"V5", "V10"}),
        ("M12-jacobian-entry", "mutation", {"jacobian": jac}, {"V7"}),
        ("M13-trial-distance", "mutation", {"trials": trials}, {"V8"}),
        ("M14-curve-point", "mutation", {"curve": cv}, {"V1"}),
        ("M15-report-digit", "mutation", {"report": rep}, {"V10"}),
        ("M16-readme-digit", "mutation", {"readme": readme}, {"V10"}),
        ("M17-source-count", "mutation", {"source": source}, {"V10"}),
        ("M18-formula-term", "mutation", {"formula_B": lambda B: 2 * B}, {"V3", "V4", "V5", "V6", "V7", "V8", "V9"}),
        ("M19-event-degrees", "mutation", {"summary": s_deg}, {"V2", "V10"}),
    ]


def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("--drill", action="store_true")
    ap.add_argument("--blind", metavar="CID")
    ap.add_argument("--log", default=None)
    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[:6]:
            print("    ", x)
    log = {"run": "RUN-006",
           "subject": "Round 05 (v0.5): the contact-event equations and the nonsmooth KKT system of the five-link",
           "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 else findings(w)}
    if log["findings"]:
        print(json.dumps(log["findings"], indent=1)[:5000])
    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) + "\n", encoding="utf-8", newline="\n")


if __name__ == "__main__":
    main()