"""RUN-002 - Round 01 (v0.1): support-skew linear programming, the finite curve library and the first
adversarial search.

Round 01 places a curve in a convex container by linear programming over a grid of support directions:
for a fixed rotation phi, the least scale s with h_gamma(theta_j - phi) + t . u(theta_j) <= s h_C(theta_j) at
every grid direction, and, for a fixed s, the least worst positive "skew". It tests a library of curves in three
containers - the disc of radius 1/2, the 30-degree sector of radius 1 and Wetzel's 30-60-90 triangle T scaled by
1/1.0048 (certified, by a cited preprint, to cover every unit arc) - and runs a five-generation adversarial search.
Everything is recomputed here with this tree's own linear programs; the package's program is not run.

Checks (each returns a list of failures; empty = pass)
  W1 containers      the three containers' areas as the report gives them; T's legs, right angle, and the inscribed
                     square of side 1/3 on its hypotenuse that defines it; the certified scale 1/1.0048
  W2 the two LPs     the report states both programs; on its 192-direction grid the scale program is a relaxation,
                     and the unit segment in the disc comes out at exactly cos(pi/192), the value the report prints
                     as the disc's and the sector's library maximum
  W3 records         every record of a curve the documents define is this tree's scanned minimum over rotations of
                     the same grid program, to 5e-7
  W4 the candidate   the adversarial curve has length 1; at its recorded rotation the grid program gives its recorded
                     scale with 3 tight directions; no rotation gives less; the 98 % test scale, the residual skew E
                     and the logged skew field all reproduce
  W5 tables          the report's top-ten table for the triangle and its three library maxima are the records'
  W6 count           the records hold as many curves as the documents say were tested
  W7 Fourier         translation adds exactly the first harmonic; rotation shifts the support signal's phase
  W8 citation        the cited preprint (arXiv:2606.14625) states T/1.0048, 599 cone models and interval certificates
  W9 readings        the names V_60, V_90, V_120, semicircle and arc_270 are not defined in the documents; exactly one
                     natural reading of each reproduces its records in all three containers

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

import numpy as np
import sympy as sp
from scipy.optimize import linprog, minimize_scalar

from corpus import TREE, members, sha256, source_bytes
import gate

RID = "round-v0.1"
LOG = TREE / "data" / "gate-logs" / "run002-round01-support-skew-lp.json"
SNAPSHOT = TREE / "data" / "external" / "round01-citations.json"
ARXIV = "2606.14625"
TOL = 5e-7                  # a record further above the scanned minimum: the rotation search stopped short
BELOW = 2e-7                # records carry LP residuals up to 1.2e-7, so one may sit this far below the minimum
READING_TOL = 1e-5          # which reading of a curve's name the records use: the right one is within 4.5e-6
                            # of every record, the other misses some record by 1.2e-5 or more
LP_OPTIONS = {"primal_feasibility_tolerance": 1e-10, "dual_feasibility_tolerance": 1e-10}
SCAN = 720                  # rotation scan per (curve, container); each local minimum is then refined
DENSE = 8192                # directions for the sector's "true" value (relaxation error below 1e-7)

CONTAINERS = ("circle_R0.5", "sector_30deg_R1", "wetzel_triangle_certified")
DEFINED = ("segment", "V_60", "V_90", "V_120", "U_equal", "semicircle", "arc_270")


# ---------------------------------------------------------------- geometry

def unit(th):
    th = np.asarray(th, dtype=float)
    return np.stack([np.cos(th), np.sin(th)], -1)


def rotate(p, phi):
    c, s = math.cos(phi), math.sin(phi)
    return p @ np.array([[c, s], [-s, c]])


def support(p, th):
    return (p @ unit(th).T).max(0)


def triangle(w):
    """T/1.0048 with its right angle at the origin; the round's own placement has the short leg a on the x-axis
    (found by W4: only that handedness reproduces the candidate's recorded scale and translation)"""
    a, b = w["legs"]
    legs = (a, b) if w["tri_orientation"] == "a-on-x" else (b, a)
    return np.array([[0.0, 0.0], [legs[0], 0.0], [0.0, legs[1]]]) / w["cert_scale_divisor"]


def container_support(w, name, th):
    th = np.asarray(th, dtype=float)
    if name == "circle_R0.5":
        return np.full(th.shape, 0.5)
    if name == "sector_30deg_R1":
        half = math.radians(15)
        d = (th + math.pi) % (2 * math.pi) - math.pi
        edge = np.maximum(np.cos(d - half), np.cos(d + half))
        return np.where(np.abs(d) <= half, 1.0, np.maximum(edge, 0.0))
    if name == "wetzel_triangle_certified":
        return support(triangle(w), th)
    raise KeyError(name)


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


def lp_scale(p, phi, th, c):
    """least s and its translation with h_p(theta - phi) + t.u(theta) <= s c(theta) on the directions th"""
    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
    # the scale is evaluated exactly at the solver's translation, so the value is attained and the solver's
    # feasibility tolerance never shows up as a smaller scale
    lhs = h + u @ np.array([tx, ty])
    pos = c > 1e-12
    s = float((lhs[pos] / c[pos]).max())
    # a direction where the container's support is 0 (the sector's back) carries no scale: it must hold outright
    back = float(lhs[~pos].max()) if (~pos).any() else -math.inf
    assert back <= 1e-9, f"translation leaves the container's back by {back}"
    # tight: a scaled direction at the scale, or a scale-free one (the triangle's legs through its right angle,
    # placed at the origin, have support 0) met with equality
    tight = int((s * c[pos] - lhs[pos] < 1e-7).sum()) + int((lhs[~pos] > -1e-7).sum())
    return s, (tx, ty), tight


def lp_skew(p, phi, th, c, s):
    """least worst skew e (>= 0) of h_p(theta - phi) + t.u(theta) - s c(theta), and the skew field there"""
    h = support(rotate(p, phi), th)
    u = unit(th)
    res = linprog([0, 0, 1], A_ub=np.column_stack([u, -np.ones_like(c)]), b_ub=s * c - h,
                  bounds=[(None, None), (None, None), (0, None)], method="highs-ds", options=LP_OPTIONS)
    assert res.status == 0, res.message
    tx, ty, e = res.x
    return e, h + u @ np.array([tx, ty]) - s * c


def global_min(f, scan=SCAN, starts=()):
    """min over phi in [0, 2pi) of f: a scan, then each local minimum - and each given start - refined in its
    bracket"""
    phis = 2 * math.pi * np.arange(scan) / scan
    vals = np.array([f(x) for x in phis])
    best = (float(vals.min()), float(phis[vals.argmin()]))
    step = 2 * math.pi / scan
    centres = [phis[i] for i in range(scan) if vals[i] <= vals[i - 1] and vals[i] <= vals[(i + 1) % scan]]
    for x0 in centres + list(starts):
        r = minimize_scalar(f, bounds=(x0 - step, x0 + step), method="bounded", options={"xatol": 1e-12})
        if r.fun < best[0]:
            best = (float(r.fun), float(r.x % (2 * math.pi)))
    return best


# ---------------------------------------------------------------- the curves the documents name

def polyline(*pts):
    return np.array(pts, dtype=float)


def v_curve(angle_deg, reading):
    """two halves of a unit segment; the name's angle read as the turn between them, or as the angle at the joint"""
    turn = math.radians(angle_deg if reading == "turn" else 180 - angle_deg)
    return polyline((0, 0), (0.5, 0), (0.5 + 0.5 * math.cos(turn), 0.5 * math.sin(turn)))


def arc(span, reading, n=64):
    """an arc of total length 1: exact (dense samples), or n chords whose total length is 1"""
    if reading == "exact":
        r, k = 1 / span, 4096
    else:
        r, k = 1 / (2 * n * math.sin(span / (2 * n))), n
    a = span * np.arange(k + 1) / k
    return r * np.column_stack([np.cos(a), np.sin(a)])


def curve(name, readings):
    if name == "segment":
        return polyline((0, 0), (1, 0))
    if name.startswith("V_"):
        return v_curve(int(name[2:]), readings["V"])
    if name == "U_equal":
        return polyline((0, 0), (1 / 3, 0), (1 / 3, 1 / 3), (0, 1 / 3))
    if name == "semicircle":
        return arc(math.pi, readings["arc"])
    if name == "arc_270":
        return arc(1.5 * math.pi, readings["arc"])
    raise KeyError(name)


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

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


def build_world():
    f = load_package()
    rows = lambda name: list(csv.DictReader(io.StringIO(f[name].decode("utf-8"))))
    snap = json.loads(SNAPSHOT.read_text(encoding="utf-8")) if SNAPSHOT.exists() else None
    return {
        "report": f["reports/ROUND_01_REPORT.md"].decode("utf-8"),
        "readme": f["README.md"].decode("utf-8"),
        "summary": json.loads(f["data/experiment_summary.json"]),
        "records": json.loads(f["data/benchmark_records.json"]),
        "history": json.loads(f["data/adversarial_history.json"]),
        "best_curve": np.array([[float(r["x"]), float(r["y"])] for r in rows("data/adversarial_best_curve.csv")]),
        "skew_field": [(float(r["theta"]), float(r["skew"])) for r in rows("data/adversarial_skew_field.csv")],
        "snapshot": snap,
        "M": 192,
        "legs": ((3 + 4 * math.sqrt(3)) / 18, (4 + math.sqrt(3)) / 6),
        "cert_scale_divisor": 1.0048,
        "readings": {"V": "turn", "arc": "chords"},
        "rot_sign": 1,
        "tri_orientation": "a-on-x",
    }


_memo = {}


def program_scale(w, name, cont, readings=None):
    """this tree's minimum over rotations of the round's grid program: scanned, every local minimum refined"""
    readings = readings or w["readings"]
    key = (name, cont, w["M"], w["legs"], w["cert_scale_divisor"], w["tri_orientation"],
           tuple(sorted(readings.items())))
    if key not in _memo:
        p, th = curve(name, readings), grid(w)
        c = container_support(w, cont, th)
        _memo[key] = global_min(lambda phi: lp_scale(p, phi, th, c)[0])[0]
    return _memo[key]


def candidate_min(w, p):
    """the scanned minimum over rotations for a chiral curve, whose basins are narrow: a scan ten times finer,
    and the recorded rotation refined as well"""
    key = ("candidate", p.tobytes(), w["M"], w["legs"], w["cert_scale_divisor"], w["tri_orientation"])
    if key not in _memo:
        th = grid(w)
        c = container_support(w, "wetzel_triangle_certified", th)
        rec = w["rot_sign"] * w["summary"]["adversarial_search"]["best_rotation"]
        _memo[key] = global_min(lambda x: lp_scale(p, x, th, c)[0], scan=10 * SCAN, starts=(rec,))
    return _memo[key]


def records_by(w):
    out = {}
    for r in w["records"]:
        out.setdefault(r["curve"], {})[r["container"]] = r
    return out


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

def w1_containers(w):
    failures = []
    a, b = w["legs"]
    areas = {"circle_R0.5": math.pi / 4, "sector_30deg_R1": math.pi / 12,
             "wetzel_triangle_certified": a * b / 2 / w["cert_scale_divisor"] ** 2}
    for label, cont in (("半徑 $1/2$ 圓盤", "circle_R0.5"), ("半徑 $1$、$30^\\circ$ 扇形", "sector_30deg_R1"),
                        ("Wetzel 三角形的 $1/1.0048$ 縮放版", "wetzel_triangle_certified")):
        m = re.search(re.escape(label) + r"\s*\|\s*([\d.]+)", w["report"])
        if not m or abs(float(m.group(1)) - areas[cont]) > 5e-10:
            failures.append(f"{cont}: report area {m.group(1) if m else None} vs {areas[cont]:.9f}")
        if abs(w["summary"]["container_areas"][cont] - areas[cont]) > 1e-15:
            failures.append(f"{cont}: summary area {w['summary']['container_areas'][cont]} vs {areas[cont]!r}")
    wt = w["summary"]["wetzel_triangle"]
    if abs(wt["a"] - a) > 1e-15 or abs(wt["b"] - b) > 1e-15:
        failures.append(f"summary legs {wt['a']}, {wt['b']} vs ({a!r}, {b!r})")
    if abs(wt["certified_scale"] - 1 / w["cert_scale_divisor"]) > 1e-15:
        failures.append(f"summary certified_scale {wt['certified_scale']} vs 1/{w['cert_scale_divisor']}")
    # T: the right triangle with legs a, b = sqrt(3) a (a 30-60-90 triangle) whose inscribed square on the
    # hypotenuse has side exactly 1/3 - symbolically, from the report's closed forms
    A, B = (3 + 4 * sp.sqrt(3)) / 18, (4 + sp.sqrt(3)) / 6
    c = sp.sqrt(A ** 2 + B ** 2)
    height = A * B / c
    square = sp.nsimplify(sp.simplify(c * height / (c + height)))
    if sp.simplify(B - sp.sqrt(3) * A) != 0:
        failures.append("b is not sqrt(3) a: T is not a 30-60-90 triangle")
    if square != sp.Rational(1, 3):
        failures.append(f"T's inscribed square on the hypotenuse has side {square}, not 1/3")
    if abs(float(A) - a) > 1e-15 or abs(float(B) - b) > 1e-15:
        failures.append("the legs in use differ from the report's closed forms")
    return failures


W2_ANCHORS = ["\\min_{t_x,t_y,s\\ge0}s", "h_\\gamma(\\theta_j-\\phi)+t\\cdotu_{\\theta_j}\\lesh_C(\\theta_j).",
              "\\min_{t_x,t_y,e\\ge0}e", "-sh_C(\\theta_j)\\lee.", "192個支撐方向"]


def w2_programs(w):
    failures = [f"anchor not found: {a}" for a in W2_ANCHORS if a not in gate.squash(w["report"])]
    seg = program_scale(w, "segment", "circle_R0.5")
    expect = math.cos(math.pi / w["M"])
    if abs(seg - expect) > 1e-9:
        failures.append(f"the unit segment in the disc on the {w['M']}-grid gives {seg!r}, not cos(pi/M) = {expect!r}")
    for key in ("circle", "sector"):
        m = re.search(r"s_\{\\max\}\^\{\\mathrm\{" + key + r"\}\}\s*=\s*(\d+\.\d+)", w["report"])
        if not m:
            failures.append(f"report: no {key} maximum")
        elif abs(float(m.group(1)) - expect) > 1e-8:
            failures.append(f"report's {key} maximum {m.group(1)} is not the grid value cos(pi/{w['M']}) = {expect:.8f}")
    return failures


def w3_records(w):
    failures = []
    recs = records_by(w)
    for name in DEFINED:
        for cont in CONTAINERS:
            r = recs.get(name, {}).get(cont)
            if r is None:
                failures.append(f"{name} / {cont}: no record")
                continue
            mine = program_scale(w, name, cont)
            if r["scale"] > mine + TOL:
                # numbers to 9 decimals: float noise in the last digits must not read as a new finding
                failures.append(f"{name} / {cont}: record {r['scale']:.9f} exceeds the scanned minimum {mine:.9f} - "
                                f"the rotation search stopped above it")
            elif r["scale"] < mine - BELOW:
                failures.append(f"{name} / {cont}: record {r['scale']:.9f} below the scanned minimum {mine:.9f}")
    return failures


def w4_candidate(w):
    failures = []
    p, s_ = w["best_curve"], w["summary"]["adversarial_search"]
    length = float(np.linalg.norm(np.diff(p, axis=0), axis=1).sum())
    if abs(length - 1) > 1e-12:
        failures.append(f"the candidate's length is {length!r}, not 1")
    th = grid(w)
    c = container_support(w, "wetzel_triangle_certified", th)
    phi = w["rot_sign"] * s_["best_rotation"]
    s, t, tight = lp_scale(p, phi, th, c)
    if abs(s - s_["best_critical_scale"]) > 1e-9:
        failures.append(f"at the recorded rotation the grid program gives {s!r}, recorded {s_['best_critical_scale']!r}")
    if tight != s_["active_angle_count"]:
        failures.append(f"{tight} tight directions at the recorded placement, recorded {s_['active_angle_count']}")
    if max(abs(a - b) for a, b in zip(t, s_["best_translation"])) > 1e-9:
        failures.append(f"the grid program's translation {tuple(round(x, 12) for x in t)} is not the recorded "
                        f"{tuple(s_['best_translation'])}")
    low, at = candidate_min(w, p)
    if low < s_["best_critical_scale"] - TOL:
        failures.append(f"rotation {at:.6f} gives {low:.9f}, below the recorded critical scale")
    if abs(s_["under_scale_test"] - 0.98 * s_["best_critical_scale"]) > 1e-15:
        failures.append(f"the test scale {s_['under_scale_test']!r} is not 98 % of {s_['best_critical_scale']!r}")
    e, field = lp_skew(p, phi, th, c, s_["under_scale_test"])
    if abs(e - s_["positive_defect"]) > 1e-9:
        failures.append(f"residual skew at the test scale {e!r}, recorded {s_['positive_defect']!r}")
    logged = w["skew_field"]
    if len(logged) != len(th) or any(abs(a - b) > 1e-12 for a, b in zip((x for x, _ in logged), th)):
        failures.append("the logged skew field is not on the program's grid")
    elif max(abs(a - b) for a, b in zip((y for _, y in logged), field)) > 1e-9:
        failures.append(f"the logged skew field differs from the recomputed one by "
                        f"{max(abs(a - b) for a, b in zip((y for _, y in logged), field)):.3e}")
    return failures


TABLE_ROW = re.compile(r"^\|\s*(\d+)\s*\|\s*(\S+)\s*\|\s*(\d+)\s*\|\s*([\d.]+)\s*\|\s*(\d+)\s*\|\s*$", re.M)


def w5_tables(w):
    failures = []
    tri = sorted((r for r in w["records"] if r["container"] == "wetzel_triangle_certified"),
                 key=lambda r: -r["scale"])
    rows = TABLE_ROW.findall(w["report"])
    if len(rows) != 10:
        failures.append(f"the report's triangle table has {len(rows)} rows, not 10")
    for (rank, name, segs, scale, act), r in zip(rows, tri):
        want = (r["curve"], r["segments"], f"{r['scale']:.8f}", len(r["active_angles"]))
        got = (name, int(segs), scale, int(act))
        if got != want:
            failures.append(f"rank {rank}: report {got} vs records {want}")
    for key, cont in (("circle", "circle_R0.5"), ("sector", "sector_30deg_R1"), ("Wetzel", "wetzel_triangle_certified")):
        m = re.search(r"s_\{\\max\}\^\{\\mathrm\{" + key + r"\}\}\s*=\s*(\d+\.\d+)", w["report"])
        top = max(r["scale"] for r in w["records"] if r["container"] == cont)
        if not m or m.group(1) != f"{top:.8f}":
            failures.append(f"report's {key} maximum {m.group(1) if m else None} vs records {top:.8f}")
    return failures


def w6_count(w):
    n = len({r["curve"] for r in w["records"]})
    stated = {"report": re.findall(r"共測試\s*(\d+)\s*條", w["report"]),
              "README": re.findall(r"測試\s*(\d+)\s*條", w["readme"]),
              "summary": [str(w["summary"]["benchmark_curve_count"])]}
    return [f"{where} says {v} curves; the records hold {n}" for where, vals in stated.items()
            for v in vals if int(v) != n]


def w7_fourier(w):
    failures = []
    th, k, tx, ty = sp.symbols("theta k t_x t_y", real=True)
    term = tx * sp.cos(th) + ty * sp.sin(th)
    for n in range(0, 7):
        for basis in (sp.cos(n * th), sp.sin(n * th)):
            val = sp.simplify(sp.integrate(term * basis, (th, 0, 2 * sp.pi)))
            if (n == 1) != (val != 0):
                failures.append(f"translation term against {basis}: {val}")
    p, grid_ = w["best_curve"], grid(w)
    for phi in (0.3, 1.1024251219368202, 4.0):
        err = np.abs(support(rotate(p, phi), grid_) - support(p, grid_ - w["rot_sign"] * phi)).max()
        if err > 1e-14:
            failures.append(f"rotation by {phi} is not a phase shift of the support signal (error {err:.2e})")
    return failures


def w8_citation(w):
    snap = w["snapshot"]
    if snap is None:
        return ["no citation snapshot - run with --fetch"]
    failures = []
    rec, abs_ = snap["record"], snap["abstract"]["abstract"]
    if "Wetzel" not in rec["title"] or "30-60-90" not in rec["title"]:
        failures.append(f"title: {rec['title']}")
    for needle in ("1.0048", "599", "interval", "0.260956"):
        if needle not in abs_:
            failures.append(f"the preprint's abstract does not state {needle}")
    a, b = w["legs"]
    area = a * b / 2 / w["cert_scale_divisor"] ** 2
    m = re.search(r"0\.260956\d*", abs_)
    if m and not f"{area:.12f}".startswith(m.group(0)):
        failures.append(f"the preprint's area {m.group(0)} vs T/1.0048 here {area:.12f}")
    if "599" not in w["report"] or "SOCP" not in w["report"]:
        failures.append("the report does not cite the 599 cone models")
    if rec["years"] != [2026]:
        failures.append(f"first version year {rec['years']}")
    return failures


def w9_readings(w):
    failures = []
    recs = records_by(w)
    for family, names, options in (("V", ("V_60", "V_90", "V_120"), ("turn", "joint")),
                                   ("arc", ("semicircle", "arc_270"), ("chords", "exact"))):
        fits = []
        for opt in options:
            readings = dict(w["readings"], **{family: opt})
            # the records stop a few millionths above their minima (W3), the wrong reading misses by hundredths
            ok = all(abs(recs[n][c]["scale"] - program_scale(w, n, c, readings)) <= READING_TOL
                     for n in names for c in CONTAINERS)
            if ok:
                fits.append(opt)
        if fits != [w["readings"][family]]:
            failures.append(f"{family}: readings that reproduce every record: {fits}; in use: {w['readings'][family]}")
    return failures


CHECKS = {"W1": w1_containers, "W2": w2_programs, "W3": w3_records, "W4": w4_candidate, "W5": w5_tables,
          "W6": w6_count, "W7": w7_fourier, "W8": w8_citation, "W9": w9_readings}
DEPENDS = {
    "W1": {"report", "summary", "legs", "cert_scale_divisor"},
    "W2": {"report", "M", "legs", "cert_scale_divisor", "readings"},
    "W3": {"records", "M", "legs", "cert_scale_divisor", "readings", "tri_orientation"},
    "W4": {"best_curve", "summary", "skew_field", "M", "legs", "cert_scale_divisor", "rot_sign", "tri_orientation"},
    "W5": {"report", "records"},
    "W6": {"report", "readme", "summary", "records"},
    "W7": {"best_curve", "M", "rot_sign"},
    "W8": {"snapshot", "legs", "cert_scale_divisor", "report"},
    "W9": {"records", "M", "legs", "cert_scale_divisor", "readings", "tri_orientation"},
}


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

def true_scales(w):
    """the continuous problem, for the curves the documents define: the disc exactly (twice the least enclosing
    radius), the triangle exactly (its three edge normals), the sector on 8,192 directions plus its edge normals"""
    exact = {"V": w["readings"]["V"], "arc": "exact"}
    out = {}
    v = triangle(w)                                     # counter-clockwise: outward normal of edge d is (dy, -dx)
    d = np.roll(v, -1, axis=0) - v
    tri_normals = np.arctan2(-d[:, 0], d[:, 1])
    dense = np.concatenate([2 * math.pi * np.arange(DENSE) / DENSE, [math.radians(105), math.radians(-105)]])
    for name in DEFINED:
        p = curve(name, exact)
        if name in ("semicircle", "arc_270"):
            span = math.pi if name == "semicircle" else 1.5 * math.pi
            disc = 2 / span   # 2r: a semicircle sits in the circle on its chord, a longer arc in its own circle
        else:
            disc = 2 * enclosing_radius(p)
        tc = container_support(w, "wetzel_triangle_certified", tri_normals)
        tri = global_min(lambda x: lp_scale(p, x, tri_normals, tc)[0])[0]
        th = grid(w)
        sc = container_support(w, "sector_30deg_R1", th)
        _, at = global_min(lambda x: lp_scale(p, x, th, sc)[0])
        dc = container_support(w, "sector_30deg_R1", dense)
        r = minimize_scalar(lambda x: lp_scale(p, x, dense, dc)[0], bounds=(at - 0.02, at + 0.02),
                            method="bounded", options={"xatol": 1e-12})
        out[name] = {"circle_R0.5": disc, "sector_30deg_R1": float(r.fun), "wetzel_triangle_certified": tri}
    return out


def enclosing_radius(p):
    """least enclosing circle of a few points, by checking every pair and triple"""
    pts = [tuple(x) for x in p]
    best = math.inf
    def covers(cx, cy, r):
        return all(math.hypot(x - cx, y - cy) <= r + 1e-12 for x, y in pts)
    for i in range(len(pts)):
        for j in range(i + 1, len(pts)):
            (x1, y1), (x2, y2) = pts[i], pts[j]
            cx, cy, r = (x1 + x2) / 2, (y1 + y2) / 2, math.hypot(x1 - x2, y1 - y2) / 2
            if r < best and covers(cx, cy, r):
                best = r
            for k in range(j + 1, len(pts)):
                x3, y3 = pts[k]
                d = 2 * (x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2))
                if abs(d) < 1e-15:
                    continue
                ux = ((x1 ** 2 + y1 ** 2) * (y2 - y3) + (x2 ** 2 + y2 ** 2) * (y3 - y1) + (x3 ** 2 + y3 ** 2) * (y1 - y2)) / d
                uy = ((x1 ** 2 + y1 ** 2) * (x3 - x2) + (x2 ** 2 + y2 ** 2) * (x1 - x3) + (x3 ** 2 + y3 ** 2) * (x2 - x1)) / d
                r = math.hypot(x1 - ux, y1 - uy)
                if r < best and covers(ux, uy, r):
                    best = r
    return best


