"""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/. Deterministic (no randomness; random_seed = null)."""
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__)))
DATA = os.path.join(ROOT, "data")
OUT = os.path.join(ROOT, "out")
Z = 1.96

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

def find(src, pattern, nth=0):
    """Return (line_no, line) of the nth line matching pattern; print it."""
    L = lines_of(src)
    hits = [(i + 1, l) for i, l in enumerate(L) if re.search(pattern, l)]
    if len(hits) <= nth:
        raise RuntimeError(f"pattern {pattern!r} not found in {src}")
    ln, line = hits[nth]
    show = line if len(line) < 400 else "..." + line[max(0, re.search(pattern, line).start() - 150): re.search(pattern, line).end() + 150] + "..."
    print(f"[{src} L{ln}] {show}")
    return ln, line

def wc_l(path):
    with open(path, "rb") as fh:
        return fh.read().count(b"\n")

studies = json.load(open(os.path.join(DATA, "studies.json")))
assert [s["id"] for s in studies] == ["S1", "S2", "S3", "S4", "S5"]
print("== Source lines relied upon ==")

items = {}
data_issues = []

# ---------------- S1 Paradis et al. 2024 ----------------
src = "arxiv-2410.12944v3"
ln_c, l_c = find(src, r"ExpCon \|")
ln_s = ln_c + 2
se_line = lines_of(src)[ln_s - 1]; print(f"[{src} L{ln_s}] {se_line}")
coefs = re.findall(r"-?\d+\.\d+", l_c); ses = re.findall(r"\((\d+\.\d+)\)", se_line)
b1, se1 = float(coefs[1]), float(ses[1])          # Model 2 column
assert (b1, se1) == (-0.24, 0.14)
ln_ci, l_ci = find(src, r"best-fit Model 2")
m = re.search(r"\[(-?\d+\.\d+),\s*(-?\d+\.\d+)\]\s*,\s*p=(\d*\.\d+)", l_ci)
ci_l, ci_u, p1 = float(m.group(1)), float(m.group(2)), float(m.group(3))
ln_n, l_n = find(src, r"TABLE III")
ln_m, l_m = find(src, r"Mean LogToT")
se1_ci = (ci_u - ci_l) / (2 * Z)
items["S1"] = dict(
    estimate_as_printed="ExpCon coefficient, Model 2 (Dev): -0.24 (marked '.'); quoted in VI-C as beta_1=-0.24, i.e. AI effect (1-exp(beta_1))*100 = 21%",
    uncertainty_as_printed="SE (0.14); 95% CI = [-0.51, 0.03]; p=0.086 (NS)",
    sample_as_printed="n = 96 (Table III caption); VI-B: N=47 (AI), N=46 (no-AI)",
    definition="OLS coefficient on experimental condition (AI vs no-AI) with developer-level covariates; dependent variable = natural log of time on task (LogToT); randomised enterprise RCT, one complex coding task",
    locator=f"{src}, Table III (L{ln_n}), row ExpCon (L{ln_c}) and SE row (L{ln_s}), column (2) 'Dev' = Model 2; CI and p-value in Section VI-C (L{ln_ci})",
    harmonisable=True, harmonised_estimate=b1, harmonised_se=se1, poolable=True,
    transformation=f"log-time coefficient: theta = b = -0.24; SE(theta) = SE(b) = 0.14 (CI-implied alternative (0.03-(-0.51))/(2*1.96) = {se1_ci:.4f})",
    reason_if_not=None)
data_issues.append({"issue": f"S1: Section VI-B (L{ln_m}) prints Mean LogToT = 4.6 for the AI group and 4.3 for the no-AI group while stating the AI group was faster (96 vs 114 min; ln 96 = 4.56, ln 114 = 4.74): internally inconsistent.",
                    "handling": "Mean LogToT values not used; theta taken from the Table III Model 2 coefficient."})
