#!/usr/bin/env python3
"""
The Fluidity Replication (plan.md, revision 1, October 14, 2026): a prospective test of two findings from Check 3 of the
pilot review, on fresh starts 115 to 134, with the pilot's model, interventions and measures unchanged.

  Hypothesis 1, the share gap: the slower half's share of the medium's extra coherence at R minus its share of the
  matched push's (Check 3's question 2, medium minus push).
  Hypothesis 2, the persistence gap: in the untouched runs, the slower half's coefficient minus the faster half's over
  the 100 units after each natural state (Check 3's question 3).

Each is judged on its own: PASS if the gap minus 3 bars is above zero and the gap is at least 0.10; FAIL if the gap plus
2 bars is below 0.10; otherwise INCONCLUSIVE, and INCONCLUSIVE whenever a validity condition fails.

medium_test.py (the pilot's code, with nudge_test.py's model and the cause test's helpers) and check3_fluidity.py are
imported unchanged, and the pilot's files are only read. The one addition is a read-only recording of each arm's exact
contributions (the slower half's, the faster half's and the free turners'), checked to change nothing.

A sandbox script kept with the project. Nothing here is served by the site.

Commands, each in the foreground with a hard limit (timeout -k 10 285 python3 -B replication.py <command>); nothing is
left running:
    checkruns   the model runs the checks need, in one piece: pilot start 103 (its base run, and moment 1 by both
                codes) and the stand-in on old start 1
    checks      the checks before the lock; writes checks.txt and checks.json (the report on made-up numbers goes to
                /tmp only)
    lock        the lock record; refused unless every check passed and no replication output exists
    go "..."    records Colson's go, in his words (used only when he gives it)
    run         one piece of the run on starts 115 to 134; refused without the lock and the go record
    report      the readings, once all 220 tasks are done
"""
import json
import shutil
import sys
import time
import warnings
from contextlib import contextmanager
from pathlib import Path

HERE = Path(__file__).resolve().parent
CODE = Path(__file__).resolve()
PILOT_DIR = HERE.parent
ROOT = HERE.parents[3]
sys.dont_write_bytecode = True                             # leave every folder as it is
sys.path.insert(0, str(PILOT_DIR))
import medium_test as mt  # noqa: E402  (first: it loads prelock_common, which sets one BLAS thread before numpy)
import check3_fluidity as c3  # noqa: E402
from medium_test import ct, nt, np, stamp, shown, sha256, write_json, plain, save_npz, run_piece, guard  # noqa: E402
from medium_test import lock_changes, machine  # noqa: E402
from prelock_common import same_environment  # noqa: E402

N, K, EPS, DT = mt.N, mt.K, mt.EPS, mt.DT
SEEDS = tuple(range(114, 134))                             # fresh starts 115 to 134
STARTS = tuple(s + 1 for s in SEEDS)
PILOT_STARTS = tuple(s + 1 for s in mt.PILOT_STARTS)       # 103 to 114: used only to check the code
CHECK_SEED = 102                                           # pilot start 103, reproduced to the last digit
MOMENTS, WINDOW, AFTER = mt.MOMENTS, mt.WINDOW, mt.AFTER
PLANNED = len(SEEDS) * MOMENTS                             # 200 planned moments
SESOI, CLEAR, SIDE = 0.10, 3.0, 2.0                        # smallest effect of interest; lower: gap - 3 bars; upper: + 2
MOST_MISSING = 20                                          # more than 1 in 10 of 200 not like for like: H1 inconclusive
TOL = 1e-12
PART_NAMES, BLOCKS, LATER = c3.PART_NAMES, c3.BLOCKS, c3.LATER
DIRECTIONS = mt.DIRECTIONS
NAN = float("nan")

WORK = ROOT / "generated-artifacts" / "medium-test" / "fluidity-replication"   # kept out of the history and the site
RUN_PARTS, CHECK_PARTS = WORK / "run", WORK / "checks"
CHECK_LOG = CHECK_PARTS / "checkruns_log.json"
PLAN, DRAFT = HERE / "plan.md", HERE / "plan-draft-1.md"
CHECKS_TXT, CHECKS_JSON = HERE / "checks.txt", HERE / "checks.json"
LOCK, GO = HERE / "replication_lock.json", HERE / "replication_go.json"
RUN_LOG = HERE / "replication_runlog.json"
RESULTS_TXT, RESULTS_JSON = HERE / "replication_results.txt", HERE / "replication_results.json"
DRY = Path("/tmp/fluidity-replication-dry")


def frozen_size():
    """The pilot's frozen size (0.5% of K), read from its freeze record."""
    return json.loads(mt.FREEZE.read_text())["chosen"]["size"]


SETTINGS = dict(
    n=N, k=K, eps=EPS, delta=nt.DELTA, dt=DT, settle=nt.SETTLE, run=mt.RUN, gap=mt.GAP, put_back_every=mt.EVERY,
    fast=nt.FAST, gap_seed=nt.LYAP_SEED, seed=mt.SEED, sham_draw=3, starts=list(STARTS), moments=MOMENTS,
    first=mt.FIRST, spacing=mt.SPACING, offsets=mt.OFFSETS, window=WINDOW, after=AFTER, size=frozen_size(),
    pulled=c3.PULLED, half=c3.HALF, scan=mt.SCAN, largest=mt.LARGEST, bisect=mt.BISECT, solved=mt.SOLVED,
    match=mt.MATCH, gap_kept=mt.KEPT, sham_band=list(mt.SHAM_BAND), resamples=mt.RESAMPLES,
    resample_seed=mt.RESAMPLE_SEED, check3_resamples=c3.RESAMPLES, check3_seed=c3.SEED, clear=CLEAR, side=SIDE,
    sesoi=SESOI, most_missing=MOST_MISSING, planned=PLANNED, blocks=[list(b) for b in BLOCKS], later=list(LATER),
    benchmark=mt.BENCHMARK, tolerance=TOL)
PLAN_NUMBERS = dict(                                       # the plan's numbers, which the settings must match
    n=4000, k=1.9, eps=0.1, delta=1.0, dt=0.025, settle=300.0, run=6000.0, gap=1e-7, put_back_every=5.0, fast=20.0,
    gap_seed=20261011, seed=20261014, sham_draw=3, starts=list(range(115, 135)), moments=10, first=2000, spacing=330,
    offsets=66, window=200, after=500, size=0.005, pulled=3872, half=1936, scan=256, largest=np.pi / 2, bisect=1e-12,
    solved=1e-6, match=1e-9, gap_kept=1e-6, sham_band=[0.95, 1.05], resamples=400, resample_seed=20261016,
    check3_resamples=400, check3_seed=20261016, clear=3.0, side=2.0, sesoi=0.10, most_missing=20, planned=200,
    blocks=[[0, 100], [100, 200], [200, 300], [300, 400], [400, 500]], later=[200, 500], benchmark=0.1602,
    tolerance=1e-12)
VERDICT_RULE = dict(
    threshold=SESOI, lower="gap - 3 bars", upper="gap + 2 bars", PASS="lower > 0 and gap >= 0.10",
    FAIL="upper < 0.10", INCONCLUSIVE="neither, or a validity condition fails",
    fail_kinds=dict(reversed="gap + 3 bars < 0", absent="lower <= 0 and gap + 3 bars >= 0",
                    real_but_smaller="lower > 0 and upper < 0.10"),
    bars=dict(share_gap="SD of the gap over 400 resamples of the starts (seed 20261016), both shares from the same "
                        "draws", persistence_gap="max(SD of the average over the same resamples, SD / sqrt(starts))"))


# ---------------------------------------------------------------------------------------------------------------------
# small helpers

def arrays_equal(a, b):
    """The same shape, type and bytes: identical to the last bit."""
    a, b = np.asarray(a), np.asarray(b)
    return bool(a.shape == b.shape and a.dtype == b.dtype and a.tobytes() == b.tobytes())


def est(x, digits=3):
    m, b = (x or {}).get("mean"), (x or {}).get("bar")
    if m is None or not np.isfinite(m):
        return "n/a"
    if b is None or not np.isfinite(b) or b == 0:
        return f"{m:+.{digits}f}"
    return f"{m:+.{digits}f} +- {b:.{digits}f} ({m / b:+.1f} bars)"


def outputs():
    """Anything the run on the fresh starts writes (none may exist before the lock)."""
    found = [p for p in sorted(RUN_PARTS.rglob("*")) if p.is_file()] if RUN_PARTS.exists() else []
    return [mt.shown_path(p) for p in found + [q for q in (RUN_LOG, RESULTS_TXT, RESULTS_JSON) if q.exists()]]


def lock_files():
    return dict(plan=PLAN, draft_1=DRAFT, runner=CODE, medium_test=Path(mt.__file__).resolve(),
                check3_fluidity=Path(c3.__file__).resolve(), nudge_test=Path(nt.__file__).resolve(),
                cause_test=Path(ct.__file__).resolve(),
                prelock_common=Path(sys.modules["prelock_common"].__file__).resolve(), pilot_freeze=mt.FREEZE,
                checks_txt=CHECKS_TXT, checks_json=CHECKS_JSON)


def pilot_records():
    """Every file of the pilot's records in medium-test/, outside this folder: they must stay exactly as they are."""
    return {str(p.relative_to(ROOT)): sha256(p) for p in sorted(PILOT_DIR.rglob("*"))
            if p.is_file() and HERE not in p.parents and "__pycache__" not in p.parts}


def pilot_data_problems():
    """The pilot's saved data against their recorded fingerprints (base runs in pilot_runlog.json, moments in
    pilot_results.json)."""
    runlog = json.loads(mt.PILOT_LOG.read_text())["parts"]
    moments = json.loads(mt.PILOT_JSON.read_text())["moment_sha256"]
    bad = []
    for s in PILOT_STARTS:
        p = mt.PILOT_PARTS / f"base_{s:03d}.npz"
        if runlog.get(p.stem, {}).get("sha256") != sha256(p):
            bad.append(p.name)
        for m in range(1, MOMENTS + 1):
            p = mt.PILOT_PARTS / f"moment_{s:03d}_{m:02d}.npz"
            if moments.get(p.name) != sha256(p):
                bad.append(p.name)
    return bad


def pilot_freeze_problems():
    """The pilot's freeze record verified: its protected files, its plan, code and calibration results; and the cause
    test's lock."""
    record = json.loads(mt.FREEZE.read_text())
    out = []
    if mt.protected_fingerprints() != record["protected"]:
        out.append("a file protected by the pilot's freeze changed")
    now = mt.current_fingerprints()
    out += [f"the pilot's {name} changed" for name in now if record["sha256"].get(name) != now[name]]
    if lock_changes():
        out.append("a file fixed by the cause test's lock changed")
    return out


# ---------------------------------------------------------------------------------------------------------------------
# the added recording (reads the state only)

def contributions(u, f, z):
    """Each group's coherence split exactly along its own rhythm into the slower half's, the faster half's and the
    free turners' parts (Check 3's measure; they add up to |z_g|): (group, 3). u (group, pulled) and f (group, free):
    one copy's state; z (group,): its coherence as the run measures it."""
    inner = mt.inner_half()
    turn = np.exp(-1j * np.angle(z))
    return np.stack([(u[:, inner].sum(-1) * turn).real / N, (u[:, ~inner].sum(-1) * turn).real / N,
                     (f.sum(-1) * turn).real / N], axis=-1)


def check3_parts(u, f):
    """Check 3's own calculation of one copy's contributions (u (group, pulled), f (group, free)): (group, 3)."""
    inner = c3.inner_half()
    return np.array([[c3.numbers(u[g], f[g], inner)[key] for key in ("c_in", "c_out", "c_free")] for g in range(2)])


