#!/usr/bin/env python3
"""
Hold feasibility (prelock_checks_plan.md, part 2), October 13, 2026. On old starts 13 to 18, ten moments each: do held
kicks at ±25% over 500 and over 100 units put each group's coherence on target at every kick, and what history and
effort do they produce? The pulse at ±25% (the cause test's dense method) runs beside them for comparison. Only the
untouched first copy is used: no gap is carried and no growth is measured. Nothing here is locked.

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

Usage (each command runs in the foreground and ends by itself; nothing is left running)
    python3 -B prelock_feasibility.py dryrun   # the machinery on old start 1 for 700 units; writes only to /tmp
    python3 -B prelock_feasibility.py run      # one piece: the 6 base runs, or the next moment tasks, 7 at a time
    python3 -B prelock_feasibility.py report   # once all 126 tasks are done: the criteria and the outputs

Each is run with a hard limit: timeout -k 10 285 python3 -B prelock_feasibility.py <command>
"""
import json
import shutil
import sys
import time
from pathlib import Path

sys.dont_write_bytecode = True
import numpy as np  # noqa: E402

from prelock_common import (HERE, WORK, SEED, N, K, EPS, DT, EVERY, ct, nt, rule, stamp, shown,  # noqa: E402
                            sha256, write_json, save_npz, run_piece, guard, lock_changes, folder_fingerprints)
from prelock_calibration import (calibrate, reproduced, revised_v1, saved_difference_text,  # noqa: E402
                                 environment_text)  # the base run in pieces; V1 as preregistered and V1-revised

STARTS = tuple(range(12, 18))                              # old starts 13 to 18 (seeds 12 to 17)
MOMENTS, SPACING, FIRST = 10, 350, 2000                    # T = 2,000 + o + 350 m
STEP = int(EVERY)                                          # a kick every 5 units
HOLDS = (500, 100)
SIZE = 0.25                                                # ±25% of each group's untouched coherence
AFTER = 500                                                # untouched after T
BLOCK = 25                                                 # the block fractions
HALF_BAND, SHAM_BAND = 0.025, (0.95, 1.05)                 # up, down: the target ratio ± 0.025; sham: its own before
ARMS = ("up", "down", "sham")
CRITERIA = dict(most_missing=3, fewest_kept=5, achieved=(0.75, 1.25), achieved_mean=(0.85, 1.15), share=0.9,
                block=(0.5, 1.5), sham_effort=(0.8, 1.2), sham_history=(-0.25, 0.25))
WINDOWS = dict(early=(0, 200), later=(200, 500))           # the coherence left after release, over (T + a, T + b]
BUDGET, START_BY = 190.0, 120.0                            # seconds into a piece: a base run still going at BUDGET
                                                           # keeps its exact state and goes on in the next piece;
                                                           # begin no task after START_BY
PARTS = WORK / "feasibility"
RUNLOG = HERE / "feasibility_runlog.json"
OUT_TXT, OUT_JSON, OUT_NPZ = (HERE / f"feasibility_results.{x}" for x in ("txt", "json", "npz"))
DRY = Path("/tmp/prelock-feasibility-dryrun")


def bands(size):
    return dict(up=(1 + size - HALF_BAND, 1 + size + HALF_BAND), down=(1 - size - HALF_BAND, 1 - size + HALF_BAND),
                sham=SHAM_BAND)


def moment_table():
    """(seed, moment, T) for every moment: T = 2,000 + o + 350 m, o = 5 x an integer from 0 to 69 drawn once per
    start from (SEED, 2, seed)."""
    rows = []
    for s in STARTS:
        o = STEP * int(np.random.default_rng((SEED, 2, s)).integers(0, 70))
        rows += [(s, m, FIRST + o + SPACING * m) for m in range(MOMENTS)]
    return rows


# ---------------------------------------------------------------------------------------------------------------------
# a kick and a hold

def one_sided(u, along, z, target):
    """The signed a closest to 0 with R(a) = |z + (1/N) sum_j u_j (e^{i a along_j} - 1)| = target, on the side that
    moves toward the target (a > 0 raises R). |a| is scanned at pi/512, 2 pi/512, ..., pi/2 (the cause test's grid);
    the first point that reaches the target and the one before are halved between until less than 1e-12 apart.
    Returns (a, 0); (nan, 1) if the target is not reached by |a| = pi/2 (infeasible); (nan, 2) if a number is not
    finite or R misses the target by more than a relative 1e-6 (a solver failure)."""
    real, imag = u.real.copy(), u.imag.copy()
    start = abs(z)
    if not (np.isfinite(target) and target > 0 and np.isfinite(start) and np.all(np.isfinite(along))):
        return np.nan, 2
    if start == target:
        return 0.0, 0
    side = 1.0 if target > start else -1.0

    def R(a):
        angle = np.multiply.outer(np.atleast_1d(a), along)
        c, s = np.cos(angle) - 1, np.sin(angle)
        cr, ci, sr, si = c @ real, c @ imag, s @ real, s @ imag
        return np.abs(z + ((cr - si) + 1j * (ci + sr)) / N)

    grid = ct.LARGEST * np.arange(1, ct.SCAN + 1) / ct.SCAN
    values = side * (R(side * grid) - target)                         # at or above 0 once the target is reached
    if not np.all(np.isfinite(values)):
        return np.nan, 2
    reached = np.flatnonzero(values >= 0)
    if not len(reached):
        return np.nan, 1
    i = int(reached[0])
    low, high = (float(grid[i - 1]) if i else 0.0), float(grid[i])
    while high - low > ct.BISECT:
        middle = (low + high) / 2
        if side * (R(side * middle)[0] - target) >= 0:
            high = middle
        else:
            low = middle
    value = R(side * high)[0]
    if not (np.isfinite(value) and abs(value / target - 1) <= ct.SOLVED):
        return np.nan, 2
    return side * high, 0