data_issues.append({"issue": "S1: Section VI-B reports N=47+46=93 developers while Table III states n = 96.",
                    "handling": "Both printed as is in sample_as_printed; not needed for theta or SE."})
data_issues.append({"issue": f"S1: the CI [-0.51, 0.03] implies SE = {se1_ci:.4f} whereas Table III prints SE = 0.14 (rounding).",
                    "handling": "Printed SE(b) = 0.14 used (first harmonisation convention; reproduces p = 0.086 under a normal approximation); sensitivity with the CI-implied SE reported as sens_* keys."})

# ---------------- S2 Becker et al. 2025 ----------------
src = "arxiv-2507.09089v2"
ln_a, l_a = find(src, r"increases completion time by 19%")
ln_r, l_r = find(src, r"19% more time on average")
ln_d, l_d = find(src, r"Developers complete 136 issues")
ln_f, l_f = find(src, r"Figure 1, which displays the confidence interval")
c2 = 0.19
items["S2"] = dict(
    estimate_as_printed="+19% (\"allowing AI actually increases completion time by 19%\"; \"they implement issues in 19% more time on average\")",
    uncertainty_as_printed="none printed in the extract (95% CI shown only in Figure 1/Figure 15, not reproduced)",
    sample_as_printed="16 developers, 246 tasks (136 issues AI-allowed, 110 AI-disallowed)",
    definition="S = E[time with AI allowed]/E[time with 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}, Abstract (L{ln_a}); Section 3 Results (L{ln_r}); estimator Sections 2.3 and D.1-D.3",
    harmonisable=True, harmonised_estimate=math.log(1 + c2), harmonised_se=None, poolable=False,
    transformation="percentage change c = +0.19: theta = ln(1 + 0.19) = ln(1.19)",
    reason_if_not="Not poolable: no SE, CI or p-value for the 19% estimate is printed in the extract; its confidence interval appears only in figures that are not reproduced.")

# ---------------- S3 Bono & Xu 2024 ----------------
src = "arxiv-2411.01067v2"
ln_t, l_t = find(src, r"^All Content \| 53\.84")
ln_cap, l_cap = find(src, r"^Table 3: Comparison of task duration")
ln_ed, l_ed = find(src, r"^Editorial note: the paper marks significance")
cells = [c.strip() for c in l_t.split("|")]
overall, treat3, ctrl3 = float(cells[1]), float(cells[2]), float(cells[3])
imp = cells[4]; c3 = float(re.match(r"(-?\d+\.\d+)%", imp).group(1)) / 100
assert (treat3, ctrl3, c3) == (44.29, 63.10, -0.2979)
items["S3"] = dict(
    estimate_as_printed="All Content: Overall 53.84 | Treatment 44.29 | Control 63.10 (minutes) | Improvement -29.79%***",
    uncertainty_as_printed="significance stars only: *** (0<p<=0.01); no SE or CI printed",
    sample_as_printed="181 subjects completed the task satisfactorily (one further subject dropped)",
    definition="Difference in mean task duration (minutes), Copilot vs control, holding accuracy constant (control duration predicted at the treatment group's accuracy, Appendix B); randomised assignment",
    locator=f"{src}, Section 4.2, Table 3 (L{ln_cap}), 'Aggregate Results' block, row 'All Content' (L{ln_t}), columns Treatment/Control/Improvement",
    harmonisable=True, harmonised_estimate=math.log(1 + c3), harmonised_se=None, poolable=False,
    transformation=f"percentage change c = -0.2979: theta = ln(1 - 0.2979) = ln(0.7021) (ratio of printed durations 44.29/63.10 would give {math.log(treat3/ctrl3):.5f})",
    reason_if_not="Not poolable: uncertainty is reported by significance stars only (***), so SE(theta) cannot be obtained.")
data_issues.append({"issue": f"S3: printed durations give 44.29/63.10 - 1 = {100*(treat3/ctrl3-1):.2f}% whereas the printed improvement is -29.79%.",
                    "handling": "Printed improvement c = -0.2979 used (theta = ln(1+c)); the ratio-based value differs by < 0.001 on the theta scale; S3 is not pooled anyway."})

