"""RUN-026 - Round 23: witness portfolio closure, strict incidence and residual-frontier routing.

Round 23 separates a strict cell x witness incidence - an edge only when a witness's whole
placement root is certified over the whole base cell - from quantitative partial edges, proves
that strict edges covering every necessity cell close the atlas (witness switching), writes the
choice of witnesses as a weighted set cover, rules that partial edges never compose, and runs a
12-cell x 5-witness pilot at lift depth 16: no edge is complete, so the strict cover is
infeasible; H19 has the least open volume on every cell, but at depth 22 B7 overtakes it on
three of the four hardest cells. Here the twelve cells are rebuilt from Round 19's segmentation,
every edge is rebuilt, the pilot's conventions are found, and the same cells are run with exact
cores and against each harmonic's limit - what its bound can reach on the cell.

Checks
  A0 anchors      the numbers and formulas below occur in Round 23 as written
  W1 theorems     sec. 3-11, 24-25: strict incidence, portfolio closure, the set cover and the
                  residual frontier, non-composition; checked on small instances
  W2 cells        sec. 12: the twelve cells are Round 19's twelve lowest outer brackets
  W3 matrix       sec. 13, 15-16: the sixty depth-16 edges; the aggregates; the routing
  W4 horizon      sec. 17-18: the depth-22 crossover on the four hardest cells; cell 0's ranking
  W5 exact        B7 with exact cores on every cell: where each tree closes
  W6 limits       each harmonic's least hull with each cell's cores: whether it can ever close
  W7 recorded     what the package records and ships
  P1 package      the matrix, the specs, the crystal and the text agree
  P2 package      VALIDATION against the archive; its checks from the data
  R1 registry     Round 23's judgments reviewed, first-used statuses defined

Usage: python code/run026_round23_portfolio_closure.py [--drill] [--blind CID]
"""
import argparse
import copy
import hashlib
import inspect
import itertools
import json
import math
from multiprocessing import Pool

import numpy as np
from scipy.optimize import minimize
from scipy.spatial import ConvexHull

import gate
import registry
import run008_round05_saturation_gate as r8
import run009_round06_witness_exchange as r9
import run011_round08_conditional_lift as r11
import run012_round09_adaptive_atlas as r12
import run013_round10_reference_dry_run as r13
import run019_round16_global_pilot as r19
import run020_round17_d3_cpg as r20
import run022_round19_necessity_lift as r22
import run023_round20_mixed_escalation as r23
import run025_round22_fourier_pareto as r25
from corpus import SOURCES, TREE, source_bytes, walk_sources
from run009_round06_witness_exchange import placed_disc
from run010_round07_minimizer_atlas import d_star, t_text
from run011_round08_conditional_lift import placed_vertices
from run012_round09_adaptive_atlas import core_pieces, pieces_support
from run013_round10_reference_dry_run import contact_points

PKG = "main.zip!/AMRAL_LUC_FC_Round_23_Package_v0.1.zip!/"
DOC = PKG + "AMRAL_Lebesgue_Finite_Closure_Round_23_Portfolio_Closure_v0.1.md"
CRYSTAL = PKG + "AMRAL_LUC_FC_State_Crystal_Round_23_v0.1.md"
MATRIX = PKG + "AMRAL_LUC_FC_Round_23_PORTFOLIO_MATRIX.json"
COVER = PKG + "AMRAL_LUC_FC_Round_23_SET_COVER_SPEC.json"
LONGRUN = PKG + "AMRAL_LUC_FC_Round_23_LONG_RUN_SPEC.json"
VALIDATION = PKG + "AMRAL_LUC_FC_Round_23_VALIDATION.json"
LOG = TREE / "data" / "gate-logs" / "run026-round23-portfolio-closure.json"
CACHE = TREE / "data" / "cache" / "run026-portfolio.json"
MANIFEST = TREE / "data" / "SOURCE_MANIFEST.v1.json"

T, TS, TM = 0.835, 0.8350002, 0.8365
WITNESSES = ("B7", "H11", "H13", "H17", "H19")
DEPTH, DEEP, CELLS, HARD = 16, 22, 12, 4
N_PILOT, N_OLD = 96, 40
EXACT_MAX, LIMIT_STARTS, SEED = 36, 14, 20260925

