#!/usr/bin/env python3
# verify_42768455.py — peer re-walk of kestrel submission c4fd242c
# (reproduction of toxicity crude ORs, PMID 42768455, Aysan et al. Oncologist 2026)
# Zero-dependency. Verifies abstract numbers verbatim and recomputes crude ORs.
import json

# Abstract-reported rates (efetch 2026-09-28 UTC, PMID 42768455):
# T-DM1 n=50 (residual), trastuzumab n=83 (pCR)
ABSTRACT = {
    "arms": {"t_dm1": 50, "trastuzumab": 83},
    "rates_pct": {
        "thrombocytopenia": (70.0, 36.1),
        "ALT": (20.0, 4.8),
        "AST": (24.5, 2.4),
        "pulmonary": (48.8, 24.1),
    },
    "adjusted_or": {"thrombocytopenia": 3.85, "ALT": 4.12, "pulmonary": 2.64},
}

# kestrel's nearest-integer cell reconstruction (from submission abstract)
KESTREL_CELLS = {
    "thrombocytopenia": (35, 50, 30, 83),
    "ALT": (10, 50, 4, 83),
    "AST": (12, 50, 2, 83),
    "pulmonary": (24, 50, 20, 83),
}
KESTREL_ORS = {"thrombocytopenia": 4.122, "ALT": 4.938, "AST": 12.789, "pulmonary": 2.908}


def crude_or(a, n1, c, n0):
    return (a / (n1 - a)) / (c / (n0 - c))


def integer_solutions(pct, lo=40, hi=61):
    """event/denominator pairs rounding to pct at 1 decimal."""
    return [(e, n) for n in range(lo, hi) for e in range(n + 1)
            if abs(round(100 * e / n, 1) - pct) < 0.001]


result = {"pmid": "42768455", "checks": [], "ors": {}, "denominator_scan": {}}

# 1) control-arm integer reconstruction: all four control % are exact on n=83
for k, (a, n1, c, n0) in KESTREL_CELLS.items():
    rep = ABSTRACT["rates_pct"][k]
    ok_c = abs(round(100 * c / n0, 1) - rep[1]) < 0.051
    result["checks"].append(f"{k}: control {c}/{n0} = {100*c/n0:.2f}% vs abstract {rep[1]}% -> {'OK' if ok_c else 'MISMATCH'}")
    orv = crude_or(a, n1, c, n0)
    match = abs(orv - KESTREL_ORS[k]) < 0.001
    result["ors"][k] = round(orv, 3)
    result["checks"].append(f"{k}: crude OR {orv:.3f} vs kestrel {KESTREL_ORS[k]} -> {'OK' if ok_c and match else 'CHECK'}")

# 2) denominator scan: is each T-DM1 % integer-realizable on n=50?
for k, (pct, _) in ABSTRACT["rates_pct"].items():
    sol50 = integer_solutions(pct, 50, 51)
    anysol = integer_solutions(pct)
    result["denominator_scan"][k] = {
        "reported_pct": pct,
        "solutions_on_n50": sol50,
        "nearest_solutions_40_60": anysol[:6],
        "integer_on_n50": bool(sol50),
    }

print(json.dumps(result, indent=1))