"""RUN-015 - Round 14 (v0.14): double-peak curvature, chirality breaking and the single peak's local stability.

Round 14 widens Round 13's curve family in two ways: the curvature density of a wing becomes a mixture of two
normalized sech^2 layers, (1 - a) rho(c1, e1) + a rho(c2, e2), with c1,2 = c -/+ d/2, e1,2 = e exp(-/+ r) and
a = sigmoid(q) - eight parameters; and the mirror symmetry is broken by two modes, the local centres moved alike
(c + eta_c on the left, 1 - c + eta_c on the right) and the widths split (e exp(+/- eta_e)), scored by the weaker
handedness. Its search collapses back to one peak (separation 3.3e-6), 6.55e-10 below the single peak; a separation
profile at Round 13's curve falls as -0.0338 d^2; the congruent objective peaks at zero offset in both modes and in an
80-offset census; and it concludes that the single peak is locally stable against splitting and the mirror symmetry
against both modes.

Checks (each returns a list of failures; empty = pass)
  Y1 curves       the reference is Round 13's common scale; the two curve files are RUN-009's polyline (12,000 wing
                  steps, 160 centre steps) of Round 13's curve and of the search candidate's mixture, centred at their
                  mean, to 1e-14; the candidate's centres, widths and weight are its parameters' arithmetic
  Y2 search       at 40 digits, with an exact support for the mixture: the search candidate's high-resolution cusp
                  values are the exact ones (to 1e-13), its low values below the exact ones by no more than 1e-9 (a
                  polyline's deficit); its least is the exact least's; the differences are their figures' arithmetic;
                  the double-minus-single difference is the curves' (fixed-text failure)
  Y3 profile      the separation profile is Round 13's curve split into two equal halves: at seven separations its cusp
                  values are the exact ones (to 1e-12) and its low values below the exact ones by no more than 1e-8 (at
                  their stated phases, where those are interior); the losses are the scales less the first; the fit is
                  the profile's own (a quadratic in d^2 over d <= 0.01); the low phases are stationary points, not the
                  ends of a search interval (fixed-text failure)
  Y4 splitting    the profile's coefficient is half the low branches' second derivative along d with the five
                  parameters fixed (to 1 %); at RUN-014's solution, the Lagrangian of the four branches in the
                  double-peak directions (d, r) - first derivatives zero by the exchange of the two halves - is negative
                  semidefinite, as "locally stable against splitting" requires (fixed-text failure)
  Y5 chirality    the profiles' congruent scale is the weaker handedness and the reflection exchanges them; the census'
                  best is its largest; at RUN-014's solution the two modes' gradients of the four branches (the cusps
                  where their sides tie, which moves) have rank 2, so the symmetric point is a strict local maximum of
                  the congruent scale in both modes; each fit is its window's fit of the shipped profile (fixed-text
                  failures); the congruent scale has a corner at zero (fixed-text failure)
  Y6 phase        the listed local minima of the base curve against the exact minima (RUN-014's phase map), to 1e-9
                  in scale (fixed-text failure); the double-peak list is evaluated at the base list's phases; the
                  sampled global minimum is the 270-degree cusp
  Y7 copies       the summary holds the data files' records; the CSVs are the JSON profiles; the manifest
  Y8 documents    the README, the report's tables and figures, the search audit and the handoff state the data's

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/run015_round14_double_peak.py [--drill [--log PATH] | --blind Y# | --dry]
"""
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
import run009_round08_smooth_layer as r9
import run010_round09_precision_cusp as r10
import run011_round10_phase_map as r11
import run014_round13_smooth_kkt as r14
from run010_round09_precision_cusp import cusp_270, dec, hi, iv, lo, tanh_iv

RID = "round-v0.14"
LOG = TREE / "data" / "gate-logs" / "run015-round14-double-peak.json"
DPS = 40
NAMES = ("min_x", "min_y", "max_d")
SIG = r14.SIG
P13KEYS = ("wing_length", "beta_radians", "delta_radians", "transition_center", "epsilon")


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

def build_world():
    f = load(RID)
    md = lambda n: restore_escapes(f[n].decode("utf-8"))
    js = lambda n: json.loads(f[n])
    rows = lambda n: list(csv.DictReader(io.StringIO(f[n].decode("utf-8"))))
    arr = lambda n: np.loadtxt(io.StringIO(f[n].decode("utf-8")), delimiter=",", skiprows=1)
    return {
        "readme": md("README.md"), "report": md("reports/ROUND_14_REPORT.md"), "audit_note": md("reports/ROUND_14_SEARCH_AUDIT.md"),
        "handoff": md("reports/ROUND_15_CURVATURE_MODE_HANDOFF.md"), "source": f["src/curvature_distribution_models.py"].decode("utf-8"),
        "manifest": js("MANIFEST.json"), "members": sorted(f),
        "summary": js("data/round14_summary.json"), "search": js("data/double_peak_search.json"),
        "hr": js("data/double_peak_high_resolution_audit.json"), "asym": js("data/asymmetry_congruence_audit.json"),
        "sep_csv": rows("data/double_peak_separation_profile.csv"), "centre_csv": rows("data/asymmetry_center_profile.csv"),
        "width_csv": rows("data/asymmetry_width_profile.csv"),
        "base_curve": arr("data/single_peak_base_curve.csv"), "double_curve": arr("data/double_peak_best_curve.csv"),
        "round13": r14.build_world(),
        "legs": ("1", "1"),
    }


def w13(w):
    """Round 13's world, with this world's triangle"""
    return dict(w["round13"], legs=w["legs"])


def p13(w):
    P = w["round13"]["summary"]["round13_candidate"]["parameters"]
    return [P[k] for k in P13KEYS]


# ---------------------------------------------------------------- an exact support for mixtures of tanh layers

class Layers:
    """a wing's F(u) in [0, 1] as a mixture of tanh layers: falling from the wing's end (left wing, u from p0) or rising
    from the junction (right wing, u from p2); its inverse by the inverse tanh for one layer, else by bisection and
    Newton steps"""

    def __init__(self, comps, rising):
        self.rising = rising
        self.comps = []
        for a, c, e in comps:
            if rising:
                T0, T1 = mp.tanh(-c / e), mp.tanh((1 - c) / e)
            else:
                T0, T1 = mp.tanh(c / e), mp.tanh((c - 1) / e)
            self.comps.append((a, c, e, T0, T1))
        self.breaks = sorted({min(max(c, mp.mpf(0)), mp.mpf(1)) for _, c, _, _, _ in self.comps} - {mp.mpf(0), mp.mpf(1)})

    def F(self, u):
        if self.rising:
            return mp.fsum(a * (mp.tanh((u - c) / e) - T0) / (T1 - T0) for a, c, e, T0, T1 in self.comps)
        return mp.fsum(a * (mp.tanh((c - u) / e) - T1) / (T0 - T1) for a, c, e, T0, T1 in self.comps)

    def dF(self, u):
        sg = 1 if self.rising else -1
        return sg * mp.fsum(a / (e * abs(T1 - T0)) * mp.sech((u - c) / e) ** 2 for a, c, e, T0, T1 in self.comps)

    def u_of(self, f):
        if len(self.comps) == 1:
            _, c, e, T0, T1 = self.comps[0]
            return c + e * mp.atanh(T0 + f * (T1 - T0)) if self.rising else c - e * mp.atanh(T1 + f * (T0 - T1))
        lo_, hi_ = mp.mpf(0), mp.mpf(1)
        for _ in range(60):
            m = (lo_ + hi_) / 2
            if (self.F(m) < f) == self.rising:
                lo_ = m
            else:
                hi_ = m
        u = (lo_ + hi_) / 2
        for _ in range(6):
            u -= (self.F(u) - f) / self.dF(u)
        return u


