"""T3-synthesis: extraction, harmonisation onto theta = ln(T_AI/T_noAI), FE/DL pooling, forest plot.
Run from workspace root:  ~/wp29/venv/bin/python out/analysis.py
Reads only data/. Writes out/results.json, out/table1.csv, out/figure1.png. No randomness."""
import json, math, os, re, csv
import numpy as np
from scipy.stats import norm
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
D = os.path.join(ROOT, "data"); OUT = os.path.join(ROOT, "out")
Z = 1.96

def lines_of(src):
    with open(os.path.join(D, "sources", src + ".txt"), encoding="utf-8") as fh:
        return fh.read().split("\n")

def find(src, pattern, start=0):
    """Return (1-based line number, line) of the first line at/after `start` matching regex."""
    L = lines_of(src)
    for i in range(start, len(L)):
        if re.search(pattern, L[i]):
            print(f"[{src} L{i+1}] {L[i].strip()[:220]}")
            return i + 1, L[i]
    raise RuntimeError(f"pattern not found in {src}: {pattern}")

studies = json.load(open(os.path.join(D, "studies.json"), encoding="utf-8"))
ids = [s["id"] for s in studies]
assert ids == ["S1", "S2", "S3", "S4", "S5"]

# n_input_rows: wc -l semantics (count of newline characters)
n_input_rows = {"studies.json": len(studies)}
for fn in ["README.md"] + ["sources/" + f for f in sorted(os.listdir(os.path.join(D, "sources")))]:
    with open(os.path.join(D, fn), "rb") as fh:
        n_input_rows[fn] = fh.read().count(b"\n")

print("=== Source lines relied on ===")
items = {}

# ---------------- S1 Paradis et al. ----------------
src = "arxiv-2410.12944v3"
ln_tab, row = find(src, r"ExpCon \|")
ln_se, _ = find(src, r"^\(", ln_tab)  # next SE line after the ExpCon row
se_line = lines_of(src)[ln_se - 1]
coefs = re.findall(r"-?\d+\.\d+", row.split("ExpCon")[0])
ses = re.findall(r"\((\d+\.\d+)\)", se_line)
b_print, se_print = coefs[1], ses[1]            # Model 2 = second column
ln_ci, ci_line = find(src, r"best-fit Model 2")
m = re.search(r"CI = \[(-?\d+\.\d+),\s*(-?\d+\.\d+)\]\s*,\s*p=(\d?\.\d+)", ci_line)
ci_l, ci_u, p1 = m.group(1), m.group(2), m.group(3)
find(src, r"TABLE III")
ln_n, nline = find(src, r"completed the task in 96 min")
b1, se1 = float(b_print), float(se_print)
se1_ci = (float(ci_u) - float(ci_l)) / (2 * Z)
items["S1"] = dict(
    estimate_as_printed=f"{b_print} (ExpCon, Model 2, LogToT; marker '.'); text: beta_1=-0.24, effect of AI = 21%",
    uncertainty_as_printed=f"SE ({se_print}); 95% CI = [{ci_l}, {ci_u}]; p={p1}",
    sample_as_printed="Table III: n = 96 (Section VI-B: N=47 AI, N=46 no-AI)",
    definition="OLS coefficient of the experimental condition (AI vs no-AI, randomised) on the natural log of time on task, controlling for developer-level factors (Model 2)",
    locator=f"{src}; Table III (L{ln_tab}-L{ln_se}), row ExpCon, column (2) 'Dev'; Section VI-C (L{ln_ci}) for CI and p",
    harmonisable=True, harmonised_estimate=b1, harmonised_se=se1, poolable=True,
    transformation=f"log-time coefficient: theta = b = {b_print}; SE(theta) = SE(b) = {se_print} (printed SE used; CI-based SE = ({ci_u} - ({ci_l}))/3.92 = {se1_ci:.4f} used only as sensitivity)",
    reason_if_not=None)