# ---------------- S4 Bono 2025 ----------------
src = "arxiv-2511.13860v1"
ln_h, l_h = find(src, r"^Metric \| Malicious")
ln_aw, l_aw = find(src, r"^Aware Time per Email")
ln_co, l_co = find(src, r"^Control Time per Email")
ln_ts, l_ts = find(src, r"^% Time Savings")
ln_p, l_p = find(src, r"^p-value \|")
ln_arms, l_arms = find(src, r"Aware \( N\{=\}56")
last = lambda l: float(l.split("|")[-1].strip())
aware4, ctrl4, sav4, p4 = last(l_aw), last(l_co), last(l_ts), last(l_p)
assert (aware4, ctrl4, sav4, p4) == (50.5, 56.9, 11.1, 0.471)
c4 = -sav4 / 100
th4 = math.log(1 + c4)
z4 = norm.ppf(1 - p4 / 2)
se4 = abs(th4) / z4
items["S4"] = dict(
    estimate_as_printed="Full Corpus: Aware Time per Email 50.5 s | Control Time per Email 56.9 s | % Time Savings 11.1",
    uncertainty_as_printed="p-value 0.471",
    sample_as_printed="Aware N=56, Control N=52 (Blind N=59; 167 analysts in total, 13 low-scoring responses filtered); 25-email queue per subject",
    definition="Mean time per email (seconds), Aware arm (agent verdicts, rationale and agent-ordered queue) vs Control arm (no AI), full queue; % time savings = reduction relative to Control; randomised three-arm RCT",
    locator=f"{src}, Section 3.3, Table 3 (L{ln_h}), rows Aware/Control Time per Email, % Time Savings, p-value (L{ln_aw}-L{ln_p}), column 'Full Corpus'",
    harmonisable=True, harmonised_estimate=th4, harmonised_se=se4, poolable=True,
    transformation=f"time saving 11.1% -> c = -0.111: theta = ln(1 - 0.111) = ln(0.889) = {th4:.5f}; approximate SE from two-sided p: SE = |theta| / z_(1-0.471/2) = {abs(th4):.5f} / {z4:.5f} = {se4:.5f} (APPROXIMATE SE)",
    reason_if_not=None)
data_issues.append({"issue": f"S4: printed Full Corpus times give 1 - 50.5/56.9 = {100*(1-aware4/ctrl4):.2f}% savings, not the printed 11.1%; the Malicious/Non-Malicious columns (2.8, 1.8, 1.9, 2.4, labelled seconds) do not match the full-corpus magnitude (50.5/56.9 s) and their implied savings differ from the printed -52.9/23.6 (ambiguous units/rounding).",
                    "handling": "Only the Full Corpus column is used, with the printed % time savings (11.1) and p-value; the ambiguous columns are not used."})
data_issues.append({"issue": "S4: only a p-value is printed; the SE is recovered as |theta|/z_(1-p/2), which assumes a normal test statistic on the log-ratio scale (the paper's test may be on another scale).",
                    "handling": "S4 pooled with an explicitly approximate SE, as allowed by the task conventions; flagged in report."})