def hold(u0, f0, T, H, R_none, key, size=SIZE, after=AFTER):
    """One hold from the untouched first copy's state at T - H (u0: (group, pulled), f0: (group, free)): up, down and
    sham side by side (rows 0, 1, 2), a kick at every multiple of 5 in [T - H, T], then untouched to T + after.
    R_none: the untouched run's R_g every time unit from t = 0, (group, time). Each kick is worked out from the row's
    own state: up and down by one_sided() to (1 +- size) R_g^none(t); the sham by the mean of their sizes with the
    cause test's paired signs (drawn from `key`, kick by kick, group A before B). If a kick is infeasible, unsolved or
    outside its band, the hold stops there and is not carried out; nothing is clipped or replaced."""
    slow, stepper = nt.model(N, K, EPS)
    count, kicks = u0.shape[-1], H // STEP + 1
    U, F = np.repeat(u0[None], 3, 0), np.repeat(f0[None], 3, 0)
    step, none = stepper(DT), np.zeros(3)
    draw = np.random.default_rng(key)
    per_unit = int(round(1.0 / DT))
    band = bands(size)
    nan = np.nan
    rec = dict(R=np.full((3, 2, H + after + 1), nan),             # R_g every unit from T - H (after the kick at kicks)
               before=np.full((3, 2, kicks), nan),                # R_g just before each kick
               amp=np.full((3, 2, kicks), nan),                   # up's and down's signed a; the sham's size
               status=np.zeros((2, 2, kicks), int), ratio=np.full((3, 2, kicks), nan),
               inband=np.zeros((3, 2, kicks), bool), sumsq=np.full((3, kicks), nan), largest=np.full((3, kicks), nan),
               kicks=kicks, done=0, carried=False, reason="")
    goals = (1 + size, 1 - size)
    for k in range(kicks):
        t = T - H + STEP * k
        z = ct.harmonic(U, F)
        R, psi = np.abs(z), np.angle(z)
        rec["before"][:, :, k] = R
        along = np.sin(psi[..., None] - np.angle(U))              # sin(psi_g - theta_j), (row, group, pulled)
        for g in range(2):
            for a in range(2):
                rec["amp"][a, g, k], rec["status"][a, g, k] = one_sided(U[a, g], along[a, g], z[a, g],
                                                                        goals[a] * R_none[g, t])
        if rec["status"][:, :, k].any():
            bad = [f"{ARMS[a]} {'AB'[g]} {'infeasible' if rec['status'][a, g, k] == 1 else 'unsolved'}"
                   for a in range(2) for g in range(2) if rec["status"][a, g, k]]
            rec["reason"] = f"kick {k + 1} at t = {t}: " + ", ".join(bad)
            return rec
        dtheta = np.zeros((3, 2, count))
        for g in range(2):
            b = 0.5 * (abs(rec["amp"][0, g, k]) + abs(rec["amp"][1, g, k]))
            signs = ct.paired_signs(along[2, g], draw)
            rec["amp"][2, g, k] = b
            dtheta[0, g] = rec["amp"][0, g, k] * along[0, g]
            dtheta[1, g] = rec["amp"][1, g, k] * along[1, g]
            dtheta[2, g] = signs * b * along[2, g]
        U *= np.exp(1j * dtheta)
        now = np.abs(ct.harmonic(U, F))
        rec["ratio"][:2, :, k] = now[:2] / R_none[:, t]
        rec["ratio"][2, :, k] = now[2] / R[2]
        for a, arm in enumerate(ARMS):
            rec["inband"][a, :, k] = (band[arm][0] <= rec["ratio"][a, :, k]) & (rec["ratio"][a, :, k] <= band[arm][1])
        rec["sumsq"][:, k] = np.square(dtheta).sum((-2, -1))
        rec["largest"][:, k] = np.abs(dtheta).max((-2, -1))
        rec["R"][:, :, STEP * k] = now
        rec["done"] = k + 1
        if not rec["inband"][:, :, k].all():
            bad = [f"{arm} {'AB'[g]} at {rec['ratio'][a, g, k]:.4f}" for a, arm in enumerate(ARMS) for g in range(2)
                   if not rec["inband"][a, g, k]]
            rec["reason"] = f"kick {k + 1} at t = {t} outside its band: " + ", ".join(bad)
            return rec
        units = STEP - 1 if k < kicks - 1 else after               # the 5th unit is the next kick's "before"
        for unit in range(1, units + 1):
            for _ in range(per_unit):
                step(U, F, none)
            rec["R"][:, :, STEP * k + unit] = np.abs(ct.harmonic(U, F))
        if k < kicks - 1:
            for _ in range(per_unit):
                step(U, F, none)
    rec["carried"] = True
    return rec