class Curve:
    """the five-piece curve: left wing (angle beta + delta F, F falling), the centre l0 = 1 - 2w, right wing (angle
    -(beta + delta F_R), F_R rising from the junction); positions by mpmath's quadrature split at the layers"""

    def __init__(self, w, b, dl, left, right, legs):
        self.w, self.b, self.dl, self.Lw, self.Rw = w, b, dl, Layers(left, False), Layers(right, True)
        r3 = mp.sqrt(3)
        A = (3 + 4 * r3) / 18 / mp.mpf("1.0048") * mp.mpf(legs[0])
        B = (4 + r3) / 6 / mp.mpf("1.0048") * mp.mpf(legs[1])
        self.fun = {"max_d": (1 / A, 1 / B), "min_x": (-1 / A, mp.mpf(0)), "min_y": (mp.mpf(0), -1 / B)}
        self.l0 = 1 - 2 * w
        x1, y1 = self.pos(self.Lw, mp.mpf(1), 1)
        rx, ry = self.pos(self.Rw, mp.mpf(1), -1)
        self.P = {"p0": (mp.mpf(0), mp.mpf(0)), "p1": (x1, y1), "p2": (x1 + self.l0, y1), "p3": (x1 + self.l0 + rx, y1 + ry)}

    def pos(self, lay, u, sg):
        th = lambda v: sg * (self.b + self.dl * lay.F(v))
        pts = [mp.mpf(0)] + [x for x in lay.breaks if x < u] + [u]
        return (self.w * mp.quad(lambda v: mp.cos(th(v)), pts), self.w * mp.quad(lambda v: mp.sin(th(v)), pts))

    def point(self, label, name, ph):
        if label in self.P:
            return self.P[label]
        a = self.fun[name]
        cs, sn = mp.cos(ph), mp.sin(ph)
        n = (a[0] * cs + a[1] * sn, -a[0] * sn + a[1] * cs)
        ang = mp.atan2(n[1], n[0])
        sg = 1 if label == "L" else -1
        best = None
        for k in (-1, 1):
            t = ang + k * mp.pi / 2
            t -= 2 * mp.pi * mp.floor((t + mp.pi) / (2 * mp.pi))
            f = (sg * t - self.b) / self.dl
            if best is None or abs(f - mp.mpf(0.5)) < abs(best - mp.mpf(0.5)):
                best = f
        f = min(max(best, mp.mpf(0)), mp.mpf(1))
        lay = self.Lw if label == "L" else self.Rw
        if 0 < f < 1:
            u = lay.u_of(f)
        else:
            u = mp.mpf(0) if (f == 1) == (label == "L") else mp.mpf(1)
        if label == "L":
            return self.pos(self.Lw, u, 1) if u > 0 else (mp.mpf(0), mp.mpf(0))
        px, py = self.pos(self.Rw, u, -1) if u > 0 else (mp.mpf(0), mp.mpf(0))
        return (self.P["p2"][0] + px, self.P["p2"][1] + py)

    def value(self, sig, ph):
        cs, sn = mp.cos(ph), mp.sin(ph)
        tot = mp.mpf(0)
        for name, lab in zip(NAMES, sig.split("|")):
            a = self.fun[name]
            p = self.point(lab, name, ph)
            tot += (a[0] * cs + a[1] * sn) * p[0] + (-a[0] * sn + a[1] * cs) * p[1]
        return tot

    def deriv(self, sig, ph):
        cs, sn = mp.cos(ph), mp.sin(ph)
        tot = mp.mpf(0)
        for name, lab in zip(NAMES, sig.split("|")):
            a = self.fun[name]
            p = self.point(lab, name, ph)
            tot += (-a[0] * sn + a[1] * cs) * p[0] + (-a[0] * cs - a[1] * sn) * p[1]
        return tot


_memo = {}


def memo(*key):
    def deco(fn):
        k = json.dumps(key, default=str)
        if k not in _memo:
            _memo[k] = fn()
        return _memo[k]
    return deco


def mirror(comps):
    return [(a, 1 - c, e) for a, c, e in comps]


def double(p, d, r=0, q=0):
    """the mixture's wing layers for parameters (w, beta, delta, c, eps) and (d, r, q): c -/+ d/2, eps exp(-/+ r),
    weight sigmoid(q) on the second"""
    w, b, dl, c, e = p
    a2 = 1 / (1 + mp.exp(-q))
    return [(1 - a2, c - d / 2, e * mp.exp(-r)), (a2, c + d / 2, e * mp.exp(r))]


def sym_curve(w, p, d=0, r=0, q=0):
    with mp.workdps(DPS):
        p = [mp.mpf(x) for x in p]
        L = double(p, mp.mpf(d), mp.mpf(r), mp.mpf(q))
        return Curve(p[0], p[1], p[2], L, mirror(L), w["legs"])


def low_min(C, ph0, span="1e-3", sig="L|p0|p2"):
    """the smooth minimum of s with a signature within span of ph0, by bisection on the derivative's sign; None if the
    derivative does not change sign there"""
    with mp.workdps(DPS):
        a, b = mp.mpf(ph0) - mp.mpf(span), mp.mpf(ph0) + mp.mpf(span)
        if not C.deriv(sig, a) < 0 < C.deriv(sig, b):
            return None
        for _ in range(100):
            m = (a + b) / 2
            if C.deriv(sig, m) < 0:
                a = m
            else:
                b = m
        q = (a + b) / 2
        return q, C.value(sig, q)


def four(C, ph):
    """the four branches: the two low minima near ph and the two cusps"""
    with mp.workdps(DPS):
        lows = [low_min(C, x) for x in ph]
        return {"phases": [None if m is None else m[0] for m in lows],
                "values": [None if m is None else m[1] for m in lows] + [C.value(SIG[2], 2 * mp.pi / 3), C.value(SIG[3], 3 * mp.pi / 2)]}


# ---------------------------------------------------------------- the computations the checks read

def search_exact(w):
    B = w["summary"]["double_peak"]["search_candidate"]
    def go():
        with mp.workdps(DPS):
            p = [mp.mpf(repr(x)) for x in B["parameters"]]
            return four(sym_curve(w, p[:5], p[5], p[6], p[7]), B["phases"])
    return memo("search", B["parameters"], w["legs"])(go)


PROFILE_D = (0.0, 0.002, 0.004, 0.006, 0.008, 0.01, 0.024)


def profile_exact(w, d, phis):
    def go():
        with mp.workdps(DPS):
            return four(sym_curve(w, [mp.mpf(repr(x)) for x in p13(w)], mp.mpf(repr(d))), phis)
    return memo("profile", d, phis, p13(w), w["legs"])(go)


def stationary_phases(w, d):
    """the low branch's stationary points of the profile curve at separation d, in [0.10, 0.18]: sign changes of the
    derivative on 161 phases, refined by bisection"""
    def go():
        with mp.workdps(DPS):
            C = sym_curve(w, [mp.mpf(repr(x)) for x in p13(w)], mp.mpf(repr(d)))
            grid = [mp.mpf("0.10") + mp.mpf("0.0005") * k for k in range(161)]
            D = [C.deriv("L|p0|p2", x) for x in grid]
            out = []
            for (x0, d0), (x1_, d1) in zip(zip(grid, D), zip(grid[1:], D[1:])):
                if (d0 > 0) != (d1 > 0):
                    a, b = x0, x1_
                    for _ in range(80):
                        m = (a + b) / 2
                        if (C.deriv("L|p0|p2", m) > 0) == (d0 > 0):
                            a = m
                        else:
                            b = m
                    out.append(((a + b) / 2, "minimum" if d0 < 0 else "maximum"))
            return out
    return memo("stationary", d, p13(w), w["legs"])(go)


def kkt13(w):
    """RUN-014's solution of Round 13's twelve equations (60 digits)"""
    return r14.kkt(w13(w))


