#!/usr/bin/env python3
"""
Cause test, exploratory follow-up, round 2: coherence history (October 12, 2026, late evening). EXPLORATORY.

Approved by Colson after round 1 (explore_cause.py), as a read-only analysis. Chosen after the cause test's decision
(FAIL, as registered) was known, and run on the same saved runs, so nothing here is a test, whatever it shows. The
decision is not revisited: it stays FAIL. The aim is to describe how the untouched runs' link between coherence and
later growth depends on the length of the history, and to set up the next discriminating experiment.

Read-only. The model is not run (no new simulation; no base run regenerated). Every file read, every file the cause
test locked, the cause test's results and round 1's outputs are fingerprinted (SHA-256) before and after, and the two
lists must match; the 16 locked files must also still match cause_lock.json. Outputs go to this folder only. A sandbox
script kept with the project; nothing here is served by the site.

Reads
    ../cause-test-2026-10-12/sim_cause4000.npz    the 42 base runs (starts 19 to 60) and 420 moments
    ../cause-test-2026-10-12/cause_test.json      the report's Phi and share left
    ../cause-test-2026-10-12/design_numbers.json  the benchmark
    ../cause-test-2026-10-12/cause_lock.json      the locked fingerprints
    explore_cause.json                            round 1's numbers, to check that this code gives the same ones
    the 18 older saved starts (1 to 18) through design_numbers.saved_sets(), as a cross-check

The list of analyses below was written before any of them was run. Nothing was added after seeing results; anything
that had to change after a first run is listed in the report under "Changes after the first run".

Definitions (untouched base runs; every quantity used at T is known at T, and no control looks past T)
  c(t)       the coherence (R_A + R_B) / 2 at every time unit t = 0, 1, ..., 6,000 of the saved run
  T          2,000, 2,005, ..., 5,500: 701 per start, as in round 1 and the benchmark
  W          the gap's growth rate over (T + 200, T + 500]; E: over (T, T + 200]
  g1 to g5   the growth rate in the five 100-unit windows up to T, as in Test 17 and the benchmark
  M_tau      the mean of c over the tau time units before T: t = T - tau, ..., T - 1. T itself is left out, so a change
             made at T moves only F_tau. Every history window ends before T, and W begins at T + 200.
  F_tau      C(T) - M_tau
  tau        5, 10, 25, 50, 100, 150, 200, 300, 400, 500, 750, 1000, 1500, 2000 (fixed now; every one is reported)
  bands      the 14 disjoint pieces between neighbouring tau: the mean of c over t = T - tau_k, ..., T - tau_(k-1) - 1,
             with tau_0 = 0
  coarse     the mean of c over lags (0, 100], (100, 500], (500, 1000] and (1000, 2000] before T, for one joint fit
M_tau is a mathematical summary of the past, not a measured physical memory. Every slope and correlation in H1, H2 and
H4 is an association in untouched runs, not an effect.

What it works out (start by start unless said otherwise; new 42 starts, and the old 18 as a cross-check)
  H1. The window-response curve, for every tau: the correlation of M_tau with W (no controls); the partial correlation
      of M_tau with W given g; the slope b in W = a + b M_tau + c1 g1 + ... + c5 g5; the partial correlation of M_tau
      with W given g and C(T). For every disjoint band: its partial correlation with W given g, and given g and C(T).
      One joint fit: W on C(T), the four coarse bands and g; each one's partial correlation given all the others, and
      its slope. Check: round 1's "C(T) with W given g and the mean over (T - 100, T]" (-0.0191) is reproduced.
  H2. Slow and fast parts, for every tau: W = a + b_M M_tau + b_F F_tau + g terms (b_M, b_F); the partial correlations
      of M_tau with W given g and F_tau, and of F_tau with W given g and M_tau (the latter equals C(T)'s, given g and
      M_tau); E = a + a_M M_tau + a_F F_tau + g terms (a_M, a_F), and F_tau's partial correlation with E given g and
      M_tau; the spreads (SD within each start) of M_tau and F_tau; natural persistence: the mean c over (T, T + 200]
      and over (T + 200, T + 500] = a + p_M M_tau + p_F F_tau (p_F: how much of a natural fast deviation is still there).
  H3. The natural fast part against the interventions (kept apart: the first is an association, the second an effect
      of the specified pulses):
      a. the induced slopes, as in round 1: the growth difference between up and down over their coherence difference
         at T, start by start, delayed (T + 200, T + 500] and early (T, T + 200], per method and pooled. A change made
         at T moves F_tau by exactly that much and leaves M_tau as it was, for every tau;
      b. at the 420 moments, the untouched continuation: W = a + b_M M_tau + b_F F_tau + g terms and the same for E,
         each start's means taken out and pooled over the moments (error bars from redrawing the starts);
      c. per start, the induced slope (pooled over the methods) minus the natural b_F (H2, every T), delayed and early;
      d. persistence: the natural p_F beside the interventions' Phi (0.511) and share left (0.297), from the report;
      e. the size of the change (about 0.0130 at T, up against sham) beside the natural spread of F_tau.
  H4. Stability and possible artifacts:
      a. independent starts: H1's M_tau with W given g, and H2's two partial correlations, in five independent blocks
         (new starts 19 to 39 and 40 to 60; old 1 to 6, 7 to 12, 13 to 18), and how many of the 60 starts are below 0;
      b. overlapping rows: only T = 2,000, 2,500, ..., 5,500 (8 per start; no two rows share a growth window, and their
         history windows up to 500 units do not overlap), pooled within starts (each start's means taken out), beside
         the same pooled measure on all 701 rows; the same three partial correlations; error bars from redrawing;
      c. autocorrelation: mismatched starts. Each start's coherence (M_tau, F_tau, C) is paired with the growth (W and
         g) of the next start in the same set (start i with i + 1, the last with the first), at the same T. The same
         measures; the mean ± error bar, the SD over starts beside the matched SD, and the count below 0. Independent
         runs share no dynamics, so this shows what autocorrelation and chance alone give;
      d. slow drift: (i) with each start's own time trend (T and T squared) added to the controls; (ii) the first half
         (T up to 3,750) and second half of each run separately; (iii) the mean c and the mean growth rate in each
         1,000-unit block of the run, averaged over the starts;
      e. a longer growth history for tau > 500: M_tau with W given g and the growth over (T - tau, T - 500], and F_tau
         with W given g, M_tau and it.
  H5. Every window, band and control above is reported, including nulls and contradictions. No window is singled out;
      the strongest window is not used for anything, and no threshold is tuned here.
  H6. What the saved data cannot resolve, stated in the report.

Fixed now, before this analysis runs, for the proposal that is to follow it (held changes against pulses): hold lengths
of 100 and 500 time units (the two history windows named in round 1's plan, before either round's results); the outcome
window stays (T + 200, T + 500]; a fresh gap is put in at T + 200, the same in every arm of a moment; the untouched arm
is part of the main comparisons; new starts from 61 on. Nothing in this analysis is used to change these; anything it
suggests in addition is added to the proposal labelled as chosen after it.

Error bars: for per-start values, as in the cause test (cause_rule.error_bar: the larger of the resampled and SD/sqrt(n),
the starts redrawn 400 times with seed 20261015, a block's own size for a block); for pooled weights and pooled partial
correlations, the SD over the same 400 redraws of the starts.

Usage
    python3 -B explore_history.py      # explore_history.txt, explore_history.json and explore_history.png, here

Needs Python 3, numpy and matplotlib, and ../cause-test-2026-10-12 and ../coherence-growth-2026-10-11 with their saved
runs, unchanged.
"""
import hashlib
import json
import sys
from pathlib import Path

