"""F4-conf-03: harmonise five printed effect sizes. Reads only data/, writes out/results.json.
Run from workspace root: ~/wp29/venv/bin/python out/analysis.py
No randomness is involved (random_seed = null)."""
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(fname):
    with open(os.path.join(SRC, fname), encoding="utf-8") as f:
        return f.read().split("\n")


def find_line(fname, prefix):
    """Return (1-based line number, text) of the unique line starting with prefix; prints it."""
    L = lines(fname)
    hits = [(i + 1, l) for i, l in enumerate(L) if l.startswith(prefix)]
    assert len(hits) == 1, (fname, prefix, hits)
    n, t = hits[0]
    print(f"[{fname}:{n}] {t}")
    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.replace("−", "-"))
    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(D, "effects.json"), encoding="utf-8"))
assert [e["id"] for e in effects] == ["E1", "E2", "E3", "E4", "E5"]

DCR = "arxiv-2609.27732v1.txt"
RES = "arxiv-2606.14621v1.txt"
MW = "arxiv-2505.04555v3.txt"
BEN = "arxiv-2510.03658v1.txt"

H = {}
checks = {}

# ---------------- E1: log points -> relative percent change ----------------
n_hdr, hdr = find_line(DCR, "| 1 km | 3 km | 5 km")
n_att, att = find_line(DCR, "ATT | -0.399")
n_pc, pc = find_line(DCR, "% change | -32.9%")
c = cells(att)
assert cells(hdr)[3] == "5 km" and c[3] == "-0.431**"
b1 = num(c[3])
e1 = 100 * (math.exp(b1) - 1)
printed_pc_5km = num(cells(pc)[3])
checks["E1_vs_printed_pct_change"] = {"computed": round(e1, 4), "printed": printed_pc_5km,
                                      "agree_to_printed_precision": round(e1, 1) == printed_pc_5km}
H["E1"] = {
    "value": round(e1, 2),
    "unit": "percent change relative to the counterfactual hospitalisation rate",
    "transformation": f"100*(exp(b)-1) with b = {b1:.3f} log points: 100*(exp(-0.431)-1) = {e1:.2f}% (printed '% change' row: -35.0%)",
    "harmonisable": True,
    "missing_quantity": None,
    "locator": f"arxiv-2609.27732v1 | Table 6 | row 'ATT' | column '5 km' (value -0.431**, extract line {n_att}); cross-check row '% change', column '5 km' (-35.0%, line {n_pc})",
}

# ---------------- E2: pp -> percent of control mean ----------------
n_r, r = find_line(RES, "Reemployment | 75.2%")
c = cells(r)
assert c[2] == "75.0%" and c[4].startswith("1.47 pp")
p_ctrl = num(c[2])          # 75.0 (%)
ate_pp = num(c[4])          # 1.47 pp
e2 = 100 * ate_pp / p_ctrl
H["E2"] = {
    "value": round(e2, 2),
    "unit": "percent change relative to the control-group mean reemployment rate",
    "transformation": f"100*Delta/p with Delta = {ate_pp:.2f} pp and p = {p_ctrl:.1f}% (control mean): 100*1.47/75.0 = {e2:.4f}%",
    "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_r})",
}

# ---------------- E3: elasticity * observed MW change ----------------
n_el, el = find_line(MW, "Employment elasticity w.r.t MW")
n_g, g = find_line(MW, "% MW changes")
n_note, note = find_line(MW, "Note: The table reports the effects of a minimum wage increase")
assert cells(el)[1] == "-0.388" and cells(g)[1] == "0.047"
eta = num(cells(el)[1])     # -0.388
gg = num(cells(g)[1])       # 0.047 (proportion)
e3 = 100 * eta * gg
H["E3"] = {
    "value": round(e3, 2),
    "unit": "percent change in employment",
    "transformation": f"100*eta*g with eta = {eta:.3f} and g = {gg:.3f} (proportional MW change, Table 2 row '% MW changes'): 100*(-0.388)*0.047 = {e3:.4f}%",
    "harmonisable": True,
    "missing_quantity": None,
    "locator": f"arxiv-2505.04555v3 | Table 2 | row 'Employment elasticity w.r.t MW' | column 'Estimate (SE)' (-0.388, line {n_el}); row '% MW changes' | column 'Estimate (SE)' (0.047, line {n_g})",
}
checks["E3_alt_with_note_4pct"] = round(100 * eta * 0.04, 4)  # sensitivity only, not the reported value

