"""RUN-016 - Round 15 (v0.15): curvature function modes, the pressure-projected Hessian, the branch-opening window
and the eight-mode candidate.

Round 15 deforms the wing's curvature density by eight bounded modes, psi_n = sqrt(2n+1) P_n(t), n = 1..8, in the CDF
coordinate t of Round 13's curve: rho(u) = rho_p(u) exp(sum a_n psi_n(t(u))) / Z, rho_p the tanh layer at the curve's own
five parameters. At a re-solved base it projects the four branches' mode gradients with the pressures, compensates each
mode to first order by the five parameters (equal heights), builds a finite-difference Hessian of the pressure-weighted
Lagrangian along those paths (seven negative eigenvalues, one positive), walks a Newton direction with the five
parameters re-equalized, finds a ninth local minimum opening near m = 2.2 and taking over at m = 3.228, and stops at
m = 3.2275: the least 0.998914480716946, 1.374194610326e-7 above Round 13's, a fifth (latent) branch 3.0046e-9 above
it. A four-branch candidate at 0.998914692247 is excluded by a hidden branch.

Checks (each returns a list of failures; empty = pass)
  Z1 basis        eight modes named; their weighted means, norms and maxima are the Legendre modes' (0, 1, sqrt(2n+1),
                  to 1e-14 and 1e-10 relative); the Gram error
  Z2 curve        the curve file is the candidate's curve: the density on 96,000 steps with F and positions by
                  cumulative trapezoid, t from Round 13's curve and rho_p at the candidate's parameters, 320 centre steps,
                  centred at its mean, to 1e-10
  Z3 base         the base curve's least is the stated base scale (to 1e-12) and the Newton line's first row (to 1e-15);
                  the pressures sum to one; the residual is the pressures times the parameter gradients, which are the
                  branches' partial derivatives at the base (60 digits), to 2e-7
  Z4 modes        the mode gradients are the branches' derivatives in the modes at the base (40 digits, step 1e-8), to
                  2e-7; the weighted gradient is the pressures times them; the compensation is the minimum-norm
                  solution of the equal-height conditions
  Z5 Hessian      symmetric; its eigenpairs reproduce it and are the stated ones; one positive eigenvalue; the file is
                  the summary's; the Newton direction, of norm 0.04, lies in the span of the six most negative
                  eigen-directions (its components along the other two below 1e-15 of its norm)
  Z6 stationary   at RUN-014's solution of Round 13's system the pressure-weighted mode gradients vanish, as a Hessian
                  test of local maximality presupposes (fixed-text failure)
  Z7 candidate    the coefficients are 3.2275 times the Newton direction; at 40 digits the four branches and the latent
                  one are the highest-resolution audit's (to 1e-13), the least is the cusp at 2pi/3, the margin and the
                  gains are their figures'; the audit's rows are their own arithmetic; the nine listed minima are the
                  curve's (phase to 1e-7, scale to 1e-13), and a dense scan of the curve file finds nine
  Z8 lines        the Newton line's and the ledger's norms are 0.04 m, their gains the scales less the base; the Newton
                  line is the first-order equal-height path p0 + Q a, a = m a_N: at m = 1 and 2 its least, the
                  270-degree cusp, is the curve's there (to 1e-13); the ninth minimum opens at m = 2.2; the ledger's
                  four-branch and full scales agree (to 1e-13) below the takeover window and part past it
  Z9 profiles     each single-mode profile row's minimum, pressure scale and spread are its branch values'; four rows'
                  cusp values are the curve's along the stated first-order equal-height path (to 1e-13)
  Z10 documents   the copies, the manifest, the README, the report's tables (sections 1, 3, 6 to 8), the audit note and
                  the handoff state the data's; the printed digits of the least, the gain and the latent margin against
                  the curve's (fixed-text failure)

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/run016_round15_functional_modes.py [--drill [--log PATH] | --blind Z# | --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 numpy.polynomial import legendre as npleg

from corpus import TREE, sha256, source_bytes
import gate
from run001_corpus_identity import restore_escapes
from run004_round03_symmetric_links import load
import run010_round09_precision_cusp as r10
import run014_round13_smooth_kkt as r14
from run010_round09_precision_cusp import dec, hi, iv, lo, tanh_iv

RID = "round-v0.15"
LOG = TREE / "data" / "gate-logs" / "run016-round15-functional-modes.json"
DPS = 40
NG = 24
SIG = r14.SIG
NAMES = ("min_x", "min_y", "max_d")
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_15_REPORT.md"), "audit_note": md("reports/ROUND_15_BRANCH_OPENING_AUDIT.md"),
        "handoff": md("reports/ROUND_16_FIVE_BRANCH_HANDOFF.md"), "source": f["src/functional_mode_model.py"].decode("utf-8"),
        "manifest": js("MANIFEST.json"), "members": sorted(f),
        "summary": js("data/round15_summary.json"), "final": js("data/final_candidate.json"), "profiles": js("data/mode_profiles.json"),
        "curve": arr("data/final_candidate_curve.csv"), "newton": rows("data/newton_line_full_phase.csv"),
        "hessian_csv": arr("data/pressure_hessian.csv"), "eigvec_csv": arr("data/pressure_hessian_eigenvectors.csv"),
        "ledger": rows("data/reoptimized_multiplier_ledger.csv"),
        "round13": r14.build_world(),
        "legs": ("1", "1"),
    }


def w13(w):
    return dict(w["round13"], legs=w["legs"])


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


def mode_names(w):
    return w["summary"]["mode_basis"]["names"]


# ---------------------------------------------------------------- the mode-deformed curve, exactly

def _gl(n):
    with mp.workdps(DPS + 10):
        xs, ws = [], []
        for i in range(1, n + 1):
            x = mp.cos(mp.pi * (i - mp.mpf(1) / 4) / (n + mp.mpf(1) / 2))
            for _ in range(100):
                p0, p1 = mp.mpf(1), x
                for k in range(2, n + 1):
                    p0, p1 = p1, ((2 * k - 1) * x * p1 - (k - 1) * p0) / k
                dp = n * (x * p1 - p0) / (x * x - 1)
                dx = p1 / dp
                x -= dx
                if abs(dx) < mp.mpf(10) ** (-(DPS + 8)):
                    break
            xs.append(x)
            ws.append(2 / ((1 - x * x) * dp * dp))
        return xs, ws


GX, GW = _gl(NG)


def legendre_all(x, n):
    P = [mp.mpf(1), x]
    for k in range(2, n + 1):
        P.append(((2 * k - 1) * x * P[-1] - (k - 1) * P[-2]) / k)
    return P


with mp.workdps(DPS + 10):
    PL = [legendre_all(x, NG) for x in GX]


class ModeWing:
    """the left wing: rho(u) = rho_p(u) exp(sum a_n sqrt(2n+1) P_n(t13(u))) / Z, F = 1 - int_0^u rho, theta = beta + delta F;
    panels a quarter of the layer's width across it, on each a 24-point Gauss-Legendre rule and the Legendre
    interpolant of rho, integrated exactly"""

    def __init__(self, p, a, q13):
        with mp.workdps(DPS + 10):
            self.w, self.b, self.dl, self.c, self.e = (mp.mpf(x) for x in p)
            self.a = [mp.mpf(x) for x in a]
            self.c13, self.e13 = mp.mpf(q13[3]), mp.mpf(q13[4])
            self.T0, self.T1 = mp.tanh(self.c / self.e), mp.tanh((self.c - 1) / self.e)
            self.U0, self.U1 = mp.tanh(self.c13 / self.e13), mp.tanh((self.c13 - 1) / self.e13)
            e = min(self.e, self.e13)
            edges = {mp.mpf(0), mp.mpf(1)}
            for k in range(-48, 49):
                for cc in (self.c, self.c13):
                    x = cc + k * e / 4
                    if 0 < x < 1:
                        edges.add(x)
            self.edges = sorted(edges)
            self.panels = []
            I = mp.mpf(0)
            for lo_, hi_ in zip(self.edges, self.edges[1:]):
                h = (hi_ - lo_) / 2
                m = (hi_ + lo_) / 2
                vals = [self.rho_raw(m + h * x) for x in GX]
                coef = [(2 * k + 1) / mp.mpf(2) * mp.fsum(GW[i] * vals[i] * PL[i][k] for i in range(NG)) for k in range(NG)]
                self.panels.append({"lo": lo_, "hi": hi_, "h": h, "m": m, "coef": coef, "I0": I})
                I += h * mp.fsum(GW[i] * vals[i] for i in range(NG))
            self.Z = I
            self.pos0 = []
            X = Y = mp.mpf(0)
            for P in self.panels:
                self.pos0.append((X, Y))
                cx, cy = self._panel_pos(P, P["hi"])
                X += cx
                Y += cy
            self.end = (X, Y)

    def rho_raw(self, u):
        rp = mp.sech((self.c - u) / self.e) ** 2 / (self.e * (self.T0 - self.T1))
        t = 1 - 2 * (mp.tanh((self.c13 - u) / self.e13) - self.U1) / (self.U0 - self.U1)
        P = legendre_all(t, 8)
        return rp * mp.exp(mp.fsum(self.a[n - 1] * mp.sqrt(2 * n + 1) * P[n] for n in range(1, 9)))

    def _panel_I(self, P, u):
        x = (u - P["m"]) / P["h"]
        L = legendre_all(x, NG)
        tot = P["coef"][0] * (x + 1)
        for k in range(1, NG):
            tot += P["coef"][k] * (L[k + 1] - L[k - 1]) / (2 * k + 1)
        return P["h"] * tot

    def panel_of(self, u):
        lo_, hi_ = 0, len(self.panels) - 1
        while lo_ < hi_:
            mid = (lo_ + hi_) // 2
            if self.panels[mid]["hi"] < u:
                lo_ = mid + 1
            else:
                hi_ = mid
        return lo_

    def _panel_pos(self, P, u):
        h = (u - P["lo"]) / 2
        m = (u + P["lo"]) / 2
        cx = cy = mp.mpf(0)
        for x, wt in zip(GX, GW):
            th = self.b + self.dl * (1 - (P["I0"] + self._panel_I(P, m + h * x)) / self.Z)
            cx += wt * mp.cos(th)
            cy += wt * mp.sin(th)
        return self.w * h * cx, self.w * h * cy

    def pos(self, u):
        if u <= 0:
            return (mp.mpf(0), mp.mpf(0))
        k = self.panel_of(u)
        X0, Y0 = self.pos0[k]
        cx, cy = self._panel_pos(self.panels[k], u)
        return X0 + cx, Y0 + cy

    def u_of(self, f):
        target = (1 - f) * self.Z
        k = 0
        while k + 1 < len(self.panels) and self.panels[k + 1]["I0"] <= target:
            k += 1
        P = self.panels[k]
        a_, b_ = P["lo"], P["hi"]
        for _ in range(200):
            m = (a_ + b_) / 2
            if P["I0"] + self._panel_I(P, m) < target:
                a_ = m
            else:
                b_ = m
            if b_ - a_ < mp.mpf(10) ** (-(DPS + 5)):
                break
        return (a_ + b_) / 2


class ModeCurve:
    def __init__(self, p, a, q13, legs):
        self.Wg = ModeWing(p, a, q13)
        with mp.workdps(DPS + 10):
            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)}
            x1, y1 = self.Wg.end
            self.l0 = 1 - 2 * self.Wg.w
            self.P = {"p0": (mp.mpf(0), mp.mpf(0)), "p1": (x1, y1), "p2": (x1 + self.l0, y1), "p3": (2 * x1 + self.l0, mp.mpf(0))}

    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]) if label == "L" else mp.atan2(n[1], -n[0])
        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 = (t - self.Wg.b) / self.Wg.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))
        u = mp.mpf(0) if f == 1 else (mp.mpf(1) if f == 0 else self.Wg.u_of(f))
        px, py = self.Wg.pos(u)
        return (px, py) if label == "L" else (2 * self.P["p1"][0] + self.l0 - px, py)

    def value(self, sig, ph):
        with mp.workdps(DPS):
            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):
        with mp.workdps(DPS):
            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


def stationary(C, sig, ph0, span="2e-3"):
    """a stationary point of s with a signature within span of ph0 (bisection on the derivative), its kind and value"""
    with mp.workdps(DPS):
        a, b = mp.mpf(ph0) - mp.mpf(span), mp.mpf(ph0) + mp.mpf(span)
        da, db = C.deriv(sig, a), C.deriv(sig, b)
        if (da > 0) == (db > 0):
            return None
        for _ in range(110):
            m = (a + b) / 2
            if (C.deriv(sig, m) > 0) == (da > 0):
                a = m
            else:
                b = m
        q = (a + b) / 2
        return q, C.value(sig, q)


_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 curve_of(w, p, a):
    """the mode-deformed curve for five parameters and eight coefficients (decimal strings or numbers)"""
    pk = [mp.nstr(mp.mpf(x), 45) if not isinstance(x, str) else x for x in p]
    ak = [mp.nstr(mp.mpf(x), 45) if not isinstance(x, str) else x for x in a]
    return memo("curve", pk, ak, p13(w), w["legs"])(lambda: ModeCurve(pk, ak, [repr(x) for x in p13(w)], w["legs"]))


def final_pa(w):
    FC = w["summary"]["final_candidate"]
    return [repr(x) for x in FC["parameters"]], [repr(FC["coefficients"][n]) for n in mode_names(w)]


ZERO = ["0"] * 8


def final_exact(w):
    """the final candidate at 40 digits: the four branches and the latent one, and the nine listed minima"""
    def go():
        p, a = final_pa(w)
        C = curve_of(w, p, a)
        L = w["summary"]["final_candidate"]["full_phase"]["local_minima"]
        out = {"listed": []}
        with mp.workdps(DPS):
            for m in L:
                ph, sg = m["phase"], m["signature"]
                if sg == "120deg":
                    out["listed"].append((2 * mp.pi / 3, C.value("p2|p1|p3", 2 * mp.pi / 3)))
                elif sg in ("p3|p3|p1", "p0|p3|p1") and abs(ph - 3 * math.pi / 2) < 1e-6:
                    out["listed"].append((3 * mp.pi / 2, C.value("p3|p3|p1", 3 * mp.pi / 2)))
                elif abs(ph - 5 * math.pi / 3) < 1e-6:
                    out["listed"].append((5 * mp.pi / 3, C.value(sg, 5 * mp.pi / 3)))
                elif sg == "p3|p1|p0":
                    # at pi/2 + alpha the support of min_x passes from the right wing's interior to its end, where s is
                    # continuously differentiable; the minimum is the smooth one just before the switch (RUN-011, RUN-014)
                    sw = mp.pi / 2 + C.Wg.b + C.Wg.dl
                    r = stationary(C, "R|p1|p0", sw - mp.mpf("5e-7"), "5e-7")
                    out["listed"].append(r if r else (None, None))
                else:
                    r = stationary(C, sg, ph)
                    out["listed"].append(r if r else (None, None))
            out["c120"], out["c270"] = out["listed"][0][1], out["listed"][1][1]
            out["B1"], out["B2"], out["latent"] = out["listed"][2][1], out["listed"][3][1], out["listed"][4][1]
            out["least"] = min(v for _, v in out["listed"] if v is not None)
        return out
    return memo("final", final_pa(w), w["summary"]["final_candidate"]["full_phase"]["local_minima"], w["legs"])(go)


def cand_low(C, ph0, span="2e-3"):
    """the smooth minimum of RUN-010's candidate with the supports L|p0|p2 within span of ph0, by bisection"""
    import run011_round10_phase_map as r11
    with mp.workdps(60):
        a, b = mp.mpf(ph0) - mp.mpf(span), mp.mpf(ph0) + mp.mpf(span)
        assert r11.deriv(C, SIG[0], a) < 0 < r11.deriv(C, SIG[0], b)
        for _ in range(150):
            m = (a + b) / 2
            if r11.deriv(C, SIG[0], m) < 0:
                a = m
            else:
                b = m
        return (a + b) / 2