def findings(w):
    recs = records_by(w)
    truth = true_scales(w)
    rows = []
    for name in DEFINED:
        for cont in CONTAINERS:
            rec = recs[name][cont]["scale"]
            rows.append({"curve": name, "container": cont, "record": rec,
                         "program_scanned_min": program_scale(w, name, cont),
                         "true": truth[name][cont], "record_minus_true": rec - truth[name][cont]})
    hist = w["history"]
    p = w["best_curve"]
    low, at = candidate_min(w, p)
    mirrored = p * np.array([1.0, -1.0])
    m_low, m_at = candidate_min(w, mirrored)
    return {
        "candidate": {"recorded": w["summary"]["adversarial_search"]["best_critical_scale"],
                      "scanned_min_over_rotations": low, "at_rotation": at,
                      "mirror_image_scanned_min": m_low, "mirror_at_rotation": m_at,
                      "congruent_critical_scale": min(low, m_low)},
        "readings": w["readings"],
        "defined_curves": rows,
        "library_maxima_true": {c: max(truth[n][c] for n in DEFINED) for c in CONTAINERS},
        # the records' arcs are 64-chord polylines of total length 1 (W9); in the disc their own true value is the
        # diameter of the circle their vertices lie on
        "arc_polylines_in_disc": {
            name: {"chord_circle_diameter": 2 / (64 * math.sin(span / 128)),
                   "record": recs[name]["circle_R0.5"]["scale"],
                   "grid": program_scale(w, name, "circle_R0.5"),
                   "record_over_true": recs[name]["circle_R0.5"]["scale"] * 64 * math.sin(span / 128) / 2,
                   "grid_over_true": program_scale(w, name, "circle_R0.5") * 64 * math.sin(span / 128) / 2,
                   "exact_arc": 2 / span}
            for name, span in (("semicircle", math.pi), ("arc_270", 1.5 * math.pi))},
        "not_rebuildable": sorted({r["curve"] for r in w["records"]} - set(DEFINED)),
        "history_vs_summary": {"history_best_by_generation": [h["best_scale"] for h in hist],
                               "summary_best": w["summary"]["adversarial_search"]["best_critical_scale"]},
        "candidate_length": float(np.linalg.norm(np.diff(w["best_curve"], axis=0), axis=1).sum()),
        "grid_relaxation_factor": math.cos(math.pi / w["M"]),
    }


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