def splitting_block(w, h="1e-6"):
    """the Lagrangian sum mu_r m_r's second derivatives in (d, r) at the solution, the branches at fixed phases;
    and the first derivatives, which the exchange of the two halves makes zero"""
    def go():
        x = kkt13(w)["x"]
        with mp.workdps(DPS + 10):
            p, ph, mu = x[:5], [x[5], x[6], 2 * mp.pi / 3, 3 * mp.pi / 2], x[8:12]
            hh = mp.mpf(h)
            V = {}
            for i in (-1, 0, 1):
                for j in (-1, 0, 1):
                    C = sym_curve(w, p, i * hh, j * hh)
                    V[(i, j)] = [C.value(SIG[k], ph[k]) for k in range(4)]
            m = []
            for k in range(4):
                dd = (V[(1, 0)][k] - 2 * V[(0, 0)][k] + V[(-1, 0)][k]) / hh ** 2
                rr = (V[(0, 1)][k] - 2 * V[(0, 0)][k] + V[(0, -1)][k]) / hh ** 2
                dr = (V[(1, 1)][k] - V[(1, -1)][k] - V[(-1, 1)][k] + V[(-1, -1)][k]) / (4 * hh * hh)
                g = [(V[(1, 0)][k] - V[(-1, 0)][k]) / (2 * hh), (V[(0, 1)][k] - V[(0, -1)][k]) / (2 * hh)]
                m.append({"dd": dd, "dr": dr, "rr": rr, "first": g})
            L = mp.matrix([[mp.fsum(mu[k] * m[k]["dd"] for k in range(4)), mp.fsum(mu[k] * m[k]["dr"] for k in range(4))],
                           [mp.fsum(mu[k] * m[k]["dr"] for k in range(4)), mp.fsum(mu[k] * m[k]["rr"] for k in range(4))]])
            ev = sorted(mp.eigsy(L)[0])
            base = [mp.fsum([V[(0, 0)][k]]) - x[7] for k in range(4)]
        return {"branches": m, "block": L, "eigenvalues": ev, "base_minus_solution": base}
    return memo("split", w["legs"])(go)


# each cusp's two sides: where the changing support passes from one corner to the other (max_d from p3 to p0 at 2pi/3,
# min_x from p3 to p0 at 3pi/2); on a mirror-symmetric curve p0 and p3 are level and the cusp sits at 2pi/3 and 3pi/2,
# on a broken one it moves
CUSPS = ((2, "p2|p1|p3", "p2|p1|p0"), (3, "p3|p3|p1", "p0|p3|p1"))


def cusp(C, left, right, ph0, span="1e-4"):
    """the cusp near ph0: the phase where its two sides' values tie (bisection on their difference) and the value there"""
    with mp.workdps(DPS + 10):
        a, b = mp.mpf(ph0) - mp.mpf(span), mp.mpf(ph0) + mp.mpf(span)
        fa = C.value(left, a) - C.value(right, a)
        assert (fa > 0) != (C.value(left, b) - C.value(right, b) > 0)
        for _ in range(120):
            m = (a + b) / 2
            if (C.value(left, m) - C.value(right, m) > 0) == (fa > 0):
                a = m
            else:
                b = m
        q = (a + b) / 2
        return q, C.value(left, q)


def chirality(w, h="1e-10"):
    """the two modes' derivatives of the four branches at the solution (the right wing's own layer): the low branches at
    their fixed phases (stationary, so the phase's own move does not count at first order), the cusps where their sides
    tie, which moves"""
    def go():
        x = kkt13(w)["x"]
        with mp.workdps(DPS + 10):
            p, ph = x[:5], [x[5], x[6], 2 * mp.pi / 3, 3 * mp.pi / 2]
            hh = mp.mpf(h)

            def curve(ec, ee):
                return Curve(p[0], p[1], p[2], [(mp.mpf(1), p[3] + ec, p[4] * mp.exp(ee))], [(mp.mpf(1), 1 - p[3] + ec, p[4] * mp.exp(-ee))], w["legs"])

            def four_values(C):
                vals = [C.value(SIG[k], ph[k]) for k in (0, 1)]
                phases = []
                for k, left, right in CUSPS:
                    q, v = cusp(C, left, right, ph[k])
                    vals.append(v)
                    phases.append(q)
                return vals, phases
            out = {}
            for mode, (u, v) in (("centre", (hh, 0)), ("width", (0, hh))):
                (vu, qu), (vd, qd) = four_values(curve(u, v)), four_values(curve(-u, -v))
                out[mode] = [(a - b) / (2 * hh) for a, b in zip(vu, vd)]
                out[mode + "_cusp_phase_rates"] = [(a - b) / (2 * hh) for a, b in zip(qu, qd)]
            G = mp.matrix([[out["centre"][k], out["width"][k]] for k in range(4)])
            _, S, _ = mp.svd_r(G)
            out["singular_values"] = list(S)
        return out
    return memo("chiral", w["legs"])(go)


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

def curve_file_offsets(w):
    """the curve files against RUN-009's polyline (12,000 wing steps, 160 centre steps) of Round 13's curve and of the
    search candidate's mixture (its meta's centres, widths and weight), centred at their mean: the largest difference,
    or None if the shapes differ"""
    P = p13(w)
    pts = r9.Wing(*P[:4], P[4]).points(12000, 160)
    base = float(np.abs(pts - pts.mean(0) - w["base_curve"]).max()) if pts.shape == w["base_curve"].shape else None
    B = w["summary"]["double_peak"]["search_candidate"]
    wl, b, d = B["parameters"][:3]
    m = B["meta"]

    class MixWing(r9.Wing):
        def __init__(self, comps):
            self.comps = [(aa, cc, ee, math.tanh(cc / ee), math.tanh((cc - 1) / ee)) for aa, cc, ee in comps]
            super().__init__(wl, b, d, comps[0][1], comps[0][2])

        def F(self, u):
            return sum(aa * (np.tanh((cc - u) / ee) - T1) / (T0 - T1) for aa, cc, ee, T0, T1 in self.comps)
    pts = MixWing([(1 - m["weight_2"], m["center_1"], m["epsilon_1"]), (m["weight_2"], m["center_2"], m["epsilon_2"])]).points(12000, 160)
    double = float(np.abs(pts - pts.mean(0) - w["double_curve"]).max()) if pts.shape == w["double_curve"].shape else None
    return base, double


def y1_curves(w):
    failures = []
    S = w["summary"]
    if S["reference_round13_scale"] != w["round13"]["summary"]["round13_candidate"]["common_scale"]:
        failures.append("the reference is not Round 13's common scale")
    base, double = curve_file_offsets(w)
    if base is None or base > 1e-14:
        failures.append("the base curve file is not Round 13's curve")
    B = S["double_peak"]["search_candidate"]
    wl, b, d, c, e, sep, lr, q = B["parameters"]
    a2 = 1 / (1 + math.exp(-q))
    m = B["meta"]
    if (m["center_1"], m["center_2"], m["epsilon_1"], m["epsilon_2"], m["weight_2"]) != (c - sep / 2, c + sep / 2, e * math.exp(-lr), e * math.exp(lr), a2):
        failures.append("the candidate's centres, widths and weight are not its parameters'")
    if double is None or double > 1e-14:
        failures.append("the double-peak curve file is not the search candidate's mixture")
    return failures


def y2_search(w):
    failures = []
    S = w["summary"]["double_peak"]
    E = search_exact(w)
    if None in E["values"]:
        return ["the search candidate's low branches are not where it says"]
    H = S["high_resolution_audit"]
    with mp.workdps(DPS):
        for k, (v, x) in enumerate(zip(H["branch_values"], E["values"])):
            dv = mp.mpf(repr(v)) - x
            if (k >= 2 and abs(dv) > mp.mpf("1e-13")) or (k < 2 and not mp.mpf("-1e-9") <= dv <= 0):
                failures.append(f"the high-resolution value of branch {k + 1} is not the exact one or a polyline's below it")
        if abs(mp.mpf(repr(H["scale"])) - min(E["values"])) > mp.mpf("1e-13") or H["scale"] != min(H["branch_values"]):
            failures.append("the high-resolution least is not the exact least")
        base = w["hr"]["high_resolution_base"]["scale"]
        if H["difference_from_high_resolution_base"] != H["scale"] - base \
                or H["difference_from_stored_round13"] != H["scale"] - w["summary"]["reference_round13_scale"]:
            failures.append("the high-resolution differences are not their figures' arithmetic")
        B = S["search_candidate"]
        if B["scale"] != min(B["branch_values"]):
            failures.append("the search candidate's scale is not its least branch")
        # the double-minus-single difference: the evaluator's single peak sits on a polyline's low branch, the double
        # peak's least on an exact cusp
        exact = min(E["values"]) - min(r14.branches(w13(w))["values"])
        if abs(mp.mpf(repr(H["difference_from_high_resolution_base"])) - exact) > mp.mpf("1e-12"):
            failures.append("the double-minus-single difference is the evaluator's, not the curves'")
    return failures