def pulse(u0, f0, T, R_none, key, size=SIZE, after=AFTER):
    """The cause test's dense method at +-size from the untouched state at T (amplitude by ct.amplitude, the sham's
    paired signs drawn from `key`, group A before B): up, down and sham side by side, then untouched to T + after."""
    slow, stepper = nt.model(N, K, EPS)
    count = u0.shape[-1]
    U, F = np.repeat(u0[None], 3, 0), np.repeat(f0[None], 3, 0)
    z = ct.harmonic(u0[None], f0[None])[0]
    R, psi = np.abs(z), np.angle(z)
    along = np.sin(psi[:, None] - np.angle(u0))
    draw = np.random.default_rng(key)
    band = bands(size)
    nan = np.nan
    rec = dict(R=np.full((3, 2, after + 1), nan), amp=np.full(2, nan), status=np.zeros(2, int),
               ratio=np.full((3, 2), nan), inband=np.zeros((3, 2), bool), sumsq=np.full(3, nan),
               largest=np.full(3, nan), carried=False, reason="")
    for g in range(2):
        rec["amp"][g], rec["status"][g] = ct.amplitude(u0[g], along[g], z[g], size * R[g])
    if rec["status"].any():
        rec["reason"] = "pulse " + ", ".join(f"{'AB'[g]} {'infeasible' if rec['status'][g] == 1 else 'unsolved'}"
                                             for g in range(2) if rec["status"][g])
        return rec
    dtheta = np.zeros((3, 2, count))
    for g in range(2):
        signs = ct.paired_signs(along[g], draw)
        move = rec["amp"][g] * along[g]
        dtheta[0, g], dtheta[1, g], dtheta[2, g] = move, -move, signs * move
    U *= np.exp(1j * dtheta)
    now = np.abs(ct.harmonic(U, F))
    rec["ratio"] = now / R
    for a, arm in enumerate(ARMS):
        rec["inband"][a] = (band[arm][0] <= rec["ratio"][a]) & (rec["ratio"][a] <= band[arm][1])
    rec["sumsq"], rec["largest"] = np.square(dtheta).sum((-2, -1)), np.abs(dtheta).max((-2, -1))
    rec["R"][:, :, 0] = now
    if not rec["inband"].all():
        rec["reason"] = "pulse outside its band: " + ", ".join(
            f"{arm} {'AB'[g]} at {rec['ratio'][a, g]:.4f}" for a, arm in enumerate(ARMS) for g in range(2)
            if not rec["inband"][a, g])
        return rec
    step, none = stepper(DT), np.zeros(3)
    per_unit = int(round(1.0 / DT))
    for unit in range(1, after + 1):
        for _ in range(per_unit):
            step(U, F, none)
        rec["R"][:, :, unit] = np.abs(ct.harmonic(U, F))
    rec["carried"] = True
    return rec


# ---------------------------------------------------------------------------------------------------------------------
# the pieces

def base_times(s):
    return sorted({t for ss, m, T in moment_table() if ss == s for t in (T - 500, T - 100, T)})


def base_work(task, deadline):
    """The base run of one start (ct.base_run step for step, with no fresh copies), in pieces if it must, keeping the
    untouched first copy's state at T - 500, T - 100 and T of every moment."""
    s = task["seed"]
    times = base_times(s)
    state, part = PARTS / f"{task['name']}.state.npz", PARTS / f"{task['name']}.npz"
    out = calibrate(s, deadline, state, moments=np.zeros(0, int), keep=times)
    if out is None:
        return dict(status="saved")
    if len(out["kept_u"]) != len(times):
        raise RuntimeError("a kept state is missing")
    same = reproduced(out, s)
    save_npz(part, start=s + 1, times=np.array(times), z=out["z"], same=np.array([same["coherence"], same["gap"],
             same["settling"]]), us=out["kept_u"], uf=out["kept_f"], spent=out["spent"], pieces=out["pieces"])
    state.unlink(missing_ok=True)
    return dict(status="done", sha256=sha256(part), error="" if all(same.values()) else f"reproduced {same}")


def moment_work(task):
    s, m, T = task["seed"], task["moment"], task["T"]
    with np.load(PARTS / f"base{s + 1:02d}.npz") as base:
        times = [int(t) for t in base["times"]]
        R_none, us, uf = np.abs(base["z"]), base["us"], base["uf"]
    out = {}
    H = 500 if task["kind"] == "hold500" else 100
    j = times.index(T - H)
    rec = hold(us[j].copy(), uf[j].copy(), T, H, R_none, (SEED, 3, s, m, H))
    out.update({f"h{H}_{key}": value for key, value in rec.items()})
    if task["kind"] == "rest":
        j = times.index(T)
        rec = pulse(us[j].copy(), uf[j].copy(), T, R_none, (SEED, 4, s, m))
        out.update({f"pulse_{key}": value for key, value in rec.items()})
    part = PARTS / f"{task['name']}.npz"
    save_npz(part, start=s + 1, moment=m, T=T, **out)
    return dict(status="done", sha256=sha256(part))


def work(task, deadline):
    return base_work(task, deadline) if task["kind"] == "base" else moment_work(task)


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


def base_tasks():
    return [dict(name=f"base{s + 1:02d}", kind="base", seed=s) for s in STARTS]