# ---------------- S2 Becker et al. ----------------
src = "arxiv-2507.09089v2"
ln_a, _ = find(src, r"increases completion time by 19%")
ln_r, rline = find(src, r"implement issues in 19% more time on average")
ln_s, _ = find(src, r"S=\\frac")
ln_c, _ = find(src, r"Developers complete 136 issues")
find(src, r"^Omitted:.*Figure 1")
c2 = float(re.search(r"(\d+)% more time", rline).group(1)) / 100
th2 = math.log(1 + c2)
items["S2"] = dict(
    estimate_as_printed="19% more time on average (AI allowed increases completion time by 19%)",
    uncertainty_as_printed="none printed (confidence interval shown only in Figure 1 / Figure 15, not reproduced in the extract)",
    sample_as_printed="16 developers, 246 tasks; 136 issues AI-allowed, 110 issues AI-disallowed",
    definition="S = E[time with AI allowed]/E[time AI disallowed] - 1, estimated as exp(beta)-1 from OLS of log implementation time on AI indicator and log forecast time (issue-level randomisation)",
    locator=f"{src}; Section 3 Results (L{ln_r}); abstract (L{ln_a}); estimator in Section 2.3 (L{ln_s}) and Appendix D.1-D.3",
    harmonisable=True, harmonised_estimate=th2, harmonised_se=None, poolable=False,
    transformation=f"percentage change: c = +0.19; theta = ln(1 + 0.19) = {th2:.6f} (equals the regression beta since S = exp(beta) - 1)",
    reason_if_not="Not poolable: no SE, CI or p-value is printed in the extract; the confidence interval appears only in a figure that is not reproduced, so SE(theta) cannot be obtained.")

# ---------------- S3 Bono & Xu ----------------
src = "arxiv-2411.01067v2"
ln_h, _ = find(src, r"\| Overall \| Treatment \| Control \| Improvement")
ln_row, row3 = find(src, r"^All Content \|", ln_h)
cells = [c.strip() for c in row3.split("|")]
ov, tr, co, imp = cells[1], cells[2], cells[3], cells[4]
find(src, r"Table 3: Comparison of task duration")
find(src, r"We had 181 subjects")
find(src, r"^Editorial note")
c3 = float(re.match(r"(-?\d+\.\d+)%", imp).group(1)) / 100
th3 = math.log(1 + c3)
items["S3"] = dict(
    estimate_as_printed=f"Treatment {tr} min, Control {co} min (Overall {ov}); Improvement {imp}",
    uncertainty_as_printed="significance stars only (***; per the paper's scheme 0<p<=0.01); no SE or CI printed",
    sample_as_printed="181 subjects completed the task satisfactorily (one subject dropped)",
    definition="Percentage difference in mean task duration (minutes), Copilot (treatment) vs control without Copilot, holding accuracy constant (control duration predicted at treatment accuracy; Appendix B), randomised",
    locator=f"{src}; Section 4.2, Table 3 (L{ln_h}-L{ln_row}), row 'Aggregate Results / All Content', columns Treatment, Control, Improvement",
    harmonisable=True, harmonised_estimate=th3, harmonised_se=None, poolable=False,
    transformation=f"percentage change: c = {c3}; theta = ln(1 + ({c3})) = {th3:.6f} (ratio of displayed means ln({tr}/{co}) = {math.log(float(tr)/float(co)):.6f}, sensitivity only)",
    reason_if_not="Not poolable: uncertainty is printed only as significance stars (p<=0.01); SE(theta) cannot be recovered from a threshold.")

# ---------------- S4 Bono ----------------
src = "arxiv-2511.13860v1"
ln_t3, _ = find(src, r"^Table 3: Comparison of Email Processing Times")
ln_aw, aw = find(src, r"^Aware Time per Email", ln_t3)
ln_co, cw = find(src, r"^Control Time per Email", ln_t3)
ln_ts, ts = find(src, r"^% Time Savings", ln_t3)
ln_p, pv = find(src, r"^p-value", ln_t3)
find(src, r"Control \( N\{=\}52")
find(src, r"Aware \( N\{=\}56")
last = lambda s: s.split("|")[-1].strip()
aw4, co4, sav4, p4 = last(aw), last(cw), last(ts), last(pv)
c4 = -float(sav4) / 100
th4 = math.log(1 + c4)
z4 = norm.ppf(1 - float(p4) / 2)
se4 = abs(th4) / z4
items["S4"] = dict(
    estimate_as_printed=f"Aware {aw4}, Control {co4} (row labels: Time per Email (s)); % Time Savings {sav4}",
    uncertainty_as_printed=f"p-value {p4}",
    sample_as_printed="Aware N=56, Control N=52 (of 167 analysts in three arms; 13 low-scoring responses filtered)",
    definition="Percentage time saving per processed email (Full Corpus, whole 25-email queue without resolve-benign protocol), Aware arm (agent verdicts, rationale, agent-ordered queue) vs Control (random order, no AI), randomised",
    locator=f"{src}; Section 3.3, Table 3 (L{ln_t3}-L{ln_p}), column 'Full Corpus', rows Aware/Control Time per Email, % Time Savings, p-value",
    harmonisable=True, harmonised_estimate=th4, harmonised_se=se4, poolable=True,
    transformation=f"percentage time saving {sav4}% -> c = {c4}; theta = ln(1 + ({c4})) = {th4:.6f}; APPROXIMATE SE from two-sided p: SE = |theta|/z_(1-{p4}/2) = {abs(th4):.6f}/{z4:.6f} = {se4:.6f} (ratio of displayed means ln({aw4}/{co4}) = {math.log(float(aw4)/float(co4)):.6f}, sensitivity only)",
    reason_if_not=None)

