"""F4-conf-03: harmonise five effect sizes from printed quantities only.
Run from workspace root: ~/wp29/venv/bin/python out/analysis.py
Reads data/ only; prints every source line relied upon; writes out/results.json."""
import json, math, os, re

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


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


def find(fn, pattern, after=None):
    """Return (line_no, text) of the first line matching regex `pattern` (after line `after`)."""
    L = lines(fn)
    for i, t in enumerate(L, 1):
        if after is not None and i <= after:
            continue
        if re.search(pattern, t):
            print(f"[{fn}:{i}] {t.strip()}")
            return i, t
    raise ValueError(f"pattern {pattern!r} not found in {fn}")


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


def num(s):
    return float(re.sub(r"[^0-9.\-]", "", s.replace("−", "-")))


H = {}

# ---- E1: Armendariz Table 6, ATT row, 5 km column (log points -> percent)
f1 = "arxiv-2609.27732v1.txt"
i6, _ = find(f1, r"^Table 6:")
ih, th = find(f1, r"^\| 1 km \| 3 km \| 5 km", after=i6)
ia, ta = find(f1, r"^ATT \|", after=i6)
ip, tp = find(f1, r"^% change \|", after=i6)
col5 = cells(th).index("5 km")  # header has a leading empty cell, same as row label position
b1_str = cells(ta)[col5]
b1 = num(b1_str)
pct_printed = cells(tp)[col5]
v1 = 100 * (math.exp(b1) - 1)
H["E1"] = {
    "value": round(v1, 2),
    "unit": "percent change in the drug-related hospitalisation rate relative to the counterfactual rate",
    "transformation": f"100*(exp(b)-1) with b = {b1_str.rstrip('*')} log points -> 100*(exp(-0.431)-1) = {v1:.2f}% (printed '% change' row: {pct_printed})",
    "harmonisable": True,
    "missing_quantity": None,
    "locator": f"arxiv-2609.27732v1 | Table 6 | row 'ATT' | column '5 km' (extract line {ia}); cross-check: row '% change', column '5 km' (line {ip})",
}

# ---- E2: RESEA Table 2, Reemployment row (pp -> percent of control mean)
f2 = "arxiv-2606.14621v1.txt"
ih2, th2 = find(f2, r"^Outcome \| Treatment \| Control")
ir2, tr2 = find(f2, r"^Reemployment \|", after=ih2)
hdr = cells(th2); row = cells(tr2)
ctrl_str = row[hdr.index("Control")]
ate_cell = row[hdr.index("Average treatment effect estimate")]
ate_str = ate_cell.split(";")[0].strip()
p0 = num(ctrl_str); d2 = num(ate_str)
v2 = 100 * d2 / p0
H["E2"] = {
    "value": round(v2, 2),
    "unit": "percent change in the reemployment probability relative to the control-group mean",
    "transformation": f"100*Delta/p with Delta = {ate_str} and p = {ctrl_str} (control mean) -> 100*1.47/75.0 = {v2:.2f}%",
    "harmonisable": True,
    "missing_quantity": None,
    "locator": f"arxiv-2606.14621v1 | Table 2 | row 'Reemployment' | columns 'Average treatment effect estimate' (Delta) and 'Control' (p) (extract line {ir2})",
}

# ---- E3: Kanayama Table 2 elasticity times printed % MW change
f3 = "arxiv-2505.04555v3.txt"
i2, _ = find(f3, r"^Table 2:")
ie, te = find(f3, r"^Employment elasticity w\.r\.t MW \|", after=i2)
ig, tg = find(f3, r"^% MW changes \|", after=i2)
find(f3, r"^Note: The table reports the effects of a minimum wage increase", after=i2)
eta_str = cells(te)[1]; g_str = cells(tg)[1]
eta = num(eta_str); g = num(g_str)
v3 = 100 * eta * g
H["E3"] = {
    "value": round(v3, 2),
    "unit": "percent change in employment implied for the observed 2023 minimum-wage increase",
    "transformation": f"100*eta*g with eta = {eta_str} and g = {g_str} (proportional MW change printed in Table 2) -> 100*(-0.388)*0.047 = {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 {ie}); row '% MW changes' | column 'Estimate (SE)' (line {ig})",
}

