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

ROOT = os.getcwd()
D = os.path.join(ROOT, "data")
SRC = os.path.join(D, "sources")


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


def find(lines, prefix, fname):
    """Return (line_no, line) of the unique line starting with prefix; print it."""
    hits = [(i + 1, l) for i, l in enumerate(lines) if l.startswith(prefix)]
    assert len(hits) == 1, (fname, prefix, hits)
    n, l = hits[0]
    print(f"[{fname}:{n}] {l}")
    return n, l


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


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


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

# n_input_rows: wc -l semantics (newline count) for text, top-level items for JSON
n_input_rows = {"effects.json": len(effects)}
for name in ["README.md"] + ["sources/" + v for v in files.values()]:
    with open(os.path.join(D, name), "rb") as f:
        n_input_rows[name] = f.read().count(b"\n")

H = {}
issues = []

# ---------------- E1: log-point ATT -> relative percent change ----------------
print("\nE1 source lines:")
_, hdr = find(L["dcr"], "| 1 km | 3 km | 5 km", files["dcr"])
n_att, att = find(L["dcr"], "ATT | -0.399", files["dcr"])
_, pct = find(L["dcr"], "% change | -32.9%", files["dcr"])
b1 = num(cells(att)[3])  # 5 km column
assert cells(hdr)[3] == "5 km"
v1 = 100 * (math.exp(b1) - 1)
printed_pct_5km = num(cells(pct)[3])
H["E1"] = {
    "value": round(v1, 2),
    "unit": "percent change relative to the counterfactual hospitalisation rate",
    "transformation": f"100*(exp(b)-1) with b = {cells(att)[3]} (log points) -> 100*(exp({b1})-1) = {v1:.4f} %; consistent with printed '% change' {cells(pct)[3]}",
    "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'",
}
assert abs(round(v1, 1) - printed_pct_5km) < 1e-9

# ---------------- E2: pp effect -> percent of control mean ----------------
print("\nE2 source lines:")
_, hdr2 = find(L["resea"], "Outcome | Treatment | Control |", files["resea"])
n_re, re_row = find(L["resea"], "Reemployment | 75.2%", files["resea"])
c = cells(re_row)
assert cells(hdr2)[2] == "Control" and cells(hdr2)[4] == "Average treatment effect estimate"
p_ctrl = num(c[2])  # 75.0 (%)
delta = num(c[4].split(";")[0])  # 1.47 pp
v2 = 100 * delta / p_ctrl
H["E2"] = {
    "value": round(v2, 2),
    "unit": "percent of the control-group mean reemployment rate",
    "transformation": f"100*Delta/p with Delta = {c[4].split(';')[0]} and p = {c[2]} (control mean) -> 100*{delta}/{p_ctrl} = {v2:.4f} %",
    "harmonisable": True,
    "missing_quantity": None,
    "locator": f"arxiv-2606.14621v1 | Table 2 | row 'Reemployment' | columns 'Average treatment effect estimate' (Delta) and 'Control' (p) (line {n_re})",
}

# ---------------- E3: elasticity x printed MW change ----------------
print("\nE3 source lines:")
n_el, el = find(L["mw"], "Employment elasticity w.r.t MW |", files["mw"])
n_g, gline = find(L["mw"], "% MW changes |", files["mw"])
_, note = find(L["mw"], "Note: The table reports the effects of a minimum wage increase", files["mw"])
eta = num(cells(el)[1])
g = num(cells(gline)[1])
v3 = 100 * eta * g
H["E3"] = {
    "value": round(v3, 2),
    "unit": "percent change in employment",
    "transformation": f"100*eta*g with eta = {cells(el)[1]} and g = {cells(gline)[1]} (proportional MW change, read as 4.7%; the table note says 'about 4%') -> 100*{eta}*{g} = {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}); row '% MW changes' | column 'Estimate (SE)' (line {n_g})",
}
issues.append({
    "issue": "arxiv-2505.04555v3 Table 2: the row '% MW changes' prints 0.047, a label suggesting percent units while the value is only consistent with a proportion (4.7%); the table note describes the change as 'about 4%'.",
    "handling": "Read 0.047 as a proportional change g = 0.047 (4.7%), the only reading compatible with the note ('about 4%'); the printed table value (not the rounded note value) is used."})

