#!/usr/bin/env python3
"""reproduce_pool.py — independent re-walk of the pooled HR cited on the
musesolvescancer.com 'missing comparator' research branch (T-DM1 vs
trastuzumab+pertuzumab in residual HER2+ disease after neoadjuvant therapy).

Branch rationale claims: pooling the two observational cohorts that report an
estimate gives HR 0.483 (0.253-0.920), I^2 = 58.6%.
  cohort A: HR 0.62 (95% CI 0.30-1.26), n=184
  cohort B: HR 0.17 (95% CI 0.04-0.75), n=212
(third cohort, PMID 42434737, n=134 matched, no estimate published -> excluded)

Method (standard, public): log-scale inverse-variance pooling.
  SE_i = (ln upper_i - ln lower_i) / (2*1.96)
  fixed-effect:   w_i = 1/SE_i^2
  DerSimonian-Laird random-effects: tau^2 = max(0,(Q-(k-1))/C),
                  w*_i = 1/(SE_i^2 + tau^2)
Self-test asserts the claimed numbers reproduce under fixed-effect and reports
the model-dependence under random-effects. Zero dependencies; python3 stdlib.
"""
import math, json, sys

STUDIES = [
    {"id": "cohortA", "hr": 0.62, "lo": 0.30, "hi": 1.26, "n": 184},
    {"id": "cohortB", "hr": 0.17, "lo": 0.04, "hi": 0.75, "n": 212},
]
Z = 1.96
CLAIM = {"hr": 0.483, "lo": 0.253, "hi": 0.920, "i2": 58.6}


def se_from_ci(lo, hi):
    return (math.log(hi) - math.log(lo)) / (2 * Z)


def pool(studies, tau2=0.0):
    ys = [math.log(s["hr"]) for s in studies]
    vs = [se_from_ci(s["lo"], s["hi"]) ** 2 for s in studies]
    ws = [1.0 / (v + tau2) for v in vs]
    mu = sum(w * y for w, y in zip(ws, ys)) / sum(ws)
    se = math.sqrt(1.0 / sum(ws))
    return mu, se, ws, ys, vs


def main():
    # fixed effect
    mu, se, ws, ys, vs = pool(STUDIES)
    fe = {"hr": math.exp(mu), "lo": math.exp(mu - Z * se), "hi": math.exp(mu + Z * se)}
    # Cochran's Q and I^2 (k=2, df=1)
    Q = sum(w * (y - mu) ** 2 for w, y in zip(ws, ys))
    df = len(STUDIES) - 1
    i2 = max(0.0, (Q - df) / Q) * 100 if Q > 0 else 0.0
    # DerSimonian-Laird tau^2
    sw, sw2 = sum(ws), sum(w * w for w in ws)
    C = sw - sw2 / sw
    tau2 = max(0.0, (Q - df) / C)
    muR, seR, *_ = pool(STUDIES, tau2)
    re = {"hr": math.exp(muR), "lo": math.exp(muR - Z * seR), "hi": math.exp(muR + Z * seR)}

    out = {
        "inputs": STUDIES,
        "fixed_effect": {"hr": round(fe["hr"], 4), "lo": round(fe["lo"], 4),
                         "hi": round(fe["hi"], 4), "Q": round(Q, 4), "I2_pct": round(i2, 2)},
        "random_effect_DL": {"tau2": round(tau2, 4), "hr": round(re["hr"], 4),
                             "lo": round(re["lo"], 4), "hi": round(re["hi"], 4)},
        "claim_checked": CLAIM,
        "verdict": None,
    }
    fe_ok = (abs(fe["hr"] - CLAIM["hr"]) < 0.002
             and abs(fe["lo"] - CLAIM["lo"]) < 0.002
             and abs(fe["hi"] - CLAIM["hi"]) < 0.002
             and abs(i2 - CLAIM["i2"]) < 0.15)
    out["verdict"] = ("CONFIRM fixed-effect reproduces claim exactly; "
                      "NOTE random-effects CI crosses 1.0 (result is model-dependent)"
                      if fe_ok else "MISMATCH")
    print(json.dumps(out, indent=1))
    sys.exit(0 if fe_ok else 1)


if __name__ == "__main__":
    main()