ANCHORS = {
    "W1": [r"J_K(H_{\mathcalF_3}(q))\geT\qquad\forallq\inC_j.", r"\text{strictedgeiscertificate-complete,notconfidence-complete}.",
           r"\Lambda(\mathcalF_3\cup\mathcalB)\geT", r"J_{K(j)}(H(q))\geT.", r"\min\sum_Kf_Ky_K+\sum_{j,K}c_{jK}x_{jK}.",
           r"\sum_Kx_{jK}\ge1\quad\forallj,", r"x_{jK}\leI_{jK}y_K.", r"\mathcalR_d=\left\{j", r"\max_KI_{jK}^{(d)}=0\right\}.",
           r"\mathcalR_d=\varnothing.", r"\text{PARTIAL}+\text{PARTIAL}\not\Rightarrow\text{STRICTCOVER}.",
           r"d^\star_{j;K,L}=\inf\{d:V_{jK}^{(d)}\leV_{jL}^{(d)}\}.", r"I_{jK}=0}$$不只是budget-zero"],
    "W2": [r"77\times5", r"$$\boxed{12}$$", r"$$\boxed{60}$$"],
    "W3": [r"d=16.", r"\mathcalR_{16}=\{1,\ldots,12\}.", r"|\mathcalR_{16}|=12.", r"$$2401.$$", r"9.42\times10^{-4}.",
           r"\approx1.93\times10^{-3}.", r"$$3749.5.$$", r"4.77\times10^{-4}.", r"7.47\times10^{-4}.", r"$$3644.2.$$",
           r"3.41\times10^{-4}.", r"6.38\times10^{-4}.", r"$$2362.5.$$", r"2.66\times10^{-4}.", r"7.83\times10^{-5}.",
           r"$$2286.3.$$", r"2.29\times10^{-4}.", r"4.96\times10^{-5}.", r"$$\boxed{12/12}$$",
           r"\text{CURRENT12-CELLSTRICTCOVERATDEPTH16:INFEASIBLE}"],
    "W4": [r"d=22.", r"cell0：$B_7$勝；", r"cell1：$B_7$勝；", r"cell2：$B_7$勝；", r"cell3：$H_{19}$勝。",
           r"B_7\text{在}3/4\text{hardestcells已於deephorizon反超}.", r"5.17\times10^{-5}", r"H_{19}:5.43\times10^{-5}",
           r"H_{17}:6.27\times10^{-5}", r"H_{13}:9.24\times10^{-5}", r"H_{11}:1.09\times10^{-4}."],
}

CLAIMS = {
    "segment_cells": 85,
    # the sixty edges against the matrix: B7 exact at 96 directions; the harmonics' pending counts
    "matrix_rebuild": {"edges_exact": 12, "pending_exact": 54, "harmonic_node_gap": [26, 64], "harmonic_pending_gap": [-2, 0],
                       "b7_all_exact": True, "b7_40_direction_matches": 0},
    "h19_deep_tolerance": 1e-12,
    # B7 with exact cores on each cell: [closing depth, nodes, open at the end, open at depth 16]
    "exact_closes": [[30, 18555, 0, 342], [31, 22009, 0, 361], [31, 26189, 0, 363], [34, 36193, 0, 379], [31, 21271, 0, 342],
                     [30, 18509, 0, 343], [33, 31209, 0, 358], [35, 45401, 0, 378], [32, 26439, 0, 364], [34, 41077, 0, 377],
                     [33, 31793, 0, 359], [30, 21895, 0, 365]],
    "limits_not_below_T": [],
    # the highest upper end, over the twelve cells, of each harmonic's least hull with the cell's cores
    "limit_worst": {"H11": 0.8342519199, "H13": 0.8339968527, "H17": 0.8334024191, "H19": 0.8333507439},
    "recorded_facts": {"ships_code": False, "states_directions": False, "states_root": False, "records_placements": False},
}

# ---------------------------------------------------------------- the lifts, sampled at N directions

def ring(N):
    th = 2 * np.pi * np.arange(N) / N
    return th, 0.5 * np.stack([np.cos(th), np.sin(th)], axis=1)


_CORES = {}


def cell_core_points(c, h, N):
    """the disc and the cell's two erosion cores at N directions"""
    kk = (tuple(c), tuple(h), N)
    if kk not in _CORES:
        k3, k5 = r19.cores(np.asarray(c, dtype=float), np.asarray(h, dtype=float))
        th, disc = ring(N)
        _CORES[kk] = np.vstack([disc] + [contact_points(p, th) for p in (k3, k5)])
    return _CORES[kk]


def root_half(name):
    """the translation root: t_7 for B7 (Round 07); for a harmonic the radius past which its
    nearest-reaching point, at least 1/2 - |a| from its centre, leaves B_{d*}"""
    if name == "B7":
        return t_text(7, d_star(TM))
    return d_star(TM) - 0.5 + r25.amp(int(name[1:]))


def box_delta(name, h):
    if name == "B7":
        return r23.delta_n(7, h)
    k = int(name[1:])
    return r25.delta("motion", k, r25.amp(k), h)


def _value(args):
    """a box's lower bound: the hull of the disc, the cell's cores and the witness's core, all
    at N directions, less the pilot's pad; pi/4 while the witness core is empty or illegal"""
    name, c, h, wc, wh, N = args
    th, _ = ring(N)
    dl = box_delta(name, wh)
    if name == "B7":
        rho = 1 - dl
        if rho <= r23.rn(7):
            return math.pi / 4
        core = contact_points(core_pieces(placed_vertices(7, wc[0], wc[1:3]), rho)[0], th)
    else:
        k = int(name[1:])
        a = r25.amp(k)
        if dl > r25.reserve(k, a):
            return math.pi / 4
        core = r25.harmonic_core(k, a, wc, dl, th)
    return float(ConvexHull(np.vstack([cell_core_points(c, h, N), core])).volume) - r19.PAD


def lift(pool, name, c, h, max_depth, N):
    """RUN-023's tree for any witness: root phi in [0, 2 pi/n), |x|, |y| <= the root half, split
    by the least child delta (ties to phi, then x); per depth (nodes, closed, open)"""
    n = 7 if name == "B7" else int(name[1:])
    t = root_half(name)
    level = [((math.pi / n, 0.0, 0.0), (math.pi / n, t, t))]
    per, nodes, closed = {}, 0, 0
    for d in range(max_depth + 1):
        vals = pool.map(_value, [(name, list(c), list(h), wc, wh, N) for wc, wh in level], chunksize=16)
        nodes += len(level)
        opened = [(wc, wh) for (wc, wh), v in zip(level, vals) if v < T]
        closed += len(level) - len(opened)
        per[d] = [nodes, closed, len(opened)]
        if not opened or d == max_depth:
            break
        level = []
        for wc, wh in opened:
            trial = []
            for ax in range(3):
                hh = list(wh)
                hh[ax] /= 2
                trial.append(box_delta(name, hh))
            ax = int(np.argmin(trial))
            hh = list(wh)
            hh[ax] /= 2
            lo, hi = list(wc), list(wc)
            lo[ax] -= hh[ax]
            hi[ax] += hh[ax]
            level += [(tuple(lo), tuple(hh)), (tuple(hi), tuple(hh))]
    return per


