"""F4-conf-03: harmonise five effect sizes from printed source extracts.
Run from workspace root: ~/wp29/venv/bin/python out/analysis.py
Reads only data/; writes out/results.json. Deterministic, no randomness."""
import json, math, os, re

DATA = "data"
SRC = os.path.join(DATA, "sources")


def lines_of(fname):
    with open(os.path.join(SRC, fname), encoding="utf-8") as f:
        return f.read().split("\n")


def find_line(fname, pattern):
    """Return (line_no, text) of the unique line matching regex pattern; print it."""
    L = lines_of(fname)
    hits = [(i + 1, l) for i, l in enumerate(L) if re.search(pattern, l)]
    assert len(hits) == 1, (fname, pattern, hits)
    n, t = hits[0]
    print(f"[{fname}:{n}] {t.strip()[:200]}")
    return n, t


def cells(text):
    return [c.strip() for c in text.split("|")]


def num(s):
    m = re.search(r"-?\d[\d,]*\.?\d*", s)
    return float(m.group(0).replace(",", ""))


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


effects = json.load(open(os.path.join(DATA, "effects.json"), encoding="utf-8"))
print(f"effects.json: {len(effects)} items")
H = {}

# ---- E1: log-point ATT -> percent change ----
f = "arxiv-2609.27732v1.txt"
n_hdr, hdr = find_line(f, r"^\| 1 km \| 3 km \| 5 km$")
n_att, att = find_line(f, r"^ATT \| -0\.399")
n_pct, pct = find_line(f, r"^% change \| -32\.9%")
col = cells(hdr).index("5 km")
b1_s = cells(att)[col].rstrip("*")
b1 = float(b1_s)
v1 = 100 * (math.exp(b1) - 1)
printed_pct = cells(pct)[col]
print(f"E1: b={b1_s}; 100*(exp(b)-1)={v1:.4f}%; printed '% change' = {printed_pct}")
H["E1"] = dict(
    value=round(v1, 2), unit="percent (relative change vs counterfactual rate)",
    transformation=f"100*(exp({b1_s}) - 1) = {v1:.4f}% (exact conversion of a log-point coefficient; matches the printed '% change' {printed_pct})",
    harmonisable=True, missing_quantity=None,
    locator=f"arxiv-2609.27732v1 | Table 6 | row 'ATT' | column '5 km' (extract line {n_att}); cross-check: row '% change', column '5 km' (line {n_pct})")

# ---- E2: pp effect -> percent of control mean ----
f = "arxiv-2606.14621v1.txt"
n_h2, h2 = find_line(f, r"^Outcome \| Treatment \| Control \|")
n_r2, r2 = find_line(f, r"^Reemployment \| 75\.2%")
hc = cells(h2); rc = cells(r2)
ate_s = rc[hc.index("Average treatment effect estimate")].split(";")[0].replace(" pp", "").strip()
ctl_s = rc[hc.index("Control")].replace("%", "").strip()
d2, p2 = float(ate_s), float(ctl_s)
v2 = 100 * d2 / p2
print(f"E2: Delta={ate_s} pp, control mean={ctl_s}% -> {v2:.4f}%")
H["E2"] = dict(
    value=round(v2, 2), unit="percent (relative to the control-group mean reemployment rate)",
    transformation=f"100*{ate_s} pp / {ctl_s}% = {v2:.4f}% (exact given the definition; denominator = raw control-group mean)",
    harmonisable=True, missing_quantity=None,
    locator=f"arxiv-2606.14621v1 | Table 2 | row 'Reemployment' | column 'Average treatment effect estimate' (1.47 pp) and column 'Control' (75.0%) (extract line {n_r2})")

# ---- E3: elasticity x %MW change -> percent change in employment ----
f = "arxiv-2505.04555v3.txt"
n_e3, e3 = find_line(f, r"^Employment elasticity w\.r\.t MW \|")
n_g3, g3 = find_line(f, r"^% MW changes \|")
find_line(f, r"about 4% minimum wage increase")
eta_s = cells(e3)[1]; g_s = cells(g3)[1]
eta, g = float(eta_s), float(g_s)
v3 = 100 * eta * g
print(f"E3: eta={eta_s}, g={g_s} -> {v3:.4f}%")
H["E3"] = dict(
    value=round(v3, 2), unit="percent (relative change in employment)",
    transformation=f"100*({eta_s})*({g_s}) = {v3:.4f}% (first-order/linear approximation; g = printed '% MW changes' = 0.047, i.e. 4.7%)",
    harmonisable=True, missing_quantity=None,
    locator=f"arxiv-2505.04555v3 | Table 2 | row 'Employment elasticity w.r.t MW' | column 'Estimate (SE)' (line {n_e3}); row '% MW changes' | column 'Estimate (SE)' (line {n_g3})")

# ---- E4: log-point ATT -> absolute change per 100,000: needs baseline level ----
f = "arxiv-2609.27732v1.txt"
n_e4, e4 = find_line(f, r"^Drug-related \(fatal\) \|")
b4_s = cells(e4)[3].rstrip("*")
b4 = float(b4_s)
rel4 = 100 * (math.exp(b4) - 1)
print(f"E4: b={b4_s}; relative change would be {rel4:.4f}%; no level of the fatal rate is printed in the extract")
# check that no baseline/mean level of the fatal rate is printed
L = lines_of(f)
lvl = [l for l in L if re.search(r"(mean|baseline|average|pre-treatment).{0,60}(rate|per 100,000)", l, re.I)
       and re.search(r"\d+(\.\d+)? (hospitali[sz]ations )?per 100,000", l)]