def moment_tasks():
    tasks = []
    for s, m, T in moment_table():
        tasks += [dict(name=f"h500_{s + 1:02d}_{m}", kind="hold500", seed=s, moment=m, T=T),
                  dict(name=f"rest_{s + 1:02d}_{m}", kind="rest", seed=s, moment=m, T=T)]
    return tasks


def run():
    PARTS.mkdir(parents=True, exist_ok=True)
    todo = [t for t in base_tasks() if not (PARTS / f"{t['name']}.npz").exists()]
    if not todo:
        todo = [t for t in moment_tasks() if not (PARTS / f"{t['name']}.npz").exists()]
    if not todo:
        print("Every task is done: run report.")
        return
    run_piece(RUNLOG, todo, job, BUDGET, START_BY)


# ---------------------------------------------------------------------------------------------------------------------
# the numbers of one hold or pulse

def hold_numbers(rec, R_none, T, H, size=SIZE, after=AFTER):
    c = rec["R"].mean(1)                                           # (arm, time from T - H)
    none = R_none[:, T - H:T + after + 1]
    c0 = none.mean(0)
    M, M0 = c[:, :H].mean(-1), c0[:H].mean()
    A = np.array([M[0] - M0, M0 - M[1], M[2] - M0]) / (size * M0)
    blocks, b0 = c[:, :H].reshape(3, -1, BLOCK).mean(-1), c0[:H].reshape(-1, BLOCK).mean(-1)
    block_A = np.stack([blocks[0] - b0, b0 - blocks[1], blocks[2] - b0]) / (size * b0)
    kick_t = T - H + STEP * np.arange(rec["kicks"])
    pre = rec["before"][:, :, 1:] / R_none[None, :, kick_t[1:]]   # just before each kick after the first, (arm, group, kick)
    goal = np.array([1 + size, 1 - size, 1.0])[:, None, None]
    ratio = rec["R"][:, :, :H] / none[None, :, :H]
    sumsq = rec["sumsq"]
    left = {}
    for name, (a, b) in WINDOWS.items():
        if b > after:                                              # (the dry run's short follow-up)
            continue
        span = slice(H + a + 1, H + b + 1)
        level = c0[span].mean()
        d = (c[:, span].mean(-1) - level) / (size * level)
        left[name] = d * np.array([1.0, -1.0, 1.0])                # down: positive while it stays below untouched
    return dict(M=M, M0=M0, A=A, block_A=block_A, dM=(M[0] - M[1]) / 2, integral=H * (M - M0),
                drift=(pre - goal).mean((1, 2)), spread=(pre - goal).std((1, 2)),
                p5=np.percentile(ratio, 5, axis=(1, 2)), p95=np.percentile(ratio, 95, axis=(1, 2)),
                sumsq=sumsq.sum(-1), first_share=sumsq[:, 0] / sumsq.sum(-1),
                maintenance=np.median(sumsq[:, 1:], -1), largest=rec["largest"].max(-1),
                sum_abs_a=np.abs(rec["amp"]).sum(-1), net_a=rec["amp"][:2].sum(-1), left=left)


def pulse_numbers(rec, R_none, T, size=SIZE, after=AFTER):
    c = rec["R"].mean(1)
    c0 = R_none[:, T:T + after + 1].mean(0)
    left = {}
    for name, (a, b) in WINDOWS.items():
        span = slice(a + 1, b + 1)
        level = c0[span].mean()
        left[name] = (c[:, span].mean(-1) - level) / (size * level) * np.array([1.0, -1.0, 1.0])
    return dict(sumsq=rec["sumsq"], ratio=rec["ratio"], left=left)


# ---------------------------------------------------------------------------------------------------------------------
# the dry run