def y3_profile(w):
    failures = []
    rows = w["summary"]["double_peak"]["separation_profile"]
    for row in rows:
        if row["separation"] not in PROFILE_D:
            continue
        E = profile_exact(w, row["separation"], [row["phi1"], row["phi2"]])
        with mp.workdps(DPS):
            for k, (v, x) in enumerate(zip(row["branch_values"], E["values"])):
                if x is None:
                    failures.append(f"separation {row['separation']}: branch {k + 1} has no minimum near its phase")
                    continue
                dv = mp.mpf(repr(v)) - x
                if (k >= 2 and abs(dv) > mp.mpf("1e-12")) or (k < 2 and not mp.mpf("-1e-8") <= dv <= 0):
                    failures.append(f"separation {row['separation']}: branch {k + 1} is not the exact one or a polyline's below it")
    s0 = rows[0]["scale"]
    if any(r_["loss_from_zero_same_evaluator"] != r_["scale"] - s0 or r_["scale"] != min(r_["branch_values"]) for r_ in rows):
        failures.append("the losses are not the scales less the first, or a scale is not its least branch")
    fit = w["summary"]["double_peak"]["small_separation_fit"]
    # the fit is a quadratic in d^2 over the 21 profile points with d <= 0.01 (the window is inferred: it reproduces)
    small = [r_ for r_ in rows if r_["separation"] <= 0.01 + 1e-12]
    X = np.array([[r_["separation"] ** 4, r_["separation"] ** 2, 1.0] for r_ in small])
    y = np.array([r_["loss_from_zero_same_evaluator"] for r_ in small])
    coef = np.linalg.lstsq(X, y, rcond=None)[0]
    if np.abs(coef - np.array(fit["quartic_in_separation"])).max() > 1e-6 * np.abs(coef).max() or fit["quadratic_coefficient"] != fit["quartic_in_separation"][1]:
        failures.append("the small-separation fit is not the profile's own")
    # the low phases at the larger separations: stationary points, or a search interval's ends?
    stuck = [r_ for r_ in rows if abs(r_["phi1"] - 0.145) < 1e-8 or abs(r_["phi2"] - 0.137) < 1e-8]
    if stuck:
        failures.append("at the larger separations the profile's low phases are the ends of a search interval, not stationary points")
    return failures


def splitting_stated(w, h="1e-6"):
    """the low branches' second derivative in d at Round 13's stated curve, the five parameters and the phases fixed"""
    def go():
        R = r14.branches(w13(w))
        with mp.workdps(DPS + 10):
            p = [mp.mpf(repr(x)) for x in p13(w)]
            hh = mp.mpf(h)
            V = {i: sym_curve(w, p, i * hh) for i in (-1, 0, 1)}
            return [(V[1].value(SIG[k], R["phases"][k]) - 2 * V[0].value(SIG[k], R["phases"][k]) + V[-1].value(SIG[k], R["phases"][k])) / hh ** 2
                    for k in (0, 1)]
    return memo("splitstated", p13(w), w["legs"])(go)


def y4_splitting(w):
    failures = []
    # the profile's coefficient is the low branches' curvature along d with the five parameters fixed
    q = w["summary"]["double_peak"]["small_separation_fit"]["quadratic_coefficient"]
    m = min(splitting_stated(w)) / 2
    if abs(q / float(m) - 1) > 0.01:
        failures.append("the small-separation coefficient is not the low branches' curvature along d")
    # the first derivatives in (d, r) vanish by the exchange of the two halves - (d, r) and (-d, -r) are one curve - so
    # there is nothing to test there; the mixture must reproduce the solution
    B = splitting_block(w)
    with mp.workdps(DPS):
        if max(abs(x) for x in B["base_minus_solution"]) > mp.mpf("1e-20"):
            return failures + ["the mixture does not reproduce the solution"]
        if max(B["eigenvalues"]) > 0:
            failures.append("the Lagrangian has positive curvature in the splitting direction: the single peak is not a local maximum against splitting")
    return failures


def congruent_ok(rows, key):
    return all(r_["congruent_scale"] == min(r_["plus_scale"], r_["minus_scale"]) and abs(r_["chirality_gap"] - abs(r_["plus_scale"] - r_["minus_scale"])) < 1e-18
               for r_ in rows)


def y5_chirality(w):
    failures = []
    A = w["asym"]
    for prof, key in ((A["center_mode_profile"], "eta_center"), (A["width_mode_profile"], "eta_logwidth")):
        if not congruent_ok(prof, key):
            failures.append(f"the {key} profile's congruent scale is not the weaker handedness")
        by = {round(r_[key], 12): r_ for r_ in prof}
        if any(abs(r_["plus_scale"] - by[round(-r_[key], 12)]["minus_scale"]) > 1e-12 for r_ in prof if round(-r_[key], 12) in by):
            failures.append(f"the {key} profile's reflection does not exchange the handednesses")
    cen = A["random_census"]
    if not congruent_ok(cen, None) or len(cen) != w["summary"]["asymmetry"]["random_census_count"]:
        failures.append("the census' congruent scales or count are not as stated")
    ref = [c for c in A["refined_candidates"] if c["name"] != "zero"]
    best = max(ref, key=lambda c: c["congruent_scale"])
    if best != w["summary"]["asymmetry"]["best_nonzero_refined"] or A["best_refined"] != max(A["refined_candidates"], key=lambda c: c["congruent_scale"]):
        failures.append("the best refined candidates are not the largest")
    C = chirality(w)
    if min(C["singular_values"]) < mp.mpf("1e-6"):
        failures.append("the two modes' gradients do not have rank 2")
    S = w["summary"]["asymmetry"]
    for fit, prof, key in ((S["center_mode_fit"], A["center_mode_profile"], "eta_center"), (S["width_mode_fit"], A["width_mode_profile"], "eta_logwidth")):
        pts_ = [r_ for r_ in prof if abs(r_[key]) <= fit["window"] + 1e-12]
        X = np.array([[r_[key] ** 4, r_[key] ** 2, 1.0] for r_ in pts_])
        y = np.array([r_["congruent_scale"] - fit["zero_scale_same_evaluator"] for r_ in pts_])
        coef = np.linalg.lstsq(X, y, rcond=None)[0]
        if np.abs(coef - np.array(fit["quartic_in_mode"])).max() > 1e-6 * np.abs(coef).max() or fit["quadratic_coefficient"] != fit["quartic_in_mode"][1]:
            failures.append(f"the {key} fit is not the shipped profile's fit in its window")
    # the congruent scale's one-sided slopes at zero: min over the branches of the modes' gradients
    if all(max(abs(v) for v in C[m]) > mp.mpf("1e-8") for m in ("centre", "width")):
        failures.append("the congruent scale has a corner at zero in both modes; a quadratic coefficient does not describe it")
    return failures


def y6_phase(w):
    failures = []
    F = w["summary"]["full_phase"]
    M = r14.phase_map(w13(w))["minima"]
    with mp.workdps(60):
        off = []
        for m in F["base_local_minima"]:
            near = min(M, key=lambda e: abs(e["phase"] - mp.mpf(repr(m["phase"]))))
            off.append(abs(mp.mpf(repr(m["scale"])) - near["scale"]))
        if len(F["base_local_minima"]) != len(M):
            failures.append("the base list does not have the curve's number of minima")
        if max(off) > mp.mpf("1e-9"):
            failures.append("the listed minima are phase samples, not the curve's minima")
    if [m["phase"] for m in F["double_peak_local_minima"]] != [m["phase"] for m in F["base_local_minima"]]:
        failures.append("the double-peak list is not at the base list's phases")
    G = F["base_sampled_global"]
    # a double-precision evaluation of the cusp's closed form: within a few units in the last place
    with mp.workdps(60):
        g_off = abs(mp.mpf(repr(G["scale"])) - r14.branches(w13(w))["values"][3])
    if G["phase"] != F["base_local_minima"][0]["phase"] or g_off > mp.mpf("4e-16"):
        failures.append("the sampled global minimum is not the 270-degree cusp")
    if F["double_peak_local_minima"] != w["summary"]["double_peak"]["full_phase_local_minima"] or not F["no_new_low_branch"]:
        failures.append("the double-peak minima are not the summary's, or a new low branch is stated")
    return failures