print(f"E4: lines printing a numeric level per 100,000: {len(lvl)}")
H["E4"] = dict(
    value=None, unit="hospitalisations per 100,000 inhabitants per year",
    transformation=(f"absolute change = R0*(exp({b4_s}) - 1), where R0 is the counterfactual (pre-treatment/baseline) fatal drug-related "
                    f"hospitalisation rate per 100,000 in the treated (5 km) regions; the relative part is exp({b4_s})-1 = {rel4/100:.4f} "
                    f"({rel4:.2f}%), but R0 is not printed in the extract, so the target cannot be reached"),
    harmonisable=False,
    missing_quantity="baseline (counterfactual / pre-treatment) level of the fatal drug-related hospitalisation rate per 100,000 inhabitants per year in the treated 5 km regions (not printed in the extract)",
    locator=f"arxiv-2609.27732v1 | Table 7 | row 'Drug-related (fatal)' | column '5 km' (extract line {n_e4})")

# ---- E5: LPM coefficient (0-1 scale) -> percentage points ----
f = "arxiv-2510.03658v1.txt"
n_e5, e5 = find_line(f, r"^Treatment \| -0\.0911")
n_m5, m5 = find_line(f, r"^Mean Dep\. Var\. \|")
b5_s = cells(e5)[4]
b5 = float(b5_s)
v5 = 100 * b5
print(f"E5: beta={b5_s} -> {v5:.4f} pp (mean dep var {cells(m5)[4]})")
H["E5"] = dict(
    value=round(v5, 2), unit="percentage points",
    transformation=f"100*({b5_s}) = {v5:.2f} pp (exact rescaling from 0-1 probability scale)",
    harmonisable=True, missing_quantity=None,
    locator=f"arxiv-2510.03658v1 | Table 3 | row 'Treatment' | column '(4)' (extract line {n_e5})")

data_issues = [
    {"issue": "arxiv-2609.27732v1: the abstract states fatal drug-related hospitalisations 'fall by 35%%', but the Table 7 5 km coefficient -0.354 log points implies 100*(exp(-0.354)-1) = %.1f%%; Section 6.1 says 'approximately 30%%'. The abstract figure appears to conflate the overall effect (-35.0%%) or log points with percent." % rel4,
     "handling": "Used the Table 7 coefficient as instructed; E4 is not harmonisable anyway (missing level). Flagged."},
    {"issue": "arxiv-2609.27732v1: the 5 km CS-DR ATT -0.431 has SE (0.210) in Table 6 but (0.218) in Table A19; Table 6/7 report 608 control regions while the A19 note says 640 never-treated comparison regions.",
     "handling": "Point estimates agree (-0.431); SEs not used (no uncertainty propagation). Flagged only."},
    {"issue": "arxiv-2609.27732v1: Section 6 text describes 'five diagnostic sub-categories' while the Figure 7 note mentions four; the Table 7 note lists four ICD sub-outcome groups and gives no definition for the 'Drug-related (fatal)' row.",
     "handling": "Did not affect the E1/E4 harmonisation; flagged."},
    {"issue": "arxiv-2505.04555v3: the Table 2 note describes the 2023 changes as 'about 4% minimum wage increase', while the table row '% MW changes' prints 0.047 (4.7%).",
     "handling": "Used the printed table value 0.047 as the target definition requires ('size of that increase as printed in the same table'); with 0.04 the result would be %.2f%%." % (100 * eta * 0.04)},
    {"issue": "arxiv-2510.03658v1: Section 4.1 refers to 'table 2' and a coefficient -0.0855 with dependent mean 0.455 (headline 18.79%%), which is not in the extract; the Table 3 column (4) '%% Change' prints -0.3685 whereas -0.0962/0.261 = %.4f (rounding difference in the last digit)." % (b5 / 0.261),
     "handling": "Used Table 3 column (4) as specified in effects.json; flagged."},
    {"issue": "arxiv-2606.14621v1: the abstract rounds the reemployment ATE to 1.5 pp versus 1.47 pp in Table 2; the control mean 75.0% is the raw (not block-adjusted) group mean.",
     "handling": "Used the Table 2 printed values 1.47 pp and 75.0% as the target definition specifies the control-group mean."},
]

n_input_rows = {"effects.json": len(effects), "README.md": wc_l(os.path.join(DATA, "README.md"))}
for fn in sorted(os.listdir(SRC)):
    n_input_rows[fn] = wc_l(os.path.join(SRC, fn))

n_h = sum(1 for v in H.values() if v["value"] is not None)
n_nh = sum(1 for v in H.values() if v["harmonisable"] is False)
res = dict(
    harmonised=H, n_harmonised=n_h, n_not_harmonisable=n_nh, task_id="F4-conf-03",
    feasibility="complete",
    feasibility_reason=None,
    data_issues=data_issues, n_input_rows=n_input_rows, n_analysed=len(effects), random_seed=None)
os.makedirs("out", exist_ok=True)
with open("out/results.json", "w", encoding="utf-8") as fh:
    json.dump(res, fh, indent=2, ensure_ascii=False)
print(json.dumps({k: (v["value"], v["harmonisable"]) for k, v in H.items()}))