def base_exact(w):
    """the base curve (no modes): its four branches at 60 digits, by RUN-010's candidate"""
    def go():
        P = [repr(x) for x in w["summary"]["reference"]["round15_base_parameters"]]
        H = r14.hp_at(w13(w), P)
        C = r10.cand(H, 60)
        with mp.workdps(60):
            ph = [cand_low(C, 0.1549), cand_low(C, 0.1239), 2 * mp.pi / 3, 3 * mp.pi / 2]
            vals = [r14.value(C, SIG[r], ph[r]) for r in range(4)]
            G = r14.grads(w13(w), [mp.mpf(x) for x in P], ph, "1e-20")
        return {"phases": ph, "values": vals, "least": min(vals), "G": G}
    return memo("base", w["summary"]["reference"]["round15_base_parameters"], w["legs"])(go)


def mode_gradients(w, p, phases, h="1e-8"):
    """the four branches' derivatives in the eight modes at fixed phases (central differences)"""
    def go():
        G = [[None] * 8 for _ in range(4)]
        with mp.workdps(DPS):
            hh = mp.mpf(h)
            for k in range(8):
                up = [h if j == k else "0" for j in range(8)]
                dn = ["-" + h if j == k else "0" for j in range(8)]
                Cu, Cd = curve_of(w, p, up), curve_of(w, p, dn)
                for r in range(4):
                    G[r][k] = (Cu.value(SIG[r], phases[r]) - Cd.value(SIG[r], phases[r])) / (2 * hh)
        return G
    return memo("modegrad", [mp.nstr(mp.mpf(x), 30) for x in p], [mp.nstr(x, 30) for x in phases], w["legs"])(go)