def recording_branch(us, uf, schedule, copies, change=None):
    """medium_test.branch, line for line, with one addition: at every sample the first copy's contributions are
    recorded (contributions(), which only reads the state). Returns branch()'s outputs and "parts" (group, 3,
    samples)."""
    U = mt.shaped(copies, "u", us[:copies] if change is None else us[:copies] * np.exp(1j * change))
    F = mt.shaped(copies, "f", uf[:copies])
    reset_every, sample_every = (int(round(x / DT)) for x in (mt.EVERY, 1.0))
    samples, logs, pulls, ends = [mt.coherence(U[:1], F[:1])[0]], [], [], []
    parts = [contributions(U[0], F[0], samples[0])]
    none = np.zeros(copies)
    for k, units in schedule:
        step = nt.model(N, k, EPS)[1](DT)
        for i in range(1, int(round(units / DT)) + 1):
            step(U, F, none)
            if i % sample_every == 0:
                samples.append(mt.coherence(U[:1], F[:1])[0])
                parts.append(contributions(U[0], F[0], samples[-1]))
                pulls.append(k)
            if copies == 2 and i % reset_every == 0:
                shift = np.angle(U[1:] * np.conj(U[:1]))
                shift -= shift.mean(axis=(-2, -1), keepdims=True)
                size = np.sqrt((shift * shift).mean(axis=(-2, -1)))
                logs.append(np.log(size / mt.GAP)[0])
                U[1:] = U[:1] * np.exp(1j * shift * (mt.GAP / size)[:, None, None])
        ends.append((U.copy(), F.copy()))
    return dict(z=np.stack(samples, -1), logs=np.array(logs), pulls=np.array(pulls), ends=ends,
                parts=np.stack(parts, -1))