# ---------------- E4: log points -> absolute change per 100,000 ----------------
n_f, f = find_line(DCR, "Drug-related (fatal) |")
assert cells(f)[3] == "-0.354**"
b4 = num(cells(f)[3])       # -0.354
rel4 = 100 * (math.exp(b4) - 1)
checks["E4_relative_pct_only"] = round(rel4, 4)
H["E4"] = {
    "value": None,
    "unit": "hospitalisations per 100,000 inhabitants per year",
    "transformation": (f"absolute change = R0*(exp(b)-1) with b = {b4:.3f}; exp(-0.354)-1 = {math.exp(b4)-1:.4f} "
                       f"(relative change {rel4:.2f}%), but R0, the counterfactual (or pre-treatment) level of the fatal "
                       "drug-related hospitalisation rate per 100,000 in treated 5 km regions, is not printed in the extract"),
    "harmonisable": False,
    "missing_quantity": ("baseline/counterfactual level R0 of the fatal drug-related hospitalisation rate per 100,000 inhabitants "
                         "per year for the treated 5 km-buffer regions (e.g. pre-opening or never-treated mean rate); "
                         "not printed in Section 6, Table 7, Section 6.1 or Table A19 of the extract"),
    "locator": f"arxiv-2609.27732v1 | Table 7 | row 'Drug-related (fatal)' | column '5 km' (-0.354**, extract line {n_f})",
}

# ---------------- E5: 0-1 probability change -> percentage points ----------------
n_hd, hd = find_line(BEN, "| (1) | (2) | (3) | (4)")
n_t, t = find_line(BEN, "Treatment | -0.0911")
n_m, m = find_line(BEN, "Mean Dep. Var.")
n_c, cc = find_line(BEN, "% Change |")
assert cells(hd)[4] == "(4)" and cells(t)[4] == "-0.0962"
b5 = num(cells(t)[4])       # -0.0962
e5 = 100 * b5
mean5 = num(cells(m)[4]); pct5 = num(cells(cc)[4])
checks["E5_pct_change_consistency"] = {"b_over_mean": round(b5 / mean5, 4), "printed": pct5}
H["E5"] = {
    "value": round(e5, 2),
    "unit": "percentage points",
    "transformation": f"100*b with b = {b5:.4f} (0-1 scale): 100*(-0.0962) = {e5:.2f} pp",
    "harmonisable": True,
    "missing_quantity": None,
    "locator": f"arxiv-2510.03658v1 | Table 3 | row 'Treatment' | column '(4)' (-0.0962, extract line {n_t})",
}

# ---------------- internal-consistency diagnostics of sources ----------------
n_a19, a19 = find_line(DCR, "ATT | -0.431** | -0.438**")
n_a19p, a19p = find_line(DCR, "% change | -35.0% | -35.5%")
find_line(DCR, "| (0.218) | (0.208)")
find_line(DCR, "| (0.412) | (0.279) | (0.210)")
find_line(DCR, "never-treated comparison regions (640 regions)")
pc_checks = []
for bt, pt in list(zip(cells(att)[1:], cells(pc)[1:])) + list(zip(cells(a19)[1:], cells(a19p)[1:])):
    bv, pv = num(bt), num(pt)
    cv = 100 * (math.exp(bv) - 1)
    pc_checks.append({"b": bv, "printed_pct": pv, "computed_pct": round(cv, 2), "match_1dp": round(cv, 1) == pv})
checks["DCR_pct_change_rows"] = pc_checks
narc_b = num(cells(find_line(DCR, "Narcotics |")[1])[3])
mental_b = num(cells(find_line(DCR, "Mental/Behavioral |")[1])[3])
checks["DCR_text_claims_5km_pct"] = {"fatal": round(rel4, 2), "narcotics": round(100 * (math.exp(narc_b) - 1), 2),
                                     "mental": round(100 * (math.exp(mental_b) - 1), 2)}
find_line(DCR, "This paper estimates the causal effect")
find_line(DCR, "The sub-outcome results point to two channels")

