"""RUN-003 - Round 02 (v0.2.2): the contact-pressure ledger, full phase scans, the topology-guided three-link
search, the exact three-normal formula, and the chirality correction.

Round 02 moves the grid program to 256 support directions and scans five curves through 360 rotation phases in the
certified triangle T/1.0048; reads the LP's dual as a ledger of contact pressures; searches a three-link family and
finds a candidate at 0.9933. Two addenda then correct it: an exact formula for the triangle (its three edge normals,
no grid) puts that candidate at 1.000837 and a re-optimized three-link at 1.007315 - above the certified scale 1 -
and a chirality correction explains both away: the triangle's certificate is for congruent copies, so a curve
counts as covered when either handedness fits, and the mirror image of the 1.007315 curve fits at 0.952939. The
corrected search maximizes the lesser handedness and ends at 0.998573794671. Everything is recomputed here with this
tree's own code; the package's program is read only after.

Checks (each returns a list of failures; empty = pass)
  Y1 ledger         at its recorded phase, the best grid candidate's 256-direction program gives its recorded scale
                    and translation; its dual has weight on exactly the three recorded directions, with the recorded
                    weights, and balances (sum lambda u = 0, sum lambda h_C = 1, s = sum lambda h_gamma); the
                    container features at those directions are the recorded ones; no phase gives less by more than 5e-7
  Y2 phase scans    each of the 1,080 logged rows of the three curves the documents define is the 256-direction
                    program's scale at its phase; for all five curves, the minimum, its phase, the maximum and the
                    barrier the summary states are the logged rows'
  Y3 the grid       the 256-direction grid misses the hypotenuse's normal, so on it the program is a relaxation of T;
                    the round's own exact values lie above its grid values
  Y4 exact formula  s(phi) = max(x/A + y/B) - min x/A - min y/B is the triangle's exact program for a polyline (it
                    equals the three-normal LP); with it, the addendum's and the correction's values reproduce
  Y5 three links    each logged three-link curve is its stated parameters: lengths summing to 1, directions
                    (alpha1, 0, -alpha3)
  Y6 consistency    the README, the report's tables and the addenda state the data's numbers
  Y7 local minima   each curve's count of strict local minima over its 360 logged phases is the summary's
  Y8 citations      the two works the report cites are what it says, and the preprint's convention is the congruent
                    copy, reflection allowed
  Y9 signatures     for the curves this tree rebuilds, the count of distinct sets of container features at the
                    dual's active directions over 360 phases is the summary's count of contact signatures

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/run003_round02_phase_scan_chirality.py --fetch      (once: the citation snapshot)
  python code/run003_round02_phase_scan_chirality.py [--drill [--log PATH] | --blind Y#]
"""
import argparse
import copy
import csv
import io
import json
import math
import pathlib
import re

import numpy as np
from scipy.optimize import linprog, minimize, minimize_scalar

from corpus import TREE, members, sha256, source_bytes
import gate
from run002_round01_support_skew_lp import LP_OPTIONS, global_min, rotate, support, unit

RID = "round-v0.2.2"
LOG = TREE / "data" / "gate-logs" / "run003-round02-phase-scan-chirality.json"
SNAPSHOT = TREE / "data" / "external" / "round02-citations.json"
CITED = ("2606.14625", "1101.5638")
CURVES = ("U_equal", "arc_270", "semicircle", "segment", "round1_adversarial")
# the scanned curves this tree can rebuild from the documents: U (Round 01's), the unit segment, Round 01's logged
# candidate. The two arcs are defined nowhere; no reading tried reproduces them (findings: arcs_closest_reading)
REBUILT = ("U_equal", "segment", "round1_adversarial")
SCAN_TOL = 1e-9             # a logged phase-scan scale against this tree's
EXACT_TOL = 1e-10           # an exact-formula value against the round's twelve decimals
FINE = 3600                 # rotation scan for an exact value, every local minimum then refined


# ---------------------------------------------------------------- the container and the programs

def legs(w):
    a, b = w["legs"]
    return a / w["cert_scale_divisor"], b / w["cert_scale_divisor"]


def triangle(w):
    """T/1.0048, right angle at the origin, the short leg on the x-axis (as RUN-002 found Round 01 placed it)"""
    A, B = legs(w)
    return np.array([[0.0, 0.0], [A, 0.0], [0.0, B]])


def grid(w):
    return 2 * math.pi * np.arange(w["M"]) / w["M"]


def lp_dual(p, phi, th, c):
    """the grid program at phase phi with its dual: (s, t, lambda); s evaluated exactly at the solver's t"""
    h = support(rotate(p, phi), th)
    u = unit(th)
    res = linprog([0, 0, 1], A_ub=np.column_stack([u, -c]), b_ub=-h,
                  bounds=[(None, None), (None, None), (0, None)], method="highs-ds", options=LP_OPTIONS)
    assert res.status == 0, res.message
    tx, ty, _ = res.x
    lhs = h + u @ np.array([tx, ty])
    pos = c > 1e-12
    s = float((lhs[pos] / c[pos]).max())
    return s, (float(tx), float(ty)), -res.ineqlin.marginals, h


def exact_scale(w, p, phi):
    """the triangle's exact program for a polyline at phase phi: the best translation puts min x on the leg x = 0
    and min y on the leg y = 0, and the hypotenuse x/A + y/B <= s binds at a vertex"""
    A, B = legs(w)
    q = rotate(p, w["rot_sign"] * phi)
    x, y = q[:, 0], q[:, 1]
    return float((x / A + y / B).max() - x.min() / A - y.min() / w["formula_B"](B))