def root_volume(name):
    n = 7 if name == "B7" else int(name[1:])
    return 2 * math.pi / n * (2 * root_half(name)) ** 2


# ---------------------------------------------------------------- a harmonic's limit on a cell

_BASE = {}


def cores_points(c, h, N):
    kk = (tuple(c), tuple(h), N)
    if kk not in _BASE:
        th = 2 * np.pi * np.arange(N) / N
        k3, k5 = r19.cores(np.asarray(c, dtype=float), np.asarray(h, dtype=float))
        pts = np.vstack([contact_points(p, th) for p in (placed_disc(), k3, k5)])
        _BASE[kk] = pts[ConvexHull(pts).vertices]
    return _BASE[kk]


def limit_inner(c, h, k, g, N):
    th = 2 * np.pi * np.arange(N) / N
    spec = {"kind": "harmonic", "k": k, "a": r25.amp(k)}
    return float(ConvexHull(np.vstack([cores_points(c, h, N), r25.witness(spec)[1](g, th)])).volume)


def limit_bracket(c, h, k, g, N):
    spec = {"kind": "harmonic", "k": k, "a": r25.amp(k)}
    sup = r25.witness(spec)[0]
    k3, k5 = r19.cores(np.asarray(c, dtype=float), np.asarray(h, dtype=float))
    sups = [lambda t, p=p: pieces_support(p, t) for p in (placed_disc(), k3, k5)] + [lambda t: sup(g, t)]
    th = 2 * np.pi * np.arange(N) / N
    return [limit_inner(c, h, k, g, N), r25.support_polygon(np.unique(np.concatenate([th, r25.kinks(th, sups)])), sups)]


def _limit_search(args):
    c, h, k, g0, N = args
    r = minimize(lambda g: limit_inner(c, h, k, g, N), g0, method="Nelder-Mead",
                 options={"xatol": 1e-10, "fatol": 1e-13, "maxiter": 3000, "maxfev": 4000})
    return [float(r.fun), [float(v) for v in r.x]]


def limit(pool, c, h, k, rng):
    st = [[rng.uniform(0, 2 * math.pi / k), rng.uniform(-0.02, 0.02), rng.uniform(-0.02, 0.02)] for _ in range(LIMIT_STARTS)]
    coarse = sorted(pool.map(_limit_search, [(c, h, k, g, 6000) for g in st]))
    fine = sorted(pool.map(_limit_search, [(c, h, k, r[1], 24000) for r in coarse[:3]]))
    return {"fine": fine, "bracket": limit_bracket(c, h, k, fine[0][1], 64000)}


# ---------------------------------------------------------------- the computation

def twelve(pool):
    """Round 19's depth-40 cells below T on the tracked depth-31 cell (RUN-023's rebuild), with
    RUN-022's brackets; the twelve lowest outer brackets"""
    c31, h31, prefix = r23.branch(31)
    _, pend = r22.refine(c31, h31, 31, 40, TS, pool)
    ex = pool.map(r19.base_area, [x[1] for x in pend], chunksize=4)
    below = [(prefix + x[0], x[1], x[2], e) for x, e in zip(pend, ex) if e < T]
    brk = pool.map(r22._bracket, [x[1] for x in below])
    seg = [{"path": p, "center": [float(v) for v in cc], "half": [float(v) for v in hh], "exact": e, "inner": b[1], "outer": b[2]}
           for (p, cc, hh, e), b in zip(below, brk)]
    return seg, sorted(seg, key=lambda x: x["outer"])[:CELLS]


def compute():
    out = {}
    rng = np.random.default_rng(SEED)
    with Pool(14) as pool:
        seg, cells = twelve(pool)
        out["segment"] = [{k: v for k, v in s.items() if k in ("path", "exact", "outer")} for s in seg]
        out["cells"] = cells
        out["pilot"] = []
        for i, cell in enumerate(cells):
            c, h = cell["center"], cell["half"]
            row = {name: lift(pool, name, c, h, DEEP if (i < HARD and name in ("B7", "H19")) or (i == 0) else DEPTH, N_PILOT)
                   for name in WITNESSES}
            row["B7-40"] = lift(pool, "B7", c, h, DEPTH, N_OLD)
            row["B7-exact"] = r23.lift(c, h, 7, True, pool, EXACT_MAX)["per"]
            out["pilot"].append(row)
        out["limits"] = [{str(k): limit(pool, cell["center"], cell["half"], k, rng) for k in (11, 13, 17, 19)} for cell in cells]
    return out


def _code_key(files):
    digest = hashlib.sha256()
    for fn in (ring, cell_core_points, root_half, box_delta, _value, lift, root_volume, cores_points, limit_inner, limit_bracket,
               _limit_search, limit, twelve, compute):
        digest.update(inspect.getsource(fn).encode())
    for mod in (r8, r9, r11, r12, r13, r19, r20, r22, r23, r25):
        digest.update(open(mod.__file__, "rb").read())
    digest.update(files)
    digest.update(repr((WITNESSES, DEPTH, DEEP, CELLS, HARD, N_PILOT, N_OLD, EXACT_MAX, LIMIT_STARTS, SEED)).encode())
    return digest.hexdigest()