def drill_cases(base):
    rec_bump = copy.deepcopy(base["records"])
    for r in rec_bump:
        if r["curve"] == "U_equal" and r["container"] == "sector_30deg_R1":
            r["scale"] += 1e-4
    rec_top = copy.deepcopy(base["records"])
    for r in rec_top:
        if r["curve"] == "random_18_5s" and r["container"] == "wetzel_triangle_certified":   # in the table, not in W3
            r["scale"] -= 2e-8
    rec_drop = [r for r in base["records"] if r["curve"] != "random_23_6s"]
    rec_shuffled = list(reversed(base["records"]))
    sm = copy.deepcopy(base["summary"])
    sm["adversarial_search"]["positive_defect"] *= 1.001
    sm_t = copy.deepcopy(base["summary"])
    sm_t["adversarial_search"]["best_translation"][0] += 1e-6
    sm_area = copy.deepcopy(base["summary"])
    sm_area["container_areas"]["sector_30deg_R1"] += 1e-9
    bc = base["best_curve"].copy()
    bc[4] += np.array([1e-3, 0.0])
    sf = list(base["skew_field"])
    sf[100] = (sf[100][0], sf[100][1] + 1e-6)
    snap = copy.deepcopy(base["snapshot"])
    snap["abstract"]["abstract"] = snap["abstract"]["abstract"].replace("1.0048", "1.0049")
    rep_table = base["report"].replace("| 2 | arc_270 | 64 | 0.92558266 | 3 |", "| 2 | arc_270 | 64 | 0.92558267 | 3 |")
    assert rep_table != base["report"], "the table plant must hit the report's raw text"
    rep_lp = base["report"].replace("s h_C(\\theta_j).", "s h_C(\\theta_{j+1}).", 1)
    rep_count = base["report"].replace("共測試 35 條", "共測試 36 條")
    a, b = base["legs"]
    return [
        ("K0-unmodified", "control", {}, set()),
        ("K1-records-reversed", "control", {"records": rec_shuffled}, set()),
        ("K2-report-crlf", "control", {"report": base["report"].replace("\n", "\r\n")}, set()),
        ("K3-report-trailing-spaces", "control", {"report": "\n".join(l + "  " for l in base["report"].split("\n"))}, set()),
        ("M1-triangle-leg", "mutation", {"legs": (a * 1.001, b)}, {"W1", "W3", "W4", "W8", "W9"}),
        ("M2-certified-divisor", "mutation", {"cert_scale_divisor": 1.0049}, {"W1", "W3", "W4", "W8", "W9"}),
        ("M3-grid-191", "mutation", {"M": 191}, {"W2", "W3", "W4", "W9"}),
        ("M4-record-bumped", "mutation", {"records": rec_bump}, {"W3"}),
        ("M5-record-below-table", "mutation", {"records": rec_top}, {"W5"}),
        ("M6-curve-dropped", "mutation", {"records": rec_drop}, {"W6"}),
        ("M7-defect-misstated", "mutation", {"summary": sm}, {"W4"}),
        ("M8-curve-point-moved", "mutation", {"best_curve": bc}, {"W4"}),
        ("M9-skew-field-entry", "mutation", {"skew_field": sf}, {"W4"}),
        ("M10-citation-scale", "mutation", {"snapshot": snap}, {"W8"}),
        ("M11-table-digit", "mutation", {"report": rep_table}, {"W5"}),
        ("M12-lp-statement", "mutation", {"report": rep_lp}, {"W2"}),
        ("M13-count-stated", "mutation", {"report": rep_count}, {"W6"}),
        ("M14-V-reading", "mutation", {"readings": {"V": "joint", "arc": "chords"}}, {"W3", "W9"}),
        ("M15-arc-reading", "mutation", {"readings": {"V": "turn", "arc": "exact"}}, {"W3", "W9"}),
        ("M16-rotation-sign", "mutation", {"rot_sign": -1}, {"W4", "W7"}),
        ("M17-sector-area", "mutation", {"summary": sm_area}, {"W1"}),
        ("M18-triangle-mirrored", "mutation", {"tri_orientation": "b-on-x"}, {"W4"}),
        ("M19-translation-misstated", "mutation", {"summary": sm_t}, {"W4"}),
    ]


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()
    log_path = pathlib.Path(args.log) if args.log else LOG
    if args.fetch:
        SNAPSHOT.parent.mkdir(parents=True, exist_ok=True)
        snap = {"record": gate.fetch_arxiv(ARXIV), "abstract": gate.fetch_arxiv_abstract(ARXIV)}
        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-002",
        "subject": "Round 01 (v0.1): support-skew linear programming, the curve library, the first adversarial search",
        "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),        # a blind run writes no log
    }
    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:
        log_path.parent.mkdir(parents=True, exist_ok=True)
        log_path.write_text(json.dumps(log, ensure_ascii=False, indent=1) + "\n", encoding="utf-8", newline="\n")


if __name__ == "__main__":
    main()