def base_mode_gradients(w):
    B = base_exact(w)
    return mode_gradients(w, [repr(x) for x in w["summary"]["reference"]["round15_base_parameters"]], B["phases"])


def kkt_mode_gradients(w):
    x = r14.kkt(w13(w))["x"]
    with mp.workdps(60):
        p = [mp.nstr(v, 45) for v in x[:5]]
        ph = [x[5], x[6], 2 * mp.pi / 3, 3 * mp.pi / 2]
    G = mode_gradients(w, p, ph)
    with mp.workdps(DPS):
        mu = x[8:12]
        return G, [mp.fsum(mu[r] * G[r][k] for r in range(4)) for k in range(8)]


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

def z1_basis(w):
    failures = []
    B = w["summary"]["mode_basis"]
    names = ["translation_like_P1", "width_split_P2", "skew_P3", "shoulder_P4", "tail_skew_P5", "triple_peak_P6", "wave_P7", "wave_P8"]
    if B["names"] != names or [m["name"] for m in B["mode_stats"]] != names:
        failures.append("the modes are not the eight named")
    for n, m in enumerate(B["mode_stats"], start=1):
        if abs(m["weighted_mean"]) > 1e-14 or abs(m["weighted_norm"] - 1) > 1e-14 or abs(m["max_abs"] / math.sqrt(2 * n + 1) - 1) > 1e-10:
            failures.append(f"mode {n}'s mean, norm or maximum is not the Legendre mode's")
    # the Gram matrix of sqrt(2n+1) P_n in t with weight dt/2 is the identity (a Gauss rule exact for degree 16)
    x, wt = npleg.leggauss(12)
    Pm = np.array([math.sqrt(2 * n + 1) * npleg.legval(x, [0] * n + [1]) for n in range(1, 9)])
    gram = (Pm * wt) @ Pm.T / 2
    if np.abs(gram - np.eye(8)).max() > 1e-14 or B["gram_max_error"] > 1e-15:
        failures.append("the Gram matrix is not the identity, or its stated error is not")
    return failures


def curve_file_offset(w, tlayer=None):
    """the curve file against the candidate's curve built on its grid; t from Round 13's curve, or from another
    layer's (w, beta, delta, c, eps) when given"""
    p, a = final_pa(w)
    wl, b, d, c, e = (float(x) for x in p)
    a = np.array([float(x) for x in a])
    n = 96000
    u = np.linspace(0, 1, n + 1)
    q = p13(w) if tlayer is None else [float(x) for x in tlayer]
    T0, T1 = math.tanh(c / e), math.tanh((c - 1) / e)
    U0, U1 = math.tanh(q[3] / q[4]), math.tanh((q[3] - 1) / q[4])
    t = 1 - 2 * (np.tanh((q[3] - u) / q[4]) - U1) / (U0 - U1)
    S = sum(a[k - 1] * math.sqrt(2 * k + 1) * npleg.legval(t, [0] * k + [1]) for k in range(1, 9))
    rho = np.exp(S) / (e * (T0 - T1)) / np.cosh((c - u) / e) ** 2
    trap = lambda y: np.concatenate([[0], np.cumsum((y[1:] + y[:-1]) / 2 / n)])
    I = trap(rho)
    th = b + d * (1 - I / I[-1])
    X, Y = wl * trap(np.cos(th)), wl * trap(np.sin(th))
    l0 = 1 - 2 * wl
    cen = np.stack([X[-1] + l0 * np.arange(1, 320) / 320, np.full(319, Y[-1])], 1)
    pts = np.concatenate([np.stack([X, Y], 1), cen, np.stack([2 * X[-1] + l0 - X[::-1], Y[::-1]], 1)])
    return float(np.abs(pts - pts.mean(0) - w["curve"]).max()) if pts.shape == w["curve"].shape else None


def z2_curve(w):
    off = curve_file_offset(w)
    return [] if off is not None and off <= 1e-10 else ["the curve file is not the candidate's curve"]


def z3_base(w):
    failures = []
    S = w["summary"]
    B = base_exact(w)
    PP = S["pressure_projection"]
    with mp.workdps(60):
        if abs(mp.mpf(repr(S["reference"]["round15_same_evaluator_base_scale"])) - B["least"]) > mp.mpf("1e-12"):
            failures.append("the base scale is not the base curve's least")
        if abs(mp.mpf(w["newton"][0]["global_scale"]) - B["least"]) > mp.mpf("1e-15") or float(w["newton"][0]["multiplier"]) != 0:
            failures.append("the Newton line's first row is not the base curve's least")
        if abs(sum(PP["pressures"]) - 1) > 1e-15:
            failures.append("the pressures do not sum to one")
        mu, G = PP["pressures"], PP["parameter_gradient_matrix"]
        if any(abs(sum(mu[r] * G[r][i] for r in range(4)) - PP["pressure_stationarity_residual"][i]) > 1e-15 for i in range(5)) \
                or abs(PP["pressure_stationarity_norm"] / math.sqrt(sum(x * x for x in PP["pressure_stationarity_residual"])) - 1) > 1e-12:
            failures.append("the stationarity residual is not the pressures times the parameter gradients")
        if max(abs(mp.mpf(repr(G[r][i])) - B["G"][r][i]) for r in range(4) for i in range(5)) > mp.mpf("2e-7"):
            failures.append("the parameter gradients are not the branches' partial derivatives at the base")
    return failures


def z4_modes(w):
    failures = []
    PP = w["summary"]["pressure_projection"]
    M = base_mode_gradients(w)
    with mp.workdps(DPS):
        if max(abs(mp.mpf(repr(PP["mode_branch_gradient_matrix"][r][k])) - M[r][k]) for r in range(4) for k in range(8)) > mp.mpf("2e-7"):
            failures.append("the mode gradients are not the branches' derivatives in the modes at the base")
    mu, Mg, G = np.array(PP["pressures"]), np.array(PP["mode_branch_gradient_matrix"]), np.array(PP["parameter_gradient_matrix"])
    if np.abs(mu @ Mg - np.array(PP["weighted_mode_gradient"])).max() > 1e-18:
        failures.append("the weighted mode gradient is not the pressures times the mode gradients")
    Q = -np.linalg.pinv(G[1:] - G[0]) @ (Mg[1:] - Mg[0])
    if np.abs(Q - np.array(PP["equalized_compensation_matrix"])).max() > 1e-14:
        failures.append("the compensation is not the minimum-norm equal-height solution")
    return failures


def newton_cosines(w):
    """the stated Newton direction against Newton steps of the stated Hessian and weighted gradient: the full step, the
    step on the seven negative eigen-directions, the step on the six most negative; its components along each
    eigen-direction, and their ratios to the Newton step's (a scaled Newton step has one ratio)"""
    Hh = w["summary"]["constrained_hessian_proxy"]
    H = np.array(Hh["matrix"])
    g = np.array(w["summary"]["pressure_projection"]["weighted_mode_gradient"])
    aN = np.array(Hh["newton_direction"])
    ev, V = np.linalg.eigh(H)
    c = V.T @ g
    step = -c / ev
    full = -np.linalg.solve(H, g)
    neg7 = V @ np.where(ev < 0, step, 0)
    neg6 = V @ np.where(np.arange(8) < 6, step, 0)
    cos = lambda x: float(abs(x @ aN) / np.linalg.norm(x) / np.linalg.norm(aN))
    comp = V.T @ aN
    ratios = comp[:6] / step[:6]
    return {"minus_H_inverse_g": cos(full), "seven_negative_directions": cos(neg7), "six_most_negative_directions": cos(neg6),
            "components_along_weakest_negative_and_positive": [float(f"{x:.3g}") for x in comp[-2:]],
            "ratio_to_newton_step_six_directions": [float(f"{x:.5g}") for x in ratios]}