def dryrun():
    """Old start 1 for 700 units, one moment at T = 600 (holds from t = 100 and 500), 50 units after release: the
    base run as the pieces run it (stopped every 5 seconds and gone on from the kept state), against ct.base_run and
    the saved run (exactly); a 100-unit hold with no change asked for (size 0), which must follow the untouched run
    (the timing of the targets and the records); then the 500- and 100-unit holds and the pulse at +-25%."""
    if DRY.exists():
        shutil.rmtree(DRY)
    DRY.mkdir(parents=True)
    s, run, T, after = 0, 700.0, 600, 50
    times = [T - 500, T - 100, T]
    clock = time.time()
    res = ct.base_run(s, keep=times, run=run)
    base_seconds = time.time() - clock
    stops, out = 0, None
    while out is None:
        out = calibrate(s, time.time() + 5.0, DRY / "base.state.npz", moments=np.zeros(0, int), run=run, keep=times)
        stops += out is None
    same = reproduced(out, s, run)
    kept = bool(len(out["kept_u"]) == len(times) and all(
        np.array_equal(out["kept_u"][j], res["states"][j][0][0]) and np.array_equal(out["kept_f"][j], res["states"][j][1][0])
        for j in range(len(times))))
    agree = bool(np.array_equal(out["z"], res["z"][0]) and np.array_equal(out["logs"], res["logs"][0]))
    R_none = np.abs(out["z"])
    state = {t: (out["kept_u"][j].copy(), out["kept_f"][j].copy()) for j, t in enumerate(times)}
    lines = [f"Dry run, {shown(stamp())}: old start 1 (seed 0), 700 units, one moment at T = 600, 50 units after",
             f"  base run in pieces: stopped {stops} times, {out['spent']:.1f} s in all (ct.base_run in one go "
             f"{base_seconds:.1f} s); reproduces the saved run exactly: {same}; the same coherence, gap and kept "
             f"states as ct.base_run to the last digit: {agree and kept}"]
    zero = hold(*state[T - 100], T, 100, R_none, (SEED, 99, 0), size=0.0, after=after)
    follow = float(np.abs(zero["R"] - R_none[None, :, T - 100:T + after + 1]).max())
    lines.append(f"  size 0, 100-unit hold: carried {zero['carried']}; largest |R - R_none| over the hold and after "
                 f"{follow:.3g}; largest kick {np.nanmax(zero['largest']):.3g}")
    ok = all(same.values()) and agree and kept and zero["carried"] and follow <= 1e-9
    for H in HOLDS:
        clock = time.time()
        rec = hold(*state[T - H], T, H, R_none, (SEED, 99, H), after=after)
        seconds = time.time() - clock
        # the hold's rows sum their free turners in a plain layout, the base run in its interleaved one: the same
        # numbers to rounding, not to the last digit
        first = float(np.abs(rec["before"][:, :, 0] - R_none[None, :, T - H]).max())
        line = (f"  {H}-unit hold at +-25%: {seconds:.1f} s; carried {rec['carried']} {rec['reason']}; all three rows "
                f"start on the untouched run (largest difference {first:.2g}, needs at most 1e-15)")
        ok = ok and rec["carried"] and first <= 1e-15
        if rec["carried"]:
            x = hold_numbers(rec, R_none, T, H, after=after)
            line += (f"\n    achieved fraction up {x['A'][0]:.3f}, down {x['A'][1]:.3f}, sham {x['A'][2]:+.3f}; "
                     f"blocks up {x['block_A'][0].min():.2f} to {x['block_A'][0].max():.2f}, down "
                     f"{x['block_A'][1].min():.2f} to {x['block_A'][1].max():.2f}"
                     f"\n    ratio just before a kick, drift from target: up {x['drift'][0]:+.4f}, down "
                     f"{x['drift'][1]:+.4f}, sham {x['drift'][2]:+.4f}; spread up {x['spread'][0]:.4f}, down "
                     f"{x['spread'][1]:.4f}"
                     f"\n    sum of squared moves: up {x['sumsq'][0]:.2f}, down {x['sumsq'][1]:.2f}, sham "
                     f"{x['sumsq'][2]:.2f}; first kick's share {np.round(x['first_share'], 3).tolist()}; largest move "
                     f"{np.round(x['largest'], 3).tolist()}")
        lines.append(line)
    clock = time.time()
    rec = pulse(*state[T], T, R_none, (SEED, 99, 1), after=after)
    seconds = time.time() - clock
    lines.append(f"  pulse at +-25%: {seconds:.1f} s; carried {rec['carried']} {rec['reason']}; ratios "
                 f"{np.round(rec['ratio'], 4).tolist()}; sum of squared moves {np.round(rec['sumsq'], 2).tolist()}")
    ok = ok and rec["carried"]
    lines.append(f"  the cause test's locked files unchanged: {not lock_changes()}")
    lines.append(f"  dry run {'passed' if ok else 'FAILED'}")
    (DRY / "dryrun.txt").write_text("\n".join(lines) + "\n")
    print("\n".join(lines))


# ---------------------------------------------------------------------------------------------------------------------
# the report

def load():
    table = moment_table()
    base, recs, sums = {}, {}, {}
    for s in STARTS:
        path = PARTS / f"base{s + 1:02d}.npz"
        sums[path.stem] = sha256(path)
        with np.load(path) as b:
            base[s] = dict(R=np.abs(b["z"]), same=b["same"].copy(), seconds=float(b["spent"]))
    for s, m, T in table:
        for name, prefixes in ((f"h500_{s + 1:02d}_{m}", ("h500_",)), (f"rest_{s + 1:02d}_{m}", ("h100_", "pulse_"))):
            path = PARTS / f"{name}.npz"
            sums[name] = sha256(path)
            with np.load(path) as p:
                assert int(p["start"]) == s + 1 and int(p["moment"]) == m and int(p["T"]) == T
                for prefix in prefixes:
                    what = "pulse" if prefix == "pulse_" else int(prefix[1:-1])
                    recs[(s, m, what)] = {key[len(prefix):]: p[key][()] if p[key].ndim == 0 else p[key].copy()
                                          for key in p.files if key.startswith(prefix)}
    return table, base, recs, sums


def natural_spreads(base):
    """The pooled within-start SD of the H-unit history mean at T = 2,000, 2,005, ..., 5,500 in these six starts."""
    out = {}
    T = np.arange(2000, 5501, STEP)
    for H in HOLDS:
        rows = []
        for s in STARTS:
            c = base[s]["R"].mean(0)
            cum = np.concatenate([[0.0], np.cumsum(c)])
            rows.append((cum[T] - cum[T - H]) / H)
        rows = np.array(rows)
        d = rows - rows.mean(1, keepdims=True)
        out[H] = dict(sd=float(np.sqrt((d * d).sum() / (d.shape[0] * (d.shape[1] - 1)))), level=float(rows.mean()))
    return out


