"""Headroom-law forecast: the weaker-model fit scored on the Claude Opus 5.5 rows (frontier report, headroom-law edition). Local re-tabulation only: no model API, VM or evaluation. Part 1 (sections A-C and SENSITIVITY) follows the synthesis paper's re-test of the law (its code: same ledger, same selection rule, same fitting code and model definitions for the models reported here); only the output path and a stdout tee differ. Given that re-test's output JSON (its path in the environment variable RETEST_JSON), the script checks every part-1 value against it and prints every mismatch. Part 2 (section D) adds: D1 profile-likelihood 95% CIs for k on the weaker rows (6) and on all 21 rows, same method as the Opus-15 CI D2 out-of-sample forecast error of the weaker-fitted models on the 15 Opus 5.5 rows (MAE, RMSE, mean |e|/s, row-wise best-model counts, per-row predictions and errors) D3 the reverse direction (Opus-fitted models scoring the 6 weaker rows) D4 leave-one-out on the Opus 15: per-row predictions and LOO MAE / RMSE D5 the 15 Opus rows and 6 weaker rows with component, instrument, control score, h, Delta, s, r = Delta / h Input (read only): forensics/ledger.csv in the paper directory: the 195-row comparison record. In the released ledger.csv these free-text cells are empty, and every numeric column and every categorical column except team is as produced: note on L009 L036 L039 L042 L043 L052 L073 L079 L083 L086 L091 L099 L100 L101 L102 L104 L105 L107 L108 L109 L110 L112 L115 L121 L122 L125 L129 L130 L138 L139 L140 L148 L150 L165 L166 L168 L176 L178 L193; settings_parity on L036 L037 L039 L052 L086 L102 L108 L110 L139 L140 L142 L143 L148 L193; instrument on L169 L170 L171 L173; contrast on L169 L170 L171 L172 L173. Agent and team names are replaced by neutral IDs or plain study names in the team column and in the instrument, contrast, settings_parity and note cells of L030 L034 L038 L044 L045 L052 L082 L103 L107-L110 L133 L134 L151 L152 L184-L191 L193. Every source cell says in plain words where the comparison was recorded (internal document names, paths and commit ids removed), as do the note of L068 and the instrument of L146. No number this script prints reads those cells. Run (from the paper directory): python3 -I evidence/headroom_forecast.py Writes headroom_forecast.json and headroom_forecast.txt beside this file. Every model has a between-comparison spread tau and the comparison's own standard error s: Delta_i ~ N(mu_i, s_i^2 + tau^2), Delta and h in points of a 0-100 scale, h = 100 - control score. M0 zero mu = 0 M1 constant mu = c M2 headroom (law) mu = k h (equivalently: a constant fraction k of the control's failures fixed) M2b headroom + shift mu = k h - b (the record's model code: 'headbias') M3 capped constant mu = min(c, h) (a constant gain that cannot exceed the headroom) M5 smooth cap mu = c (1 - exp(-h / c)), c > 0 (slope 1 at small h, saturating at c) """ import io import json import math import os import sys import numpy as np import pandas as pd from scipy import optimize, stats HERE = os.path.dirname(os.path.abspath(__file__)) OUT = os.path.join(HERE, "headroom_forecast.json") TXT = os.path.join(HERE, "headroom_forecast.txt") PAPER = os.path.dirname(HERE) ORIG_JSON = os.environ.get("RETEST_JSON", "") # optional: the re-test's output JSON class _Tee: def __init__(self, *streams): self.streams = streams def write(self, s): for st in self.streams: st.write(s) def flush(self): for st in self.streams: st.flush() _BUF = io.StringIO() sys.stdout = _Tee(sys.__stdout__, _BUF) # ============================================================================================================ # PART 1: the synthesis paper's re-test code, from here to the end of SENSITIVITY # ============================================================================================================ LEDGER = os.path.join(PAPER, "forensics", "ledger.csv") RES = {} L = pd.read_csv(LEDGER, keep_default_na=False) for c in ["est", "lo", "hi", "p", "n", "se", "z", "ours", "ctrl", "mx", "headroom_pts"]: L[c] = pd.to_numeric(L[c], errors="coerce") FLAGS = L["flags"] LAYER = {"layer_vs_bare", "layer_vs_prior", "layer_vs_rival", "carried_w2s", "carried_s2w", "carried_self", "system_vs_single", "trained_vs_frontier"} WEAK = {"haiku", "qwen", "gpt-oss-20b", "opus-4.8"} # ================================================================== B. selection: the rule, criterion by criterion CRIT = [ ("a component compared with a control (ledger type is a layer comparison)", lambda d: d.ctype.isin(LAYER)), ("not a structural row (control cannot do the task at all) and not a verdict traced to a fixed harness bug", lambda d: ~d["flags"].str.contains("structural") & ~d["flags"].str.contains("harness_bug")), ("the control is a no-harness control (bare, told-bare, diligence instruction, fresh/alone), not a rival system", lambda d: ~d.ctype.isin(["layer_vs_rival", "system_vs_single"])), ("a bounded metric with a recorded control score (headroom computable on a 0-100 scale)", lambda d: d.ctrl.notna() & d.mx.notna() & d.headroom_pts.notna()), ("a standard error is recorded", lambda d: d.se.notna()), ("the harnessed model is Claude Opus 5.5", lambda d: d.model == "opus-5.5"), ] mask = pd.Series(True, index=L.index) seq = [] for name, f in CRIT: keep = f(L) removed = int((mask & ~keep).sum()) mask &= keep seq.append(dict(criterion=name, removed=removed, remaining=int(mask.sum()))) only = [] for j, (name, f) in enumerate(CRIT): others = pd.Series(True, index=L.index) for i, (_, g) in enumerate(CRIT): if i != j: others &= g(L) fails_only = L[others & ~f(L)] only.append(dict(criterion=name, rows_failing_only_this=fails_only.id.tolist())) RES["selection"] = dict(total=int(len(L)), sequential=seq, failing_only_one=only) print("=" * 100) print("B. SELECTION from the 195-row ledger") for s in seq: print(f" {s['removed']:3d} removed, {s['remaining']:3d} remain: {s['criterion']}") for o in only: print(f" fails only '{o['criterion'][:60]}...': {o['rows_failing_only_this']}") def prep(d): d = d.copy() d["y"] = d.est * 100 / d.mx d["s"] = d.se * 100 / d.mx d["h"] = d.headroom_pts return d A = prep(L[mask]) # the 15 Opus 5.5 rows keep_all_but_model = pd.Series(True, index=L.index) for name, f in CRIT[:-1]: keep_all_but_model &= f(L) W = prep(L[keep_all_but_model & L.model.isin(WEAK)]) # the weaker-model rows print("\nOpus 5.5 rows:", A.id.tolist()) print("weaker-model rows:", W[["id", "model"]].values.tolist()) # ================================================================== models def mean(kind, th, h): if kind == "M0": return np.zeros_like(h) if kind == "M1": return np.full_like(h, th[0]) if kind == "M2": return th[0] * h if kind == "M2b": return th[0] * h - th[1] if kind == "M3": return np.minimum(th[0], h) if kind == "M5": c = math.exp(th[0]) return c * (1 - np.exp(-h / c)) raise ValueError(kind) NPAR = {"M0": 0, "M1": 1, "M2": 1, "M2b": 2, "M3": 1, "M5": 1} STARTS = {"M0": [[]], "M1": [[0.0], [5.0]], "M2": [[0.3], [0.6]], "M2b": [[0.3, 0.0], [0.6, 1.0]], "M3": [[3.0], [8.0], [0.0]], "M5": [[math.log(3)], [math.log(10)], [math.log(30)]]} def nll(p, d, kind, tau_fixed=None): th, tau = (p, tau_fixed) if tau_fixed is not None else (p[:-1], math.exp(p[-1])) sd = np.sqrt(d.s.values ** 2 + tau ** 2) return -stats.norm.logpdf(d.y.values, mean(kind, th, d.h.values), sd).sum() def fit(d, kind, tau_fixed=None): best = None for st in STARTS[kind]: for lt in ([0.0, 1.5, 2.5] if tau_fixed is None else [None]): x0 = list(st) + ([lt] if lt is not None else []) if not x0: continue r = optimize.minimize(nll, x0, args=(d, kind, tau_fixed), method="Nelder-Mead", options=dict(maxiter=40000, xatol=1e-9, fatol=1e-11)) if best is None or r.fun < best.fun: best = r if best is None: # M0 with fixed tau: no free parameter return [], tau_fixed, nll([], d, kind, tau_fixed) th = list(best.x[:-1]) if tau_fixed is None else list(best.x) tau = math.exp(best.x[-1]) if tau_fixed is None else tau_fixed return th, tau, float(best.fun) def lpd(d, kind, th, tau): return stats.norm.logpdf(d.y.values, mean(kind, th, d.h.values), np.sqrt(d.s.values ** 2 + tau ** 2)) KINDS = ["M0", "M1", "M2", "M2b", "M3", "M5"] def table(d, label): n = len(d) out = {} for kind in KINDS: th, tau, f = fit(d, kind) k = NPAR[kind] + 1 loo = 0.0 for i in range(n): tr = d.drop(d.index[i]) th_i, tau_i, _ = fit(tr, kind) loo += float(lpd(d.iloc[[i]], kind, th_i, tau_i)[0]) out[kind] = dict(params=[round(float(x), 4) for x in th], tau=round(tau, 3), nll=round(f, 3), aic=round(2 * k + 2 * f, 2), bic=round(k * math.log(n) + 2 * f, 2), loo_lpd=round(loo, 3)) print(f"\n{label} (n={n})") print(f" {'model':5s} {'params':>22s} {'tau':>7s} {'AIC':>8s} {'BIC':>8s} {'LOO lpd':>9s}") for kind in KINDS: o = out[kind] shown = o["params"] if kind != "M5" else [round(math.exp(o["params"][0]), 3)] print(f" {kind:5s} {str(shown):>22s} {o['tau']:7.3f} {o['aic']:8.2f} {o['bic']:8.2f} {o['loo_lpd']:9.3f}") return out # ================================================================== A. reproduce the published fit def profile_k(d): th, tau, f = fit(d, "M2") grid = np.linspace(-0.5, 1.5, 2001) prof = np.array([optimize.minimize(lambda p, g=g: nll([g, p[0]], d, "M2"), [math.log(max(tau, 1e-3))], method="Nelder-Mead").fun for g in grid]) ok = grid[prof - f <= 1.92] return th[0], float(ok.min()), float(ok.max()) print("\n" + "=" * 100) print("A. REPRODUCTION") k, klo, khi = profile_k(A) repro = dict(k=round(k, 3), k_ci=[round(klo, 3), round(khi, 3)]) for kind, lab in (("M0", "zero"), ("M1", "constant"), ("M2", "headroom")): th, tau, f = fit(A, kind) repro[lab] = dict(aic=round(2 * (NPAR[kind] + 1) + 2 * f, 2), tau=round(tau, 2), params=[round(x, 3) for x in th]) thw, tauw, _ = fit(W, "M2") repro["weaker_k"] = round(thw[0], 3) repro["weaker_tau"] = round(tauw, 2) th21, tau21, _ = fit(pd.concat([A, W]), "M2") repro["all21_k"] = round(th21[0], 3) # the published Bayes factor: k from the weaker rows, both models scored with the weaker rows' tau lp_hr = float(lpd(A, "M2", thw, tauw).sum()) lp_zero = float(lpd(A, "M0", [], tauw).sum()) repro["published_bf"] = dict(lpd_head=round(lp_hr, 2), lpd_zero=round(lp_zero, 2), bf=round(math.exp(lp_hr - lp_zero), 1), note="both scored with tau fitted on the weaker rows under the headroom model") RES["reproduction"] = repro print(json.dumps(repro, indent=1)) # ================================================================== C. stronger nulls on the 15 rows, and on 21 print("\n" + "=" * 100) print("C. MODEL COMPARISON") RES["models_opus15"] = table(A, "Opus 5.5 rows") RES["models_weaker6"] = table(W, "weaker-model rows") RES["models_all21"] = table(pd.concat([A, W]), "all 21 rows") # cross-model: parameters (and tau) fitted on the weaker rows, then the Opus rows scored def cross(src, dst, label): res = {} for kind in KINDS: th, tau, _ = fit(src, kind) res[kind] = dict(params=[round(float(x), 4) for x in th], tau=round(tau, 3), lpd_own_tau=round(float(lpd(dst, kind, th, tau).sum()), 3), lpd_head_tau=round(float(lpd(dst, kind, th, tauw if src is W else tau).sum()), 3)) # mean structure transferred, tau re-estimated on the destination (semi out-of-sample) if kind != "M0": f_t = optimize.minimize(lambda p: -lpd(dst, kind, th, math.exp(p[0])).sum(), [1.0], method="Nelder-Mead") res[kind]["lpd_tau_refit"] = round(-float(f_t.fun), 3) else: f_t = optimize.minimize(lambda p: -lpd(dst, "M0", [], math.exp(p[0])).sum(), [1.0], method="Nelder-Mead") res[kind]["lpd_tau_refit"] = round(-float(f_t.fun), 3) print(f"\n{label}") print(f" {'model':5s} {'params (src)':>22s} {'tau src':>8s} {'lpd own tau':>12s} {'lpd tau_w':>10s} {'lpd tau refit':>14s}") for kind in KINDS: o = res[kind] print(f" {kind:5s} {str(o['params']):>22s} {o['tau']:8.3f} {o['lpd_own_tau']:12.3f} {o['lpd_head_tau']:10.3f} " f"{o['lpd_tau_refit']:14.3f}") bfs = {} for other in KINDS: if other == "M2": continue bfs[other] = dict( own_tau=round(math.exp(res["M2"]["lpd_own_tau"] - res[other]["lpd_own_tau"]), 2), shared_tau_w=round(math.exp(res["M2"]["lpd_head_tau"] - res[other]["lpd_head_tau"]), 2), tau_refit=round(math.exp(res["M2"]["lpd_tau_refit"] - res[other]["lpd_tau_refit"]), 2)) print(" Bayes factors of M2 (headroom) against each model, on the destination rows:") for other, v in bfs.items(): print(f" M2:{other} own tau {v['own_tau']:>10} shared tau_w {v['shared_tau_w']:>10} tau refit {v['tau_refit']:>10}") return dict(models=res, bf_m2_vs=bfs) RES["cross_weaker_to_opus"] = cross(W, A, "weaker rows -> Opus 5.5 rows") RES["cross_opus_to_weaker"] = cross(A, W, "Opus 5.5 rows -> weaker rows") # which rows discriminate: per-row log-likelihood differences at each model's full-data MLE (Opus 15) fits = {kind: fit(A, kind) for kind in KINDS} rowtab = A[["id", "h", "y", "s"]].copy() for kind in ["M1", "M2", "M3"]: th, tau, _ = fits[kind] rowtab[f"mu_{kind}"] = mean(kind, th, A.h.values) rowtab[f"ll_{kind}"] = lpd(A, kind, th, tau) for other in ["M1", "M3"]: rowtab[f"M2-{other}"] = rowtab["ll_M2"] - rowtab[f"ll_{other}"] rowtab = rowtab.sort_values("h") print("\nPer-row log-likelihood advantage of M2 at each model's MLE (Opus rows, sorted by headroom):") print(rowtab[["id", "h", "y", "s", "mu_M1", "mu_M2", "mu_M3", "M2-M1", "M2-M3"]].round(2).to_string(index=False)) RES["per_row"] = rowtab.round(4).to_dict(orient="records") # ================================================================== sensitivity: rows excluded only by the SE criterion print("\n" + "=" * 100) print("SENSITIVITY: Opus rows excluded only because no SE is recorded") se_only = [r for r in only if r["criterion"].startswith("a standard error")][0]["rows_failing_only_this"] Z = prep(L[L.id.isin(se_only)]) print(Z[["id", "model", "instrument", "y", "h", "flags"]].to_string(index=False)) sens = {} opus_se_only = Z[Z.model == "opus-5.5"] for s_imp in (0.5, 1.0, 2.0): Zi = opus_se_only.copy() Zi["s"] = s_imp D = pd.concat([A, Zi]) row = {} for kind in ["M1", "M2", "M3"]: th, tau, f = fit(D, kind) row[kind] = dict(params=[round(float(x), 3) for x in th], tau=round(tau, 3), aic=round(2 * (NPAR[kind] + 1) + 2 * f, 2)) sens[f"imputed_s_{s_imp}"] = dict(rows=Zi.id.tolist(), **row) print(f" imputed s = {s_imp}: k {row['M2']['params']} AIC M1 {row['M1']['aic']} M2 {row['M2']['aic']} " f"M3 {row['M3']['aic']}") RES["sensitivity_se_only_rows"] = sens # ============================================================================================================ # PART 2: additions for the headroom-law edition # ============================================================================================================ print("\n" + "=" * 100) print("D. HEADROOM-LAW FORECAST ADDITIONS") CMP = ["M0", "M1", "M2", "M3"] # the forecast comparison set the paper reports FKINDS = ["M0", "M1", "M2", "M2b", "M3"] # models scored out of sample in D2-D4 NAMES = {"M0": "zero", "M1": "constant", "M2": "headroom", "M2b": "headroom+shift", "M3": "capped constant", "M5": "smooth cap"} D21 = pd.concat([A, W]) def _clean(v): if isinstance(v, (float, np.floating)): return None if not np.isfinite(v) else float(v) if isinstance(v, (np.integer,)): return int(v) return v # ---------------------------------------------------------------- D1 profile-likelihood CIs for k def profile_k_wide(d, lo=-0.5, hi=1.5): """profile_k above, on a grid that is widened by 1.0 each side whenever the interval touches its edge.""" th, tau, f = fit(d, "M2") while True: grid = np.linspace(lo, hi, int(round((hi - lo) / 0.001)) + 1) prof = np.array([optimize.minimize(lambda p, g=g: nll([g, p[0]], d, "M2"), [math.log(max(tau, 1e-3))], method="Nelder-Mead").fun for g in grid]) inside = prof - f <= 1.92 ok = grid[inside] if ok.min() > grid[0] and ok.max() < grid[-1]: break lo, hi = lo - 1.0, hi + 1.0 idx = np.flatnonzero(inside) contiguous = bool(np.all(np.diff(idx) == 1)) return th[0], float(ok.min()), float(ok.max()), float(tau), [round(float(lo), 3), round(float(hi), 3)], contiguous print("\nD1. Profile-likelihood 95% CI for k (tau profiled out; cut at min nll + 1.92; grid step 0.001)") D1 = {} for lab, d in (("opus15", A), ("weaker6", W), ("all21", D21)): k_, lo_, hi_, tau_, gr, contig = profile_k_wide(d) D1[lab] = dict(n=int(len(d)), k=round(k_, 3), k_ci95=[round(lo_, 3), round(hi_, 3)], tau=round(tau_, 3), grid=gr, interval_contiguous=contig) print(f" {lab:8s} n={len(d):2d} k = {k_:.3f} 95% CI [{lo_:.3f}, {hi_:.3f}] tau {tau_:.3f} grid {gr}" f" contiguous={contig}") D1["opus15_matches_reproduction"] = D1["opus15"]["k_ci95"] == RES["reproduction"]["k_ci"] D1["method"] = ("profile likelihood of M2 (mu = k h) with log tau profiled out by Nelder-Mead at each grid value of " "k; the 95% set is {k : nll_profile(k) - nll_min <= 1.92}; identical to profile_k in part 1, grid " "widened only if the set touches the grid edge") print(f" Opus-15 interval identical to the part-1 reproduction: {D1['opus15_matches_reproduction']}") RES["D1_profile_ci"] = D1 # ---------------------------------------------------------------- D2/D3 out-of-sample forecast errors def forecast(src, dst, label): rows = dst[["id", "model", "h", "y", "s"]].copy().reset_index(drop=True) summ = {} for kind in FKINDS: th, tau, _ = fit(src, kind) mu = mean(kind, th, dst.h.values) e = dst.y.values - mu rows[f"mu_{kind}"] = mu rows[f"err_{kind}"] = e shown = [round(float(x), 4) for x in th] summ[kind] = dict(name=NAMES[kind], params_src=shown, tau_src=round(float(tau), 3), mae=round(float(np.mean(np.abs(e))), 3), rmse=round(float(np.sqrt(np.mean(e ** 2))), 3), mean_abs_err_over_s=round(float(np.mean(np.abs(e) / dst.s.values)), 3), mean_abs_err_over_sqrt_s2_tau2=round(float(np.mean(np.abs(e) / np.sqrt(dst.s.values ** 2 + tau ** 2))), 3), mean_err_bias=round(float(np.mean(e)), 3), max_abs_err=round(float(np.max(np.abs(e))), 3)) def winners(cands): E = np.abs(np.column_stack([rows[f"err_{k}"].values for k in cands])) best = E.min(axis=1) strict = {k: 0 for k in cands} tied = {k: 0 for k in cands} per = [] for i in range(len(rows)): w = [cands[j] for j in range(len(cands)) if E[i, j] <= best[i] + 1e-9] if len(w) == 1: strict[w[0]] += 1 else: for k in w: tied[k] += 1 per.append("/".join(w)) return dict(candidates=cands, strictly_best_rows=strict, tied_best_rows=tied, n_rows_with_tie=sum( 1 for p in per if "/" in p)), per win_main, per_main = winners(CMP) rows["best_of_M0_M1_M2_M3"] = per_main pair = {} for other in FKINDS: if other == "M2": continue a2, ao = np.abs(rows["err_M2"].values), np.abs(rows[f"err_{other}"].values) pair[other] = dict(m2_closer=int((a2 < ao - 1e-9).sum()), other_closer=int((ao < a2 - 1e-9).sum()), tie=int((np.abs(a2 - ao) <= 1e-9).sum())) print(f"\n{label}: parameters (and tau) fitted on the source rows, mean predicted on the destination rows") print(f" {'model':5s} {'name':18s} {'params(src)':>16s} {'MAE':>7s} {'RMSE':>7s} {'|e|/s':>7s} " f"{'|e|/sd':>7s} {'bias':>7s} {'max|e|':>7s} {'strict best':>11s}") for kind in FKINDS: o = summ[kind] sb = win_main["strictly_best_rows"].get(kind, "-") print(f" {kind:5s} {o['name']:18s} {str(o['params_src']):>16s} {o['mae']:7.3f} {o['rmse']:7.3f} " f"{o['mean_abs_err_over_s']:7.3f} {o['mean_abs_err_over_sqrt_s2_tau2']:7.3f} {o['mean_err_bias']:7.3f} " f"{o['max_abs_err']:7.3f} {str(sb):>11s}") print(f" row-wise best among {CMP}: strict {win_main['strictly_best_rows']}, tied {win_main['tied_best_rows']}") print(" M2 closer than each other model (rows): " + ", ".join(f"{o}: {v['m2_closer']}-{v['other_closer']} (ties {v['tie']})" for o, v in pair.items())) show = ["id", "model", "h", "y", "s"] + [f"mu_{k}" for k in CMP] + [f"err_{k}" for k in CMP] + \ ["best_of_M0_M1_M2_M3"] print(rows.sort_values("h")[show].round(2).to_string(index=False)) recs = [{k: _clean(v) for k, v in r.items()} for r in rows.round(4).to_dict(orient="records")] return dict(summary=summ, rowwise_best=win_main, m2_vs_other_rowwise=pair, per_row=recs) RES["D2_forecast_weaker_to_opus"] = forecast(W, A, "D2. Weaker rows (6) -> Opus 5.5 rows (15)") RES["D3_forecast_opus_to_weaker"] = forecast(A, W, "D3. Opus 5.5 rows (15) -> weaker rows (6)") # ---------------------------------------------------------------- D4 leave-one-out on the Opus 15 print("\nD4. Leave-one-out on the Opus 5.5 rows: each row predicted by each model fitted on the other 14") loo = A[["id", "h", "y", "s"]].copy().reset_index(drop=True) loo_sum = {} for kind in FKINDS: preds = [] for i in range(len(A)): th_i, tau_i, _ = fit(A.drop(A.index[i]), kind) preds.append(float(mean(kind, th_i, A.h.values[i:i + 1])[0])) preds = np.array(preds) e = A.y.values - preds th_f, tau_f, _ = fit(A, kind) e_in = A.y.values - mean(kind, th_f, A.h.values) loo[f"loo_mu_{kind}"] = preds loo[f"loo_err_{kind}"] = e loo_sum[kind] = dict(name=NAMES[kind], loo_mae=round(float(np.mean(np.abs(e))), 3), loo_rmse=round(float(np.sqrt(np.mean(e ** 2))), 3), loo_mean_abs_err_over_s=round(float(np.mean(np.abs(e) / A.s.values)), 3), insample_mae=round(float(np.mean(np.abs(e_in))), 3), insample_rmse=round(float(np.sqrt(np.mean(e_in ** 2))), 3), loo_lpd=RES["models_opus15"][kind]["loo_lpd"]) E = np.abs(np.column_stack([loo[f"loo_err_{k}"].values for k in CMP])) best = E.min(axis=1) loo_best = {k: 0 for k in CMP} per = [] for i in range(len(loo)): w = [CMP[j] for j in range(len(CMP)) if E[i, j] <= best[i] + 1e-9] if len(w) == 1: loo_best[w[0]] += 1 per.append("/".join(w)) loo["loo_best_of_M0_M1_M2_M3"] = per print(f" {'model':5s} {'name':18s} {'LOO MAE':>8s} {'LOO RMSE':>9s} {'LOO |e|/s':>10s} {'in-sample MAE':>14s} {'LOO lpd':>8s}") for kind in FKINDS: o = loo_sum[kind] print(f" {kind:5s} {o['name']:18s} {o['loo_mae']:8.3f} {o['loo_rmse']:9.3f} {o['loo_mean_abs_err_over_s']:10.3f} " f"{o['insample_mae']:14.3f} {o['loo_lpd']:8.3f}") print(f" LOO row-wise strictly best among {CMP}: {loo_best}") print(loo.sort_values("h")[["id", "h", "y", "s"] + [f"loo_mu_{k}" for k in ["M0", "M1", "M2", "M3"]] + ["loo_best_of_M0_M1_M2_M3"]].round(2).to_string(index=False)) RES["D4_loo_opus15"] = dict(summary=loo_sum, rowwise_strictly_best=loo_best, per_row=[{k: _clean(v) for k, v in r.items()} for r in loo.round(4).to_dict(orient="records")]) # ---------------------------------------------------------------- D5 the rows themselves def rowtable(d): out = [] for _, r in d.sort_values("h").iterrows(): r = r.to_dict() # plain dict: Series attributes such as .flags shadow column names scale = 100.0 / r["mx"] n = _clean(r["n"]) out.append({ "id": r["id"], "model": r["model"], "team": r["team"], "outcome": r["outcome"], "component": r["contrast"].split(" - ")[0].strip(), "contrast": r["contrast"], "ctype": r["ctype"], "instrument": r["instrument"], "control_score": _clean(round(r["ctrl"] * scale, 4)), "harnessed_score": _clean(round(r["ours"] * scale, 4)), "h": _clean(round(r["h"], 4)), "delta": _clean(round(r["y"], 4)), "s": _clean(round(r["s"], 4)), "r_delta_over_h": _clean(round(r["y"] / r["h"], 4)) if r["h"] > 0 else None, "ci_lo": _clean(round(r["lo"] * scale, 4)), "ci_hi": _clean(round(r["hi"] * scale, 4)), "p": _clean(r["p"]), "n": int(n) if n is not None else None, "se_method": r["se_method"], "verdict": r["verdict"], "prereg": r["prereg"], "flags": r["flags"], "source": r["source"], "note": r["note"], "metric_scale_max": _clean(r["mx"])}) return out RES["D5_rows"] = dict(opus15=rowtable(A), weaker6=rowtable(W), note="sorted by h; control/harnessed/delta/s/h on a 0-100 scale (x 100/mx); r = delta/h " "(None where h = 0); component = left side of the ledger contrast") print("\nD5. Rows (0-100 scale, sorted by h)") for lab in ("opus15", "weaker6"): print(f" {lab}:") for r in RES["D5_rows"][lab]: rr = " n/a" if r["r_delta_over_h"] is None else f"{r['r_delta_over_h']:5.2f}" print(f" {r['id']} {r['model']:9s} ctrl {r['control_score']:6.1f} h {r['h']:5.1f} delta {r['delta']:6.2f} " f"s {r['s']:5.2f} r {rr} n {r['n']} {r['component'][:28]:28s} | {r['instrument'][:50]}") # ---------------------------------------------------------------- headline block for the figure and tables cw = RES["cross_weaker_to_opus"]["bf_m2_vs"] f2 = RES["D2_forecast_weaker_to_opus"]["summary"] RES["headline"] = dict( opus15=dict(k=RES["reproduction"]["k"], k_ci95=RES["reproduction"]["k_ci"], aic=dict(headroom=RES["models_opus15"]["M2"]["aic"], constant=RES["models_opus15"]["M1"]["aic"], zero=RES["models_opus15"]["M0"]["aic"], capped=RES["models_opus15"]["M3"]["aic"], smooth_cap=RES["models_opus15"]["M5"]["aic"])), weaker6=dict(k=RES["reproduction"]["weaker_k"], k_ci95=D1["weaker6"]["k_ci95"], tau=RES["reproduction"]["weaker_tau"]), all21=dict(k=RES["reproduction"]["all21_k"], k_ci95=D1["all21"]["k_ci95"]), forecast_bf_own_tau=dict(vs_zero=cw["M0"]["own_tau"], vs_constant=cw["M1"]["own_tau"], vs_capped=cw["M3"]["own_tau"], vs_smooth_cap=cw["M5"]["own_tau"]), forecast_bf_shared_tau_published=RES["reproduction"]["published_bf"]["bf"], forecast_error_weaker_to_opus={NAMES[k]: dict(mae=f2[k]["mae"], rmse=f2[k]["rmse"], mean_abs_err_over_s=f2[k]["mean_abs_err_over_s"], strictly_best_rows=RES["D2_forecast_weaker_to_opus"]["rowwise_best"][ "strictly_best_rows"].get(k)) for k in CMP}, ) print("\nHEADLINE") print(json.dumps(RES["headline"], indent=1)) # ---------------------------------------------------------------- reproduction check against the re-test's JSON print("\n" + "=" * 100) print("CHECK: every part-1 value against the synthesis paper's re-test of the law") orig = {} if ORIG_JSON and os.path.isfile(ORIG_JSON): with open(ORIG_JSON) as fh: orig = json.load(fh) mine = json.loads(json.dumps({k: RES.get(k) for k in orig}, default=float)) MISM = [] NLEAF = [0] def _cmp(a, b, path): if isinstance(a, dict): if not isinstance(b, dict): MISM.append(dict(path=path, committed=str(a)[:80], recomputed=str(b)[:80])) return for k in b: # every value computed here; the re-test's other models are not if k not in a: MISM.append(dict(path=f"{path}.{k}", committed=a.get(k), recomputed=b.get(k))) else: _cmp(a[k], b[k], f"{path}.{k}") elif isinstance(a, list): if not isinstance(b, list) or len(a) != len(b): MISM.append(dict(path=path, committed=a, recomputed=b)) return for i, (x, y) in enumerate(zip(a, b)): _cmp(x, y, f"{path}[{i}]") else: NLEAF[0] += 1 if a != b: MISM.append(dict(path=path, committed=a, recomputed=b)) if orig: _cmp(orig, mine, "") RES["reproduction_check"] = dict(source="the synthesis paper's re-test of the law", leaves_compared=NLEAF[0], mismatches=MISM, exact=len(MISM) == 0) print(f" leaves compared: {NLEAF[0]}; mismatches: {len(MISM)}") else: print(" the re-test's output JSON is not given (RETEST_JSON); check skipped") for m in MISM: print(f" {m['path']}: committed {m['committed']} recomputed {m['recomputed']}") with open(OUT, "w") as fh: json.dump(RES, fh, indent=1, default=float) print("\nwrote", os.path.relpath(OUT, PAPER)) sys.stdout = sys.__stdout__ with open(TXT, "w") as fh: fh.write(_BUF.getvalue()) print("wrote", os.path.relpath(TXT, PAPER))