def golden(f, a, b, iters=120):
    """golden-section minimum of f on [a, b]: unlike a parabolic step it cannot stall at a kink, and 120 steps
    shrink the bracket below 1e-25"""
    g = (math.sqrt(5) - 1) / 2
    c, d = b - g * (b - a), a + g * (b - a)
    fc, fd = f(c), f(d)
    for _ in range(iters):
        if fc < fd:
            b, d, fd = d, c, fc
            c = b - g * (b - a)
            fc = f(c)
        else:
            a, c, fc = c, d, fd
            d = a + g * (b - a)
            fd = f(d)
    x = (a + b) / 2
    return f(x), x


def exact_min(w, p, scan=FINE):
    """the least exact scale over all phases: a scan of 3,600 phases, every local minimum refined by golden section
    (the exact scale is a maximum of smooth pieces, so its minima are often kinks)"""
    key = ("exact", p.tobytes(), scan, w["legs"], w["cert_scale_divisor"], w["rot_sign"], id(w["formula_B"]))
    if key not in _memo:
        f = lambda x: exact_scale(w, p, x)
        phis = 2 * math.pi * np.arange(scan) / scan
        vals = np.array([f(x) for x in phis])
        step = 2 * math.pi / scan
        best = (float(vals.min()), float(phis[vals.argmin()]))
        for i in range(scan):
            if vals[i] <= vals[i - 1] and vals[i] <= vals[(i + 1) % scan]:
                v, x = golden(f, phis[i] - step, phis[i] + step)
                if v < best[0]:
                    best = (float(v), float(x % (2 * math.pi)))
        _memo[key] = best
    return _memo[key]


def mirror(p):
    return p * np.array([1.0, -1.0])


def three_link(l1, l2, l3, a1, a3):
    """directions (alpha1, 0, -alpha3), lengths (l1, l2, l3), from the origin"""
    pts = [np.zeros(2)]
    for length, ang in ((l1, a1), (l2, 0.0), (l3, -a3)):
        pts.append(pts[-1] + length * np.array([math.cos(ang), math.sin(ang)]))
    return np.array(pts)


def u_family(alpha, l):
    """the symmetric family (alpha, 0, -alpha), (l, 1 - 2l, l)"""
    return three_link(l, 1 - 2 * l, l, alpha, alpha)


_memo = {}


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

def load_package():
    return {m.path.split("/", 1)[1]: b for m, b in members(RID, source_bytes(RID), keep=lambda m: True)}


def round01_candidate():
    f = {m.path.split("/", 1)[1]: b for m, b in members("round-v0.1", source_bytes("round-v0.1"), keep=lambda m: True)}
    rows = list(csv.DictReader(io.StringIO(f["data/adversarial_best_curve.csv"].decode("utf-8"))))
    return np.array([[float(r["x"]), float(r["y"])] for r in rows])


def build_world():
    f = load_package()
    text = lambda n: f[n].decode("utf-8")
    table = lambda n: [{k: float(v) for k, v in r.items()} for r in csv.DictReader(io.StringIO(text(n)))]
    snap = json.loads(SNAPSHOT.read_text(encoding="utf-8")) if SNAPSHOT.exists() else None
    # the curves the documents define; Round 02's two arcs are not among them (see REBUILT)
    curves = {"U_equal": np.array([[0, 0], [1 / 3, 0], [1 / 3, 1 / 3], [0, 1 / 3]], dtype=float),
              "segment": np.array([[0, 0], [1, 0]], dtype=float),
              "round1_adversarial": round01_candidate()}
    return {
        "readme": text("README.md"), "report": text("reports/ROUND_02_REPORT.md"),
        "addendum": text("reports/ROUND_02_EXACT_ADDENDUM.md"),
        "correction": text("reports/ROUND_02_CHIRALITY_CORRECTION.md"),
        "summary": json.loads(f["data/round02_summary.json"]),
        "ledger": json.loads(f["data/best_contact_stress_ledger.json"]),
        "chirality": json.loads(f["data/chirality_corrected_summary.json"]),
        "exact_json": json.loads(f["data/exact_triangle_validation.json"]),
        "scans": {c: table(f"data/phase_scan_{c}.csv") for c in CURVES},
        "logged_curves": {k: np.array([[r["x"], r["y"]] for r in table(f"data/{k}.csv")])
                          for k in ("best_topology_guided_curve", "best_exact_three_link_curve",
                                    "best_congruence_three_link_curve")},
        "surface": [[float(x) for x in r.split(",")] for r in text("data/u_family_surface.csv").strip().split("\n")],
        "curves": curves,
        "snapshot": snap,
        "M": 256,
        "legs": ((3 + 4 * math.sqrt(3)) / 18, (4 + math.sqrt(3)) / 6),
        "cert_scale_divisor": 1.0048,
        "rot_sign": 1,
        "formula_B": lambda B: B,
    }


def container_c(w, th):
    return support(triangle(w), th)


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

def feature(w, theta):
    """which part of T supports direction theta: an edge whose outward normal it is, else a vertex"""
    v = triangle(w)
    d = np.roll(v, -1, axis=0) - v
    normals = np.arctan2(-d[:, 0], d[:, 1]) % (2 * math.pi)
    names = {0: "edge:0-1", 1: "edge:1-2", 2: "edge:0-2"}
    for i, n in enumerate(normals):
        if abs((theta - n + math.pi) % (2 * math.pi) - math.pi) < 1e-12:
            return names[i]
    return f"vertex:{int(np.argmax(v @ unit(np.array([theta]))[0]))}"


