"""Harmonise five printed effect sizes (F4-conf-03). Reads only data/, writes out/results.json.
Run from workspace root: ~/wp29/venv/bin/python out/analysis.py"""
import json, math, os, re, csv

ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
D = os.path.join(ROOT, "data")
SRC = os.path.join(D, "sources")
OUT = os.path.join(ROOT, "out")

def read_lines(path):
    with open(path, encoding="utf-8") as f:
        return f.read().split("\n")

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

def find_row(lines, pattern, fname):
    """Return (line_number, cells) of the unique line matching regex pattern; print it."""
    hits = [(i + 1, l) for i, l in enumerate(lines) if re.match(pattern, l)]
    assert len(hits) == 1, f"pattern {pattern!r} matched {len(hits)} lines in {fname}"
    n, l = hits[0]
    print(f"[{fname}:{n}] {l}")
    return n, [c.strip() for c in l.split("|")]

def num(cell):
    """Parse a printed number, stripping significance stars, %, 'pp' and thousands separators."""
    m = re.match(r"^\s*([-+]?\d[\d,]*\.?\d*)", cell.replace("*", ""))
    assert m, cell
    return float(m.group(1).replace(",", ""))

effects = json.load(open(os.path.join(D, "effects.json"), encoding="utf-8"))
files = {
    "arxiv-2609.27732v1": "arxiv-2609.27732v1.txt",
    "arxiv-2606.14621v1": "arxiv-2606.14621v1.txt",
    "arxiv-2505.04555v3": "arxiv-2505.04555v3.txt",
    "arxiv-2510.03658v1": "arxiv-2510.03658v1.txt",
}
L = {k: read_lines(os.path.join(SRC, v)) for k, v in files.items()}
H = {}

# ---- E1: DCR overall ATT, Table 6, 5 km -> relative % change ------------------------
s = "arxiv-2609.27732v1"
_, hdr = find_row(L[s], r"^\| 1 km \| 3 km \| 5 km$", files[s])
n_att, att = find_row(L[s], r"^ATT \| -0\.399", files[s])
n_pct, pct = find_row(L[s], r"^% change \| -32\.9%", files[s])
col = hdr.index("5 km")  # header has a leading empty cell, same layout as data rows
b1 = num(att[col])
v1 = 100 * (math.exp(b1) - 1)
H["E1"] = dict(value=round(v1, 4), unit="percent (relative change vs. counterfactual hospitalisation rate)",
               transformation=f"100*(exp(b)-1) with b = {att[col].replace('*','')} (printed ATT, log points) -> 100*(exp({b1})-1) = {v1:.4f}%",
               harmonisable=True, missing_quantity=None,
               locator=f"arxiv-2609.27732v1 | Table 6 | row 'ATT' | column '5 km' (line {n_att}); cross-check: row '% change' | column '5 km' = {pct[col]} (line {n_pct})")
print(f"E1: b={b1} -> {v1:.4f}% (printed % change {pct[col]})")

# ---- E4: DCR fatal, Table 7, 5 km -> absolute change per 100,000 --------------------
_, hdr7 = find_row(L[s], r"^Outcome \| 1 km \| 3 km \| 5 km$", files[s])
n_fat, fat = find_row(L[s], r"^Drug-related \(fatal\) \|", files[s])
b4 = num(fat[hdr7.index("5 km")])
rel4 = 100 * (math.exp(b4) - 1)
# Search the extract for any printed baseline/mean level of the fatal rate.
lvl = [l for l in L[s] if re.search(r"(mean|baseline|pre-treatment level|average rate)", l, re.I)
       and re.search(r"fatal", l, re.I) and re.search(r"\d+(\.\d+)? (per|hospitali[sz]ations per)", l, re.I)]
print(f"E4: b={b4}; relative change would be {rel4:.4f}%; lines printing a fatal-rate level: {len(lvl)}")
H["E4"] = dict(value=None, unit="hospitalisations per 100,000 inhabitants per year",
               transformation=(f"absolute change = R0*(exp(b)-1) with b = {b4} (printed ATT, log points); "
                               f"exp(b)-1 = {rel4/100:.4f} is computable, but R0 (counterfactual/baseline fatal rate) is not printed"),
               harmonisable=False,
               missing_quantity=("baseline (counterfactual, e.g. pre-opening treated-region mean) level of the fatal drug-related "
                                 "hospitalisation rate per 100,000 inhabitants per year in the 5 km treated regions; not printed in the extract"),
               locator=f"arxiv-2609.27732v1 | Table 7 | row 'Drug-related (fatal)' | column '5 km' (line {n_fat})")