# ---------------- E4: log-point ATT on fatal hosp. -> absolute change per 100k ----------------
print("\nE4 source lines:")
_, hdr4 = find(L["dcr"], "Outcome | 1 km | 3 km | 5 km", files["dcr"])
n_f, fat = find(L["dcr"], "Drug-related (fatal) |", files["dcr"])
b4 = num(cells(fat)[3])
rel4 = 100 * (math.exp(b4) - 1)
# search the extract for any printed level of the outcome (mean/baseline rate); none expected
level_hits = [l for l in L["dcr"] if re.search(r"(mean|baseline rate|pre-treatment (mean|level)|average rate)", l, re.I)
              and re.search(r"fatal", l, re.I)]
print(f"Lines printing a level of the fatal-hospitalisation rate: {len(level_hits)}")
H["E4"] = {
    "value": None,
    "unit": "hospitalisations per 100,000 inhabitants per year",
    "transformation": f"Absolute change = R0*(exp(b)-1) with b = {cells(fat)[3]} (log points), i.e. relative change 100*(exp({b4})-1) = {rel4:.4f} %; R0 (counterfactual/baseline fatal rate per 100,000 in treated regions) is not printed in the extract, so the absolute change cannot be computed",
    "harmonisable": False,
    "missing_quantity": "baseline (counterfactual, e.g. pre-opening or control-group mean) level of the fatal drug-related hospitalisation rate per 100,000 inhabitants per year in the 5 km treated regions",
    "locator": f"arxiv-2609.27732v1 | Table 7 | row 'Drug-related (fatal)' | column '5 km' (line {n_f})",
}
assert not level_hits

# ---------------- E5: LPM coefficient -> percentage points ----------------
print("\nE5 source lines:")
_, hdr5 = find(L["cs"], "| (1) | (2) | (3) | (4)", files["cs"])
n_t, tr = find(L["cs"], "Treatment |", files["cs"])
_, mean5 = find(L["cs"], "Mean Dep. Var. |", files["cs"])
_, pc5 = find(L["cs"], "% Change |", files["cs"])
b5 = num(cells(tr)[4])
v5 = 100 * b5
H["E5"] = {
    "value": round(v5, 2),
    "unit": "percentage points",
    "transformation": f"100*beta with beta = {cells(tr)[4]} (change in probability, 0-1 scale) -> 100*({b5}) = {v5:.4f} pp",
    "harmonisable": True,
    "missing_quantity": None,
    "locator": f"arxiv-2510.03658v1 | Table 3 | row 'Treatment' | column '(4)' (line {n_t})",
}

# ---------------- consistency checks on the sources (recorded as data issues) ----------------
print("\nSource-consistency checks:")
_, abs_dcr = find(L["dcr"], "This paper estimates the causal effect of drug consumption room", files["dcr"])
issues.append({
    "issue": f"arxiv-2609.27732v1: the abstract states that fatal outcomes 'fall by 35%', but Table 7 (5 km) prints {cells(fat)[3]} log points, i.e. 100*(exp({b4})-1) = {rel4:.1f}%; Section 6.1 says 'approximately 30%', consistent with Table 7. The 35% matches the overall ATT (Table 6), not the fatal outcome.",
    "handling": "Table 7 is used as the authoritative printed estimate for E4; the abstract figure is not used."})
n_a19, a19 = find(L["dcr"], "ATT | -0.431** | -0.438**", files["dcr"])
a19_se = L["dcr"][n_a19 + 1]
t6_se = L["dcr"][n_att + 1]
print(f"[{files['dcr']}:{n_att+2}] {t6_se}")
print(f"[{files['dcr']}:{n_a19+2}] {a19_se}")
issues.append({
    "issue": f"arxiv-2609.27732v1: the 5 km CS-DR ATT -0.431 has SE {cells(t6_se)[3]} in Table 6 but {cells(a19_se)[1]} in Table A19; the Table A19 note mentions 640 never-treated comparison regions while the table body and Tables 6-7 print 608.",
    "handling": "Does not affect the point estimate used in E1 (identical -0.431 in both tables); Table 6 is used as the cited location. No uncertainty is propagated."})