def y1_ledger(w):
    failures = []
    L = w["ledger"]
    p = w["logged_curves"]["best_topology_guided_curve"]
    th = grid(w)
    c = container_c(w, th)
    s, t, lam, h = lp_dual(p, L["phase"], th, c)
    if abs(s - L["critical_scale"]) > 1e-9:
        failures.append(f"at the recorded phase the 256-direction program gives {s:.12f}, recorded {L['critical_scale']:.12f}")
    if max(abs(a - b) for a, b in zip(t, L["translation"])) > 1e-9:
        failures.append(f"translation {tuple(round(x, 12) for x in t)} vs recorded {tuple(L['translation'])}")
    active = np.flatnonzero(lam > 1e-9)
    got = [round(math.degrees(th[i]), 6) for i in active]
    want = [round(x, 6) for x in L["contact_centers_degrees"]]
    if got != want:
        failures.append(f"dual weight on directions {got}, recorded {want}")
    else:
        for i, wt in zip(active, L["contact_weights"]):
            if abs(lam[i] - wt) > 1e-7:
                failures.append(f"weight at {math.degrees(th[i]):.5f} deg is {lam[i]:.9f}, recorded {wt:.9f}")
    u = unit(th)
    eq = float(np.linalg.norm(lam @ u))
    nm = abs(float(lam @ c) - 1)
    gap = abs(float(lam @ h) - s)
    if max(eq, nm, gap) > 1e-9:
        failures.append(f"the dual does not balance: |sum lambda u| {eq:.1e}, |sum lambda c - 1| {nm:.1e}, gap {gap:.1e}")
    feats = [feature(w, th[i]) for i in active]
    if feats != L["container_features"]:
        failures.append(f"container features {feats}, recorded {L['container_features']}")
    # and no phase gives less (RUN-002's W4 clause, same tolerance): the recorded scale is a minimum over phases
    low, at = grid_min(w, p)
    if low < L["critical_scale"] - 5e-7:
        failures.append(f"phase {at:.6f} gives {low:.9f}, below the recorded critical scale")
    return failures


def grid_min(w, p, scan=720):
    """the 256-direction program's minimum over phases: a scan, every local minimum refined by golden section"""
    key = ("gridmin", p.tobytes(), scan, w["M"], w["legs"], w["cert_scale_divisor"], w["rot_sign"])
    if key not in _memo:
        th = grid(w)
        c = container_c(w, th)
        f = lambda x: lp_dual(p, w["rot_sign"] * x, th, c)[0]
        phis = 2 * math.pi * np.arange(scan) / scan
        vals = np.array([f(x) for x in phis])
        step = 2 * math.pi / scan
        best = (float(vals.min()), float(phis[vals.argmin()]))
        for i in range(scan):
            if vals[i] <= vals[i - 1] and vals[i] <= vals[(i + 1) % scan]:
                v, x = golden(f, phis[i] - step, phis[i] + step, iters=70)
                if v < best[0]:
                    best = (float(v), float(x % (2 * math.pi)))
        _memo[key] = best
    return _memo[key]


def scan_rows(w, name):
    key = ("scan", name, w["M"], w["legs"], w["cert_scale_divisor"], w["rot_sign"])
    if key not in _memo:
        p, th = w["curves"][name], grid(w)
        c = container_c(w, th)
        _memo[key] = [lp_dual(p, w["rot_sign"] * (math.pi * k / 180), th, c)[0] for k in range(360)]
    return _memo[key]


def y2_scans(w):
    failures = []
    for name in CURVES:
        rows = w["scans"][name]
        if len(rows) != 360:
            failures.append(f"{name}: {len(rows)} rows, not 360")
            continue
        if name in REBUILT:
            mine = scan_rows(w, name)
            off = [k for k, (r, m) in enumerate(zip(rows, mine)) if abs(r["scale"] - m) > SCAN_TOL
                   or abs(r["phase"] - math.pi * k / 180) > 1e-12]
            if off:
                k = off[0]
                failures.append(f"{name}: {len(off)} of 360 rows differ, first at phase {k} deg: logged "
                                f"{rows[k]['scale']:.12f}, this tree {mine[k]:.12f}")
        s = w["summary"]["phase_scan"]["curves"][name]
        logged = [r["scale"] for r in rows]
        k = int(np.argmin(logged))
        stated = {"minimum_scale": min(logged), "minimum_phase": rows[k]["phase"], "maximum_scale": max(logged),
                  "phase_barrier": max(logged) - min(logged)}
        for key, v in stated.items():
            if abs(s[key] - v) > 1e-12:
                failures.append(f"{name}: summary {key} {s[key]!r}, the logged rows give {v!r}")
    return failures


def y3_grid(w):
    failures = []
    th = grid(w)
    hyp = math.radians(30)                      # the hypotenuse's outward normal in this placement
    if np.min(np.abs((th - hyp + math.pi) % (2 * math.pi) - math.pi)) < 1e-12:
        failures.append("the 256-direction grid holds the hypotenuse's normal: the program would be exact")
    # on the grid the program is a relaxation of T: its value never exceeds the exact one, at any phase
    for name in ("U_equal", "segment", "round1_adversarial"):
        p = w["curves"][name]
        worst = max(scan_rows(w, name)[k] - exact_scale(w, p, math.pi * k / 180) for k in range(360))
        if worst > 1e-12:
            failures.append(f"{name}: the grid value exceeds the exact one by {worst:.2e} at some phase")
    s_grid = w["summary"]["phase_scan"]["curves"]["U_equal"]["minimum_scale"]
    s_exact = exact_min(w, w["curves"]["U_equal"])[0]
    if not s_exact > s_grid:
        failures.append(f"U: exact {s_exact} not above the grid's {s_grid}")
    return failures


