1"""Quantitative figures (vector PDF) generated from the frozen analysis dataset. Every figure shows n and uncertainty where relevant."""2import os, sys, json3sys.path.insert(0, os.path.dirname(__file__))4from common import *5import matplotlib; matplotlib.use("Agg")6import matplotlib.pyplot as plt7FIG = os.path.join(REPO, "results", "figures"); os.makedirs(FIG, exist_ok=True)8COL = {"A": "#7f7f7f", "B": "#4c72b0", "C": "#dd8452", "D": "#55a868"}9plt.rcParams.update({"font.size": 9, "axes.spines.top": False, "axes.spines.right": False, "pdf.fonttype": 42})1011def fig_primary(split):12 at = pd.read_csv(os.path.join(EST, f"{split}_arm_table.csv")); ct = pd.read_csv(os.path.join(EST, f"{split}_contrasts.csv")); ct = ct[ct["endpoint"] == "Y"]13 df = load_runs(split)14 fig, ax = plt.subplots(1, 2, figsize=(7.4, 3.1), gridspec_kw={"width_ratios": [1, 1.4]})15 # left: task-weighted rates with cluster-bootstrap CI (per arm)16 for i, a in enumerate(ARMS):17 d = df[df["treatment"] == a].groupby("task_id")["Y"].mean()18 est, lo, hi = cluster_bootstrap(d - 0, 3000)19 ax[0].bar(i, est, color=COL[a], width=0.65); ax[0].errorbar(i, est, yerr=[[est - lo], [hi - est]], color="k", capsize=3, lw=1)20 ax[0].text(i, min(hi + 0.03, 1.02), f"{len(d)} tasks\n{len(df[df['treatment']==a])} runs", ha="center", va="bottom", fontsize=6.5)21 ax[0].set_xticks(range(4)); ax[0].set_xticklabels([f"Arm {a}" for a in ARMS]); ax[0].set_ylim(0, 1.2); ax[0].set_ylabel("Verified completion $\\hat p_a$ (task-weighted)")22 ax[0].set_title("(a) Primary endpoint by arm", loc="left")23 # right: paired contrasts24 cts = ct.reset_index(drop=True)25 for i, r in cts.iterrows():26 ax[1].errorbar(r["estimate"], i, xerr=[[r["estimate"] - r["ci_low"]], [r["ci_high"] - r["estimate"]]], fmt="o", color="k" if r["role"] == "primary" else "#555", capsize=3)27 ax[1].text(0.99, i + 0.32, f"p={r['perm_p']:.3f}; tasks +/=/−: {int(r['n_tasks_favoring_first'])}/{int(r['n_tasks_tied'])}/{int(r['n_tasks_favoring_second'])}", fontsize=6, va="center", ha="right", transform=ax[1].get_yaxis_transform(), color="#444")28 ax[1].axvline(0, color="#999", lw=0.8); ax[1].set_yticks(range(len(cts))); ax[1].set_yticklabels([f"{c} ({r})" for c, r in zip(cts["contrast"], cts["role"])]); ax[1].invert_yaxis()29 ax[1].set_xlabel("Paired task-level difference (95% task-cluster bootstrap CI)"); ax[1].set_title("(b) Prespecified contrasts", loc="left")30 fig.tight_layout(); fig.savefig(os.path.join(FIG, f"fig_primary_{split}.pdf")); plt.close(fig)3132def fig_quality(split):33 df = load_runs(split)34 comps = ["q_correctness", "q_integrity", "q_reproducibility", "q_reporting", "q_baseline"]; labels = ["correctness", "integrity", "reproducibility", "reporting", "baseline"]35 fig, ax = plt.subplots(1, 2, figsize=(7.2, 2.8))36 w = 0.237 for i, a in enumerate(ARMS):38 d = df[df["treatment"] == a]; means = [d.groupby("task_id")[c].mean().mean() for c in comps]39 ax[0].bar(np.arange(len(comps)) + (i - 1.5) * w, means, width=w, color=COL[a], label=f"Arm {a} (n={len(d)})")40 ax[0].set_xticks(range(len(comps))); ax[0].set_xticklabels(labels, rotation=20); ax[0].set_ylim(0, 1.05); ax[0].set_ylabel("Mean component score (task-weighted)"); ax[0].legend(fontsize=6, frameon=False); ax[0].set_title("(a) Quality components", loc="left")41 # violations by type and arm42 cr = pd.read_csv(os.path.join(DATA, "criteria.csv")); cr = cr[cr["kind"] == "violation"].merge(df[["run_id", "treatment"]], on="run_id")43 if len(cr):44 tab = cr.groupby(["passed_or_type", "treatment"])["run_id"].nunique().unstack(fill_value=0).reindex(columns=ARMS, fill_value=0)45 n_runs = df.groupby("treatment").size().reindex(ARMS)46 rate = tab / n_runs47 rate.plot.bar(ax=ax[1], color=[COL[a] for a in ARMS], width=0.8); ax[1].set_ylabel("Share of runs with ≥1 violation of type"); ax[1].set_xlabel(""); ax[1].legend(fontsize=6, frameon=False, title=None)48 ax[1].set_xticklabels(ax[1].get_xticklabels(), rotation=20)49 ax[1].set_title("(b) Critical integrity violations by type", loc="left")50 fig.tight_layout(); fig.savefig(os.path.join(FIG, f"fig_quality_{split}.pdf")); plt.close(fig)5152def fig_cost(split):53 df = load_runs(split)54 fig, ax = plt.subplots(1, 2, figsize=(7.2, 2.8))55 for a in ARMS:56 d = df[df["treatment"] == a]57 ax[0].scatter(d["spent_usd"], d["spent_wall_s"] / 60, s=14, color=COL[a], alpha=0.7, label=f"Arm {a}", marker="o" if True else "x")58 ax[1].scatter(d["spent_usd"].mean(), task_weighted_rate(d, a), s=60, color=COL[a], edgecolor="k", zorder=3)59 ax[1].annotate(f"{a}", (d["spent_usd"].mean(), task_weighted_rate(d, a)), xytext=(5, 4), textcoords="offset points", fontsize=8)60 ax[1].errorbar(d["spent_usd"].mean(), task_weighted_rate(d, a), xerr=d["spent_usd"].std() / np.sqrt(max(1, len(d))), color=COL[a], capsize=2, lw=0.8)61 ax[0].set_xlabel("Inference cost per run (USD, client estimate)"); ax[0].set_ylabel("Wall time per run (min)"); ax[0].legend(fontsize=6, frameon=False); ax[0].set_title("(a) Cost and time of every run", loc="left")62 ax[1].set_xlabel("Mean cost per run (USD) ± s.e."); ax[1].set_ylabel("Intervention-free verified completion"); ax[1].set_ylim(0, 1.05); ax[1].set_title("(b) Quality–cost frontier", loc="left")63 fig.tight_layout(); fig.savefig(os.path.join(FIG, f"fig_cost_{split}.pdf")); plt.close(fig)6465def fig_family(split):66 fam = pd.read_csv(os.path.join(EST, f"{split}_family_table.csv"))67 fams = sorted(fam["family"].unique()); fig, ax = plt.subplots(figsize=(7.2, 2.6)); w = 0.268 for i, a in enumerate(ARMS):69 d = fam[fam["treatment"] == a].set_index("family").reindex(fams)70 ax.bar(np.arange(len(fams)) + (i - 1.5) * w, d["Y"], width=w, color=COL[a], label=f"Arm {a}")71 for j, f in enumerate(fams):72 if i == 0: ax.text(j, 1.02, f"{int(d.loc[f, 'n_tasks'])} tasks", ha="center", fontsize=7)73 ax.set_xticks(range(len(fams))); ax.set_xticklabels(fams); ax.set_ylim(0, 1.12); ax.set_ylabel("Completion (task-weighted)"); ax.legend(fontsize=7, frameon=False, ncol=4); ax.set_title("Task-family heterogeneity of the primary endpoint", loc="left")74 fig.tight_layout(); fig.savefig(os.path.join(FIG, f"fig_family_{split}.pdf")); plt.close(fig)7576def fig_errors(split):77 df = load_runs(split); cr = pd.read_csv(os.path.join(DATA, "criteria.csv")).merge(df[["run_id", "treatment", "task_id"]], on="run_id")78 fails = cr[(cr["kind"] == "criterion") & (cr["passed_or_type"].astype(str) == "0")]79 if fails.empty: return80 fails["error"] = fails["task_id"] + ":" + fails["id"]81 top = fails.groupby("error")["run_id"].nunique().sort_values(ascending=False).head(25)82 tab = fails[fails["error"].isin(top.index)].groupby(["error", "treatment"])["run_id"].nunique().unstack(fill_value=0).reindex(columns=ARMS, fill_value=0).loc[top.index]83 fig, ax = plt.subplots(figsize=(7.2, 0.22 * len(tab) + 1.2)); tab.plot.barh(ax=ax, stacked=True, color=[COL[a] for a in ARMS], width=0.8)84 ax.invert_yaxis(); ax.set_xlabel("Number of runs failing the criterion (denominator: runs of that task per arm)"); ax.legend(fontsize=7, frameon=False); ax.set_title("Most frequent failed scientific criteria (task:criterion)", loc="left")85 fig.tight_layout(); fig.savefig(os.path.join(FIG, f"fig_errors_{split}.pdf")); plt.close(fig)8687if __name__ == "__main__":88 split = sys.argv[1] if len(sys.argv) > 1 else "confirmatory"89 for f in (fig_primary, fig_quality, fig_cost, fig_family, fig_errors):90 try: f(split); print("ok", f.__name__)91 except Exception as e: print("skip", f.__name__, e)92