def replication_moment(task, parts, window=WINDOW, after=AFTER):
    """medium_test.pilot_moment, line for line, with the arms run by recording_branch: every array the pilot saves,
    plus each arm's contributions at every time unit (parts_medium (direction, group, part, 701 from T); parts_push and
    parts_sham (direction, group, part, 501 from R), not-a-number where not built), the untouched state's at T and R,
    and two checks of the recording: against Check 3's own calculation from the states at R, and adding up."""
    s, m, f = task["seed"], task["moment"], task["size"]
    T = task.get("T", mt.moment_times(s)[m])
    R = T + window
    with np.load(parts / f"base_{s + 1:03d}.npz") as base:
        us_T, uf_T = mt.kept_state(base, T)
        us_R, uf_R = mt.kept_state(base, R)
        z_none = base["z"][:, T:R + after + 1]
        logs_none = base["logs"][T // 5:(R + after) // 5]
    medium = [recording_branch(us_T, uf_T, [(K * (1 + sign * f), window), (K, after)], 2) for sign in (1, -1)]
    records, changes, pushes, shams = [], [], [], []
    for d in range(2):
        rec, dtheta = mt.matched(us_R, uf_R, np.abs(medium[d]["z"][:, window]), (3, s, m, d))
        records.append(rec)
        changes.append(dtheta)
        pushes.append(recording_branch(us_R, uf_R, [(K, after)], 2, change=dtheta[0]) if rec["built"] else None)
        shams.append(recording_branch(us_R, uf_R, [(K, after)], 2, change=dtheta[1]) if rec["built"] else None)
    states = [(us_R[0], uf_R[0])] + [(arm["ends"][0][0][0], arm["ends"][0][1][0]) for arm in medium]
    for b in range(2):                                     # the pushes', then the shams', states at R
        states += [(us_R[0] * np.exp(1j * changes[d][b]), uf_R[0]) if records[d]["built"] else None for d in range(2)]
    desc = {key: np.full((len(mt.STATES), 2), np.nan) for key in ("R", "inner", "outer", "z2", "z3", "z4")}
    desc.update(C=np.full(len(mt.STATES), np.nan), ab=np.full(len(mt.STATES), np.nan))
    for i, state in enumerate(states):
        if state is not None:
            numbers = mt.describe(*state)
            for key in mt.DESC_KEYS:
                desc[key][i] = numbers[key]
    apart = [mt.distance(states[1 + d][0], states[3 + d][0]) if records[d]["built"] else np.nan for d in range(2)]
    delivered = np.array([f * K * np.abs(arm["z"][:, 1:window + 1]).sum(-1) for arm in medium])
    restored = np.array([bool(len(arm["pulls"]) == window + after and np.all(arm["pulls"][window:] == K)
                              and np.all(arm["pulls"][:window] == K * (1 + sign * f)))
                         for arm, sign in zip(medium, (1, -1))])
    blank_z, blank_logs = np.full((2, after + 1), np.nan + 0j), np.full(after // 5, np.nan)
    blank_parts = np.full((2, 3, after + 1), np.nan)
    arrays = {f"rec_{key}": np.array([r[key] for r in records]) for key in mt.REC_KEYS}
    # the recording, checked: against Check 3's calculation from the states at R, and adding up at every time
    pairs = [(medium[d]["parts"][..., window], states[1 + d]) for d in range(2)]
    pairs += [(pushes[d]["parts"][..., 0], states[3 + d]) for d in range(2) if pushes[d]]
    pairs += [(shams[d]["parts"][..., 0], states[5 + d]) for d in range(2) if shams[d]]
    parts_check = max(float(np.max(np.abs(recorded - check3_parts(*state)))) for recorded, state in pairs)
    arms = medium + [p for p in pushes + shams if p]
    adds_up = max(float(np.max(np.abs(arm["parts"].sum(-2) - np.abs(arm["z"])))) for arm in arms)
    part = parts / f"{task['name']}.npz"
    save_npz(part, start=s + 1, moment=m + 1, T=T, R=R, size=f, z_none=z_none, logs_none=logs_none,
             z_medium=np.stack([arm["z"] for arm in medium]), logs_medium=np.stack([arm["logs"] for arm in medium]),
             z_push=np.stack([p["z"] if p else blank_z for p in pushes]),
             logs_push=np.stack([p["logs"] if p else blank_logs for p in pushes]),
             z_sham=np.stack([p["z"] if p else blank_z for p in shams]),
             logs_sham=np.stack([p["logs"] if p else blank_logs for p in shams]),
             built=np.array([r["built"] for r in records]), reason=np.array([r["reason"] for r in records]),
             distance=np.array(apart), delivered=delivered, restored=restored,
             **arrays, **{f"desc_{key}": value for key, value in desc.items()},
             parts_medium=np.stack([arm["parts"] for arm in medium]),
             parts_push=np.stack([p["parts"] if p else blank_parts for p in pushes]),
             parts_sham=np.stack([p["parts"] if p else blank_parts for p in shams]),
             parts_none_T=check3_parts(us_T[0], uf_T[0]), parts_none_R=check3_parts(us_R[0], uf_R[0]),
             parts_check=parts_check, parts_adds_up=adds_up)
    medium_finite = all(np.isfinite(arm["z"]).all() and np.isfinite(arm["logs"]).all() for arm in medium)
    problems = [] if restored.all() and medium_finite else [f"restored {restored.tolist()}, finite {medium_finite}"]
    if not (parts_check <= TOL and adds_up <= TOL):
        problems.append("the recorded parts don't add up or don't match Check 3's calculation")
    return part, "; ".join(problems)


# ---------------------------------------------------------------------------------------------------------------------
# the stand-in on old start 1 (equality and mechanics only)

def same_arm(a, b):
    """Every output of the unchanged branch and of the recording one, to the last digit."""
    return bool(arrays_equal(a["z"], b["z"]) and arrays_equal(a["logs"], b["logs"])
                and arrays_equal(a["pulls"], b["pulls"]) and len(a["ends"]) == len(b["ends"])
                and all(arrays_equal(x[0], y[0]) and arrays_equal(x[1], y[1]) for x, y in zip(a["ends"], b["ends"])))


def largest_misfit(arm):
    return float(np.max(np.abs(arm["parts"].sum(-2) - np.abs(arm["z"]))))


def standin(task, parts):
    """Old start 1 (seed 0), short runs, as in the pilot's checks: a 400-unit base run; a zero change and a changed
    pull, each by the unchanged branch and the recording one; a push and its sham, each both ways; and the pilot's and
    the replication's moment tasks on a stand-in moment (T = 100, the pull changed 5% for 20 units, 50 after)."""
    s, run, T0, short, after = mt.CHECK_START, 400, 100, 20, 50
    base = mt.base_run(s, np.inf, parts / "never.state.npz", (T0, T0 + short), run)
    save_npz(parts / "base_001.npz", start=1, z=base["z"], logs=base["logs"], kept_t=base["kept_t"],
             kept_u=base["kept_u"], kept_f=base["kept_f"])
    kept = dict(kept_t=base["kept_t"], kept_u=base["kept_u"], kept_f=base["kept_f"])
    us, uf = mt.kept_state(kept, T0)
    out = {}
    zero_old, zero_new = (fn(us, uf, [(K, 50), (K, 50)], 2) for fn in (mt.branch, recording_branch))
    out["zero change"] = dict(
        reproduces_untouched=bool(arrays_equal(zero_new["z"], base["z"][:, T0:T0 + 101])
                                  and arrays_equal(zero_new["logs"], base["logs"][T0 // 5:T0 // 5 + 20])),
        recording_arm_identical=same_arm(zero_old, zero_new), parts_add_up=largest_misfit(zero_new) <= TOL)
    changed_old, changed_new = (fn(us, uf, [(K * 1.05, 50), (K, 50)], 2) for fn in (mt.branch, recording_branch))
    out["pull changed 5%, then restored"] = dict(
        recording_arm_identical=same_arm(changed_old, changed_new),
        pull_as_scheduled=bool(np.all(changed_new["pulls"][:50] == K * 1.05) and np.all(changed_new["pulls"][50:] == K)),
        change_took_effect=not arrays_equal(changed_new["z"][:, 1:51], base["z"][:, T0 + 1:T0 + 51]),
        parts_add_up=largest_misfit(changed_new) <= TOL)
    usR, ufR = mt.kept_state(kept, T0 + short)
    R0 = np.abs(mt.coherence(usR[:1], ufR[:1])[0])
    for name, factor, code in (("push raising", mt.MARKS["raise_low"], 0), ("push lowering", mt.MARKS["lower_high"], 1)):
        rec, dtheta = mt.matched(usR, ufR, factor * R0, (9, s, 0, code))
        mech = dict(built=bool(rec["built"]), on_target=bool(np.all(np.abs(rec["push_ratio"] - 1) <= mt.SOLVED)),
                    sham_in_band=bool(np.all((mt.SHAM_BAND[0] <= rec["sham_ratio"])
                                             & (rec["sham_ratio"] <= mt.SHAM_BAND[1]))),
                    same_sizes=bool(rec["same"] and rec["match"] <= mt.MATCH),
                    gap_kept=bool(np.all(rec["gaprel"] < mt.KEPT)))
        if rec["built"]:
            for b, arm in ((0, "push"), (1, "sham")):
                old = mt.branch(usR, ufR, [(K, 20)], 2, change=dtheta[b])
                new = recording_branch(usR, ufR, [(K, 20)], 2, change=dtheta[b])
                mech[f"{arm}_recording_arm_identical"] = same_arm(old, new)
                mech[f"{arm}_parts_add_up"] = largest_misfit(new) <= TOL
        out[name] = mech
    new_part, new_problem = replication_moment(dict(name="standin_moment_new", seed=s, moment=0, size=0.05, T=T0),
                                               parts, short, after)
    old_part, old_problem = mt.pilot_moment(dict(name="standin_moment_old", seed=s, moment=0, size=0.05, T=T0),
                                            parts, short, after)
    with np.load(new_part) as a, np.load(old_part) as b:
        out["moment task on the stand-in"] = dict(
            pilot_arrays_identical=all(key in a.files and arrays_equal(a[key], b[key]) for key in b.files),
            pilot_arrays=len(b.files), added=sorted(set(a.files) - set(b.files)),
            shapes=bool(a["parts_medium"].shape == (2, 2, 3, short + after + 1)
                        and a["parts_push"].shape == (2, 2, 3, after + 1)
                        and a["parts_sham"].shape == (2, 2, 3, after + 1)),
            restored=bool(a["restored"].all()), problem=new_problem or old_problem or "none",
            pushes_built=[bool(x) for x in a["built"]], parts_match_check3=float(a["parts_check"]) <= TOL,
            parts_add_up=float(a["parts_adds_up"]) <= TOL)
    path = parts / "standin.json"
    write_json(path, out)
    return path, ""


# ---------------------------------------------------------------------------------------------------------------------
# the tasks

def work(task, deadline):
    parts = Path(task["parts"])
    parts.mkdir(parents=True, exist_ok=True)
    if task["kind"] == "base":
        return mt.base_task(task, deadline, parts)
    run = dict(pilot=mt.pilot_moment, replication=replication_moment, standin=standin)[task["kind"]]
    part, problem = run(task, parts)
    return dict(status="done", sha256=sha256(part), error=problem)


def job(args):
    return guard(work, args)


def checkrun_tasks(size):
    return [dict(name="base_103", kind="base", seed=CHECK_SEED, keep=mt.keep_times(CHECK_SEED), logs=True,
                 parts=str(CHECK_PARTS / "rerun")),
            dict(name="pilot_code_moment_103_01", kind="pilot", seed=CHECK_SEED, moment=0, size=size,
                 parts=str(CHECK_PARTS / "pilot-code")),
            dict(name="new_code_moment_103_01", kind="replication", seed=CHECK_SEED, moment=0, size=size,
                 parts=str(CHECK_PARTS / "new-code")),
            dict(name="standin", kind="standin", parts=str(CHECK_PARTS / "standin"))]


def run_tasks(size):
    base = [dict(name=f"base_{s + 1:03d}", kind="base", seed=s, keep=mt.keep_times(s), logs=True, parts=str(RUN_PARTS))
            for s in SEEDS]
    moments = [dict(name=f"moment_{s + 1:03d}_{m + 1:02d}", kind="replication", seed=s, moment=m, size=size,
                    parts=str(RUN_PARTS)) for s in SEEDS for m in range(MOMENTS)]
    return base, moments


# ---------------------------------------------------------------------------------------------------------------------
# Check 3's calculations for any list of starts

def share_n(a, b, n):
    """check3_fluidity.share for n starts; None (undefined) if a start has no value."""
    a, b = np.asarray(a, float), np.asarray(b, float)
    if not (len(a) == n and np.isfinite(a).all() and np.isfinite(b).all()):
        return None
    idx = np.random.default_rng(c3.SEED).integers(0, len(a), (c3.RESAMPLES, len(a)))
    return float(a.mean() / (a.mean() + b.mean())), a[idx].mean(-1) / (a[idx].mean(-1) + b[idx].mean(-1))


def as_dict_n(s):
    return dict(mean=NAN, bar=NAN) if s is None else c3.as_dict(s)


def compared_n(x, y):
    return dict(mean=NAN, bar=NAN, clear=False) if x is None or y is None else c3.compared(x, y)


def ratio_n(numerator, denominator):
    a, b = np.asarray(numerator, float), np.asarray(denominator, float)
    if int((np.isfinite(a) & np.isfinite(b)).sum()) < 2:
        return dict(mean=NAN, bar=NAN, n=int((np.isfinite(a) & np.isfinite(b)).sum()), clear=False)
    return c3.ratio_summary(a, b)


def fluidity_numbers(starts, base_file, moment_file, base_hash=None, moment_hash=None, report=None):
    """Check 3's calculations (check3_fluidity.main), step for step and with its own functions, for any list of
    starts. The differences: the starts and file locations are arguments; fingerprints are checked when given; the two
    checks against the pilot's report run only when it is given; and a share with a start missing is left undefined
    rather than stopping. Pointed at the pilot's files, it must give check3_fluidity.json's numbers exactly."""
    inner = c3.inner_half()
    if inner.size != c3.PULLED or int(inner.sum()) != c3.HALF:
        sys.exit("The slower half isn't 1,936 of 3,872 pulled oscillators. Stopping.")
    S_, M_ = len(starts), len(c3.MOMENTS)
    nat = {key: np.full((S_, 2, 2 * M_), np.nan) for key in ("R", "c_in", "c_out", "c_free", "m_in", "m_out")}
    later = np.full((S_, 2, 2 * M_, len(BLOCKS)), np.nan)
    late = np.full((S_, 2, 2 * M_), np.nan)
    at_R = np.zeros((S_, 2 * M_), bool)
    rows = []
    worst = dict(desc=0.0, base_z=0.0, z_none=0.0, parts=0.0)
    for i, s in enumerate(starts):
        path = base_file(s)
        if base_hash is not None and base_hash(path) != sha256(path):
            sys.exit(f"{path.name} differs from its recorded fingerprint. Stopping.")
        with np.load(path) as b:
            z, kept_t, kept_u, kept_f = b["z"], b["kept_t"], b["kept_u"], b["kept_f"]
        moment_rows = []
        for m in c3.MOMENTS:
            path = moment_file(s, m)
            if moment_hash is not None and moment_hash(path) != sha256(path):
                sys.exit(f"{path.name} differs from its recorded fingerprint. Stopping.")
            with np.load(path) as p:
                moment_rows.append({key: p[key] for key in ("T", "R", "z_none", "z_medium", "z_push", "built", "desc_R",
                                                            "desc_inner", "desc_outer")})
        times = [int(t) for t in kept_t]
        if times != sorted({int(r["T"]) for r in moment_rows} | {int(r["R"]) for r in moment_rows}):
            sys.exit(f"base_{s:03d}'s saved states aren't the moments' T and R. Stopping.")
        for j, t in enumerate(times):
            at_R[i, j] = t in {int(r["R"]) for r in moment_rows}
            for g in range(2):
                values = c3.numbers(kept_u[j, 0, g], kept_f[j, 0, g], inner)
                for key in nat:
                    nat[key][i, g, j] = values[key]
                worst["base_z"] = max(worst["base_z"], abs(values["R"] - abs(z[g, t])))
                worst["parts"] = max(worst["parts"], abs(values["c_in"] + values["c_out"] + values["c_free"]
                                                         - values["R"]))
                coherence = np.abs(z[g])
                for k, (lo, hi) in enumerate(BLOCKS):
                    later[i, g, j, k] = coherence[t + lo + 1:t + hi + 1].mean()
                late[i, g, j] = coherence[t + LATER[0] + 1:t + LATER[1] + 1].mean()
        for r in moment_rows:
            j, T = times.index(int(r["R"])), int(r["T"])
            for g in range(2):
                worst["desc"] = max(worst["desc"], abs(nat["R"][i, g, j] - r["desc_R"][0, g]),
                                    abs(nat["m_in"][i, g, j] - r["desc_inner"][0, g]),
                                    abs(nat["m_out"][i, g, j] - r["desc_outer"][0, g]))
            worst["z_none"] = max(worst["z_none"], float(np.max(np.abs(r["z_none"] - z[:, T:T + c3.W + c3.AFTER + 1]))))
        rows += moment_rows

    def grid(fn):
        return np.array([[fn(rows[i * M_ + m]) for m in range(M_)] for i in range(S_)], float)

    def C(zz):
        return np.abs(zz).mean(-2)

    built = np.array([[rows[i * M_ + m]["built"] for m in range(M_)] for i in range(S_)], bool)
    mask = np.where(built.all(-1), 1.0, np.nan)
    checks = {
        "untouched states at R from the base runs against the report's descriptions (coherence, slower half, faster "
        "half)": worst["desc"],
        "coherence from the saved states against the base runs' record at T and R": worst["base_z"],
        "the moment files' untouched coherence against the base runs": worst["z_none"],
        "the three contributions adding up to the group's coherence": worst["parts"]}
    if report is not None:
        worst_state = 0.0
        for i_state, state in enumerate(c3.STATES):
            m_ = mask if i_state >= 3 else 1.0
            for key in ("R", "inner", "outer"):
                mine = float(np.mean(c3.by_start(grid(lambda r, i_state=i_state, key=key:
                                                      np.mean(np.abs(r[f"desc_{key}"][i_state]))) * m_)))
                worst_state = max(worst_state, abs(mine - report["descriptions"]["states"][state][key]))
        checks["the report's state descriptions at R (coherence, slower half, faster half; seven states)"] = worst_state
    diffR = grid(lambda r: C(r["z_medium"][0])[c3.W] - C(r["z_medium"][1])[c3.W])

    def medium_block(lo, hi):
        return grid(lambda r: (C(r["z_medium"][0]) - C(r["z_medium"][1]))[c3.W + lo + 1:c3.W + hi + 1].mean())

    def push_block(lo, hi):
        return grid(lambda r: (C(r["z_push"][0]) - C(r["z_push"][1]))[lo + 1:hi + 1].mean())

    persistence = dict(medium=ratio_n(c3.by_start(medium_block(*LATER) * mask), c3.by_start(diffR * mask)),
                       push=ratio_n(c3.by_start(push_block(*LATER) * mask), c3.by_start(diffR * mask)))
    if report is not None:
        checks["the report's persistence over (R + 200, R + 500], medium and push"] = max(
            abs(persistence[k][q] - report["persistence"][k][q]) for k in ("medium", "push") for q in ("mean", "bar"))

    # question 1: where natural order sits (exact contributions)
    keys = ("c_in", "c_out", "c_free")

    def shares_of(select):
        return np.array([[c3.within_slope(nat[k][i][:, select(i)], nat["R"][i][:, select(i)]) for k in keys]
                         for i in range(S_)])

    every, only_R, only_T = (lambda i: slice(None)), (lambda i: at_R[i]), (lambda i: ~at_R[i])
    a = shares_of(every)
    q1 = c3.summary(a[:, 0] - c3.EVEN)
    verdict1 = ("more in the slower half than an even spread" if q1["clear"] and q1["mean"] > 0
                else "less in the slower half than an even spread" if q1["clear"] else "not clear")
    one = dict(
        shares={name: c3.summary(a[:, k]) for k, name in enumerate(PART_NAMES)}, slower_minus_even=q1,
        sum_of_shares_largest_departure_from_1=float(np.max(np.abs(a.sum(-1) - 1))),
        at_R_alone={name: c3.summary(v) for name, v in zip(PART_NAMES, shares_of(only_R).T)},
        at_T_alone={name: c3.summary(v) for name, v in zip(PART_NAMES, shares_of(only_T).T)},
        made_up_of={name: ratio_n(nat[k].reshape(S_, -1).mean(-1), nat["R"].reshape(S_, -1).mean(-1))
                    for name, k in zip(PART_NAMES, keys)},
        average_R=c3.summary(nat["R"].reshape(S_, -1).mean(-1)),
        spread_of_R=c3.summary([float(np.sqrt(np.mean((nat["R"][i] - nat["R"][i].mean(-1, keepdims=True)) ** 2)))
                                for i in range(S_)]),
        random_floor=dict(half=float(np.sqrt(np.pi / (4 * c3.HALF))), group=float(np.sqrt(np.pi / (4 * c3.N)))))

    # question 2: the same relationship across routes (each half's own coherence)
    def slopes_of(select):
        return np.array([[c3.within_slope(nat[k][i][:, select(i)], nat["R"][i][:, select(i)]) for k in ("m_in", "m_out")]
                         for i in range(S_)])

    b = slopes_of(every)

    def half_split(key, x, y):
        return grid(lambda r: float(np.mean(r[f"desc_{key}"][x] - r[f"desc_{key}"][y])) / 2) * mask

    d = {route: {key: c3.by_start(half_split(key, x, y)) for key in ("inner", "outer", "R")}
         for route, (x, y) in (("medium", (1, 2)), ("push", (3, 4)))}
    sh = dict(natural=share_n(b[:, 0], b[:, 1], S_), medium=share_n(d["medium"]["inner"], d["medium"]["outer"], S_),
              push=share_n(d["push"]["inner"], d["push"]["outer"], S_))
    nm, np_, mp = compared_n(sh["natural"], sh["medium"]), compared_n(sh["natural"], sh["push"]), \
        compared_n(sh["medium"], sh["push"])
    if np_["clear"] and not nm["clear"]:
        verdict2 = "like the medium"
    elif nm["clear"] and not np_["clear"]:
        verdict2 = "like the push"
    elif nm["clear"] and np_["clear"]:
        low, high = sorted((sh["medium"][0], sh["push"][0]))
        verdict2 = ("between the two" if low < sh["natural"][0] < high else
                    "outside the two (above both)" if sh["natural"][0] > high else "outside the two (below both)")
    else:
        verdict2 = "not clear"

    def against_untouched(key, arm, lowered):
        if lowered:
            return grid(lambda r: float(np.mean(r[f"desc_{key}"][0] - r[f"desc_{key}"][arm]))) * mask
        return grid(lambda r: float(np.mean(r[f"desc_{key}"][arm] - r[f"desc_{key}"][0]))) * mask

    separately, separate_shares = {}, {}
    for name, arm, lowered in (("medium raised", 1, False), ("medium lowered", 2, True), ("push raised", 3, False),
                               ("push lowered", 4, True)):
        x, y = c3.by_start(against_untouched("inner", arm, lowered)), c3.by_start(against_untouched("outer", arm, lowered))
        separate_shares[name] = share_n(x, y, S_)
        separately[name] = dict(share=as_dict_n(separate_shares[name]), slower=c3.summary(x), faster=c3.summary(y))
    all_moments = share_n(c3.by_start(grid(lambda r: float(np.mean(r["desc_inner"][1] - r["desc_inner"][2])) / 2)),
                          c3.by_start(grid(lambda r: float(np.mean(r["desc_outer"][1] - r["desc_outer"][2])) / 2)), S_)
    natural_R, natural_T = share_n(*slopes_of(only_R).T, S_), share_n(*slopes_of(only_T).T, S_)
    two = dict(
        share={k: as_dict_n(v) for k, v in sh.items()},
        natural_minus_medium=nm, natural_minus_push=np_, medium_minus_push=mp,
        per_unit=dict(natural=dict(slower=c3.summary(b[:, 0]), faster=c3.summary(b[:, 1])),
                      **{route: dict(slower=ratio_n(d[route]["inner"], d[route]["R"]),
                                     faster=ratio_n(d[route]["outer"], d[route]["R"]),
                                     group_change=c3.summary(d[route]["R"])) for route in ("medium", "push")}),
        natural_at_R_alone=as_dict_n(natural_R), natural_at_T_alone=as_dict_n(natural_T),
        separately=separately, medium_all_moments=as_dict_n(all_moments))

    # question 3: which order lasts (exact contributions at t, the group's coherence later)
    beta = np.array([[c3.within_joint(later[i, :, :, k], [nat[key][i] for key in keys]) for k in range(len(BLOCKS))]
                     for i in range(S_)])                    # (start, block, part)
    q3 = c3.summary(beta[:, 0, 0] - beta[:, 0, 1])
    verdict3 = ("the slower half's order lasts longer" if q3["clear"] and q3["mean"] > 0
                else "the faster half's order lasts longer" if q3["clear"] else "not clear")
    by_block = [dict(block=f"(t+{lo}, t+{hi}]", **{name: c3.summary(beta[:, k, p]) for p, name in enumerate(PART_NAMES)},
                     slower_minus_faster=c3.summary(beta[:, k, 0] - beta[:, k, 1])) for k, (lo, hi) in enumerate(BLOCKS)]
    still_clear = [row["block"] for row in by_block if row["slower half"]["clear"] and row["slower half"]["mean"] > 0]
    survival = [dict(block=f"({lo}, {hi}]",
                     natural=c3.summary([c3.within_slope(later[i, :, :, k], nat["R"][i]) for i in range(S_)]),
                     medium=ratio_n(c3.by_start(medium_block(lo, hi) * mask), c3.by_start(diffR * mask)),
                     push=ratio_n(c3.by_start(push_block(lo, hi) * mask), c3.by_start(diffR * mask)))
                for k, (lo, hi) in enumerate(BLOCKS)]
    survival_late = dict(natural=c3.summary([c3.within_slope(late[i], nat["R"][i]) for i in range(S_)]), **persistence)
    three = dict(main=q3, by_block=by_block, slower_half_still_clear=still_clear, survival=survival,
                 survival_over_200_500=survival_late)
    result = dict(verdicts=dict(where=verdict1, routes=verdict2, lasts=verdict3), reproduced=checks,
                  question_1=one, question_2=two, question_3=three)
    return result, dict(built=built, mask=mask, d=d, sh=sh, beta=beta, all_moments=all_moments, nat=nat)


def as_json(x):
    """As check3_fluidity wrote its record (numpy values as lists, not-a-number kept), read back."""
    return json.loads(json.dumps(x, default=lambda o: o.tolist() if hasattr(o, "tolist") else str(o)))


def number_misfit(a, b):
    """How many numbers two records share, and their largest difference (inf if their structure differs)."""
    if isinstance(a, dict) and isinstance(b, dict):
        if set(a) != set(b):
            return 0, float("inf")
        pairs = [number_misfit(a[k], b[k]) for k in a]
    elif isinstance(a, list) and isinstance(b, list):
        if len(a) != len(b):
            return 0, float("inf")
        pairs = [number_misfit(x, y) for x, y in zip(a, b)]
    elif isinstance(a, bool) or isinstance(b, bool) or isinstance(a, str) or isinstance(b, str) or a is None or b is None:
        return 0, (0.0 if a == b else float("inf"))
    else:
        both_nan = a != a and b != b
        return 1, (0.0 if both_nan or a == b else abs(a - b))
    return sum(p[0] for p in pairs), max([p[1] for p in pairs], default=0.0)


# ---------------------------------------------------------------------------------------------------------------------
# the secondary readings, by the pilot's own analysis function

@contextmanager
def starts_for_analysis(seeds, most_missing):
    """medium_test.analyse reads its list of starts and its 1-in-10 cap from two module settings. They are set for the
    call and put back afterwards; the function itself is unchanged."""
    saved = mt.PILOT_STARTS, mt.MOST_MISSING
    mt.PILOT_STARTS, mt.MOST_MISSING = tuple(seeds), most_missing
    try:
        yield
    finally:
        mt.PILOT_STARTS, mt.MOST_MISSING = saved


def secondary_readings(rows, size, seeds=SEEDS, most_missing=MOST_MISSING):
    with starts_for_analysis(seeds, most_missing):
        return mt.analyse(rows, size)


SECONDARY_ITEMS = (
    ("the medium's S (later growth against untouched)", "later", "S_medium"),
    ("the medium's P", "later", "P_medium"),
    ("D_raised", "later", "D_raised"),
    ("D_lowered", "later", "D_lowered"),
    ("the pushes against their shams: S", "later", "S_push_sham"),
    ("the pushes against their shams: P", "later", "P_push_sham"),
    ("the pushes against untouched: S", "later", "S_push"),
    ("the pushes against untouched: P", "later", "P_push"),
    ("the shams against untouched, raised", "later", "sham_raised"),
    ("the shams against untouched, lowered", "later", "sham_lowered"),
    ("medium minus push, like for like: S difference", "later", "S_difference"),
    ("medium minus push, like for like: P difference", "later", "P_difference"),
    ("persistence over (R + 200, R + 500]: medium", "persistence", "medium"),
    ("persistence: push", "persistence", "push"),
    ("persistence: difference (medium minus push)", "persistence", "difference"),
    ("lasting: the medium's last block", "lasting", "medium_last"),
    ("lasting: the medium's decline", "lasting", "medium_decline"),
    ("lasting: the push's decline", "lasting", "push_decline"))


def beside(new, old):
    m, p = new.get("mean"), (old or {}).get("mean")
    if m is None or p is None or not (np.isfinite(m) and np.isfinite(p)):
        return "n/a"
    same = 'same' if np.sign(m) == np.sign(p) else 'opposite'
    return f"{same} sign and clear" if new.get("clear") else f"{same} sign, not clear"


# ---------------------------------------------------------------------------------------------------------------------
# the verdicts

def verdict(gap, bar, problems=()):
    """The plan's rule for one primary gap. PASS: gap - 3 bars > 0 and gap >= 0.10. FAIL: gap + 2 bars < 0.10, with
    its kind. INCONCLUSIVE: neither, no estimate, or a validity condition failed (never negative)."""
    gap, bar = float(gap), float(bar)
    out = dict(gap=gap, bar=bar, lower=gap - CLEAR * bar, upper=gap + SIDE * bar, threshold=SESOI,
               problems=list(problems))
    if problems:
        out.update(verdict="INCONCLUSIVE", kind="a validity condition failed: " + "; ".join(problems))
    elif not (np.isfinite(gap) and np.isfinite(bar)):
        out.update(verdict="INCONCLUSIVE", kind="no estimate")
    elif out["lower"] > 0 and gap >= SESOI:
        out.update(verdict="PASS", kind="clearly above zero and at least 0.10")
    elif out["upper"] < SESOI:
        out.update(verdict="FAIL", kind="reversed" if gap + CLEAR * bar < 0 else
                   "real but smaller than 0.10" if out["lower"] > 0 else "absent")
    else:
        out.update(verdict="INCONCLUSIVE", kind="the bounds can't tell it from 0.10")
    return out


def h1_problems(missing, starts_without, extras, exact, failed=()):
    """Hypothesis 1's validity conditions (plan.md); any one makes it INCONCLUSIVE."""
    out = []
    if missing > MOST_MISSING:
        out.append(f"{missing} of the {PLANNED} planned moments aren't like for like (at most {MOST_MISSING})")
    if starts_without:
        out.append("no like-for-like moment in start " + ", ".join(str(s) for s in starts_without))
    for route, value in extras.items():
        if not (np.isfinite(value) and value > 0):
            out.append(f"the {route}'s extra isn't above zero")
    if np.isfinite(exact["mean"]) and np.isfinite(exact["bar"]) and exact["mean"] + CLEAR * exact["bar"] < 0:
        out.append("the exact-contribution version is clearly opposite in sign")
    return out + list(failed)


VERDICT_CASES = (  # (gap, bar, problems, verdict, kind): the plan's examples, then each validity condition
    (0.25, 0.03, (), "PASS", None), (0.12, 0.03, (), "PASS", None), (0.08, 0.03, (), "INCONCLUSIVE", None),
    (0.06, 0.03, (), "INCONCLUSIVE", None), (0.02, 0.03, (), "FAIL", "absent"), (-0.15, 0.03, (), "FAIL", "reversed"),
    (0.06, 0.01, (), "FAIL", "real but smaller than 0.10"), (0.10, 0.03, (), "PASS", None),
    (NAN, NAN, (), "INCONCLUSIVE", "no estimate"))
CONDITION_CASES = (  # (missing, starts without, extras, exact gap, failed checks): each must give INCONCLUSIVE
    (21, [], dict(medium=0.01, push=0.01), dict(mean=0.2, bar=0.02), ()),
    (3, [117], dict(medium=0.01, push=0.01), dict(mean=0.2, bar=0.02), ()),
    (3, [], dict(medium=0.01, push=-0.001), dict(mean=0.2, bar=0.02), ()),
    (3, [], dict(medium=0.01, push=NAN), dict(mean=0.2, bar=0.02), ()),
    (3, [], dict(medium=0.01, push=0.01), dict(mean=-0.2, bar=0.02), ()),
    (3, [], dict(medium=0.01, push=0.01), dict(mean=0.2, bar=0.02), ("a failed check",)))


def verdict_problems():
    """The verdict code against the plan's examples and each validity condition: the cases it gets wrong."""
    wrong = []
    for gap, bar, problems, want, kind in VERDICT_CASES:
        got = verdict(gap, bar, problems)
        if got["verdict"] != want or (kind is not None and got["kind"] != kind):
            wrong.append(f"gap {gap}, bar {bar}: {got['verdict']} ({got['kind']}), planned {want}")
    if h1_problems(20, [], dict(medium=0.01, push=0.01), dict(mean=0.2, bar=0.02)):
        wrong.append("20 moments not like for like should still be valid")
    for case in CONDITION_CASES:
        problems = h1_problems(*case)
        if not problems or verdict(0.25, 0.03, problems)["verdict"] != "INCONCLUSIVE":
            wrong.append(f"condition {case} didn't give INCONCLUSIVE")
    return wrong


# ---------------------------------------------------------------------------------------------------------------------
# the exploratory reading: does slower-half concentration explain the route difference in persistence?

def composition(rows, mask):
    """plan.md, "Exploratory". Each route's extra (half of raised minus lowered) split exactly into the three parts at
    every time unit after R, like for like. Per block: the route difference in persistence (the gate), each part's
    retention, the push re-weighted to the medium's composition at R, and the share of the difference it closes."""
    S_, M_ = mask.shape

    def C(zz):
        return np.abs(zz).mean(-2)

    def grid(fn):
        return np.array([[fn(rows[i * M_ + m]) for m in range(M_)] for i in range(S_)], float)

    diffR = grid(lambda r: C(r["z_medium"][0])[WINDOW] - C(r["z_medium"][1])[WINDOW])

    def gate(lo, hi):
        wm = grid(lambda r: (C(r["z_medium"][0]) - C(r["z_medium"][1]))[WINDOW + lo + 1:WINDOW + hi + 1].mean())
        wp = grid(lambda r: (C(r["z_push"][0]) - C(r["z_push"][1]))[lo + 1:hi + 1].mean())
        return ratio_n(c3.by_start((wm - wp) * mask), c3.by_start(diffR * mask))

    def extras(key, offset):
        out = np.full((S_, M_, 3, AFTER + 1), np.nan)
        for i in range(S_):
            for m in range(M_):
                p = rows[i * M_ + m][key][..., offset:offset + AFTER + 1]   # (direction, group, part, time)
                out[i, m] = (p[0] - p[1]).mean(0) / 2
        return c3.by_start(out * mask[:, :, None, None])                    # (start, part, time)

    E = dict(medium=extras("parts_medium", WINDOW), push=extras("parts_push", 0))
    keep = np.isfinite(E["medium"]).all(axis=(1, 2)) & np.isfinite(E["push"]).all(axis=(1, 2))
    E = {route: value[keep] for route, value in E.items()}
    n = int(keep.sum())
    draws = np.random.default_rng(c3.SEED).integers(0, n, (c3.RESAMPLES, n)) if n >= 2 else []

    def numbers(lo, hi, pick):
        at = {r: E[r][pick, :, 0].mean(0) for r in E}
        over = {r: E[r][pick, :, lo + 1:hi + 1].mean(-1).mean(0) for r in E}
        total = {r: float(over[r].sum() / at[r].sum()) for r in E}
        kept = {r: over[r] / at[r] for r in E}
        weights = at["medium"] / at["medium"].sum()
        reweighted = float((weights * kept["push"]).sum())
        return dict(medium=total["medium"], push=total["push"], reweighted_push=reweighted,
                    closed=(reweighted - total["push"]) / (total["medium"] - total["push"]),
                    slower_kept_medium=float(kept["medium"][0]), slower_kept_push=float(kept["push"][0]),
                    slower_kept_difference=float(kept["medium"][0] - kept["push"][0]),
                    kept_medium=kept["medium"].tolist(), kept_push=kept["push"].tolist(),
                    medium_shares_at_R=weights.tolist(), push_shares_at_R=(at["push"] / at["push"].sum()).tolist())

    out = []
    with warnings.catch_warnings(), np.errstate(divide="ignore", invalid="ignore"):
        warnings.simplefilter("ignore", category=RuntimeWarning)
        for lo, hi in list(BLOCKS) + [tuple(LATER)]:
            g = gate(lo, hi)
            if n < 2:
                out.append(dict(block=f"(R+{lo}, R+{hi}]", route_difference=g, reading="no estimate", starts=n))
                continue
            point = numbers(lo, hi, slice(None))
            boot = [numbers(lo, hi, dr) for dr in draws]
            bars = {key: float(np.std([bb[key] for bb in boot], ddof=1))
                    for key in ("medium", "push", "reweighted_push", "closed", "slower_kept_medium", "slower_kept_push",
                                "slower_kept_difference")}
            if not g["clear"]:
                reading = "nothing clear to explain"
            elif point["slower_kept_difference"] > CLEAR * bars["slower_kept_difference"]:
                reading = "composition alone doesn't account for it"
            elif point["closed"] > 0.5:
                reading = "consistent with composition"
            else:
                reading = "not clear"
            out.append(dict(block=f"(R+{lo}, R+{hi}]", route_difference=g, starts=n, **point, bars=bars,
                            reading=reading))
    return out


# ---------------------------------------------------------------------------------------------------------------------
# the readings

def replication_readings(parts_dir, size, title, base_hash=None, moment_hash=None, dry=False):
    """Every reading fixed in plan.md (revision 1), from the fresh starts 115 to 134 only. Returns the record and the
    report's lines. dry: made-up numbers, for trying the code (their consistency checks fail, as they should, and are
    left out of the verdicts)."""
    numbers, inside = fluidity_numbers(STARTS, lambda s: parts_dir / f"base_{s:03d}.npz",
                                       lambda s, m: parts_dir / f"moment_{s:03d}_{m:02d}.npz", base_hash, moment_hash)
    failed = [f"{key} (largest difference {value:.1e})" for key, value in numbers["reproduced"].items()
              if not value <= TOL]
    rows = []
    for s in STARTS:
        for m in c3.MOMENTS:
            with np.load(parts_dir / f"moment_{s:03d}_{m:02d}.npz") as p:
                rows.append({key: p[key] for key in p.files})
    recorded = dict(check=max(float(r["parts_check"]) for r in rows),
                    adds_up=max(float(r["parts_adds_up"]) for r in rows))
    if not (recorded["check"] <= TOL and recorded["adds_up"] <= TOL):
        failed.append("the recorded parts against Check 3's calculation, or adding up")
    finite = all(np.isfinite(r["z_none"]).all() and np.isfinite(r["logs_none"]).all()
                 and np.isfinite(r["parts_medium"]).all()
                 and all(np.isfinite(r[key][d]).all() for key in ("z_push", "logs_push", "z_sham", "logs_sham",
                                                                  "parts_push", "parts_sham")
                         for d in range(2) if r["built"][d]) for r in rows)
    if not finite:
        failed.append("every number finite (untouched, medium, every built push and sham, and their recordings)")
    secondary = secondary_readings(rows, size)
    if not (secondary["checks"]["restored"] and secondary["checks"]["medium_finite"]):
        failed.append("the pull restored at R, and every medium number finite")
    if not secondary["checks"]["built_checks"]:
        failed.append("every built push and sham within its tolerances")
    counted = [] if dry else failed

    S_, M_ = len(STARTS), MOMENTS
    built, mask = inside["built"], inside["mask"]
    both = built.all(-1)
    missing = int((~both).sum())
    without = [STARTS[i] for i in range(S_) if not both[i].any()]
    d, sh = inside["d"], inside["sh"]
    with warnings.catch_warnings():
        warnings.simplefilter("ignore", category=RuntimeWarning)
        extras = {route: float(np.nanmean(d[route]["inner"]) + np.nanmean(d[route]["outer"]))
                  for route in ("medium", "push")}

    def grid(fn):
        return np.array([[fn(rows[i * M_ + m]) for m in range(M_)] for i in range(S_)], float)

    # beside hypothesis 1: the exact contributions, like for like
    exact_parts = {route: [c3.by_start(grid(lambda r, key=key, at=at, p=p:
                                            float(np.mean(r[key][0, :, p, at] - r[key][1, :, p, at])) / 2) * mask)
                           for p in range(3)]
                   for route, key, at in (("medium", "parts_medium", WINDOW), ("push", "parts_push", 0))}
    exact_sh = {route: share_n(v[0], v[1], S_) for route, v in exact_parts.items()}
    exact_gap = compared_n(exact_sh["medium"], exact_sh["push"])
    # beside hypothesis 1: all 200 planned moments, the push counted like the medium where not like for like
    medium_all = {key: grid(lambda r, key=key: float(np.mean(r[f"desc_{key}"][1] - r[f"desc_{key}"][2])) / 2)
                  for key in ("inner", "outer")}
    push_all = {key: np.where(both, grid(lambda r, key=key: float(np.mean(r[f"desc_{key}"][3] - r[f"desc_{key}"][4]))
                                         / 2), medium_all[key]) for key in ("inner", "outer")}
    conservative_sh = dict(medium=share_n(c3.by_start(medium_all["inner"]), c3.by_start(medium_all["outer"]), S_),
                           push=share_n(c3.by_start(push_all["inner"]), c3.by_start(push_all["outer"]), S_))
    conservative_gap = compared_n(conservative_sh["medium"], conservative_sh["push"])

    gap1 = numbers["question_2"]["medium_minus_push"]
    h1 = verdict(gap1["mean"], gap1["bar"], h1_problems(missing, without, extras, exact_gap, counted))
    gap2 = numbers["question_3"]["main"]
    h2 = verdict(gap2["mean"], gap2["bar"], list(counted))

    C = lambda zz: np.abs(zz).mean(-2)  # noqa: E731
    dC = grid(lambda r: C(r["z_medium"][0])[WINDOW] - C(r["z_medium"][1])[WINDOW])
    selection = dict(like_for_like=float(np.mean(dC[both])) if both.any() else NAN, n_like_for_like=int(both.sum()),
                     others=float(np.mean(dC[~both])) if (~both).any() else NAN, n_others=int((~both).sum()))
    comp = composition(rows, mask)
    pilot = json.loads(mt.PILOT_JSON.read_text())
    secondary_beside = [dict(name=name, replication=secondary[group][key], pilot=pilot[group][key],
                             described=beside(secondary[group][key], pilot[group][key]))
                        for name, group, key in SECONDARY_ITEMS]
    result = dict(
        title=title, starts=list(STARTS), size=size,
        hypothesis_1=dict(**h1, medium_share=as_dict_n(sh["medium"]), push_share=as_dict_n(sh["push"]),
                          like_for_like=PLANNED - missing, not_like_for_like=missing, starts_without=without,
                          extras=extras),
        hypothesis_2=dict(**h2, slower=numbers["question_3"]["by_block"][0]["slower half"],
                          faster=numbers["question_3"]["by_block"][0]["faster half"]),
        beside_hypothesis_1=dict(exact=dict(medium_share=as_dict_n(exact_sh["medium"]),
                                            push_share=as_dict_n(exact_sh["push"]), gap=exact_gap),
                                 conservative=dict(medium_share=as_dict_n(conservative_sh["medium"]),
                                                   push_share=as_dict_n(conservative_sh["push"]), gap=conservative_gap),
                                 selection=selection, not_built=secondary["checks"]["not_built"]),
        failed_checks=failed, recorded_parts=recorded, check3=numbers, secondary=secondary,
        secondary_beside_pilot=secondary_beside, exploratory=comp)

    def sh_text(s):
        return "n/a" if s is None else f"{s[0]:.3f} +- {float(np.std(s[1], ddof=1)):.3f}"

    def gap_text(v):
        return (f"{v['gap']:+.3f} +- {v['bar']:.3f}; lower bound {v['lower']:+.3f}, upper bound {v['upper']:+.3f}"
                if np.isfinite(v["gap"]) and np.isfinite(v["bar"]) else "n/a")

    q1, q2, q3 = numbers["question_1"], numbers["question_2"], numbers["question_3"]
    lines = [
        title,
        "A toy model: structure inside the model, not evidence about cells. Histories, not memory: differences between "
        "routes are between intervention histories and the states they leave.",
        f"Fresh starts 115 to 134: {S_} starts, {M_} moments each, {PLANNED} planned moments. Every number comes from "
        "these starts; the pilot's values appear only beside them, for comparison.",
        "Reproduced before reading on: " + "; ".join(f"{key} (largest difference {value:.1e})"
                                                     for key, value in numbers["reproduced"].items())
        + f"; the recorded parts against Check 3's calculation (largest difference {recorded['check']:.1e}) and adding"
          f" up (largest {recorded['adds_up']:.1e}).",
        f"Failed checks: {'; '.join(failed) if failed else 'none'}"
        + (" (made-up numbers: expected, and left out of the verdicts)" if dry and failed else ""),
        "",
        "Primary verdicts. Each hypothesis is judged on its own, with no combined verdict. Smallest effect of interest "
        "0.10; lower bound = gap - 3 bars; upper bound = gap + 2 bars. PASS: lower bound above zero and gap at least "
        "0.10. FAIL: upper bound below 0.10. Otherwise INCONCLUSIVE.",
        f"  Hypothesis 1, coherence distribution (the share gap): {h1['verdict']} ({h1['kind']})",
        f"    slower half's share of the extra at R: medium {sh_text(sh['medium'])}; push {sh_text(sh['push'])}",
        f"    gap (medium minus push): {gap_text(h1)}",
        f"    like-for-like moments: {PLANNED - missing} of {PLANNED} ({missing} not; more than {MOST_MISSING} would make "
        f"it inconclusive); starts without any: {', '.join(map(str, without)) or 'none'}",
        f"  Hypothesis 2, persistence (the persistence gap): {h2['verdict']} ({h2['kind']})",
        f"    coefficients over (t, t + 100]: slower half {est(q3['by_block'][0]['slower half'])}; faster half "
        f"{est(q3['by_block'][0]['faster half'])}",
        f"    gap (slower minus faster): {gap_text(h2)}",
        "",
        "Beside hypothesis 1 (reported; they don't decide its verdict)",
        f"  from the exact contributions: medium {sh_text(exact_sh['medium'])}; push {sh_text(exact_sh['push'])}; gap "
        f"{est(exact_gap)}; clearly opposite in sign: "
        f"{'yes' if np.isfinite(exact_gap['mean']) and exact_gap['mean'] + CLEAR * exact_gap['bar'] < 0 else 'no'}",
        f"  conservative, all {PLANNED} planned moments (the push counted like the medium where not like for like): "
        f"medium {sh_text(conservative_sh['medium'])}; push {sh_text(conservative_sh['push'])}; gap "
        f"{est(conservative_gap)}",
        f"  the medium's raised-minus-lowered coherence difference at R, averaged over moments: like-for-like moments "
        f"{selection['like_for_like']:+.5f} ({selection['n_like_for_like']}); others "
        + (f"{selection['others']:+.5f} ({selection['n_others']})" if selection["n_others"] else "none"),
        "  directions not built: " + (f"{len(secondary['checks']['not_built'])}, with reasons:"
                                      if secondary["checks"]["not_built"] else "none"),
        *(f"    {x}" for x in secondary["checks"]["not_built"]),
        "",
        "Secondary: delayed perturbation growth, by the pilot's own analysis function (descriptions, never PASS or "
        "FAIL), each beside the pilot's value",
        *(f"  {row['name']}: {est(row['replication'], 5)}; pilot {est(row['pilot'], 5)}: {row['described']}"
          for row in secondary_beside),
        f"  shapes: medium {secondary['shapes']['medium']} (pilot {pilot['shapes']['medium']}); pushes against shams "
        f"{secondary['shapes']['push_sham']} (pilot {pilot['shapes']['push_sham']}); pushes against untouched "
        f"{secondary['shapes']['push_untouched']} (pilot {pilot['shapes']['push_untouched']})",
        f"  the medium's S ruled out at half the benchmark: {'yes' if secondary['readings']['medium_ruled_out'] else 'no'}"
        f" (pilot {'yes' if pilot['readings']['medium_ruled_out'] else 'no'})",
        f"  readings: {', '.join(secondary['patterns']) or 'none'} (pilot: {', '.join(pilot['patterns']) or 'none'})",
        *(["  More than 1 moment in 10 lacked a matched push or sham: the medium-vs-push readings are inconclusive."]
          if secondary["readings"]["push_inconclusive"] else []),
        "",
        "Exploratory: does slower-half concentration explain the route difference in persistence? (exploratory, not a "
        "test; the re-weighting assumes each half's extra fades the same way whatever the other half holds)",
        *(f"  {row['block']}: route difference {est(row['route_difference'])}"
          + ("" if "medium" not in row else
             f"; persistence medium {row['medium']:.3f}, push {row['push']:.3f}, push re-weighted "
             f"{row['reweighted_push']:.3f} +- {row['bars']['reweighted_push']:.3f}; share of the difference closed "
             f"{row['closed']:+.2f} +- {row['bars']['closed']:.2f}; slower half kept: medium "
             f"{row['slower_kept_medium']:.3f}, push {row['slower_kept_push']:.3f} (difference "
             f"{row['slower_kept_difference']:+.3f} +- {row['bars']['slower_kept_difference']:.3f})")
          + f": {row['reading']}" for row in comp),
        "",
        "Also reported (descriptions)",
        f"  where natural order sits (Check 3's question 1): {numbers['verdicts']['where']}; shares: "
        + "; ".join(f"{name} {est(q1['shares'][name])}" for name in PART_NAMES)
        + f"; slower half minus an even spread ({c3.EVEN:.3f}): {est(q1['slower_minus_even'])}",
        f"  the natural share beside the routes (question 2): natural {sh_text(sh['natural'])}; medium "
        f"{sh_text(sh['medium'])}; push {sh_text(sh['push'])}: {numbers['verdicts']['routes']}; from states at R alone "
        f"{est(q2['natural_at_R_alone'])}, at T alone {est(q2['natural_at_T_alone'])}",
        "  hypothesis 2's coefficients, block by block (slower half; faster half; free turners; slower minus faster):",
        *(f"    {row['block']}: {est(row['slower half'])}; {est(row['faster half'])}; {est(row['free turners'])}; "
          f"{est(row['slower_minus_faster'])}" for row in q3["by_block"]),
        "  share of a coherence excess still there (natural slope; medium; push):",
        *(f"    {row['block']} after t or R: {est(row['natural'])}; {est(row['medium'])}; {est(row['push'])}"
          for row in q3["survival"]),
        "  arms separately against untouched (slower half's share of the change): "
        + "; ".join(f"{name} {est(v['share'])}" for name, v in q2["separately"].items()),
        f"  the medium over all {PLANNED} moments: slower half's share {est(q2['medium_all_moments'])}",
        "  states at R (coherence; slower half; faster half; |z2|; |z3|; |z4|; A-B timing):",
        *(f"    {state}: " + "; ".join(f"{secondary['descriptions']['states'][state][key]:.4f}"
                                       for key in ("R", "inner", "outer", "z2", "z3", "z4", "ab"))
          for state in mt.STATES),
        "  phase distance between the medium's and the push's states at R: "
        + "; ".join(f"{name} {est(v)}" for name, v in secondary["descriptions"]["distance"].items()),
        "",
        "Limits",
        "  A toy model at one pull, size and window, with one push method. Natural states within a start share one "
        "trajectory; the start is the unit throughout.",
        "  A half's own coherence responds only weakly to small amounts of order near the random-phase floor; the exact "
        "contributions beside hypothesis 1 don't have that limit.",
        "  The exploratory reading isn't a controlled comparison: no arm changes the composition with everything else "
        "fixed, and the routes also differ in A-B timing and phase arrangement at R."]
    return result, lines


def made_up(folder):
    """Files shaped like the run's, with made-up numbers (no model run), for trying the report only: base runs with
    random states and the pilot's made-up moments with made-up recordings (some directions not built)."""
    rng = np.random.default_rng(11)
    with starts_for_analysis(SEEDS, MOST_MISSING):
        rows = mt.made_up_rows(frozen_size())
    samples = int(mt.RUN) + 1
    for s in SEEDS:
        times = mt.keep_times(s)
        u = np.exp(1j * rng.vonmises(0.0, 0.08, (len(times), 2, 2, N)))
        z = 0.03 * (1 + 0.3 * rng.standard_normal((2, samples))) * np.exp(1j * rng.uniform(-np.pi, np.pi, (2, samples)))
        save_npz(folder / f"base_{s + 1:03d}.npz", start=s + 1, z=z, logs=0.07 + 0.02 * rng.standard_normal(samples // 5),
                 kept_t=np.array(times, float), kept_u=u[..., :c3.PULLED], kept_f=u[..., c3.PULLED:],
                 trace=np.zeros((2, 61), np.float32))
    for row in rows:
        def parts(length, built=None):
            p = 0.01 * (1 + rng.standard_normal((2, 2, 3, length)))
            if built is not None:
                p[~built] = np.nan
            return p
        row.update(parts_medium=parts(WINDOW + AFTER + 1), parts_push=parts(AFTER + 1, row["built"]),
                   parts_sham=parts(AFTER + 1, row["built"]), parts_none_T=rng.uniform(0, 0.02, (2, 3)),
                   parts_none_R=rng.uniform(0, 0.02, (2, 3)), parts_check=0.0, parts_adds_up=0.0)
        save_npz(folder / f"moment_{int(row['start']):03d}_{int(row['moment']):02d}.npz", **row)


# ---------------------------------------------------------------------------------------------------------------------
# the commands

def checkruns():
    if LOCK.exists():
        sys.exit("Already locked: the checks before the lock are done.")
    found = outputs()
    if found:
        sys.exit(f"Replication output exists ({found[0]} ...): the checks before the lock can't run now.")
    problems = pilot_freeze_problems()
    if problems:
        sys.exit(f"Stopping: {'; '.join(problems)}.")
    if CHECK_PARTS.exists():
        shutil.rmtree(CHECK_PARTS)                         # a fresh set of check runs each time
    for sub in ("rerun", "pilot-code", "new-code", "standin"):
        (CHECK_PARTS / sub).mkdir(parents=True)
    saved = mt.PILOT_PARTS / f"base_{CHECK_SEED + 1:03d}.npz"
    for sub in ("pilot-code", "new-code"):                 # both codes read the pilot's saved base run (read only)
        (CHECK_PARTS / sub / saved.name).symlink_to(saved)
    run_piece(CHECK_LOG, checkrun_tasks(frozen_size()), job, 265.0, 30.0)
    log = json.loads(CHECK_LOG.read_text())
    log["pieces"][-1].update(plan_sha256=sha256(PLAN), code_sha256=sha256(CODE))
    write_json(CHECK_LOG, log)


def checks():
    if LOCK.exists():
        sys.exit("Already locked: the checks before the lock are done.")
    began = time.time()
    size = frozen_size()
    lines = [f"Checks before the lock: Fluidity Replication (plan.md, revision 1), {shown(stamp())}.",
             "No fresh start was touched. Equality and mechanics only: no outcome of any intervention is read, shown or"
             " kept. The pilot's own numbers are used only to check the code.", ""]
    results, details = {}, {}

    def note(key, ok, text, detail=None):
        results[key] = bool(ok)
        details[key] = detail
        lines.append(f"{text}: {'yes' if ok else 'NO'}")

    now = json.loads(json.dumps(plain(SETTINGS)))
    planned = json.loads(json.dumps(plain(PLAN_NUMBERS)))
    wrong = sorted({key for key in planned if now.get(key) != planned[key]} | {key for key in now if key not in planned})
    halves = bool(np.array_equal(mt.inner_half(), c3.inner_half()))
    note("0 settings", not wrong and halves,
         "0. The runner's settings match the plan's numbers (tolerances, starts, moment rule, seeds, size, verdict rule),"
         " and the pilot's code and Check 3 define the same slower half" + (f" [differ: {', '.join(wrong)}]" if wrong
                                                                             else ""), dict(differ=wrong, halves=halves))

    log = json.loads(CHECK_LOG.read_text()) if CHECK_LOG.exists() else None
    names = [t["name"] for t in checkrun_tasks(size)]
    done = bool(log) and all(name in log["parts"] for name in names)
    under = bool(log) and all(p.get("plan_sha256") == sha256(PLAN) and p.get("code_sha256") == sha256(CODE)
                              for p in log["pieces"])
    seconds = sum(p.get("seconds", 0) for p in log["pieces"]) if log else 0
    note("1 check runs", done and under, f"1. The check runs finished ({len(log['pieces']) if log else 0} piece, "
         f"{seconds:.0f} s), all under this plan and this code", dict(done=done, under=under))
    if not done:
        lines += ["", "All checks passed: NO (the check runs are missing: run checkruns first)"]
        CHECKS_TXT.write_text("\n".join(lines) + "\n")
        write_json(CHECKS_JSON, dict(made=stamp(), all_passed=False, results=results, details=details))
        sys.exit("\n".join(lines))

    runlog = json.loads(mt.PILOT_LOG.read_text())["parts"]
    report = json.loads(mt.PILOT_JSON.read_text())
    saved_base = mt.PILOT_PARTS / "base_103.npz"
    with np.load(saved_base) as a, np.load(CHECK_PARTS / "rerun" / "base_103.npz") as b:
        same_base = {key: arrays_equal(a[key], b[key]) for key in ("start", "z", "logs", "trace", "kept_t", "kept_u",
                                                                   "kept_f")}
        unbroken, base_seconds = int(b["pieces"]) == 1, float(b["spent"])
    saved_ok = runlog["base_103"]["sha256"] == sha256(saved_base)
    note("1a base run", all(same_base.values()) and unbroken and saved_ok,
         f"   a. Pilot start 103's base run, re-run unbroken now ({base_seconds:.0f} s, one piece), against the saved "
         f"file (which matches its recorded fingerprint), every array ({', '.join(same_base)}): identical to the last "
         "digit", dict(arrays=same_base, unbroken=unbroken, saved_matches_fingerprint=saved_ok))

    saved_moment = mt.PILOT_PARTS / "moment_103_01.npz"
    with np.load(saved_moment) as a, np.load(CHECK_PARTS / "pilot-code" / "pilot_code_moment_103_01.npz") as b, \
            np.load(CHECK_PARTS / "new-code" / "new_code_moment_103_01.npz") as c:
        keys = a.files
        old_saved = all(k in b.files and arrays_equal(a[k], b[k]) for k in keys) and sorted(b.files) == sorted(keys)
        new_saved = all(k in c.files and arrays_equal(a[k], c[k]) for k in keys)
        new_old = all(arrays_equal(b[k], c[k]) for k in keys)
        added = sorted(set(c.files) - set(keys))
        T = int(a["T"])
        built = [bool(x) for x in c["built"]]
        recording = dict(check=float(c["parts_check"]), adds_up=float(c["parts_adds_up"]),
                         shapes=bool(c["parts_medium"].shape == (2, 2, 3, WINDOW + AFTER + 1)
                                     and c["parts_push"].shape == (2, 2, 3, AFTER + 1)
                                     and c["parts_sham"].shape == (2, 2, 3, AFTER + 1)
                                     and c["parts_none_T"].shape == (2, 3) and c["parts_none_R"].shape == (2, 3)),
                         finite=bool(np.isfinite(c["parts_medium"]).all()
                                     and all(np.isfinite(c[k][d]).all() for k in ("parts_push", "parts_sham")
                                             for d in range(2) if built[d])))
    saved_moment_ok = report["moment_sha256"]["moment_103_01.npz"] == sha256(saved_moment)
    note("1b moment", old_saved and new_saved and new_old and saved_moment_ok,
         f"   b. Its moment 1 (T = {T}), run now by the pilot's own code and by the new runner from the saved base run, "
         f"against the saved file, every array the pilot saved ({len(keys)}): pilot code = saved "
         f"{'yes' if old_saved else 'NO'}; new runner = saved {'yes' if new_saved else 'NO'}; new runner = pilot code "
         f"{'yes' if new_old else 'NO'}; identical to the last digit",
         dict(pilot_code_equals_saved=old_saved, new_runner_equals_saved=new_saved, new_equals_pilot_code=new_old,
              arrays=len(keys), added=added, saved_matches_fingerprint=saved_moment_ok, pushes_built=built))

    stand = json.loads((CHECK_PARTS / "standin" / "standin.json").read_text())
    arms_identical = all(v for part in stand.values() for k, v in part.items() if k.endswith("recording_arm_identical"))
    arms_counted = sum(1 for part in stand.values() for k in part if k.endswith("recording_arm_identical"))
    note("2a non-invasive", new_old and arms_identical and arms_counted == 6,
         f"2. The added recording.\n   a. Non-invasive: the new runner's arrays are identical to the pilot code's (1b), "
         f"and on old start 1 each arm with the recording matches the same arm without it to the last digit "
         f"({arms_counted} arms: a zero change, a changed pull, two pushes, two shams)",
         dict(moment=new_old, standin=arms_identical, arms=arms_counted))
    stand_adds = all(v for part in stand.values() for k, v in part.items() if k.endswith("parts_add_up"))
    note("2b adds up", recording["adds_up"] <= TOL and stand_adds and recording["shapes"] and recording["finite"],
         f"   b. Adds up: in every arm of moment 1, at every recorded time, the three parts add up to the group's "
         f"coherence (largest difference {recording['adds_up']:.1e}; at most 1e-12), likewise on old start 1; every "
         f"recorded array has its planned shape and is finite where built",
         dict(largest=recording["adds_up"], standin=stand_adds, shapes=recording["shapes"], finite=recording["finite"]))
    stand_match = stand["moment task on the stand-in"]["parts_match_check3"]
    note("2c correct", recording["check"] <= TOL and stand_match,
         f"   c. Correct: at R, the recorded parts against Check 3's own calculation from the states (largest "
         f"difference {recording['check']:.1e}; at most 1e-12), likewise on the stand-in",
         dict(largest=recording["check"], standin=stand_match))

    zero, changed = stand["zero change"], stand["pull changed 5%, then restored"]
    pushes = [stand["push raising"], stand["push lowering"]]
    moment = stand["moment task on the stand-in"]
    mech = dict(zero_reproduces_untouched=zero["reproduces_untouched"], pull_as_scheduled=changed["pull_as_scheduled"],
                change_took_effect=changed["change_took_effect"],
                pushes_built=all(p["built"] for p in pushes), on_target=all(p["on_target"] for p in pushes),
                sham_in_band=all(p["sham_in_band"] for p in pushes), same_sizes=all(p["same_sizes"] for p in pushes),
                gap_kept=all(p["gap_kept"] for p in pushes),
                moment_task_identical_to_pilot_code=moment["pilot_arrays_identical"], moment_shapes=moment["shapes"],
                pull_restored_after_R=moment["restored"], moment_problem=moment["problem"] == "none")
    note("3 mechanics", all(mech.values()),
         "3. Mechanics on old start 1 (short runs, as in the pilot's checks): a zero change reproduces the untouched run "
         "exactly; the pull is exactly as scheduled and back at exactly 1.9 after the change; both pushes built, on "
         "target, shams in band, sizes the same, gap kept; the replication's moment task on a stand-in moment gives "
         "every array the pilot's does, identically, with the pull restored at R", mech)

    began_report = time.time()
    hashes = report["moment_sha256"]
    mine, _ = fluidity_numbers(PILOT_STARTS, lambda s: mt.PILOT_PARTS / f"base_{s:03d}.npz",
                               lambda s, m: mt.PILOT_PARTS / f"moment_{s:03d}_{m:02d}.npz",
                               lambda p: runlog.get(p.stem, {}).get("sha256"), lambda p: hashes.get(p.name), report)
    mine = as_json(mine)
    recorded3 = json.loads((PILOT_DIR / "check3_fluidity.json").read_text())
    compared_keys = ("verdicts", "reproduced", "question_1", "question_2", "question_3")
    count3, largest3 = number_misfit({k: mine[k] for k in compared_keys}, {k: recorded3[k] for k in compared_keys})
    note("4a check 3", count3 > 0 and largest3 == 0.0,
         f"4. The report code reproduces the pilot's records.\n   a. Check 3's calculations, pointed at the pilot's "
         f"files: check3_fluidity.json's verdicts and numbers exactly ({count3} numbers; largest difference "
         f"{largest3:.1e})", dict(numbers=count3, largest=largest3))
    pilot_rows = mt.load_pilot(mt.PILOT_PARTS)
    res = as_json(plain(secondary_readings(pilot_rows, size, mt.PILOT_STARTS, 12)))
    put_back = mt.PILOT_STARTS == tuple(range(102, 114)) and mt.MOST_MISSING == 12
    count2, largest2 = number_misfit(res, {k: report[k] for k in res})
    note("4b analysis", count2 > 0 and largest2 == 0.0 and put_back,
         f"   b. The pilot's analysis function, called the way the runner calls it (its starts and cap set for the call "
         f"and put back after{'' if put_back else ' [NOT put back]'}): pilot_results.json's readings exactly ({count2} "
         f"numbers; largest difference {largest2:.1e})",
         dict(numbers=count2, largest=largest2, put_back=put_back))
    report_seconds = time.time() - began_report

    wrong_verdicts = verdict_problems()
    note("5 verdicts", not wrong_verdicts,
         f"5. The verdict code gives the planned verdict in each of the plan's examples and each validity condition "
         f"({len(VERDICT_CASES) + len(CONDITION_CASES) + 1} cases)" + (f" [wrong: {wrong_verdicts}]" if wrong_verdicts
                                                                        else ""), wrong_verdicts)

    began_dry = time.time()
    if DRY.exists():
        shutil.rmtree(DRY)
    (DRY / "parts").mkdir(parents=True)
    made_up(DRY / "parts")
    dry, dry_lines = replication_readings(DRY / "parts", size, "DRY RUN ON MADE-UP NUMBERS (no model run; nothing here "
                                          "is a result).", dry=True)
    (DRY / "dryreport.txt").write_text("\n".join(dry_lines) + "\n")
    write_json(DRY / "dryreport.json", dry)
    shutil.rmtree(DRY / "parts")
    dry_ok = bool(dry_lines) and dry["hypothesis_1"]["verdict"] in ("PASS", "FAIL", "INCONCLUSIVE") \
        and dry["hypothesis_2"]["verdict"] in ("PASS", "FAIL", "INCONCLUSIVE")
    note("6 dry report", dry_ok,
         f"6. The whole report on made-up numbers (written to {DRY} only; {len(dry_lines)} lines; "
         f"{dry['hypothesis_1']['not_like_for_like']} of {PLANNED} moments not like for like; "
         f"{time.time() - began_dry:.0f} s; meaningless numbers): ran end to end",
         dict(lines=len(dry_lines), not_like_for_like=dry["hypothesis_1"]["not_like_for_like"]))

    freeze = pilot_freeze_problems()
    data = pilot_data_problems()
    records = pilot_records()
    note("7 nothing earlier changed", not freeze and not data,
         f"7. Nothing earlier changed: the pilot's freeze record verifies (its protected files, plan, code and "
         f"calibration results); the pilot's saved data match their fingerprints ({len(PILOT_STARTS) * (MOMENTS + 1)} "
         f"files); the cause test's lock holds; the pilot's records ({len(records)} files) fingerprinted for the lock"
         + (f" [{'; '.join(freeze + data)}]" if freeze or data else ""), dict(freeze=freeze, data=data))
    found = outputs()
    note("8 no output", not found, "8. No replication output exists" + (f" [{found[0]} ...]" if found else ""), found)

    all_passed = all(results.values())
    lines += ["", f"All checks passed: {'yes' if all_passed else 'NO'}",
              f"Time: the check runs took {seconds:.0f} s in one piece; these comparisons {time.time() - began:.0f} s "
              f"(the report code on the pilot's files {report_seconds:.0f} s).",
              f"Run under plan.md SHA-256 {sha256(PLAN)} and replication.py SHA-256 {sha256(CODE)}."]
    write_json(CHECKS_JSON, dict(made=stamp(), all_passed=all_passed, results=results, details=details,
                                 sha256=dict(plan=sha256(PLAN), runner=sha256(CODE)), settings=SETTINGS,
                                 machine=machine(), check_runs_machine=log.get("machine"), pilot_records=records))
    CHECKS_TXT.write_text("\n".join(lines) + "\n")
    print("\n".join(lines))


def lock():
    if LOCK.exists():
        made = json.loads(LOCK.read_text())["made"]
        sys.exit(f"Already locked ({shown(made)}). A new lock needs Colson's approval and is labelled as a change.")
    found = outputs()
    if found:
        sys.exit(f"Replication output exists ({found[0]} ...): the lock is refused.")
    if not CHECKS_JSON.exists():
        sys.exit("No check results: run the checks first.")
    record = json.loads(CHECKS_JSON.read_text())
    if not record["all_passed"]:
        sys.exit("A check failed: the lock is refused.")
    if record["sha256"] != dict(plan=sha256(PLAN), runner=sha256(CODE)) or \
            record["settings"] != json.loads(json.dumps(plain(SETTINGS))):
        sys.exit("The plan, the code or the settings changed after the checks ran: run the checks again.")
    problems = pilot_freeze_problems() + pilot_data_problems()
    if pilot_records() != record["pilot_records"]:
        problems.append("the pilot's records changed after the checks ran")
    if problems:
        sys.exit(f"Stopping: {'; '.join(problems)}.")
    fingerprints = {name: sha256(p) for name, p in lock_files().items()}
    lock_record = dict(
        made=stamp(),
        what=("The lock of the Fluidity Replication (plan.md, revision 1), written after every check passed and before "
              "any fresh start ran. The run checks every fingerprint below before each piece, and the report before it "
              "reads anything. Nothing here changes afterwards; any change would come to Colson first and be a new, "
              "dated lock labelled as a change, with this record kept."),
        sha256=fingerprints, pilot_records=record["pilot_records"], settings=SETTINGS,
        starts=list(STARTS), moments={str(s + 1): dict(T=mt.moment_times(s), R=[T + WINDOW for T in mt.moment_times(s)])
                                      for s in SEEDS},
        verdict_rule=VERDICT_RULE, checks=dict(made=record["made"], all_passed=True, results=record["results"]),
        machine=machine(), check_runs_machine=record.get("check_runs_machine"), outputs_at_lock=found,
        pilot_freeze_verified=True, cause_lock_unchanged=True)
    write_json(LOCK, lock_record)
    print(f"Locked {shown(lock_record['made'])}: replication_lock.json (SHA-256 {sha256(LOCK)}).")
    for name, value in fingerprints.items():
        print(f"  {name}: {value}")
    print("The run on starts 115 to 134 waits for Colson's separate go.")


INCONCLUSIVE_STOP = (" Under the plan both hypotheses are INCONCLUSIVE: nothing is reported, and this comes to "
                     "Colson.")


def locked(need_go=True, at_report=False):
    """The lock record, after checking that everything it fixes is as it was, and (for the run and the report) that
    Colson's go is recorded. At the report, a mismatch makes both hypotheses INCONCLUSIVE (plan.md)."""
    after = INCONCLUSIVE_STOP if at_report else ""
    if not LOCK.exists():
        sys.exit("No lock record: nothing on a fresh start can run before the lock.")
    record = json.loads(LOCK.read_text())
    changed = [name for name, p in lock_files().items() if record["sha256"].get(name) != sha256(p)]
    if record["settings"] != json.loads(json.dumps(plain(SETTINGS))):
        changed.append("settings")
    if pilot_records() != record["pilot_records"]:
        changed.append("the pilot's records")
    changed += pilot_freeze_problems()
    if changed:
        sys.exit(f"Changed since the lock ({shown(record['made'])}): {', '.join(changed)}. Stopping.{after}")
    if need_go:
        if not GO.exists():
            sys.exit("No record of Colson's go (replication_go.json): the run waits for it.")
        if json.loads(GO.read_text()).get("lock_sha256") != sha256(LOCK):
            sys.exit(f"The go record doesn't refer to this lock. Stopping.{after}")
    return record


def go(words):
    if GO.exists():
        sys.exit("Colson's go is already recorded.")
    if not words.strip():
        sys.exit("The go needs Colson's words.")
    record = locked(need_go=False)
    write_json(GO, dict(made=stamp(), words=words, lock_made=record["made"], lock_sha256=sha256(LOCK)))
    print(f"Colson's go recorded {shown(stamp())} in replication_go.json, for the lock of {shown(record['made'])}.")


def run():
    record = locked()
    now = machine()
    if not all(same_environment(env, now) for env in (record["machine"], record["check_runs_machine"]) if env):
        sys.exit("The numerical environment differs from the one the checks ran on. Stopping.")
    base, moments = run_tasks(record["settings"]["size"])
    if RUN_LOG.exists():
        failures = {}
        for piece in json.loads(RUN_LOG.read_text())["pieces"]:
            for name, r in piece["results"].items():
                failures[name] = failures.get(name, 0) + (r["status"] == "failed")
        twice = sorted(name for name, count in failures.items() if count >= 2)
        if twice:
            sys.exit(f"Failed twice: {', '.join(twice)}. Stopping; this comes to Colson without any numbers read.")
    todo = mt.unfinished(base)
    ready = [t for t in mt.unfinished(moments) if (RUN_PARTS / f"base_{t['seed'] + 1:03d}.npz").exists()]
    if not todo and not ready:
        print("All 220 tasks are done: run report.")
        return
    run_piece(RUN_LOG, todo + ready, job, mt.BUDGET, mt.START_BY)
    log = json.loads(RUN_LOG.read_text())
    log["pieces"][-1].update(plan_sha256=sha256(PLAN), code_sha256=sha256(CODE))
    write_json(RUN_LOG, log)
    left = len(mt.unfinished(base)) + len(mt.unfinished(moments))
    print(f"{220 - left} of 220 tasks done.")


def report():
    if RESULTS_JSON.exists():
        made = json.loads(RESULTS_JSON.read_text())["made"]
        sys.exit(f"The report was made {shown(made)}; it is made once.")
    record = locked(at_report=True)
    if not RUN_LOG.exists():
        sys.exit("No run log: nothing has run.")
    log = json.loads(RUN_LOG.read_text())
    base, moments = run_tasks(record["settings"]["size"])
    if mt.unfinished(base) or mt.unfinished(moments):
        sys.exit("The run is incomplete: the report needs all 220 tasks.")
    hashes = {name: value["sha256"] for name, value in log["parts"].items()}
    bad = [t["name"] for t in base + moments if hashes.get(t["name"]) != sha256(RUN_PARTS / f"{t['name']}.npz")]
    if bad:
        sys.exit(f"Files differ from their logged fingerprints ({', '.join(bad[:3])} ...). Stopping.{INCONCLUSIVE_STOP}")
    if any(p.get("plan_sha256") != record["sha256"]["plan"] or p.get("code_sha256") != record["sha256"]["runner"]
           for p in log["pieces"]):
        sys.exit(f"A piece ran under a plan or code other than the locked ones. Stopping.{INCONCLUSIVE_STOP}")
    title = f"Fluidity Replication: results ({shown(stamp())})."
    first, lines = replication_readings(RUN_PARTS, record["settings"]["size"], title, lambda p: hashes.get(p.stem),
                                        lambda p: hashes.get(p.stem))
    second, _ = replication_readings(RUN_PARTS, record["settings"]["size"], title, lambda p: hashes.get(p.stem),
                                     lambda p: hashes.get(p.stem))
    if json.dumps(plain(first), sort_keys=True) != json.dumps(plain(second), sort_keys=True):
        sys.exit(f"The report gave different numbers on a second computation. Stopping.{INCONCLUSIVE_STOP}")
    seconds = sum(p.get("seconds", 0) for p in log["pieces"])
    go_record = json.loads(GO.read_text())
    lines += ["", f"Compute: {len(log['pieces'])} pieces, {seconds / 60:.1f} minutes of foreground time. The report was "
                  "computed twice, with identical numbers.",
              f"Plan, code, model files, the pilot's freeze and records as locked ({shown(record['made'])}); Colson's go "
              f"recorded {shown(go_record['made'])}; the cause test's locked files unchanged: yes."]
    first.update(made=stamp(), lock=dict(made=record["made"], sha256=sha256(LOCK)), go=go_record,
                 part_sha256={t["name"]: hashes[t["name"]] for t in base + moments})
    write_json(RESULTS_JSON, first)
    RESULTS_TXT.write_text("\n".join(lines) + "\n")
    print("\n".join(lines))


if __name__ == "__main__":
    commands = dict(checkruns=checkruns, checks=checks, lock=lock, run=run, report=report)
    if len(sys.argv) >= 2 and sys.argv[1] == "go":
        go(" ".join(sys.argv[2:]))
    elif len(sys.argv) == 2 and sys.argv[1] in commands:
        commands[sys.argv[1]]()
    else:
        sys.exit(__doc__)