"""Agent G — preregistered estimation from data/validated/runs_<tag>.parquet.
Outputs: analysis/tables/*.csv (+ .tex) and analysis/derived/estimates_<tag>.json. All numbers in the manuscript come from here.
Usage: .venv/bin/python analysis/econometrics/estimate.py --tag pilot [--q 70] [--B 5000] [--perm 20000]"""
from __future__ import annotations
import argparse, json, os, warnings
import numpy as np, pandas as pd
warnings.filterwarnings("ignore")
REPO = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", ".."))
TAB = os.path.join(REPO, "analysis", "tables"); DER = os.path.join(REPO, "analysis", "derived")
ARMS = ["A0", "A1", "A2", "A3", "A4"]
CONTRASTS = [("A3", "A1", "primary"), ("A1", "A0", "secondary"), ("A4", "A3", "secondary"), ("A2", "A1", "secondary")]

def load(tag, q):
    df = pd.read_parquet(os.path.join(REPO, "data", "validated", f"runs_{tag}.parquet"))
    df = df[df["status"].isin(["collected", "evaluated"])].copy()   # metrics are computed before the run is marked evaluated
    df["S"] = ((df["det_completion"].astype(bool)) & (~df["critical_error"].astype(bool)) & (df["quality"].astype(float) >= q)).astype(float)
    df.loc[df["quality"].isna(), "S"] = np.nan
    df["D"] = df["det_completion"].astype(float); df["E"] = df["critical_error"].astype(float)
    df["C"] = df["total_cost_usd"].astype(float); df["T"] = df["wall_total_s"].astype(float)
    df["censored"] = df["timed_out"].astype(bool)
    # Killed sessions emit no `result` event, so their marginal cost is unobserved (cost_usd = 0). Two cases:
    # (a) an A4 run whose recovery session resumed the same conversation: the client reports a CUMULATIVE total_cost_usd for the
    #     resumed session (WP29 gotcha), stored in recovery_cost_usd → that total already covers both sessions, no imputation;
    # (b) a killed session without resume (A3 interrupted control): impute from wall time with a per-model slope fitted by least
    #     squares on the completed pilot runs (cost ≈ β·wall_s; R² = 0.93 for Opus; token counts streamed per event are not usable
    #     because per-event output usage is partial — see data section). Flagged in cost_imputed_flag.
    df["cost_observed"] = df["cost_usd"].fillna(0).astype(float) + df["recovery_cost_usd"].fillna(0).astype(float); df["cost_imputed_flag"] = False; df["cost_cumulative_resume"] = False
    try:
        ref = pd.read_parquet(os.path.join(REPO, "data", "validated", "runs_pilot.parquet")); ref = ref[ref["cost_usd"].fillna(0) > 0]
    except Exception: ref = df[df["cost_usd"].fillna(0) > 0]  # noqa: BLE001
    for mdl, g in df.groupby("model"):
        fit = ref[ref["model"] == mdl]
        resumed = g[(g["cost_usd"].fillna(0) == 0) & (g["recovery_cost_usd"].fillna(0) > 0)]
        df.loc[resumed.index, "C"] = resumed["recovery_cost_usd"].astype(float); df.loc[resumed.index, "cost_cumulative_resume"] = True
        miss = g[(g["cost_usd"].fillna(0) == 0) & (g["recovery_cost_usd"].fillna(0) == 0) & g["censored"]]
        if len(fit) >= 5 and len(miss):
            w = fit["wall_s"].astype(float).values; beta = float((w * fit["cost_usd"].values).sum() / (w * w).sum())
            df.loc[miss.index, "C"] = beta * miss["wall_s"].astype(float); df.loc[miss.index, "cost_imputed_flag"] = True
    return df

def task_means(df, col, arm):
    return df[df.arm == arm].groupby("task_id")[col].mean()

def weighted_rate(df, col, arm, weights=None):
    tm = task_means(df, col, arm).dropna()
    if weights is None: return float(tm.mean())
    w = weights.reindex(tm.index).fillna(0); return float((tm * w).sum() / w.sum())

def cluster_bootstrap(df, col, a, b, B, rng, weights=None):
    ta = task_means(df, col, a); tb = task_means(df, col, b); tasks = sorted(set(ta.index) & set(tb.index))
    d = np.array([ta[t] - tb[t] for t in tasks]); w = np.array([weights[t] if weights is not None else 1.0 for t in tasks])
    est = float((d * w).sum() / w.sum()); boots = []
    for _ in range(B):
        idx = rng.integers(0, len(tasks), len(tasks)); boots.append((d[idx] * w[idx]).sum() / w[idx].sum())
    lo, hi = np.percentile(boots, [2.5, 97.5]); return est, float(lo), float(hi), len(tasks)