mism = [x for x in pc_checks if not x["match_1dp"]]
data_issues = [
    {"issue": "arxiv-2609.27732v1: the abstract states fatal drug-related hospitalisations 'fall by 35%%', but Section 6.1 says 'approximately 30%% at 5 km' and the Table 7 5 km coefficient -0.354 implies 100*(exp(-0.354)-1) = %.1f%%." % rel4,
     "handling": "Internal inconsistency of the source. E4 uses only the Table 7 coefficient as specified; the 35% figure in the abstract is not used."},
    {"issue": "arxiv-2609.27732v1: the standard error of the same 5 km CS-DR ATT (-0.431) is printed as (0.210) in Table 6 but (0.218) in Table A19.",
     "handling": "Does not affect point estimates (no uncertainty propagation required); recorded only."},
    {"issue": "arxiv-2609.27732v1: the Table A19 note says 640 never-treated comparison regions while the tables (and the Table 6 note) report 608 control regions.",
     "handling": "Recorded; does not affect the harmonised values."},
]
for x in mism:
    data_issues.append({"issue": "arxiv-2609.27732v1: printed '%% change' %.1f%% for b = %.3f does not equal 100*(exp(b)-1) = %.2f%% at one decimal." % (x["printed_pct"], x["b"], x["computed_pct"]),
                        "handling": "Minor discrepancy (at most 0.1 percentage point; plausibly from the printed coefficient being rounded to 3 decimals) in a cell not used for any harmonised item (the E1 cell matches); recorded only."})
data_issues += [
    {"issue": "arxiv-2609.27732v1: no level (mean) of any hospitalisation rate is printed in the extract, so log-point effects cannot be converted to absolute changes per 100,000.",
     "handling": "E4 declared not harmonisable (value null) with the missing quantity named; no external value substituted."},
    {"issue": "arxiv-2505.04555v3: Table 2 prints '% MW changes' = 0.047 while the table note says 'about 4% minimum wage increase'. Despite the '%' label, 0.047 is read as a proportion (4.7%), consistent with the note's order of magnitude.",
     "handling": "Used the table value g = 0.047 as the target definition requires ('as printed in the same table'); with g = 0.04 the implied change would be %.2f%% (sensitivity only, not reported as the value)." % checks["E3_alt_with_note_4pct"]},
    {"issue": "arxiv-2510.03658v1: the abstract and the first paragraph of Section 4.1 report an effect of 0.0855 (18.79%%, mean 0.455) attributed to 'table 2', which is not in the extract; Table 3 column (4) prints -0.0962 with Mean Dep. Var. 0.261 and %% Change -0.3685 (computed b/mean = %.4f). Footnote markers are rendered as duplicated text (e.g. 'percent^9^9 9')." % (b5 / mean5),
     "handling": "E5 uses Table 3 column (4) as specified by effects.json; the 0.0855 figure is not used. The rendering defect does not affect any number used."},
    {"issue": "arxiv-2510.03658v1: the source calls b/mean a 'percentage change' (a relative change), which differs from the percentage-point target definition.",
     "handling": "E5 reports the absolute change in percentage points (100*b); the relative figure is not used as the harmonised value."},
    {"issue": "arxiv-2606.14621v1: the Table 2 control mean (75.0%) is a raw aggregate group mean, whereas the ATE (1.47 pp) adjusts for randomisation week (the paper notes that raw group comparisons are biased under blocking). Also a rendering defect 'p ¡ 0.001' in the text; the abstract rounds the effect to 1.5 pp.",
     "handling": "Used the printed control mean as p, as the target definition specifies (percent of the control-group mean); flagged as approximate. Used 1.47 pp from Table 2, not the rounded 1.5."},
]

n_input_rows = {"effects.json": len(effects), "README.md": wc_l(os.path.join(D, "README.md"))}
for fn in sorted(os.listdir(SRC)):
    n_input_rows["sources/" + 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)
assert n_h + n_nh == 5

results = {
    "task_id": "F4-conf-03",
    "harmonised": H,
    "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,
    "baseline_printed_estimates": {"E1": b1, "E2": ate_pp, "E3": eta, "E4": b4, "E5": b5},
    "consistency_checks": checks,
}
with open(os.path.join(ROOT, "out", "results.json"), "w", encoding="utf-8") as fo:
    json.dump(results, fo, indent=2, ensure_ascii=False)
print("\nHarmonised values:")
for k, v in H.items():
    print(f"  {k}: value={v['value']} unit={v['unit']} harmonisable={v['harmonisable']}")
print(f"n_harmonised={n_h} n_not_harmonisable={n_nh} feasibility=complete")