"""RUN-012 - Round 11 (v0.11): exact contact boundaries, derivative boxes, root boxes, boundary neighbourhoods and the
special 120/270-degree enclosures - a semi-validated certificate.

Round 11 writes Round 10's 17 switches as exact angles (pi/2 - alpha, beta - pi/3, ...), encloses s(3pi/2) and
s(2pi/3) - s(3pi/2) by a composite midpoint rule with a second-derivative bound and a float padding, covers each of the
18 intervals' interiors with 579 adaptive derivative boxes of one sign, isolates 12 smooth stationary points in root
boxes whose second-derivative boxes exclude zero, and treats each switch's neighbourhood of radius 2e-6 by one-sided
derivative intervals or a lower bound of the scale. It calls the whole a semi-validated certificate: analytic error
bounds and float padding, without directed rounding throughout.

Checks (each returns a list of failures; empty = pass)
  U1 boundaries   each of the 19 boundary formulas, at 60 digits, is the stated phase (to 1e-15) and degrees, and
                  satisfies its switch's own condition (to 1e-40); the largest difference from Round 10's switches
  U2 special      the s270, s120 and s120 - s270 intervals contain the exact values (200 digits); the stated error
                  bounds are the composite-midpoint bound M2 / (24 N^2) with the stated constants; the centres are the
                  midpoint sums
  U3 event        the gain interval is the s270 interval less Round 09's event branch B1, exactly; positive
  U4 boxes        579 boxes; each box's derivative interval contains the exact derivative at 13 points of the box (to
                  1e-12) and has the stated sign; per-interval counts and signs; the boxes, root boxes and guards tile
                  each interval
  U5 roots        each root box holds an exact stationary point; its second-derivative interval contains s'' at four
                  points inside it (60 digits, to 1e-10);
                  its scale interval contains s over its scale box; the margins; the boxes are all the stationary
                  points (fixed-text failure)
  U6 boundaries'  each neighbourhood's scale intervals contain s over that side; its one-sided derivative intervals
      neighbour-  contain the derivative over that side, the switch included (60 digits; fixed-text failures, one per
      hoods       boundary); the classifications are the derivative's signs; the safe lower bounds, margins and flags
                  (a neighbourhood with a special certificate is excluded by that certificate, not by its margin)
  U7 copies       the summary and the CSVs hold the data files' figures
  U8 documents    the README, the report's tables and figures, the audit note, the handoff and the manifest
  U9 method       the verifier's constants are Round 09's decimals; its bounds are the stated M1, M2

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

import mpmath as mp
import numpy as np

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

RID = "round-v0.11"
LOG = TREE / "data" / "gate-logs" / "run012-round11-box-certificate.json"
NAMES = ("min_x", "min_y", "max_d")


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

def build_world():
    f, f10 = load(RID), load("round-v0.10")
    md = lambda n: restore_escapes(f[n].decode("utf-8"))
    rows = lambda n: list(csv.DictReader(io.StringIO(f[n].decode("utf-8"))))
    js = lambda n: json.loads(f[n])
    return {
        "readme": md("README.md"), "report": md("reports/ROUND_11_REPORT.md"),
        "audit_note": md("reports/ROUND_11_CERTIFICATE_AUDIT.md"), "handoff": md("reports/ROUND_12_BALL_ARITHMETIC_HANDOFF.md"),
        "method": md("src/METHOD.md"), "verifier": f["src/special_branch_verify.py"].decode("utf-8"),
        "manifest": js("MANIFEST.json"), "members": sorted(f),
        "summary": js("data/round11_summary.json"),
        "boundaries": js("data/exact_contact_boundaries.json"),
        "special": js("data/special_branch_interval_certificate.json"),
        "boxes": js("data/phase_derivative_box_certificate.json"),
        "roots": js("data/stationary_root_boxes.json"),
        "neighbourhoods": js("data/boundary_neighborhood_audit.json"),
        "box_csv": rows("data/derivative_interval_summary.csv"),
        "root_csv": rows("data/root_box_summary.csv"),
        "round10": js("data/round10_summary.json") if "data/round10_summary.json" in f else json.loads(f10["data/round10_summary.json"]),
        "round09": json.loads(load("round-v0.9")["data/round09_summary.json"]),
        "round05": json.loads(load("round-v0.5")["data/round05_summary.json"])["event_solution"]["parameters"],
        "legs": ("1", "1"),
    }


def w11(w):
    """the world RUN-011's tools read: Round 10's map and the curve"""
    return {"summary": w["round10"], "round05": w["round05"], "legs": w["legs"]}


def consts(w):
    D = w["round09"]["definitions"]["candidate"]
    return {k: mp.mpf(D[k]) for k in ("epsilon", "wing_length", "beta", "delta", "transition_center")}


# ---------------------------------------------------------------- the exact support in double precision, vectorized

class Model:
    """the candidate in double precision (RUN-009's wing: positions exact to rounding) with the exact support"""

    def __init__(self, w):
        P = w["round10"]["curve_parameters"]
        self.wing = Wing(P["wing_length"], P["beta"], P["delta"], P["transition_center"], P["epsilon"])
        fa, fb = (float(x) for x in w["legs"])
        self.A = (3 + 4 * math.sqrt(3)) / 18 / 1.0048 * fa
        self.B = (4 + math.sqrt(3)) / 6 / 1.0048 * fb
        x1, y1 = self.wing.P1
        l0 = self.wing.l0
        self.P = {"p0": (0.0, 0.0), "p1": (x1, y1), "p2": (x1 + l0, y1), "p3": (2 * x1 + l0, 0.0)}
        self.a = {"max_d": (1 / self.A, 1 / self.B), "min_x": (-1 / self.A, 0.0), "min_y": (0.0, -1 / self.B)}

    def point(self, label, name, ph):
        if label in self.P:
            return np.full_like(ph, self.P[label][0]), np.full_like(ph, self.P[label][1])
        a0, a1 = self.a[name]
        c, s = np.cos(ph), np.sin(ph)
        n0, n1 = a0 * c + a1 * s, -a0 * s + a1 * c
        ang = np.arctan2(n1, n0) if label == "L" else np.arctan2(n1, -n0)
        W = self.wing
        best = None
        for k in (-1, 1):
            t = np.remainder(ang + k * math.pi / 2 + math.pi, 2 * math.pi) - math.pi
            f = (t - W.b) / W.d
            best = f if best is None else np.where(np.abs(f - 0.5) < np.abs(best - 0.5), f, best)
        f = np.clip(best, 0.0, 1.0)
        u = np.clip(W.c - W.eps * np.arctanh(np.clip(W.T1 + f * (W.T0 - W.T1), -1 + 1e-16, 1 - 1e-16)), 0.0, 1.0)
        p = W.pos(u)
        if label == "L":
            return p[:, 0], p[:, 1]
        return 2 * self.P["p1"][0] + W.l0 - p[:, 0], p[:, 1]

    def deriv(self, signature, ph):
        ph = np.atleast_1d(np.asarray(ph, dtype=float))
        c, s = np.cos(ph), np.sin(ph)
        tot = np.zeros_like(ph)
        for name, lab in zip(NAMES, signature.split("|")):
            a0, a1 = self.a[name]
            px, py = self.point(lab, name, ph)
            tot += (-a0 * s + a1 * c) * px + (-a0 * c - a1 * s) * py
        return tot

    def scale(self, ph):
        """s at phases, the exact support: the corners and each wing's interior point"""
        ph = np.atleast_1d(np.asarray(ph, dtype=float))
        c, s = np.cos(ph), np.sin(ph)
        tot = np.zeros_like(ph)
        for name in NAMES:
            a0, a1 = self.a[name]
            n0, n1 = a0 * c + a1 * s, -a0 * s + a1 * c
            vals = [self.P[k][0] * n0 + self.P[k][1] * n1 for k in self.P]
            for lab in ("L", "R"):
                px, py = self.point(lab, name, ph)
                vals.append(px * n0 + py * n1)
            tot += np.max(np.stack(vals, 1), 1)
        return tot


_models = {}


def model(w):
    key = (json.dumps(w["round10"]["curve_parameters"], sort_keys=True), tuple(w["legs"]))
    if key not in _models:
        _models[key] = Model(w)
    return _models[key]


def stationary(w):
    """this run's stationary points, from RUN-011 (exact derivative, 40 digits)"""
    return r11.all_stationary(w11(w))


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

def alpha_beta(w):
    K = consts(w)
    return K["beta"] + K["delta"], K["beta"]


def formula_value(text, al, be):
    expr = text.replace("\\alpha", "(AL)").replace("\\beta", "(BE)").replace("\\pi", "(PI)").replace("{", "(").replace("}", ")")
    expr = re.sub(r"(\d)\(", r"\1*(", expr)
    return eval(expr, {"__builtins__": {}}, {"AL": al, "BE": be, "PI": mp.pi})


def switch_condition(C, sw, ph):
    """how far a phase is from being the switch: for two corners, their gap under the changing functional; for a
    wing's interior support and a corner, the distance of the wing's perpendicular point from the wing's end (F = 0
    or 1, the tangent at beta or alpha)"""
    L, R = sw["left_signature"].split("|"), sw["right_signature"].split("|")
    k = [i for i in range(3) if L[i] != R[i]][0]
    name, a, b = NAMES[k], L[k], R[k]
    with mp.workdps(60):
        A_ = C.a[name]
        cs, sn = mp.cos(ph), mp.sin(ph)
        n = (A_[0] * cs + A_[1] * sn, -A_[0] * sn + A_[1] * cs)
        if a in C.P and b in C.P:
            return abs(n[0] * (C.P[a][0] - C.P[b][0]) + n[1] * (C.P[a][1] - C.P[b][1]))
        lab = a if a in ("L", "R") else b
        ang = mp.atan2(n[1], n[0]) if lab == "L" else mp.atan2(n[1], -n[0])
        best = mp.inf
        for kk in (-1, 1):
            tt = ang + kk * mp.pi / 2
            tt -= 2 * mp.pi * mp.floor((tt + mp.pi) / (2 * mp.pi))
            f = (tt - C.b) / C.d
            best = min(best, abs(f), abs(f - 1))
        return best


def u1_boundaries(w):
    failures = []
    B = w["boundaries"]
    sw = w["round10"]["contact_switches"]
    with mp.workdps(60):
        al, be = alpha_beta(w)
        vals = []
        for b in B:
            v = formula_value(b["formula"], al, be)
            vals.append(v)
            if abs(mp.mpf(b["phase"]) - v) > 1e-15 or abs(b["degrees"] - float(v * 180 / mp.pi)) > 1e-9:
                failures.append(f"boundary {b['index']}: {b['formula']} is not its phase")
        if len(B) != 19 or B[0]["phase"] != 0 or abs(vals[-1] - 2 * mp.pi) > 1e-40:
            failures.append("the boundaries are not 19 from 0 to 2pi")
        C = r11.cand(w11(w), 60)
        for i, s in enumerate(sw):
            if switch_condition(C, s, vals[i + 1]) > mp.mpf(10) ** -40:
                failures.append(f"boundary {i + 1} is not the switch's condition")
    diff = max(abs(b["phase"] - s["phase"]) for b, s in zip(B[1:-1], sw))
    if abs(diff - w["summary"]["exact_contact_partition"]["maximum_difference_from_round10_numeric_edges"]) > 1e-20:
        failures.append("the largest difference from Round 10's switches is not as stated")
    return failures


def midpoint_sums(w, N):
    """the composite midpoint sums of cos(theta), sin(theta) and q(theta) over N parts, in double precision"""
    P = w["round10"]["curve_parameters"]
    eps, wl, b, d, c = P["epsilon"], P["wing_length"], P["beta"], P["delta"], P["transition_center"]
    raw = lambda u: 0.5 * (1 - np.tanh((u - c) / eps))
    r0, r1 = float(raw(0.0)), float(raw(1.0))
    u = (np.arange(N) + 0.5) / N
    th = b + d * (raw(u) - r1) / (r0 - r1)
    qa, qb = 2 / math.sqrt(3) - 1, 1 / math.sqrt(3)
    Cc = d * 0.5 / eps / (r0 - r1)
    M2 = Cc ** 2 + 2 * Cc / eps * 2 / (3 * math.sqrt(3))
    return {"cos": float(np.sum(np.cos(th)) / N), "sin": float(np.sum(np.sin(th)) / N),
            "q": float(np.sum(qa * np.sin(th) - qb * np.cos(th)) / N), "M1": Cc, "M2": M2, "Mq": math.hypot(qa, qb)}


def u2_special(w):
    failures = []
    H = r11.hp(w11(w))
    s = r10.s270(H)
    with mp.workdps(100):
        s120 = r11.cand(w11(w), 100).s(2 * mp.pi / 3)
        d = s120 - s
    SP = w["special"]
    for label, iv_, true in (("s270", SP["s270_interval"], s), ("s120", SP["s120_interval"], s120),
                             ("s120 - s270", SP["s120_minus_s270_interval"], d)):
        with mp.workdps(100):
            if not mp.mpf(iv_[0]) <= true <= mp.mpf(iv_[1]):
                failures.append(f"the {label} interval does not contain the exact value")
    N = SP["quadrature_partition_count"]
    m = midpoint_sums(w, N)
    if abs(SP["coordinate_midpoint_error_bound"] - (m["M2"] / (24 * N * N) + 2e-15)) > 1e-25 \
            or abs(SP["q_difference_midpoint_error_bound"] - (m["Mq"] * m["M2"] / (24 * N * N) + 2e-15)) > 1e-25:
        failures.append("the error bounds are not M2/(24 N^2) with the stated padding")
    P = w["round10"]["curve_parameters"]
    A, B = r10.legs_mp(H)
    A, B = float(A), float(B)
    l0 = 1 - 2 * P["wing_length"]
    c270 = P["wing_length"] * m["sin"] / A + (P["wing_length"] * m["cos"] + l0) / B
    cdiff = l0 * (1 / (2 * A) - 1 / B) + P["wing_length"] / A * m["q"]
    if abs(sum(SP["s270_interval"]) / 2 - c270) > 1e-14 or abs(sum(SP["s120_minus_s270_interval"]) / 2 - cdiff) > 1e-14:
        failures.append("the intervals are not centred on the midpoint sums")
    return failures


def u3_event(w):
    failures = []
    SP = w["special"]
    B1 = w["round09"]["arbitrary_precision"]["event_branch_values"][0]
    if SP["event_control_branch_B1"] != B1:
        failures.append("the event branch is not Round 09's B1")
    with mp.workdps(130):
        for k in (0, 1):
            # the bound as its shortest decimal, less B1; the stated string cuts it to its own decimals
            want = mp.mpf(repr(SP["s270_interval"][k])) - mp.mpf(B1)
            txt = SP["smooth_minus_event_interval"][k]
            n = len(txt.split(".")[1])
            if not 0 <= want - mp.mpf(txt) < mp.mpf(10) ** -n:
                failures.append(f"the gain interval's end {k} is not the s270 bound less B1")
        if not mp.mpf(SP["smooth_minus_event_interval"][0]) > 0:
            failures.append("the gain's lower end is not positive")
    return failures


def box_samples(w, b):
    a, z = b["start"], b["end"]
    fr = np.array([j / 8 for j in range(9)] + [1e-4, 1e-8, 1 - 1e-4, 1 - 1e-8])
    return model(w).deriv(b["signature"], a + fr * (z - a))


def u4_boxes(w):
    failures = []
    BX = w["boxes"]
    if len(BX["boxes"]) != 579 or BX["total_certified_subboxes"] != len(BX["boxes"]) or BX["total_unresolved_subboxes"] != 0:
        failures.append("the boxes are not 579, all resolved")
    for k, b in enumerate(BX["boxes"]):
        v = box_samples(w, b)
        L, H = b["derivative_interval"]
        if v.min() < L - 1e-12 or v.max() > H + 1e-12:
            failures.append(f"box {k}: the interval does not contain the derivative")
        if (b["sign"] == "negative" and not H < 0) or (b["sign"] == "positive" and not L > 0):
            failures.append(f"box {k}: the sign is not its interval's")
    roots = w["roots"]
    for S in BX["summary"]:
        i = S["interval_index"]
        mine = sorted([(b["start"], b["end"]) for b in BX["boxes"] if b["interval_index"] == i] +
                      [tuple(r["uniqueness_box"]) for r in roots if r["interval_index"] == i])
        n = sum(1 for b in BX["boxes"] if b["interval_index"] == i)
        signs = sorted({b["sign"] for b in BX["boxes"] if b["interval_index"] == i})
        if n != S["certified_subbox_count"] or signs != sorted(S["certified_signs"]):
            failures.append(f"interval {i}: the counts or signs are not its boxes'")
        g = S["guard"]
        ok = abs(mine[0][0] - (S["start"] + g)) < 1e-15 and abs(mine[-1][1] - (S["end"] - g)) < 1e-15 \
            and all(x[1] >= y[0] - 1e-15 for x, y in zip(mine, mine[1:]))
        if not ok:
            failures.append(f"interval {i}: the boxes do not tile it within its guards")
    return failures


def u5_roots(w):
    failures = []
    M = model(w)
    stat = stationary(w)
    s270_hi = w["special"]["s270_interval"][1]
    for j, r in enumerate(w["roots"]):
        a, z = r["uniqueness_box"]
        inside = [e for e in stat if a <= e["phase"] <= z]
        if len(inside) != 1 or inside[0]["interval"] != r["interval_index"] or inside[0]["kind"] != r["kind"]:
            failures.append(f"root box {j} does not hold one exact stationary point of its kind")
            continue
        L, H = r["second_derivative_interval"]
        d2 = root_second(w, j)
        if min(d2) < L - 1e-10 or max(d2) > H + 1e-10 or (L <= 0 <= H):
            failures.append(f"root box {j}: the second-derivative interval does not contain s'' or contains 0")
        sb = np.linspace(r["scale_evaluation_box"][0], r["scale_evaluation_box"][1], 5)
        sv = M.scale(sb)
        if sv.min() < r["scale_interval"][0] - 1e-13 or sv.max() > r["scale_interval"][1] + 1e-13:
            failures.append(f"root box {j}: the scale interval does not contain s")
        if abs(r["scale_lower_minus_s270_upper"] - (r["scale_interval"][0] - s270_hi)) > 1e-18:
            failures.append(f"root box {j}: the margin")
    if len(stat) != len(w["roots"]):
        failures.append("the root boxes are not all the stationary points")
    return failures


def boundary_phase(w, k):
    with mp.workdps(60):
        al, be = alpha_beta(w)
        return formula_value(w["boundaries"][k]["formula"], al, be)


def root_second(w, j):
    """s'' at four points inside root box j (2 % and 98 % of the way, the middle, the root), by a central difference of
    the 60-digit derivative with h = 1e-25"""
    def go():
        r = w["roots"][j]
        C = r11.cand(w11(w), 60)
        stat = stationary(w)
        with mp.workdps(60):
            a, z = mp.mpf(repr(r["uniqueness_box"][0])), mp.mpf(repr(r["uniqueness_box"][1]))
            q = next(e["phase"] for e in stat if a <= e["phase"] <= z)
            h = mp.mpf("1e-25")
            pts = [a + (z - a) * mp.mpf("0.02"), (a + z) / 2, a + (z - a) * mp.mpf("0.98"), mp.mpf(q)]
            return [float((r11.deriv(C, r["signature"], x + h) - r11.deriv(C, r["signature"], x - h)) / (2 * h)) for x in pts]
    return r10.memo("second", r11.hp(w11(w)), j, json.dumps(w["roots"][j]["uniqueness_box"]), w["roots"][j]["signature"])(go)


def side_hp(w, nb, side):
    """the derivative on one side of a switch at 60 digits: at the switch itself (that side's one-sided value), at
    2e-10 rad from it, and at offsets down to 1e-12 of the radius"""
    def go():
        C = r11.cand(w11(w), 60)
        with mp.workdps(60):
            ph = boundary_phase(w, nb["boundary_index"])
            r = mp.mpf(repr(nb["neighborhood_radius"]))
            sgn = -1 if side == "left" else 1
            offs = [mp.mpf(0), mp.mpf("2e-10")] + [r * mp.mpf(f) for f in ("1", "0.5", "0.1", "1e-2", "1e-3", "1e-4", "1e-6", "1e-8", "1e-10", "1e-12")]
            return [(o, r11.deriv(C, nb[f"{side}_signature"], ph + sgn * o)) for o in offs]
    return r10.memo("side", r11.hp(w11(w)), nb["boundary_index"], side, nb[f"{side}_signature"], nb["neighborhood_radius"])(go)


def side_samples(w, nb, side):
    ph, r = nb["phase"], nb["neighborhood_radius"]
    sgn = -1 if side == "left" else 1
    offs = r * np.array([1, 0.5, 0.1, 1e-2, 1e-3, 1e-4, 1e-6, 1e-8, 1e-10, 1e-12])
    M = model(w)
    sig = nb[f"{side}_signature"]
    return M.deriv(sig, ph + sgn * offs), M.scale(ph + sgn * offs)


def u6_neighbourhoods(w):
    failures = []
    s270_hi = w["special"]["s270_interval"][1]
    for nb in w["neighbourhoods"]:
        k = nb["boundary_index"]
        enclose = True
        signs = {}
        for side in ("left", "right"):
            _, s = side_samples(w, nb, side)
            d = [v for _, v in side_hp(w, nb, side)]
            DL, DH = nb[f"{side}_derivative_interval"]
            SL, SH = nb[f"{side}_scale_interval"]
            if s.min() < SL - 1e-12 or s.max() > SH + 1e-12:
                failures.append(f"boundary {k}: the {side} scale interval does not contain s")
            if min(d) < DL - 1e-12 or max(d) > DH + 1e-12:
                enclose = False
            signs[side] = "negative" if max(d) < 0 else ("positive" if min(d) > 0 else "mixed")
        if not enclose:
            failures.append(f"boundary {k}: a one-sided derivative interval does not contain the derivative")
        want = {("negative", "negative"): "decreasing_through_switch", ("positive", "positive"): "increasing_through_switch",
                ("negative", "positive"): "cusp_minimum"}.get((signs["left"], signs["right"]), "unresolved_derivative")
        if nb["classification"] != want:
            failures.append(f"boundary {k}: the classification is not the derivative's")
        lo_ = min(nb["left_scale_interval"][0], nb["right_scale_interval"][0])
        if nb["special_certificate"] is None and (abs(nb["safe_lower_bound"] - lo_) > 0
                                                  or abs(nb["safe_lower_minus_s270_upper"] - (lo_ - s270_hi)) > 1e-18):
            failures.append(f"boundary {k}: the safe lower bound or margin")
        if nb["special_certificate"] is None and nb["scale_exclusion_passed"] != (nb["safe_lower_minus_s270_upper"] > 0):
            failures.append(f"boundary {k}: the scale exclusion flag")
    return failures


def u7_copies(w):
    failures = []
    S = w["summary"]
    if S["special_branch_certificate"] != w["special"]:
        failures.append("the summary's special certificate is not the data file")
    D, R, NB = S["derivative_boxes"], S["stationary_roots"], S["boundary_neighborhoods"]
    BX, roots, nbs = w["boxes"], w["roots"], w["neighbourhoods"]
    if (D["total_certified_subboxes"], D["total_unresolved_subboxes"], D["interval_count"]) != (BX["total_certified_subboxes"], BX["total_unresolved_subboxes"], len(BX["summary"])):
        failures.append("the summary's box counts are not the data's")
    if R["root_count"] != len(roots) or R["smooth_minimum_count"] != sum(r["kind"] == "smooth_minimum" for r in roots) \
            or R["smooth_maximum_count"] != sum(r["kind"] == "smooth_maximum" for r in roots) \
            or R["all_second_derivative_boxes_exclude_zero"] != all(not (r["second_derivative_interval"][0] <= 0 <= r["second_derivative_interval"][1]) for r in roots):
        failures.append("the summary's root counts are not the data's")
    if NB["boundary_count"] != len(nbs) or NB["side_derivative_signs_fully_resolved"] != sum(n["derivative_side_signs_resolved"] for n in nbs) \
            or NB["scale_exclusion_passed_count"] != sum(n["scale_exclusion_passed"] for n in nbs) \
            or NB["unresolved_derivative_but_scale_safe"] != [n["boundary_index"] for n in nbs if not n["derivative_side_signs_resolved"]] \
            or NB["globally_safe_count"] != sum(n["global_exclusion_safe"] for n in nbs):
        failures.append("the summary's boundary counts are not the data's")
    if S["global_candidate"]["s270_interval"] != w["special"]["s270_interval"] or S["global_candidate"]["event_control_comparison_interval"] != w["special"]["smooth_minus_event_interval"]:
        failures.append("the summary's global candidate is not the special certificate")
    for row, s in zip(w["box_csv"], BX["summary"]):
        same = (int(row["interval_index"]) == s["interval_index"] and float(row["start"]) == s["start"] and float(row["end"]) == s["end"]
                and row["signature"] == s["signature"] and float(row["guard"]) == s["guard"] and int(row["root_count"]) == s["root_count"]
                and int(row["certified_subbox_count"]) == s["certified_subbox_count"] and row["status"] == s["status"])
        if not same:
            failures.append(f"the box CSV's row {row['interval_index']} is not the summary's")
    for row, r in zip(w["root_csv"], roots):
        if int(row["interval_index"]) != r["interval_index"] or row["kind"] != r["kind"] or float(row["root_center"]) != r["root_center"] \
                or json.loads(row["second_derivative_interval"]) != r["second_derivative_interval"] or json.loads(row["scale_interval"]) != r["scale_interval"]:
            failures.append(f"the root CSV's row for {r['root_center']:.6f} is not the data's")
    if len(w["box_csv"]) != 18 or len(w["root_csv"]) != len(roots):
        failures.append("the CSVs' row counts")
    return failures


def u8_documents(w):
    failures = []
    rp = w["report"].replace("\r\n", "\n")
    rs = gate.squash(rp)
    rd = "\n".join(l.rstrip() for l in w["readme"].replace("\r\n", "\n").split("\n"))
    SP, BX, roots, nbs = w["special"], w["boxes"], w["roots"], w["neighbourhoods"]
    # section 2: the boundary table
    rows = re.findall(r"^\|\s*(\d+)\s*\|\s*\$([^$]+)\$\s*\|\s*([\d.]+)\s*\|\s*([\d.]+)°\s*\|", rp, re.M)
    if len(rows) != 19:
        failures.append(f"section 2: {len(rows)} rows")
    for (k, fm, ph, dg), b in zip(rows, w["boundaries"]):
        if fm != b["formula"] or not r10.rounds_to(ph, b["phase"]) or not r10.rounds_to(dg, b["degrees"]):
            failures.append(f"section 2 row {k} is not the data's")
    m = re.search(r"(\d\.\d+)e-14\.", rs)
    if not m or not r10.rounds_to(m.group(1) + "e-14", w["summary"]["exact_contact_partition"]["maximum_difference_from_round10_numeric_edges"]):
        failures.append("section 2's largest difference is not the data's")
    # section 3: the special intervals
    m = re.search(r"s_\{270\}\\in\[([\d.]+),\\,([\d.]+)\]", rs)
    if not m or not (r10.rounds_to(m.group(1), SP["s270_interval"][0]) and r10.rounds_to(m.group(2), SP["s270_interval"][1])):
        failures.append("section 3's s270 interval is not the data's")
    m = re.search(r"s_\{120\}-s_\{270\}\\in\[([\d.e-]+),\\,([\d.e-]+)\]", rs)
    if not m or not (r10.rounds_to(m.group(1), SP["s120_minus_s270_interval"][0]) and r10.rounds_to(m.group(2), SP["s120_minus_s270_interval"][1])):
        failures.append("section 3's difference interval is not the data's")
    # section 4: the gain interval
    m = re.search(r"\[(0\.0000105819\d+),\\,(0\.0000105819\d+)\]", rs)
    if not m or [m.group(1), m.group(2)] != SP["smooth_minus_event_interval"]:
        failures.append("section 4's interval is not the data's")
    # section 5: the box table
    S = {s["interval_index"]: s for s in BX["summary"]}
    rows = re.findall(r"^\|\s*(\d+)\s*\|\s*`([^`]+)`\s*\|\s*(\d+)\s*\|\s*(\d+)\s*\|\s*([a-z, ]+)\s*\|\s*(\d+)\s*\|", rp, re.M)
    if len(rows) != 18:
        failures.append(f"section 5: {len(rows)} rows")
    for k, sg, rc, n, signs, un in rows:
        s = S[int(k)]
        if sg != s["signature"] or int(rc) != s["root_count"] or int(n) != s["certified_subbox_count"] \
                or [x.strip() for x in signs.split(",")] != s["certified_signs"] or int(un) != s["unresolved_subbox_count"]:
            failures.append(f"section 5 row {k} is not the data's")
    # section 6: the root table
    rows = re.findall(r"^\|\s*(\d+)\s*\|\s*(\d+)\s*\|\s*(smooth_\w+)\s*\|\s*([\d.]+)\s*\|\s*\[([-\d.e+]+),\s*([-\d.e+]+)\]\s*\|\s*([\d.e+-]+)\s*\|", rp, re.M)
    if len(rows) != len(roots):
        failures.append(f"section 6: {len(rows)} rows")
    for (j, i, kd, ctr, d2l, d2h, mar), r in zip(rows, roots):
        if int(i) != r["interval_index"] or kd != r["kind"] or not r10.rounds_to(ctr, r["root_center"]) \
                or not r10.rounds_to(d2l, r["second_derivative_interval"][0]) or not r10.rounds_to(d2h, r["second_derivative_interval"][1]) \
                or not r10.rounds_to(mar, r["scale_lower_minus_s270_upper"]):
            failures.append(f"section 6 row {j} is not the data's")
    # section 7: the boundary table
    rows = re.findall(r"^\|\s*(\d+)\s*\|\s*\$([^$]+)\$\s*\|\s*(\w+)\s*\|\s*(True|False)\s*\|\s*(-?[\d.e+-]+)\s*\|", rp, re.M)
    if len(rows) != 17:
        failures.append(f"section 7: {len(rows)} rows")
    for (k, fm, cl, res, mar), n in zip(rows, nbs):
        if fm != n["formula"] or cl != n["classification"] or (res == "True") != n["derivative_side_signs_resolved"] \
                or not r10.rounds_to(mar, n["safe_lower_minus_s270_upper"]):
            failures.append(f"section 7 row {k} is not the data's")
    # the README and the counts in the text
    for label, v in (("精確接觸區間", 18), ("自適應導數盒", len(BX["boxes"])), ("未解析導數盒", BX["total_unresolved_subboxes"]),
                     ("光滑駐點根盒", len(roots))):
        m = re.search(label + r"：(\d+)", rd)
        if not m or int(m.group(1)) != v:
            failures.append(f"the README's {label} is not the data's")
    m = re.search(r"\$s_\{120\}-s_\{270\}\$ 下界：([\d.e-]+)", rd)
    if not m or not r10.rounds_to(m.group(1), SP["s120_minus_s270_interval"][0]):
        failures.append("the README's difference bound is not the data's")
    m = re.search(r"平滑候選減事件控制下界：([\d.]+)", rd)
    if not m or m.group(1) != SP["smooth_minus_event_interval"][0]:
        failures.append("the README's gain bound is not the data's")
    if not re.search(r"共有：\$\$12\$\$個光滑駐點", rs) or "12\\text{個唯一駐點根盒}" not in rs:
        failures.append("the report's count of stationary points is not 12")
    # the audit note, the handoff, the manifest
    an = w["audit_note"]
    if f"[{SP['s270_interval'][0]!r}, {SP['s270_interval'][1]!r}]" not in an or "自適應子盒：579" not in an or "光滑駐點：12" not in an:
        failures.append("the audit note's figures are not the data's")
    if "N(X_r)\\subset\\operatorname{int}X_r" not in gate.squash(w["handoff"]):
        failures.append("the handoff does not state its root criterion")
    M = w["manifest"]
    if M["derivative_subbox_count"] != len(BX["boxes"]) or M["stationary_root_count"] != len(roots) \
            or M["s120_minus_s270_lower"] != SP["s120_minus_s270_interval"][0] or M["smooth_minus_event_lower"] != SP["smooth_minus_event_interval"][0] \
            or sorted(M["files"]) != sorted(set(w["members"]) - {"MANIFEST.json"}):
        failures.append("the manifest is not the package's")
    return failures


def u9_method(w):
    failures = []
    V = w["verifier"]
    D = w["round09"]["definitions"]["candidate"]
    for name, k in (("EPS", "epsilon"), ("W", "wing_length"), ("BETA", "beta"), ("DELTA", "delta"), ("CENTER", "transition_center")):
        m = re.search(rf"^{name}=([\d.]+)$", V, re.M)
        if not m or m.group(1) != D[k]:
            failures.append(f"the verifier's {name} is not Round 09's")
    m = midpoint_sums(w, 262144)
    K = consts(w)
    with mp.workdps(40):
        den = mp.tanh(K["transition_center"] / K["epsilon"]) - mp.tanh((K["transition_center"] - 1) / K["epsilon"])
        M1 = K["delta"] / (K["epsilon"] * den)
        M2 = M1 ** 2 + K["delta"] * (4 / (3 * mp.sqrt(3))) / (K["epsilon"] ** 2 * den)
    if abs(m["M1"] - float(M1)) > 1e-12 or abs(m["M2"] - float(M2)) > 1e-10:
        failures.append("the verifier's bounds are not |theta'| and |theta'|^2 + |theta''|")
    for frag in ("M_COORD_2=M_THETA_1**2+M_THETA_2", "E_xy=M_COORD_2/(24*N*N)+2e-15", "q=qa*np.sin(th)-qb*np.cos(th)"):
        if frag not in V:
            failures.append(f"the verifier lacks {frag}")
    if "|f''(u)|" not in w["method"] or "\\|\\theta'\\|_\\infty^2" not in w["method"]:
        failures.append("the method note does not state the bound")
    return failures


CHECKS = {"U1": u1_boundaries, "U2": u2_special, "U3": u3_event, "U4": u4_boxes, "U5": u5_roots, "U6": u6_neighbourhoods,
          "U7": u7_copies, "U8": u8_documents, "U9": u9_method}
CURVE = {"round10", "round09", "round05", "legs"}
DEPENDS = {
    "U1": CURVE | {"boundaries", "summary"},
    "U2": CURVE | {"special"},
    "U3": {"special", "round09"},
    "U4": CURVE | {"boxes", "roots"},
    "U5": CURVE | {"roots", "special"},
    "U6": CURVE | {"neighbourhoods", "special"},
    "U7": {"summary", "special", "boxes", "roots", "neighbourhoods", "box_csv", "root_csv"},
    "U8": {"readme", "report", "audit_note", "handoff", "manifest", "members", "boundaries", "special", "boxes", "roots", "neighbourhoods", "summary"},
    "U9": {"verifier", "method", "round09", "round10"},
}


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

def findings(w):
    out = {}
    M = model(w)
    nbrows = []
    for nb in w["neighbourhoods"]:
        row = {"boundary": nb["boundary_index"], "formula": nb["formula"], "classification": nb["classification"],
               "scale_exclusion_passed": nb["scale_exclusion_passed"]}
        row["radius"] = nb["neighborhood_radius"]
        for side in ("left", "right"):
            hp_ = side_hp(w, nb, side)
            d = [float(v) for _, v in hp_]
            DL, DH = nb[f"{side}_derivative_interval"]
            at0, at2 = d[0], d[1]
            # the interval's end nearest the switch, against the derivative 2e-10 rad from the switch
            near = DL if (DL - at2) ** 2 < (DH - at2) ** 2 else DH
            row[side] = {"at_switch": at0, "at_2e-10": at2, "range_60_digits": [min(d), max(d)], "stated": [DL, DH],
                         "encloses": bool(min(d) >= DL - 1e-12 and max(d) <= DH + 1e-12),
                         "shortfall": max(DL - min(d), max(d) - DH, 0.0),
                         "near_end_minus_derivative_at_2e-10": near - at2, "move_within_2e-10": abs(at0 - at2)}
        nbrows.append(row)
    out["neighbourhoods"] = nbrows
    out["neighbourhoods_not_enclosing"] = [r["boundary"] for r in nbrows if not (r["left"]["encloses"] and r["right"]["encloses"])]
    out["largest_shortfall"] = max(max(r["left"]["shortfall"], r["right"]["shortfall"]) for r in nbrows)
    # the boundary formulas' residue under their switches' conditions, at 60 digits
    C = r11.cand(w11(w), 60)
    with mp.workdps(60):
        out["boundary_condition_residue"] = mp.nstr(max(switch_condition(C, s, boundary_phase(w, i + 1))
                                                       for i, s in enumerate(w["round10"]["contact_switches"])), 3)
    out["smooth_minimum_count"] = sum(1 for e in stationary(w) if e["kind"] == "smooth_minimum")
    out["root_second_derivative_margin"] = min(min(min(v) - r["second_derivative_interval"][0], r["second_derivative_interval"][1] - max(v))
                                               for r, v in ((r, root_second(w, j)) for j, r in enumerate(w["roots"])))
    stat = stationary(w)
    out["stationary_count"] = len(stat)
    nb10 = next(n for n in w["neighbourhoods"] if n["boundary_index"] == 10)
    extra = [e for e in stat if not any(r["uniqueness_box"][0] <= e["phase"] <= r["uniqueness_box"][1] for r in w["roots"])]
    with mp.workdps(40):
        out["stationary_outside_root_boxes"] = [{"phase": mp.nstr(e["phase"], 20), "interval": e["interval"],
                                                 "before_boundary_10": mp.nstr(mp.mpf(nb10["phase"]) - e["phase"], 5),
                                                 "within_boundary_10_neighbourhood": bool(abs(mp.mpf(nb10["phase"]) - e["phase"]) <= nb10["neighborhood_radius"])}
                                                for e in extra]
    # the largest gaps between a box's derivative interval and the derivative sampled in it
    slack = [min(box_samples(w, b).min() - b["derivative_interval"][0], b["derivative_interval"][1] - box_samples(w, b).max())
             for b in w["boxes"]["boxes"]]
    out["boxes_least_slack"] = float(min(slack))
    # the special enclosures against this run's (RUN-011's certificate)
    H = r11.hp(w11(w))
    with mp.workdps(60):
        s = r10.s270(H)
        out["s270_minus_interval_ends"] = [mp.nstr(s - mp.mpf(repr(x)), 5) for x in w["special"]["s270_interval"]]
        s120 = r11.cand(w11(w), 60).s(2 * mp.pi / 3)
        out["s120_minus_s270_exact"] = mp.nstr(s120 - s, 15)
    m = midpoint_sums(w, w["special"]["quadrature_partition_count"])
    out["recomputed_bounds"] = {"M1": m["M1"], "M2": m["M2"], "coordinate": m["M2"] / (24 * 262144 ** 2) + 2e-15}
    return out


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

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

    def rrep(old, new):
        assert rp.count(old) == 1, old
        return rp.replace(old, new)
    # a boundary's formula, with section 2's cell alike
    bf = edit("boundaries", lambda B: B[3].__setitem__("formula", "\\beta-\\pi/4"))
    rep_bf = rrep("| 3 | $\\beta-\\pi/3$ | 0.358185243642795 |", "| 3 | $\\beta-\\pi/4$ | 0.358185243642795 |")
    # a boundary's phase moved by 1e-14, with section 2's cell alike
    bp = edit("boundaries", lambda B: B[5].__setitem__("phase", B[5]["phase"] + 1e-14))
    rep_bp = rrep("| 5 | $\\pi/2$ | 1.570796326794897 |", f"| 5 | $\\pi/2$ | {bp[5]['phase']:.15f} |")
    # the special certificate and the summary's copy alike: an error bound, the s120 interval, the event branch
    def sp(fn):
        return edit("special", fn), edit("summary", lambda s: fn(s["special_branch_certificate"]))
    eb = sp(lambda c: c.__setitem__("coordinate_midpoint_error_bound", c["coordinate_midpoint_error_bound"] * 1.01))
    s120 = sp(lambda c: c["s120_interval"].__setitem__(0, c["s120_interval"][0] + 3e-11))
    b1 = sp(lambda c: c.__setitem__("event_control_branch_B1", c["event_control_branch_B1"][:30] + str((int(c["event_control_branch_B1"][30]) + 1) % 10) + c["event_control_branch_B1"][31:]))
    # a box's interval shrunk below the derivative it holds (its sign kept)
    bx = edit("boxes", lambda B: B["boxes"][100]["derivative_interval"].__setitem__(1, 8.0e-05))
    # an interval's box count, in the summary, the CSV and section 5 alike
    bc = edit("boxes", lambda B: B["summary"][4].__setitem__("certified_subbox_count", 44))
    bcsv = edit("box_csv", lambda R: R[4].__setitem__("certified_subbox_count", "44"))
    rep_bc = rrep("| 4 | `p1|p0|p3` | 1 | 43 |", "| 4 | `p1|p0|p3` | 1 | 44 |")
    # a root box's second-derivative interval made to hold 0, with the CSV, section 6 and the summary's flag alike
    rb = edit("roots", lambda R: R[3]["second_derivative_interval"].__setitem__(1, 0.5))
    rcsv = edit("root_csv", lambda R: R[3].__setitem__("second_derivative_interval", json.dumps(rb[3]["second_derivative_interval"])))
    rep_rb = rrep("| 4 | 4 | smooth_maximum | 0.994829121177570 | [-1.247879e+00, -1.247879e+00] |",
                  "| 4 | 4 | smooth_maximum | 0.994829121177570 | [-1.247879e+00, 5.000000e-01] |")
    rsum = edit("summary", lambda s: s["stationary_roots"].__setitem__("all_second_derivative_boxes_exclude_zero", False))
    # a neighbourhood's right scale interval raised above s there (its safe lower bound, the left one, unchanged)
    nbs = edit("neighbourhoods", lambda N: N[3]["right_scale_interval"].__setitem__(0, 1.0406501169233402 + 1e-9))
    # a neighbourhood's classification, with section 7's cell alike
    nbc = edit("neighbourhoods", lambda N: N[4].__setitem__("classification", "increasing_through_switch"))
    rep_nbc = rrep("| 5 | $\\pi/2$ | decreasing_through_switch |", "| 5 | $\\pi/2$ | increasing_through_switch |")
    cnt = edit("summary", lambda s: s["boundary_neighborhoods"].__setitem__("scale_exclusion_passed_count", 15))
    rep_sign = rrep("| 3 | `p1|p0|R` | 0 | 1 | positive | 0 |", "| 3 | `p1|p0|R` | 0 | 1 | negative | 0 |")
    readme = base["readme"].replace("- 自適應導數盒：579", "- 自適應導數盒：578")
    handoff = base["handoff"].replace("N(X_r)\\subset\\operatorname{int}X_r.", "N(X_r)\\subseteq X_r.")
    ver = base["verifier"].replace("EPS=0.037", "EPS=0.038")
    meth = base["method"].replace("\\|\\theta'\\|_\\infty^2", "\\|\\theta'\\|_\\infty")
    for x, y in ((readme, base["readme"]), (handoff, base["handoff"]), (ver, base["verifier"]), (meth, base["method"])):
        assert x != y, "a plant must hit the raw text"
    return [
        ("K0-unmodified", "control", {}, set()),
        ("K1-report-crlf", "control", {"report": base["report"].replace("\n", "\r\n")}, set()),
        ("K2-readme-trailing-spaces", "control", {"readme": "\n".join(l + " " for l in base["readme"].split("\n"))}, set()),
        ("K3-boxes-copied", "control", {"boxes": copy.deepcopy(base["boxes"])}, set()),
        ("M1-triangle-leg", "mutation", {"legs": ("1.0005", "1")}, {"U1", "U2", "U4", "U5", "U6"}),
        ("M2-boundary-formula", "mutation", {"boundaries": bf, "report": rep_bf}, {"U1"}),
        ("M3-boundary-phase", "mutation", {"boundaries": bp, "report": rep_bp}, {"U1"}),
        ("M4-error-bound", "mutation", {"special": eb[0], "summary": eb[1]}, {"U2"}),
        ("M5-s120-interval", "mutation", {"special": s120[0], "summary": s120[1]}, {"U2"}),
        ("M6-event-branch", "mutation", {"special": b1[0], "summary": b1[1]}, {"U3"}),
        ("M7-box-interval", "mutation", {"boxes": bx}, {"U4"}),
        ("M8-box-count", "mutation", {"boxes": bc, "box_csv": bcsv, "report": rep_bc}, {"U4"}),
        ("M9-root-second-derivative", "mutation", {"roots": rb, "root_csv": rcsv, "report": rep_rb, "summary": rsum}, {"U5"}),
        ("M10-neighbourhood-scale", "mutation", {"neighbourhoods": nbs}, {"U6"}),
        ("M11-neighbourhood-class", "mutation", {"neighbourhoods": nbc, "report": rep_nbc}, {"U6"}),
        ("M12-summary-count", "mutation", {"summary": cnt}, {"U7"}),
        ("M13-report-sign", "mutation", {"report": rep_sign}, {"U8"}),
        ("M14-readme-count", "mutation", {"readme": readme}, {"U8"}),
        ("M15-handoff-criterion", "mutation", {"handoff": handoff}, {"U8"}),
        ("M16-verifier-constant", "mutation", {"verifier": ver}, {"U9"}),
        ("M17-method-bound", "mutation", {"method": meth}, {"U9"}),
    ]


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[:10]:
            print("    ", x)
    log = {"run": "RUN-012",
           "subject": "Round 11 (v0.11): exact boundaries, derivative boxes, root boxes, boundary neighbourhoods, special enclosures",
           "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, ensure_ascii=False, default=str)[:9000])
    if args.drill:
        drill = gate.run_drill(build_world, CHECKS, DEPENDS, drill_cases(build_world()))
        log["drill"] = gate.summarize(drill, CHECKS)
        print(gate.report_drill(log["drill"]))
        for d in drill:
            if not d["ok"]:
                print("   DRILL FAIL", d)
    if args.blind:
        seen = sum(1 for d in log["drill"]["cases"] if args.blind in d["expected"] and args.blind not in d["tripped"])
        print(f"blind {args.blind}: reported blind in {seen} cases (log not written)")
        return
    if args.drill:
        path = pathlib.Path(args.log) if args.log else LOG
        path.parent.mkdir(parents=True, exist_ok=True)
        path.write_text(json.dumps(log, ensure_ascii=False, indent=1, default=str) + "\n", encoding="utf-8", newline="\n")


if __name__ == "__main__":
    main()