import numpy as np

HERE = Path(__file__).resolve().parent
CAUSE = HERE.parent / "cause-test-2026-10-12"
CG = HERE.parent / "coherence-growth-2026-10-11"
NUDGE = HERE.parent / "nudge-test-2026-10-10"
sys.dont_write_bytecode = True                             # leave every folder as it is
sys.path.insert(0, str(CAUSE))
sys.path.insert(1, str(CG))
import cause_rule as rule  # noqa: E402  (the decision rule's error bars; unchanged)
import design_numbers as dn  # noqa: E402  (the 18 older saved starts; unchanged)
import earlier_growth_check as egc  # noqa: E402  (the earlier growth and partial correlations, as in Test 17)

SAVED = CAUSE / "sim_cause4000.npz"
REPORT, DESIGN, LOCK = CAUSE / "cause_test.json", CAUSE / "design_numbers.json", CAUSE / "cause_lock.json"
ROUND1 = HERE / "explore_cause.json"
LOCKED = dict(                                             # cause_lock.json's names for the files it fixed
    plan=CAUSE / "cause_test_prereg.md", cause_test=CAUSE / "cause_test.py", cause_rule=CAUSE / "cause_rule.py",
    design_numbers=CAUSE / "design_numbers.py", design_text=CAUSE / "design_numbers.txt",
    design_json=DESIGN, moment_table=CAUSE / "cause_moments.csv", targets_text=CAUSE / "cause_targets.txt",
    targets_json=CAUSE / "cause_targets.json", nudge_test=NUDGE / "nudge_test.py",
    coherence_growth=CG / "coherence_growth.py", earlier_growth_check=CG / "earlier_growth_check.py",
    test16_results=NUDGE / "lyapunov_long.json", earlier_runs=NUDGE / "sim_lyapunov4000_long.npz",
    fresh_runs=CG / "sim_fresh4000.npz", replication_runs=CG / "sim_replication4000.npz")
RESULTS = [SAVED, REPORT, LOCK, CAUSE / "cause_test.txt", CAUSE / "cause_runlog.json",
           HERE.parent.parent / "screenshots" / "nudge-test-cause-test.png",
           *(HERE / f"explore_cause.{x}" for x in ("py", "txt", "json", "png"))]
EVERY = 5                                                  # time units per interval of the gap
FIRST, LAST = 2000, 5500                                   # base-run moments T
EARLY, LATER = (0, 200), (200, 500)                        # immediate and delayed windows after T
TAUS = (5, 10, 25, 50, 100, 150, 200, 300, 400, 500, 750, 1000, 1500, 2000)
EDGES = (0,) + TAUS
BANDS = tuple(zip(EDGES[:-1], EDGES[1:]))                  # lags (a, b] before T
COARSE = ((0, 100), (100, 500), (500, 1000), (1000, 2000))
LONG = tuple(tau for tau in TAUS if tau > 500)
SPLIT, THIN, BLOCK = 3750, 500, 1000                       # the time halves; non-overlapping rows; drift blocks
BLOCKS = (("new 19-39", "new", slice(0, 21)), ("new 40-60", "new", slice(21, 42)), ("old 1-6", "old", slice(0, 6)),
          ("old 7-12", "old", slice(6, 12)), ("old 13-18", "old", slice(12, 18)))
TRIO = ("r_M|g", "r_M|g,F", "r_F|g,M")                     # the three partial correlations checked throughout
OUT_TXT, OUT_JSON, FIGURE = (HERE / f"explore_history.{x}" for x in ("txt", "json", "png"))


def fingerprints(paths):
    return {str(p.resolve().relative_to(HERE.parent.parent)): hashlib.sha256(p.read_bytes()).hexdigest()
            for p in paths}


def fit(y, columns):
    """Least squares of y on a constant and the columns: the weights (constant first)."""
    design = np.column_stack([np.ones(len(y)), *columns])
    return np.linalg.lstsq(design, y, rcond=None)[0]


def summary(values, picks, keep=False):
    values = np.asarray(values, float)
    bar = float(rule.error_bar(values, picks))
    mean = float(values.mean())
    out = dict(mean=mean, error_bar=bar, bars=mean / bar if bar > 0 else float("nan"),
               below=int((values < 0).sum()), n=int(values.size), sd=float(values.std(ddof=1)))
    if keep:
        out["values"] = values.tolist()
    return out


def within(y, X, picks):
    """y: (start, row); X: (start, row, k). Each start's own means taken out, least squares pooled over the rows; the
    weights and their SD over the 400 redraws of the starts."""
    yd = y - y.mean(1, keepdims=True)
    Xd = X - X.mean(1, keepdims=True)

    def solve(rows):
        return np.linalg.lstsq(Xd[rows].reshape(-1, X.shape[-1]), yd[rows].reshape(-1), rcond=None)[0]

    return dict(weights=solve(np.arange(len(y))).tolist(), bars=np.array([solve(p) for p in picks]).std(0).tolist())


def products(x, y, Z):
    """Per start, the sums of products of [x, y, Z] after each start's own means are taken out.
    x, y: (start, row); Z: (start, row, k)."""
    A = np.concatenate([x[..., None], y[..., None], Z], -1)
    A = A - A.mean(1, keepdims=True)
    return np.einsum("srk,srl->skl", A, A)


def from_products(S):
    """The pooled partial correlation of x and y given Z, from the summed products."""
    held, bx, by = S[2:, 2:], S[2:, 0], S[2:, 1]
    sol = np.linalg.solve(held, np.column_stack([bx, by]))
    sxy, sxx, syy = S[0, 1] - bx @ sol[:, 1], S[0, 0] - bx @ sol[:, 0], S[1, 1] - by @ sol[:, 1]
    return float(sxy / np.sqrt(sxx * syy))