def y4_exact(w):
    failures = []
    A, B = legs(w)
    normals = np.array([math.radians(d) for d in (270, 30, 180)])
    c3 = support(triangle(w), normals)
    rng = np.random.default_rng(20261019)
    for _ in range(25):                          # the formula is the three-normal LP, for random polylines
        p = rng.normal(size=(int(rng.integers(2, 9)), 2)) * 0.3
        phi = float(rng.uniform(0, 2 * math.pi))
        lp = lp_dual(p, w["rot_sign"] * phi, normals, c3)[0]
        if abs(lp - exact_scale(w, p, phi)) > 1e-10:
            failures.append(f"formula {exact_scale(w, p, phi)!r} vs the three-normal LP {lp!r}")
            break
    ej, ch = w["exact_json"], w["chirality"]
    if max(abs(x - y) for x, y in zip(ej["exact_triangle_method"]["triangle_legs"], (A, B))) > 1e-15:
        failures.append(f"the addendum's legs {ej['exact_triangle_method']['triangle_legs']} are not T/1.0048's")
    # each stated value, recomputed as the minimum over all phases of the exact formula
    tg = w["logged_curves"]["best_topology_guided_curve"]
    ex = w["logged_curves"]["best_exact_three_link_curve"]
    cg = w["logged_curves"]["best_congruence_three_link_curve"]
    claims = [("U_equal", w["curves"]["U_equal"], ej["u_equal"]["exact_scale"]),
              ("the grid candidate", tg, ej["round02_grid_candidate"]["exact_scale"]),
              ("the optimized three-link", ex, ej["optimized_three_link"]["exact_scale"]),
              ("its mirror image", mirror(ex), ch["false_alarm_candidate"]["reflected"]["scale"]),
              ("the corrected three-link", cg, ch["corrected_three_link_search"]["orientation_preserving"]["scale"]),
              ("its mirror image", mirror(cg), ch["corrected_three_link_search"]["reflected"]["scale"])]
    for label, p, stated in claims:
        got = exact_min(w, p)[0]
        if abs(got - stated) > EXACT_TOL:
            failures.append(f"{label}: exact minimum {got:.12f}, stated {stated:.12f}")
    cs = ch["corrected_three_link_search"]
    if abs(cs["congruent_scale"] - min(cs["orientation_preserving"]["scale"], cs["reflected"]["scale"])) > 1e-15:
        failures.append("the corrected congruent scale is not the lesser handedness")
    if abs(cs["improvement_over_u_equal_congruent"] - (cs["congruent_scale"] - ch["u_equal"]["congruent_scale"])) > 1e-15:
        failures.append("the improvement over U is not the difference it names")
    if abs(cs["distance_to_certified_scale_one"] - (1 - cs["congruent_scale"])) > 1e-15:
        failures.append("the distance to 1 is not 1 - s_cong")
    for label, block in (("false alarm", ch["false_alarm_candidate"]), ("corrected", cs), ("U", ch["u_equal"])):
        gap = abs(block["orientation_preserving"]["scale"] - block["reflected"]["scale"])
        if abs(block["chirality_gap"] - gap) > 1e-15:
            failures.append(f"{label}: chirality gap {block['chirality_gap']!r} is not |s+ - s-| = {gap!r}")
    o3 = ej["optimized_three_link"]
    if abs(o3["distance_to_certified_scale_one"] - (1 - o3["exact_scale"])) > 1e-15:
        failures.append("the addendum's distance to 1 is not 1 - s_3link")
    return failures


def params_of(p):
    d = np.diff(p, axis=0)
    lengths = np.linalg.norm(d, axis=1)
    ang = np.arctan2(d[:, 1], d[:, 0])
    return lengths, ang


def y5_three_links(w):
    failures = []
    stated = {"best_topology_guided_curve": w["summary"]["topology_guided_search"]["best_parameters"],
              "best_exact_three_link_curve": w["exact_json"]["optimized_three_link"]["parameters"],
              "best_congruence_three_link_curve": w["chirality"]["corrected_three_link_search"]["parameters"]}
    for name, prm in stated.items():
        lengths, ang = params_of(w["logged_curves"][name])
        a1 = prm.get("alpha1", prm.get("alpha1_radians"))
        a3 = prm.get("alpha3", prm.get("alpha3_radians"))
        want_l = np.array([prm["l1"], prm["l2"], prm["l3"]])
        want_a = np.array([a1, 0.0, -a3])
        if abs(lengths.sum() - 1) > 1e-12 or abs(want_l.sum() - 1) > 1e-12:
            failures.append(f"{name}: lengths sum to {lengths.sum():.15f} (stated {want_l.sum():.15f})")
        if np.abs(lengths - want_l).max() > 1e-12 or np.abs((ang - want_a + math.pi) % (2 * math.pi) - math.pi).max() > 1e-12:
            failures.append(f"{name}: the curve's lengths {np.round(lengths, 10).tolist()} and directions "
                            f"{np.round(np.degrees(ang), 8).tolist()} are not its stated parameters")
    return failures


def num(pattern, text):
    m = re.search(pattern, text)
    return float(m.group(1).rstrip(".")) if m else None      # a sentence's full stop is not the number's