_, a19pct = find(L["dcr"], "% change | -35.0% | -35.5%", files["dcr"])
mism = []
for b, p in zip(cells(a19)[1:], cells(a19pct)[1:]):
    bb, pp = num(b), num(p)
    implied = 100 * (math.exp(bb) - 1)
    if abs(round(implied, 1) - pp) > 1e-9:
        mism.append(f"{b} -> {implied:.2f}% vs printed {p}")
print("A19 % change mismatches:", mism)
if mism:
    issues.append({"issue": "arxiv-2609.27732v1 Table A19: printed '% change' differs in the last decimal from 100*(exp(ATT)-1) for: " + "; ".join(mism) + ".",
                   "handling": "Minor rounding inconsistency in a robustness table not used for any item; recorded only."})
ratio = b5 / num(cells(mean5)[4])
print(f"Tossou col (4): beta/mean = {ratio:.4f} vs printed % Change {cells(pc5)[4]}")
issues.append({
    "issue": f"arxiv-2510.03658v1: Table 3 column (4) prints '% Change' {cells(pc5)[4]} while {b5}/{num(cells(mean5)[4])} = {ratio:.4f} (last-digit difference, presumably from an unrounded mean). The abstract's '0.0855 (a 18.79 percentage change)' refers to a 'table 2' (mean 0.455) that is not in the extract.",
    "handling": "E5 uses the Table 3 column (4) coefficient as printed and converts only the scale (x100 -> pp); the abstract/Table 2 numbers and the '% Change' row are not used."})
issues.append({
    "issue": "Rendering defects: arxiv-2510.03658v1 footnote markers are rendered inline (e.g. 'percent^9^9 9') and footnote 11 is interleaved with main text; arxiv-2606.14621v1 renders 'p < 0.001' as 'p ¡ 0.001'.",
    "handling": "No numbers used in E1-E5 are affected; recorded only."})
issues.append({
    "issue": "arxiv-2606.14621v1 Table 2: the control-group mean (75.0%) is a raw aggregate mean, whereas the ATE (1.47 pp) adjusts for randomisation week (blocks); the abstract rounds the ATE to 1.5 pp.",
    "handling": "Following the target definition ('percent of the control-group mean'), the printed control mean 75.0% is used as p and the Table 2 ATE 1.47 pp (not the rounded 1.5) as Delta; the denominator mismatch is noted as a limitation."})

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)
results = {
    "task_id": "F4-conf-03",
    "harmonised": H,
    "n_harmonised": n_h,
    "n_not_harmonisable": n_nh,
    "feasibility": "complete",
    "feasibility_reason": None,
    "data_issues": issues,
    "n_input_rows": n_input_rows,
    "n_analysed": len(effects),
    "random_seed": None,
    # P2 baseline: the printed estimates on their reported definitions, untransformed
    "baseline_printed_values": {"E1": b1, "E2": delta, "E3": eta, "E4": b4, "E5": b5},
    "baseline_printed_units": {"E1": "log points", "E2": "percentage points", "E3": "elasticity (unitless)",
                               "E4": "log points", "E5": "probability change (0-1 scale)"},
}
assert [e["id"] for e in effects] == list(H.keys())
assert n_h + n_nh == len(effects)
with open(os.path.join(ROOT, "out", "results.json"), "w", encoding="utf-8") as f:
    json.dump(results, f, indent=2, ensure_ascii=False)
print("\nSummary:")
for k, v in H.items():
    print(f"  {k}: value={v['value']} [{v['unit']}] harmonisable={v['harmonisable']}")
print(f"n_harmonised={n_h} n_not_harmonisable={n_nh} feasibility=complete")