def summarise(values):
    values = np.asarray(values, float)
    return dict(mean=float(values.mean()), sd=float(values.std(ddof=1)) if values.size > 1 else 0.0,
                low=float(values.min()), high=float(values.max()))


def analyse(table, base, recs, v1_revised):
    natural = natural_spreads(base)
    out = dict(holds={}, pulse={}, natural={str(H): v for H, v in natural.items()},
               base_reproduced=bool(all(base[s]["same"].all() for s in STARTS)),             # as preregistered, kept
               base_reproduced_count=int(sum(bool(base[s]["same"].all()) for s in STARTS)),
               base_revised=v1_revised, reading=bool(v1_revised["met"]))                     # amendment 1
    values = {}
    for H in HOLDS:
        carried = np.array([[bool(recs[(s, m, H)]["carried"]) for m in range(MOMENTS)] for s in STARTS])
        reasons = [f"start {s + 1}, moment {m}: {recs[(s, m, H)]['reason']}" for s, m, T in table
                   if not recs[(s, m, H)]["carried"]]
        rows = [(s, m, T) for s, m, T in table if recs[(s, m, H)]["carried"]]
        nums = [hold_numbers(recs[(s, m, H)], base[s]["R"], T, H) for s, m, T in rows]
        missing, fewest = int(carried.size - carried.sum()), int(carried.sum(1).min())
        F1 = missing <= CRITERIA["most_missing"] and fewest >= CRITERIA["fewest_kept"]
        res = dict(carried=int(carried.sum()), missing=missing, fewest_per_start=fewest, reasons=reasons, F1=bool(F1))
        if nums:
            A = np.array([x["A"] for x in nums])                   # (moment, arm)
            lo, hi = CRITERIA["achieved"]
            inside = (A[:, 0] >= lo) & (A[:, 0] <= hi) & (A[:, 1] >= lo) & (A[:, 1] <= hi)
            mlo, mhi = CRITERIA["achieved_mean"]
            F2 = bool(inside.mean() >= CRITERIA["share"] and mlo <= A[:, 0].mean() <= mhi
                      and mlo <= A[:, 1].mean() <= mhi)
            blo, bhi = CRITERIA["block"]
            held = np.array([bool(((x["block_A"][:2] >= blo) & (x["block_A"][:2] <= bhi)).all()) for x in nums])
            F3 = bool(held.mean() >= CRITERIA["share"])
            sumsq = np.array([x["sumsq"] for x in nums])           # (moment, arm)
            sham_effort = sumsq[:, 2] / sumsq[:, :2].mean(1)
            elo, ehi = CRITERIA["sham_effort"]
            hlo, hhi = CRITERIA["sham_history"]
            F4a = bool(((sham_effort >= elo) & (sham_effort <= ehi)).all())
            F4b = bool(hlo <= A[:, 2].mean() <= hhi)
            pulses = {(s, m): pulse_numbers(recs[(s, m, "pulse")], base[s]["R"], T) for s, m, T in table
                      if recs[(s, m, "pulse")]["carried"]}
            versus = np.array([x["sumsq"][:2] / pulses[(s, m)]["sumsq"][:2]
                               for (s, m, T), x in zip(rows, nums) if (s, m) in pulses])
            planned = SIZE * natural[H]["level"]
            res.update(
                F2=F2, F3=F3, F4a=F4a, F4b=F4b, feasible=bool(F1 and F2 and F3),
                achieved={arm: summarise(A[:, a]) for a, arm in enumerate(ARMS)}, inside_share=float(inside.mean()),
                held_share=float(held.mean()),
                block_range={arm: [float(min(x["block_A"][a].min() for x in nums)),
                                   float(max(x["block_A"][a].max() for x in nums))] for a, arm in enumerate(ARMS[:2])},
                dM=summarise([x["dM"] for x in nums]),
                dM_in_natural_sd=float(np.mean([x["dM"] for x in nums]) / natural[H]["sd"]),
                planned_dM=planned, planned_in_natural_sd=planned / natural[H]["sd"],
                integral={arm: summarise([x["integral"][a] for x in nums]) for a, arm in enumerate(ARMS)},
                drift={arm: summarise([x["drift"][a] for x in nums]) for a, arm in enumerate(ARMS)},
                spread={arm: summarise([x["spread"][a] for x in nums]) for a, arm in enumerate(ARMS)},
                ratio_5_95={arm: [float(np.mean([x["p5"][a] for x in nums])), float(np.mean([x["p95"][a] for x in nums]))]
                            for a, arm in enumerate(ARMS)},
                kicks=H // STEP + 1,
                sumsq={arm: summarise(sumsq[:, a]) for a, arm in enumerate(ARMS)},
                sham_effort=summarise(sham_effort),
                versus_pulse={arm: summarise(versus[:, a]) for a, arm in enumerate(ARMS[:2])} if len(versus) else {},
                first_share={arm: summarise([x["first_share"][a] for x in nums]) for a, arm in enumerate(ARMS)},
                maintenance={arm: summarise([x["maintenance"][a] for x in nums]) for a, arm in enumerate(ARMS)},
                largest={arm: float(max(x["largest"][a] for x in nums)) for a, arm in enumerate(ARMS)},
                sum_abs_a={arm: summarise([x["sum_abs_a"][a].sum() for x in nums]) for a, arm in enumerate(ARMS)},
                net_a={arm: summarise([x["net_a"][a].sum() for x in nums]) for a, arm in enumerate(ARMS[:2])},
                left={name: {arm: summarise([x["left"][name][a] for x in nums]) for a, arm in enumerate(ARMS)}
                      for name in WINDOWS})
            values[f"A_{H}"] = A
            values[f"dM_{H}"] = np.array([x["dM"] for x in nums])
            values[f"sumsq_{H}"] = sumsq
        else:
            res.update(F2=False, F3=False, F4a=False, F4b=False, feasible=False)
        values[f"carried_{H}"] = carried
        out["holds"][str(H)] = res
    pulse_carried = np.array([[bool(recs[(s, m, "pulse")]["carried"]) for m in range(MOMENTS)] for s in STARTS])
    pn = [pulse_numbers(recs[(s, m, "pulse")], base[s]["R"], T) for s, m, T in table
          if recs[(s, m, "pulse")]["carried"]]
    out["pulse"] = dict(
        carried=int(pulse_carried.sum()),
        reasons=[f"start {s + 1}, moment {m}: {recs[(s, m, 'pulse')]['reason']}" for s, m, T in table
                 if not recs[(s, m, "pulse")]["carried"]],
        sumsq={arm: summarise([x["sumsq"][a] for x in pn]) for a, arm in enumerate(ARMS)} if pn else {},
        left={name: {arm: summarise([x["left"][name][a] for x in pn]) for a, arm in enumerate(ARMS)}
              for name in WINDOWS} if pn else {})
    values["carried_pulse"] = pulse_carried
    values["moments"] = np.array(table)
    return out, values