_COMPUTED = {}


def computed(key_bytes):
    key = _code_key(key_bytes)
    if key not in _COMPUTED:
        if CACHE.exists():
            cached = json.loads(CACHE.read_text(encoding="utf-8"))
            if cached.get("key") == key:
                _COMPUTED[key] = cached["results"]
        if key not in _COMPUTED:
            _COMPUTED[key] = json.loads(json.dumps(compute()))
            CACHE.parent.mkdir(parents=True, exist_ok=True)
            CACHE.write_text(json.dumps({"key": key, "results": _COMPUTED[key]}) + "\n", encoding="utf-8")
    return _COMPUTED[key]


# ---------------------------------------------------------------- checks

def near(a, b, tol):
    return abs(a - b) <= tol


def at(per, d):
    return r23.at(per, d)


def cover_optimum(I, f, c):
    """the strict weighted set cover by enumeration of the witness set: None when infeasible"""
    m, n = len(I), len(I[0])
    best = None
    for y in itertools.product((0, 1), repeat=n):
        cost = sum(fk for fk, yk in zip(f, y) if yk)
        for j in range(m):
            opts = [c[j][k] for k in range(n) if I[j][k] and y[k]]
            if not opts:
                cost = None
                break
            cost += min(opts)
        if cost is not None and (best is None or cost < best):
            best = cost
    return best


def w1_theorems(w):
    """sec. 5-11: portfolio closure is monotonicity of the hull; the set cover is feasible exactly
    when the residual frontier is empty; checked on random instances"""
    out = []
    rng = np.random.default_rng(SEED)
    worst = np.inf
    for _ in range(24):
        q = [rng.uniform(-0.01, 0.01), rng.uniform(-0.01, 0.01), rng.uniform(0, 1.2), rng.uniform(-0.02, 0.02), rng.uniform(-0.02, 0.02)]
        g1 = [rng.uniform(0, 1), rng.uniform(-0.02, 0.02), rng.uniform(-0.02, 0.02)]
        g2 = [rng.uniform(0, 1), rng.uniform(-0.02, 0.02), rng.uniform(-0.02, 0.02)]
        base = [p for body in r25.base_bodies(["cell", q]) for p in body]
        one = r19.exact_area(base + r9.placed_reuleaux(7, g1[0], g1[1:3]))
        two = r19.exact_area(base + r9.placed_reuleaux(11, g2[0], g2[1:3]))
        both = r19.exact_area(base + r9.placed_reuleaux(7, g1[0], g1[1:3]) + r9.placed_reuleaux(11, g2[0], g2[1:3]))
        worst = min(worst, both - max(one, two))
    frontier_ok = True
    for _ in range(200):
        I = (rng.random((6, 5)) < rng.uniform(0.05, 0.5)).astype(int).tolist()
        f, c = rng.uniform(1, 5, 5).tolist(), rng.uniform(0, 3, (6, 5)).tolist()
        feasible = cover_optimum(I, f, c) is not None
        frontier = [j for j in range(6) if max(I[j]) == 0]
        frontier_ok &= feasible == (not frontier)
    if worst < -1e-12 or not frontier_ok:
        out.append(f"monotonicity worst {worst:.2e}; frontier criterion {frontier_ok}")
    cov = w["cover"]
    if cov["strict_weighted_set_cover"]["constraints"] != ["for every cell j: sum_K x_jK >= 1", "x_jK <= I[j,K]*y_K"] \
            or "PARTIAL edges never count as cover" not in json.dumps(w["longrun"]):
        out.append("the set-cover spec's constraints or the long-run rule")
    w["w1"] = {"monotonicity_least_gain": worst, "frontier_criterion_holds": frontier_ok}
    return out


def w2_cells(w):
    """sec. 12: the twelve cells are the twelve lowest outer brackets of Round 19's depth-40 cells"""
    out, comp, rows = [], w["computed"], w["matrix"]["rows"]
    seg = comp["segment"]
    outer_below = sum(s["outer"] < T for s in seg)
    worst = 0.0
    for r, cell in zip(rows, comp["cells"]):
        worst = max(worst, abs(r["marker_outer"] - cell["outer"]))
        if r["path"] != cell["path"]:
            out.append(f"cell {r['cell_id']}: path {r['path']} against {cell['path']}")
    if [r["cell_id"] for r in rows] != list(range(CELLS)) or worst > 1e-13 \
            or outer_below != w["longrun"]["input"]["expected_outer-below-target_cells"] or len(seg) != w["segment_cells"]:
        out.append(f"the cells: {len(seg)} below T, {outer_below} with outer below T, marker difference {worst:.1e}")
    w["w2"] = {"cells_below_T": len(seg), "outer_below_T": outer_below, "marker_worst": worst}
    return out