def z5_hessian(w):
    failures = []
    Hh = w["summary"]["constrained_hessian_proxy"]
    H = np.array(Hh["matrix"])
    ev = np.linalg.eigvalsh(H)
    E = w["eigvec_csv"]
    if np.abs(H - H.T).max() > 0 or np.abs(ev - np.array(Hh["eigenvalues"])).max() > 1e-18 \
            or np.abs(E @ np.diag(Hh["eigenvalues"]) @ E.T - H).max() > 1e-18:
        failures.append("the Hessian is not symmetric, or its eigenpairs are not the stated ones")
    if Hh["positive_eigenvalue_count"] != int((ev > 0).sum()) or Hh["weakest_negative_eigenvalue"] != max(x for x in Hh["eigenvalues"] if x < 0):
        failures.append("the positive count or the weakest negative eigenvalue is not the spectrum's")
    if not np.array_equal(w["hessian_csv"], H):
        failures.append("the Hessian file is not the summary's")
    aN = np.array(Hh["newton_direction"])
    if abs(np.linalg.norm(aN) - Hh["newton_direction_norm"]) > 1e-12:
        failures.append("the Newton direction's norm is not stated")
    ev, V = np.linalg.eigh(H)
    comp = V.T @ aN
    if np.abs(comp[-2:]).max() > 1e-15 * np.linalg.norm(aN):
        failures.append("the Newton direction is not in the span of the six most negative eigen-directions")
    return failures


def z6_stationary(w):
    _, wg = kkt_mode_gradients(w)
    with mp.workdps(DPS):
        if max(abs(x) for x in wg) > mp.mpf("1e-10"):
            return ["at the five-parameter solution the pressure-weighted mode gradients do not vanish: it is not stationary in the modes"]
    return []


def z7_candidate(w):
    failures = []
    S = w["summary"]
    FC = S["final_candidate"]
    Hh = S["constrained_hessian_proxy"]
    a = np.array([FC["coefficients"][n] for n in mode_names(w)])
    if np.abs(a - FC["multiplier"] * np.array(Hh["newton_direction"])).max() > 1e-18 or abs(np.linalg.norm(a) - FC["coefficient_norm"]) > 1e-12:
        failures.append("the coefficients are not the multiplier times the Newton direction")
    E = final_exact(w)
    top = FC["resolution_audit"][-1]
    with mp.workdps(DPS):
        for k, v in zip(("B1", "B2", "c120", "c270"), top["four_branch_values"]):
            if E[k] is None or abs(mp.mpf(repr(v)) - E[k]) > mp.mpf("1e-13"):
                failures.append(f"the highest-resolution {k} value is not the curve's")
        if E["latent"] is None or abs(mp.mpf(repr(top["latent_scale"])) - E["latent"]) > mp.mpf("1e-13"):
            failures.append("the latent branch is not the curve's")
        if E["least"] != E["c120"] or abs(mp.mpf(repr(FC["highest_resolution_global_scale"])) - E["least"]) > mp.mpf("1e-13"):
            failures.append("the least is not the 120-degree cusp's")
        if abs(mp.mpf(repr(FC["highest_resolution_latent_margin"])) - (E["latent"] - E["least"])) > mp.mpf("1e-13"):
            failures.append("the latent margin is not the curve's")
        r13 = w["round13"]["summary"]["round13_candidate"]["common_scale"]
        for r in FC["resolution_audit"]:
            fb = r["four_branch_values"]
            if r["control_scale"] != min(fb) or r["control_spread"] != max(fb) - min(fb) or r["latent_minus_control"] != r["latent_scale"] - r["control_scale"] \
                    or r["global_scale"] != min(r["control_scale"], r["latent_scale"]) or r["gain_over_stored_round13"] != r["global_scale"] - r13:
                failures.append(f"the {r['u_count']}-point audit row is not its own arithmetic")
        if FC["highest_resolution_global_scale"] != top["global_scale"] or FC["highest_resolution_gain_over_round13"] != top["gain_over_stored_round13"] \
                or FC["highest_resolution_latent_margin"] != top["latent_minus_control"]:
            failures.append("the highest-resolution figures are not the last audit row's")
        for m, (q, v) in zip(FC["full_phase"]["local_minima"], E["listed"]):
            if q is None or abs(mp.mpf(repr(m["phase"])) - q) > mp.mpf("1e-7") or abs(mp.mpf(repr(m["scale"])) - v) > mp.mpf("1e-13"):
                failures.append(f"the listed minimum at {m['phase']:.6f} is not the curve's")
    if dense_minima(w) != FC["full_phase"]["sampled_local_minimum_count"]:
        failures.append("a dense scan of the curve does not find the stated number of local minima")
    return failures


def dense_minima(w, nph=65536):
    """the local minima of the scale on a phase grid, from the curve file's polyline (double precision): every support
    is a vertex of its convex hull"""
    def go():
        from scipy.spatial import ConvexHull
        P = w["curve"]
        V = P[ConvexHull(P).vertices]
        r3 = math.sqrt(3)
        A = (3 + 4 * r3) / 18 / 1.0048 * float(w["legs"][0])
        B = (4 + r3) / 6 / 1.0048 * float(w["legs"][1])
        ph = 2 * np.pi * np.arange(nph) / nph
        s = np.zeros(nph)
        for a0, a1 in ((1 / A, 1 / B), (-1 / A, 0.0), (0.0, -1 / B)):
            for k in range(0, nph, 1024):
                q = ph[k:k + 1024]
                n0, n1 = a0 * np.cos(q) + a1 * np.sin(q), -a0 * np.sin(q) + a1 * np.cos(q)
                s[k:k + 1024] += (np.outer(n0, V[:, 0]) + np.outer(n1, V[:, 1])).max(1)
        return int(((s < np.roll(s, 1)) & (s <= np.roll(s, -1))).sum())
    return memo("dense", w["legs"], float(w["curve"].sum()))(go)


def newton_row_cusp(w, m):
    """the 270-degree cusp of the curve at multiplier m on the first-order path: p = p0 + Q a, a = m a_N"""
    def go():
        S = w["summary"]
        aN = S["constrained_hessian_proxy"]["newton_direction"]
        Q = S["pressure_projection"]["equalized_compensation_matrix"]
        p0 = S["reference"]["round15_base_parameters"]
        a = [m * x for x in aN]
        p = [repr(float(p0[i] + sum(Q[i][k] * a[k] for k in range(8)))) for i in range(5)]
        C = curve_of(w, p, [repr(float(x)) for x in a])
        return C.value(SIG[3], 3 * mp.pi / 2)
    return memo("newtonrow", m, S_key(w), w["legs"])(go)


def S_key(w):
    S = w["summary"]
    return [S["constrained_hessian_proxy"]["newton_direction"], S["pressure_projection"]["equalized_compensation_matrix"],
            S["reference"]["round15_base_parameters"]]


def z8_lines(w):
    failures = []
    base = w["summary"]["reference"]["round15_same_evaluator_base_scale"]
    for rows, gk, nk in ((w["newton"], "global_scale", "local_minimum_count"), (w["ledger"], "full_global_scale", "full_local_minimum_count")):
        for r in rows:
            m = float(r["multiplier"])
            if abs(float(r["coefficient_norm"]) - 0.04 * m) > 1e-12 or abs(float(r["gain_over_base"]) - (float(r[gk]) - base)) > 1e-18:
                failures.append(f"line row m = {m}: the norm or gain is not its own arithmetic")
    with mp.workdps(DPS):
        for m in ("1.0", "2.0"):
            r = [x for x in w["newton"] if x["multiplier"] == m][0]
            if abs(float(r["global_phase"]) - 3 * math.pi / 2) > 1e-12 \
                    or abs(mp.mpf(r["global_scale"]) - newton_row_cusp(w, float(m))) > mp.mpf("1e-13"):
                failures.append(f"the Newton line at m = {m} is not the first-order path's 270-degree cusp")
    first9 = min(float(r["multiplier"]) for r in w["newton"] if r["local_minimum_count"] == "9")
    if not 2.0 < first9 <= 2.2 or any(r["local_minimum_count"] == "9" for r in w["newton"] if float(r["multiplier"]) < first9):
        failures.append("the ninth minimum does not open at m = 2.2 on the Newton line")
    win = w["summary"]["branch_opening"]["takeover_window"]
    for r in w["ledger"]:
        m, full, four = float(r["multiplier"]), float(r["full_global_scale"]), float(r["four_branch_scale"])
        if m <= win[0] and abs(full - four) > 1e-13:
            failures.append(f"ledger row m = {m}: below the window the full scale is not the four-branch scale")
        if m >= win[1] and not full < four - 1e-12:
            failures.append(f"ledger row m = {m}: past the window the full scale is not below the four-branch scale")
    return failures


