"""RUN-005 - Round 04 (v0.4.1): the phase-branch ledger and the four-branch equalization of the five-link.

Round 04 takes Round 03's five-link (RUN-004), tracks the four lowest local minima of its phase landscape - B1 and B2
at low phases, B3 at 120 deg, B4 at 270 deg - as functions m_r(eta) of eta = (l1, l2, beta, delta), and solves the
epigraph problem max z subject to m_r(eta) >= z (SLSQP). It reports the equalized five-link at 0.998903750476, the
four tracked values level to 2.734812376559e-12, a finite-difference sensitivity matrix G_rj = dm_r/deta_j, a
"branch pressure" vector mu on the simplex with sum_r mu_r grad m_r ~ 0, contact ledgers per branch, and a
resolution table. Its own cusp audit corrects the level claim: the tracker stopped beside the 270-deg cusp, the true
global scale is 0.998903750476358 at phi ~ 3pi/2, 1.301910645957e-08 below the tracked values, and the four true
branch heights lie in a band of about 1e-8.

Checks (each returns a list of failures; empty = pass)
  E1 parameters     the logged curve is the "after" parameters (l3 = 1 - 2 l1 - 2 l2, directions (alpha, beta, 0,
                    -beta, -alpha)), of length 1, its own mirror image traversed backwards; the "before" parameters
                    are Round 03's five-link
  E2 global         the least exact scale over all phases is each curve's stated global scale, to 1e-10, and the
                    after curve's is at the audit's phase, to 1e-8
  E3 minima         the after curve's local minima over phases are the nine stated, each value to 1e-10 and phase
                    to 1e-8, and there are no others
  E4 resolution     each row: the strict local minima on that phase grid number as stated, and its global scale is
                    the least exact scale, to 1e-12
  E5 tracker        each tracked branch's phase (after) lies within 1e-6 rad of a true local minimum, and each tracked
                    value (before and after) lies above that minimum by at most 2e-8; the stated spread is the
                    tracked values' spread
  E6 derived        the improvement, the distance to 1 and the README's statement of the audit's tracker
                    overestimate are the differences they name
  E7 sensitivity    the matrix in the summary is the csv's; each entry is this run's central difference of the true
                    branch minimum in eta = (l1, l2, beta, delta), delta = alpha - beta, to 2e-4, and to 1e-6
  E8 pressure       mu lies on the simplex and mu G is the stated residual, with the stated norm; this run's own G
                    admits a least-norm simplex combination within 1e-3 of mu
  E9 contacts       at each branch's true local minimum the vertices on the leg x = 0, the leg y = 0 and the
                    hypotenuse are the recorded ones
  E10 consistency   the README, the report and the audit state the data's values, parameters, tables and ledgers,
                    and METHOD.md states the report's epigraph problem and pressure condition

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

import numpy as np
from scipy.optimize import minimize

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, golden, mirror
from run004_round03_symmetric_links import (branch_minima, contacts, five_link, grid_values, load, same_B,
                                            seven_link, strict_minima)

RID = "round-v0.4.1"
LOG = TREE / "data" / "gate-logs" / "run005-round04-branch-equalization.json"
RUN004_LOG = TREE / "data" / "gate-logs" / "run004-round03-symmetric-links.json"
LABELS = ("B1", "B2", "B3", "B4")
FD_STEP = 1e-6
FINE_G = 1e-6          # far above this run's derivatives' stability across steps, far below the printed rows' gaps


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

def build_world():
    f, f3 = load(RID), load("round-v0.3")
    md = lambda n: restore_escapes(f[n].decode("utf-8"))
    rows = list(csv.DictReader(io.StringIO(f["data/equalized_5link_curve.csv"].decode("utf-8"))))
    return {
        "readme": md("README.md"), "report": md("reports/ROUND_04_REPORT.md"), "audit": md("reports/ROUND_04_CUSP_AUDIT.md"),
        "method": md("src/METHOD.md"),
        "summary": json.loads(f["data/round04_summary.json"]),
        "csv_matrix": [[float(v) for v in r.values()] for r in csv.DictReader(io.StringIO(
            f["data/branch_sensitivity_matrix.csv"].decode("utf-8")))],
        "curve": np.array([[float(r["x"]), float(r["y"])] for r in rows]),
        "round03": json.loads(f3["data/round03_summary.json"])["best_5link"]["parameters"],
        "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 curve_of(prm):
    return five_link(prm["l1"], prm["l2"], prm["l3"], prm["outer_angle_degrees"], prm["inner_angle_degrees"])


def eta_of(prm):
    b = math.radians(prm["inner_angle_degrees"])
    return np.array([prm["l1"], prm["l2"], b, math.radians(prm["outer_angle_degrees"]) - b])


def from_eta(e):
    l1, l2, b, d = e
    return five_link(l1, l2, 1 - 2 * l1 - 2 * l2, math.degrees(b + d), math.degrees(b))


def branch_value(w, p, phase, half=0.02, n=801):
    """the true local minimum of the exact scale nearest a phase: a scan of +-half, then golden section"""
    xs = np.linspace(phase - half, phase + half, n)
    v = np.array([exact_scale(w, p, x) for x in xs])
    i = int(v.argmin())
    return golden(lambda t: exact_scale(w, p, t), xs[max(i - 1, 0)], xs[min(i + 1, n - 1)])


def all_minima(w, p):
    """every local minimum over phases: an 8,192-phase scan, each refined by golden section, as (phase, value)"""
    return sorted((x * math.pi, v) for v, x in branch_minima(w, p, n=8192, k=10 ** 6))


def nearest(minima, phase):
    return min(minima, key=lambda m: abs(math.remainder(m[0] - phase, 2 * math.pi)))


def sensitivity(w, prm, phases, h=FD_STEP):
    e0 = eta_of(prm)
    G = np.zeros((4, 4))
    for j in range(4):
        a, b = e0.copy(), e0.copy()
        a[j] += h
        b[j] -= h
        pa, pb = from_eta(a), from_eta(b)
        for r, ph in enumerate(phases):
            G[r, j] = (branch_value(w, pa, ph)[0] - branch_value(w, pb, ph)[0]) / (2 * h)
    return G


def simplex_least_norm(G):
    res = minimize(lambda m: float(np.sum((m @ G) ** 2)), np.ones(len(G)) / len(G), method="SLSQP",
                   bounds=[(0, 1)] * len(G), constraints=[{"type": "eq", "fun": lambda m: m.sum() - 1}],
                   options={"ftol": 1e-20, "maxiter": 500})
    return res.x


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

def e1_parameters(w):
    failures = []
    S = w["summary"]
    after = S["after"]["parameters"]
    p = w["curve"]
    q = curve_of(after)
    if p.shape != q.shape or np.abs((p - p[0]) - q).max() > 1e-12:
        failures.append("the logged curve is not the after parameters")
    for tag, prm in (("after", after), ("before", S["before"]["parameters"])):
        if abs(prm["l3"] - (1 - 2 * prm["l1"] - 2 * prm["l2"])) > 1e-12:
            failures.append(f"{tag}: l3 is not 1 - 2 l1 - 2 l2")
    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")
    for k, v in w["round03"].items():
        if abs(S["before"]["parameters"][k] - v) > 1e-12:
            failures.append(f"before {k} {S['before']['parameters'][k]!r} is not Round 03's {v!r}")
    return failures


def e2_global(w):
    failures = []
    S = w["summary"]
    s, at = exact_min(w, w["curve"])
    if abs(s - S["after"]["global_scale"]) > 1e-10:
        failures.append(f"after: least exact scale {s:.13f}, stated {S['after']['global_scale']:.13f}")
    stated_at = min(S["all_final_local_minima"], key=lambda m: m["scale"])["phase"]
    if abs(at - stated_at) > 1e-8:
        failures.append(f"after: least scale at phase {at:.10f}, the audit's {stated_at:.10f}")
    sb = exact_min(w, curve_of(S["before"]["parameters"]))[0]
    if abs(sb - S["before"]["global_scale"]) > 1e-10:
        failures.append(f"before: least exact scale {sb:.13f}, stated {S['before']['global_scale']:.13f}")
    return failures


def e3_minima(w):
    failures = []
    got = all_minima(w, w["curve"])
    stated = sorted((m["phase"], m["scale"]) for m in w["summary"]["all_final_local_minima"])
    if len(got) != len(stated):
        failures.append(f"{len(got)} local minima over phases, stated {len(stated)}")
    for ph, v in stated:
        gp, gv = nearest(got, ph)
        if abs(gv - v) > 1e-10 or abs(math.remainder(gp - ph, 2 * math.pi)) > 1e-8:
            failures.append(f"the minimum stated at ({ph:.10f}, {v:.13f}) is ({gp:.10f}, {gv:.13f})")
    return failures


def e4_resolution(w):
    failures = []
    s = exact_min(w, w["curve"])[0]
    for row in w["summary"]["resolution_stability"]:
        n = strict_minima(grid_values(w, w["curve"], row["phase_grid"])).size
        if n != row["local_minimum_count"]:
            failures.append(f"{row['phase_grid']} phases: {n} strict local minima, stated {row['local_minimum_count']}")
        if abs(row["global_scale"] - s) > 1e-12:
            failures.append(f"{row['phase_grid']} phases: global scale {row['global_scale']!r}, the least exact scale {s!r}")
    return failures


def e5_tracker(w):
    failures = []
    S = w["summary"]
    phases = S["after"]["tracked_branch_phases"]
    for tag, prm in (("after", S["after"]["parameters"]), ("before", S["before"]["parameters"])):
        mins = all_minima(w, curve_of(prm))
        for lab, ph, v in zip(LABELS, phases, S[tag]["tracked_branch_values"]):
            gp, gv = nearest(mins, ph)
            if tag == "after" and abs(math.remainder(gp - ph, 2 * math.pi)) > 1e-6:
                failures.append(f"after {lab}: tracked at phase {ph:.10f}, the nearest true minimum at {gp:.10f}")
            if not 0 <= v - gv <= 2e-8:
                failures.append(f"{tag} {lab}: tracked {v:.13f}, the true minimum {gv:.13f} ({v - gv:+.3e})")
    tv = S["after"]["tracked_branch_values"]
    if abs(S["after"]["branch_spread"] - (max(tv) - min(tv))) > 1e-18:
        failures.append(f"the stated spread {S['after']['branch_spread']!r} is not the tracked values' {max(tv) - min(tv)!r}")
    return failures


def e6_derived(w):
    failures = []
    S = w["summary"]
    a, b = S["after"]["global_scale"], S["before"]["global_scale"]
    for label, got, want in (("improvement", S["after"]["improvement"], a - b),
                             ("distance to one", S["after"]["distance_to_one"], 1 - a)):
        if abs(got - want) > 1e-15:
            failures.append(f"{label}: stated {got!r}, the difference is {want!r}")
    over = num(r"追蹤器高估量：(\d+\.\d+e-?\d+)", gate.squash(w["readme"]))
    want = min(S["after"]["tracked_branch_values"]) - a
    if over is None or abs(over - want) > 6e-21:
        failures.append(f"the tracker's overestimate is stated {over}, the least tracked value minus the global scale is {want!r}")
    return failures


def e7_sensitivity(w):
    failures = []
    S = w["summary"]
    G = np.array(S["sensitivity_matrix"])
    if S["sensitivity_columns"] != ["l1", "l2", "beta_radians", "delta_radians"]:
        failures.append(f"columns {S['sensitivity_columns']}")
    if G.shape != (4, 4) or np.abs(G - np.array(w["csv_matrix"])).max() > 1e-15:
        failures.append("the summary's matrix is not branch_sensitivity_matrix.csv")
    mine = sensitivity(w, S["after"]["parameters"], S["after"]["tracked_branch_phases"])
    for r, lab in enumerate(LABELS):
        d = np.abs(mine[r] - G[r]).max()
        if d > 2e-4:
            failures.append(f"{lab}: the stated row differs from this run's derivatives by {d:.2e}")
        if d > FINE_G:
            # fixed text: a plant must make a new failure, not reprint this one with another number
            failures.append(f"{lab}: the stated row is not this run's derivatives to {FINE_G:.0e}")
    return failures


def e8_pressure(w):
    failures = []
    L = w["summary"]["branch_pressure_ledger"]
    mu, G = np.array(L["weights"]), np.array(w["summary"]["sensitivity_matrix"])
    if (mu < 0).any() or abs(mu.sum() - 1) > 1e-12:
        failures.append(f"mu {mu.tolist()} is not on the simplex")
    res = mu @ G
    if np.abs(res - np.array(L["stationarity_residual"])).max() > 1e-15:
        failures.append(f"mu G is {res.tolist()}, stated {L['stationarity_residual']}")
    if abs(np.linalg.norm(L["stationarity_residual"]) - L["stationarity_norm"]) > 1e-15:
        failures.append(f"the stated norm {L['stationarity_norm']!r} is not the residual's")
    mine = simplex_least_norm(sensitivity(w, w["summary"]["after"]["parameters"], w["summary"]["after"]["tracked_branch_phases"]))
    if np.abs(mine - mu).max() > 1e-3:
        failures.append(f"this run's G gives weights {np.round(mine, 6).tolist()}, stated {np.round(mu, 6).tolist()}")
    return failures


def e9_contacts(w):
    failures = []
    mins = all_minima(w, w["curve"])
    for led in w["summary"]["branch_contact_ledgers"]:
        ph, _ = nearest(mins, led["phase"])
        got = list(contacts(w, w["curve"], ph))
        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}")
    return failures


def num(pattern, text):
    m = re.search(pattern, text)
    return float(m.group(1).rstrip(".")) if m else None


def table(text, header):
    """the rows of the markdown table whose header line contains `header`, as lists of stripped cells"""
    lines = text.splitlines()
    i = next(k for k, l in enumerate(lines) if l.startswith("|") and header in l)
    out = []
    for l in lines[i + 2:]:
        if not l.startswith("|"):
            break
        out.append([c.strip() for c in l.strip().strip("|").split("|")])
    return out


def e10_consistency(w):
    failures = []
    S = w["summary"]
    rd, rp, au = gate.squash(w["readme"]), gate.squash(w["report"]), gate.squash(w["audit"])
    A, B = S["after"], S["before"]

    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}")
    tv = A["tracked_branch_values"]
    over = min(tv) - A["global_scale"]
    same("README round 3", num(r"第3輪五連桿：(\d+\.\d+)", rd), B["global_scale"], 6e-13)
    same("README round 4", num(r"第4輪等高五連桿：(\d+\.\d+)", rd), A["global_scale"], 6e-13)
    same("README improvement", num(r"提升：(\d+\.\d+e-?\d+)", rd), A["improvement"], 6e-17)
    same("README distance", num(r"距離認證尺度1：(\d+\.\d+e-?\d+)", rd), A["distance_to_one"], 6e-16)
    same("README spread", num(r"四分支等高誤差：(\d+\.\d+e-?\d+)", rd), A["branch_spread"], 6e-25)
    same("README tracked spread", num(r"追蹤epigraph分支差距：(\d+\.\d+e-?\d+)", rd), A["branch_spread"], 6e-25)
    same("README global", num(r"完整相位全域尺度：(\d+\.\d+)", rd), A["global_scale"], 6e-16)
    same("README overestimate", num(r"追蹤器高估量：(\d+\.\d+e-?\d+)", rd), over, 6e-21)
    same("report s5(3)", num(r"s_5\^\{\(3\)\}=(\d+\.\d+)", rp), B["global_scale"], 6e-13)
    same("report s5(4)", num(r"s_5\^\{\(4\)\}=(\d+\.\d+)", rp), A["global_scale"], 6e-13)
    same("report improvement", num(r"\\Deltas=(\d+\.\d+e-?\d+)", rp), A["improvement"], 6e-17)
    same("report distance", num(r"1-s_5\^\{\(4\)\}=(\d+\.\d+e-?\d+)", rp), A["distance_to_one"], 6e-16)
    same("report spread", num(r"\\max_rm_r-\\min_rm_r=(\d+\.\d+e-?\d+)", rp), A["branch_spread"], 6e-25)
    same("report residual norm", num(r"\\right\\\|=(\d+\.\d+e-?\d+)", rp), S["branch_pressure_ledger"]["stationarity_norm"], 6e-20)
    same("report conclusion s5", num(r"結論.*?s_5=(\d+\.\d+)", rp), A["global_scale"], 6e-13)
    for tag, prm in (("第3輪", B["parameters"]), ("第4輪", A["parameters"])):
        part = rp[rp.index(tag + "：", rp.index("參數重新分配")):]
        m = re.search(r"\(l_1,l_2,l_3\)=\(([\d.]+),([\d.]+),([\d.]+)\)", part)
        got = [float(x) for x in m.groups()] if m else None
        if not got or max(abs(a - b) for a, b in zip(got, (prm["l1"], prm["l2"], prm["l3"]))) > 6e-11:
            failures.append(f"report {tag} lengths {got}")
        m = re.search(r"\(\\beta,\\alpha\)=\(([\d.]+)\^\\circ,([\d.]+)\^\\circ\)", part)
        got = [float(x) for x in m.groups()] if m else None
        if not got or max(abs(a - b) for a, b in zip(got, (prm["inner_angle_degrees"], prm["outer_angle_degrees"]))) > 6e-9:
            failures.append(f"report {tag} angles {got}")
    rows = table(w["report"], "相位區域")
    for i, row in enumerate(rows):
        for got, want, tol in ((row[2], B["tracked_branch_values"][i], 6e-13), (row[3], tv[i], 6e-13),
                               (row[4], A["tracked_branch_phases"][i], 6e-10)):
            same(f"report branch table {row[1]}", float(got), want, tol)
    for row, led in zip(table(w["report"], "左邊界接觸"), S["branch_contact_ledgers"]):
        got = [json.loads(c) for c in row[2:5]]
        want = [led["left_contact_vertices"], led["bottom_contact_vertices"], led["hypotenuse_contact_vertices"]]
        if row[1] != led["label"] or got != want:
            failures.append(f"report contact table {row[1]}: {got} vs the data's {want}")
    for row, g in zip(table(w["report"], "\\partial_{l_1}"), S["sensitivity_matrix"]):
        if max(abs(float(c) - x) for c, x in zip(row[1:5], g)) > 6e-9:
            failures.append(f"report sensitivity row {row[0]} is not the data's")
    for row, mu in zip(table(w["report"], "分支壓力"), S["branch_pressure_ledger"]["weights"]):
        same(f"report pressure {row[1]}", float(row[2]), mu, 6e-11)
    m = re.search(r"\$\$([\d.]+),\\quad([\d.]+),\\quad([\d.]+),\\quad([\d.]+)\.\$\$", au)
    if not m or max(abs(float(a) - b) for a, b in zip(m.groups(), tv)) > 6e-16:
        failures.append(f"the audit's tracked values {m.groups() if m else None} vs the data's {tv}")
    same("audit s*", num(r"s_\\ast=(\d+\.\d+)", au), A["global_scale"], 6e-16)
    same("audit phase", num(r"\\phi=(\d+\.\d+)", au), min(S["all_final_local_minima"], key=lambda x: x["scale"])["phase"], 6e-16)
    same("audit overestimate", num(r"之差：\$\$(\d+\.\d+e-?\d+)", au), over, 6e-21)
    for formula in ("\\max_{\\eta,z}z", "m_r(\\eta)\\gez,\\qquadr=1,\\ldots,4", "\\mu_r\\ge0", "\\sum_r\\mu_r\\nabla_\\etam_r\\approx0"):
        for label, text in (("METHOD.md", gate.squash(w["method"])), ("the report", rp)):
            if formula not in text:
                failures.append(f"{label} does not state {formula}")
    return failures


CHECKS = {"E1": e1_parameters, "E2": e2_global, "E3": e3_minima, "E4": e4_resolution, "E5": e5_tracker,
          "E6": e6_derived, "E7": e7_sensitivity, "E8": e8_pressure, "E9": e9_contacts, "E10": e10_consistency}
GEOM = {"legs", "cert_scale_divisor", "rot_sign", "formula_B"}
DEPENDS = {
    "E1": {"summary", "curve", "round03"},
    "E2": {"summary", "curve"} | GEOM,
    "E3": {"summary", "curve"} | GEOM,
    "E4": {"summary", "curve"} | GEOM,
    "E5": {"summary"} | GEOM,
    "E6": {"summary", "readme"},
    "E7": {"summary", "csv_matrix"} | GEOM,
    "E8": {"summary"} | GEOM,
    "E9": {"summary", "curve"} | GEOM,
    "E10": {"readme", "report", "audit", "method", "summary"},
}


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

def special_values(w, e):
    """the exact scale at the four phases where the branches' contact identities switch: link 1 parallel to the leg
    x = 0 (pi/2 - alpha), link 2 parallel to it (pi/2 - beta), the symmetry axis parallel to the hypotenuse's normal
    (2 pi/3) and to the leg x = 0's (3 pi/2)"""
    l1, l2, b, d = e
    p = from_eta(e)
    return np.array([exact_scale(w, p, math.pi / 2 - (b + d)), exact_scale(w, p, math.pi / 2 - b),
                     exact_scale(w, p, 2 * math.pi / 3), exact_scale(w, p, 3 * math.pi / 2)])