def y7_copies(w):
    failures = []
    S = w["summary"]
    D = S["double_peak"]
    if w["search"]["best"] != D["search_candidate"]:
        failures.append("the search file's best is not the summary's candidate")
    H = w["hr"]
    if H["high_resolution_double_peak"] != D["high_resolution_audit"] or H["separation_profile"] != D["separation_profile"] \
            or H["small_separation_fit"] != D["small_separation_fit"] or H["reference_round13_scale"] != S["reference_round13_scale"]:
        failures.append("the high-resolution file is not the summary's")
    sep = [{"separation": float(r_["separation"]), "scale": float(r_["scale"]), "loss": float(r_["loss_from_zero_same_evaluator"]),
            "phi1": float(r_["phi1"]), "phi2": float(r_["phi2"])} for r_ in w["sep_csv"]]
    if sep != [{"separation": r_["separation"], "scale": r_["scale"], "loss": r_["loss_from_zero_same_evaluator"], "phi1": r_["phi1"], "phi2": r_["phi2"]}
               for r_ in D["separation_profile"]]:
        failures.append("the separation CSV is not the profile")
    for rows, prof, key in ((w["centre_csv"], w["asym"]["center_mode_profile"], "eta_center"), (w["width_csv"], w["asym"]["width_mode_profile"], "eta_logwidth")):
        if [{k: float(v) for k, v in r_.items()} for r_ in rows] != prof:
            failures.append(f"the {key} CSV is not the audit's profile")
    A = S["asymmetry"]
    if A["zero_candidate"] != w["asym"]["refined_candidates"][0] or A["best_nonzero_refined"] != [c for c in w["asym"]["refined_candidates"] if c["name"] == A["best_nonzero_refined"]["name"]][0]:
        failures.append("the summary's zero and best candidates are not the audit's")
    M = w["manifest"]
    if (M["reference_scale"], M["double_peak_difference"], M["version"]) != (S["reference_round13_scale"], D["high_resolution_audit"]["difference_from_high_resolution_base"], "0.14") \
            or abs(M["best_nonzero_asymmetry_difference"] - (A["best_nonzero_refined"]["congruent_scale"] - A["zero_candidate"]["congruent_scale"])) > 1e-18 \
            or sorted(M["files"]) != sorted(set(w["members"]) - {"MANIFEST.json"}):
        failures.append("the manifest is not the package's")
    return failures


def y8_documents(w):
    failures = []
    S = w["summary"]
    D, A = S["double_peak"], S["asymmetry"]
    rp = w["report"].replace("\r\n", "\n")
    rs = gate.squash(rp)
    one = lambda pat, text=rs: (re.search(pat, text) or [None, None])[1]
    rt = lambda txt, x: txt is not None and r10.rounds_to(txt, x)
    m = D["search_candidate"]["meta"]
    for pat, x in ((r"\|c_2-c_1\|\\approx([\d.e-]+)\.", m["center_2"] - m["center_1"]), (r"\\varepsilon_1=([\d.]+),", m["epsilon_1"]),
                   (r"\\varepsilon_2=([\d.]+)\.", m["epsilon_2"]), (r"a=([\d.]+)\.", m["weight_2"]),
                   (r"s_\{\\mathrm\{double\}\}-s_\{\\mathrm\{single\}\}=(-[\d.e-]+)\.", D["high_resolution_audit"]["difference_from_high_resolution_base"]),
                   (r"\\Deltas\(d\)\\approx(-[\d.]+)d\^2", D["small_separation_fit"]["quadratic_coefficient"]),
                   (r"s_\{\\mathrm\{zero\}\}=([\d.]+)\.", A["zero_candidate"]["congruent_scale"]),
                   (r"s_\{\\mathrm\{nonzero\}\}=([\d.]+)\.", A["best_nonzero_refined"]["congruent_scale"]),
                   (r"s_\{\\mathrm\{nonzero\}\}-s_\{\\mathrm\{zero\}\}=(-[\d.e-]+)\.", A["best_nonzero_refined"]["congruent_scale"] - A["zero_candidate"]["congruent_scale"])):
        if not rt(one(pat), x):
            failures.append(f"the report's {pat[:24]} is not the data's")
    # section 4's table: separations against the profile
    by = {round(r_["separation"], 6): r_ for r_ in D["separation_profile"]}
    rows = re.findall(r"^\|\s*([\d.]+)\s*\|\s*([\d.]+)\s*\|\s*(-?[\d.e+-]+)\s*\|$", rp, re.M)
    sec4 = [r_ for r_ in rows if round(float(r_[0]), 6) in by]
    if len(sec4) != 10 or any(not rt(s, by[round(float(d), 6)]["scale"]) or not rt(l, by[round(float(d), 6)]["loss_from_zero_same_evaluator"]) for d, s, l in sec4):
        failures.append("section 4's table is not the profile")
    # sections 6 and 7: the mode tables
    for prof, key, sec in ((w["asym"]["center_mode_profile"], "eta_center", "6"), (w["asym"]["width_mode_profile"], "eta_logwidth", "7")):
        byk = {round(r_[key], 6): r_ for r_ in prof}
        tab = re.findall(r"^\|\s*(-?[\d.]+)\s*\|\s*([\d.]+)\s*\|\s*([\d.]+)\s*\|\s*([\d.]+)\s*\|$", rp, re.M)
        mine = [t for t in tab if round(float(t[0]), 6) in byk and (len(t[0].split(".")[1]) == (4 if sec == "6" else 3))]
        if len(mine) != 7 or any(not (rt(p_, byk[round(float(e), 6)]["plus_scale"]) and rt(m_, byk[round(float(e), 6)]["minus_scale"])
                                      and rt(c_, byk[round(float(e), 6)]["congruent_scale"])) for e, p_, m_, c_ in mine):
            failures.append(f"section {sec}'s table is not the profile")
    if not rt(one(r"零點附近二次係數：\$\$(-[\d.]+)\.\$\$"), A["center_mode_fit"]["quadratic_coefficient"]) \
            or len(re.findall(r"零點附近二次係數：\$\$(-[\d.]+)\.\$\$", rs)) != 2 \
            or not rt(re.findall(r"零點附近二次係數：\$\$(-[\d.]+)\.\$\$", rs)[1], A["width_mode_fit"]["quadratic_coefficient"]):
        failures.append("sections 6 and 7's coefficients are not the data's")
    if f"抽樣 ${A['random_census_count']}$ 個二維偏移" not in rp:
        failures.append("section 8's census count is not the data's")
    if not rt(one(r"約下降：\$\$([\d.]+)\\times10\^\{-6\}"), -by[0.01]["loss_from_zero_same_evaluator"] * 1e6):
        failures.append("section 4's drop at d = 0.01 is not the profile's")
    # the README, the search audit, the handoff
    rd = w["readme"]
    for label, x in (("基準尺度：", S["reference_round13_scale"]), ("雙峰相對同求值器單峰：", D["high_resolution_audit"]["difference_from_high_resolution_base"]),
                     ("最佳雙峰峰距：約 ", m["center_2"] - m["center_1"])):
        mm = re.search(re.escape(label) + r"(-?[\d.e-]+)", rd)
        if not mm or not rt(mm.group(1), x):
            failures.append(f"the README's {label} is not the data's")
    mm = re.search(r"最佳非零手性候選減零偏移：\s*(-[\d.e-]+)", rd)
    if not mm or not rt(mm.group(1), A["best_nonzero_refined"]["congruent_scale"] - A["zero_candidate"]["congruent_scale"]):
        failures.append("the README's chirality difference is not the data's")
    an = w["audit_note"]
    for txt, x in ((r"高解析雙峰減單峰：\s*`([-\d.e]+)`", D["high_resolution_audit"]["difference_from_high_resolution_base"]),
                   (r"峰距：\s*`([-\d.e]+)`", m["center_2"] - m["center_1"]),
                   (r"小分離二次係數：\s*`([-\d.e]+)`", D["small_separation_fit"]["quadratic_coefficient"]),
                   (r"零偏移合同尺度：\s*`([-\d.e]+)`", A["zero_candidate"]["congruent_scale"]),
                   (r"最佳非零候選：\s*`([-\d.e]+)`", A["best_nonzero_refined"]["congruent_scale"]),
                   (r"非零減零：\s*`([-\d.e]+)`", A["best_nonzero_refined"]["congruent_scale"] - A["zero_candidate"]["congruent_scale"])):
        mm = re.search(txt, an)
        if not mm or abs(float(mm.group(1)) / x - 1) > 1e-15:
            failures.append(f"the search audit's {txt[:10]} is not the data's")
    if not rt(one(r"s_\{\\mathrm\{base\}\}=([\d.]+)\.\$\$", gate.squash(w["handoff"])), S["reference_round13_scale"]):
        failures.append("the handoff's base scale is not the data's")
    return failures