# ---------------- S5 Kreitmeir & Raschky 2023 ----------------
src = "arxiv-2304.09339v1"
ln_hd, l_hd = find(src, r"^\| Events \| Release \| Output")
ln_tp, l_tp = find(src, r"TreatPost \| 0\.115")
ln_tps = ln_tp + 2; l_tps = lines_of(src)[ln_tps - 1]; print(f"[{src} L{ln_tps}] {l_tps}")
ln_N, l_N = find(src, r"^N \| 18332")
ln_leg, l_leg = find(src, r"\*\* p<0\.5")
c5 = [c.strip() for c in l_tp.split("|")][3]; s5 = [c.strip() for c in l_tps.split("|")][3]; n5 = [c.strip() for c in l_N.split("|")][3]
assert (c5, s5, n5) == ("-0.278", "(0.295)", "18332")
items["S5"] = dict(
    estimate_as_printed="Treat x Post (Output, column 3, panel A, users with affiliation): -0.278",
    uncertainty_as_printed="SE (0.295), no significance stars",
    sample_as_printed="N = 18332 (user-day observations)",
    definition="Difference-in-differences coefficient on an aggregate GitHub Output count measure (Push, Pull Requests, PR comments, commit comments, Create, Issues) for Italian vs Austrian/French users, pre (27/28.03.2023) vs post (03/04.04.2023) the ChatGPT ban; user FE, user-clustered SE",
    locator=f"{src}, Section 4, Table 1 (header L{ln_hd}), panel A, row Treat x Post (L{ln_tp}) and SE row (L{ln_tps}), column (3) Output; N row (L{ln_N})",
    harmonisable=False, harmonised_estimate=None, harmonised_se=None, poolable=False,
    transformation="none applied (not harmonisable)",
    reason_if_not="Not harmonisable: the printed effect is on output counts (activity volume) around a withdrawal of access in an observational DiD, not a comparison of task-completion times with and without AI on the same tasks.")
data_issues.append({"issue": f"S5: Table 1 notes print the significance legend as '** p<0.5' (L{ln_leg}), presumably a typo for p<0.05.",
                    "handling": "Irrelevant to the analysis (S5 not harmonisable); recorded."})

# ---------------- pooling ----------------
def pool(th, se):
    th, se = np.asarray(th, float), np.asarray(se, float)
    w = 1 / se**2
    fe = float(np.sum(w * th) / np.sum(w)); fe_se = float(math.sqrt(1 / np.sum(w)))
    k = len(th); df = k - 1
    Q = float(np.sum(w * (th - fe) ** 2))
    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 * th) / np.sum(wr)); re_se = float(math.sqrt(1 / np.sum(wr)))
    return dict(fe=fe, fe_se=fe_se, re=re_, re_se=re_se, Q=Q, df=df, C=C, tau2=tau2, i2=i2,
                w_fe=(w / w.sum()).tolist(), w_re=(wr / wr.sum()).tolist())

pool_ids = [k for k in items if items[k]["poolable"]]
P = pool([items[k]["harmonised_estimate"] for k in pool_ids], [items[k]["harmonised_se"] for k in pool_ids])
# sensitivity: S1 with CI-implied SE
Ps = pool([items[k]["harmonised_estimate"] for k in pool_ids],
          [se1_ci if k == "S1" else items[k]["harmonised_se"] for k in pool_ids])
harm_ids = [k for k in items if items[k]["harmonised_estimate"] is not None]
i2 = P["i2"]
het = "low" if i2 < 25 else ("moderate" if i2 < 50 else "substantial")
re_lo, re_hi = P["re"] - Z * P["re_se"], P["re"] + Z * P["re_se"]

n_input_rows = {"studies.json": len(studies)}
for s in studies:
    fn = s["source"] + ".txt"
    n_input_rows[fn] = wc_l(os.path.join(DATA, "sources", fn))