# ---- E2: RESEA reemployment ATE (pp) -> % of control mean ---------------------------
s = "arxiv-2606.14621v1"
_, hdr2 = find_row(L[s], r"^Outcome \| Treatment \| Control \| Treatment-control difference", files[s])
n_re, re_row = find_row(L[s], r"^Reemployment \| \d+\.\d%", files[s])
d2 = num(re_row[hdr2.index("Average treatment effect estimate")].split(";")[0])
p2 = num(re_row[hdr2.index("Control")])
v2 = 100 * d2 / p2
H["E2"] = dict(value=round(v2, 4), unit="percent (relative change vs. control-group mean reemployment rate)",
               transformation=f"100*Delta/p with Delta = {d2} pp (printed ATE) and p = {p2}% (printed control mean) -> 100*{d2}/{p2} = {v2:.4f}%",
               harmonisable=True, missing_quantity=None,
               locator=(f"arxiv-2606.14621v1 | Table 2 | row 'Reemployment' | column 'Average treatment effect estimate' (1.47 pp) "
                        f"and column 'Control' (75.0%) (line {n_re})"))
print(f"E2: {d2}/{p2} -> {v2:.4f}%")

# ---- E3: MW employment elasticity * printed % MW change ----------------------------
s = "arxiv-2505.04555v3"
n_el, el = find_row(L[s], r"^Employment elasticity w\.r\.t MW \|", files[s])
n_g, g = find_row(L[s], r"^% MW changes \|", files[s])
n_note, _ = find_row(L[s], r"^Note: The table reports the effects of a minimum wage increase", files[s])
eta, gg = num(el[1]), num(g[1])
v3 = 100 * eta * gg
H["E3"] = dict(value=round(v3, 4), unit="percent (relative change in employment)",
               transformation=f"100*eta*g with eta = {eta} (printed elasticity) and g = {gg} (printed '% MW changes', proportional) -> 100*{eta}*{gg} = {v3:.4f}%",
               harmonisable=True, missing_quantity=None,
               locator=(f"arxiv-2505.04555v3 | Table 2 | row 'Employment elasticity w.r.t MW' | column 'Estimate (SE)' (line {n_el}); "
                        f"row '% MW changes' | column 'Estimate (SE)' (line {n_g})"))
print(f"E3: {eta}*{gg} -> {v3:.4f}%")

# ---- E5: LPM coefficient -> percentage points --------------------------------------
s = "arxiv-2510.03658v1"
_, hdr5 = find_row(L[s], r"^\| \(1\) \| \(2\) \| \(3\) \| \(4\)$", files[s])
n_t, tr = find_row(L[s], r"^Treatment \| ", files[s])
n_m, mean = find_row(L[s], r"^Mean Dep\. Var\. \|", files[s])
n_pc, pc = find_row(L[s], r"^% Change \|", files[s])
c4 = hdr5.index("(4)")
b5 = num(tr[c4])
v5 = 100 * b5
H["E5"] = dict(value=round(v5, 4), unit="percentage points",
               transformation=f"100*beta with beta = {b5} (printed LPM coefficient, 0-1 scale) -> {v5:.4f} pp",
               harmonisable=True, missing_quantity=None,
               locator=f"arxiv-2510.03658v1 | Table 3 | row 'Treatment' | column '(4)' (line {n_t})")
print(f"E5: {b5} -> {v5:.4f} pp (check: beta/mean = {b5/num(mean[c4]):.5f}, printed % Change {pc[c4]})")

# ---- consistency diagnostics used for data_issues ----------------------------------
s = "arxiv-2609.27732v1"
_, a19 = find_row(L[s], r"^ATT \| -0\.431\*\* \| -0\.438", files[s])
_, a19se = find_row(L[s], r"^\| \(0\.218\)", files[s])
_, t6se = find_row(L[s], r"^\| \(0\.412\)", files[s])
_, a19n = find_row(L[s], r"^Note: All specifications use the 5 km buffer", files[s])
abs_line = [l for l in L[s] if "which fall by 35%" in l]
mech_line = [l for l in L[s] if "by approximately 30% at 5 km" in l]
print(f"check E4 text: abstract '35%' present={bool(abs_line)}, Sec 6.1 '30%' present={bool(mech_line)}, exp(-0.354)-1={rel4:.2f}%")