PROFILE_CHECK = (("shoulder_P4", -0.04), ("shoulder_P4", 0.04), ("wave_P8", -0.04), ("wave_P8", 0.04))


def profile_cusps(w, name, amp):
    """a profile row's cusp values: the mode at its amplitude along the first-order equal-height path, the five
    parameters moved by the stated compensation (the rows' cusps are not the fixed-parameter curve's: they differ from
    it by the compensation's change, to first order)"""
    def go():
        k = mode_names(w).index(name)
        Q = w["summary"]["pressure_projection"]["equalized_compensation_matrix"]
        p0 = w["summary"]["reference"]["round15_base_parameters"]
        p = [repr(float(p0[i] + Q[i][k] * amp)) for i in range(5)]
        a = [repr(amp) if j == k else "0" for j in range(8)]
        C = curve_of(w, p, a)
        return [C.value(SIG[2], 2 * mp.pi / 3), C.value(SIG[3], 3 * mp.pi / 2)]
    return memo("profile", name, amp, w["summary"]["pressure_projection"]["equalized_compensation_matrix"], w["legs"])(go)


def z9_profiles(w):
    failures = []
    mu = w["summary"]["pressure_projection"]["pressures"]
    for name, rows in w["profiles"].items():
        for r in rows:
            b = r["branch_values"]
            if r["minimum_scale"] != min(b) or abs(r["branch_spread"] - (max(b) - min(b))) > 1e-18 \
                    or abs(r["pressure_scale"] - sum(m * x for m, x in zip(mu, b))) > 2e-16:
                failures.append(f"{name} at {r['amplitude']}: the minimum, spread or pressure scale is not its branches'")
    with mp.workdps(DPS):
        for name, amp in PROFILE_CHECK:
            r = [x for x in w["profiles"][name] if x["amplitude"] == amp][0]
            ex = profile_cusps(w, name, amp)
            if any(abs(mp.mpf(repr(r["branch_values"][2 + k])) - ex[k]) > mp.mpf("1e-13") for k in (0, 1)):
                failures.append(f"{name} at {amp}: the cusp values are not the curve's")
    return failures


def z10_documents(w):
    failures = []
    S = w["summary"]
    FC = S["final_candidate"]
    if w["final"] != FC or w["profiles"] is None or S["final_candidate_round13_gain"] != FC["highest_resolution_gain_over_round13"]:
        failures.append("the candidate file is not the summary's")
    M = w["manifest"]
    if (M["final_scale"], M["gain_over_round13"], M["latent_branch_margin"], M["version"]) != (
            FC["highest_resolution_global_scale"], FC["highest_resolution_gain_over_round13"], FC["highest_resolution_latent_margin"], "0.15") \
            or sorted(M["files"]) != sorted(set(w["members"]) - {"MANIFEST.json"}):
        failures.append("the manifest is not the package's")
    rp = w["report"].replace("\r\n", "\n")
    rs = gate.squash(rp)
    rt = lambda txt, x: txt is not None and r10.rounds_to(txt, x)
    one = lambda pat, text=rs: (re.search(pat, text) or [None, None])[1]
    # section 1's mode table
    rows = re.findall(r"^\|\s*(\d)\s*\|\s*`(\w+)`\s*\|\s*([\d.]+)\s*\|\s*(-?[\d.e+-]+)\s*\|$", rp, re.M)
    mg = S["pressure_projection"]["weighted_mode_gradient"]
    if len(rows) != 8 or any(nm != S["mode_basis"]["names"][int(k) - 1] or not rt(amx, S["mode_basis"]["mode_stats"][int(k) - 1]["max_abs"])
                             or not rt(g, mg[int(k) - 1]) for k, nm, amx, g in rows):
        failures.append("section 1's table is not the data's")
    if not rt(one(r"八模態Gram矩陣最大誤差：\$\$([\d.e-]+)\.\$\$"), S["mode_basis"]["gram_max_error"]):
        failures.append("section 1's Gram error is not the data's")
    ev = re.findall(r"^\|\s*(\d)\s*\|\s*(-?[\d.]+e-\d+)\s*\|$", rp, re.M)
    if len(ev) != 8 or any(not rt(v, S["constrained_hessian_proxy"]["eigenvalues"][int(k) - 1]) for k, v in ev):
        failures.append("section 3's eigenvalues are not the data's")
    inv = S["invalidated_candidate"]
    if not rt(one(r"表面尺度：\$\$([\d.]+)\.\$\$"), inv["four_branch_candidate_scale"]) or not rt(one(r"\\phi\\approx([\d.]+),"), inv["hidden_branch_phase"]) \
            or not rt(one(r"s\\approx([\d.]+)\.\$\$"), inv["full_phase_scale"]):
        failures.append("section 4's excluded candidate is not the data's")
    for k, sym in zip(range(5), ("w", "\\beta", "\\delta", "c", "\\varepsilon")):
        if not rt(one(r"\$\$" + re.escape(sym) + r"=([\d.]+)[,.]\$\$"), FC["parameters"][k]):
            failures.append(f"section 6's {sym} is not the data's")
    co = re.findall(r"^\|\s*`(\w+)`\s*\|\s*(-?[\d.]+)\s*\|$", rp, re.M)
    if len(co) != 8 or any(not rt(v, FC["coefficients"][n]) for n, v in co):
        failures.append("section 6's coefficients are not the data's")
    if not rt(one(r"m=([\d.]+)\.\$\$"), FC["multiplier"]) or not rt(one(r"\\\|a\\\|_2=([\d.]+)\.\$\$"), FC["coefficient_norm"]):
        failures.append("section 6's multiplier or norm is not the data's")
    res = re.findall(r"^\|\s*(\d+)\s*\|\s*([\d.]+)\s*\|\s*([\d.e-]+)\s*\|\s*([\d.e-]+)\s*\|$", rp, re.M)
    by = {r["u_count"]: r for r in FC["resolution_audit"]}
    if len(res) != 5 or any(not (rt(g, by[int(n)]["global_scale"]) and rt(l, by[int(n)]["latent_minus_control"]) and rt(d, by[int(n)]["gain_over_stored_round13"]))
                            for n, g, l, d in res):
        failures.append("section 7's table is not the data's")
    mins = re.findall(r"^\|\s*(\d)\s*\|\s*([\d.]+)\s*\|\s*([\d.]+)\s*\|\s*([\d.e+-]+)\s*\|\s*`([^`]+)`\s*\|$", rp, re.M)
    L = FC["full_phase"]["local_minima"]
    if len(mins) != 9 or any(not (rt(ph, L[int(k) - 1]["phase"]) and rt(sc, L[int(k) - 1]["scale"]) and rt(ab, L[int(k) - 1]["scale"] - L[0]["scale"]))
                             or sg != L[int(k) - 1]["signature"] for k, ph, sc, ab, sg in mins):
        failures.append("section 8's minima table is not the data's")
    for pat, x in ((r"\\boxed\{s_\{15\}=([\d.]+)\}", FC["highest_resolution_global_scale"]),
                   (r"\\boxed\{s_\{15\}-s_\{13\}=([\d.e-]+)\}", FC["highest_resolution_gain_over_round13"])):
        if not rt(one(pat), x):
            failures.append(f"the report's {pat[8:22]} is not the data's")
    rd = w["readme"]
    for label, x in (("第 15 輪有限模態候選：", FC["highest_resolution_global_scale"]), ("提升：", FC["highest_resolution_gain_over_round13"]),
                     ("第五分支餘量：", FC["highest_resolution_latent_margin"]), ("第 13 輪尺度：", w["round13"]["summary"]["round13_candidate"]["common_scale"])):
        m = re.search(re.escape(label) + r"([\d.e-]+)", rd)
        if not m or not rt(m.group(1), x):
            failures.append(f"the README's {label} is not the data's")
    if f"局部極小：{FC['full_phase']['sampled_local_minimum_count']}" not in rd or f"相位點：{FC['full_phase']['phase_count']}" not in rd:
        failures.append("the README's counts are not the data's")
    an = w["audit_note"]
    m = re.search(r"```json\n(.*?)\n```", an, re.S)
    if not m or json.loads(m.group(1)) != FC:
        failures.append("the audit note's candidate is not the data's")
    m = re.search(r"four-branch scale = ([\d.e+-]+)\nfull-phase hidden branch = ([\d.e+-]+)\nhidden phase = ([\d.e+-]+)", an)
    if not m or (float(m.group(1)), float(m.group(2)), float(m.group(3))) != (inv["four_branch_candidate_scale"], inv["full_phase_scale"], inv["hidden_branch_phase"]):
        failures.append("the audit note's excluded candidate is not the data's")
    hs = gate.squash(w["handoff"])
    if not rt(one(r"\\phi_5\\approx([\d.]+)\.\$\$", hs), S["branch_opening"]["latent_branch_phase"]) or not rt(one(r"s_\{15\}=([\d.]+)\.\$\$", hs), FC["highest_resolution_global_scale"]):
        failures.append("the handoff's figures are not the data's")
    # the printed digits against the curve's
    E = final_exact(w)
    with mp.workdps(DPS):
        r13 = min(r14.branches(w13(w))["values"])
        figs = ((one(r"\\boxed\{s_\{15\}=([\d.]+)\}"), E["least"]), (one(r"\\boxed\{s_\{15\}-s_\{13\}=([\d.e-]+)\}"), E["least"] - r13),
                (one(r"\\boxed\{([\d.e-]+)\.\}"), E["latent"] - E["least"]))
        if any(txt is None or not r10.rounds_to(txt, float(x)) for txt, x in figs):
            failures.append("the report prints the least, the gain and the latent margin past the digits the curve holds")
    return failures