# ---------------- S5 Kreitmeir & Raschky ----------------
src = "arxiv-2304.09339v1"
ln_tp, tp = find(src, r"TreatPost \|")
ln_tse, tse = find(src, r"^\| \(", ln_tp)
ln_N, nl = find(src, r"^N \|", ln_tp)
find(src, r"aggregate Output measure")
b5 = [c.strip() for c in tp.split("|")][3]
s5 = [c.strip() for c in tse.split("|")][3]
n5 = [c.strip() for c in nl.split("|")][3]
items["S5"] = dict(
    estimate_as_printed=f"{b5} (Treat x Post, Output, panel A, column (3))",
    uncertainty_as_printed=f"SE {s5} (clustered at user level)",
    sample_as_printed=f"N = {n5} (user-day observations, users with affiliation)",
    definition="Difference-in-differences coefficient (Italy vs Austria/France, post vs pre ChatGPT ban) on an aggregate count of GitHub Output events (Push, Pull Requests, comments, Create, Issues)",
    locator=f"{src}; Section 4, Table 1 (L{ln_tp}-L{ln_N}), panel A 'Pre 27/28.03 - Post 03/04.04', row Treat x Post, column (3) Output (users with affiliation)",
    harmonisable=False, harmonised_estimate=None, harmonised_se=None, poolable=False,
    transformation="none: no transformation maps an effect on output counts onto ln(T_AI/T_noAI)",
    reason_if_not="Not harmonisable: the effect is on output counts (GitHub events) around the withdrawal of ChatGPT access in a non-randomised difference-in-differences design, not a comparison of task-completion times with and without AI on the same tasks.")

# ---------------- Pooling ----------------
def pool(est, se):
    est, se = np.asarray(est, float), np.asarray(se, float)
    w = 1 / se**2
    fe = float(np.sum(w * est) / np.sum(w)); fe_se = float(math.sqrt(1 / np.sum(w)))
    q = float(np.sum(w * (est - fe) ** 2)); df = len(est) - 1
    C = float(np.sum(w) - np.sum(w**2) / np.sum(w))
    tau2 = max(0.0, (q - df) / C) if C > 0 else 0.0
    i2 = max(0.0, (q - df) / q) * 100 if q > 0 else 0.0
    wr = 1 / (se**2 + tau2)
    re_ = float(np.sum(wr * est) / np.sum(wr)); re_se = float(math.sqrt(1 / np.sum(wr)))
    return dict(fe=fe, fe_se=fe_se, q=q, df=df, tau2=tau2, i2=i2, re=re_, re_se=re_se)

pool_ids = [k for k in ids if items[k]["poolable"]]
P = pool([items[k]["harmonised_estimate"] for k in pool_ids], [items[k]["harmonised_se"] for k in pool_ids])
interp = "low" if P["i2"] < 25 else ("moderate" if P["i2"] < 50 else "substantial")
re_lo, re_hi = P["re"] - Z * P["re_se"], P["re"] + Z * P["re_se"]

# sensitivity (not headline)
S_ci = pool([b1, th4], [se1_ci, se4])
th4_ratio = math.log(float(aw4) / float(co4)); se4_ratio = abs(th4_ratio) / z4
S_ratio = pool([b1, th4_ratio], [se1, se4_ratio])
harm = [items[k]["harmonised_estimate"] for k in ids if items[k]["harmonisable"]]