data_issues = [
    {"issue": "arxiv-2609.27732v1: the abstract says fatal drug-related hospitalisations 'fall by 35%%', but Table 7 (5 km) prints -0.354 log points, i.e. 100*(exp(-0.354)-1) = %.1f%%, matching Section 6.1 ('approximately 30%% at 5 km'); 35%% is the overall (Table 6) effect." % rel4,
     "handling": "Table 7 cell used as the printed estimate for E4; abstract figure not used. E4 is not harmonisable for a separate reason (missing baseline level)."},
    {"issue": "arxiv-2609.27732v1: the same CS-DR 5 km ATT (-0.431) is printed with SE (0.210) in Table 6 but (0.218) in Appendix Table A19; A19 note says 640 never-treated regions while its table and Table 6 say 608.",
     "handling": "Only the point estimate (identical in both tables) is used; SEs and region counts are not needed and were not used."},
    {"issue": "arxiv-2505.04555v3: Table 2 prints '% MW changes' = 0.047 (4.7%), while the table note describes the 2023 increases as 'about 4%'.",
     "handling": "Per the target definition, the table value 0.047 is used for E3; with 'about 4%%' the implied change would be 100*%s*0.04 = %.4f%% (sensitivity only, not the reported value)." % (eta, 100 * eta * 0.04)},
    {"issue": "arxiv-2505.04555v3: Table 2 prints the same SE (0.091) for 'Affected employment' and 'Emp. elasticity w.r.t affected wage' (possible duplication); point estimates are internally consistent (-0.268/0.366 = -0.732).",
     "handling": "Not used; no uncertainty is propagated."},
    {"issue": "arxiv-2510.03658v1: Table 3 prints '%% Change' = -0.3685 for column (4), whereas -0.0962/0.261 = %.5f (rounds to -0.3686); the text also attributes the 0.0855 estimate to 'table 2' and calls an 18.79%% relative change a 'percentage change' while interpreting it as lives saved per 100 deaths." % (b5 / num(mean[c4])),
     "handling": "Not material for E5, whose target (percentage points) needs only the printed coefficient -0.0962; noted for transparency."},
    {"issue": "arxiv-2606.14621v1: the ATE (1.47 pp) is block(week)-adjusted while the 75.0% control mean is the raw aggregate control-group mean; the abstract rounds the ATE to 1.5 pp.",
     "handling": "Used the printed control mean 75.0% as the counterfactual level p, as the target definition specifies (percent of the control-group mean); 1.47 (table) used rather than 1.5 (abstract)."},
]

order = ["E1", "E2", "E3", "E4", "E5"]
assert [e["id"] for e in effects] == order
harm = {k: H[k] for k in order}
n_h = sum(1 for k in order if harm[k]["value"] is not None)
n_nh = sum(1 for k in order if harm[k]["harmonisable"] is False)
n_input_rows = {"effects.json": len(effects), "README.md": wc_l(os.path.join(D, "README.md"))}
for v in files.values():
    n_input_rows[v] = wc_l(os.path.join(SRC, v))

res = {
    "task_id": "F4-conf-03",
    "harmonised": harm,
    "n_harmonised": n_h,
    "n_not_harmonisable": n_nh,
    "feasibility": "complete",
    "feasibility_reason": None,
    "data_issues": data_issues,
    "n_input_rows": n_input_rows,
    "n_analysed": len(effects),
    "random_seed": None,
    # P2 baseline: naive readings without the proper transformation (for comparison only)
    "baseline_naive": {
        "E1_logpoints_times_100_pct": round(100 * b1, 4),
        "E2_printed_pp": d2,
        "E3_printed_elasticity": eta,
        "E4_logpoints_times_100_pct": round(100 * b4, 4),
        "E5_printed_coefficient": b5,
        "note": "Naive values = printed numbers read at face value (log points x100 treated as percent; no rescaling). Not harmonised estimates.",
    },
}
os.makedirs(OUT, exist_ok=True)
with open(os.path.join(OUT, "results.json"), "w", encoding="utf-8") as f:
    json.dump(res, f, indent=2, ensure_ascii=False)
with open(os.path.join(OUT, "exclusions.csv"), "w", newline="", encoding="utf-8") as f:
    w = csv.writer(f)
    w.writerow(["step", "rule", "n_before", "n_after", "n_removed", "rationale"])
    w.writerow([1, "keep all items listed in data/effects.json", len(effects), len(effects), 0,
                "No item excluded; E4 is analysed and reported as not harmonisable rather than removed."])
print(f"n_harmonised={n_h} n_not_harmonisable={n_nh} feasibility=complete -> out/results.json written")