def pooled_partial(x, y, Z, picks):
    """Pooled within starts: each start's means taken out, then the partial correlation over all rows; the error bar is
    the SD over the 400 redraws of the starts (a start drawn twice counts twice)."""
    S = products(x, y, Z)
    counts = np.stack([np.bincount(p, minlength=len(S)) for p in picks])
    redrawn = np.array([from_products(np.tensordot(cnt, S, 1)) for cnt in counts])
    value = from_products(S.sum(0))
    bar = float(redrawn.std())
    return dict(mean=value, error_bar=bar, bars=value / bar if bar > 0 else float("nan"))


# ---------------------------------------------------------------------------------------------------------------------
# the untouched runs

def history_start(z, logs):
    """One start's series at T = 2,000, 2,005, ..., 5,500. z: (group, time), every time unit from t = 0; logs:
    (interval,), interval j covering (5j, 5j + 5]."""
    c = np.abs(z.astype(np.complex128)).mean(0)
    T = np.arange(FIRST, LAST + 1, EVERY)
    js = T // EVERY - 1                                    # the interval ending at T
    total = np.concatenate([[0.0], np.cumsum(logs)])

    def growth(a, b):                                      # the growth rate over (T + a, T + b]
        return (total[js + 1 + b // EVERY] - total[js + 1 + a // EVERY]) / (b - a)

    cum = np.concatenate([[0.0], np.cumsum(c)])

    def back(a, b):                                        # the mean c over t = T - b, ..., T - a - 1
        return (cum[T - a] - cum[T - b]) / (b - a)

    def ahead(a, b):                                       # the mean c over (T + a, T + b]
        return (cum[T + b + 1] - cum[T + a + 1]) / (b - a)

    C = c[T]
    M = np.stack([back(0, tau) for tau in TAUS])
    return dict(T=T, C=C, M=M, F=C - M, bands=np.stack([back(a, b) for a, b in BANDS]),
                coarse=np.stack([back(a, b) for a, b in COARSE]), W=growth(*LATER), E=growth(*EARLY),
                g=egc.before(logs, js, egc.MAIN), G_long=np.stack([growth(-tau, -500) for tau in LONG]),
                early_C=ahead(*EARLY), later_C=ahead(*LATER), past100_incl=ahead(-100, 0), c=c, logs=logs)


def trend(T):
    s = (T - 3750.0) / 1750.0                              # each start's own time, scaled to about -1 to 1
    return np.column_stack([s, s * s])


def numbers(f):
    """One start's numbers for H1, H2, H4d and H4e."""
    W, E, g, C = f["W"], f["E"], f["g"], f["C"]
    tt = trend(f["T"])
    halves = dict(first=f["T"] <= SPLIT, second=f["T"] > SPLIT)
    out = {}
    for k, tau in enumerate(TAUS):
        M, F = f["M"][k], f["F"][k]
        gM, gF = np.column_stack([g, M]), np.column_stack([g, F])
        out[f"r_plain:{tau}"] = egc.corr(M, W)
        out[f"r_M|g:{tau}"] = egc.partial(M, W, g)
        out[f"b_M|g:{tau}"] = fit(W, [M, *g.T])[1]
        out[f"r_M|g,C:{tau}"] = egc.partial(M, W, np.column_stack([g, C]))
        out[f"r_M|g,F:{tau}"] = egc.partial(M, W, gF)
        out[f"r_F|g,M:{tau}"] = egc.partial(F, W, gM)
        w = fit(W, [M, F, *g.T])
        out[f"b_M:{tau}"], out[f"b_F:{tau}"] = w[1], w[2]
        w = fit(E, [M, F, *g.T])
        out[f"a_M:{tau}"], out[f"a_F:{tau}"] = w[1], w[2]
        out[f"rE_F|g,M:{tau}"] = egc.partial(F, E, gM)
        for name in ("early_C", "later_C"):
            w = fit(f[name], [M, F])
            out[f"p_M {name}:{tau}"], out[f"p_F {name}:{tau}"] = w[1], w[2]
        out[f"sd_M:{tau}"], out[f"sd_F:{tau}"], out[f"mean_M:{tau}"] = M.std(), F.std(), M.mean()
        out[f"r_M|g trend:{tau}"] = egc.partial(M, W, np.column_stack([g, tt]))
        out[f"r_M|g,F trend:{tau}"] = egc.partial(M, W, np.column_stack([gF, tt]))
        out[f"r_F|g,M trend:{tau}"] = egc.partial(F, W, np.column_stack([gM, tt]))
        for half, rows in halves.items():
            out[f"r_M|g {half}:{tau}"] = egc.partial(M[rows], W[rows], g[rows])
            out[f"r_M|g,F {half}:{tau}"] = egc.partial(M[rows], W[rows], gF[rows])
            out[f"r_F|g,M {half}:{tau}"] = egc.partial(F[rows], W[rows], gM[rows])
        if tau in LONG:
            G = f["G_long"][LONG.index(tau)]
            out[f"r_M|g,G:{tau}"] = egc.partial(M, W, np.column_stack([g, G]))
            out[f"r_F|g,M,G:{tau}"] = egc.partial(F, W, np.column_stack([gM, G]))
        band = f["bands"][k]
        out[f"band r|g:{tau}"] = egc.partial(band, W, g)
        out[f"band r|g,C:{tau}"] = egc.partial(band, W, np.column_stack([g, C]))
    X = [C, *f["coarse"]]
    w = fit(W, [*X, *g.T])
    for j, name in enumerate(joint_names()):
        others = [x for i, x in enumerate(X) if i != j]
        out[f"joint r:{name}"] = egc.partial(X[j], W, np.column_stack([*others, g]))
        out[f"joint b:{name}"] = w[1 + j]
    out["link"] = egc.partial(C, W, g)
    out["beta_held"] = fit(W, [C, *g.T])[1]
    out["sd_C"], out["mean_C"] = C.std(), C.mean()
    out["check past100"] = egc.partial(C, W, np.column_stack([g, f["past100_incl"]]))
    for b in range(6):                                     # H4d(iii): the run in 1,000-unit blocks
        out[f"drift c:{b}"] = f["c"][BLOCK * b + 1:BLOCK * (b + 1) + 1].mean()     # t in (1000b, 1000(b + 1)]
        span = BLOCK // EVERY
        out[f"drift growth:{b}"] = f["logs"][span * b:span * (b + 1)].sum() / BLOCK
    return out


def joint_names():
    return ["C(T)"] + [f"lags ({a}, {b}]" for a, b in COARSE]


def pick(table, tau):
    """A per-tau entry, whether the keys are numbers (in memory) or text (read back from JSON)."""
    return table[str(tau)] if str(tau) in table else table[tau]


def mismatched(fc, fw):
    """Coherence from one start (fc), growth from another (fw), at the same T."""
    out = {}
    for k, tau in enumerate(TAUS):
        M, F = fc["M"][k], fc["F"][k]
        out[f"r_M|g:{tau}"] = egc.partial(M, fw["W"], fw["g"])
        out[f"r_M|g,F:{tau}"] = egc.partial(M, fw["W"], np.column_stack([fw["g"], F]))
        out[f"r_F|g,M:{tau}"] = egc.partial(F, fw["W"], np.column_stack([fw["g"], M]))
    out["link"] = egc.partial(fc["C"], fw["W"], fw["g"])
    return out


def pooled_trio(series, picks, rows):
    """H4b: the three partial correlations pooled within starts, for every tau, on the given rows of every start."""
    W = np.stack([f["W"][rows] for f in series])
    g = np.stack([f["g"][rows] for f in series])
    out = {}
    for k, tau in enumerate(TAUS):
        M = np.stack([f["M"][k][rows] for f in series])
        F = np.stack([f["F"][k][rows] for f in series])
        out[tau] = {"r_M|g": pooled_partial(M, W, g, picks),
                    "r_M|g,F": pooled_partial(M, W, np.concatenate([g, F[..., None]], -1), picks),
                    "r_F|g,M": pooled_partial(F, W, np.concatenate([g, M[..., None]], -1), picks)}
    return out


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

def main():
    lock = json.loads(LOCK.read_text())["sha256"]
    inputs = [*LOCKED.values(), *RESULTS, DESIGN]
    inputs = list(dict.fromkeys(p.resolve() for p in inputs))
    before = fingerprints(inputs)
    lock_ok = all(hashlib.sha256(path.read_bytes()).hexdigest() == lock[name] for name, path in LOCKED.items())
    assert lock_ok and set(lock) == set(LOCKED), "a locked file no longer matches cause_lock.json"
    d = np.load(SAVED)
    report = json.loads(REPORT.read_text())
    round1 = json.loads(ROUND1.read_text())
    bench = float(json.loads(DESIGN.read_text())["benchmark"]["beta"])
    starts = [int(s) for s in d["starts"]]
    n = len(starts)
    picks = rule.draws(n)
    assert np.asarray(d["carried"], bool).all(), "this script assumes every moment was carried out (all 420 were)"
    Tm = np.asarray(d["T"], int)
    per = Tm.shape[1]
    z_base, logs_base = np.asarray(d["z"]), np.asarray(d["logs"], float)
    res = dict(exploratory=("round 2, approved by Colson after round 1; chosen after the cause test's decision (FAIL) was"
                            " known and run on the same saved runs: not a test; the decision stays FAIL as registered"),
               starts=starts, taus=list(TAUS), bands=[list(b) for b in BANDS], changes_after_first_run=[])

    # the untouched runs, start by start
    series = dict(new=[history_start(z_base[i], logs_base[i]) for i in range(n)],
                  old=[history_start(z[j], logs[j]) for _, z, logs in dn.saved_sets() for j in range(len(z))])
    per_start = {key: [numbers(f) for f in fs] for key, fs in series.items()}
    draws = dict(new=picks, old=rule.draws(len(series["old"])))
    keep = {f"{m}:{tau}" for m in (*TRIO, "b_M", "b_F") for tau in TAUS}

    def both(key):
        return {name: summary([x[key] for x in per_start[name]], draws[name], key in keep) for name in ("new", "old")}

    res["per_start"] = {key: both(key) for key in per_start["new"][0]}
    S = res["per_start"]
    assert abs(S["beta_held"]["old"]["mean"] + bench) < 1e-9, "the benchmark is not reproduced"
    assert abs(S["check past100"]["new"]["mean"] - round1["A"]["past100:C|g,F"]["new"]["mean"]) < 1e-12, \
        "round 1's (T - 100, T] value is not reproduced"

    # H3a: the interventions (as in round 1)
    bl = np.asarray(d["branch_logs"], float)

    def rate(x, window):
        low, high = (v // EVERY for v in window)
        return x[..., low:high].sum(-1) / (window[1] - window[0])

    C1v = np.asarray(d["C1"], float).reshape(n, per, 2, 3)
    G = {name: rate(bl, w).reshape(n, per, 2, 3) for name, w in (("early", EARLY), ("delayed", LATER))}
    spread = (C1v[..., 0] - C1v[..., 2]) - (C1v[..., 1] - C1v[..., 2])      # (C_up - C_down) at T, against the sham
    induced = {name: (g_[..., 0] - g_[..., 1]).mean(1) / spread.mean(1) for name, g_ in G.items()}   # (start, method)
    res["induced"] = {name: dict(dense=summary(v[:, 0], picks), sparse=summary(v[:, 1], picks),
                                 pooled=summary(v.mean(-1), picks)) for name, v in induced.items()}
    res["change_at_T"] = summary((spread / 2).mean((1, 2)), picks)                # (C_up - C_down) / 2 at T
    for name, key in (("delayed", "C(T) -> delayed growth"), ("early", "C(T) -> early growth")):
        assert abs(res["induced"][name]["pooled"]["mean"] - round1["slopes"][key]["induced"]["pooled"]["mean"]) < 1e-12
    res["persistence_interventions"] = dict(phi=report["phi"]["pooled"], left_over=report["left_over"]["pooled"])

    # H3b: the 420 moments, untouched continuation
    c_base = np.abs(z_base.astype(np.complex128)).mean(1)
    assert np.allclose(np.asarray(d["C0"], float), c_base[np.arange(n)[:, None], Tm], atol=1e-6)
    cum = np.concatenate([np.zeros((n, 1)), np.cumsum(c_base, 1)], 1)
    rows = np.arange(n)[:, None]
    first = Tm // EVERY

    def none_rate(window):
        low, high = (v // EVERY for v in window)
        return np.array([[logs_base[i, first[i, m] + low:first[i, m] + high].sum() for m in range(per)]
                         for i in range(n)]) / (window[1] - window[0])

    W_none, E_none = none_rate(LATER), none_rate(EARLY)
    g_m = np.stack([egc.before(logs_base[i], Tm[i] // EVERY - 1, egc.MAIN) for i in range(n)])
    C_m = c_base[rows, Tm]
    H3b = {}
    for tau in TAUS:
        M_m = (cum[rows, Tm] - cum[rows, Tm - tau]) / tau
        X = np.concatenate([np.stack([M_m, C_m - M_m], -1), g_m], -1)
        late, early = within(W_none, X, picks), within(E_none, X, picks)
        H3b[tau] = dict(b_M=dict(weight=late["weights"][0], bar=late["bars"][0]),
                        b_F=dict(weight=late["weights"][1], bar=late["bars"][1]),
                        a_M=dict(weight=early["weights"][0], bar=early["bars"][0]),
                        a_F=dict(weight=early["weights"][1], bar=early["bars"][1]))
    res["moments"] = H3b

    # H3c: induced minus natural, start by start (new starts only)
    ind_late, ind_early = induced["delayed"].mean(-1), induced["early"].mean(-1)
    res["induced_minus_natural"] = {
        tau: dict(delayed=summary(ind_late - np.array([x[f"b_F:{tau}"] for x in per_start["new"]]), picks),
                  early=summary(ind_early - np.array([x[f"a_F:{tau}"] for x in per_start["new"]]), picks))
        for tau in TAUS}

    # H4a: independent blocks
    res["blocks"] = {}
    for label, name, part in BLOCKS:
        block = per_start[name][part]
        bp = rule.draws(len(block))
        res["blocks"][label] = {f"{m}:{tau}": summary([x[f"{m}:{tau}"] for x in block], bp)
                                for m in TRIO for tau in TAUS}
    res["below_of_60"] = {f"{m}:{tau}": int(sum(x[f"{m}:{tau}"] < 0 for name in ("new", "old") for x in per_start[name]))
                          for m in (*TRIO, "b_F") for tau in TAUS}

    # H4b: non-overlapping rows against all rows, pooled within starts
    T_all = series["new"][0]["T"]
    thin = np.isin(T_all, np.arange(FIRST, LAST + 1, THIN))
    res["pooled"] = {name: dict(all=pooled_trio(series[name], draws[name], np.ones(len(T_all), bool)),
                                thin=pooled_trio(series[name], draws[name], thin)) for name in ("new", "old")}
    res["thin_rows_per_start"] = int(thin.sum())

    # H4c: mismatched starts
    res["mismatched"] = {}
    for name in ("new", "old"):
        fs = series[name]
        mis = [mismatched(fs[i], fs[(i + 1) % len(fs)]) for i in range(len(fs))]
        res["mismatched"][name] = {key: summary([x[key] for x in mis], draws[name]) for key in mis[0]}

    after = fingerprints(inputs)
    res["fingerprints"] = dict(inputs=before, unchanged=before == after, lock_matches=lock_ok)
    assert before == after, "an input changed while this script ran"
    res["checks"] = dict(benchmark_old_18=S["beta_held"]["old"]["mean"], round1_past100=S["check past100"]["new"]["mean"],
                         induced_delayed=res["induced"]["delayed"]["pooled"]["mean"],
                         induced_early=res["induced"]["early"]["pooled"]["mean"])
    OUT_JSON.write_text(json.dumps(res, indent=1))
    lines = write_report(res)
    OUT_TXT.write_text("\n".join(lines) + "\n")
    figure(res)
    print(f"\nWrote {OUT_TXT.name}, {OUT_JSON.name} and {FIGURE.name}")


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

def write_report(res):
    lines = []
    S = res["per_start"]

    def say(text=""):
        print(text)
        lines.append(text)

    def c(s, digits=4, below=False):                       # a cell: mean ± error bar [starts below 0 / n]
        text = f"{s['mean']:+.{digits}f} ± {s['error_bar']:.{digits}f}"
        return text + (f" [{s['below']}/{s['n']}]" if below else "")

    def w(x, digits=4):
        return f"{x['weight']:+.{digits}f} ± {x['bar']:.{digits}f}"

    def m(s, digits=3):                                    # a mean only
        return f"{s['mean']:+.{digits}f}"

    say("Cause test, exploratory follow-up, round 2 (coherence history): October 12, 2026. N = 4,000 per group,"
        " K = 1.90, eps = 0.1")
    say("EXPLORATORY: approved by Colson after round 1; chosen after the cause test's decision was known and run on the")
    say("same saved runs, so not a test. The decision stands: FAIL as registered. Nothing here changes it.")
    say(f"Read-only: the model was not run. Every input fingerprinted before and after: "
        f"{'unchanged' if res['fingerprints']['unchanged'] else 'CHANGED'}; the 16 locked files "
        f"{'still match' if res['fingerprints']['lock_matches'] else 'DO NOT match'} cause_lock.json.")
    say("Data: the 42 new starts (19 to 60) and their 420 moments; the 18 older starts (1 to 18) as a cross-check.")
    say("Cells: mean ± error bar; [k/n] = starts below 0. Error bars from the starts, as in the cause test.")
    ch = res["checks"]
    say(f"Checks: the benchmark on the old 18 ({ch['benchmark_old_18']:+.4f}); round 1's C(T) with W given g and the"
        f" mean over (T - 100, T] ({ch['round1_past100']:+.4f});")
    say(f"  the induced slopes ({ch['induced_delayed']:+.4f} delayed, {ch['induced_early']:+.4f} early).")
    say("Definitions: M_tau = mean C over the tau units before T (T left out); F_tau = C(T) - M_tau; W = growth over")
    say("  (T + 200, T + 500]; E = growth over (T, T + 200]; g = growth in the five 100-unit windows up to T.")
    say("  Every history window ends before T and W starts at T + 200; no control looks past T.")
    say("  H1, H2 and H4 are associations in untouched runs, not effects. M_tau is a summary of the past, not a measured")
    say("  physical memory.")
    say()
    say("H1. THE WINDOW-RESPONSE CURVE (start by start; new 42 starts unless marked old)")
    say(f"  {'tau':>5}  {'M with W, plain':<17}{'M with W, given g':<29}{'slope b, given g':<22}"
        f"{'M with W, given g and C(T)':<29}old: given g | given g, C(T)")
    for tau in TAUS:
        say(f"  {tau:>5}  {m(S[f'r_plain:{tau}']['new'], 4):<17}{c(S[f'r_M|g:{tau}']['new'], 4, True):<29}"
            f"{c(S[f'b_M|g:{tau}']['new'], 4):<22}{c(S[f'r_M|g,C:{tau}']['new'], 4, True):<29}"
            f"{m(S[f'r_M|g:{tau}']['old'])} | {m(S[f'r_M|g,C:{tau}']['old'])}")
    say(f"  for comparison, C(T) with W given g: {c(S['link']['new'])} (old {m(S['link']['old'])})")
    say("  Disjoint bands (mean C over the lags (a, b] before T), partial correlation with W")
    say(f"  {'lags':>13}  {'given g':<29}{'given g and C(T)':<29}old: given g | given g, C(T)")
    for (a, b), tau in zip(BANDS, TAUS):
        say(f"  {f'({a}, {b}]':>13}  {c(S[f'band r|g:{tau}']['new'], 4, True):<29}"
            f"{c(S[f'band r|g,C:{tau}']['new'], 4, True):<29}{m(S[f'band r|g:{tau}']['old'])} | "
            f"{m(S[f'band r|g,C:{tau}']['old'])}")
    say("  Joint fit: W on C(T), the four coarse bands and g; each one's partial correlation given all the others")
    say(f"  {'':<20}{'partial correlation':<29}{'slope':<22}old: partial correlation")
    for name in joint_names():
        say(f"  {name:<20}{c(S[f'joint r:{name}']['new'], 4, True):<29}{c(S[f'joint b:{name}']['new'], 4):<22}"
            f"{m(S[f'joint r:{name}']['old'])}")
    say()
    say("H2. SLOW AND FAST PARTS, C(T) = M_tau + F_tau (W = a + b_M M + b_F F + g terms; slopes are growth per unit C)")
    say(f"  {'tau':>5}  {'mean M':<8}{'SD M':<8}{'SD F':<8}{'b_M':<20}{'b_F':<20}{'M with W, given g, F':<29}"
        f"{'F with W, given g, M':<29}old: M|g,F | F|g,M")
    for tau in TAUS:
        say(f"  {tau:>5}  {S[f'mean_M:{tau}']['new']['mean']:<8.4f}{S[f'sd_M:{tau}']['new']['mean']:<8.4f}"
            f"{S[f'sd_F:{tau}']['new']['mean']:<8.4f}"
            f"{c(S[f'b_M:{tau}']['new'], 3):<20}{c(S[f'b_F:{tau}']['new'], 3):<20}"
            f"{c(S[f'r_M|g,F:{tau}']['new'], 4, True):<29}{c(S[f'r_F|g,M:{tau}']['new'], 4, True):<29}"
            f"{m(S[f'r_M|g,F:{tau}']['old'])} | {m(S[f'r_F|g,M:{tau}']['old'])}")
    say(f"  SD of C(T) within a start: {S['sd_C']['new']['mean']:.4f}; mean C: {S['mean_C']['new']['mean']:.4f}")
    say("  Early growth and natural persistence (E = a + a_M M + a_F F + g terms; later mean C = a + p_M M + p_F F)")
    say(f"  {'tau':>5}  {'a_M':<20}{'a_F':<20}{'F with E, given g, M':<22}{'p_M (T, T+200]':<18}{'p_F (T, T+200]':<18}"
        f"{'p_M (T+200, T+500]':<20}p_F (T+200, T+500]")
    for tau in TAUS:
        say(f"  {tau:>5}  {c(S[f'a_M:{tau}']['new'], 3):<20}{c(S[f'a_F:{tau}']['new'], 3):<20}"
            f"{c(S[f'rE_F|g,M:{tau}']['new'], 4):<22}{c(S[f'p_M early_C:{tau}']['new'], 3):<18}"
            f"{c(S[f'p_F early_C:{tau}']['new'], 3):<18}{c(S[f'p_M later_C:{tau}']['new'], 3):<20}"
            f"{c(S[f'p_F later_C:{tau}']['new'], 3)}")
    say()
    say("H3. THE NATURAL FAST PART AGAINST THE INTERVENTIONS")
    say("  Interventions (the cause test's arms; effects of the specified one-moment pulses, up against down):")
    for name, span in (("delayed", "(T + 200, T + 500]"), ("early", "(T, T + 200]")):
        v = res["induced"][name]
        say(f"    induced slope, growth over {span} per unit change at T: dense {c(v['dense'])},"
            f" sparse {c(v['sparse'])}, pooled {c(v['pooled'])}")
    pi = res["persistence_interventions"]
    say(f"    change at T, (C_up - C_down) / 2: {c(res['change_at_T'], 5)}; persistence (report): {pi['phi']:.3f} over"
        f" (T, T + 200], {pi['left_over']:.3f} over (T + 200, T + 500]")
    say("    M_tau leaves out T, so each pulse moved F_tau by the change at T and left M_tau unchanged, for every tau.")
    say("  Observational (untouched runs; associations, not effects), beside the induced slopes:")
    say(f"  {'tau':>5}  {'b_F, every T':<20}{'b_F, 420 moments':<20}{'induced - b_F':<22}{'a_F, every T':<20}"
        f"{'a_F, 420 moments':<20}{'induced early - a_F':<22}{'p_F (T, T+200]':<16}{'p_F (T+200, T+500]':<20}"
        "SD F / change")
    for tau in TAUS:
        mo, dif = pick(res["moments"], tau), pick(res["induced_minus_natural"], tau)
        say(f"  {tau:>5}  {c(S[f'b_F:{tau}']['new'], 3):<20}{w(mo['b_F'], 3):<20}{c(dif['delayed'], 3, True):<22}"
            f"{c(S[f'a_F:{tau}']['new'], 3):<20}{w(mo['a_F'], 3):<20}{c(dif['early'], 3, True):<22}"
            f"{m(S[f'p_F early_C:{tau}']['new']):<16}{m(S[f'p_F later_C:{tau}']['new']):<20}"
            f"{S[f'sd_F:{tau}']['new']['mean'] / res['change_at_T']['mean']:.2f}")
    say("  At the 420 moments, the slow part (b_M, untouched continuation), for comparison:")
    say("  " + "   ".join(f"{tau}: {w(pick(res['moments'], tau)['b_M'], 3)}" for tau in TAUS))
    say()
    say("H4. STABILITY AND POSSIBLE ARTIFACTS")
    say("  a. Independent blocks of starts (mean ± error bar within each block); starts below 0 of all 60")
    labels = [label for label, _, _ in BLOCKS]
    for measure, title in (("r_M|g", "M with W, given g"), ("r_M|g,F", "M with W, given g and F"),
                           ("r_F|g,M", "F with W, given g and M")):
        say(f"    {title}")
        say(f"    {'tau':>5}  " + "".join(f"{label:<22}" for label in labels) + "below 0")
        for tau in TAUS:
            say(f"    {tau:>5}  " + "".join(f"{c(res['blocks'][label][f'{measure}:{tau}'], 3):<22}" for label in labels)
                + f"{res['below_of_60'][f'{measure}:{tau}']}/60")
    say(f"  b. Non-overlapping rows (T = 2,000, 2,500, ..., 5,500: {res['thin_rows_per_start']} per start), pooled"
        " within starts, against the same pooled measure on all 701 rows")
    for name, label in (("new", "new 42"), ("old", "old 18")):
        say(f"    {label}")
        say(f"    {'tau':>5}  " + "".join(f"{f'{x}, all rows':<24}{f'{x}, 8 rows':<24}" for x in TRIO))
        for tau in TAUS:
            P = res["pooled"][name]
            a, t = pick(P["all"], tau), pick(P["thin"], tau)
            say(f"    {tau:>5}  " + "".join(f"{c(a[x], 3):<24}{c(t[x], 3):<24}" for x in TRIO))
    say("  c. Mismatched starts: coherence of one start, growth of the next (an autocorrelation-only reference)")
    for name, label in (("new", "new 42"), ("old", "old 18")):
        say(f"    {label}: {'':<1}measure: matched mean (SD over starts) | mismatched mean ± error bar (SD) [below 0]")
        mis = res["mismatched"][name]
        say(f"      C(T), given g: {m(S['link'][name], 4)} ({S['link'][name]['sd']:.3f}) | "
            f"{c(mis['link'], 4, True)} ({mis['link']['sd']:.3f})")
        for tau in TAUS:
            cells = []
            for x in TRIO:
                cells.append(f"{x} {m(S[f'{x}:{tau}'][name], 3)} ({S[f'{x}:{tau}'][name]['sd']:.3f}) | "
                             f"{c(mis[f'{x}:{tau}'], 3, True)} ({mis[f'{x}:{tau}']['sd']:.3f})")
            say(f"      {tau:>5}  " + "   ".join(cells))
    say("  d. Slow drift (new 42; old means after the bar)")
    say(f"    {'tau':>5}  " + "".join(f"{label:<20}" for label in ("M|g", "M|g + trend", "M|g 1st half",
                                                                  "M|g 2nd half", "F|g,M", "F|g,M + trend",
                                                                  "F|g,M 1st half", "F|g,M 2nd half")) + "old: M|g + trend, F|g,M + trend")
    for tau in TAUS:
        keys = (f"r_M|g:{tau}", f"r_M|g trend:{tau}", f"r_M|g first:{tau}", f"r_M|g second:{tau}",
                f"r_F|g,M:{tau}", f"r_F|g,M trend:{tau}", f"r_F|g,M first:{tau}", f"r_F|g,M second:{tau}")
        say(f"    {tau:>5}  " + "".join(f"{c(S[k]['new'], 3):<20}" for k in keys)
            + f"{m(S[f'r_M|g trend:{tau}']['old'])}, {m(S[f'r_F|g,M trend:{tau}']['old'])}")
    say("    M with W given g and F, with the trend: " + "  ".join(
        f"{tau}: {m(S[f'r_M|g,F trend:{tau}']['new'])}" for tau in TAUS))
    say("    The run in 1,000-unit blocks (mean C; mean growth rate), averaged over starts:")
    for b in range(6):
        say(f"      ({BLOCK * b}, {BLOCK * (b + 1)}]: new C {c(S[f'drift c:{b}']['new'], 5)}, growth"
            f" {c(S[f'drift growth:{b}']['new'], 5)}; old C {m(S[f'drift c:{b}']['old'], 5)}, growth"
            f" {m(S[f'drift growth:{b}']['old'], 5)}")
    say("  e. Longer growth history (tau > 500): adding the growth over (T - tau, T - 500] to the controls")
    say(f"    {'tau':>5}  {'M|g':<22}{'M|g + longer growth':<26}{'F|g,M':<22}{'F|g,M + longer growth':<26}old: last two")
    for tau in LONG:
        say(f"    {tau:>5}  {c(S[f'r_M|g:{tau}']['new'], 4):<22}{c(S[f'r_M|g,G:{tau}']['new'], 4):<26}"
            f"{c(S[f'r_F|g,M:{tau}']['new'], 4):<22}{c(S[f'r_F|g,M,G:{tau}']['new'], 4):<26}"
            f"{m(S[f'r_M|g,G:{tau}']['old'])}, {m(S[f'r_F|g,M,G:{tau}']['old'])}")
    say()
    say("H5. Every window, band and control that was run is reported above, nulls and contradictions included. No")
    say("    window is singled out: the strongest one is not used for anything, and no threshold was tuned here.")
    say()
    say("H6. WHAT THE SAVED DATA CANNOT RESOLVE")
    for text in (
            "1. Whether coherence history is a cause. Every history comparison here is between moments that differ"
            " naturally. No saved run changed coherence over a stretch of time (the cause test changed it at one"
            " instant), so 'history drives later growth' cannot be told apart from 'a slow hidden state drives both'.",
            "2. What M_tau stands for physically. The saved runs keep only each group's z = R e^(i psi) at every time"
            " unit, not the oscillators' phases, so no frequency-resolved measure can be computed (for example, how"
            " the slow, nearly locked oscillators are arranged).",
            "3. The higher harmonics over time and their phases: only |z2|, |z3|, |z4| at the moments were saved.",
            "4. The gap's direction: only its growth was saved, so whether the carried gap's alignment carries part of"
            " the link, or of the cause test's delayed effects, cannot be checked. That needs a fresh gap.",
            "5. Anything later than 500 units after a moment: the branches stop at T + 500.",
            "6. Dose against level: each pulse was +-50% of the level at that moment, so the two cannot be separated.",
            "7. Other settings: one coupling (K = 1.90, just below the threshold), one N and one eps. Whether the role of"
            " history depends on the distance from the threshold is unknown here.",
            "8. The window-response curve describes these runs. Its shape, or the window where it is strongest, is not"
            " an estimate of a physical memory time."):
        words, line = text.split(), "   "
        for word in words:
            if len(line) + len(word) + 1 > 118:
                say(line)
                line = "      "
            line += " " + word
        say(line)
    say()
    say("Changes after the first run: " + ("; ".join(res["changes_after_first_run"]) or "none"))
    return lines


# ---------------------------------------------------------------------------------------------------------------------
# the figure

def figure(res):
    import matplotlib
    matplotlib.use("Agg")
    import matplotlib.pyplot as plt

    S = res["per_start"]
    taus = np.array(TAUS, float)

    def series(key, name="new"):
        return (np.array([S[f"{key}:{tau}"][name]["mean"] for tau in TAUS]),
                np.array([S[f"{key}:{tau}"][name]["error_bar"] for tau in TAUS]))

    fig, axes = plt.subplots(2, 2, figsize=(16, 11.5))
    fig.subplots_adjust(left=0.07, right=0.985, top=0.885, bottom=0.11, hspace=0.36, wspace=0.2)
    fig.suptitle("Cause test, exploratory follow-up, round 2: coherence history (October 12, 2026). EXPLORATORY:"
                 " chosen after the FAIL was known; not a test\nUntouched runs, N = 4,000 per group, K = 1.90; 42 new"
                 " starts (filled) and 18 old (open). The decision stays FAIL as registered.", fontsize=13)

    ax = axes[0, 0]
    styles = (("r_M|g", "#222222", "D", "history mean M, given g"),
              ("r_M|g,F", "#2b8cbe", "o", "M, given g and the fast part F"),
              ("r_F|g,M", "#e66101", "s", "fast part F = C(T) - M, given g and M"))
    for key, color, marker, label in styles:
        y, e = series(key)
        ax.errorbar(taus, y, yerr=e, fmt=marker + "-", color=color, ms=6, lw=1.4, capsize=3, label=label)
        ax.plot(taus, series(key, "old")[0], marker, mfc="white", mec=color, mew=1.5, ms=6)
    mis = res["mismatched"]["new"]
    my = np.array([mis[f"r_M|g:{tau}"]["mean"] for tau in TAUS])
    me = np.array([mis[f"r_M|g:{tau}"]["error_bar"] for tau in TAUS])
    ax.fill_between(taus, my - 2 * me, my + 2 * me, color="0.75", alpha=0.6,
                    label="M given g, mismatched starts (± 2 error bars)")
    ax.axhline(S["link"]["new"]["mean"], color="0.3", ls=":", lw=1, label="C(T), given g (no history)")
    ax.axhline(0, color="0.5", lw=0.9)
    ax.set_xscale("log")
    ax.set(xlabel="history window tau (time units before T)", ylabel="partial correlation with W",
           title="The window-response curve: every window, none selected")
    ax.legend(fontsize=8, loc="lower left")

    ax = axes[0, 1]
    mids = np.array([np.sqrt(max(a, 1) * b) for a, b in BANDS])
    for key, color, marker, label in (("band r|g", "#222222", "D", "band, given g"),
                                      ("band r|g,C", "#7570b3", "^", "band, given g and C(T)")):
        y, e = series(key)
        ax.errorbar(mids, y, yerr=e, fmt=marker + "-", color=color, ms=6, lw=1.4, capsize=3, label=label)
        ax.plot(mids, series(key, "old")[0], marker, mfc="white", mec=color, mew=1.5, ms=6)
    ax.axhline(0, color="0.5", lw=0.9)
    ax.set_xscale("log")
    ax.set(xlabel="lag band before T (plotted at its geometric middle)", ylabel="partial correlation with W",
           title="Disjoint history bands (no band overlaps another)")
    ax.legend(fontsize=8.5, loc="lower left")

    ax = axes[1, 0]
    for key, color, marker, label in (("b_M", "#222222", "D", "b_M: slow part (every T)"),
                                      ("b_F", "#e66101", "s", "b_F: fast part (every T)")):
        y, e = series(key)
        ax.errorbar(taus, y, yerr=e, fmt=marker + "-", color=color, ms=6, lw=1.4, capsize=3, label=label)
    yb = np.array([pick(res["moments"], tau)["b_F"]["weight"] for tau in TAUS])
    eb = np.array([pick(res["moments"], tau)["b_F"]["bar"] for tau in TAUS])
    ax.errorbar(taus * 1.08, yb, yerr=eb, fmt="s", mfc="white", mec="#e66101", ecolor="#e66101", ms=6, capsize=2,
                label="b_F at the 420 moments (untouched)")
    ind = res["induced"]["delayed"]
    lo, hi = ind["pooled"]["mean"] - ind["pooled"]["error_bar"], ind["pooled"]["mean"] + ind["pooled"]["error_bar"]
    ax.axhspan(lo, hi, color="#d7301f", alpha=0.25,
               label=f"intervention effect, pooled ({ind['pooled']['mean']:+.3f}; a pulse moves F only)")
    for method, ls in (("dense", "--"), ("sparse", ":")):
        ax.axhline(ind[method]["mean"], color="#d7301f", ls=ls, lw=1.2, label=f"intervention effect, {method}")
    ax.axhline(0, color="0.5", lw=0.9)
    ax.set_xscale("log")
    ax.set(xlabel="history window tau", ylabel="growth over (T + 200, T + 500] per unit coherence",
           title="Slopes: observational (untouched runs) against the pulses' effect")
    ax.legend(fontsize=7.8, loc="lower left", ncol=2)

    ax = axes[1, 1]
    colors = ("#1b9e77", "#66c2a5", "#7570b3", "#9e9ac8", "#cbc9e2")
    for (label, _, _), color in zip(BLOCKS, colors):
        y = np.array([res["blocks"][label][f"r_M|g:{tau}"]["mean"] for tau in TAUS])
        ax.plot(taus, y, "-", color=color, lw=1.6, label=f"starts {label}")
    for key, ls, label in (("r_M|g trend", "--", "new 42, time trend held"),
                           ("r_M|g first", "-.", "new 42, first half of the run"),
                           ("r_M|g second", ":", "new 42, second half of the run")):
        ax.plot(taus, series(key)[0], ls, color="k", lw=1.3, label=label)
    thin = res["pooled"]["new"]["thin"]
    ty = np.array([pick(thin, tau)["r_M|g"]["mean"] for tau in TAUS])
    te = np.array([pick(thin, tau)["r_M|g"]["error_bar"] for tau in TAUS])
    ax.errorbar(taus, ty, yerr=te, fmt="s", color="#d95f02", ms=5, capsize=2,
                label="new 42, non-overlapping rows, pooled")
    ax.axhline(0, color="0.5", lw=0.9)
    ax.set_xscale("log")
    ax.set(xlabel="history window tau", ylabel="partial correlation of M with W, given g",
           title="Stability: independent blocks, non-overlapping rows, drift checks")
    ax.legend(fontsize=7.8, loc="lower left", ncol=2)

    for y, text in ((0.05, "Exploratory: every analysis was listed before it was run, on saved results only (no new"
                           " simulation). Error bars from the starts, as in the cause test."),
                    (0.025, "Observational slopes are associations in untouched runs; only the red band is an effect of"
                            " an intervention. A toy model: structure for choosing the next question, not proof.")):
        fig.text(0.07, y, text, fontsize=9.5, color="0.25")
    fig.savefig(FIGURE, dpi=110)


if __name__ == "__main__":
    main()