CHECKS = {"Z1": z1_basis, "Z2": z2_curve, "Z3": z3_base, "Z4": z4_modes, "Z5": z5_hessian, "Z6": z6_stationary,
          "Z7": z7_candidate, "Z8": z8_lines, "Z9": z9_profiles, "Z10": z10_documents}
CURVE = {"summary", "round13", "legs"}
DEPENDS = {
    "Z1": {"summary"},
    "Z2": {"summary", "round13", "curve"},
    "Z3": CURVE | {"newton"},
    "Z4": CURVE,
    "Z5": {"summary", "hessian_csv", "eigvec_csv"},
    "Z6": {"round13", "legs"},
    "Z7": CURVE | {"curve"},
    "Z8": CURVE | {"newton", "ledger"},
    "Z9": CURVE | {"profiles"},
    "Z10": CURVE | {"final", "profiles", "manifest", "members", "report", "readme", "audit_note", "handoff"},
}


# ---------------------------------------------------------------- the certificate: interval boxes for the mode curve

TSN = 5


class TS:
    """a truncated Taylor series with interval coefficients: on an interval base point, c[k] encloses f^(k)/k! over it"""

    def __init__(self, c):
        self.c = list(c) + [iv.mpf(0)] * (TSN - len(c))

    def __add__(self, o):
        return TS([a + b for a, b in zip(self.c, o.c)]) if isinstance(o, TS) else TS([self.c[0] + o] + self.c[1:])
    __radd__ = __add__

    def __neg__(self):
        return TS([-a for a in self.c])

    def __sub__(self, o):
        return self + (-o)

    def __rsub__(self, o):
        return (-self) + o

    def __mul__(self, o):
        if not isinstance(o, TS):
            return TS([a * o for a in self.c])
        return TS([sum((self.c[i] * o.c[k - i] for i in range(k + 1)), iv.mpf(0)) for k in range(TSN)])
    __rmul__ = __mul__


def ts_exp(f):
    g = [iv.exp(f.c[0])]
    for k in range(1, TSN):
        g.append(sum((j * f.c[j] * g[k - j] for j in range(1, k + 1)), iv.mpf(0)) / k)
    return TS(g)


def ts_tanh(f):
    y = [tanh_iv(f.c[0])]
    q = [1 - y[0] * y[0]]
    for k in range(1, TSN):
        y.append(sum((j * f.c[j] * q[k - j] for j in range(1, k + 1)), iv.mpf(0)) / k)
        q.append(-sum((y[i] * y[k - i] for i in range(k + 1)), iv.mpf(0)))
    return TS(y)


def legendre_monomials(n):
    P = [[mp.mpf(1)], [mp.mpf(0), mp.mpf(1)]]
    for k in range(2, n + 1):
        a = [mp.mpf(0)] + [(2 * k - 1) * x / k for x in P[k - 1]]
        b = [(k - 1) * x / k for x in P[k - 2]] + [mp.mpf(0), mp.mpf(0)]
        P.append([x - y for x, y in zip(a, b)])
    return P