def y6_consistency(w):
    failures = []
    S, ch, ej, tg = w["summary"], w["chirality"], w["exact_json"], w["summary"]["topology_guided_search"]
    def same(label, got, want, digits):
        if got is None or abs(got - want) > 0.6 * 10 ** -digits:
            failures.append(f"{label}: stated {got}, the data {want:.{digits}f}")
    same("README U", num(r"等長 U 形臨界尺度：([\d.]+)", w["readme"]), tg["u_equal_fine_scale"], 10)
    same("README best", num(r"最佳三連桿臨界尺度：([\d.]+)", w["readme"]), tg["best_fine_scale"], 10)
    same("README grid candidate exact", num(r"舊網格候選精確尺度：([\d.]+)", w["readme"]), ej["round02_grid_candidate"]["exact_scale"], 12)
    same("README re-optimized", num(r"精確再優化三連桿：([\d.]+)", w["readme"]), ej["optimized_three_link"]["exact_scale"], 12)
    same("README corrected", num(r"修正後三連桿合同難度：([\d.]+)", w["readme"]), ch["corrected_three_link_search"]["congruent_scale"], 12)
    rows = re.findall(r"^\|\s*(\w+)\s*\|\s*([\d.]+)\s*\|\s*([\d.]+)\s*\|\s*(\d+)\s*\|\s*(\d+)\s*\|\s*(\d+)\s*\|\s*$",
                      w["report"], re.M)
    if len(rows) != 5:
        failures.append(f"the report's phase table has {len(rows)} rows, not 5")
    for name, mn, bar, lmin, jumps, sigs in rows:
        s = S["phase_scan"]["curves"].get(name)
        if s is None:
            failures.append(f"the report's table names {name}, not in the summary")
            continue
        got = (round(float(mn), 9), round(float(bar), 9), int(lmin), int(jumps), int(sigs))
        want = (round(s["minimum_scale"], 9), round(s["phase_barrier"], 9), s["local_minimum_count"],
                s["jump_count"], s["distinct_container_contact_signatures"])
        if got != want:
            failures.append(f"{name}: report {got} vs summary {want}")
    for label, pat, want, d in (("report U", r"s_U=([\d.]+)", tg["u_equal_fine_scale"], 10),
                                ("report best", r"s_\{\\mathrm\{topo\}\}=([\d.]+)", tg["best_fine_scale"], 10),
                                ("addendum U", r"s_U=([\d.]+)", ej["u_equal"]["exact_scale"], 12),
                                ("addendum re-optimized", r"s_\{\\mathrm\{3link\}\}=([\d.]+)", ej["optimized_three_link"]["exact_scale"], 12),
                                ("correction s_-", r"s_-=([\d.]+)", ch["false_alarm_candidate"]["reflected"]["scale"], 12),
                                ("correction U", r"s_\{\\mathrm\{cong\}\}\(U\)=([\d.]+)", ch["u_equal"]["congruent_scale"], 12)):
        text = w["addendum"] if label.startswith("addendum") else w["correction"] if label.startswith("correction") else w["report"]
        same(label, num(pat, gate.squash(text)), want, d)
    fa = ch["false_alarm_candidate"]
    rd = gate.squash(w["readme"])
    same("README s+", num(r"s_\+=(\d+\.\d+)", rd), fa["orientation_preserving"]["scale"], 9)
    same("README s-", num(r"s_-=(\d+\.\d+)", rd), fa["reflected"]["scale"], 9)
    same("README s_cong", num(r"s_\{cong\}=(\d+\.\d+)", rd), fa["congruent_scale"], 9)
    chis = [float(x) for x in re.findall(r"\\chi=(\d+\.\d+)", gate.squash(w["correction"]))]
    if [round(x, 12) for x in chis] != [round(fa["chirality_gap"], 12), round(ch["corrected_three_link_search"]["chirality_gap"], 12)]:
        failures.append(f"the correction's chi values {chis} vs the data's gaps")
    d3 = re.search(r"1-s_\{\\mathrm\{3link\}\}=(-?\d+\.\d+e-?\d+)", gate.squash(w["addendum"]))
    if not d3 or abs(float(d3.group(1)) - ej["optimized_three_link"]["distance_to_certified_scale_one"]) > 1e-14:
        failures.append(f"the addendum's 1 - s_3link {d3.group(1) if d3 else None} vs the data")
    ledger = re.findall(r"^\|\s*(\d)\s*\|\s*([\d.]+)\s*\|\s*([\d.]+)\s*\|\s*(\S+)\s*\|\s*$", w["report"], re.M)
    L = w["ledger"]
    for (k, deg, wt, feat), d0, w0, f0 in zip(ledger, L["contact_centers_degrees"], L["contact_weights"], L["container_features"]):
        if abs(float(deg) - d0) > 5e-6 or abs(float(wt) - w0) > 5e-9 or feat != f0:
            failures.append(f"ledger row {k}: report ({deg}, {wt}, {feat}) vs data ({d0}, {w0}, {f0})")
    if len(ledger) != 3:
        failures.append(f"the report's ledger has {len(ledger)} rows")
    return failures


def y7_local_minima(w):
    failures = []
    for name in CURVES:
        v = [r["scale"] for r in w["scans"][name]]
        if len(v) != 360:
            failures.append(f"{name}: {len(v)} logged phases, not 360")
            continue
        n = sum(1 for k in range(360) if v[k] < v[k - 1] and v[k] < v[(k + 1) % 360])
        stated = w["summary"]["phase_scan"]["curves"][name]["local_minimum_count"]
        if n != stated:
            failures.append(f"{name}: {n} strict local minima over the 360 logged phases, stated {stated}")
    return failures


def y8_citations(w):
    snap = w["snapshot"]
    if snap is None:
        return ["no citation snapshot - run with --fetch"]
    failures = []
    refs = {int(n): (a, t, rest) for n, a, t, rest in
            re.findall(r"^(\d+)\.\s+(.+?),\s+\*(.+?)\*,\s+(.+?)\.?\s*$", w["report"], re.M)}
    for n, aid in enumerate(CITED, 1):
        rec = snap["records"][aid]
        authors, title, rest = refs.get(n, ("", "", ""))
        if gate.norm_title(title) != gate.norm_title(rec["title"]):
            failures.append(f"ref {n}: title {title!r} vs arXiv {rec['title']!r}")
        fams = [re.sub(r"^(?:[A-Z]\.\s*)+", "", a).strip() for a in re.split(r",\s*(?:and\s+)?|\s+and\s+", authors) if a.strip()]
        if fams != rec["families"]:
            failures.append(f"ref {n}: authors {fams} vs arXiv {rec['families']}")
        if f"arXiv:{aid}" not in rest:
            failures.append(f"ref {n}: arXiv id not {aid}")
        year = re.findall(r",\s*((?:19|20)\d{2})$", rest)       # a year token, not the tail of an arXiv id
        if year and int(year[0]) not in rec["years"]:
            failures.append(f"ref {n}: year {year[0]} vs {rec['years']}")
    facts = snap["convention"]["facts"]
    if not facts.get("congruent"):
        failures.append("the preprint does not state its cover as containing a congruent copy of every unit arc")
    if not facts.get("reflection") or "rotations and reflections" not in facts["reflection"]:
        failures.append("the preprint does not state that its congruent-copy convention includes reflections")
    return failures