def closed_gradient(w, e, h=1e-7):
    """central differences of the four closed forms (smooth in eta where the contact identities hold)"""
    G = np.zeros((4, 4))
    for j in range(4):
        a, b = e.copy(), e.copy()
        a[j] += h
        b[j] -= h
        G[:, j] = (special_values(w, a) - special_values(w, b)) / (2 * h)
    return G


def one_sided(w, p, at, e=1e-6):
    s = exact_scale(w, p, at)
    return [(s - exact_scale(w, p, at - e)) / e, (exact_scale(w, p, at + e) - s) / e]


def findings(w):
    S = w["summary"]
    A = S["after"]
    after, before = w["curve"], curve_of(S["before"]["parameters"])
    e0 = eta_of(A["parameters"])
    alpha, beta = e0[2] + e0[3], e0[2]
    special = [math.pi / 2 - alpha, math.pi / 2 - beta, 2 * math.pi / 3, 3 * math.pi / 2]
    mins = all_minima(w, after)
    true4 = [nearest(mins, ph) for ph in A["tracked_branch_phases"]]
    steps = {h: sensitivity(w, A["parameters"], A["tracked_branch_phases"], h=h) for h in (1e-5, 1e-6, 1e-7)}
    mine = steps[FD_STEP]
    G = np.array(S["sensitivity_matrix"])
    mu_mine = simplex_least_norm(mine)
    # the epigraph problem with each branch at its exact phase: smooth in eta
    v0 = special_values(w, e0)

    def epigraph(x0):
        cons = [{"type": "ineq", "fun": (lambda x, r=r: special_values(w, x[:4])[r] - x[4])} for r in range(4)]
        return minimize(lambda x: -x[4], x0, method="SLSQP", constraints=cons, options={"ftol": 1e-16, "maxiter": 400})
    res = epigraph(np.append(e0, v0.min()))
    xs = res.x[:4]
    eq = from_eta(xs)
    veq = special_values(w, xs)
    seq = exact_min(w, eq)
    # at the levelled point: the gradients' rank, the multiplier, and second-order probes
    Gs = closed_gradient(w, xs)
    sv = np.linalg.svd(Gs, compute_uv=False)
    null = np.linalg.svd(Gs)[2][-1]
    mu_s = simplex_least_norm(Gs)
    rng = np.random.default_rng(0)
    dirs = rng.normal(size=(200, 4))
    dirs /= np.linalg.norm(dirs, axis=1)[:, None]
    restarts = []
    for k in range(3):
        x0 = xs + 1e-3 * rng.normal(size=4)
        rr = epigraph(np.append(x0, special_values(w, x0).min()))
        restarts.append(float(special_values(w, rr.x[:4]).min() - veq.min()))
    run004 = json.loads(RUN004_LOG.read_text(encoding="utf-8"))["findings"]["seven_link_search"]
    x7 = max(run004["runs"], key=lambda r: r["value"])["parameters"]
    s7 = exact_min(w, seven_link(x7[0], x7[1], x7[2], 1 - 2 * (x7[0] + x7[1] + x7[2]), x7[3], x7[4], x7[5]))[0]
    v05 = num(r"事件共同尺度：(\d+\.\d+)", gate.squash(load("round-v0.5")["README.md"].decode("utf-8")))
    stated = S["all_final_local_minima"]
    s_after, at_after = exact_min(w, after)
    grid_rows = {}
    for row in S["resolution_stability"]:
        n = row["phase_grid"]
        vals = grid_values(w, after, n)
        grid_rows[n] = {"contains_3pi_2_under_2pi_k_over_n": n % 4 == 0, "grid_minimum": float(vals.min()),
                        "argmin_k": int(vals.argmin()), "three_quarters_n": 3 * n // 4,
                        "stated_minus_least": row["global_scale"] - s_after}
    stated_col = [row["global_scale"] for row in S["resolution_stability"]]
    return {
        "global": [s_after, at_after],
        "global_phase_minus_3pi_2": at_after - 3 * math.pi / 2,
        "audit_phase_minus_3pi_2": min(stated, key=lambda m: m["scale"])["phase"] - 3 * math.pi / 2,
        "curve_minus_parameters": float(np.abs((after - after[0]) - curve_of(A["parameters"])).max()),
        "minima_max_deviation": {
            "value": max(abs(nearest(mins, m["phase"])[1] - m["scale"]) for m in stated),
            "phase": max(abs(math.remainder(nearest(mins, m["phase"])[0] - m["phase"], 2 * math.pi)) for m in stated)},
        "branch_minima_at_special_phases": {
            lab: {"true_phase": tp, "special_phase": sp, "difference": math.remainder(tp - sp, 2 * math.pi),
                  "stated_phase_difference": math.remainder(nearest([(m["phase"], m["scale"]) for m in stated], tp)[0] - sp, 2 * math.pi),
                  "true_value": tv, "value_at_special_phase": float(v0[i]), "one_sided_slopes": one_sided(w, after, tp)}
            for i, (lab, (tp, tv), sp) in enumerate(zip(LABELS, true4, special))},
        "true_branch_spread": max(v for _, v in true4) - min(v for _, v in true4),
        "tracked_phase_minus_true": [math.remainder(ph - tp, 2 * math.pi) for ph, (tp, _) in zip(A["tracked_branch_phases"], true4)],
        "tracked_minus_true": {"after": [t - v for t, (_, v) in zip(A["tracked_branch_values"], true4)],
                               "before": [t - nearest(all_minima(w, before), ph)[1] for t, ph in
                                          zip(S["before"]["tracked_branch_values"], A["tracked_branch_phases"])]},
        "resolution_rows": grid_rows,
        "resolution_stated_column_scatter": max(stated_col) - min(stated_col),
        "sensitivity_this_run": mine.tolist(),
        "sensitivity_step_stability": max(float(np.abs(steps[h] - mine).max()) for h in steps),
        "sensitivity_closed_form_vs_true_minima": float(np.abs(closed_gradient(w, e0) - mine).max()),
        "sensitivity_row_max_difference": dict(zip(LABELS, np.abs(mine - G).max(axis=1).tolist())),
        "pressure_this_run": {"weights": mu_mine.tolist(), "residual_norm": float(np.linalg.norm(mu_mine @ mine)),
                              "max_weight_difference": float(np.abs(mu_mine - np.array(S["branch_pressure_ledger"]["weights"])).max()),
                              "stated_weights_on_this_G": float(np.linalg.norm(np.array(S["branch_pressure_ledger"]["weights"]) @ mine))},
        "exact_phase_equalization": {
            "slsqp_status": int(res.status), "slsqp_message": str(res.message),
            "eta": xs.tolist(), "alpha_degrees": math.degrees(xs[2] + xs[3]),
            "beta_degrees": math.degrees(xs[2]), "l3": 1 - 2 * xs[0] - 2 * xs[1],
            "eta_distance_from_round": float(np.linalg.norm(xs - e0)),
            "branch_values": veq.tolist(), "spread": float(veq.max() - veq.min()),
            "least_scale_3600": list(seq), "least_scale_65536": exact_min(w, eq, scan=65536)[0],
            "above_round": seq[0] - A["global_scale"],
            "branches_8192": branch_minima(w, eq, k=6),
            "gradient_singular_values": sv.tolist(),
            "least_norm_multiplier": mu_s.tolist(), "multiplier_residual": float(np.linalg.norm(mu_s @ Gs)),
            "null_direction_change_at_1e-4": [float(special_values(w, xs + s * 1e-4 * null).min() - veq.min()) for s in (1, -1)],
            "random_directions_best_change_at_1e-5": float(max(special_values(w, xs + 1e-5 * d).min() - veq.min() for d in dirs)),
            "restarts_from_1e-3_minus_value": restarts,
            "round_v05_event_scale": v05, "round_v05_minus_this": None if v05 is None else v05 - seq[0]},
        "run004_seven_link": {"value": s7, "above_round": s7 - A["global_scale"], "above_this_run": s7 - seq[0]},
    }


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