# ---- E4: Armendariz Table 7 fatal, 5 km: absolute change needs a level
i7, _ = find(f1, r"^Table 7:")
ih7, th7 = find(f1, r"^Outcome \| 1 km \| 3 km \| 5 km", after=i7)
if7, tf7 = find(f1, r"^Drug-related \(fatal\) \|", after=i7)
b4_str = cells(tf7)[cells(th7).index("5 km")]
b4 = num(b4_str)
rel4 = 100 * (math.exp(b4) - 1)
print(f"E4 intermediate (relative change only, not the target): 100*(exp({b4})-1) = {rel4:.2f}%")
# Check whether any baseline/mean level of the fatal rate is printed in the extract
lvl = [t for t in lines(f1) if re.search(r"(?i)(mean|baseline|pre-treatment|average).*(rate|per 100,000)", t) and re.search(r"\d+\.\d+\s*(per|hospitali)", t)]
print(f"E4: lines printing a level of the fatal hospitalisation rate: {len(lvl)}")
H["E4"] = {
    "value": None,
    "unit": "hospitalisations per 100,000 inhabitants per year",
    "transformation": f"Would require R0*(exp(b)-1) with b = {b4_str.rstrip('*')} (relative change {rel4:.2f}%); R0 is not printed in the extract",
    "harmonisable": False,
    "missing_quantity": "the counterfactual (or pre-opening baseline) level of the fatal drug-related hospitalisation rate per 100,000 inhabitants per year for the treated regions within 5 km (e.g. a pre-treatment mean of the outcome in levels); the extract prints only log-point ATTs and no outcome level",
    "locator": f"arxiv-2609.27732v1 | Table 7 | row 'Drug-related (fatal)' | column '5 km' (extract line {if7})",
}

# ---- E5: Tossou Table 3, Treatment row, column (4): 0-1 -> pp
f5 = "arxiv-2510.03658v1.txt"
i3, _ = find(f5, r"^Table 3:")
ih5, th5 = find(f5, r"^\| \(1\) \| \(2\) \| \(3\) \| \(4\)", after=i3)
it5, tt5 = find(f5, r"^Treatment \|", after=i3)
im5, tm5 = find(f5, r"^Mean Dep\. Var\. \|", after=i3)
c4 = cells(th5).index("(4)")
b5_str = cells(tt5)[c4]
b5 = num(b5_str)
v5 = 100 * b5
print(f"E5 cross-check: relative change = b/mean = {b5}/{num(cells(tm5)[c4])} = {b5/num(cells(tm5)[c4]):.4f}")
H["E5"] = {
    "value": round(v5, 2),
    "unit": "percentage points (change in the probability of infant death or stillbirth)",
    "transformation": f"100*beta with beta = {b5_str} (0-1 scale) -> {v5:.2f} pp (exact rescaling)",
    "harmonisable": True,
    "missing_quantity": None,
    "locator": f"arxiv-2510.03658v1 | Table 3 | row 'Treatment' | column '(4)' (extract line {it5})",
}