def surface_probe(w):
    """The U-family surface (14 x 21) is not rebuildable: its grid and its evaluation are stated nowhere. Its
    corners - alpha 65 or 115 degrees, l 0.18 or 0.42, rows running over l - reproduce exactly on a 128-direction
    grid with 36 phases; its interior does not on evenly spaced axes. Logged as a finding, not checked."""
    S = np.array(w["surface"])
    th = 2 * math.pi * np.arange(128) / 128
    c = container_c(w, th)
    out = {}
    alphas, ls = np.linspace(65, 115, S.shape[1]), np.linspace(0.18, 0.42, S.shape[0])
    for i, j in ((0, 0), (0, S.shape[1] - 1), (S.shape[0] - 1, 0), (S.shape[0] - 1, S.shape[1] - 1), (8, 6), (5, 10)):
        p = u_family(math.radians(alphas[j]), ls[i])
        v = min(lp_dual(p, 2 * math.pi * k / 36, th, c)[0] for k in range(36))
        out[f"{i},{j}"] = {"alpha_deg_if_even": float(alphas[j]), "l_if_even": float(ls[i]), "logged": float(S[i, j]),
                           "grid128_phases36": v, "difference": v - float(S[i, j])}
    return out


def signatures(w, name):
    """the distinct sets of container features (Y1's feature()) at the dual's active directions, over 360 phases"""
    key = ("signatures", name, w["M"], w["legs"], w["cert_scale_divisor"], w["rot_sign"])
    if key not in _memo:
        p, th = w["curves"][name], grid(w)
        c = container_c(w, th)
        seen = set()
        for k in range(360):
            lam = lp_dual(p, w["rot_sign"] * math.pi * k / 180, th, c)[2]
            seen.add(tuple(sorted({feature(w, th[i]) for i in np.flatnonzero(lam > 1e-9)})))
        _memo[key] = seen
    return _memo[key]


def y9_signatures(w):
    """the summary's count of container contact signatures, for the curves this tree rebuilds"""
    failures = []
    for name in REBUILT:
        n = len(signatures(w, name))
        stated = w["summary"]["phase_scan"]["curves"][name]["distinct_container_contact_signatures"]
        if n != stated:
            failures.append(f"{name}: {n} distinct sets of contact features over 360 phases, stated {stated}")
    return failures


CHECKS = {"Y1": y1_ledger, "Y2": y2_scans, "Y3": y3_grid, "Y4": y4_exact, "Y5": y5_three_links,
          "Y6": y6_consistency, "Y7": y7_local_minima, "Y8": y8_citations, "Y9": y9_signatures}
DEPENDS = {
    "Y1": {"ledger", "logged_curves", "M", "legs", "cert_scale_divisor", "rot_sign"},
    "Y2": {"scans", "summary", "curves", "M", "legs", "cert_scale_divisor", "rot_sign"},
    "Y3": {"summary", "curves", "M", "legs", "cert_scale_divisor", "rot_sign", "formula_B"},
    "Y4": {"exact_json", "chirality", "logged_curves", "curves", "legs", "cert_scale_divisor", "rot_sign", "formula_B"},
    "Y5": {"logged_curves", "summary", "exact_json", "chirality"},
    "Y6": {"readme", "report", "addendum", "correction", "summary", "chirality", "exact_json", "ledger"},
    "Y7": {"scans", "summary"},
    "Y8": {"snapshot", "report"},
    "Y9": {"summary", "curves", "M", "legs", "cert_scale_divisor", "rot_sign"},
}


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

def findings(w):
    out = {"handedness": {}, "family_search": None, "arcs_closest_reading": {}, "u_surface": surface_probe(w)}
    th = grid(w)
    c = container_c(w, th)
    for name, span in (("semicircle", math.pi), ("arc_270", 1.5 * math.pi)):
        b = -span / 2 + span * np.arange(2049) / 2048          # the unit arc, centred on the x-axis
        p = (1 / span) * np.column_stack([np.cos(b), np.sin(b)])
        mine = [lp_dual(p, w["rot_sign"] * math.pi * k / 180, th, c)[0] for k in range(360)]
        logged = [r["scale"] for r in w["scans"][name]]
        out["arcs_closest_reading"][name] = {"reading": "the unit arc centred on the x-axis, 2,048 chords",
                                             "max_difference_over_360_phases": max(abs(a - b) for a, b in zip(mine, logged))}
    for name, p in list(w["curves"].items()) + [(k, v) for k, v in w["logged_curves"].items()]:
        plus, at_p = exact_min(w, p)
        minus, at_m = exact_min(w, mirror(p))
        out["handedness"][name] = {"exact_orientation_preserving": plus, "exact_mirror": minus,
                                   "congruent": min(plus, minus), "chirality_gap": abs(plus - minus),
                                   "grid_256_scan_min": min(scan_rows(w, name)) if name in w["curves"] else None}
    # this tree's own search of the three-link family for the largest congruent critical scale
    def neg_cong(x):
        l1, l3, a1, a3 = x
        l2 = 1 - l1 - l3
        if min(l1, l2, l3) <= 0.01 or not (0.1 < a1 < math.pi - 0.1 and 0.1 < a3 < math.pi - 0.1):
            return 0.0
        p = three_link(l1, l2, l3, a1, a3)
        return -min(exact_min(w, p, scan=360)[0], exact_min(w, mirror(p), scan=360)[0])
    prm = w["chirality"]["corrected_three_link_search"]["parameters"]
    starts = [(prm["l1"], prm["l3"], prm["alpha1_radians"], prm["alpha3_radians"]),
              (1 / 3, 1 / 3, math.pi / 2, math.pi / 2), (0.3, 0.3, math.radians(80), math.radians(80)),
              (0.35, 0.3, math.radians(85), math.radians(75))]
    best = None
    for x0 in starts:
        r = minimize(neg_cong, x0, method="Nelder-Mead", options={"xatol": 1e-9, "fatol": 1e-12, "maxiter": 400})
        if best is None or r.fun < best[0]:
            best = (float(r.fun), [float(v) for v in r.x])
    l1, l3, a1, a3 = best[1]
    p = three_link(l1, 1 - l1 - l3, l3, a1, a3)
    plus, minus = exact_min(w, p)[0], exact_min(w, mirror(p))[0]
    out["family_search"] = {"best_congruent": min(plus, minus), "orientation_preserving": plus, "mirror": minus,
                            "l1": l1, "l2": 1 - l1 - l3, "l3": l3,
                            "alpha1_degrees": math.degrees(a1), "alpha3_degrees": math.degrees(a3),
                            "round_s_cong": w["chirality"]["corrected_three_link_search"]["congruent_scale"]}
    # the search ends at equal thirds with equal turns, where the maximum is a sharp kink: refine the turn there
    def thirds(alpha_deg):
        q = three_link(1 / 3, 1 / 3, 1 / 3, math.radians(alpha_deg), math.radians(alpha_deg))
        return min(exact_min(w, q)[0], exact_min(w, mirror(q))[0])
    v, a = golden(lambda x: -thirds(x), 82.40, 82.47, iters=60)
    out["family_search"]["equal_thirds_refined"] = {"alpha_degrees": a, "congruent": -v,
                                                    "above_round": -v - out["family_search"]["round_s_cong"],
                                                    "slope_per_degree_below": (-v - thirds(a - 2e-4)) / 2e-4,
                                                    "slope_per_degree_above": (-v - thirds(a + 2e-4)) / 2e-4}
    out.update(more_findings(w))
    return out