def mode_boxes(Pdec, coeffs, q13, K, block=64):
    """interval boxes for the points at u = j/K on both wings of the mode curve:
    I = int rho~ by Simpson's rule on each step, within h^5/2880 max|rho~''''| over its block (interval Taylor arithmetic
    of order 4 gives the derivative enclosures over each block of 64 steps); F = 1 - I/Z, theta = beta + delta F;
    X, Y by trapezoid steps, within w h^3/12 (theta'^2 + |theta''|) over the block, theta' = -delta rho~/Z and
    theta'' = -delta rho~'/Z"""
    iv.dps = 30
    with mp.workdps(50):
        mono = legendre_monomials(8)
        sc = [mp.mpf(0)] * 9
        for n in range(1, 9):
            for j, v in enumerate(mono[n]):
                sc[j] += mp.sqrt(2 * n + 1) * v * mp.mpf(coeffs[n - 1])
        S0 = [iv.mpf([x - mp.mpf("1e-40") * (1 + abs(x)), x + mp.mpf("1e-40") * (1 + abs(x))]) for x in sc]
    wl, b, d, c, e = (iv.mpf(x) for x in Pdec)
    c13, e13 = iv.mpf(q13[3]), iv.mpf(q13[4])
    T0, T1 = tanh_iv(c / e), tanh_iv((c - 1) / e)
    den = T0 - T1
    U0, U1 = tanh_iv(c13 / e13), tanh_iv((c13 - 1) / e13)
    den13 = U0 - U1

    def rho_ts(u):
        x = TS([u, iv.mpf(1)])
        th = ts_tanh(((-x) + c) * (1 / e))
        rp = (1 - th * th) * (1 / (e * den))
        th13 = ts_tanh(((-x) + c13) * (1 / e13))
        t = 1 - (th13 - U1) * (2 / den13)
        Sx = TS([S0[-1]])
        for co in reversed(S0[:-1]):
            Sx = Sx * t + co
        return rp * ts_exp(Sx)

    def rho_pt(u):
        z = tanh_iv((c - u) / e)
        t = 1 - 2 * (tanh_iv((c13 - u) / e13) - U1) / den13
        Sv = iv.mpf(0)
        for co in reversed(S0):
            Sv = Sv * t + co
        return (1 - z * z) / (e * den) * iv.exp(Sv)
    h = iv.mpf(1) / K
    mag = lambda x: max(abs(lo(x)), abs(hi(x)))
    blk = []
    for k in range(K // block):
        r = rho_ts(iv.mpf([mp.mpf(k * block) / K, mp.mpf((k + 1) * block) / K]))
        blk.append((mag(r.c[0]), mag(r.c[1]), hi(24 * iv.mpf(mag(r.c[4])))))
    I = [iv.mpf(0)]
    r_prev = rho_pt(iv.mpf(0))
    for j in range(K):
        rm = rho_pt(iv.mpf(2 * j + 1) / (2 * K))
        rn = rho_pt(iv.mpf(j + 1) / K)
        e_ = hi(h ** 5 / 2880 * blk[j // block][2])
        I.append(I[-1] + h * (r_prev + 4 * rm + rn) / 6 + iv.mpf([-e_, e_]))
        r_prev = rn
    Z = I[-1]
    Zlo = lo(Z)
    assert Zlo > 0, "the density's integral is not enclosed away from zero"
    cs = [iv.cos(b + d * (1 - Ij / Z)) for Ij in I]
    sn = [iv.sin(b + d * (1 - Ij / Z)) for Ij in I]
    X, Y = [iv.mpf(0)], [iv.mpf(0)]
    for j in range(K):
        r0, r1, _ = blk[j // block]
        e_ = hi(wl * h ** 3 / 12 * ((iv.mpf(hi(d)) * r0 / iv.mpf(Zlo)) ** 2 + iv.mpf(hi(d)) * r1 / iv.mpf(Zlo)))
        X.append(X[-1] + wl * h * (cs[j] + cs[j + 1]) / 2 + iv.mpf([-e_, e_]))
        Y.append(Y[-1] + wl * h * (sn[j] + sn[j + 1]) / 2 + iv.mpf([-e_, e_]))
    x1, y1 = X[-1], Y[-1]
    l0 = 1 - 2 * wl
    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, "Z_width": float(Z.delta)}


def certify_final(w):
    """the final candidate's least at least its exact least less 1.5e-10, in interval arithmetic: RUN-014's least bound,
    reading the mode curve's boxes in place of RUN-010's"""
    p, a = final_pa(w)
    q13 = [repr(x) for x in p13(w)]
    E = final_exact(w)
    with mp.workdps(40):
        T = E["least"] - mp.mpf("1.5e-10")
    H = r14.hp_at(w13(w), p)
    saved = r10.certificate_boxes
    r10.certificate_boxes = lambda H_, K: mode_boxes(p, a, q13, K)
    try:
        out = r14.least_bounds(H, T)
    finally:
        r10.certificate_boxes = saved
    return out


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

def findings(w):
    out = {}
    S = w["summary"]
    FC = S["final_candidate"]
    out["curve_file_largest_difference"] = curve_file_offset(w)
    E = final_exact(w)
    B = base_exact(w)
    with mp.workdps(DPS):
        r13 = min(r14.branches(w13(w))["values"])
        out["final_exact"] = {k: mp.nstr(E[k], 20) for k in ("c120", "c270", "B1", "B2", "latent", "least")}
        out["final_least_minus_stated"] = {"highest_resolution": mp.nstr(mp.mpf(repr(FC["highest_resolution_global_scale"])) - E["least"], 3),
                                           "by_resolution": {r["u_count"]: mp.nstr(mp.mpf(repr(r["global_scale"])) - E["least"], 3) for r in FC["resolution_audit"]}}
        out["final_branches_above_least"] = {k: mp.nstr(E[k] - E["least"], 6) for k in ("c270", "B1", "B2", "latent")}
        out["latent_margin_exact"] = mp.nstr(E["latent"] - E["least"], 8)
        out["gain_over_round13_exact"] = mp.nstr(E["least"] - r13, 12)
        out["listed_minima_exact"] = [{"phase": mp.nstr(q, 12), "scale": mp.nstr(v, 16)} if q is not None else None for q, v in E["listed"]]
        out["listed_minima_stated_minus_exact"] = [(mp.nstr(mp.mpf(repr(m["phase"])) - q, 3), mp.nstr(mp.mpf(repr(m["scale"])) - v, 3)) if q is not None else None
                                                  for m, (q, v) in zip(FC["full_phase"]["local_minima"], E["listed"])]
        out["dense_scan_local_minima"] = dense_minima(w)
        out["base_exact"] = {"values": [mp.nstr(v, 17) for v in B["values"]], "least": mp.nstr(B["least"], 17),
                             "stated_minus_exact": mp.nstr(mp.mpf(repr(S["reference"]["round15_same_evaluator_base_scale"])) - B["least"], 3),
                             "low_branches_above_cusps": mp.nstr(min(B["values"][:2]) - min(B["values"][2:]), 4)}
        Gs = S["pressure_projection"]["parameter_gradient_matrix"]
        out["base_parameter_gradients_stated_minus_exact_largest"] = {
            "low_branches": mp.nstr(max(abs(mp.mpf(repr(Gs[r][i])) - B["G"][r][i]) for r in (0, 1) for i in range(5)), 3),
            "cusps": mp.nstr(max(abs(mp.mpf(repr(Gs[r][i])) - B["G"][r][i]) for r in (2, 3) for i in range(5)), 3)}
        z0 = [x for x in w["profiles"]["shoulder_P4"] if x["amplitude"] == 0.0][0]["branch_values"]
        out["profile_zero_row_minus_base_exact"] = [mp.nstr(mp.mpf(repr(z0[r])) - B["values"][r], 3) for r in range(4)]
        M = base_mode_gradients(w)
        Ms = S["pressure_projection"]["mode_branch_gradient_matrix"]
        out["base_mode_gradients_stated_minus_exact_largest"] = mp.nstr(max(abs(mp.mpf(repr(Ms[r][k])) - M[r][k]) for r in range(4) for k in range(8)), 3)
        mu = S["pressure_projection"]["pressures"]
        out["base_weighted_mode_gradient_exact"] = [mp.nstr(mp.fsum(mp.mpf(repr(mu[r])) * M[r][k] for r in range(4)), 6) for k in range(8)]
        _, wg = kkt_mode_gradients(w)
        out["solution_weighted_mode_gradient"] = dict(zip(mode_names(w), [mp.nstr(x, 6) for x in wg]))
        x = r14.kkt(w13(w))["x"]
        out["base_minus_solution_parameters"] = [mp.nstr(mp.mpf(repr(v)) - x[k], 3) for k, v in enumerate(S["reference"]["round15_base_parameters"])]
        out["base_pressures_minus_solution"] = [mp.nstr(mp.mpf(repr(v)) - x[8 + k], 3) for k, v in enumerate(mu)]
    out["newton_direction_cosines"] = newton_cosines(w)
    led = [(float(r["multiplier"]), float(r["full_global_scale"]) - float(r["four_branch_scale"])) for r in w["ledger"]]
    out["ledger_full_minus_four_branch"] = [(m, float(f"{d:.4g}")) for m, d in led if abs(d) > 1e-13]
    with mp.workdps(DPS):
        out["profile_cusps_stated_minus_exact"] = {f"{n} {a}": [mp.nstr(mp.mpf(repr([x for x in w["profiles"][n] if x["amplitude"] == a][0]["branch_values"][2 + k])) - v, 3)
                                                               for k, v in enumerate(profile_cusps(w, n, a))] for n, a in PROFILE_CHECK}
    with mp.workdps(DPS):
        a = [mp.mpf(repr(FC["coefficients"][n])) for n in mode_names(w)]
        gs = [mp.mpf(repr(x)) for x in S["pressure_projection"]["weighted_mode_gradient"]]
        H = S["constrained_hessian_proxy"]["matrix"]
        out["first_order_term_at_solution"] = mp.nstr(mp.fsum(x * y for x, y in zip(wg, a)), 6)
        out["first_order_term_at_base"] = mp.nstr(mp.fsum(x * y for x, y in zip(gs, a)), 6)
        out["second_order_term_stated_hessian"] = mp.nstr(mp.fsum(a[i] * mp.mpf(repr(H[i][j])) * a[j] for i in range(8) for j in range(8)) / 2, 6)
        out["newton_line_rows_rebuilt"] = {m: mp.nstr(mp.mpf([x for x in w["newton"] if x["multiplier"] == m][0]["global_scale"]) - newton_row_cusp(w, float(m)), 3)
                                           for m in ("1.0", "2.0")}
    out["curve_file_t_from_other_layers"] = {k: curve_file_offset(w, v) for k, v in (
        ("round15_base", S["reference"]["round15_base_parameters"]),
        ("run014_solution", [float(x) for x in r14.kkt(w13(w))["x"][:5]]),
        ("candidate_own", FC["parameters"]))}
    cert = certify_final(w)
    lo_, up_ = cert.pop("_lower"), cert.pop("_upper")
    lo13, hi13 = r14.s270_enclosure(r14.hp(w13(w)))
    # a certified lower bound on Round 13's least (its least is the B2 low branch, below its s(3pi/2)): RUN-014's bound
    with mp.workdps(40):
        T13 = mp.mpf(repr(w["round13"]["summary"]["round13_candidate"]["common_scale"])) - mp.mpf("1.5e-10")
    c13 = r14.least_bounds(r14.hp(w13(w)), T13)
    with mp.workdps(40):
        cert["round13_least_at_most"] = dec(hi13, 18, up=True)
        cert["round13_least_at_least"] = dec(T13, 18) if c13["certified"] else None
        cert["gain_over_round13_at_least"] = dec(lo_ - hi13, 18)
        cert["gain_over_round13_at_most"] = dec(up_ - T13, 18, up=True) if c13["certified"] else None
    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 audit_json(fn):
        """the audit note's JSON block, planted as the candidate file is"""
        an = base["audit_note"]
        m = re.search(r"```json\n(.*?)\n```", an, re.S)
        FC = json.loads(m.group(1))
        fn(FC)
        return an[:m.start(1)] + json.dumps(FC, indent=2, ensure_ascii=False) + an[m.end(1):]

    def cand_plant(fn):
        """a record of the final candidate, in the summary, the candidate file and the audit note alike"""
        return {"summary": edit("summary", lambda S: fn(S["final_candidate"])), "final": edit("final", fn), "audit_note": audit_json(fn)}
    stats = edit("summary", lambda S: S["mode_basis"]["mode_stats"][2].__setitem__("max_abs", 2.6457513136551471))
    gram = edit("summary", lambda S: S["mode_basis"].__setitem__("gram_max_error", 2.220446049250313e-14))
    gram_r = rrep(rp, "2.220e-16.", "2.220e-14.")
    curve = base["curve"].copy()
    curve[50000, 0] += 1e-9
    rev = base["curve"][::-1].copy()
    newton = copy.deepcopy(base["newton"])
    newton[0]["global_scale"] = "0.9989143432437673"
    newton[0]["gain_over_base"] = repr(0.9989143432437673 - base["summary"]["reference"]["round15_same_evaluator_base_scale"])
    norm = edit("summary", lambda S: S["pressure_projection"].__setitem__("pressure_stationarity_norm", 1.3544257117e-07))

    def mode_g(S):
        PP = S["pressure_projection"]
        PP["mode_branch_gradient_matrix"][0][3] += 5e-7
        mu, M, G = np.array(PP["pressures"]), np.array(PP["mode_branch_gradient_matrix"]), np.array(PP["parameter_gradient_matrix"])
        PP["weighted_mode_gradient"] = [float(x) for x in mu @ M]
        PP["equalized_compensation_matrix"] = [[float(x) for x in row] for row in -np.linalg.pinv(G[1:] - G[0]) @ (M[1:] - M[0])]
    mg = edit("summary", mode_g)
    old_g4 = base["summary"]["pressure_projection"]["weighted_mode_gradient"][3]
    mg_r = rrep(rp, f"| {old_g4:.6e} |", f"| {mg['pressure_projection']['weighted_mode_gradient'][3]:.6e} |")
    comp = edit("summary", lambda S: S["pressure_projection"]["equalized_compensation_matrix"][2].__setitem__(4, -3.067015583242591e-04))
    hcsv = base["hessian_csv"].copy()
    hcsv[2, 5] += 1e-9
    pos = edit("summary", lambda S: S["constrained_hessian_proxy"].__setitem__("positive_eigenvalue_count", 2))
    coef = cand_plant(lambda FC: FC["coefficients"].__setitem__("shoulder_P4", -0.0752235430442821))
    coef["report"] = rrep(rp, "| `shoulder_P4` | -0.075223543043 |", "| `shoulder_P4` | -0.075223543044 |")
    spread = cand_plant(lambda FC: FC["resolution_audit"][0].__setitem__("control_spread", 1.0237176883148559e-09))
    minim = cand_plant(lambda FC: FC["full_phase"]["local_minima"][5].__setitem__("scale", 1.029348172728573))
    minim["report"] = rrep(rp, "| 1.736148508754 | 1.029348172727573 |", "| 1.736148508754 | 1.029348172728573 |")
    count = cand_plant(lambda FC: FC["full_phase"].__setitem__("sampled_local_minimum_count", 10))
    count["readme"] = rrep(base["readme"], "局部極小：9", "局部極小：10")
    nrow = copy.deepcopy(base["newton"])
    nrow[5]["coefficient_norm"] = "0.0041"
    led = copy.deepcopy(base["ledger"])
    k2 = [i for i, r in enumerate(led) if r["multiplier"] == "2.0"][0]
    led[k2]["four_branch_scale"] = "0.9989144327404801"
    win = edit("summary", lambda S: S["branch_opening"].__setitem__("takeover_window", [3.2275, 3.226]))

    def prof(P):
        r = [x for x in P["shoulder_P4"] if x["amplitude"] == 0.04][0]
        r["branch_values"][2] += 1e-12
        mu = base["summary"]["pressure_projection"]["pressures"]
        r["pressure_scale"] = sum(m * x for m, x in zip(mu, r["branch_values"]))
    pc = edit("profiles", prof)
    pa = edit("profiles", lambda P: P["width_split_P2"][3].__setitem__("minimum_scale", max(P["width_split_P2"][3]["branch_values"])))
    readme = rrep(base["readme"], "提升：1.374194610326e-07", "提升：1.374194610336e-07")
    rep7 = rrep(rp, "| 12001 | 0.998914480716965 | 3.005e-09 |", "| 12001 | 0.998914480716965 | 3.015e-09 |")
    an = rrep(base["audit_note"], "hidden phase = 1.3908150396590241e-01", "hidden phase = 1.3908150396690241e-01")
    hand = rrep(base["handoff"], "\\phi_5\\approx0.13842465.", "\\phi_5\\approx0.13852465.")
    files = edit("manifest", lambda M: M["files"].remove("figures/05_density_comparison.png"))
    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-profiles-copied", "control", {"profiles": copy.deepcopy(base["profiles"])}, set()),
        ("M1-triangle-leg", "mutation", {"legs": ("1.000000005", "1")}, {"Z3", "Z7", "Z8", "Z9"}),
        ("M2-mode-maximum", "mutation", {"summary": stats}, {"Z1"}),
        ("M3-gram-error", "mutation", {"summary": gram, "report": gram_r}, {"Z1"}),
        ("M4-curve-point", "mutation", {"curve": curve}, {"Z2"}),
        ("M5-curve-reversed", "mutation", {"curve": rev}, {"Z2"}),
        ("M6-newton-first-row", "mutation", {"newton": newton}, {"Z3"}),
        ("M7-stationarity-norm", "mutation", {"summary": norm}, {"Z3"}),
        ("M8-mode-gradient", "mutation", {"summary": mg, "report": mg_r}, {"Z4", "Z8", "Z9"}),
        ("M9-compensation", "mutation", {"summary": comp}, {"Z4", "Z8"}),
        ("M10-hessian-file", "mutation", {"hessian_csv": hcsv}, {"Z5"}),
        ("M11-positive-count", "mutation", {"summary": pos}, {"Z5"}),
        ("M12-coefficient", "mutation", coef, {"Z7"}),
        ("M13-audit-spread", "mutation", spread, {"Z7"}),
        ("M14-listed-minimum", "mutation", minim, {"Z7"}),
        ("M15-minimum-count", "mutation", count, {"Z7"}),
        ("M16-newton-norm", "mutation", {"newton": nrow}, {"Z8"}),
        ("M17-ledger-four-branch", "mutation", {"ledger": led}, {"Z8"}),
        ("M18-takeover-window", "mutation", {"summary": win}, {"Z8"}),
        ("M19-profile-cusp", "mutation", {"profiles": pc}, {"Z9"}),
        ("M20-profile-minimum", "mutation", {"profiles": pa}, {"Z9"}),
        ("M21-readme-gain", "mutation", {"readme": readme}, {"Z10"}),
        ("M22-report-resolution-row", "mutation", {"report": rep7}, {"Z10"}),
        ("M23-audit-hidden-phase", "mutation", {"audit_note": an}, {"Z10"}),
        ("M24-handoff-phase", "mutation", {"handoff": hand}, {"Z10"}),
        ("M25-manifest-files", "mutation", {"manifest": files}, {"Z10"}),
    ]


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-016",
           "subject": "Round 15 (v0.15): curvature function modes, the pressure-projected Hessian, the branch-opening window and the eight-mode candidate",
           "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)[:15000])
    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()