data_issues = [
 {"issue": "S1 (Paradis): Table III states n = 96, while Section VI-B reports N=47 (AI) + N=46 (no-AI) = 93 completed tasks.", "handling": "Both sample statements recorded in sample_as_printed; coefficient and SE taken from Table III as designated; no effect on theta."},
 {"issue": "S1 (Paradis): Section VI-B says the AI group (Mean LogToT=4.6) was faster than the no-AI group (Mean LogToT=4.3), which contradicts the raw means (96 vs 114 min); the log means appear swapped.", "handling": "Not used; harmonisation relies on the Model 2 regression coefficient designated in studies.json."},
 {"issue": "S1 (Paradis): printed SE 0.14 and printed 95% CI [-0.51, 0.03] imply slightly different SEs (0.14 vs 0.1378); '21% faster' denotes a 21% reduction in time (1-exp(beta)).", "handling": "Printed SE used (convention for a coefficient with SE); CI-based SE reported as sensitivity_ci_se_s1_* keys."},
 {"issue": "S2 (Becker): the confidence interval for the 19% estimate is shown only in Figure 1/Figure 15, which are not reproduced in the extract.", "handling": "No SE imputed; S2 harmonised (theta = ln 1.19) but not poolable."},
 {"issue": "S3 (Bono & Xu): uncertainty only as significance stars with a non-standard scheme (* 0.05<p<=0.1, ** 0.01<p<=0.05, *** p<=0.01); displayed minutes give 44.29/63.10-1 = -29.81% vs printed -29.79%; control subjects could use other AI tools (25% did).", "handling": "Printed -29.79% used for theta; S3 not poolable; contamination noted as a limitation."},
 {"issue": "S4 (Bono): printed % Time Savings 11.1 does not equal 1-50.5/56.9 = 11.25% (Malicious and Non-Malicious columns also do not reproduce from displayed times); row labels say seconds per email, but Full Corpus values (50.5, 56.9) are inconsistent with the per-email Malicious/Non-Malicious values (2.8/1.9, 1.8/2.4), so the unit of the Full Corpus column is ambiguous.", "handling": "Printed % saving (11.1) used for theta (unit-free); ratio of displayed means used as sensitivity (sensitivity_ratio_s4_* keys); SE recovered approximately from p = 0.471."},
 {"issue": "S4 (Bono): the Aware arm differs from Control by both agent verdicts/rationales and agent-prioritised queue ordering; only a p-value is printed.", "handling": "Treated as AI vs no-AI; SE flagged as approximate."},
 {"issue": "S5 (Kreitmeir & Raschky): table notes print '** p<0.5' (presumably 0.05); row label rendered as '\\times TreatPost'; date line under the author block is a rendering artefact; whether Output is logged is ambiguous in Section 3.1.", "handling": "Values copied verbatim; S5 not harmonisable regardless of these issues."},
 {"issue": "Stray 'p<', 'p=' and footnote-number lines in the extracts.", "handling": "Ignored as rendering artefacts per data/README.md."},
]