def w3_matrix(w):
    """sec. 13, 15-16: the sixty depth-16 edges, their volumes, the aggregates, the routing"""
    out, comp, m = [], w["computed"], w["matrix"]
    exact_edges, pending_edges, node_gap, pend_gap, vol_worst = 0, 0, [], [], 0.0
    routes = []
    for r, row in zip(m["rows"], comp["pilot"]):
        vols = {}
        for name in WITNESSES:
            q = r["quantitative_edges"][name]
            nodes, _, opened = at(row[name], DEPTH)
            root = root_volume(name)
            vols[name] = opened / 2 ** DEPTH * root
            exact_edges += (nodes, opened) == (q["nodes"], q["pending_leaves"])
            pending_edges += opened == q["pending_leaves"]
            if name != "B7":
                node_gap.append(nodes - q["nodes"])
                pend_gap.append(opened - q["pending_leaves"])
            vol_worst = max(vol_worst, abs(q["unresolved_volume"] - q["pending_leaves"] / 2 ** DEPTH * q["root_volume"]) / q["unresolved_volume"],
                            abs(q["root_volume"] - root) / root)
            if q["status"] != "PARTIAL" or r["strict_incidence"][name] or opened == 0:
                out.append(f"cell {r['cell_id']} {name}: status {q['status']}, {opened} open at depth {DEPTH}")
        routes.append(min(vols, key=vols.get))
        if r["depth16_proxy_best_by_unresolved_volume"] != routes[-1]:
            out.append(f"cell {r['cell_id']}: the package routes to {r['depth16_proxy_best_by_unresolved_volume']}, this tree to {routes[-1]}")
    b7_exact = all((at(row["B7"], DEPTH)[0], at(row["B7"], DEPTH)[2]) == (r["quantitative_edges"]["B7"]["nodes"], r["quantitative_edges"]["B7"]["pending_leaves"])
                   for r, row in zip(m["rows"], comp["pilot"]))
    old = [at(row["B7-40"], DEPTH) for row in comp["pilot"]]
    old_matches = sum((o[0], o[2]) == (r["quantitative_edges"]["B7"]["nodes"], r["quantitative_edges"]["B7"]["pending_leaves"])
                      for o, r in zip(old, m["rows"]))
    got = {"edges_exact": exact_edges, "pending_exact": pending_edges, "harmonic_node_gap": [min(node_gap), max(node_gap)],
           "harmonic_pending_gap": [min(pend_gap), max(pend_gap)], "b7_all_exact": b7_exact, "b7_40_direction_matches": old_matches}
    if got != w["matrix_rebuild"] or vol_worst > 1e-12:
        out.append(f"the rebuild: {got}, volume worst {vol_worst:.1e}")
    ag = m["aggregate"]
    for name in WITNESSES:
        qs = [r["quantitative_edges"][name] for r in m["rows"]]
        mine = {"mean_nodes": sum(q["nodes"] for q in qs) / CELLS, "mean_pending": sum(q["pending_leaves"] for q in qs) / CELLS,
                "mean_unresolved_volume": sum(q["unresolved_volume"] for q in qs) / CELLS,
                "median_unresolved_volume": float(np.median([q["unresolved_volume"] for q in qs])),
                "mean_seconds": sum(q["seconds"] for q in qs) / CELLS, "strict_complete_cells": 0}
        if any(not near(mine[k], ag[name][k], 1e-12 * max(1.0, abs(ag[name][k]))) for k in mine):
            out.append(f"{name}: the aggregate against the rows")
    counts = {name: routes.count(name) for name in WITNESSES}
    if counts != m["proxy_winner_counts"] or m["strict_complete_edge_count"] != 0 or m["strict_set_cover_feasible"] is not False:
        out.append(f"routing {counts}")
    w["w3"] = {**got, "routes": counts, "b7_40_first": list(old[0])}
    return out


def w4_horizon(w):
    """sec. 17-18: depth 22 on the four hardest cells, and cell 0's whole portfolio"""
    out, comp, m = [], w["computed"], w["matrix"]
    rows = []
    for rec, row in zip(m["depth22_crossover_first4"], comp["pilot"]):
        v = {name: at(row[name], DEEP)[2] / 2 ** DEEP * root_volume(name) for name in ("B7", "H19")}
        winner = min(v, key=v.get)
        rows.append({"cell": rec["cell"], "B7": v["B7"], "H19": v["H19"], "winner": winner, "package": rec["winner"],
                     "B7_package": rec["B7"], "H19_package": rec["H19"]})
        if winner != rec["winner"] or not near(v["B7"], rec["B7"], 1e-12 * rec["B7"]) or not near(v["H19"], rec["H19"], w["h19_deep_tolerance"] * rec["H19"]):
            out.append(f"cell {rec['cell']}: B7 {v['B7']!r}, H19 {v['H19']!r}, package {rec['B7']!r}, {rec['H19']!r}")
    full = m["cell0_depth22_full_portfolio"]
    mine = {name: list(at(comp["pilot"][0][name], DEEP)) for name in WITNESSES}
    rank = sorted(WITNESSES, key=lambda n: mine[n][2] / 2 ** DEEP * root_volume(n))
    if rank != sorted(WITNESSES, key=lambda n: full[n]["unresolved_volume"]) or (mine["B7"][0], mine["B7"][2]) != (full["B7"]["nodes"], full["B7"]["pending"]):
        out.append(f"cell 0 at depth {DEEP}: {mine}")
    w["w4"] = {"crossover": rows, "cell0": mine, "cell0_rank": rank}
    return out


def w5_exact(w):
    """B7 with exact cores on every cell: every tree closes, and where"""
    out, comp = [], w["computed"]
    closes = []
    for row in comp["pilot"]:
        per = row["B7-exact"]
        last = max(per, key=int)
        closes.append([int(last), per[last][0], per[last][2], at(per, DEPTH)[2]])
    if closes != w["exact_closes"] or any(c[2] != 0 or c[3] == 0 for c in closes):
        out.append(f"B7 with exact cores: {closes}")
    w["w5"] = {"closes": closes}
    return out