def more_findings(w):
    """what the report's text rests on beyond the checks: every reading tried, with its numbers"""
    from run002_round01_support_skew_lp import global_min as brent_min    # RUN-002's refinement, for comparison
    th = grid(w)
    c = container_c(w, th)
    tg = w["summary"]["topology_guided_search"]
    u_min, u_at = grid_min(w, w["curves"]["U_equal"])
    b_min, b_at = grid_min(w, w["logged_curves"]["best_topology_guided_curve"])
    out = {"stated_vs_grid_minimum": {
        "U": {"stated_fine": tg["u_equal_fine_scale"], "scan_minimum": w["summary"]["phase_scan"]["curves"]["U_equal"]["minimum_scale"],
              "grid_minimum": u_min, "stated_above": tg["u_equal_fine_scale"] - u_min},
        "best_curve": {"stated": tg["best_fine_scale"], "grid_minimum": b_min, "at_phase": b_at,
                       "recorded_phase": w["ledger"]["phase"], "stated_above": tg["best_fine_scale"] - b_min}}}
    # the arcs: every reading tried, and its largest difference from the logged scan
    readings = {}
    for name, span in (("semicircle", math.pi), ("arc_270", 1.5 * math.pi)):
        logged = [r["scale"] for r in w["scans"][name]]
        def diff(p, phases=range(360)):
            return max(abs(lp_dual(p, w["rot_sign"] * math.pi * k / 180, th, c)[0] - logged[k]) for k in phases)
        def arc(n, radius, offset=0.0):
            b = -span / 2 + offset + span * np.arange(n + 1) / n
            return radius * np.column_stack([np.cos(b), np.sin(b)])
        rows = {}
        for n in (64, 128, 256):
            rows[f"{n} chords of total length 1"] = diff(arc(n, 1 / (2 * n * math.sin(span / (2 * n)))))
            rows[f"{n} chords inscribed in the unit arc"] = diff(arc(n, 1 / span))
        rows["the unit arc (2,048 chords)"] = diff(arc(2048, 1 / span))
        every10 = range(0, 360, 10)
        rot = minimize_scalar(lambda o: diff(arc(2048, 1 / span, o), every10), bounds=(-0.06, 0.06), method="bounded",
                              options={"xatol": 1e-6})
        rad = minimize_scalar(lambda r: diff(arc(2048, r / span), every10), bounds=(0.9, 1.1), method="bounded",
                              options={"xatol": 1e-9})
        rows["the unit arc turned by a fitted angle (fitted on every 10th phase)"] = {"degrees": math.degrees(rot.x), "max_difference": rot.fun}
        rows["the unit arc with a fitted radius (fitted on every 10th phase)"] = {"radius_ratio": rad.x, "max_difference": rad.fun}
        readings[name] = rows
    out["arc_readings"] = readings
    # the surface's corners on the round's own 256-direction grid, 36 and 360 phases
    S = np.array(w["surface"])
    corners = {}
    for i, j in ((0, 0), (0, S.shape[1] - 1), (S.shape[0] - 1, 0), (S.shape[0] - 1, S.shape[1] - 1)):
        p = u_family(math.radians(65 + 50 * j / (S.shape[1] - 1)), 0.18 + 0.24 * i / (S.shape[0] - 1))
        corners[f"{i},{j}"] = {str(n): min(lp_dual(p, 2 * math.pi * k / n, th, c)[0] for k in range(n)) - float(S[i, j])
                               for n in (36, 360)}
    out["u_surface_corners_on_256"] = corners
    # the jumps, read as changes of the dual's active set from one logged phase to the next
    jumps = {}
    for name in REBUILT:
        act = [tuple(np.flatnonzero(lp_dual(w["curves"][name], w["rot_sign"] * math.pi * k / 180, th, c)[2] > 1e-9))
               for k in range(360)]
        jumps[name] = {"active_set_changes": sum(1 for k in range(360) if act[k] != act[k - 1]),
                       "stated_jump_count": w["summary"]["phase_scan"]["curves"][name]["jump_count"]}
    out["jumps_as_active_set_changes"] = jumps
    # RUN-002's refinement (SciPy's bounded Brent) against golden section, on the twelve handedness minima
    gaps = {}
    for name, p in list(w["curves"].items()) + list(w["logged_curves"].items()):
        for label, q in (("+", p), ("-", mirror(p))):
            b = brent_min(lambda x: exact_scale(w, q, x), scan=FINE)[0]
            gaps[f"{name} {label}"] = b - exact_min(w, q)[0]
    out["brent_minus_golden"] = gaps
    return out


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