def permutation_p(df, col, a, b, n_perm, rng):
    ta = task_means(df, col, a); tb = task_means(df, col, b); tasks = sorted(set(ta.index) & set(tb.index))
    d = np.array([ta[t] - tb[t] for t in tasks]); obs = d.mean(); cnt = 0
    for _ in range(n_perm):
        s = rng.choice([-1, 1], len(d)); cnt += abs((d * s).mean()) >= abs(obs) - 1e-12
    return float(cnt / n_perm)

def holm(pvals):
    idx = np.argsort(pvals); m = len(pvals); adj = np.empty(m); run = 0
    for r, i in enumerate(idx): run = max(run, (m - r) * pvals[i]); adj[i] = min(1.0, run)
    return adj.tolist()

def mixed_logit(df, col):
    try:
        import statsmodels.formula.api as smf
        d = df.dropna(subset=[col]).copy(); d["y"] = d[col].astype(int)
        md = smf.mixedlm if False else None
        import statsmodels.api as sm
        X = pd.get_dummies(d["arm"], drop_first=False).astype(float); X = X[[a for a in ARMS if a in X.columns]]
        m = sm.GLM(d["y"], X, family=sm.families.Binomial()).fit(cov_type="cluster", cov_kwds={"groups": pd.factorize(d["task_id"])[0]})
        out = {a: {"logit": float(m.params[a]), "se_cluster": float(m.bse[a]), "p": float(1 / (1 + np.exp(-m.params[a])))} for a in X.columns}
        return {"model": "cluster-robust logit (arm dummies, no intercept), clustered by scenario", "params": out}
    except Exception as e:  # noqa: BLE001
        return {"error": repr(e)}

def survival(df):
    try:
        from lifelines import KaplanMeierFitter
        out = {}
        for a in ARMS:
            d = df[df.arm == a]
            if len(d) == 0: continue
            km = KaplanMeierFitter(); km.fit(d["T"], event_observed=~d["censored"])
            out[a] = {"median_s": float(km.median_survival_time_), "n": int(len(d)), "censored": int(d["censored"].sum()),
                      "p_done_600": float(1 - km.survival_function_at_times(600).iloc[0]), "p_done_1200": float(1 - km.survival_function_at_times(1200).iloc[0])}
        return out
    except Exception as e:  # noqa: BLE001
        return {"error": repr(e)}

def alf_curves(df, budgets, weights=None, h_allow=0.0, q_col="S"):
    """ALF_o(B,H,q) = Σ_j w_oj P(S=1, C≤B, H≤H). H = 5 min per unanswered question (assumption) + logged interventions."""
    rows = []
    for occ, g in df.groupby("occupation"):
        for a in ARMS:
            ga = g[g.arm == a]
            if len(ga) == 0: continue
            for B in budgets:
                ok = (ga[q_col] == 1) & (ga["C"] <= B) & (ga["H_min"] <= h_allow)
                tm = ok.groupby(ga["task_id"]).mean()
                w = (weights.reindex(tm.index).fillna(0) if weights is not None else pd.Series(1.0, index=tm.index))
                rows.append({"occupation": occ, "arm": a, "B": B, "H": h_allow, "ALF": float((tm * w).sum() / w.sum())})
    return pd.DataFrame(rows)

def agreement(jdg):
    """Krippendorff's alpha (interval) across judges on the per-run quality scores + pairwise Pearson."""
    piv = jdg.pivot_table(index="run_id", columns="judge_idx", values="score")
    try:
        vals = piv.values.astype(float); n_units = vals.shape[0]
        obs = []; allv = vals[~np.isnan(vals)]
        De = np.var(allv) * 2 if len(allv) > 1 else np.nan
        Do = 0.0; m = 0
        for row in vals:
            r = row[~np.isnan(row)]
            if len(r) > 1:
                Do += np.sum((r[:, None] - r[None, :]) ** 2) / (len(r) - 1); m += len(r)
        Do = Do / m if m else np.nan
        alpha = 1 - Do / De if De and not np.isnan(De) else np.nan
    except Exception: alpha = np.nan  # noqa: BLE001
    corr = piv.corr().round(3).to_dict() if piv.shape[1] > 1 else {}
    return {"krippendorff_alpha_interval": float(alpha) if alpha == alpha else None, "pairwise_pearson": corr, "n_runs": int(piv.shape[0])}