def w6_limits(w):
    """each harmonic's least hull with each cell's erosion cores: an upper end below T means no
    lift depth can close that cell with that witness"""
    out, comp = [], w["computed"]
    table, above = [], []
    for i, lim in enumerate(comp["limits"]):
        row = {}
        for k in (11, 13, 17, 19):
            b, g = lim[str(k)]["bracket"], lim[str(k)]["fine"][0][1]
            inside = max(abs(g[1]), abs(g[2])) <= root_half(f"H{k}")      # phi counts modulo the period 2 pi/k
            row[f"H{k}"] = b
            if not b[1] < T or not inside:
                above.append((i, k))
        table.append(row)
    worst = {f"H{k}": max(r[f"H{k}"][1] for r in table) for k in (11, 13, 17, 19)}
    if above != w["limits_not_below_T"] or any(not near(worst[n], v, 5e-10) for n, v in w["limit_worst"].items()):
        out.append(f"limits not below T: {above}; highest upper ends {worst}")
    w["w6"] = {"table": table, "highest_upper_end": worst}
    return out


def w7_recorded(w):
    """what the package records and ships"""
    import re
    blob = json.dumps(w["matrix"]) + json.dumps(w["cover"]) + json.dumps(w["longrun"])
    facts = {"ships_code": any(p.endswith(".py") for p in w["package_files"]),
             "states_directions": bool(re.search(r"direction|sample|resolution", blob)),
             "states_root": "root_half" in blob or "translation_root" in blob,
             "records_placements": "placement" in blob and "phi" in blob}
    w["w7"] = facts
    return [] if facts == w["recorded_facts"] else [f"recorded: {facts}"]


def p1_package(w):
    out, text, m = [], gate.squash(w["text"]), w["matrix"]
    ag = m["aggregate"]
    for name, nodes, vol, margin in (("B7", "2401.", "9.42", "1.93"), ("H11", "3749.5.", "4.77", "7.47"), ("H13", "3644.2.", "3.41", "6.38"),
                                     ("H17", "2362.5.", "2.66", "7.83"), ("H19", "2286.3.", "2.29", "4.96")):
        a = ag[name]
        if not near(a["mean_nodes"], float(nodes.rstrip(".")), 0.5) or f"$${nodes}$$" not in text \
                or f"{a['mean_unresolved_volume']:.2e}"[:4] != vol or f"{a['forcing_margin_search']:.2e}"[:4] != margin:
            out.append(f"sec. 15, {name}: the aggregate against the text")
    full = m["cell0_depth22_full_portfolio"]
    for name, s in (("B7", "5.17"), ("H19", "5.43"), ("H17", "6.27"), ("H13", "9.24"), ("H11", "1.09")):
        if f"{full[name]['unresolved_volume']:.2e}"[:4] != s:
            out.append(f"sec. 18, {name}")
    winners = [r["winner"] for r in m["depth22_crossover_first4"]]
    if winners != ["B7", "B7", "B7", "H19"] or m["proxy_winner_counts"]["H19"] != 12 or m["lift_depth"] != DEPTH \
            or m["cell_count"] != CELLS or m["edge_count"] != CELLS * len(WITNESSES) or m["witnesses"] != list(WITNESSES):
        out.append("sec. 12-17 against the matrix")
    pr = w["cover"]["pilot_result"]
    if pr["cells"] != CELLS or pr["strict_edges"] != 0 or pr["strict_set_cover_feasible"] is not False or w["longrun"]["input"]["portfolio"] != list(WITNESSES):
        out.append("the specs' pilot result or portfolio")
    crystal = gate.squash(w["crystal"])
    for frag in ("B7~9.42e-4", "H11~4.77e-4", "H13~3.41e-4", "H17~2.66e-4", "H19~2.29e-4", "-B7beatsH19on3/4cells.", "-B7~5.17e-5",
                 "-H19~5.43e-5", "-H17~6.27e-5", "-H13~9.24e-5", "-H11~1.09e-4", r"\boxed{0}.", r"\boxed{\text{INFEASIBLE}}."):
        if frag not in crystal:
            out.append(f"crystal: {frag}")
    w["p1"] = {"winners": winners}
    return out


def p2_package(w):
    import re
    out, v = [], w["validation"]
    leaves = {row[0]: row for row in w["manifest"]["leaves"]}
    for key, leaf_id, text in (("round23", DOC, w["text"]), ("state_crystal", CRYSTAL, w["crystal"])):
        rec, (_, size, digest) = v[key], leaves[leaf_id]
        n = text.count("$$")
        if rec["sha256"] != digest or rec["bytes"] != size or not rec["utf8"] \
                or rec["double_dollar_count"] != n or rec["double_dollar_even"] != (n % 2 == 0) \
                or rec["contains_unicode_escape_sequence"] != bool(re.search(r"\\u[0-9a-fA-F]{4}", text)) \
                or rec["contains_bracket_math_delimiter"] != bool(re.search(r"\\[\[\]()]", text)):
            out.append(f"VALIDATION {key} against the archive and the document")
    for key, leaf_id in (("portfolio_matrix_sha256", MATRIX), ("set_cover_spec_sha256", COVER), ("long_run_spec_sha256", LONGRUN)):
        if v[key] != leaves[leaf_id][2]:
            out.append(f"VALIDATION {key} against the archive")
    m = w["matrix"]
    kc = {"matrix_has_12_cells": len(m["rows"]) == 12,
          "matrix_has_60_edges": sum(len(r["quantitative_edges"]) for r in m["rows"]) == 60,
          "strict_edge_count_zero": not any(any(r["strict_incidence"].values()) for r in m["rows"]),
          "strict_set_cover_infeasible": m["strict_set_cover_feasible"] is False,
          "H19_shallow_proxy_wins_all_12": all(r["depth16_proxy_best_by_unresolved_volume"] == "H19" for r in m["rows"]),
          "B7_depth22_wins_3_of_first4_vs_H19": sum(r["winner"] == "B7" for r in m["depth22_crossover_first4"]) == 3}
    if kc != v["checks"] or v["status"] != "PASS":
        out.append(f"VALIDATION checks {kc}")
    w["p2"] = kc
    return out