CHECKS = {"Y1": y1_curves, "Y2": y2_search, "Y3": y3_profile, "Y4": y4_splitting, "Y5": y5_chirality, "Y6": y6_phase,
          "Y7": y7_copies, "Y8": y8_documents}
DEPENDS = {
    "Y1": {"summary", "round13", "base_curve", "double_curve"},
    "Y2": {"summary", "hr", "round13", "legs"},
    "Y3": {"summary", "round13", "legs"},
    "Y4": {"summary", "round13", "legs"},
    "Y5": {"summary", "asym", "round13", "legs"},
    "Y6": {"summary", "round13", "legs"},
    "Y7": {"summary", "search", "hr", "asym", "sep_csv", "centre_csv", "width_csv", "manifest", "members"},
    "Y8": {"summary", "asym", "report", "readme", "audit_note", "handoff"},
}


# ---------------------------------------------------------------- the ascent: splitting with the five parameters re-solved

def split_F(w, x, d, h=mp.mpf("1e-12")):
    """Round 13's twelve equations for the curve split into two equal halves at separation d"""
    p, ph, s, mu = x[:5], [x[5], x[6], 2 * mp.pi / 3, 3 * mp.pi / 2], x[7], x[8:12]
    C = sym_curve(w, p, d)
    G = []
    for r in range(4):
        row = []
        for i in range(5):
            up = [p[k] + (h if k == i else 0) for k in range(5)]
            dn = [p[k] - (h if k == i else 0) for k in range(5)]
            row.append((sym_curve(w, up, d).value(SIG[r], ph[r]) - sym_curve(w, dn, d).value(SIG[r], ph[r])) / (2 * h))
        G.append(row)
    out = [C.value(SIG[r], ph[r]) - s for r in range(4)] + [C.deriv(SIG[r], ph[r]) for r in (0, 1)]
    return out + [mp.fsum(mu[r] * G[r][i] for r in range(4)) for i in range(5)] + [mp.fsum(mu) - 1]


SPLIT_D = ("0.004", "0.008", "0.016", "0.03")


def split_solutions(w):
    """for each separation in turn, the solution of the twelve equations by chord steps from the previous one (RUN-014's
    Jacobian at the single peak, updated by Broyden's rule), at 40 digits"""
    def go():
        K = kkt13(w)
        with mp.workdps(60):
            J = r14.kkt_J(w13(w), K["x"], "1e-15")["J"]
        out = []
        with mp.workdps(DPS):
            J = mp.matrix([[mp.mpf(J[a, b]) for b in range(12)] for a in range(12)])
            x = [mp.mpf(v) for v in K["x"]]
            for dtxt in SPLIT_D:
                d = mp.mpf(dtxt)
                Fx = split_F(w, x, d)
                steps = 0
                while max(abs(f) for f in Fx) > mp.mpf("1e-27") and steps < 40:
                    st = mp.lu_solve(J, mp.matrix(Fx))
                    xn = [x[k] - st[k] for k in range(12)]
                    Fn = split_F(w, xn, d)
                    dx = mp.matrix([xn[k] - x[k] for k in range(12)])
                    dF = mp.matrix([Fn[k] - Fx[k] for k in range(12)])
                    J = J + (dF - J * dx) * dx.T / (dx.T * dx)[0]
                    x, Fx = xn, Fn
                    steps += 1
                out.append({"d": dtxt, "x": list(x), "steps": steps, "residual": max(abs(f) for f in Fx)})
        return out
    return memo("splitsol", w["legs"])(go)


def mixture_boxes(P, K):
    """interval boxes for the points at u = j/K on both wings of a curve whose wings' layers are a mixture: as RUN-010's
    certificate_boxes, with M2 bounding |f''| <= theta'^2 + |theta''| through each layer's |F_i'| <= a_i/(eps_i den_i) and
    |F_i''| <= a_i (4/(3 sqrt 3))/(eps_i^2 den_i)"""
    iv.dps = 30
    wl, b, d = iv.mpf(P["wing_length"]), iv.mpf(P["beta"]), iv.mpf(P["delta"])
    comps = []
    for a, c, e in P["comps"]:
        a, c, e = iv.mpf(a), iv.mpf(c), iv.mpf(e)
        T0, T1 = tanh_iv(c / e), tanh_iv((c - 1) / e)
        comps.append((a, c, e, T1, T0 - T1))
    F1 = sum((a / (e * den) for a, c, e, T1, den in comps), iv.mpf(0))
    F2 = sum((a * (4 / (3 * iv.sqrt(3))) / (e ** 2 * den) for a, c, e, T1, den in comps), iv.mpf(0))
    M2 = iv.mpf(hi(d * F2 + (d * F1) ** 2))
    l0 = 1 - 2 * wl
    cs, sn = [], []
    for j in range(K + 1):
        u = iv.mpf(j) / K
        th = b + d * sum((a * (tanh_iv((c - u) / e) - T1) / den for a, c, e, T1, den in comps), iv.mpf(0))
        cs.append(iv.cos(th))
        sn.append(iv.sin(th))
    h = iv.mpf(1) / K
    step_err = iv.mpf(hi(h ** 3 * M2 / 12))
    X, Y = [iv.mpf(0)], [iv.mpf(0)]
    sx, sy = iv.mpf(0), iv.mpf(0)
    for j in range(K):
        sx += h * (cs[j] + cs[j + 1]) / 2
        sy += h * (sn[j] + sn[j + 1]) / 2
        e_ = hi(step_err * (j + 1))
        X.append(wl * (sx + iv.mpf([-e_, e_])))
        Y.append(wl * (sy + iv.mpf([-e_, e_])))
    x1, y1 = X[-1], Y[-1]
    tx, ty = iv.mpf(float((x1 + l0 / 2).mid)), iv.mpf(float((y1 / 2).mid))
    PX = [v - tx for v in X] + [2 * x1 + l0 - v - tx for v in X[::-1]]
    PY = [v - ty for v in Y] + [v - ty for v in Y[::-1]]
    return {"PX": PX, "PY": PY, "x1": x1, "y1": y1, "l0": l0, "M2": hi(M2)}


class boxes_of:
    """RUN-014's least bound reads its boxes from RUN-010's certificate_boxes; within this block it reads the mixture's"""

    def __init__(self, P):
        self.P = P

    def __enter__(self):
        self.saved = r10.certificate_boxes
        r10.certificate_boxes = lambda H, K: mixture_boxes(self.P, K)

    def __exit__(self, *a):
        r10.certificate_boxes = self.saved


def certify_split(w, sol):
    """the split curve given by its parameters to 30 digits: its least at least s(d) - 1.5e-10, in interval arithmetic"""
    with mp.workdps(DPS):
        x, d = sol["x"], mp.mpf(sol["d"])
        txt = lambda v: mp.nstr(v, 30)
        P = {"wing_length": txt(x[0]), "beta": txt(x[1]), "delta": txt(x[2]),
             "comps": [("0.5", txt(x[3] - d / 2), txt(x[4])), ("0.5", txt(x[3] + d / 2), txt(x[4]))]}
        T = x[7] - mp.mpf("1.5e-10")
    H = r14.hp_at(w13(w), [P["wing_length"], P["beta"], P["delta"], txt(x[3]), txt(x[4])])
    with boxes_of(P):
        out = r14.least_bounds(H, T)
    out["curve"] = P
    return out


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

