"""Extra descriptives for the manuscript's empirical sections (Agent H). Everything is recomputed from data/validated/*.csv|parquet,
analysis/tables/*.csv and the archived first attempts of amendment A-010; nothing is typed by hand.
Output: analysis/tables/table14_extra_descriptives.csv (long format: block, key, arm/group, metric, value) and a printed summary.
Usage: .venv/bin/python analysis/robustness/extra_descriptives.py"""
from __future__ import annotations
import json, os, sys
import numpy as np, pandas as pd
REPO = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", ".."))
VAL = os.path.join(REPO, "data", "validated"); TAB = os.path.join(REPO, "analysis", "tables"); RAW = os.path.join(REPO, "data", "raw", "runs")
sys.path.insert(0, os.path.join(REPO, "analysis", "econometrics"))
from estimate import load, cluster_bootstrap, permutation_p, weighted_rate, ARMS
ROWS = []

def add(block, key, group, metric, value):
    ROWS.append({"block": block, "key": key, "group": group, "metric": metric, "value": value})

def tokens(df):
    # same supervision proxy as estimate.py main(): 5 min per unanswered question (assumption) + logged interventions
    df["H_min"] = df["human_intervention_min"].astype(float) + 5.0 * df["n_questions"].fillna(0).astype(float)
    tok = df["tokens"].apply(lambda x: json.loads(x) if isinstance(x, str) else (x or {}))
    for k in ("input", "output", "cache_read", "cache_create"):
        df["tok_" + k] = tok.apply(lambda t: t.get(k, 0)).astype(float)
    df["tok_total"] = df[["tok_input", "tok_output", "tok_cache_read", "tok_cache_create"]].sum(axis=1)
    return df