res = {
 "task_id": "T3-synthesis",
 "target_definition": "theta = ln(T_AI / T_noAI): natural log of the ratio of mean task-completion time with access to a generative-AI assistant to mean completion time without it, for the same tasks, from a randomised comparison; negative values mean faster with AI.",
 "items": items,
 "k_extracted": sum(1 for k in ids if items[k]["estimate_as_printed"]),
 "k_harmonised": sum(1 for k in ids if items[k]["harmonised_estimate"] is not None),
 "k_pooled": len(pool_ids),
 "pooled_study_ids": pool_ids,
 "fe_estimate": P["fe"], "fe_se": P["fe_se"], "fe_ci_low": P["fe"] - Z * P["fe_se"], "fe_ci_high": P["fe"] + Z * P["fe_se"],
 "re_estimate": P["re"], "re_se": P["re_se"], "re_ci_low": re_lo, "re_ci_high": re_hi,
 "q_stat": P["q"], "q_df": P["df"], "i2": P["i2"], "tau2": P["tau2"],
 "heterogeneity_interpretation": interp,
 "pooled_significant_5pct": bool(re_lo > 0 or re_hi < 0),
 "pooled_pct_change_re": 100 * (math.exp(P["re"]) - 1),
 "pooled_pct_change_re_ci_low": 100 * (math.exp(re_lo) - 1),
 "pooled_pct_change_re_ci_high": 100 * (math.exp(re_hi) - 1),
 "baseline_unweighted_mean_theta_harmonised": float(np.mean(harm)),
 "baseline_unweighted_mean_theta_pooled": float(np.mean([items[k]["harmonised_estimate"] for k in pool_ids])),
 "sensitivity_ci_se_s1_re_estimate": S_ci["re"], "sensitivity_ci_se_s1_re_se": S_ci["re_se"],
 "sensitivity_ratio_s4_theta": th4_ratio, "sensitivity_ratio_s4_re_estimate": S_ratio["re"], "sensitivity_ratio_s4_re_se": S_ratio["re_se"],
 "feasibility": "complete",
 "feasibility_reason": None,
 "data_issues": data_issues,
 "n_input_rows": n_input_rows,
 "n_analysed": len(ids),
 "random_seed": None,
}
os.makedirs(OUT, exist_ok=True)
with open(os.path.join(OUT, "results.json"), "w") as fh:
    json.dump(res, fh, indent=1)

cols = ["study", "estimate_as_printed", "uncertainty_as_printed", "sample_as_printed", "definition", "locator",
        "harmonised_estimate", "harmonised_se", "harmonisable", "poolable", "transformation", "reason_if_not"]
with open(os.path.join(OUT, "table1.csv"), "w", newline="") as fh:
    w = csv.writer(fh); w.writerow(cols)
    for k in ids:
        it = items[k]
        w.writerow([k] + [("" if it[c] is None else it[c]) for c in cols[1:]])

# ---------------- Forest plot ----------------
short = {s["id"]: s["short"].split(",")[0] for s in studies}
hk = [k for k in ids if items[k]["harmonisable"]]
fig, ax = plt.subplots(figsize=(8, 4.2))
ys = list(range(len(hk) + 1, 1, -1))
for y, k in zip(ys, hk):
    t, s = items[k]["harmonised_estimate"], items[k]["harmonised_se"]
    if s is not None:
        ax.errorbar(t, y, xerr=Z * s, fmt="s", color="k", capsize=3)
        lab = f"{t:.3f} [{t-Z*s:.3f}, {t+Z*s:.3f}]" + (" (approx. SE)" if k == "S4" else "")
    else:
        ax.plot(t, y, "o", mfc="white", mec="k")
        lab = f"{t:.3f} (no SE; not pooled)"
    ax.text(0.45, y, lab, va="center", fontsize=8)
ax.fill([re_lo, P["re"], re_hi, P["re"]], [1, 1.25, 1, 0.75], color="tab:blue")
ax.text(0.45, 1, f"RE (k={len(pool_ids)}): {P['re']:.3f} [{re_lo:.3f}, {re_hi:.3f}]", va="center", fontsize=8)
ax.axvline(0, color="grey", ls="--", lw=0.8)
ax.set_yticks(ys + [1]); ax.set_yticklabels([f"{k} {short[k]}" for k in hk] + ["Pooled DL random effects"])
ax.set_xlim(-0.9, 1.3); ax.set_xlabel("theta = ln(T_AI / T_noAI)   (negative = faster with AI)")
ax.set_title("Forest plot: harmonised effects of AI access on task-completion time\n(S5 not harmonisable; S2, S3 shown without intervals)", fontsize=9)
fig.tight_layout(); fig.savefig(os.path.join(OUT, "figure1.png"), dpi=150); plt.close(fig)

print("=== Results ===")
for k in ids:
    print(k, "theta =", items[k]["harmonised_estimate"], "SE =", items[k]["harmonised_se"], "poolable =", items[k]["poolable"])
print(f"FE {res['fe_estimate']:.4f} (SE {res['fe_se']:.4f}); RE {res['re_estimate']:.4f} (SE {res['re_se']:.4f}) "
      f"[{re_lo:.4f}, {re_hi:.4f}]; Q={P['q']:.4f}, I2={P['i2']:.2f}%, tau2={P['tau2']:.4f}; pct change RE {res['pooled_pct_change_re']:.2f}%")