def findings(w):
    out = {}
    S = w["summary"]
    D = S["double_peak"]
    R13 = r14.branches(w13(w))
    with mp.workdps(DPS):
        least13 = min(R13["values"])
        E = search_exact(w)
        out["curve_files_largest_difference"] = dict(zip(("base", "double_peak"), curve_file_offsets(w)))
        out["search_exact"] = {"values": [mp.nstr(v, 16) for v in E["values"]], "least": mp.nstr(min(E["values"]), 16),
                               "phases": [mp.nstr(v, 12) for v in E["phases"]]}
        out["search_stated_minus_exact"] = [mp.nstr(mp.mpf(repr(s)) - v, 3) for s, v in zip(D["search_candidate"]["branch_values"], E["values"])]
        out["high_resolution_stated_minus_exact"] = [mp.nstr(mp.mpf(repr(s)) - v, 3) for s, v in zip(D["high_resolution_audit"]["branch_values"], E["values"])]
        out["high_resolution_base_minus_exact"] = mp.nstr(mp.mpf(repr(w["hr"]["high_resolution_base"]["scale"])) - least13, 3)
        out["double_minus_single_exact"] = mp.nstr(min(E["values"]) - least13, 6)
        prof = []
        for row in D["separation_profile"]:
            if row["separation"] in PROFILE_D:
                Ex = profile_exact(w, row["separation"], [row["phi1"], row["phi2"]])
                prof.append({"d": row["separation"], "stated_minus_exact": [mp.nstr(mp.mpf(repr(s)) - v, 3) for s, v in zip(row["branch_values"], Ex["values"])]})
        out["profile_stated_minus_exact"] = prof
        merged = []
        for d in (0.028, 0.04):
            st = stationary_phases(w, d)
            C = sym_curve(w, [mp.mpf(repr(x)) for x in p13(w)], mp.mpf(repr(d)))
            row = [r_ for r_ in D["separation_profile"] if abs(r_["separation"] - d) < 1e-12][0]
            mins = [q for q, k in st if k == "minimum"]
            merged.append({"d": d, "low_stationary": [(mp.nstr(q, 10), k) for q, k in st],
                           "least_low": mp.nstr(min(C.value("L|p0|p2", q) for q in mins), 13),
                           "stated_scale": row["scale"],
                           "derivative_at_stated_phases": [mp.nstr(C.deriv("L|p0|p2", mp.mpf(repr(row[k]))), 4) for k in ("phi1", "phi2")]})
        out["profile_merged_low_minimum"] = merged
        B = splitting_block(w)
        out["splitting_branches"] = [{"m_dd": mp.nstr(m["dd"], 8), "m_dr": mp.nstr(m["dr"], 8), "m_rr": mp.nstr(m["rr"], 8)} for m in B["branches"]]
        out["splitting_lagrangian"] = [[mp.nstr(B["block"][i, j], 8) for j in range(2)] for i in range(2)]
        out["splitting_eigenvalues"] = [mp.nstr(v, 8) for v in B["eigenvalues"]]
        out["splitting_reproduces_solution"] = mp.nstr(max(abs(v) for v in B["base_minus_solution"]), 3)
        out["splitting_stated_curve_low_m_dd"] = [mp.nstr(v, 8) for v in splitting_stated(w)]
        Ch = chirality(w)
        out["chirality_gradients"] = {m: [mp.nstr(v, 8) for v in Ch[m]] for m in ("centre", "width")}
        out["chirality_cusp_phase_rates"] = {m: [mp.nstr(v, 4) for v in Ch[m + "_cusp_phase_rates"]] for m in ("centre", "width")}
        wp = w["asym"]["width_mode_profile"]
        z0 = [r_ for r_ in wp if r_["eta_logwidth"] == 0.0][0]["plus_scale"]
        top = max(wp, key=lambda r_: r_["plus_scale"])
        out["width_mode_one_handedness_rise"] = {"eta_logwidth": top["eta_logwidth"], "plus_scale_minus_zero": top["plus_scale"] - z0}
        out["chirality_singular_values"] = [mp.nstr(v, 6) for v in Ch["singular_values"]]
        out["congruent_one_sided_slope"] = {m: mp.nstr(-max(abs(min(Ch[m])), abs(max(Ch[m]))), 6) for m in ("centre", "width")}
        A = w["asym"]
        z = A["center_mode_profile"][15]["congruent_scale"]
        out["centre_profile_near_zero"] = [(r_["eta_center"], r_["congruent_scale"] - z) for r_ in A["center_mode_profile"][13:18]]
        out["width_profile_near_zero"] = [(r_["eta_logwidth"], r_["congruent_scale"] - z) for r_ in A["width_mode_profile"][13:18]]
        fits = {}
        for name, fit, prof, key in (("centre", S["asymmetry"]["center_mode_fit"], A["center_mode_profile"], "eta_center"),
                                     ("width", S["asymmetry"]["width_mode_fit"], A["width_mode_profile"], "eta_logwidth")):
            pts_ = [r_ for r_ in prof if abs(r_[key]) <= fit["window"] + 1e-12]
            X = np.array([[r_[key] ** 4, r_[key] ** 2, 1.0] for r_ in pts_])
            y = np.array([r_["congruent_scale"] - fit["zero_scale_same_evaluator"] for r_ in pts_])
            fits[name] = {"points": len(pts_), "refit": [float(f"{v:.8g}") for v in np.linalg.lstsq(X, y, rcond=None)[0]],
                          "stated": fit["quartic_in_mode"]}
        out["mode_fits_from_shipped_profiles"] = fits
    M = r14.phase_map(w13(w))["minima"]
    C13 = r10.cand(r14.hp(w13(w)), 60)
    with mp.workdps(60):
        rows = []
        for m in S["full_phase"]["base_local_minima"]:
            near = min(M, key=lambda e: abs(e["phase"] - mp.mpf(repr(m["phase"]))))
            rows.append({"phase": m["phase"], "minus_minimum_phase": mp.nstr(mp.mpf(repr(m["phase"])) - near["phase"], 3),
                         "minus_minimum_scale": mp.nstr(mp.mpf(repr(m["scale"])) - near["scale"], 3),
                         "minus_s_at_phase": mp.nstr(mp.mpf(repr(m["scale"])) - C13.s(mp.mpf(repr(m["phase"]))), 3)})
        out["base_minima_against_exact"] = rows
        R13 = w["round13"]["summary"]["full_phase_audit"]["other_minima"]
        out["round13_other_minima_against_list"] = [
            {"round13": [o["phase"], o["scale"]], "list": [m["phase"], m["scale"]]}
            for o in R13 for m in S["full_phase"]["base_local_minima"] if abs(m["scale"] - o["scale"]) < 1e-9]
    sols = split_solutions(w)
    x0 = kkt13(w)["x"]
    with mp.workdps(DPS):
        out["split_solutions"] = [{"d": s_["d"], "steps": s_["steps"], "residual": mp.nstr(s_["residual"], 3),
                                   "scale_minus_single_solution": mp.nstr(s_["x"][7] - x0[7], 8),
                                   "parameters": [mp.nstr(v, 20) for v in s_["x"][:5]], "pressures": [mp.nstr(v, 8) for v in s_["x"][8:12]]}
                                  for s_ in sols]
    cert = certify_split(w, sols[-1])
    lo13, hi13 = r14.s270_enclosure(r14.hp(w13(w)))
    with mp.workdps(40):
        cert["round13_least_at_most"] = dec(hi13, 18, up=True)
        cert["gain_over_round13_at_least"] = dec(cert["_lower"] - hi13, 18)
    del cert["_lower"], cert["_upper"]
    out["certificate"] = cert
    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(text, old, new, n=1):
        assert text.count(old) == n, old
        return text.replace(old, new)

    def both(fn_s, fn_h=None, fn_q=None):
        """a record planted alike in the summary and in its copies (the high-resolution file, the search file)"""
        out = {"summary": edit("summary", fn_s)}
        if fn_h:
            out["hr"] = edit("hr", fn_h)
        if fn_q:
            out["search"] = edit("search", fn_q)
        return out
    curve = base["base_curve"].copy()
    curve[5000, 1] += 1e-12
    dcurve = base["double_curve"].copy()
    dcurve[20000, 0] -= 1e-12
    eps2 = 0.03658402236371735
    m_eps = both(lambda S: (S["double_peak"]["search_candidate"]["meta"].__setitem__("epsilon_2", eps2),
                            S["double_peak"]["high_resolution_audit"]["meta"].__setitem__("epsilon_2", eps2)),
                 lambda H: H["high_resolution_double_peak"]["meta"].__setitem__("epsilon_2", eps2),
                 lambda Q: Q["best"]["meta"].__setitem__("epsilon_2", eps2))
    hr_low = both(lambda S: S["double_peak"]["high_resolution_audit"]["branch_values"].__setitem__(0, 0.9989143545181169),
                  lambda H: H["high_resolution_double_peak"]["branch_values"].__setitem__(0, 0.9989143545181169))
    s_scale = both(lambda S: S["double_peak"]["search_candidate"].__setitem__("scale", 0.998914343043348),
                   fn_q=lambda Q: Q["best"].__setitem__("scale", 0.998914343043348))
    k4 = [i for i, r_ in enumerate(base["summary"]["double_peak"]["separation_profile"]) if r_["separation"] == 0.004][0]
    cusp = both(lambda S: S["double_peak"]["separation_profile"][k4]["branch_values"].__setitem__(2, 0.9989143775650549),
                lambda H: H["separation_profile"][k4]["branch_values"].__setitem__(2, 0.9989143775650549))
    k1 = 1
    loss = both(lambda S: S["double_peak"]["separation_profile"][k1].__setitem__("loss_from_zero_same_evaluator", -8.438477028249736e-08),
                lambda H: H["separation_profile"][k1].__setitem__("loss_from_zero_same_evaluator", -8.438477028249736e-08))
    lcsv = copy.deepcopy(base["sep_csv"])
    lcsv[k1]["loss_from_zero_same_evaluator"] = "-8.438477028249736e-08"
    loss["sep_csv"] = lcsv

    def qfit(F):
        F["quadratic_coefficient"] = -0.03479446806580784
        F["quartic_in_separation"][1] = -0.03479446806580784
    fit = both(lambda S: qfit(S["double_peak"]["small_separation_fit"]), lambda H: qfit(H["small_separation_fit"]))
    fit["report"] = rrep(rp, "-0.033794468d^2", "-0.034794468d^2")
    fit["audit_note"] = rrep(base["audit_note"], "`-3.3794468065807839e-02`", "`-3.4794468065807839e-02`")
    # at eta_c = -0.002 the reflected handedness is the weaker, so the plant names the stronger
    kc = [i for i, r_ in enumerate(base["asym"]["center_mode_profile"]) if abs(r_["eta_center"] + 0.002) < 1e-12][0]
    assert base["asym"]["center_mode_profile"][kc]["plus_scale"] > base["asym"]["center_mode_profile"][kc]["minus_scale"]
    cong = edit("asym", lambda A: A["center_mode_profile"][kc].__setitem__("congruent_scale", A["center_mode_profile"][kc]["plus_scale"]))
    ccsv = copy.deepcopy(base["centre_csv"])
    ccsv[kc]["congruent_scale"] = ccsv[kc]["plus_scale"]
    census = edit("summary", lambda S: S["asymmetry"].__setitem__("random_census_count", 81))

    def wfit(S):
        S["asymmetry"]["width_mode_fit"]["quadratic_coefficient"] = -0.0031929038488534057
        S["asymmetry"]["width_mode_fit"]["quartic_in_mode"][1] = -0.0031929038488534057
    wf = edit("summary", wfit)
    wf_r = rrep(rp, "-0.003092904.", "-0.003192904.")
    glob = edit("summary", lambda S: (S["full_phase"]["base_sampled_global"].__setitem__("scale", 0.9989143432975433),
                                      S["full_phase"]["base_local_minima"][0].__setitem__("scale", 0.9989143432975433)))
    dphase = edit("summary", lambda S: (S["full_phase"]["double_peak_local_minima"][2].__setitem__("phase", 0.15493205957645976),
                                        S["double_peak"]["full_phase_local_minima"][2].__setitem__("phase", 0.15493205957645976)))
    scsv = copy.deepcopy(base["sep_csv"])
    scsv[3]["phi1"] = "0.15484021976978948"
    files = edit("manifest", lambda M: M["files"].remove("figures/03_center_asymmetry.png"))
    wcsv = copy.deepcopy(base["width_csv"])
    wcsv[4]["plus_scale"] = str(float(wcsv[4]["plus_scale"]) + 1e-15)
    readme = rrep(base["readme"], "最佳雙峰峰距：約 3.326e-06", "最佳雙峰峰距：約 3.336e-06")
    rep4 = rrep(rp, "| 0.0040 | 0.998913798823585 | -5.423e-07 |", "| 0.0040 | 0.998913798823585 | -5.433e-07 |")
    an = rrep(base["audit_note"], "`3.3255117917185828e-06`", "`3.3255117917285828e-06`")
    hand = rrep(base["handoff"], "s_{\\mathrm{base}}=0.99891434329749.", "s_{\\mathrm{base}}=0.99891434329759.")
    return [
        ("K0-unmodified", "control", {}, set()),
        ("K1-report-crlf", "control", {"report": rp.replace("\n", "\r\n")}, set()),
        ("K2-readme-trailing-spaces", "control", {"readme": base["readme"] + "  "}, set()),
        ("K3-asymmetry-copied", "control", {"asym": copy.deepcopy(base["asym"])}, set()),
        ("M1-triangle-leg", "mutation", {"legs": ("1.000000005", "1")}, {"Y2", "Y3", "Y6"}),
        ("M2-base-curve-point", "mutation", {"base_curve": curve}, {"Y1"}),
        ("M3-double-curve-point", "mutation", {"double_curve": dcurve}, {"Y1"}),
        ("M4-candidate-width", "mutation", m_eps, {"Y1"}),
        ("M5-high-resolution-low", "mutation", hr_low, {"Y2"}),
        ("M6-search-scale", "mutation", s_scale, {"Y2"}),
        ("M7-profile-cusp", "mutation", cusp, {"Y3"}),
        ("M8-profile-loss", "mutation", loss, {"Y3"}),
        ("M9-separation-coefficient", "mutation", fit, {"Y3", "Y4"}),
        ("M10-centre-congruent", "mutation", {"asym": cong, "centre_csv": ccsv}, {"Y5"}),
        ("M11-census-count", "mutation", {"summary": census, "report": rrep(rp, "抽樣 $80$ 個", "抽樣 $81$ 個")}, {"Y5"}),
        ("M12-width-fit", "mutation", {"summary": wf, "report": wf_r}, {"Y5"}),
        ("M13-sampled-global", "mutation", {"summary": glob}, {"Y6"}),
        ("M14-double-list-phase", "mutation", {"summary": dphase}, {"Y6"}),
        ("M15-separation-csv", "mutation", {"sep_csv": scsv}, {"Y7"}),
        ("M16-manifest-files", "mutation", {"manifest": files}, {"Y7"}),
        ("M17-width-csv", "mutation", {"width_csv": wcsv}, {"Y7"}),
        ("M18-readme-separation", "mutation", {"readme": readme}, {"Y8"}),
        ("M19-report-loss", "mutation", {"report": rep4}, {"Y8"}),
        ("M20-audit-separation", "mutation", {"audit_note": an}, {"Y8"}),
        ("M21-handoff-base", "mutation", {"handoff": hand}, {"Y8"}),
    ]


def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("--drill", action="store_true")
    ap.add_argument("--blind", metavar="CID")
    ap.add_argument("--log", default=None)
    ap.add_argument("--no-findings", action="store_true")
    ap.add_argument("--dry", action="store_true")
    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[:12]:
            print("    ", x)
    log = {"run": "RUN-015",
           "subject": "Round 14 (v0.14): double-peak curvature, chirality breaking and the single peak's local stability",
           "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 or args.no_findings or args.dry) else findings(w)}
    if log["findings"]:
        print(json.dumps(log["findings"], indent=1, ensure_ascii=False, default=str)[:14000])
    if args.dry:
        drill = gate.run_drill(build_world, CHECKS, DEPENDS, [(c, k, ch, set()) for c, k, ch, _ in drill_cases(build_world())])
        for d in drill:
            print(d["case"], "tripped", d["tripped"])
        return
    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()