results = {
    "task_id": "T3-synthesis",
    "target_definition": "theta = ln(T_AI / T_noAI): natural log of the ratio of mean task-completion time with generative-AI assistance to mean completion time without it, for the same tasks, from a randomised comparison (negative = faster with AI).",
    "items": items,
    "k_extracted": sum(1 for k in items if items[k]["estimate_as_printed"]),
    "k_harmonised": len(harm_ids),
    "k_pooled": len(pool_ids),
    "pooled_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"], "q_pvalue": float(1 - __import__("scipy.stats").stats.chi2.cdf(P["Q"], P["df"])),
    "i2": i2, "tau2": P["tau2"],
    "heterogeneity_interpretation": het,
    "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": [100 * (math.exp(re_lo) - 1), 100 * (math.exp(re_hi) - 1)],
    "weights_fe": dict(zip(pool_ids, P["w_fe"])),
    "weights_re": dict(zip(pool_ids, P["w_re"])),
    "baseline_unweighted_mean_theta_harmonised": float(np.mean([items[k]["harmonised_estimate"] for k in harm_ids])),
    "baseline_unweighted_mean_theta_poolable": float(np.mean([items[k]["harmonised_estimate"] for k in pool_ids])),
    "baseline_vote_count_faster_with_ai": sum(1 for k in harm_ids if items[k]["harmonised_estimate"] < 0),
    "sens_re_estimate_s1_ci_se": Ps["re"], "sens_re_se_s1_ci_se": Ps["re_se"],
    "sens_re_ci_low_s1_ci_se": Ps["re"] - Z * Ps["re_se"], "sens_re_ci_high_s1_ci_se": Ps["re"] + Z * Ps["re_se"],
    "feasibility": "complete",
    "feasibility_reason": None,
    "data_issues": data_issues,
    "n_input_rows": n_input_rows,
    "n_analysed": len(items),
    "random_seed": None,
}

with open(os.path.join(OUT, "results.json"), "w") as fh:
    json.dump(results, 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, it in items.items():
        w.writerow([k] + [("" if it[c] is None else it[c]) for c in cols[1:]])

# ---------------- forest plot ----------------
labels = {s["id"]: s["short"].split(",")[0] for s in studies}
fig, ax = plt.subplots(figsize=(8, 4.2))
rows = harm_ids[::-1]
for y, k in enumerate(rows, start=1):
    th, se = items[k]["harmonised_estimate"], items[k]["harmonised_se"]
    if se is not None:
        ax.errorbar(th, y, xerr=Z * se, fmt="s", color="black", capsize=3)
        txt = f"{th:+.3f} [{th-Z*se:+.3f}, {th+Z*se:+.3f}]" + (" (approx. SE)" if k == "S4" else "")
    else:
        ax.plot(th, y, "o", mfc="white", mec="grey")
        txt = f"{th:+.3f} (no SE: not pooled)"
    ax.text(0.62, y, txt, va="center", fontsize=8)
ax.errorbar(P["re"], 0, xerr=Z * P["re_se"], fmt="D", color="tab:red", capsize=3)
ax.text(0.62, 0, f"{P['re']:+.3f} [{re_lo:+.3f}, {re_hi:+.3f}]", va="center", fontsize=8, color="tab:red")
ax.set_yticks(range(0, len(rows) + 1))
ax.set_yticklabels([f"RE pooled (k={len(pool_ids)}: {', '.join(pool_ids)})"] + [f"{k} {labels[k]}" for k in rows])
ax.axvline(0, color="grey", ls="--", lw=0.8)
ax.set_xlim(-0.75, 1.25)
ax.set_xlabel("theta = ln(T_AI / T_noAI)   (negative = faster with AI)")
ax.set_title("Figure 1. Harmonised effects of AI access on task-completion time (95% CI)", fontsize=10)
fig.tight_layout()
fig.savefig(os.path.join(OUT, "figure1.png"), dpi=150)

print("== Results ==")
for k in items:
    it = items[k]
    print(k, "theta=", it["harmonised_estimate"], "se=", it["harmonised_se"], "poolable=", it["poolable"])
print(f"FE {results['fe_estimate']:.4f} (SE {results['fe_se']:.4f}) [{results['fe_ci_low']:.4f}, {results['fe_ci_high']:.4f}]")
print(f"RE {results['re_estimate']:.4f} (SE {results['re_se']:.4f}) [{re_lo:.4f}, {re_hi:.4f}]; Q={P['Q']:.4f} I2={i2:.2f} tau2={P['tau2']:.4f} ({het})")
print(f"pooled % change (RE) = {results['pooled_pct_change_re']:.2f}%; significant at 5%: {results['pooled_significant_5pct']}")