def load_inputs():
    blobs = {n: source_bytes(n) for n in SOURCES}
    want = {DOC, CRYSTAL, MATRIX, COVER, LONGRUN, VALIDATION}
    got = {leaf.id: d for leaf, d in walk_sources(blobs, keep=lambda l: l.id in want) if d is not None}
    return {"text": got[DOC].decode("utf-8"), "crystal": got[CRYSTAL].decode("utf-8"), "matrix": json.loads(got[MATRIX]),
            "cover": json.loads(got[COVER]), "longrun": json.loads(got[LONGRUN]), "validation": json.loads(got[VALIDATION]),
            "key_bytes": got[MATRIX], "manifest": json.loads(MANIFEST.read_text(encoding="utf-8")),
            "package_files": sorted(row[0][len(PKG):] for row in json.loads(MANIFEST.read_text(encoding="utf-8"))["leaves"]
                                    if row[0].startswith(PKG))}


R23_DEFINITIONS = {
    "PARTIAL": ("in-line", "sec. 4, 9: a lift tree that still has open leaves at its budget - a quantitative edge for the "
                           "scheduler, never an incidence"),
    "DISQUALIFIED": ("in-line", "sec. 25: an edge whose witness has a valid counterexample placement on the cell (Round 20's "
                                "negative filter) - out of single-witness cover until the semantics change"),
    "INFEASIBLE": ("in-line", "sec. 13: no strict cover exists at the budget - some cell has no complete edge"),
}
R23_REVIEWS = [
    ("CLOSED", "STRICT PORTFOLIO INCIDENCE", "analytic", "PASS",
     "an edge only for a whole certified root over the whole cell; theorem 5.1 is monotonicity of the hull, which adding a "
     "witness never lowers (W1)"),
    ("CLOSED", "SET-COVER LEGALITY", "analytic", "PASS",
     "the facility-location program is feasible exactly when the residual frontier is empty - by enumeration on 200 random "
     "instances (W1)"),
    ("DISQUALIFIED", "25. Counterexample effect", "analytic", "PASS",
     "Round 20's negative filter (RUN-023 W3): one placement below T bars a witness as sole closer"),
    ("INFEASIBLE", "CURRENT 12-CELL STRICT COVER AT DEPTH 16", "computation", "PASS",
     "no tree closes by depth 16, with the pilot's cores or exact ones (W3, W5); with exact cores B7 closes every cell by lift "
     "depth 35, so the cover exists at that budget - and no harmonic can supply an edge on these cells at any depth (W6)"),
    ("INFEASIBLE", "", "computation", "PASS", "sec. 0, the same verdict as sec. 13 (W3, W5)"),
    ("INFEASIBLE", "", "computation", "PASS", "sec. 13, the same verdict as the boxed one (W3, W5)"),
    ("STILL COMPUTE-DEFERRED", "a_{Leb}≥0.8350", "n/a", "n/a", "sec. 0's boxed summary; read from census v5 on (RUN-003)"),
]


def r1_registry(w):
    return registry.completeness(w["registry"], "R23", w["definitions"], w["reviews"])


CHECKS = {"A0": gate.anchors_check(ANCHORS), "W1": w1_theorems, "W2": w2_cells, "W3": w3_matrix, "W4": w4_horizon,
          "W5": w5_exact, "W6": w6_limits, "W7": w7_recorded, "P1": p1_package, "P2": p2_package, "R1": r1_registry}
DEPENDS = {"A0": {"text"}, "W1": {"cover", "longrun"},
           "W2": {"computed", "matrix", "longrun", "segment_cells"},
           "W3": {"computed", "matrix", "matrix_rebuild"},
           "W4": {"computed", "matrix", "h19_deep_tolerance"},
           "W5": {"computed", "exact_closes"},
           "W6": {"computed", "limits_not_below_T", "limit_worst"},
           "W7": {"matrix", "cover", "longrun", "package_files", "recorded_facts"},
           "P1": {"text", "matrix", "cover", "longrun", "crystal"},
           "P2": {"validation", "manifest", "text", "crystal", "matrix"},
           "R1": {"registry", "definitions", "reviews"}}