def drill_cases(base):
    def edit(key, fn):
        v = copy.deepcopy(base[key])
        fn(v)
        return v
    s_global = edit("summary", lambda s: s["after"].__setitem__("global_scale", s["after"]["global_scale"] + 1e-9))
    s_track = edit("summary", lambda s: s["after"]["tracked_branch_values"].__setitem__(
        2, s["after"]["tracked_branch_values"][2] - 1e-8))
    s_min = edit("summary", lambda s: s["all_final_local_minima"][2].__setitem__(
        "phase", s["all_final_local_minima"][2]["phase"] + 1e-6))
    s_res = edit("summary", lambda s: s["resolution_stability"][0].__setitem__("local_minimum_count", 9))
    s_g = edit("summary", lambda s: s["sensitivity_matrix"][3].__setitem__(2, s["sensitivity_matrix"][3][2] + 1e-3))
    s_mu = edit("summary", lambda s: s["branch_pressure_ledger"]["weights"].__setitem__(
        0, s["branch_pressure_ledger"]["weights"][0] + 2e-3))
    s_led = edit("summary", lambda s: s["branch_contact_ledgers"][3].__setitem__("bottom_contact_vertices", [0]))
    s_ang = edit("summary", lambda s: s["after"]["parameters"].__setitem__(
        "inner_angle_degrees", s["after"]["parameters"]["inner_angle_degrees"] + 1e-6))
    s_bef = edit("summary", lambda s: s["before"]["parameters"].__setitem__(
        "l2", s["before"]["parameters"]["l2"] + 1e-9))
    s_imp = edit("summary", lambda s: s["after"].__setitem__("improvement", 1.4937880820e-04))
    cv = base["curve"].copy()
    cv[2, 1] += 1e-6
    rep = base["report"].replace("s_5^{(4)}=0.998903750476", "s_5^{(4)}=0.998903750477")
    readme = base["readme"].replace("追蹤器高估量：1.301910645957e-08", "追蹤器高估量：1.301910645958e-08")
    audit = base["audit"].replace("\\phi=4.712388980384675", "\\phi=4.712388980384685")
    method = base["method"].replace("m_r(\\eta)\\ge z", "m_r(\\eta)\\le z")
    for a, b in ((rep, base["report"]), (readme, base["readme"]), (audit, base["audit"]), (method, base["method"])):
        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)}, {"E2", "E3", "E4", "E5", "E7"}),
        ("M2-global-scale", "mutation", {"summary": s_global}, {"E2", "E6", "E10"}),
        ("M3-tracked-value", "mutation", {"summary": s_track}, {"E5", "E6", "E10"}),
        ("M4-minimum-phase", "mutation", {"summary": s_min}, {"E3"}),
        ("M5-resolution-count", "mutation", {"summary": s_res}, {"E4"}),
        ("M6-sensitivity-entry", "mutation", {"summary": s_g}, {"E7", "E8", "E10"}),
        ("M7-pressure-weight", "mutation", {"summary": s_mu}, {"E8", "E10"}),
        ("M8-contact-ledger", "mutation", {"summary": s_led}, {"E9", "E10"}),
        ("M9-after-angle", "mutation", {"summary": s_ang}, {"E1", "E10"}),
        ("M10-before-length", "mutation", {"summary": s_bef}, {"E1", "E2", "E5", "E10"}),
        ("M11-improvement", "mutation", {"summary": s_imp}, {"E6", "E10"}),
        ("M12-curve-point", "mutation", {"curve": cv}, {"E1", "E2", "E3", "E4"}),
        ("M13-report-digit", "mutation", {"report": rep}, {"E10"}),
        ("M14-readme-digit", "mutation", {"readme": readme}, {"E6", "E10"}),
        ("M15-audit-phase", "mutation", {"audit": audit}, {"E10"}),
        ("M16-method-inequality", "mutation", {"method": method}, {"E10"}),
        ("M17-formula-term", "mutation", {"formula_B": lambda B: 2 * B}, {"E2", "E3", "E4", "E5", "E7", "E8", "E9"}),
    ]


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-005",
           "subject": "Round 04 (v0.4.1): the phase-branch ledger and the four-branch equalization 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()