def report_lines(res, log, sums):
    yes = {True: "met", False: "NOT met"}
    v = res["base_revised"]
    lines = ["Hold feasibility (prelock_checks_plan.md, part 2): old starts 13 to 18, ten moments each, +-25%",
             f"made {shown(res['made'])}; not a locked plan; no new start; no gap and no growth measured", "",
             "Reading: " + ("given below, with the revised exactness check (adopted after the preregistered one "
                            "failed in the calibration; prelock_amendment.md)" if res["reading"] else
                            "NO READING: the revised exactness check is not met; the numbers below are reported as "
                            "they stand"),
             f"Base runs, as preregistered (each reproduces its saved run exactly): {yes[res['base_reproduced']]} "
             f"({res['base_reproduced_count']} of 6)",
             f"   first difference from the saved run, by start (gap; coherence): "
             f"{saved_difference_text(v['first_saved_difference'])}",
             f"Base runs, revised (not the preregistered check): each agrees to the last digit with the cause test's "
             f"unchanged base-run code run in one go on the same machine (coherence at full precision and every kept "
             f"state), in one numerical environment throughout: {yes[v['met']]} ({sum(v['agree'].values())} of 6 "
             f"agree; one environment: {v['environment_same']})",
             f"   the environment: {environment_text(v)}",
             "   bitwise agreement on one machine shows reproducibility, not numerical accuracy (prelock_amendment.md)",
             "Seen before this check (prelock_amendment.md; dryruns/feasibility_dryrun.txt): the dry run's holds on old "
             "start 1 at T = 600, outside these starts and moments: 500-unit fractions up 0.976, down 0.962, sham "
             "+0.218, one up block at 0.15; 100-unit fractions 0.927 and 0.907, sham +0.031."]
    for H in HOLDS:
        h = res["holds"][str(H)]
        lines += ["", f"{H}-unit hold ({H // STEP + 1} kicks): {'FEASIBLE' if h['feasible'] else 'NOT FEASIBLE'}",
                  f"  F1 kicks on target: {h['carried']} of 60 moments carried out, fewest in a start "
                  f"{h['fewest_per_start']} (needs at most 3 left out, at least 5 per start): {yes[h['F1']]}"]
        for reason in h["reasons"][:6]:
            lines.append(f"     left out: {reason}")
        if "achieved" not in h:
            continue
        a = h["achieved"]
        lines += [f"  F2 history achieved: both fractions in 0.75-1.25 in {100 * h['inside_share']:.0f}% of moments "
                  f"(needs 90%); mean up {a['up']['mean']:.3f} (SD {a['up']['sd']:.3f}), down {a['down']['mean']:.3f} "
                  f"(SD {a['down']['sd']:.3f}) (each needs 0.85-1.15): {yes[h['F2']]}",
                  f"  F3 held throughout: every 25-unit block in 0.5-1.5 in {100 * h['held_share']:.0f}% of moments "
                  f"(needs 90%); blocks ranged up {h['block_range']['up'][0]:.2f} to {h['block_range']['up'][1]:.2f}, "
                  f"down {h['block_range']['down'][0]:.2f} to {h['block_range']['down'][1]:.2f}: {yes[h['F3']]}",
                  f"  F4 sham (reported apart): effort against up and down {h['sham_effort']['low']:.3f} to "
                  f"{h['sham_effort']['high']:.3f} (needs 0.8-1.2): {yes[h['F4a']]}; mean history fraction "
                  f"{a['sham']['mean']:+.3f} (needs within +-0.25): {yes[h['F4b']]}",
                  "  Diagnostics:",
                  f"    history change dM = {h['dM']['mean']:.5f} (SD {h['dM']['sd']:.5f}) = "
                  f"{h['dM_in_natural_sd']:.2f} natural SDs of the {H}-unit mean (planned {h['planned_dM']:.5f}, "
                  f"{h['planned_in_natural_sd']:.2f} SDs)",
                  f"    time-integrated history H x (M - M_none): up {h['integral']['up']['mean']:+.3f}, down "
                  f"{h['integral']['down']['mean']:+.3f}, sham {h['integral']['sham']['mean']:+.4f}",
                  f"    just before each kick, ratio to untouched minus target: up {h['drift']['up']['mean']:+.4f} "
                  f"(spread {h['spread']['up']['mean']:.4f}), down {h['drift']['down']['mean']:+.4f} (spread "
                  f"{h['spread']['down']['mean']:.4f}), sham {h['drift']['sham']['mean']:+.4f} (spread "
                  f"{h['spread']['sham']['mean']:.4f})",
                  f"    ratio to untouched over the hold, 5th to 95th percentile: up {h['ratio_5_95']['up'][0]:.3f} to "
                  f"{h['ratio_5_95']['up'][1]:.3f}, down {h['ratio_5_95']['down'][0]:.3f} to "
                  f"{h['ratio_5_95']['down'][1]:.3f}, sham {h['ratio_5_95']['sham'][0]:.3f} to "
                  f"{h['ratio_5_95']['sham'][1]:.3f}",
                  f"    effort, sum of squared moves over the hold: up {h['sumsq']['up']['mean']:.2f}, down "
                  f"{h['sumsq']['down']['mean']:.2f}, sham {h['sumsq']['sham']['mean']:.2f}; against one pulse: up "
                  f"{h['versus_pulse'].get('up', {}).get('mean', float('nan')):.2f}x, down "
                  f"{h['versus_pulse'].get('down', {}).get('mean', float('nan')):.2f}x",
                  f"    first kick's share: up {h['first_share']['up']['mean']:.3f}, down "
                  f"{h['first_share']['down']['mean']:.3f}; median maintenance kick: up "
                  f"{h['maintenance']['up']['mean']:.4f}, down {h['maintenance']['down']['mean']:.4f}; largest single "
                  f"move: up {h['largest']['up']:.3f}, down {h['largest']['down']:.3f}, sham {h['largest']['sham']:.3f}",
                  f"    sum of |a| (both groups): up {h['sum_abs_a']['up']['mean']:.3f}, down "
                  f"{h['sum_abs_a']['down']['mean']:.3f}; net signed push: up {h['net_a']['up']['mean']:+.3f}, down "
                  f"{h['net_a']['down']['mean']:+.3f}",
                  "    coherence left after release, per unit of the planned change (down's sign turned): "
                  + "; ".join(f"{name} " + ", ".join(f"{arm} {h['left'][name][arm]['mean']:+.3f}" for arm in ARMS)
                              for name in WINDOWS)]
    p = res["pulse"]
    lines += ["", f"Pulse at +-25%: {p['carried']} of 60 carried out"]
    for reason in p["reasons"][:6]:
        lines.append(f"     left out: {reason}")
    if p["sumsq"]:
        lines += [f"  sum of squared moves: {p['sumsq']['up']['mean']:.2f} (each arm the same)",
                  "  coherence left after release, per unit of the planned change: "
                  + "; ".join(f"{name} " + ", ".join(f"{arm} {p['left'][name][arm]['mean']:+.3f}" for arm in ARMS)
                              for name in WINDOWS)]
    lines += ["", f"Pieces: {len(log['pieces'])}; parts matching the run log: "
              f"{sum(log['parts'].get(name, {}).get('sha256') == value for name, value in sums.items())} of {len(sums)}"]
    return lines