# ---- data issues (detected by inspection; key checks recomputed here)
iA19, tA19 = find(f1, r"^\| \(0\.218\)")
find(f1, r"never-treated comparison regions \(640 regions\)")
find(f1, r"fatal outcomes, which fall by 35%")
find(f1, r"by approximately 30% at 5 km")
find(f3, r"about 4% minimum wage increase")
find(f2, r"increased reemployment by 1\.5 percentage points")
find(f5, r"by 0\.0855 \(a 18\.79 percentage change\)")
issues = [
    {"issue": "Armendariz: the same CS-DR 5 km ATT (-0.431) has standard error (0.210) in Table 6 but (0.218) in Table A19.",
     "handling": "Point estimate identical, so E1 is unaffected; no uncertainty propagated. Recorded only."},
    {"issue": "Armendariz: Tables 6/7 and the Figure 6 note state 608 never-treated control regions; the Table A19 note states 640.",
     "handling": "Does not enter any transformation; recorded only."},
    {"issue": f"Armendariz: abstract says fatal hospitalisations 'fall by 35%', but the Table 7 5 km estimate -0.354 implies 100*(exp(-0.354)-1) = {rel4:.1f}%, matching Section 6.1 ('approximately 30% at 5 km'). The abstract figure appears to repeat the overall effect.",
     "handling": "Used the Table 7 coefficient as instructed; E4 target (absolute change) is in any case not reachable."},
    {"issue": "Armendariz: Figure 7 note mentions 'four diagnostic sub-categories' and the Table 7 note lists ICD codes for four, while five outcomes (incl. fatal) are reported; the definition of the fatal outcome is not given in the extract.",
     "handling": "Recorded; does not affect the harmonised values."},
    {"issue": "Kanayama: Table 2 note describes the 2023 increase as 'about 4%' while the table row '% MW changes' prints 0.047.",
     "handling": "Used the printed table value 0.047 as the task requires (target definition points to the table). With 0.04 the implied change would be about -1.55%; the 0.047-based -1.82% is consistent with missing+excess jobs (-0.03+0.012 = -0.018) and with the abstract's ~2% decline."},
    {"issue": "RESEA: abstract rounds the reemployment ATE to 1.5 pp; Table 2 prints 1.47 pp. The 'Control' column is the raw aggregate control mean, not a block-adjusted counterfactual (the paper notes raw means are confounded by randomisation week).",
     "handling": "Used 1.47 pp and the printed control mean 75.0% as the target definition specifies (percent of control-group mean)."},
    {"issue": "Tossou: abstract/Section 4.1 headline (0.0855, 18.79 percentage change, mean 0.455) refers to a 'table 2' that is not in the extract and calls a relative change a 'percentage change'; Table 3 '% Change' is a relative change (b/mean, -0.0962/0.261 = -0.3686 vs printed -0.3685, rounding), not percentage points. Footnote markers are garbled by rendering ('^9^9 9').",
     "handling": "E5 uses Table 3 column (4) only; the target is percentage points (100*beta), distinct from the printed relative change."},
    {"issue": "RESEA extract: rendering defect 'p ¡ 0.001' (for p < 0.001) in Section 3 text.",
     "handling": "Irrelevant to the harmonised values; recorded."},
]

n_rows = {}
for fn in sorted(os.listdir(SRC)):
    with open(os.path.join(SRC, fn), encoding="utf-8") as f:
        n_rows["sources/" + fn] = f.read().count("\n")
with open(os.path.join(D, "README.md"), encoding="utf-8") as f:
    n_rows["README.md"] = f.read().count("\n")
with open(os.path.join(D, "effects.json"), encoding="utf-8") as f:
    eff = json.load(f)
n_rows["effects.json"] = len(eff)
assert [e["id"] for e in eff] == ["E1", "E2", "E3", "E4", "E5"]

res = {
    "task_id": "F4-conf-03",
    "harmonised": H,
    "n_harmonised": sum(v["value"] is not None for v in H.values()),
    "n_not_harmonisable": sum(not v["harmonisable"] for v in H.values()),
    "feasibility": "complete",
    "feasibility_reason": None,
    "data_issues": issues,
    "n_input_rows": n_rows,
    "n_analysed": len(H),
    "random_seed": None,
}
with open(os.path.join(ROOT, "out", "results.json"), "w", encoding="utf-8") as f:
    json.dump(res, f, indent=2, ensure_ascii=False)
print(json.dumps({k: H[k]["value"] for k in H}))