def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("--drill", action="store_true")
    ap.add_argument("--blind", metavar="CID")
    args = ap.parse_args()
    inputs = load_inputs()
    results_computed = computed(inputs["key_bytes"])

    def build():
        return {**{k: v for k, v in inputs.items() if k != "key_bytes"}, "computed": results_computed,
                "registry": registry.load(), "definitions": dict(R23_DEFINITIONS), "reviews": list(R23_REVIEWS),
                **copy.deepcopy(CLAIMS)}

    if args.blind:
        CHECKS[args.blind] = lambda w: []
        args.drill = True
    w = build()
    results = {c: fn(w) for c, fn in CHECKS.items()}
    for c, f in results.items():
        print(f"{c} {'PASS' if not f else 'FAIL'}", *(f"\n     {x}" for x in f[:6]))
    for key in ("w1", "w2", "w3", "w4", "w5", "w6", "w7", "p1", "p2"):
        print(key, json.dumps(w.get(key), default=float, ensure_ascii=False)[:1200])
    if not args.drill:
        return
    m, text = inputs["matrix"], inputs["text"]

    def edit(fn):
        mm = copy.deepcopy(m)
        fn(mm)
        return mm

    cases = [
        ("K0-unmodified", "control", {}, set()),
        ("K1-doc-crlf", "control", {"text": text.replace("\n", "\r\n")}, set()),
        ("K2-matrix-keys-sorted", "control", {"matrix": json.loads(json.dumps(m, sort_keys=True))}, set()),
        ("K3-cover-keys-sorted", "control", {"cover": json.loads(json.dumps(inputs["cover"], sort_keys=True))}, set()),
        ("M1-doc-winner", "mutation", {"text": text.replace("cell 3：$H_{19}$ 勝。", "cell 3：$B_7$ 勝。")}, {"A0"}),
        ("M2-cover-constraint", "mutation",
         {"cover": {**inputs["cover"], "strict_weighted_set_cover": {**inputs["cover"]["strict_weighted_set_cover"],
                                                                      "constraints": ["for every cell j: sum_K x_jK >= 0", "x_jK <= I[j,K]*y_K"]}}},
         {"W1"}),
        ("M3-marker-outer", "mutation", {"matrix": edit(lambda mm: mm["rows"][5].update(marker_outer=mm["rows"][5]["marker_outer"] + 1e-12))}, {"W2"}),
        ("M4-edge-pending", "mutation",
         {"matrix": edit(lambda mm: mm["rows"][3]["quantitative_edges"]["B7"].update(
             pending_leaves=mm["rows"][3]["quantitative_edges"]["B7"]["pending_leaves"] + 1))}, {"W3"}),
        ("M5-rebuild-claim", "mutation", {"matrix_rebuild": {**CLAIMS["matrix_rebuild"], "edges_exact": CLAIMS["matrix_rebuild"]["edges_exact"] + 1}},
         {"W3"}),
        ("M6-deep-B7-volume", "mutation",
         {"matrix": edit(lambda mm: mm["depth22_crossover_first4"][1].update(B7=mm["depth22_crossover_first4"][1]["B7"] + 1e-9))}, {"W4"}),
        ("M7-exact-claim", "mutation", {"exact_closes": [[c[0] + 1] + c[1:] if i == 7 else c for i, c in enumerate(CLAIMS["exact_closes"])]},
         {"W5"}),
        ("M8-limit-claim", "mutation", {"limits_not_below_T": [[0, 11]]}, {"W6"}),
        ("M9-limit-value", "mutation", {"limit_worst": {**CLAIMS["limit_worst"], "H19": CLAIMS["limit_worst"]["H19"] + 1e-9}}, {"W6"}),
        ("M10-ships-code", "mutation", {"recorded_facts": {**CLAIMS["recorded_facts"], "ships_code": True}}, {"W7"}),
        ("M11-margin", "mutation",
         {"matrix": edit(lambda mm: mm["aggregate"]["H19"].update(forcing_margin_search=5.96e-5))}, {"P1"}),
        ("M12-crystal-value", "mutation", {"crystal": inputs["crystal"].replace("H17 ~2.66e-4", "H17 ~2.67e-4")}, {"P1"}),
        ("M13-validation-digest", "mutation", {"validation": {**inputs["validation"], "portfolio_matrix_sha256": "0" * 64}}, {"P2"}),
        ("M14-validation-check", "mutation",
         {"validation": {**inputs["validation"], "checks": {**inputs["validation"]["checks"], "matrix_has_60_edges": False}}}, {"P2"}),
        ("M15-review-dropped", "mutation", {"reviews": list(R23_REVIEWS[1:])}, {"R1"}),
    ]
    if args.blind:
        cases = [c for c in cases if args.blind in c[3]]
    d = gate.run_drill(build, CHECKS, DEPENDS, cases)
    summary = gate.summarize(d, CHECKS)
    print(gate.report_drill(summary))
    for x in d:
        if not x["ok"]:
            print("   DRILL FAIL", x)
    if args.blind:
        n = sum(1 for x in d if args.blind in x["expected"] and args.blind not in x["tripped"])
        print(f"blind {args.blind}: drill reported it blind in {n} cases (nothing written)")
        return
    log = {"run": "RUN-026", "subject": "Round 23: witness portfolio closure, strict incidence and residual-frontier routing",
           "document": DOC, "checks": {c: {"pass": not f, "failures": f} for c, f in results.items()},
           **{k: w.get(k) for k in ("w1", "w2", "w3", "w4", "w5", "w6", "w7", "p1", "p2")}, "drill": summary}
    LOG.write_text(json.dumps(log, ensure_ascii=False, indent=1, default=float) + "\n", encoding="utf-8")
    if not results["R1"]:
        registry.record("R23", "RUN-026", R23_DEFINITIONS, R23_REVIEWS)
        print(f"registry: {len(R23_DEFINITIONS)} definitions, {len(R23_REVIEWS)} reviews recorded for R23")


if __name__ == "__main__":
    main()