def main():
    pilot = tokens(load("pilot", 70)); sens = tokens(load("sensitivity", 70)); rec = tokens(load("recovery", 70))
    # ---- A. process descriptives by arm (pilot)
    for arm, d in pilot.groupby("arm"):
        for m in ("tok_input", "tok_output", "tok_cache_read", "tok_cache_create", "tok_total", "num_turns", "tool_time_s", "first_tool_s", "retries", "n_artifacts"):
            add("A_process", "mean", arm, m, float(d[m].mean()))
        for qn, qv in (("p10", 0.10), ("p50", 0.50), ("p90", 0.90), ("max", 1.0), ("min", 0.0)):
            add("A_process", qn, arm, "wall_s", float(d["T"].quantile(qv))); add("A_process", qn, arm, "cost_usd", float(d["C"].quantile(qv)))
            add("A_process", qn, arm, "tool_calls", float(d["tool_calls"].quantile(qv))); add("A_process", qn, arm, "quality", float(d["quality"].quantile(qv)))
        add("A_process", "sd", arm, "wall_s", float(d["T"].std())); add("A_process", "n", arm, "runs", int(len(d)))
    # ---- B. node × arm (Assumption: no node effect) — runs and acceptance
    for (node, arm), d in pilot.groupby(["node", "arm"]):
        add("B_node", node, arm, "n", int(len(d))); add("B_node", node, arm, "acceptance", float(d["S"].mean())); add("B_node", node, arm, "quality", float(d["quality"].mean()))
    for node, d in pilot.groupby("node"):
        add("B_node", node, "all", "n", int(len(d))); add("B_node", node, "all", "acceptance", float(d["S"].mean()))
        dd = d[d.arm != "A0"]; add("B_node", node, "A1-A4", "acceptance", float(dd["S"].mean())); add("B_node", node, "A1-A4", "n", int(len(dd)))
    # ---- C. questions / supervision proxy
    for arm, d in pilot.groupby("arm"):
        nq = d["n_questions"].fillna(0).astype(float)
        add("C_questions", "share_zero_questions", arm, "share", float((nq == 0).mean())); add("C_questions", "share_le_3", arm, "share", float((nq <= 3).mean()))
        add("C_questions", "mean", arm, "n_questions", float(nq.mean())); add("C_questions", "median", arm, "n_questions", float(nq.median())); add("C_questions", "max", arm, "n_questions", float(nq.max()))
        add("C_questions", "mean", arm, "H_min", float(d["H_min"].mean()))
    # ---- D. ALF with the human-time allowance unbounded (budget-only frontier) and pooled over occupations
    budgets = [0.25, 0.5, 0.75, 1.0, 1.5, 2.0, 3.0, 4.0]
    for arm in ARMS:
        g = pilot[pilot.arm == arm]
        for B in budgets:
            ok = ((g["S"] == 1) & (g["C"] <= B)).groupby(g["task_id"]).mean()
            add("D_alf_Hinf", f"B={B}", arm, "ALF_pooled_equal", float(ok.mean()))
            for H in (0.0, 15.0, 30.0):
                okh = ((g["S"] == 1) & (g["C"] <= B) & (g["H_min"] <= H)).groupby(g["task_id"]).mean()
                add("D_alf_H", f"B={B};H={H:g}", arm, "ALF_pooled_equal", float(okh.mean()))
        for occ, go in g.groupby("occupation"):
            for B in (1.0, 2.0, 4.0):
                ok = ((go["S"] == 1) & (go["C"] <= B)).groupby(go["task_id"]).mean()
                add("D_alf_Hinf_occ", f"B={B}", f"{occ}|{arm}", "ALF", float(ok.mean()))
    # ---- E. A-008 sensitivity: the 8 re-judged runs (RES-01, arms A1–A4, repeats 1–2, judged after the sanitisation fix)
    rejudged = [f"RES-01_{a}_r{r}" for a in ("A1", "A2", "A3", "A4") for r in (1, 2)]
    jd = pd.read_parquet(os.path.join(VAL, "judgments_pilot.parquet")); first = jd.groupby("run_id")["judged_utc"].min()
    add("E_rejudged", "n_rejudged_runs_present", "pilot", "n", int(sum(r in first.index for r in rejudged)))
    for r in rejudged: add("E_rejudged", r, "pilot", "first_judged_utc", str(first.get(r)))
    rng = np.random.default_rng(42)
    for label, dsub in (("full", pilot), ("excl_8_rejudged", pilot[~pilot.run_id.isin(rejudged)]), ("excl_RES-01", pilot[pilot.task_id != "RES-01"])):
        for arm in ARMS: add("E_rejudged", label, arm, "acceptance_taskw", weighted_rate(dsub, "S", arm))
        for x, y in (("A3", "A1"), ("A1", "A0"), ("A4", "A3"), ("A2", "A1")):
            est, lo, hi, J = cluster_bootstrap(dsub, "S", x, y, 5000, rng); p = permutation_p(dsub, "S", x, y, 20000, rng)
            add("E_rejudged", label, f"{x}-{y}", "estimate", est); add("E_rejudged", label, f"{x}-{y}", "ci_lo", lo); add("E_rejudged", label, f"{x}-{y}", "ci_hi", hi); add("E_rejudged", label, f"{x}-{y}", "perm_p", p); add("E_rejudged", label, f"{x}-{y}", "J", J)
    res = pilot[pilot.task_id == "RES-01"]
    for _, r in res.iterrows(): add("E_rejudged", r["run_id"], r["arm"], "quality", float(r["quality"])); add("E_rejudged", r["run_id"], r["arm"], "S", float(r["S"]))
    # ---- F. A-010: the three lost collections (first attempts archived as *.attempt1) and their reruns
    for rid in ("RE-01_A2_r1", "RE-01_A2_r2", "LEGAL-01_A2_r1"):
        p1 = os.path.join(RAW, rid, "raw_stream.jsonl.attempt1")
        ev = [json.loads(l) for l in open(p1) if l.strip()]
        ntool = sum(1 for e in ev if e.get("type") == "assistant" for c in (e.get("message", {}).get("content") or []) if isinstance(c, dict) and c.get("type") == "tool_use")
        resv = [e for e in ev if e.get("type") == "result"]
        add("F_a010", rid, "attempt1", "events", len(ev)); add("F_a010", rid, "attempt1", "tool_use", ntool)
        add("F_a010", rid, "attempt1", "result_subtype", resv[0].get("subtype") if resv else None); add("F_a010", rid, "attempt1", "cost_usd", float(resv[0].get("total_cost_usd")) if resv else None)
        add("F_a010", rid, "attempt1", "num_turns", resv[0].get("num_turns") if resv else None); add("F_a010", rid, "attempt1", "duration_s", float(resv[0].get("duration_ms")) / 1000 if resv else None)
        rc1 = json.load(open(os.path.join(RAW, rid, "receipt.json.attempt1"))); add("F_a010", rid, "attempt1", "det_completion", bool(rc1.get("det_completion"))); add("F_a010", rid, "attempt1", "det_score", float(rc1.get("det_score")))
        j1 = json.load(open(os.path.join(RAW, rid, "judgments.json.attempt1")))
        add("F_a010", rid, "attempt1", "quality_mean", float(np.mean([x.get("score") for x in j1])) if isinstance(j1, list) and j1 else None)
        r2 = pilot[pilot.run_id == rid].iloc[0]
        for m in ("D", "quality", "S", "C", "tool_calls", "T", "node"): add("F_a010", rid, "rerun", m, r2[m] if m == "node" else float(r2[m]))
    # ---- G. recovery stage: cost rule of estimate.py (A4 resumed → client's cumulative total; A3 killed → wall-time slope fitted on pilot)
    ref = pd.read_parquet(os.path.join(VAL, "runs_pilot.parquet")); ref = ref[ref["cost_usd"].fillna(0) > 0]
    w = ref["wall_s"].astype(float).values; y = ref["cost_usd"].astype(float).values; beta = float((w * y).sum() / (w * w).sum()); pred = beta * w
    add("G_recovery_cost", "wall_slope_fit", "claude-opus-5-5", "beta_usd_per_s", beta); add("G_recovery_cost", "wall_slope_fit", "claude-opus-5-5", "n_runs", int(len(ref)))
    add("G_recovery_cost", "wall_slope_fit", "claude-opus-5-5", "R2", float(1 - ((y - pred) ** 2).sum() / ((y - y.mean()) ** 2).sum())); add("G_recovery_cost", "wall_slope_fit", "claude-opus-5-5", "mean_abs_err_usd", float(np.abs(y - pred).mean()))
    for arm, d in rec.groupby("arm"):
        add("G_recovery_cost", "cost_rule", arm, "n_imputed_from_wall", int(d["cost_imputed_flag"].sum())); add("G_recovery_cost", "cost_rule", arm, "n_cumulative_resume", int(d.get("cost_cumulative_resume", pd.Series(False, index=d.index)).sum()))
        add("G_recovery_cost", "total", arm, "mean_usd", float(d["C"].mean())); add("G_recovery_cost", "total", arm, "min_usd", float(d["C"].min())); add("G_recovery_cost", "total", arm, "max_usd", float(d["C"].max()))
        add("G_recovery_cost", "wall_total_s", arm, "mean", float(d["wall_total_s"].mean())); add("G_recovery_cost", "wall_total_s", arm, "median", float(d["wall_total_s"].median()))
        add("G_recovery_cost", "wall_total_s", arm, "max", float(d["wall_total_s"].max())); add("G_recovery_cost", "wall_first_s", arm, "mean", float(d["wall_s"].mean())); add("G_recovery_cost", "wall_first_s", arm, "max", float(d["wall_s"].max()))
        add("G_recovery_cost", "n_artifacts", arm, "mean", float(d["n_artifacts"].mean())); add("G_recovery_cost", "n_artifacts", arm, "share_zero", float((d["n_artifacts"] == 0).mean()))
        add("G_recovery_cost", "tool_calls", arm, "mean", float(d["tool_calls"].mean())); add("G_recovery_cost", "n_questions", arm, "mean", float(d["n_questions"].fillna(0).mean()))
        add("G_recovery_cost", "quality_missing", arm, "n", int(d["quality"].isna().sum())); add("G_recovery_cost", "S", arm, "sum", float(d["S"].sum())); add("G_recovery_cost", "S", arm, "n_nonmissing", int(d["S"].notna().sum()))
        add("G_recovery_cost", "critical", arm, "n", int(d["E"].sum())); add("G_recovery_cost", "det", arm, "n", int(d["D"].sum()))
    rr = pd.read_csv(os.path.join(VAL, "receipts_recovery.csv"))
    for _, x in rr.iterrows():
        v = json.loads(x["violations"]) if isinstance(x["violations"], str) else []
        for i in v: add("G_recovery_cost", "violation", x["run_id"], i["id"], i.get("description", "")[:120])
    rs = pd.read_csv(os.path.join(VAL, "receipts_sensitivity.csv"))
    for _, x in rs.iterrows():
        v = json.loads(x["violations"]) if isinstance(x["violations"], str) else []
        for i in v: add("J_sensitivity", "violation", x["run_id"], i["id"], i.get("description", "")[:120])
    # ---- H. spend by stage (validated data): worker cost (observed or imputed) and judge cost
    for tag, d in (("pilot", pilot), ("sensitivity", sens), ("recovery", rec)):
        j = pd.read_parquet(os.path.join(VAL, f"judgments_{tag}.parquet"))
        add("H_spend", tag, "worker", "sum_usd", float(d["C"].sum())); add("H_spend", tag, "worker", "n_runs", int(len(d))); add("H_spend", tag, "worker", "n_cost_imputed", int(d["cost_imputed_flag"].sum()))
        add("H_spend", tag, "judges", "sum_usd", float(j["cost_usd"].sum())); add("H_spend", tag, "judges", "n_judgments", int(len(j))); add("H_spend", tag, "judges", "mean_per_run_usd", float(j.groupby("run_id")["cost_usd"].sum().mean()))
        for jm, jj in j.groupby("judge_model"): add("H_spend", tag, f"judge:{jm}", "mean_cost_usd", float(jj["cost_usd"].mean()))
        add("H_spend", tag, "runs", "n_judges_0", int((d["n_judges"] == 0).sum())); add("H_spend", tag, "runs", "n_judges_3", int((d["n_judges"] == 3).sum()))
        add("H_spend", tag, "runs", "timed_out", int(d["censored"].sum())); add("H_spend", tag, "runs", "result_subtypes", ";".join(sorted(d["result_subtype"].astype(str).unique())))
        add("H_spend", tag, "runs", "nodes", ";".join(f"{k}:{v}" for k, v in d["node"].value_counts().sort_index().items()))
    # ---- I. judge-vote agreement (acceptable votes) vs D and vs Q>=70, and judge means by arm × model
    pilot["acc_votes"] = pilot["judge_acceptable_votes"].apply(lambda s: sum(1 for v in json.loads(s) if v) if isinstance(s, str) else np.nan)
    for (dv, av), n in pilot.groupby(["D", "acc_votes"]).size().items(): add("I_judges", f"D={int(dv)}", f"votes={int(av)}", "n", int(n))
    for (qv, av), n in pilot.groupby([(pilot["quality"] >= 70), "acc_votes"]).size().items(): add("I_judges", f"Q>=70={bool(qv)}", f"votes={int(av)}", "n", int(n))
    add("I_judges", "unanimous", "all", "share", float(pilot["acc_votes"].isin([0, 3]).mean()))
    jd2 = jd.merge(pilot[["run_id", "arm"]], on="run_id")
    for (arm, jm), g in jd2.groupby(["arm", "judge_model"]): add("I_judges", arm, jm, "mean_score", float(g["score"].mean()))
    for (arm, ji), g in jd2.groupby(["arm", "judge_idx"]): add("I_judges", arm, f"judge{ji}", "mean_score", float(g["score"].mean()))
    piv = jd.pivot_table(index="run_id", columns="judge_idx", values="score")
    add("I_judges", "opus_minus_sonnet_mean", "all", "points", float((piv[0] - piv[[1, 2]].mean(axis=1)).mean()))
    add("I_judges", "abs_diff_opus_sonnet", "all", "median_points", float((piv[0] - piv[[1, 2]].mean(axis=1)).abs().median()))
    add("I_judges", "abs_diff_opus_sonnet", "all", "p90_points", float((piv[0] - piv[[1, 2]].mean(axis=1)).abs().quantile(0.9)))
    add("I_judges", "abs_diff_sonnet_sonnet", "all", "median_points", float((piv[1] - piv[2]).abs().median()))
    add("I_judges", "runs_with_judge_disagreement_on_q", "all", "n", int(((piv >= 70).sum(axis=1).isin([1, 2])).sum()))
    # ---- J. model sensitivity: process by model (A3)
    a3 = pilot[pilot.arm == "A3"].copy(); a3["model_tag"] = "opus"; sens["model_tag"] = sens["model"].map({"claude-sonnet-5-5": "sonnet", "claude-fable-5-1": "fable"})
    for tag, d in pd.concat([a3, sens]).groupby("model_tag"):
        for m in ("num_turns", "tool_errors", "app_transitions", "tok_output", "tok_cache_read", "n_questions"): add("J_sensitivity", tag, "mean", m, float(d[m].fillna(0).mean()))
        for qn, qv in (("p10", 0.1), ("p50", 0.5), ("p90", 0.9), ("max", 1.0)): add("J_sensitivity", tag, qn, "wall_s", float(d["T"].quantile(qv))); add("J_sensitivity", tag, qn, "cost_usd", float(d["C"].quantile(qv)))
        add("J_sensitivity", tag, "abstention_rate", "share", float(d["abstention_status"].isin(["needs_clarification", "abstained"]).mean()))
        add("J_sensitivity", tag, "n", "runs", int(len(d)))
    # ---- K. q-gate anatomy: deterministic passes failing the gate, by scenario
    gate_fail = pilot[(pilot["D"] == 1) & (pilot["E"] == 0) & (pilot["S"] == 0)]
    for t, g in gate_fail.groupby("task_id"): add("K_gate", t, "det_pass_quality_below_q", "n", int(len(g))); add("K_gate", t, "det_pass_quality_below_q", "quality_min", float(g["quality"].min())); add("K_gate", t, "det_pass_quality_below_q", "quality_max", float(g["quality"].max()))
    det_fail_hi = pilot[(pilot["D"] == 0) & (pilot["quality"] >= 70)]
    for _, r in det_fail_hi.iterrows(): add("K_gate", r["run_id"], "det_fail_quality_above_q", "quality", float(r["quality"]))
    adm = pilot[pilot.task_id == "ADM-01"]
    for arm, g in adm.groupby("arm"): add("K_gate", "ADM-01", arm, "quality_mean", float(g["quality"].mean())); add("K_gate", "ADM-01", arm, "n_questions_mean", float(g["n_questions"].fillna(0).mean()))
    # judges' per-criterion means for ADM-01 (which rubric criteria pull the score below q)
    jadm = jd[jd.run_id.str.startswith("ADM-01") & ~jd.run_id.str.contains("_A0_")]
    crit = {}
    for _, x in jadm.iterrows():
        sc = json.loads(x["scores"]).get("scores", {}) if isinstance(x["scores"], str) else {}
        for k, v in sc.items(): crit.setdefault(k, []).append(v)
    for k, v in crit.items(): add("K_gate", "ADM-01_A1-A4", k, "mean_criterion_score_0to4", float(np.mean(v)))
    # ---- L. dependency: within-arm correlation of acceptance with depth and with tool calls (A1–A4 pooled and per arm)
    tasks = {t: json.load(open(os.path.join(REPO, "benchmark", "scenarios", t, "task.json"))) for t in pilot.task_id.unique()}
    pilot["depth"] = pilot.task_id.map({t: v.get("designed_dependency_depth") for t, v in tasks.items()})
    for arm in ("A1", "A2", "A3", "A4"):
        d = pilot[pilot.arm == arm]; tm = d.groupby("task_id").agg(S=("S", "mean"), depth=("depth", "first"), tc=("tool_calls", "mean"))
        add("L_dependency", arm, "scenario_level", "spearman_S_depth", float(tm[["S", "depth"]].corr(method="spearman").iloc[0, 1]))
        add("L_dependency", arm, "scenario_level", "spearman_S_toolcalls", float(tm[["S", "tc"]].corr(method="spearman").iloc[0, 1]))
        add("L_dependency", arm, "scenario_level", "spearman_toolcalls_depth", float(tm[["tc", "depth"]].corr(method="spearman").iloc[0, 1]))
        add("L_dependency", arm, "run_level", "mean_toolcalls_accepted", float(d[d.S == 1]["tool_calls"].mean())); add("L_dependency", arm, "run_level", "mean_toolcalls_rejected", float(d[d.S == 0]["tool_calls"].mean()) if (d.S == 0).any() else np.nan)
    shell = pilot[pilot.arm != "A0"]; tm = shell.groupby("task_id").agg(S=("S", "mean"), depth=("depth", "first"), tc=("tool_calls", "mean"))
    add("L_dependency", "A1-A4", "scenario_level", "spearman_S_depth", float(tm[["S", "depth"]].corr(method="spearman").iloc[0, 1])); add("L_dependency", "A1-A4", "scenario_level", "spearman_S_toolcalls", float(tm[["S", "tc"]].corr(method="spearman").iloc[0, 1]))
    for dcls, g in shell.groupby(shell["depth"].map(lambda x: "standard(21-24)" if x <= 24 else ("complex(31-34)" if x <= 34 else "workflow(44)"))):
        add("L_dependency", "A1-A4", dcls, "acceptance", float(g["S"].mean())); add("L_dependency", "A1-A4", dcls, "n", int(len(g))); add("L_dependency", "A1-A4", dcls, "tool_calls_mean", float(g["tool_calls"].mean()))
    # ---- M. tool families per run (pilot) and denied calls context
    fam = pd.read_csv(os.path.join(TAB, "audit_families_pilot.csv"))
    for arm, g in fam.groupby("arm"):
        tot = g["calls"].sum(); add("M_families", arm, "total", "calls_per_run", float(tot / 30))
        for _, r in g.iterrows(): add("M_families", arm, r["family"], "share", float(r["calls"] / tot))

    # ---- N. cost accounting of the recovery stage, re-derived from the raw streams (per-API-call usage deduplicated by message id)
    #       Prices are not taken from a list: they are fitted on the 150 pilot runs' result usage (input, output, cache read, cache creation).
    def dedup_usage(path):
        ev = [json.loads(l) for l in open(path) if l.strip()]; by = {}
        for e in ev:
            if e.get("type") != "assistant": continue
            m = e.get("message", {}); u = m.get("usage"); mid = m.get("id")
            if u and mid: by[mid] = u
        res = [e for e in ev if e.get("type") == "result"]
        return {"calls": len(by), "inp": sum(u.get("input_tokens", 0) for u in by.values()), "cr": sum(u.get("cache_read_input_tokens", 0) for u in by.values()),
                "cc": sum(u.get("cache_creation_input_tokens", 0) for u in by.values()), "res_usage": res[0].get("usage") if res else None, "res_cost": res[0].get("total_cost_usd") if res else None}
    pil = pd.read_parquet(os.path.join(VAL, "runs_pilot.parquet")); U = []; calls = []
    for rid in pil.run_id:
        d = dedup_usage(os.path.join(RAW, rid, "raw_stream.jsonl")); u = d["res_usage"]
        U.append([u["input_tokens"], u["output_tokens"], u["cache_read_input_tokens"], u["cache_creation_input_tokens"], d["res_cost"], d["calls"]])
    U = np.array(U, dtype=float); price, *_ = np.linalg.lstsq(U[:, :4], U[:, 4], rcond=None); predc = U[:, :4] @ price
    for k, b in zip(("input", "output", "cache_read", "cache_create"), price): add("N_recovery_cost_v2", "fitted_price_usd_per_Mtoken", "claude-opus-5-5", k, float(b * 1e6))
    add("N_recovery_cost_v2", "fit_on_result_usage", "claude-opus-5-5", "R2", float(1 - ((U[:, 4] - predc) ** 2).sum() / ((U[:, 4] - U[:, 4].mean()) ** 2).sum()))
    add("N_recovery_cost_v2", "fit_on_result_usage", "claude-opus-5-5", "max_abs_err_usd", float(np.abs(U[:, 4] - predc).max()))
    pil_arm = pil["arm"].values
    for arm in ("A3", "A4"):
        sel = pil_arm == arm; add("N_recovery_cost_v2", "pilot_output_tokens_per_call", arm, "mean", float((U[sel, 1] / U[sel, 5]).mean()))
        add("N_recovery_cost_v2", "pilot_cost_per_wall_second", arm, "usd_per_s", float(pil.loc[sel, "cost_usd"].sum() / pil.loc[sel, "wall_s"].sum()))
        add("N_recovery_cost_v2", "pilot_output_share_of_cost", arm, "share", float((U[sel, 1] * price[1]).sum() / U[sel, 4].sum()))
    opc = {arm: (U[pil_arm == arm, 1] / U[pil_arm == arm, 5]).mean() for arm in ("A3", "A4")}
    cps = {arm: pil.loc[pil_arm == arm, "cost_usd"].sum() / pil.loc[pil_arm == arm, "wall_s"].sum() for arm in ("A3", "A4")}
    recv = pd.read_csv(os.path.join(VAL, "runs_recovery.csv")); nrows = []
    for _, r in recv.iterrows():
        d1 = dedup_usage(os.path.join(RAW, r["run_id"], "raw_stream.jsonl"))
        lb = (d1["inp"] * price[0] + d1["cr"] * price[2] + d1["cc"] * price[3]); est = lb + d1["calls"] * opc[r["arm"]] * price[1]
        row = {"run_id": r["run_id"], "arm": r["arm"], "first_calls": d1["calls"], "first_cost_lower_bound": lb, "first_cost_est_tokens": est, "first_cost_est_wall": cps[r["arm"]] * r["wall_s"], "reported_recovery_cost": r["recovery_cost_usd"]}
        p2 = os.path.join(RAW, r["run_id"], "raw_stream_recovery.jsonl")
        if os.path.exists(p2):
            d2 = dedup_usage(p2); u2 = d2["res_usage"] or {}
            row["recovery_only_cost_from_usage"] = u2.get("input_tokens", 0) * price[0] + u2.get("output_tokens", 0) * price[1] + u2.get("cache_read_input_tokens", 0) * price[2] + u2.get("cache_creation_input_tokens", 0) * price[3]
            row["recovery_result_cost"] = d2["res_cost"]; row["recovery_calls"] = d2["calls"]
        nrows.append(row)
    N = pd.DataFrame(nrows)
    for arm, g in N.groupby("arm"):
        for m in ("first_calls", "first_cost_lower_bound", "first_cost_est_tokens", "first_cost_est_wall", "reported_recovery_cost", "recovery_only_cost_from_usage", "recovery_result_cost", "recovery_calls"):
            if m in g and g[m].notna().any(): add("N_recovery_cost_v2", arm, m, "mean", float(g[m].mean()))
        if arm == "A4":
            tot = g["first_cost_est_tokens"] + g["recovery_only_cost_from_usage"]; add("N_recovery_cost_v2", arm, "total_est_tokens(first+recovery_only)", "mean", float(tot.mean()))
            add("N_recovery_cost_v2", arm, "reported_minus_recovery_only", "mean", float((g["reported_recovery_cost"] - g["recovery_only_cost_from_usage"]).mean()))
            add("N_recovery_cost_v2", arm, "reported_minus_recovery_only", "min", float((g["reported_recovery_cost"] - g["recovery_only_cost_from_usage"]).min()))
            add("N_recovery_cost_v2", arm, "corr_reported_vs_total_est", "pearson", float(np.corrcoef(g["reported_recovery_cost"], tot)[0, 1]))
            add("N_recovery_cost_v2", arm, "corr_reported_vs_recovery_only", "pearson", float(np.corrcoef(g["reported_recovery_cost"], g["recovery_only_cost_from_usage"])[0, 1]))
            S4 = rec[rec.arm == "A4"].set_index("run_id").loc[g["run_id"], "S"]
            add("N_recovery_cost_v2", arm, "cost_per_accepted_if_reported_is_cumulative", "usd", float(g["reported_recovery_cost"].sum() / S4.sum()))
            add("N_recovery_cost_v2", arm, "cost_per_accepted_tokens_est", "usd", float(tot.sum() / S4.sum()))
        if arm == "A3":
            S3 = rec[rec.arm == "A3"].set_index("run_id").loc[g["run_id"], "S"]
            add("N_recovery_cost_v2", arm, "cost_per_accepted_tokens_est", "usd", float(g["first_cost_est_tokens"].sum() / S3.sum()) if S3.sum() else np.nan)
    for _, r in N.iterrows(): add("N_recovery_cost_v2", r["run_id"], r["arm"], "first_cost_est_tokens", float(r["first_cost_est_tokens"])); 
    for _, r in N[N.arm == "A4"].iterrows(): add("N_recovery_cost_v2", r["run_id"], r["arm"], "recovery_only_cost_from_usage", float(r["recovery_only_cost_from_usage"])); add("N_recovery_cost_v2", r["run_id"], r["arm"], "reported_recovery_cost", float(r["reported_recovery_cost"]))

    # ---- O. session token totals from the client's per-model usage (result event), A0 tool-error anatomy, criteria counts per scenario
    def mu(x):
        try: d = json.loads(x)
        except Exception: return {}  # noqa: BLE001
        out = {}
        for m, v in d.items():
            for k in ("inputTokens", "outputTokens", "cacheReadInputTokens", "cacheCreationInputTokens", "thinkingTokens", "costUSD"): out[k] = out.get(k, 0) + (v.get(k, 0) or 0)
        return out
    for tag, d in (("pilot", pilot), ("sensitivity", sens)):
        m = pd.DataFrame(list(d["model_usage"].apply(mu))); m["grp"] = (d["arm"] if tag == "pilot" else d["model"]).values
        for g, gg in m.groupby("grp"):
            for k in ("inputTokens", "outputTokens", "cacheReadInputTokens", "cacheCreationInputTokens", "thinkingTokens"): add("O_tokens", tag, g, k + "_mean", float(gg[k].mean()))
            add("O_tokens", tag, g, "n_with_usage", int(gg["outputTokens"].notna().sum()))
    evp = pd.read_parquet(os.path.join(VAL, "events_pilot.parquet"))
    ev0 = evp[evp.run_id.isin(pilot[pilot.arm == "A0"].run_id)]
    add("O_tokens", "events_columns", "pilot", "cols", ";".join(evp.columns))
    if "error_class" in evp.columns:
        err = ev0[ev0["error_class"].notna() & (ev0["error_class"].astype(str) != "")]
        for (tool, ec), n in err.groupby(["tool_name", "error_class"]).size().sort_values(ascending=False).head(8).items(): add("O_a0_errors", str(tool), str(ec), "n", int(n))
        add("O_a0_errors", "total", "A0", "n", int(len(err)))
        pref = err["summary"].astype(str).str.replace("tool_result:error:", "", regex=False).str[:40]
        for k, n in pref.value_counts().head(6).items(): add("O_a0_errors", "summary_prefix", k, "n", int(n))
        for arm in ("A1", "A2", "A3", "A4"):
            e1 = evp[evp.run_id.isin(pilot[pilot.arm == arm].run_id)]; e1 = e1[e1["error_class"].notna() & (e1["error_class"].astype(str) != "")]
            for (tool, ec), n in e1.groupby(["tool_name", "error_class"]).size().sort_values(ascending=False).head(4).items(): add("O_a0_errors", f"{arm}:{tool}", str(ec), "n", int(n))
    rcp = pd.read_parquet(os.path.join(VAL, "receipts_pilot.parquet"))
    for t, g in rcp.assign(task_id=rcp.run_id.str.split("_").str[0]).groupby("task_id"):
        c = json.loads(g.iloc[0]["criteria"]) if isinstance(g.iloc[0]["criteria"], str) else g.iloc[0]["criteria"]
        req = [x for x in c if x.get("required", True)] if isinstance(c, list) else c
        add("O_criteria", t, "required_criteria", "n", int(len(req)) if isinstance(req, list) else None); add("O_criteria", t, "evaluator_version", "v", str(g.iloc[0]["evaluator_version"]))
        add("O_criteria", t, "oracle_sha256_prefix", "v", str(g.iloc[0]["oracle_sha256"])[:12]); add("O_criteria", t, "n_distinct_oracle_hashes", "n", int(g["oracle_sha256"].nunique()))
    out = pd.DataFrame(ROWS); out.to_csv(os.path.join(TAB, "table14_extra_descriptives.csv"), index=False)
    pd.set_option("display.width", 220); pd.set_option("display.max_rows", 2000); pd.set_option("display.max_colwidth", 60)
    print(out.to_string())

if __name__ == "__main__":
    main()