#!/usr/bin/env python3 """ Cause test (October 12, 2026): if coherence is raised or lowered at a moment, does the growth of a tiny gap 200 to 500 time units later change as the correlation in the untouched runs predicts? The plan, fixed before the run, is cause_test_prereg.md. This file carries it out with nudge_test.py's model and step code (imported unchanged) and makes the decision with cause_rule.py (unchanged). 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 cause_test.py check # the checks before the lock, on old start 1; prints them, writes nothing python3 -B cause_test.py targets # the planned changes at every multiple of 5 from t = 2,000 to 5,495 in old # starts 13 to 18 (Test 17's replication starts): immediate changes only, # no branch is run on; writes cause_targets.txt and cause_targets.json python3 -B cause_test.py dryrun # both kinds of task on old start 1 for a short stretch, the pieces, the # assembly, and the report on made-up numbers; writes only to /tmp python3 -B cause_test.py freeze # the checks again, then the lock; refused once any new data exist python3 -B cause_test.py run # one piece: the next unfinished tasks (7 base runs or 14 moment tasks) python3 -B cause_test.py report # once all 462 tasks are done: assemble, check, decide, write the outputs Each is run with a hard time limit under 5 minutes (timeout -k 10 285 ...). A piece cut off leaves no part for a task it did not finish, and the next piece runs that task again from the beginning. Needs Python 3, numpy, scipy and matplotlib; cause_rule.py and design_numbers.py here; ../coherence-growth-2026-10-11 and ../nudge-test-2026-10-10 with their saved runs, unchanged. """ import csv import hashlib import json import os import platform import shutil import sys import time from datetime import datetime, timezone from multiprocessing import Pool from pathlib import Path from zoneinfo import ZoneInfo os.environ["OPENBLAS_NUM_THREADS"] = os.environ["OMP_NUM_THREADS"] = "1" # one thread per worker (set before numpy) import numpy as np # noqa: E402 HERE = Path(__file__).resolve().parent ROOT = HERE.parent.parent CG = HERE.parent / "coherence-growth-2026-10-11" sys.dont_write_bytecode = True # leave every folder as it is sys.path.insert(0, str(HERE)) sys.path.insert(1, str(CG)) import cause_rule as rule # noqa: E402 (the decision rule, fixed with the plan) import coherence_growth as cg # noqa: E402 (the model's settings and the saved runs; unchanged) import design_numbers as dn # noqa: E402 (the plan's numbers and moment table; unchanged) import earlier_growth_check as egc # noqa: E402 (Test 17's measure, for the base runs; unchanged) nt = cg.nt N, K, EPS, DT = cg.N, cg.K, nt.EPS, nt.DT NEW = dn.NEW # starts 19 to 60 (seeds 18 to 59) RUN, EVERY, GAP = cg.RUN, nt.LYAP_EVERY, nt.LYAP_GAP # 6,000 time units; the gap put back to 1e-7 every 5 MOMENTS, DRAWS = dn.MOMENTS, dn.SEED # ten moments per start; every draw from 20261012 FRACTION, SPARSE = dn.FRACTION, dn.SPARSE # coherence x 1.5 or x 0.5; sparse: a tenth of the pulled ARMS, METHODS = rule.ARMS, rule.METHODS # a moment's six branches: dense up, down, sham; sparse ... GROUPS = ("A", "B") BANDS = dict(up=(1.45, 1.55), down=(0.45, 0.55), sham=(0.95, 1.05)) # R_g after / R_g before, each group SCAN, LARGEST, BISECT, SOLVED = 256, np.pi / 2, 1e-12, 1e-6 # amplitude: scan to pi/2, halve to 1e-12; |h/target - 1| MATCH, KEPT = 1e-9, 1e-6 # M1: sums of squared moves alike; the gap as it was AFTER = 500.0 # each branch runs this long after T WINDOWS = dict(later=dn.LATER, early=(0.0, 200.0), whole=(0.0, 500.0)) # growth after T; "later" is the outcome HALF, BENCHMARK, NATURAL, FLOOR = dn.HALF, 0.1602, 0.363, 0.181 # the threshold; natural persistence; M2 WORKERS, PIECE = 7, dict(base=7, moment=14) # workers; tasks per piece (each about 3 minutes) CHECK_RUN, CHECK_AT = cg.CHECK_RUN, 25 # the checks: start 1 for 50 time units, branches from t = 25 TARGET_STARTS, TARGET_SPAN, TARGET_KEY = tuple(range(12, 18)), (2000, 5495), 100 # old starts 13 to 18 HALIFAX = ZoneInfo("America/Halifax") PLAN, LOCK, RUNLOG = HERE / "cause_test_prereg.md", HERE / "cause_lock.json", HERE / "cause_runlog.json" PARTS, SAVED = HERE / "cause_parts", HERE / "sim_cause4000.npz" OUT_TXT, OUT_JSON = HERE / "cause_test.txt", HERE / "cause_test.json" TARGETS_TXT, TARGETS_JSON = HERE / "cause_targets.txt", HERE / "cause_targets.json" TEST16 = nt.HERE / "lyapunov_long.json" FIGURE = ROOT / "screenshots" / "nudge-test-cause-test.png" DRY = Path("/tmp/cause-dryrun") PAGE_SHA = "b06730e7772be1f719cadae7924cdca3c2e0ebb89efb17104969363431f782e7" # nudge_test.py as the page carries it PLAN_SHA = "1c44034e64f4ec74d5afba7080cf1df1e96262680583f785447b153f3ce61112" # the plan as Colson approves it; the freeze refuses any other PLAN_FILES = dict( # the plan's table of files fixed with it cause_rule=(HERE / "cause_rule.py", "5fdd1c90df371675d1113a5eca799abdc2c639c00adbeaeb0339c05b02f487ca"), design_numbers=(HERE / "design_numbers.py", "103241b2fb361649f3e53cea900fedb1ee475c11200e8e8e8feeb2c2af4fecd3"), design_text=(dn.OUT_TXT, "76d68062720f21e7c5f62ae14092e6c135ae207b036564ccb30e1a9770edca39"), design_json=(dn.OUT_JSON, "39bc39762c551631327d5120b69a6b5daf27f7a6c5cdb5f06fe59a22b281af28"), moment_table=(dn.MOMENT_TABLE, "1a1d72113ab035b14b9e80e413e7030291915d8d1a365ad1f658f95c277709a8"), targets_text=(TARGETS_TXT, "5e3f8a8dba33a11bad8a0693e65803cfbfdb3bd83bfbaa7e868d41b02882d30c"), targets_json=(TARGETS_JSON, "225058b7c11159c412ced9b2f1d8c7ea48d9e48d2805b913be4bdc2278acbe19"), nudge_test=(nt.HERE / "nudge_test.py", PAGE_SHA)) LOCKED = dict( # every file the result depends on plan=PLAN, cause_test=Path(__file__).resolve(), **{name: path for name, (path, _) in PLAN_FILES.items()}, coherence_growth=CG / "coherence_growth.py", earlier_growth_check=CG / "earlier_growth_check.py", test16_results=TEST16, earlier_runs=cg.EARLIER, fresh_runs=cg.SAVED, replication_runs=dn.REPLICATION) SETTINGS = dict( n=N, k=K, eps=EPS, delta=nt.DELTA, dt=DT, settle=nt.SETTLE, run=RUN, starts=[s + 1 for s in NEW], gap=GAP, put_back_every=EVERY, gap_seed=nt.LYAP_SEED, moments=MOMENTS, spacing=dn.SPACING, first=dn.FIRST, draws=DRAWS, fraction=FRACTION, sparse=SPARSE, bands={arm: list(band) for arm, band in BANDS.items()}, scan=SCAN, largest=LARGEST, bisect=BISECT, solved=SOLVED, match=MATCH, gap_kept=KEPT, after=AFTER, later=list(WINDOWS["later"]), benchmark=BENCHMARK, half=HALF, natural=NATURAL, floor=FLOOR, most_missing=rule.MOST_MISSING, fewest_kept=rule.FEWEST_KEPT, rule=rule.RULE, side=rule.SIDE, resamples=rule.RESAMPLES, resample_seed=rule.RESAMPLE_SEED, pieces=PIECE, workers=WORKERS) FROZEN = dict( # the plan's numbers, which the code must match n=4000, k=1.9, eps=0.1, delta=1.0, dt=0.025, settle=300.0, run=6000.0, starts=list(range(19, 61)), gap=1e-7, put_back_every=5.0, gap_seed=20261011, moments=10, spacing=350.0, first=2000.0, draws=20261012, fraction=0.5, sparse=0.1, bands=dict(up=[1.45, 1.55], down=[0.45, 0.55], sham=[0.95, 1.05]), scan=256, largest=np.pi / 2, bisect=1e-12, solved=1e-6, match=1e-9, gap_kept=1e-6, after=500.0, later=[200.0, 500.0], benchmark=0.1602, half=0.5, natural=0.363, floor=0.181, most_missing=20, fewest_kept=5, rule=3.0, side=2.0, resamples=400, resample_seed=20261015, pieces=dict(base=7, moment=14), workers=7) MEANING = dict( PASS=("in this toy model, the interventions changed later growth as predicted: raising coherence at a moment slowed" " the growth of tiny gaps 200 to 500 time units later, and lowering it sped that growth up, each against a" " matched sham, in both methods, by at least half the size the correlation predicts. This supports a causal" " effect of these interventions within the model. It does not show that coherence is the only route (the" " changes also rearrange the phases in other ways, recorded below), and it says nothing beyond the model."), SMALL=("a real effect of the interventions in the predicted direction, but under half the benchmark: most of the" " correlation would come from something else."), PARTIAL=("raising and lowering differ as predicted, but up or down is not clearly apart from its sham, or one method" " does not show it."), REVERSED="the interventions move later growth the opposite way to the correlation.", FAIL=("no effect of the interventions as large as half the benchmark. The correlation in the untouched runs is then" " most likely carried by something else, such as a slower feature of the state that drives both."), INCONCLUSIVE=("not a negative result: the change did not happen as planned, the branches were not matched, the" " change faded too fast, or the data cannot tell a benchmark-sized effect from none.")) # --------------------------------------------------------------------------------------------------------------------- # small helpers def stamp(t=None): moment = datetime.fromtimestamp(time.time() if t is None else t, timezone.utc) return dict(utc=moment.isoformat(timespec="seconds"), halifax=moment.astimezone(HALIFAX).isoformat(timespec="seconds")) def shown(moment): return moment["halifax"][:16].replace("T", " ") + " Halifax" def sha256(path): return hashlib.sha256(Path(path).read_bytes()).hexdigest() def write_json(path, data): """Written whole or not at all.""" partial = path.with_name(path.name + ".partial") partial.write_text(json.dumps(plain(data), indent=1) + "\n") os.replace(partial, path) def plain(x): """Numbers and arrays as JSON holds them (not-a-number as null).""" if isinstance(x, dict): return {str(key): plain(value) for key, value in x.items()} if isinstance(x, (list, tuple)): return [plain(value) for value in x] if isinstance(x, np.ndarray): return plain(x.tolist()) if isinstance(x, np.bool_): return bool(x) if isinstance(x, np.integer): return int(x) if isinstance(x, (np.floating, float)): return float(x) if np.isfinite(x) else None if isinstance(x, np.str_): return str(x) return x def fingerprints(): return {name: sha256(path) for name, path in LOCKED.items()} def frozen(): """The settings in the code are the plan's, and the plan's rounded numbers are design_numbers.json's.""" now = json.loads(json.dumps(SETTINGS)) wrong = [key for key, value in FROZEN.items() if now.get(key) != value] design = json.loads(dn.OUT_JSON.read_text()) for name, written, worked_out, digits in ( ("benchmark", BENCHMARK, design["benchmark"]["beta"], 4), ("natural", NATURAL, design["persistence"]["0-200"]["mean"], 3), ("floor", FLOOR, design["benchmark"]["persistence_floor"], 3)): if round(worked_out, digits) != written: wrong.append(f"{name} ({written} in the plan, {worked_out} worked out)") if wrong: sys.exit(f"The settings differ from the plan: {', '.join(wrong)}") def unchanged(): """The lock, after checking that nothing it covers has changed since.""" if not LOCK.exists(): sys.exit("No lock yet: run `python3 -B cause_test.py freeze` first, before any new start is run.") lock = json.loads(LOCK.read_text()) now = fingerprints() changed = [name for name in LOCKED if lock["sha256"].get(name) != now[name]] if lock["settings"] != json.loads(json.dumps(SETTINGS)): changed.append("settings") if changed: sys.exit(f"Changed since the lock ({shown(lock['locked'])}): {', '.join(changed)}. Stopping.") return lock def wrap(x): """An angle (or a change of angle) in (-pi, pi].""" return np.angle(np.exp(1j * np.asarray(x))) def rms(x): return float(np.sqrt(np.mean(np.square(x)))) def harmonic(us, uf, k=1): """Each group's z_k = (1/N) sum_j e^{i k theta_j} over all N oscillators, (row, group). For k = 1 it is the coherence exactly as nudge_test measures it.""" if k == 1: return (us.sum(-1) + uf.sum(-1)) / N return ((us ** k).sum(-1) + (uf ** k).sum(-1)) / N def gap_of(first, second): """The gap as the base run measures it: the phase differences of the two copies over both groups' pulled oscillators, with the overall timing lined up (the mean removed).""" shift = np.angle(second * np.conj(first)) return shift - shift.mean() def saved_run(s): """Old start s (by seed) as saved: its untouched copy's coherence (group, time), gap growth (interval) and settling (group, 61).""" if s < 6: with np.load(cg.EARLIER) as old: i = [round(float(x), 2) for x in old["ks"]].index(K) return dict(z=old["z"][i, s], logs=old["logs"][i, s], trace=old["trace"][i, s]) for path in (cg.SAVED, dn.REPLICATION): with np.load(path) as data: starts = [int(x) for x in data["starts"]] if s in starts: j = starts.index(s) return dict(z=data["z"][j], logs=data["logs"][j], trace=data["trace"][j]) raise KeyError(f"start {s + 1} is not saved") # --------------------------------------------------------------------------------------------------------------------- # the model runs def base_run(s, keep=(), run=RUN, visit=None): """Start s (seed), step for step as nudge_test's lyapunov runs one start without a deadline: settling, the second copy with the tiny random gap, then `run` time units with the gap measured and put back every 5. The full state (both copies, just after the gap was put back) is kept at each time in `keep`. visit(t, us, uf), if given, is called just after every put-back and must leave the state as it is (the target check).""" clock = time.time() slow, stepper = nt.model(N, K, EPS) run_steps, reset_every, sample_every = (int(round(x / DT)) for x in (run, EVERY, 1.0)) def coherence(us, uf): return (us.sum(-1) + uf.sum(-1)) / N us, uf = nt.starting_phases(N, slow, s, 1) settle_steps, trace_every = int(round(nt.SETTLE / DT)), int(round(5.0 / DT)) trace = [np.abs(coherence(us, uf))] step, none = stepper(DT), np.zeros(1) for j in range(1, settle_steps + 1): step(us, uf, none) if j % trace_every == 0: trace.append(np.abs(coherence(us, uf))) trace = np.stack(trace, -1).astype(np.float32) offset = np.stack([np.random.default_rng((nt.LYAP_SEED, s)).standard_normal(us.shape[1:])]) offset -= offset.mean(axis=(-2, -1), keepdims=True) offset *= GAP / np.sqrt((offset * offset).mean(axis=(-2, -1), keepdims=True)) us, uf = np.concatenate([us, us * np.exp(1j * offset)]), np.concatenate([uf, uf]) samples, logs, states = [coherence(us[:1], uf[:1])], [], [] keep_steps = {int(round(t / DT)) for t in keep} step, none = stepper(DT), np.zeros(2) for i in range(1, run_steps + 1): step(us, uf, none) if i % sample_every == 0: samples.append(coherence(us[:1], uf[:1])) if i % reset_every == 0: shift = np.angle(us[1:] * np.conj(us[:1])) shift -= shift.mean(axis=(-2, -1), keepdims=True) size = np.sqrt((shift * shift).mean(axis=(-2, -1))) logs.append(np.log(size / GAP)) us[1:] = us[:1] * np.exp(1j * shift * (GAP / size)[:, None, None]) if i in keep_steps: states.append((us.copy(), uf.copy())) if visit is not None: visit(i // sample_every, us, uf) return dict(z=np.stack(samples, -1), logs=np.stack(logs, -1), trace=trace, states=states, seconds=time.time() - clock) def branches(us, uf, dtheta, length): """Every branch of a moment from its saved state (us: (copy, group, pulled), uf: (copy, group, free)), one after another. Branch b has the change dtheta[b] applied to the same oscillators of both copies; then it runs exactly as the base run does, with the two copies as its two rows (the base run's own shapes), for `length` time units. Returns each branch's coherence every time unit from T, (branch, group, time), and its gap growth in each 5-unit interval, (branch, interval).""" slow, stepper = nt.model(N, K, EPS) steps, reset_every, sample_every = (int(round(x / DT)) for x in (length, EVERY, 1.0)) paths, growths = [], [] for change in dtheta: U, F = us * np.exp(1j * change), uf.copy() # the same change on both copies step, none = stepper(DT), np.zeros(2) samples, logs = [harmonic(U[:1], F[:1])[0]], [] for i in range(1, steps + 1): step(U, F, none) if i % sample_every == 0: samples.append(harmonic(U[:1], F[:1])[0]) if 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 / GAP)[0]) U[1:] = U[:1] * np.exp(1j * shift * (GAP / size)[:, None, None]) paths.append(np.stack(samples, -1)) growths.append(np.array(logs)) return np.stack(paths), np.stack(growths) # --------------------------------------------------------------------------------------------------------------------- # the planned change at a moment def paired_signs(along, draw): """The sham's signs for one group: the oscillators ranked by sin^2(psi - theta), taken in pairs of neighbours; a draw decides which one of each pair moves as in up (+1) and which as in down (-1). If the count is odd, the last one gets a random sign.""" order = np.argsort(along * along, kind="stable") count = len(order) pairs = count // 2 bits = draw.integers(0, 2, pairs + count % 2) first = np.where(bits[:pairs] == 0, 1.0, -1.0) signs = np.empty(count) signs[order[0:2 * pairs:2]] = first signs[order[1:2 * pairs:2]] = -first if count % 2: signs[order[-1]] = 1.0 if bits[-1] == 0 else -1.0 return signs def amplitude(u, along, z, target): """The smallest a in (0, pi/2] with h(a) = [R(+a) - R(-a)] / 2 = target, where R(+-a) = |z + (1/N) sum_j u_j (e^{+-i a along_j} - 1)| is the group's coherence after moving the chosen oscillators (u: their e^{i theta}) by +-a along_j. h is worked out at a = pi/512, 2 pi/512, ..., pi/2; the first value where it reaches the target and the one before are then halved between until less than 1e-12 apart. Returns (a, 0); (nan, 1) if h never reaches the target (infeasible); (nan, 2) if a number is not finite or h(a) misses the target by more than a relative 1e-6 (a solver failure).""" real, imag = u.real.copy(), u.imag.copy() def h(a): angle = np.multiply.outer(np.atleast_1d(a), along) # (value of a, oscillator) c, s = np.cos(angle) - 1, np.sin(angle) cr, ci, sr, si = c @ real, c @ imag, s @ real, s @ imag plus = np.abs(z + ((cr - si) + 1j * (ci + sr)) / N) # sum_j u_j (e^{+i a along_j} - 1) / N minus = np.abs(z + ((cr + si) + 1j * (ci - sr)) / N) # sum_j u_j (e^{-i a along_j} - 1) / N return (plus - minus) / 2 if not (np.isfinite(target) and target > 0 and np.all(np.isfinite(along))): return np.nan, 2 grid = LARGEST * np.arange(1, SCAN + 1) / SCAN values = h(grid) if not np.all(np.isfinite(values)): return np.nan, 2 reached = np.flatnonzero(values >= target) 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 > BISECT: middle = (low + high) / 2 if h(middle)[0] >= target: high = middle else: low = middle value = h(high)[0] if not (np.isfinite(value) and abs(value / target - 1) <= SOLVED): return np.nan, 2 return high, 0 def state_numbers(u, f): """The first copy's state: each group's R_g, psi_g, |z_2|, |z_3|, |z_4|; C; psi_A - psi_B.""" zs = [harmonic(u[None], f[None], k)[0] for k in (1, 2, 3, 4)] R, psi = np.abs(zs[0]), np.angle(zs[0]) return dict(R=R, psi=psi, z2=np.abs(zs[1]), z3=np.abs(zs[2]), z4=np.abs(zs[3]), C=float(R.mean()), ab=float(wrap(psi[0] - psi[1]))) def record_template(): """Every immediate record of a moment, before it is filled in (not a number where nothing was built).""" nan = np.nan return dict( R0=np.full(2, nan), psi0=np.full(2, nan), z2_0=np.full(2, nan), z3_0=np.full(2, nan), z4_0=np.full(2, nan), C0=nan, ab0=nan, gap0=nan, R1=np.full((6, 2), nan), psi1=np.full((6, 2), nan), z2_1=np.full((6, 2), nan), z3_1=np.full((6, 2), nan), z4_1=np.full((6, 2), nan), ratio=np.full((6, 2), nan), inband=np.zeros((6, 2), bool), C1=np.full(6, nan), ab1=np.full(6, nan), sumsq=np.full(6, nan), moved=np.full(6, nan), largest=np.full(6, nan), rms=np.full(6, nan), gap1=np.full(6, nan), gaprel=np.full(6, nan), amp=np.full((2, 2), nan), status=np.zeros((2, 2), int), built=np.zeros(2, bool), same=np.zeros(2, bool), match=np.full(2, nan), carried=False, code_error=False, reason="") RECORD_KEYS = tuple(record_template()) def planned(us, uf, key): """The six changes of one moment from its saved state (both copies: us (copy, group, pulled), uf (copy, group, free)), worked out from the first copy, with every check made before any branch runs. key: (20261012, seed of the start, moment); the draws are key + (0,) for the dense sham's signs, (1,) for the sparse subsets, (2,) for the sparse sham's signs, group A before B. For each method and group: the amplitude by amplitude(); up moves each chosen oscillator by +a sin(psi - theta), down by -a sin(psi - theta), sham by the same size with paired signs. Each branch must take each group's coherence to its band (up 1.45 to 1.55 times R_g, down 0.45 to 0.55, sham 0.95 to 1.05). The moment is carried out only if every target was feasible and solved, every branch is inside its bands, the arms of each method move the same oscillators by the same sizes (sums of squares within a relative 1e-9) and the gap is as it was (within a relative 1e-6). Otherwise none of its branches is run, the reason is kept, and nothing is clipped or replaced. Returns the records and the changes, (branch, group, pulled).""" u, f = us[0], uf[0] count = u.shape[-1] z = harmonic(u[None], f[None])[0] R, psi = np.abs(z), np.angle(z) along = np.sin(psi[:, None] - np.angle(u)) # sin(psi_g - theta_j), (group, pulled) draw = [np.random.default_rng((*key, use)) for use in range(3)] dense_signs = [paired_signs(along[g], draw[0]) for g in range(2)] chosen = [np.sort(draw[1].choice(count, int(round(SPARSE * count)), replace=False)) for g in range(2)] sparse_signs = [paired_signs(along[g, chosen[g]], draw[2]) for g in range(2)] sets = dict(dense=[np.arange(count)] * 2, sparse=chosen) signs = dict(dense=dense_signs, sparse=sparse_signs) rec = record_template() dtheta = np.zeros((6, 2, count)) reasons = [] before = state_numbers(u, f) rec.update(R0=before["R"], psi0=before["psi"], z2_0=before["z2"], z3_0=before["z3"], z4_0=before["z4"], C0=before["C"], ab0=before["ab"]) d0 = gap_of(us[0], us[1]) rec["gap0"] = rms(d0) for k, method in enumerate(METHODS): for g in range(2): idx = sets[method][g] rec["amp"][k, g], rec["status"][k, g] = amplitude(u[g, idx], along[g, idx], z[g], FRACTION * R[g]) if rec["status"][k, g]: what = "target infeasible" if rec["status"][k, g] == 1 else "solver failure" reasons.append(f"{method} {what}, group {GROUPS[g]} (R = {R[g]:.4f})") if (rec["status"][k] != 0).any(): continue rec["built"][k] = True for g in range(2): idx = sets[method][g] move = rec["amp"][k, g] * along[g, idx] dtheta[3 * k, g, idx] = move dtheta[3 * k + 1, g, idx] = -move dtheta[3 * k + 2, g, idx] = signs[method][g] * move moved = sum(len(sets[method][g]) for g in range(2)) for a, arm in enumerate(ARMS): b = 3 * k + a turn = np.exp(1j * dtheta[b]) after = state_numbers(u * turn, f) rec["R1"][b], rec["psi1"][b], rec["C1"][b], rec["ab1"][b] = after["R"], after["psi"], after["C"], after["ab"] rec["z2_1"][b], rec["z3_1"][b], rec["z4_1"][b] = after["z2"], after["z3"], after["z4"] rec["ratio"][b] = after["R"] / R rec["inband"][b] = (BANDS[arm][0] <= rec["ratio"][b]) & (rec["ratio"][b] <= BANDS[arm][1]) for g in np.flatnonzero(~rec["inband"][b]): reasons.append(f"{method} {arm}, group {GROUPS[g]} at {rec['ratio'][b, g]:.4f} x R") d1 = gap_of(us[0] * turn, us[1] * turn) rec["gap1"][b], rec["gaprel"][b] = rms(d1), rms(d1 - d0) / rms(d0) rec["sumsq"][b] = float(np.square(dtheta[b]).sum()) rec["moved"][b], rec["largest"][b] = moved, float(np.abs(dtheta[b]).max()) rec["rms"][b] = float(np.sqrt(rec["sumsq"][b] / moved)) sizes = np.abs(dtheta[3 * k:3 * k + 3]) rec["same"][k] = np.array_equal(sizes[0], sizes[1]) and np.array_equal(sizes[0], sizes[2]) sums = rec["sumsq"][3 * k:3 * k + 3] rec["match"][k] = float(np.abs(sums[1:] - sums[0]).max() / sums[0]) if not (rec["same"][k] and rec["match"][k] <= MATCH): rec["code_error"] = True reasons.append(f"{method} arms not matched") if not np.all(rec["gaprel"][3 * k:3 * k + 3] < KEPT): rec["code_error"] = True reasons.append(f"{method} gap changed by the change") rec["carried"] = bool(rec["built"].all() and rec["inband"].all() and not rec["code_error"]) rec["reason"] = "; ".join(reasons) return rec, dtheta # --------------------------------------------------------------------------------------------------------------------- # tasks, pieces and parts def tasks(parts=PARTS): """All 462 tasks in the fixed order: the 42 base runs, then the 420 moments, start by start, moment by moment.""" table = dn.moment_table() keep = {} for start, m, T in table: keep.setdefault(start - 1, []).append(T) base = [dict(name=f"base_{s + 1:02d}", kind="base", seed=s, keep=tuple(keep[s]), run=RUN, parts=str(parts)) for s in NEW] moments = [dict(name=f"moment_{start:02d}_{m:02d}", kind="moment", seed=start - 1, moment=m - 1, T=T, after=AFTER, parts=str(parts)) for start, m, T in table] return base + moments def base_task(task): res = base_run(task["seed"], keep=task["keep"], run=task["run"]) if len(res["states"]) != len(task["keep"]): raise RuntimeError("a moment's state was not kept") return dict(start=task["seed"] + 1, moments=np.array(task["keep"]), z=res["z"].astype(np.complex64), logs=res["logs"], trace=res["trace"], us=np.stack([u for u, _ in res["states"]]), uf=np.stack([f for _, f in res["states"]]), seconds=res["seconds"]) def moment_task(task): clock = time.time() with np.load(Path(task["parts"]) / f"base_{task['seed'] + 1:02d}.npz") as base: where = [int(t) for t in base["moments"]].index(task["T"]) us, uf = base["us"][where], base["uf"][where] rec, dtheta = planned(us, uf, (DRAWS, task["seed"], task["moment"])) length = int(round(task["after"])) z = np.full((6, 2, length + 1), np.nan, np.complex64) logs = np.full((6, int(round(task["after"] / EVERY))), np.nan) if rec["carried"]: z, logs = branches(us, uf, dtheta, task["after"]) z = z.astype(np.complex64) return dict(start=task["seed"] + 1, moment=task["moment"] + 1, T=task["T"], branch_z=z, branch_logs=logs, seconds=time.time() - clock, **{key: np.asarray(value) for key, value in rec.items()}) def work(task): """One task from beginning to end. Its part is written whole or not at all.""" clock = time.time() parts = Path(task["parts"]) final, partial = parts / f"{task['name']}.npz", parts / f"{task['name']}.partial.npz" try: out = base_task(task) if task["kind"] == "base" else moment_task(task) np.savez(partial, **out) os.replace(partial, final) return dict(name=task["name"], ok=True, sha256=sha256(final), seconds=time.time() - clock) except Exception as error: # noqa: BLE001 (a failed task leaves nothing; the next piece runs it again) partial.unlink(missing_ok=True) return dict(name=task["name"], ok=False, error=f"{type(error).__name__}: {error}", seconds=time.time() - clock) def piece(todo, parts, runlog, header): """One piece: the next unfinished tasks of one kind in the fixed order, on WORKERS workers. Parts already logged must be there unchanged; anything else in the folder is from a task that did not finish and is removed. Prints counts and times only.""" parts.mkdir(exist_ok=True) (parts / ".gitignore").write_text("*\n") log = json.loads(runlog.read_text()) if runlog.exists() else dict(header, pieces=[], parts={}) for name, entry in log["parts"].items(): path = parts / f"{name}.npz" if not path.exists() or sha256(path) != entry["sha256"]: sys.exit(f"{path.name} is missing or differs from the run log. Stopping.") for path in parts.iterdir(): if path.name != ".gitignore" and path.name.removesuffix(".npz") not in log["parts"]: path.unlink() left = [task for task in todo if task["name"] not in log["parts"]] if not left: print(f" nothing left to run: all {len(todo)} tasks are done", flush=True) return log kind = left[0]["kind"] now = [task for task in left if task["kind"] == kind][:PIECE[kind]] entry = dict(piece=len(log["pieces"]) + 1, began=stamp(), kind=kind, tasks=[task["name"] for task in now], done=[], failed={}) log["pieces"].append(entry) write_json(runlog, log) clock = time.time() with Pool(min(WORKERS, len(now))) as pool: for res in pool.imap_unordered(work, now): if res["ok"]: log["parts"][res["name"]] = dict(sha256=res["sha256"], piece=entry["piece"], finished=stamp(), seconds=round(res["seconds"], 1)) entry["done"].append(res["name"]) print(f" done: {res['name']} ({res['seconds']:.0f} s)", flush=True) else: entry["failed"][res["name"]] = res["error"] print(f" FAILED: {res['name']} ({res['seconds']:.0f} s): {res['error']}", flush=True) write_json(runlog, log) entry["ended"] = stamp() write_json(runlog, log) done = len(log["parts"]) print(f" piece {entry['piece']}: {len(entry['done'])} of {len(now)} {kind} tasks done in" f" {time.time() - clock:.0f} s; {done} of {len(todo)} tasks done, {len(todo) - done} left", flush=True) return log def assemble(parts, starts, table, out, run, after): """Every part into one file. The saved full states stay behind: the base runs regenerate them exactly.""" base = [] for s in starts: with np.load(parts / f"base_{s:02d}.npz") as p: assert int(p["start"]) == s base.append({key: p[key] for key in ("z", "logs", "trace", "seconds")}) per = len(table) // len(starts) moments = [] for start, m, T in table: with np.load(parts / f"moment_{start:02d}_{m:02d}.npz") as p: assert (int(p["start"]), int(p["moment"]), int(p["T"])) == (start, m, T) moments.append({key: p[key] for key in p.files}) def stacked(key): return np.stack([x[key] for x in moments]).reshape((len(starts), per) + moments[0][key].shape) data = dict(starts=np.array(starts), T=np.array([T for _, _, T in table]).reshape(len(starts), per), z=np.concatenate([b["z"] for b in base]), logs=np.concatenate([b["logs"] for b in base]), trace=np.concatenate([b["trace"] for b in base]), base_seconds=np.array([float(b["seconds"]) for b in base]), moment_seconds=stacked("seconds"), n=N, k=K, dt=DT, settle=nt.SETTLE, run=run, every=EVERY, gap=GAP, seed=nt.LYAP_SEED, draws=DRAWS, after=after) for key in RECORD_KEYS + ("branch_z", "branch_logs"): data[key] = stacked(key) np.savez_compressed(out, **data) with np.load(out) as saved: for key, value in data.items(): same = (np.array_equal(saved[key], value, equal_nan=True) if np.asarray(value).dtype.kind in "fc" else np.array_equal(saved[key], value)) if not same: raise RuntimeError(f"{out.name}: {key} did not save as assembled") return data # --------------------------------------------------------------------------------------------------------------------- # the analysis def analyse(d, lock=None, log=None, dry=False): """Every number of the plan from the assembled data, the decision by cause_rule.py, and the report's lines.""" starts = [int(s) for s in d["starts"]] n, picks = len(starts), rule.draws(len(starts)) carried = np.asarray(d["carried"], bool) # (start, moment) per = carried.shape[1] T = np.asarray(d["T"], int) built = np.asarray(d["built"], bool) # (start, moment, method) logs = np.where(carried[..., None, None], np.asarray(d["branch_logs"], float), np.nan) def rate(x, window): low, high = (int(round(v / EVERY)) for v in window) return x[..., low:high].sum(-1) / (window[1] - window[0]) G = {name: rate(logs, window).reshape(n, per, 2, 3) for name, window in WINDOWS.items()} D = {name: rule.per_start(G[name]) for name in WINDOWS} # (D_up, D_down), each (start, method) d_up, d_down = D["later"] # M1 count_ok = rule.enough_moments(carried) missing = int((~carried).sum()) same, match = np.asarray(d["same"], bool), np.asarray(d["match"], float) gaprel = np.asarray(d["gaprel"], float).reshape(n, per, 2, 3) inband = np.asarray(d["inband"], bool).reshape(n, per, 2, 3, 2) m1_checks = {} for k, method in enumerate(METHODS): b = built[..., k] m1_checks[method] = dict( matched=bool(np.all(same[..., k][b]) and np.all(match[..., k][b] <= MATCH)), gap=bool(np.all(gaprel[..., k, :][b] < KEPT)), bands=bool(np.all(inband[..., k, :, :][carried])), largest_match=float(np.max(match[..., k][b], initial=0.0)), largest_gap=float(np.max(gaprel[..., k, :][b], initial=0.0))) m1 = {method: bool(count_ok and m1_checks[method]["matched"] and m1_checks[method]["gap"] and m1_checks[method]["bands"]) for method in METHODS} m1["pooled"] = m1["dense"] and m1["sparse"] # M2 and the coherence paths paths = np.abs(np.asarray(d["branch_z"]).astype(np.complex128)).mean(-2) # C, (start, moment, branch, time) paths = np.where(carried[..., None, None], paths, np.nan).reshape(n, per, 2, 3, -1) up_down = paths[..., 0, :] - paths[..., 1, :] # (start, moment, method, time) held, at_t = up_down[..., 1:201].mean(-1), up_down[..., 0] phi = dict(pooled=float(held[carried].mean() / at_t[carried].mean()), **{method: float(held[..., k][carried].mean() / at_t[..., k][carried].mean()) for k, method in enumerate(METHODS)}) m2 = {key: bool(value >= FLOOR) for key, value in phi.items()} # the threshold C0 = np.asarray(d["C0"], float) C1 = np.asarray(d["C1"], float).reshape(n, per, 2, 3) size = (np.abs(C1[..., 0] - C0[..., None]) + np.abs(C1[..., 1] - C0[..., None])) / 2 # (start, moment, method) sizes = dict(pooled=float(size[carried].mean()), **{method: float(size[..., k][carried].mean()) for k, method in enumerate(METHODS)}) delta = {key: HALF * BENCHMARK * value for key, value in sizes.items()} # the decision (out, numbers), each = rule.decide(d_up, d_down, delta, picks, m1, m2) decision = str(np.asarray(out).item()) readings = {method: str(np.asarray(each[method][0]).item()) for method in METHODS} def summary(values): values = np.asarray(values, float) bar = float(rule.error_bar(values, picks)) mean = float(values.mean()) return dict(mean=mean, error_bar=bar, textbook=float(values.std(ddof=1) / np.sqrt(len(values))), resampled=float(values[picks].mean(-1).std()), bars=mean / bar if bar > 0 else float("nan"), above=int((values > 0).sum()), below=int((values < 0).sum()), starts=values.tolist()) def four(up, down): out = {} for key in ("pooled",) + METHODS: u = up.mean(-1) if key == "pooled" else up[:, METHODS.index(key)] w = down.mean(-1) if key == "pooled" else down[:, METHODS.index(key)] out[key] = dict(S=summary((w - u) / 2), up=summary(u), down=summary(w), P=summary((u + w) / 2)) return out primary = four(d_up, d_down) assert abs(primary["pooled"]["S"]["mean"] - float(numbers["s"])) < 1e-15 assert abs(primary["pooled"]["S"]["error_bar"] - float(numbers["s_bar"])) < 1e-15 other = {name: four(*D[name]) for name in ("early", "whole")} s_k = (d_down - d_up) / 2 methods_differ = summary(s_k[:, 0] - s_k[:, 1]) def clear_reading(s): if s["mean"] > rule.RULE * s["error_bar"]: return "speeds later growth" if s["mean"] < -rule.RULE * s["error_bar"]: return "slows later growth" return "no clear effect" # sham against none: the base run's own growth over the same window base_logs = np.asarray(d["logs"], float) first = T // int(round(EVERY)) def untouched(window): low, high = (int(round(v / EVERY)) for v in window) values = np.array([[base_logs[i, first[i, m] + low:first[i, m] + high].sum() for m in range(per)] for i in range(n)]) / (window[1] - window[0]) return np.where(carried, values, np.nan) none = {name: untouched(window) for name, window in WINDOWS.items()} sham_none = {name: rule.moment_average(G[name][..., 2] - none[name][..., None]) for name in WINDOWS} sham_vs_none = {name: {method: summary(sham_none[name][:, k]) for k, method in enumerate(METHODS)} for name in WINDOWS} # the causal slope s_pool = primary["pooled"]["S"] slope = dict(value=-s_pool["mean"] / sizes["pooled"], error_bar=s_pool["error_bar"] / sizes["pooled"]) slope.update(scaled=slope["value"] * NATURAL / phi["pooled"], scaled_bar=slope["error_bar"] * NATURAL / phi["pooled"]) def start_mean(x): # over each start's moments carried out kept = ~np.isnan(x) with np.errstate(invalid="ignore", divide="ignore"): return np.where(kept, x, 0.0).sum(1) / kept.sum(1) vs_sham = {arm: paths[..., a, :] - paths[..., 2, :] for a, arm in ((0, "up"), (1, "down"))} coherence_paths = {arm: start_mean(v).mean(0) for arm, v in vs_sham.items()} # (method, time) phi_sham = {method: {arm: float(v[..., k, 1:201].mean(-1)[carried].mean() / v[..., k, 0][carried].mean()) for arm, v in vs_sham.items()} for k, method in enumerate(METHODS)} late = up_down[..., 201:501].mean(-1) left_over = dict(pooled=float(late[carried].mean() / at_t[carried].mean()), **{method: float(late[..., k][carried].mean() / at_t[..., k][carried].mean()) for k, method in enumerate(METHODS)}) cumulative = np.cumsum(logs, -1).reshape(n, per, 2, 3, -1) growth_paths = {arm: start_mean(cumulative[..., a, :] - cumulative[..., 2, :]).mean(0) for a, arm in ((0, "up"), (1, "down"))} # (method, interval) # what the changes did to the state (immediate; reported, not used) R0, R1 = np.asarray(d["R0"], float), np.asarray(d["R1"], float) psi0, psi1 = np.asarray(d["psi0"], float), np.asarray(d["psi1"], float) C1flat = np.asarray(d["C1"], float) quantities = dict( ratio_A=R1[..., 0] / R0[..., None, 0], ratio_B=R1[..., 1] / R0[..., None, 1], dC=C1flat - C0[..., None], relC=(C1flat - C0[..., None]) / C0[..., None], dpsi_A=wrap(psi1[..., 0] - psi0[..., None, 0]), dpsi_B=wrap(psi1[..., 1] - psi0[..., None, 1]), dab=wrap(np.asarray(d["ab1"], float) - np.asarray(d["ab0"], float)[..., None])) for order in (2, 3, 4): before, after = np.asarray(d[f"z{order}_0"], float), np.asarray(d[f"z{order}_1"], float) for g, group in enumerate(GROUPS): quantities[f"dz{order}_{group}"] = after[..., g] - before[..., None, g] state = {} for name, x in quantities.items(): x = np.where(carried[..., None], x, np.nan).reshape(n, per, 2, 3) state[name] = {} for k, method in enumerate(METHODS): entry = {} for a, arm in enumerate(ARMS): values = x[..., k, a][carried] entry[arm] = dict(summary(start_mean(x[..., k, a][..., None])[:, 0]), smallest=float(values.min()), largest=float(values.max())) for a, arm in ((0, "up"), (1, "down")): entry[f"{arm}_minus_sham"] = summary(start_mean((x[..., k, a] - x[..., k, 2])[..., None])[:, 0]) state[name][method] = entry levels = {} for order in (2, 3, 4): before = np.asarray(d[f"z{order}_0"], float)[carried] # (moment carried, group) after = np.asarray(d[f"z{order}_1"], float)[carried].reshape(-1, 2, 3, 2) levels[f"z{order}"] = dict(before=before.mean(0).tolist(), after={method: {arm: after[:, k, a].mean(0).tolist() for a, arm in enumerate(ARMS)} for k, method in enumerate(METHODS)}) moves = {} for name in ("sumsq", "moved", "largest", "rms"): x = np.asarray(d[name], float)[carried].reshape(-1, 2, 3) moves[name] = {method: {arm: dict(mean=float(x[:, k, a].mean()), largest=float(x[:, k, a].max())) for a, arm in enumerate(ARMS)} for k, method in enumerate(METHODS)} not_carried = [dict(start=starts[i], moment=int(m) + 1, T=int(T[i, m]), reason=str(d["reason"][i, m])) for i, m in zip(*np.nonzero(~carried))] # how much the moments of a start go together def icc(x): kept = ~np.isnan(x) counts = kept.sum(1) total = counts.sum() with np.errstate(invalid="ignore", divide="ignore"): means = np.where(kept, x, 0.0).sum(1) / counts grand = x[kept].mean() between = (counts * (means - grand) ** 2).sum() / (n - 1) within = np.square(np.where(kept, x - means[:, None], 0.0)).sum() / (total - n) n0 = (total - (counts ** 2).sum() / total) / (n - 1) value = (between - within) / (between + (n0 - 1) * within) return dict(icc=float(value), design_effect=float(1 + (n0 - 1) * value), n0=float(n0)) g_later = G["later"] du, dd = g_later[..., 0] - g_later[..., 2], g_later[..., 1] - g_later[..., 2] sm = (dd - du) / 2 iccs = dict(pooled=dict(S=icc(sm.mean(-1)), up=icc(du.mean(-1)), down=icc(dd.mean(-1))), **{method: dict(S=icc(sm[..., k]), up=icc(du[..., k]), down=icc(dd[..., k])) for k, method in enumerate(METHODS)}) # the base runs z_base = np.asarray(d["z"]) measured = cg.measures(z_base, base_logs, np.random.default_rng(cg.SEED)) test16 = json.loads(TEST16.read_text()) combined = float(np.hypot(measured["long_rate_se"], test16["rate_se"])) long_rate = dict(new=measured["long_rate"], new_se=measured["long_rate_se"], test16=test16["rate"], test16_se=test16["rate_se"], agree=bool(abs(measured["long_rate"] - test16["rate"]) <= 3 * combined)) checked = egc.checked(z_base, base_logs, np.random.default_rng(cg.SEED)) main_values = np.asarray(checked["main"]["starts"], float) main_textbook = float(main_values.std(ddof=1) / np.sqrt(len(main_values))) main_bar = max(float(checked["main"]["se"]), main_textbook) test17 = dict(mean=float(checked["main"]["mean"]), resampled=float(checked["main"]["se"]), textbook=main_textbook, error_bar=main_bar, bars=float(checked["main"]["mean"]) / main_bar, below=int((main_values < 0).sum())) # why, for the readings that need it def m1_failures(method_list): out = [] if rule.MOST_MISSING * missing > carried.size: out.append(f"{missing} of {carried.size} moments not carried out (at most" f" {carried.size // rule.MOST_MISSING} allowed)") if int(carried.sum(1).min()) < rule.FEWEST_KEPT: out.append(f"a start kept only {int(carried.sum(1).min())} of its {per} moments") for method in method_list: for key, words in (("matched", "arms not matched"), ("gap", "gap changed"), ("bands", "a change outside its band")): if not m1_checks[method][key]: out.append(f"{method}: {words}") return out def why(reading, nums, method=None): key = method or "pooled" if reading == "INCONCLUSIVE": if not m1[key]: return "M1 failed: " + "; ".join(m1_failures(METHODS if method is None else (method,))) upper = float(nums["s"]) + rule.SIDE * float(nums["s_bar"]) if upper >= delta[key]: return (f"an effect of half the benchmark cannot be ruled out (S + 2 error bars = {upper:+.5f}," f" at least delta = {delta[key]:.5f})") return f"the change faded too fast (M2: Phi = {phi[key]:.3f}, below {FLOOR})" if reading == "PARTIAL": failed = [] if not float(nums["up"]) < -rule.SIDE * float(nums["up_bar"]): failed.append("up is not clearly below its sham") if not float(nums["down"]) > rule.SIDE * float(nums["down_bar"]): failed.append("down is not clearly above its sham") if method is None: for k, name in enumerate(METHODS): s = primary[name]["S"] if not s["mean"] > rule.SIDE * s["error_bar"]: failed.append(f"the {name} method's S is not clearly above 0") return "; ".join(failed) return "" result = dict( dry=dry, question=("If coherence is raised or lowered at a chosen moment, does the growth of a tiny gap 200 to" " 500 time units later change in the direction the correlation predicts, compared with a" " matched sham?"), decision=decision, meaning=MEANING[decision], why=why(decision, numbers), readings={method: dict(reading=readings[method], why=why(readings[method], each[method][1], method), numbers=plain(each[method][1])) for method in METHODS}, numbers=plain(numbers), m1=m1, m1_checks=m1_checks, moments_carried=int(carried.sum()), moments=int(carried.size), fewest_in_a_start=int(carried.sum(1).min()), m2=m2, phi=phi, sizes=sizes, delta=delta, primary=primary, other_windows=other, methods_differ=methods_differ, P_readings={key: clear_reading(primary[key]["P"]) for key in ("pooled",) + METHODS}, sham_vs_none=sham_vs_none, slope=slope, phi_vs_sham=phi_sham, left_over=left_over, coherence_paths={arm: v.tolist() for arm, v in coherence_paths.items()}, growth_paths={arm: v.tolist() for arm, v in growth_paths.items()}, state=state, levels=levels, moves=moves, not_carried=not_carried, icc=iccs, long_rate=long_rate, test17=test17, starts=starts, per_start=[dict(start=starts[i], S=primary["pooled"]["S"]["starts"][i], up=primary["pooled"]["up"]["starts"][i], down=primary["pooled"]["down"]["starts"][i], S_dense=primary["dense"]["S"]["starts"][i], S_sparse=primary["sparse"]["S"]["starts"][i], carried=int(carried[i].sum())) for i in range(n)], lock=None if lock is None else dict(locked=lock["locked"], plan_sha256=lock.get("plan_sha256")), runlog=None if log is None else dict(pieces=len(log["pieces"]), parts=len(log["parts"]), began=log["pieces"][0]["began"], ended=log["pieces"][-1].get("ended")), deviations=[]) return plain(result), report_lines(result) def report_lines(r): """The report, in words and numbers.""" def pm(s, digits=5): return (f"{s['mean']:+.{digits}f} (error bar {s['error_bar']:.{digits}f}; {s['bars']:+.1f} error bars;" f" above 0 in {s['above']} of {len(s['starts'])})") n = len(r["starts"]) lines = [f"Cause test: October 12, 2026. N = {N:,} per group, K = {K:.2f}, eps = {EPS:g}; starts" f" {r['starts'][0]} to {r['starts'][-1]} ({n}), {MOMENTS} moments each."] if r["dry"]: lines.append("DRY RUN ON MADE-UP NUMBERS: NOT A RESULT.") if r["lock"]: lines.append(f"Pre-registered: {PLAN.name} (SHA-256 {r['lock']['plan_sha256']}), locked" f" {shown(r['lock']['locked'])}.") if r["runlog"]: lines.append(f"Run: {r['runlog']['parts']} tasks in {r['runlog']['pieces']} pieces," f" {shown(r['runlog']['began'])} to {shown(r['runlog']['ended'])}.") lines += ["A toy model (two groups of phase oscillators): this shows structure in the model, not proof beyond it.", "", f"Question: {r['question']}", "", f"DECISION: {r['decision']}: {r['meaning']}"] if r["why"]: lines.append(f" Why: {r['why']}.") lines += ["", "Manipulation checks", f" moments carried out: {r['moments_carried']} of {r['moments']} (fewest in a start:" f" {r['fewest_in_a_start']} of {MOMENTS}); at most {r['moments'] // rule.MOST_MISSING} may be missing"] for method in METHODS: c = r["m1_checks"][method] lines.append(f" {method}: every branch carried out inside its band {c['bands']}; same oscillators and sizes" f" {c['matched']} (largest relative difference in sum of squared moves {c['largest_match']:.1e});" f" gap as it was {c['gap']} (largest {c['largest_gap']:.1e}); M1 {'holds' if r['m1'][method] else 'FAILS'}") lines.append(f" M2: Phi = {r['phi']['pooled']:.3f} pooled (dense {r['phi']['dense']:.3f}, sparse" f" {r['phi']['sparse']:.3f}); floor {FLOOR}: {'holds' if r['m2']['pooled'] else 'FAILS'}") p = r["primary"]["pooled"] lines += ["", "Primary measure: later growth over (T + 200, T + 500], each arm minus its own sham; each start's" f" moments averaged, then the two methods, then the {n} starts (error bar: the larger of resampled and" f" SD/sqrt({n}))", f" S = (D_down - D_up)/2 {pm(p['S'])}", f" D_up (up - sham) {pm(p['up'])}", f" D_down (down - sham) {pm(p['down'])}", f" delta = 1/2 x {BENCHMARK} x |dC| = {r['delta']['pooled']:.5f} (|dC| = {r['sizes']['pooled']:.5f})", f" clear (S > 3 error bars): {bool(r['numbers']['clear'])}; sides: {bool(r['numbers']['sides'])};" f" methods agree: {bool(r['numbers']['methods_agree'])}; excluded (S + 2 error bars < delta):" f" {p['S']['mean'] + rule.SIDE * p['S']['error_bar'] < r['delta']['pooled']}", "", "Each method's own reading (beside the decision, not the decision)"] for method in METHODS: q, rd = r["primary"][method], r["readings"][method] lines.append(f" {method}: {rd['reading']}" + (f" ({rd['why']})" if rd["why"] else "")) lines.append(f" S {pm(q['S'])}; D_up {pm(q['up'])}; D_down {pm(q['down'])}; delta {r['delta'][method]:.5f};" f" M1 {'holds' if r['m1'][method] else 'fails'}; Phi {r['phi'][method]:.3f}") lines.append(f" dense S - sparse S: {pm(r['methods_differ'])}") lines += ["", "Also reported (fixed before the run, not used for the decision)", " P = (D_up + D_down)/2, a change of this size either way: " + "; ".join( f"{key} {pm(r['primary'][key]['P'])}: {r['P_readings'][key]}" for key in ("pooled",) + METHODS), " sham against none (G_sham - G_none): " + "; ".join( f"{method} {pm(r['sham_vs_none']['later'][method])}" for method in METHODS), f" causal slope -S/|dC|: {r['slope']['value']:+.4f} (error bar {r['slope']['error_bar']:.4f}); benchmark" f" -{BENCHMARK}; scaled by {NATURAL}/Phi: {r['slope']['scaled']:+.4f} ({r['slope']['scaled_bar']:.4f})"] for name, words in (("early", "(T, T + 200]"), ("whole", "(T, T + 500]")): o = r["other_windows"][name] lines.append(f" growth over {words}: " + "; ".join( f"{key} S {o[key]['S']['mean']:+.5f} ({o[key]['S']['error_bar']:.5f}), D_up {o[key]['up']['mean']:+.5f}" f" ({o[key]['up']['error_bar']:.5f}), D_down {o[key]['down']['mean']:+.5f} ({o[key]['down']['error_bar']:.5f})" for key in ("pooled",) + METHODS)) lines.append(" coherence paths: Phi against sham: " + "; ".join( f"{method} up {r['phi_vs_sham'][method]['up']:.3f}, down {r['phi_vs_sham'][method]['down']:.3f}" for method in METHODS) + "; share of the up - down difference left over (T + 200, T + 500]: " + ", ".join( f"{key} {r['left_over'][key]:.3f}" for key in ("pooled",) + METHODS)) lines += ["", "What the changes did to the state, just after the change (moments carried out; average over the" " starts, with its error bar, and the range over moments). Reported, not used for the decision."] names = (("ratio_A", "R_A after / before", 4), ("ratio_B", "R_B after / before", 4), ("dC", "change in C", 5), ("relC", "change in C / C", 4), ("dpsi_A", "change in psi_A", 5), ("dpsi_B", "change in psi_B", 5), ("dab", "change in psi_A - psi_B", 5), ("dz2_A", "change in |z2| A", 6), ("dz2_B", "change in |z2| B", 6), ("dz3_A", "change in |z3| A", 6), ("dz3_B", "change in |z3| B", 6), ("dz4_A", "change in |z4| A", 6), ("dz4_B", "change in |z4| B", 6)) for method in METHODS: lines.append(f" {method}") for key, words, digits in names: e = r["state"][key][method] cells = "; ".join(f"{arm} {e[arm]['mean']:+.{digits}f} [{e[arm]['smallest']:+.{digits}f}," f" {e[arm]['largest']:+.{digits}f}]" for arm in ARMS) diffs = "; ".join(f"{arm} - sham {e[f'{arm}_minus_sham']['mean']:+.{digits}f}" f" ({e[f'{arm}_minus_sham']['error_bar']:.{digits}f})" for arm in ("up", "down")) lines.append(f" {words}: {cells} | {diffs}") for key in ("z2", "z3", "z4"): lv = r["levels"][key] lines.append(f" |{key}| before (A, B): {lv['before'][0]:.5f}, {lv['before'][1]:.5f}; after: " + "; ".join( f"{method} {arm} {lv['after'][method][arm][0]:.5f}, {lv['after'][method][arm][1]:.5f}" for method in METHODS for arm in ARMS)) lines += ["", "The moves themselves (matched by construction; kept apart from their effects above)"] for method in METHODS: mv = r["moves"] lines.append(f" {method}: " + "; ".join( f"{arm} sum of squares {mv['sumsq'][method][arm]['mean']:.4g} (largest {mv['sumsq'][method][arm]['largest']:.4g})," f" moved {mv['moved'][method][arm]['mean']:.0f}, largest move {mv['largest'][method][arm]['mean']:.4f}" f" (largest {mv['largest'][method][arm]['largest']:.4f}), rms {mv['rms'][method][arm]['mean']:.4f}" for arm in ARMS)) lines += ["", f"Moments not carried out: {len(r['not_carried'])}"] lines += [f" start {x['start']}, moment {x['moment']} (T = {x['T']:,}): {x['reason']}" for x in r["not_carried"]] lines += ["", "How much the moments of a start go together (intraclass correlation; design effect)"] for key in ("pooled",) + METHODS: lines.append(f" {key}: " + "; ".join(f"{name} {r['icc'][key][name]['icc']:+.3f}" f" ({r['icc'][key][name]['design_effect']:.2f})" for name in ("S", "up", "down"))) lr, t17 = r["long_rate"], r["test17"] lines += ["", f"From the base runs (the {n} untouched new starts)", f" long-run rate, t = 3,000 to 6,000: {lr['new']:+.4f} ({lr['new_se']:.4f}); Test 16: {lr['test16']:+.4f}" f" ({lr['test16_se']:.4f}): {'agree' if lr['agree'] else 'DIFFER'} (within 3 combined error bars)", f" Test 17's conditional coefficient on these starts (its own code; not part of this decision):" f" {t17['mean']:+.4f} (error bar {t17['error_bar']:.4f}; {t17['bars']:+.1f} error bars; below 0 in" f" {t17['below']} of {n})", "", "Per start: S, D_up, D_down (pooled); S dense, S sparse; moments carried out"] lines += [f" {x['start']:2d}: {x['S']:+.5f} {x['up']:+.5f} {x['down']:+.5f} | {x['S_dense']:+.5f}" f" {x['S_sparse']:+.5f} | {x['carried']}" for x in r["per_start"]] lines += ["", "Deviations from the plan: " + ("; ".join(r["deviations"]) if r["deviations"] else "none")] return lines def figure(r, path): """Four panels: S start by start; the arms against their shams; the coherence paths; the growth paths.""" import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt fig, axes = plt.subplots(2, 2, figsize=(13, 9.6)) fig.subplots_adjust(left=0.07, right=0.98, top=0.88, bottom=0.15, hspace=0.42, wspace=0.22) title = f"Cause test: {r['decision']}" + (" (DRY RUN ON MADE-UP NUMBERS)" if r["dry"] else "") fig.text(0.06, 0.955, title, fontsize=15, weight="bold") fig.text(0.06, 0.925, f"Starts {r['starts'][0]} to {r['starts'][-1]}; later growth over (T + 200, T + 500], each" " arm against its own sham; N = 4,000 per group, K = 1.90", fontsize=10, color="0.3") colors = dict(pooled="#1f3b73", dense="#c0392b", sparse="#2e8b57") ax = axes[0, 0] x = np.array(r["starts"]) for key, marker, size in (("dense", "^", 18), ("sparse", "v", 18), ("pooled", "o", 30)): ax.scatter(x, r["primary"][key]["S"]["starts"], s=size, marker=marker, color=colors[key], label=key, alpha=1.0 if key == "pooled" else 0.55, zorder=3 if key == "pooled" else 2) s = r["primary"]["pooled"]["S"] ax.axhspan(s["mean"] - s["error_bar"], s["mean"] + s["error_bar"], color=colors["pooled"], alpha=0.12) ax.axhline(s["mean"], color=colors["pooled"], lw=1.2) ax.axhline(r["delta"]["pooled"], color="#d68910", lw=1.2, ls="--", label=f"delta = {r['delta']['pooled']:.5f}") ax.axhline(0, color="0.4", lw=0.8) ax.set_title("S = (D_down - D_up)/2, start by start", fontsize=11) ax.set_xlabel("start") ax.legend(fontsize=8, ncol=2, loc="best") ax = axes[0, 1] names = (("up", "D_up"), ("down", "D_down"), ("S", "S"), ("P", "P")) for j, key in enumerate(("pooled",) + METHODS): for i, (name, label) in enumerate(names): e = r["primary"][key][name] ax.errorbar(i + (j - 1) * 0.22, e["mean"], yerr=e["error_bar"], fmt="o", color=colors[key], capsize=3, label=key if i == 0 else None) ax.axhline(0, color="0.4", lw=0.8) ax.set_xticks(range(len(names)), [label for _, label in names]) ax.set_title("The arms against their shams (mean and error bar over the starts)", fontsize=11) ax.legend(fontsize=8) ax = axes[1, 0] for k, method in enumerate(METHODS): for arm, style in (("up", "-"), ("down", "--")): path_ = np.array(r["coherence_paths"][arm][k]) ax.plot(np.arange(len(path_)), path_, style, color=colors[method], lw=1.3, label=f"{method} {arm} - sham") ax.axhline(0, color="0.4", lw=0.8) ax.axvspan(200, 500, color="0.85", alpha=0.4) ax.set_title("Coherence C after the change, against the sham", fontsize=11) ax.set_xlabel("time after T") ax.legend(fontsize=8, ncol=2) ax = axes[1, 1] for k, method in enumerate(METHODS): for arm, style in (("up", "-"), ("down", "--")): path_ = np.array(r["growth_paths"][arm][k]) ax.plot(EVERY * np.arange(1, len(path_) + 1), path_, style, color=colors[method], lw=1.3, label=f"{method} {arm} - sham") ax.axhline(0, color="0.4", lw=0.8) ax.axvspan(200, 500, color="0.85", alpha=0.4) ax.set_title("Accumulated log growth of the gap, against the sham", fontsize=11) ax.set_xlabel("time after T (shaded: the outcome's window)") ax.legend(fontsize=8, ncol=2) fig.text(0.06, 0.06, f"Decision: {r['decision']}. " + (r["why"] + "." if r["why"] else ""), fontsize=10, color="0.2") fig.text(0.06, 0.03, "Error bars from the starts (the larger of resampled and SD/sqrt(n)). Fixed and locked before" " the run. A toy model: this shows structure in the model, not proof beyond it.", fontsize=9, color="0.35") path.parent.mkdir(exist_ok=True) fig.savefig(path, dpi=115) plt.close(fig) # --------------------------------------------------------------------------------------------------------------------- # the checks before the lock def checks(): """The checks before the lock, on old start 1 only. Returns every result; writes nothing.""" out = dict(page=sha256(PLAN_FILES["nudge_test"][0]) == PAGE_SHA, files={name: bool(path.exists() and sha256(path) == sha) for name, (path, sha) in PLAN_FILES.items()}) res = base_run(0, keep=(CHECK_AT,), run=CHECK_RUN) old = saved_run(0) samples, resets = int(CHECK_RUN) + 1, int(round(CHECK_RUN / EVERY)) out["base"] = dict(settling=bool(np.array_equal(res["trace"][0], old["trace"])), coherence=bool(np.array_equal(res["z"].astype(np.complex64)[0], old["z"][:, :samples])), gap=bool(np.array_equal(res["logs"][0], old["logs"][:resets])), seconds=res["seconds"]) out["base"]["repeats"] = out["base"]["settling"] and out["base"]["coherence"] and out["base"]["gap"] us, uf = res["states"][0] clock = time.time() zb, lb = branches(us, uf, np.zeros((6, 2, us.shape[-1])), CHECK_RUN - CHECK_AT) spent = time.time() - clock ref_z, ref_logs = res["z"][0][:, CHECK_AT:], res["logs"][0][int(CHECK_AT / EVERY):] out["unchanged"] = dict(coherence=float(np.abs(zb - ref_z[None]).max()), gap=float(np.abs(lb - ref_logs[None]).max()), seconds=spent, per_task=spent * AFTER / (CHECK_RUN - CHECK_AT)) out["unchanged"]["holds"] = out["unchanged"]["coherence"] <= 1e-9 and out["unchanged"]["gap"] <= 1e-6 clock = time.time() rec, dtheta = planned(us, uf, (DRAWS, 0, 0)) out["planned"] = dict(carried=rec["carried"], reason=rec["reason"], R=rec["R0"], ratio=rec["ratio"], amp=rec["amp"], status=rec["status"], same=rec["same"], match=rec["match"], gaprel=float(np.nanmax(rec["gaprel"])), largest=rec["largest"], seconds=time.time() - clock) with open(dn.MOMENT_TABLE, newline="") as handle: rows = [tuple(int(x) for x in row) for row in list(csv.reader(handle))[1:]] out["table"] = rows == dn.moment_table() out["new_data"] = [path.name for path in (SAVED, RUNLOG, PARTS) if path.exists()] return out def show_checks(c): print("Checks before the lock (old start 1 only; no new start is run)") print(f" nudge_test.py is the copy the page carries: {c['page']}") print(" files in the plan's table: " + ", ".join(f"{name} {ok}" for name, ok in c["files"].items())) b = c["base"] print(f" the base-run code against the saved long run, start 1 for {CHECK_RUN:.0f} time units: settling" f" {b['settling']}, coherence {b['coherence']}, gap {b['gap']} ({b['seconds']:.1f} s)") u = c["unchanged"] print(f" six branches with no change, run as in a moment task, from t = {CHECK_AT} to {CHECK_RUN:.0f}" f" against the base run: coherence within {u['coherence']:.1e} (limit 1e-9), log of the gap within" f" {u['gap']:.1e} (limit 1e-6): {u['holds']} ({u['seconds']:.1f} s for {CHECK_RUN - CHECK_AT:.0f} time" f" units; about {u['per_task']:.0f} s of branches per moment task on one worker)") p = c["planned"] ratio = np.asarray(p["ratio"]) print(f" the planned changes at t = {CHECK_AT} (R_A {p['R'][0]:.4f}, R_B {p['R'][1]:.4f}): carried out" f" {p['carried']}" + (f" ({p['reason']})" if p["reason"] else "")) for k, method in enumerate(METHODS): print(f" {method}: amplitudes {p['amp'][k][0]:.5f}, {p['amp'][k][1]:.5f} (status {p['status'][k].tolist()});" " R after / before (A, B): " + "; ".join(f"{arm} {ratio[3 * k + a, 0]:.4f}, {ratio[3 * k + a, 1]:.4f}" for a, arm in enumerate(ARMS)) + f"; same {p['same'][k]}; sums of squares differ by {p['match'][k]:.1e}; largest move" f" {p['largest'][3 * k]:.4f}") print(f" gap after / before: largest relative change {p['gaprel']:.1e} (limit 1e-6) ({p['seconds']:.2f} s)") print(f" the regenerated moment table equals cause_moments.csv: {c['table']}") print(f" new data: {', '.join(c['new_data']) if c['new_data'] else 'none'}") def check(): frozen() print(f"Settings match the plan; plan SHA-256 {sha256(PLAN)}" + (" (the approved version)" if sha256(PLAN) == PLAN_SHA else " (not yet the approved version)")) show_checks(checks()) def target_sweep(s): """One old start: its base run to t = 5,495, with the planned changes worked out (not run on) at every multiple of 5 from t = 2,000, from the draws (20261012, seed, 100 + index, use), which the run never uses.""" rows = [] def visit(t, us, uf): if TARGET_SPAN[0] <= t <= TARGET_SPAN[1]: rec, _ = planned(us, uf, (DRAWS, s, TARGET_KEY + (t - TARGET_SPAN[0]) // int(EVERY))) rows.append(rec) res = base_run(s, run=float(TARGET_SPAN[1]), visit=visit) saved = saved_run(s) samples, resets = TARGET_SPAN[1] + 1, TARGET_SPAN[1] // int(EVERY) repeats = dict(coherence=bool(np.array_equal(res["z"].astype(np.complex64)[0], saved["z"][:, :samples])), gap=bool(np.array_equal(res["logs"][0], saved["logs"][:resets]))) records = {key: np.stack([np.asarray(row[key]) for row in rows]) for key in RECORD_KEYS} return dict(start=s + 1, records=records, repeats=repeats, seconds=res["seconds"]) def target_summary(found, seconds): rec = {key: np.concatenate([f["records"][key] for f in found]) for key in RECORD_KEYS} count = len(rec["carried"]) R0, status, built = rec["R0"], rec["status"], rec["built"].astype(bool) ratio, inband, carried = rec["ratio"], rec["inband"].astype(bool), rec["carried"].astype(bool) centiles = {f"{q:g}": float(np.percentile(R0, q)) for q in (1, 50, 99)} centiles["largest"] = float(R0.max()) targets = {method: dict(solved=int((status[:, k] == 0).sum()), infeasible=int((status[:, k] == 1).sum()), failures=int((status[:, k] == 2).sum()), infeasible_R=sorted(float(x) for x in R0[status[:, k] == 1])) for k, method in enumerate(METHODS)} achieved = {} for k, method in enumerate(METHODS): b = built[:, k] for a, arm in enumerate(ARMS): x = ratio[b, 3 * k + a] achieved[f"{method} {arm}"] = dict(smallest=float(x.min()), mean=float(x.mean()), largest=float(x.max()), outside=int((~inband[b, 3 * k + a]).sum())) kinds = dict(infeasible=int((status == 1).any((1, 2)).sum()), failure=int((status == 2).any((1, 2)).sum()), band=int((built.all(1) & ~inband.all((1, 2))).sum()), code=int(rec["code_error"].sum())) examples = sorted({str(x) for x in rec["reason"][~carried]})[:12] matching = dict(same=bool(np.all(rec["same"][built])), largest=float(np.max(rec["match"][built], initial=0.0))) gap = float(np.nanmax(rec["gaprel"])) amp = rec["amp"] a_over_R = {method: dict(median=float(np.nanmedian(amp[:, k] / R0)), largest=float(np.nanmax(amp[:, k] / R0))) for k, method in enumerate(METHODS)} moves = {method: {f"{q:g}": float(np.nanpercentile(rec["largest"][:, 3 * k], q)) for q in (50, 99, 100)} for k, method in enumerate(METHODS)} quantities = dict(ratio_A=rec["R1"][..., 0] / R0[:, None, 0], ratio_B=rec["R1"][..., 1] / R0[:, None, 1], relC=(rec["C1"] - rec["C0"][:, None]) / rec["C0"][:, None], dpsi_A=wrap(rec["psi1"][..., 0] - rec["psi0"][:, None, 0]), dpsi_B=wrap(rec["psi1"][..., 1] - rec["psi0"][:, None, 1]), dab=wrap(rec["ab1"] - rec["ab0"][:, None])) for order in (2, 3, 4): for g, group in enumerate(GROUPS): quantities[f"dz{order}_{group}"] = rec[f"z{order}_1"][..., g] - rec[f"z{order}_0"][:, None, g] state = {} for name, x in quantities.items(): x = x[carried] state[name] = {} for k, method in enumerate(METHODS): entry = {arm: dict(mean=float(x[:, 3 * k + a].mean()), sd=float(x[:, 3 * k + a].std())) for a, arm in enumerate(ARMS)} for a, arm in ((0, "up"), (1, "down")): diff = x[:, 3 * k + a] - x[:, 3 * k + 2] entry[f"{arm}_minus_sham"] = dict(mean=float(diff.mean()), sd=float(diff.std())) state[name][method] = entry result = dict( written=stamp(), states=count, starts=[f["start"] for f in found], span=list(TARGET_SPAN), key=TARGET_KEY, repeats={str(f["start"]): f["repeats"] for f in found}, seconds=round(seconds, 1), base_seconds={str(f["start"]): round(f["seconds"], 1) for f in found}, coherence=centiles, targets=targets, achieved=achieved, not_carried=int((~carried).sum()), kinds=kinds, examples=examples, matching=matching, gap=gap, a_over_R=a_over_R, largest_moves=moves, state=state) exact = sum(f["repeats"]["coherence"] and f["repeats"]["gap"] for f in found) lines = [f"Target check (written {shown(result['written'])}, before the lock): the planned changes at every" f" multiple of 5 from t = {TARGET_SPAN[0]:,} to {TARGET_SPAN[1]:,} in old starts {found[0]['start']} to" f" {found[-1]['start']} ({count:,} states, each just after the gap was put back, as at the run's moments).", "Immediate changes only: no branch was run on and no growth was looked at. Draws from (20261012, seed," f" {TARGET_KEY} + state, use), which the run never uses.", "", f"Base runs: the code repeated the saved runs exactly (coherence and gap, every value to t =" f" {TARGET_SPAN[1]:,}) in {exact} of {len(found)} starts.", f"Each group's coherence R_g at these states: 1% {centiles['1']:.4f}, median {centiles['50']:.4f}," f" 99% {centiles['99']:.4f}, largest {centiles['largest']:.4f}.", "", f"Targets (each method: {2 * count:,} group targets)"] for method in METHODS: t = targets[method] lines.append(f" {method}: feasible and solved {t['solved']:,}; infeasible {t['infeasible']}; solver failures" f" {t['failures']}" + (f" (infeasible at R = {', '.join(f'{x:.4f}' for x in t['infeasible_R'][:8])})" if t["infeasible_R"] else "")) lines.append("Achieved R_g after / R_g before, both groups (smallest, average, largest; band; outside the band)") for name, x in achieved.items(): arm = name.split()[1] lines.append(f" {name:12s} {x['smallest']:.4f} {x['mean']:.4f} {x['largest']:.4f} {BANDS[arm][0]:.2f} to" f" {BANDS[arm][1]:.2f} {x['outside']}") lines += [f"Moments that would not be carried out: {result['not_carried']} of {count:,}" f" ({100 * result['not_carried'] / count:.2f}%); the run allows at most 21 of 420 (5%). By kind:" f" infeasible {kinds['infeasible']}, solver failure {kinds['failure']}, outside a band {kinds['band']}," f" code check {kinds['code']}"] lines += [f" e.g. {x}" for x in examples] lines += [f"Matching: the same oscillators and the same size of move in all three arms at every state:" f" {matching['same']}; largest relative difference in the sum of squared moves: {matching['largest']:.1e}" f" (limit 1e-9)", f"Gap kept: largest |d_after - d_before| / |d_before|: {gap:.1e} (limit 1e-6)", "Amplitude / R_g: " + "; ".join(f"{m} median {a_over_R[m]['median']:.3f}, largest {a_over_R[m]['largest']:.3f}" for m in METHODS), "Largest move |dtheta| (radians): " + "; ".join( f"{m} median {moves[m]['50']:.4f}, 99% {moves[m]['99']:.4f}, largest {moves[m]['100']:.4f}" for m in METHODS), "", "What the changes did to the state (average over the states carried out, with the SD over states;" " reported, not used for any decision)"] names = (("ratio_A", "R_A after / before", 4), ("ratio_B", "R_B after / before", 4), ("relC", "change in C / C", 4), ("dpsi_A", "change in psi_A", 5), ("dpsi_B", "change in psi_B", 5), ("dab", "change in psi_A - psi_B", 5), ("dz2_A", "change in |z2| A", 6), ("dz2_B", "change in |z2| B", 6), ("dz3_A", "change in |z3| A", 6), ("dz3_B", "change in |z3| B", 6), ("dz4_A", "change in |z4| A", 6), ("dz4_B", "change in |z4| B", 6)) for method in METHODS: lines.append(f" {method}") for key, words, digits in names: e = state[key][method] lines.append(f" {words}: " + "; ".join(f"{arm} {e[arm]['mean']:+.{digits}f} ({e[arm]['sd']:.{digits}f})" for arm in ARMS) + " | " + "; ".join(f"{arm} - sham {e[f'{arm}_minus_sham']['mean']:+.{digits}f}" f" ({e[f'{arm}_minus_sham']['sd']:.{digits}f})" for arm in ("up", "down"))) lines.append(f"\nTime: {seconds:.0f} s on {len(found)} workers.") return result, lines def targets(): frozen() print(f"Target check: old starts {TARGET_STARTS[0] + 1} to {TARGET_STARTS[-1] + 1} on {len(TARGET_STARTS)} workers" " (immediate changes only)...", flush=True) clock = time.time() with Pool(len(TARGET_STARTS)) as pool: found = pool.map(target_sweep, TARGET_STARTS) found.sort(key=lambda f: f["start"]) result, lines = target_summary(found, time.time() - clock) write_json(TARGETS_JSON, result) partial = TARGETS_TXT.with_name(TARGETS_TXT.name + ".partial") partial.write_text("\n".join(lines) + "\n") os.replace(partial, TARGETS_TXT) print("\n".join(lines)) print(f"\nWrote {TARGETS_TXT.name} and {TARGETS_JSON.name}") # --------------------------------------------------------------------------------------------------------------------- # the dry run def made_up(moment, rng): """A data set in the shape of the real one, for trying the report: the saved old starts 1 to 18 (repeated) as the base runs, the dry run's moment records at every moment, three moments not carried out, and random numbers for the growth and the coherence paths. Nothing in it is a result.""" n, per, steps = len(NEW), MOMENTS, int(round(AFTER / EVERY)) old = [saved_run(s) for s in range(18)] pick = np.arange(n) % 18 d = dict(starts=np.array([s + 1 for s in NEW]), T=np.array([T for _, _, T in dn.moment_table()]).reshape(n, per), z=np.stack([old[i]["z"] for i in pick]), logs=np.stack([old[i]["logs"] for i in pick])) for key in RECORD_KEYS: value = np.asarray(moment[key]) d[key] = np.broadcast_to(value, (n, per) + value.shape).copy() carried = np.ones((n, per), bool) carried[2, 3] = carried[17, 0] = carried[30, 9] = False d["carried"], d["reason"] = carried, np.where(carried, "", "made up: not carried out") logs = rng.normal(0.014, 0.02, (n, per, 6, steps)) logs[~carried] = np.nan t = np.arange(int(AFTER) + 1) C0, C1 = d["C0"][..., None, None], d["C1"][..., None] level = C0 + (C1 - C0) * np.exp(-t / 150.0) + rng.normal(0, 0.002, (n, per, 6, len(t))) paths = np.repeat(level[..., None, :], 2, -2).astype(np.complex64) paths[~carried] = np.nan d["branch_logs"], d["branch_z"] = logs, paths return d def dryrun(): """Both kinds of task on old start 1 for a short stretch, through the same piece, part and assembly code as the run; the not-carried-out path on a made-up state; the report on made-up numbers. Writes only to /tmp.""" frozen() shutil.rmtree(DRY, ignore_errors=True) DRY.mkdir(parents=True) parts, runlog = DRY / "parts", DRY / "runlog.json" todo = [dict(name="base_01", kind="base", seed=0, keep=(CHECK_AT,), run=CHECK_RUN, parts=str(parts)), dict(name="moment_01_01", kind="moment", seed=0, moment=0, T=CHECK_AT, after=CHECK_RUN - CHECK_AT, parts=str(parts))] header = dict(dry=True, written=stamp()) print(f"Dry run on old start 1 ({DRY}; nothing in the project is written)") piece(todo, parts, runlog, header) (parts / "moment_01_01.partial.npz").write_bytes(b"left by a piece cut off") piece(todo, parts, runlog, header) stray_gone = not (parts / "moment_01_01.partial.npz").exists() (parts / "moment_01_09.npz").write_bytes(b"not in the run log") piece(todo, parts, runlog, header) print(f" a half-written part was removed before the next piece: {stray_gone}; a part not in the run log was" f" removed: {not (parts / 'moment_01_09.npz').exists()}") data = assemble(parts, [1], [(1, 1, CHECK_AT)], DRY / "sim_dry.npz", run=CHECK_RUN, after=CHECK_RUN - CHECK_AT) with np.load(parts / "moment_01_01.npz") as p: moment = {key: p[key] for key in p.files} with np.load(parts / "base_01.npz") as p: base_same = bool(np.array_equal(p["z"], data["z"]) and np.array_equal(p["logs"], data["logs"])) branch_same = bool(np.array_equal(moment["branch_logs"], data["branch_logs"][0, 0], equal_nan=True)) first_sample = np.abs(moment["branch_z"][:, :, 0].astype(np.complex128)).mean(-1) print(f" assembled into {data['branch_z'].shape} coherence paths and {data['branch_logs'].shape} gap growths;" f" base run and branches as in their parts: {base_same and branch_same}") print(f" the moment was carried out: {bool(moment['carried'])}; each branch's first coherence sample equals its" f" record of C just after the change within {np.abs(first_sample - moment['C1']).max():.1e}" f" (stored to single precision); task {float(moment['seconds']):.1f} s for {CHECK_RUN - CHECK_AT:.0f} time" f" units, about {float(moment['seconds']) * AFTER / (CHECK_RUN - CHECK_AT):.0f} s for a full moment task") with np.load(parts / "base_01.npz") as p: us, uf = p["us"][0].copy(), p["uf"][0].copy() rng = np.random.default_rng(7) tight = np.exp(1j * 0.1 * rng.standard_normal(1500)) us[1, 0, :1500] = tight * (us[1, 0, :1500] * np.conj(us[0, 0, :1500])) # the gap as it was us[0, 0, :1500] = tight rec, _ = planned(us, uf, (DRAWS, 0, 99)) print(f" a made-up state with group A's coherence at {rec['R0'][0]:.3f}: carried out {rec['carried']}; reason:" f" {rec['reason']}") print(f" a target that is not a number: {amplitude(us[0, 0], np.zeros(us.shape[-1]), 0j, np.nan)}" " (not a number, 2: a solver failure)") made = made_up(moment, np.random.default_rng(1)) result, lines = analyse(made, dry=True) (DRY / "cause_test.txt").write_text("\n".join(lines) + "\n") write_json(DRY / "cause_test.json", result) figure(result, DRY / "cause_test.png") print(f" the report on made-up numbers: {DRY / 'cause_test.txt'} ({len(lines)} lines), cause_test.json and" " cause_test.png (not a result; not shown)") # --------------------------------------------------------------------------------------------------------------------- # the lock, the run and the report def freeze(): """The checks, then the lock. Refused once any new data exist: the plan can only be fixed before the run.""" found = [path.name for path in (SAVED, RUNLOG, PARTS) if path.exists()] if found: sys.exit(f"New data already exist ({', '.join(found)}): the plan can no longer be locked.") frozen() if not PLAN_SHA or sha256(PLAN) != PLAN_SHA: sys.exit("The plan is not the version Colson approved. Stopping.") c = checks() show_checks(c) problems = [name for name, ok in (("page", c["page"]), ("files", all(c["files"].values())), ("base run", c["base"]["repeats"]), ("no-change branches", c["unchanged"]["holds"]), ("planned changes", c["planned"]["carried"]), ("moment table", c["table"]), ("new data", not c["new_data"])) if not ok] if problems: sys.exit(f"Checks failed: {', '.join(problems)}. Stopping.") import scipy import matplotlib lock = dict(locked=stamp(), plan=PLAN.name, plan_sha256=PLAN_SHA, settings=SETTINGS, sha256=fingerprints(), checks=c, environment=dict(python=platform.python_version(), numpy=np.__version__, scipy=scipy.__version__, matplotlib=matplotlib.__version__), new_data="none: no saved run, no run log, no parts") write_json(LOCK, lock) print(f"Locked {shown(lock['locked'])}: {LOCK.name}") def run(): """One piece of the run. Prints counts and times only.""" lock = unchanged() frozen() if SAVED.exists(): print(f" nothing left to run: the report has assembled {SAVED.name}") return piece(tasks(), PARTS, RUNLOG, dict(locked=lock["locked"], plan_sha256=lock["plan_sha256"])) def report(): """Once all 462 tasks are done: assemble, check, decide and write the outputs; then remove the parts.""" lock = unchanged() frozen() log = json.loads(RUNLOG.read_text()) todo = tasks() if SAVED.exists(): with np.load(SAVED) as saved: data = {key: saved[key] for key in saved.files} else: missing = [task["name"] for task in todo if task["name"] not in log["parts"]] if missing: sys.exit(f"{len(missing)} of {len(todo)} tasks are not done yet: run the pieces first.") for task in todo: if sha256(PARTS / f"{task['name']}.npz") != log["parts"][task["name"]]["sha256"]: sys.exit(f"{task['name']}.npz differs from the run log. Stopping.") data = assemble(PARTS, [s + 1 for s in NEW], dn.moment_table(), SAVED, run=RUN, after=AFTER) shutil.rmtree(PARTS) assert [int(s) for s in data["starts"]] == [s + 1 for s in NEW] assert all(float(data[key]) == float(value) for key, value in ( ("n", N), ("k", K), ("dt", DT), ("settle", nt.SETTLE), ("run", RUN), ("every", EVERY), ("gap", GAP), ("seed", nt.LYAP_SEED), ("draws", DRAWS), ("after", AFTER))) result, lines = analyse(data, lock, log) write_json(OUT_JSON, result) OUT_TXT.write_text("\n".join(lines) + "\n") figure(result, FIGURE) print("\n".join(lines)) print(f"\nWrote {OUT_TXT.name}, {OUT_JSON.name} and {FIGURE.relative_to(ROOT)}") if __name__ == "__main__": commands = dict(check=check, targets=targets, dryrun=dryrun, freeze=freeze, run=run, report=report) command = sys.argv[1] if len(sys.argv) > 1 else "" if command not in commands: sys.exit(__doc__) commands[command]()