def report():
    names = [t["name"] for t in base_tasks() + moment_tasks()]
    missing = [name for name in names if not (PARTS / f"{name}.npz").exists()]
    if missing:
        sys.exit(f"{len(missing)} of {len(names)} tasks are not done yet: run the pieces first.")
    log = json.loads(RUNLOG.read_text())
    changed = lock_changes()
    now = folder_fingerprints()
    moved = sorted(name for name in set(now) | set(log["folders_before"])
                   if now.get(name) != log["folders_before"].get(name))
    table, base, recs, sums = load()
    res, values = analyse(table, base, recs, revised_v1(log, "feasibility", STARTS))
    res.update(made=stamp(), criteria=CRITERIA, size=SIZE, holds_order=list(HOLDS),
               unchanged=dict(lock=not changed, folders=not moved, changed=changed + moved))
    lines = report_lines(res, log, sums)
    lines.append(f"The cause test's locked files unchanged: {not changed}; its folder and the exploratory folder "
                 f"unchanged: {not moved}")
    write_json(OUT_JSON, res)
    save_npz(OUT_NPZ, **values)
    OUT_TXT.write_text("\n".join(lines) + "\n")
    print("\n".join(lines))


if __name__ == "__main__":
    commands = dict(dryrun=dryrun, run=run, report=report)
    command = sys.argv[1] if len(sys.argv) > 1 else ""
    if command not in commands:
        sys.exit(__doc__)
    commands[command]()