def main():
    ap = argparse.ArgumentParser(); ap.add_argument("--tag", default="pilot"); ap.add_argument("--q", type=float, default=70); ap.add_argument("--B", type=int, default=5000); ap.add_argument("--perm", type=int, default=20000); ap.add_argument("--seed", type=int, default=42)
    a = ap.parse_args(); rng = np.random.default_rng(a.seed); os.makedirs(TAB, exist_ok=True); os.makedirs(DER, exist_ok=True)
    df = load(a.tag, a.q)
    # human-time proxy: 5 min per question in COMPLETION.json (assumption) + logged interventions
    def nq(x):
        try: return len(json.loads(x)) if isinstance(x, str) else 0
        except Exception: return 0  # noqa: BLE001
    df["H_min"] = df["human_intervention_min"].astype(float) + 5.0 * df.get("n_questions", pd.Series(0, index=df.index)).fillna(0).astype(float) if "n_questions" in df else df["human_intervention_min"].astype(float)
    # task weights: equal (primary) and professional minutes (sensitivity)
    tasks = {t: json.load(open(os.path.join(REPO, "benchmark", "scenarios", t, "task.json"))) for t in df.task_id.unique()}
    w_min = pd.Series({t: float(v.get("professional_minutes_estimate", 60)) for t, v in tasks.items()})
    # --- Table 1: outcomes by arm
    rows = []
    for arm in ARMS:
        d = df[df.arm == arm]
        if len(d) == 0: continue
        S = d["S"]; acc = S.sum(); rows.append({"arm": arm, "n_runs": len(d), "n_tasks": d.task_id.nunique(), "acceptance_taskw": weighted_rate(df, "S", arm), "acceptance_minw": weighted_rate(df, "S", arm, w_min),
            "det_completion": weighted_rate(df, "D", arm), "critical_error_rate": float(d["E"].mean()), "quality_mean": float(d["quality"].mean()), "quality_sd": float(d["quality"].std()),
            "timeout_rate": float(d["censored"].mean()), "cost_mean": float(d["C"].mean()), "cost_sd": float(d["C"].std()), "cost_per_accepted": float(d["C"].sum() / acc) if acc else np.nan,
            "wall_median_s": float(d["T"].median()), "tool_calls_mean": float(d["tool_calls"].mean()), "tool_errors_mean": float(d["tool_errors"].mean()), "app_transitions_mean": float(d["app_transitions"].mean()),
            "abstention_rate": float((d["abstention_status"].isin(["needs_clarification", "abstained"])).mean())})
    t1 = pd.DataFrame(rows); t1.to_csv(os.path.join(TAB, f"table1_outcomes_by_arm_{a.tag}.csv"), index=False)
    # --- Table 2: contrasts
    crows = []; pv = []
    for x, y, kind in CONTRASTS:
        for col in ("S", "D", "E", "quality"):
            if col == "quality":
                est, lo, hi, J = cluster_bootstrap(df, "quality", x, y, a.B, rng)
                p = permutation_p(df, "quality", x, y, a.perm, rng)
            else:
                est, lo, hi, J = cluster_bootstrap(df, col, x, y, a.B, rng); p = permutation_p(df, col, x, y, a.perm, rng)
            est_w = cluster_bootstrap(df, col, x, y, 200, rng, w_min)[0]
            crows.append({"contrast": f"{x}-{y}", "kind": kind, "outcome": col, "estimate": est, "ci_lo": lo, "ci_hi": hi, "perm_p": p, "J": J, "estimate_minutes_weighted": est_w})
            if col == "S" and kind == "secondary": pv.append(p)
    t2 = pd.DataFrame(crows)
    if pv:
        adj = holm(pv); k = 0
        t2["holm_p"] = np.nan
        for i, r in t2.iterrows():
            if r["outcome"] == "S" and r["kind"] == "secondary": t2.loc[i, "holm_p"] = adj[k]; k += 1
    t2.to_csv(os.path.join(TAB, f"table2_contrasts_{a.tag}.csv"), index=False)
    # --- Table 3: by occupation × arm
    t3 = df.groupby(["occupation", "arm"]).agg(n=("S", "size"), acceptance=("S", "mean"), det=("D", "mean"), quality=("quality", "mean"), critical=("E", "mean"), cost=("C", "mean"), tool_calls=("tool_calls", "mean")).reset_index()
    t3.to_csv(os.path.join(TAB, f"table3_occupation_arm_{a.tag}.csv"), index=False)
    # --- dependency: acceptance vs designed depth and vs tool calls (within arm)
    df["depth"] = df.task_id.map({t: v.get("designed_dependency_depth") for t, v in tasks.items()})
    dep = df.groupby(["arm", "depth"]).agg(n=("S", "size"), acceptance=("S", "mean"), tool_calls=("tool_calls", "mean")).reset_index(); dep.to_csv(os.path.join(TAB, f"table4_dependency_{a.tag}.csv"), index=False)
    # --- recovery index (A4)
    a4 = df[df.arm == "A4"]; fail_first = a4[a4["selfcheck_pass_first"] == 0]
    ari = {"n_A4": int(len(a4)), "n_failed_first_selfcheck": int(len(fail_first)), "n_recovery_attempted": int(a4["recovery_attempted"].astype(bool).sum()),
           "n_recovered": int((fail_first["selfcheck_pass_final"] == 1).sum()), "ARI": float((fail_first["selfcheck_pass_final"] == 1).mean()) if len(fail_first) else None,
           "acceptance_given_recovery": float(fail_first["S"].mean()) if len(fail_first) else None, "mean_recovery_cost": float(a4["recovery_cost_usd"].fillna(0).mean())}
    # --- ALF curves
    budgets = [0.25, 0.5, 0.75, 1, 1.5, 2, 3, 4]
    alf0 = alf_curves(df, budgets, None, 0.0); alf0["weights"] = "equal"; alf1 = alf_curves(df, budgets, w_min, 0.0); alf1["weights"] = "minutes"
    alfH = alf_curves(df, budgets, None, 15.0); alfH["weights"] = "equal_H15"
    alf = pd.concat([alf0, alf1, alfH]); alf.to_csv(os.path.join(TAB, f"table5_alf_{a.tag}.csv"), index=False)
    # --- sensitivity to q
    qs = []
    for q in (50, 60, 70, 80, 90):
        d2 = df.copy(); d2["S"] = ((d2["D"] == 1) & (d2["E"] == 0) & (d2["quality"] >= q)).astype(float)
        for arm in ARMS:
            if (d2.arm == arm).any(): qs.append({"q": q, "arm": arm, "acceptance": weighted_rate(d2, "S", arm)})
    pd.DataFrame(qs).to_csv(os.path.join(TAB, f"table6_q_sensitivity_{a.tag}.csv"), index=False)
    # --- judges agreement
    jdg = pd.read_parquet(os.path.join(REPO, "data", "validated", f"judgments_{a.tag}.parquet"))
    agr = agreement(jdg)
    jm = jdg.groupby("judge_model")["score"].agg(["mean", "std", "count"]).reset_index(); jm.to_csv(os.path.join(TAB, f"table7_judges_{a.tag}.csv"), index=False)
    # det vs judge-acceptable agreement
    df["judge_majority_acceptable"] = df["judge_acceptable_votes"].apply(lambda s: (sum(1 for v in json.loads(s) if v) >= 2) if isinstance(s, str) else None)
    cross = pd.crosstab(df["D"], df["judge_majority_acceptable"]); cross.to_csv(os.path.join(TAB, f"table8_det_vs_judges_{a.tag}.csv"))
    out = {"tag": a.tag, "q": a.q, "n_runs": int(len(df)), "n_tasks": int(df.task_id.nunique()), "total_cost_usd": float(df["C"].sum()), "by_arm": t1.to_dict("records"), "contrasts": t2.to_dict("records"),
           "mixed_logit_S": mixed_logit(df, "S"), "survival": survival(df), "recovery": ari, "agreement": agr, "models": sorted(df["model"].dropna().unique().tolist()),
           "claude_code_versions": sorted(df["claude_code_version"].dropna().unique().tolist())}
    json.dump(out, open(os.path.join(DER, f"estimates_{a.tag}.json"), "w"), indent=1, default=str)
    print(t1.round(3).to_string()); print(t2.round(3).to_string()); print(json.dumps(ari, indent=1)); print(json.dumps(agr, indent=1, default=str))

if __name__ == "__main__":
    main()