def drill_cases(base):
    def edit(key, fn):
        v = copy.deepcopy(base[key])
        fn(v)
        return v
    scans_bumped = edit("scans", lambda s: s["U_equal"][200].__setitem__("scale", s["U_equal"][200]["scale"] + 1e-6))
    scans_dropped = edit("scans", lambda s: s["segment"].pop(359))
    ledger_w = edit("ledger", lambda L: L["contact_weights"].__setitem__(1, L["contact_weights"][1] + 1e-5))
    ledger_f = edit("ledger", lambda L: L["container_features"].__setitem__(0, "edge:1-2"))
    summ = edit("summary", lambda s: s["phase_scan"]["curves"]["arc_270"].__setitem__("phase_barrier", 0.1275896627))   # visible at 9 places
    ch_gap = edit("chirality", lambda c: c["false_alarm_candidate"].__setitem__("chirality_gap", 0.054376064834))
    summ_sig = edit("summary", lambda s: s["phase_scan"]["curves"]["segment"].__setitem__("distinct_container_contact_signatures", 5))
    summ_lm = edit("summary", lambda s: s["phase_scan"]["curves"]["semicircle"].__setitem__("local_minimum_count", 4))
    ch = edit("chirality", lambda c: c["corrected_three_link_search"]["reflected"].__setitem__("scale", 0.99857379467))
    ej = edit("exact_json", lambda e: e["optimized_three_link"]["parameters"].__setitem__("l1", 0.2862216772))
    lc = edit("logged_curves", lambda c: c["best_exact_three_link_curve"].__setitem__((2, 0), c["best_exact_three_link_curve"][2, 0] + 1e-4))
    snap = edit("snapshot", lambda s: s["convention"]["facts"].__setitem__("reflection", None))
    rep_tab = base["report"].replace("| semicircle | 0.913857419 |", "| semicircle | 0.913857418 |")
    rep_ref = base["report"].replace("Lower Bound for Convex Hull Area", "Lower Bounds for Convex Hull Area")
    readme = base["readme"].replace("最佳三連桿臨界尺度：0.9933108210", "最佳三連桿臨界尺度：0.9933108211")
    for a, b in ((rep_tab, base["report"]), (rep_ref, base["report"]), (readme, base["readme"])):
        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-scans-copied", "control", {"scans": copy.deepcopy(base["scans"])}, 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)}, {"Y1", "Y2", "Y4"}),
        ("M2-grid-255", "mutation", {"M": 255}, {"Y1", "Y2", "Y9"}),   # 180 and 270 degrees fall off the grid
        ("M3-scan-row", "mutation", {"scans": scans_bumped}, {"Y2"}),
        ("M4-scan-row-dropped", "mutation", {"scans": scans_dropped}, {"Y2", "Y7"}),
        ("M5-ledger-weight", "mutation", {"ledger": ledger_w}, {"Y1", "Y6"}),
        ("M6-ledger-feature", "mutation", {"ledger": ledger_f}, {"Y1", "Y6"}),
        ("M7-barrier", "mutation", {"summary": summ}, {"Y2", "Y6"}),
        ("M8-local-minima", "mutation", {"summary": summ_lm}, {"Y6", "Y7"}),
        ("M9-mirror-scale", "mutation", {"chirality": ch}, {"Y4"}),
        ("M10-parameter", "mutation", {"exact_json": ej}, {"Y5"}),
        ("M11-curve-point", "mutation", {"logged_curves": lc}, {"Y4", "Y5"}),
        ("M13-convention", "mutation", {"snapshot": snap}, {"Y8"}),
        ("M14-table-digit", "mutation", {"report": rep_tab}, {"Y6"}),
        ("M15-reference-title", "mutation", {"report": rep_ref}, {"Y8"}),
        ("M16-readme-digit", "mutation", {"readme": readme}, {"Y6"}),
        ("M17-rotation-sign", "mutation", {"rot_sign": -1}, {"Y2"}),
        ("M18-formula-term", "mutation", {"formula_B": lambda B: 2 * B}, {"Y3", "Y4"}),
        ("M19-signature-count", "mutation", {"summary": summ_sig}, {"Y6", "Y9"}),
        ("M20-chirality-gap", "mutation", {"chirality": ch_gap}, {"Y4", "Y6"}),
    ]


def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("--fetch", action="store_true")
    ap.add_argument("--drill", action="store_true")
    ap.add_argument("--blind", metavar="CID")
    ap.add_argument("--log", default=None, help="write the gate log here instead (a second, cold run to compare)")
    args = ap.parse_args()
    if args.fetch:
        SNAPSHOT.parent.mkdir(parents=True, exist_ok=True)
        snap = {"records": {aid: gate.fetch_arxiv(aid) for aid in CITED},
                "convention": gate.fetch_arxiv_facts(CITED[0], 1, {
                    "congruent": r"contain a congruent copy of every unit-length arc",
                    # the convention itself, not any sentence that happens to say "reflected"
                    "reflection": r"Other orientations are obtained by rotations and reflections, in the congruent-copy "
                                  r"convention used for the cover"})}
        SNAPSHOT.write_text(json.dumps(snap, ensure_ascii=False, indent=1) + "\n", encoding="utf-8", newline="\n")
        print("snapshot written")
        return
    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[:8]:
            print("    ", x)
    log = {
        "run": "RUN-003",
        "subject": "Round 02 (v0.2.2): the contact ledger, phase scans, the three-link search, the exact formula, chirality",
        "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 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()