"""
test_calibration_vs_cohorts.py
==============================

Calibration test for the Nivel-1 ypN0 logistic predictor against six
published cohorts. For each cohort we:

1.  Encode the reported case-mix (subtype proportion, modal cN, dominant
    regimen, median age, IO/HER2 agent flags).
2.  Call ``LogisticPredictor.predict(...)`` for every subtype slice.
3.  Weight predictions by case-mix proportion -> aggregate predicted ypN0.
4.  Compare against the published observed ypN0 rate.
5.  Compute calibration metrics (mean absolute error, slope, intercept)
    and write a Markdown report at ``tests/calibration_report.md``.

Case-mix proportions are extracted from the cohort source files in
``02_ensayos_clinicos/`` and ``05_tecnicas_quirurgicas_TAD_SLNB/``.
Where the publication does NOT report exact proportions (most cases),
we use plausible defaults from contemporaneous breast-cancer
epidemiology and the trial's eligibility criteria, and flag the
assumption in the ``case_mix_notes`` field.

NOTE on Montagna 2024: the cohort was *selected* by ypN0 (all 1144
patients are ypN0 by inclusion). It is therefore NOT a calibration
target for the predictor's marginal ypN0 rate. We compute predictions
anyway and flag it as ``EXCLUDED FROM CALIBRATION`` so future readers
do not over-interpret the comparison.
"""

from __future__ import annotations

import math
import os
import sys
from pathlib import Path
from typing import Optional

# ---------------------------------------------------------------------------
# Predictor import with fallback
# ---------------------------------------------------------------------------
THIS_DIR = Path(__file__).resolve().parent
ROOT = THIS_DIR.parent
if str(ROOT) not in sys.path:
    sys.path.insert(0, str(ROOT))

_FALLBACK = False
try:
    from nivel1_calculadora.predictor import LogisticPredictor  # type: ignore
except Exception as exc:                                          # noqa: BLE001
    _FALLBACK = True
    print(f"[warn] could not import LogisticPredictor ({exc!r}); using fallback.")

    class LogisticPredictor:                                     # type: ignore[no-redef]
        """Minimal fallback replicating the documented coefficients.

        Used only if Agent B's module is unavailable. Coefficients are
        copied from predictor_data_dictionary.md / predictor_effect_sizes.csv
        so the test script is self-contained.
        """

        _SUBTYPE_OR = {
            "HR+/HER2-": 1.0,
            "TNBC": 6.91,
            "HR+/HER2+": 14.82,
            "HR-/HER2+": 40.0,
        }
        _INTERCEPT = math.log(0.10 / 0.90)
        _AGE_OR = 0.99
        _CN23_OR = 0.65
        _PERTU_OR = 1.5
        _CARBO_OR = 1.4
        _PEMBRO_OR = 1.42
        _KI67_HI_OR = 3.0
        _GEN_HI_OR = 2.0
        _GEN_LO_OR = 0.5

        def predict(
            self,
            subtype: str,
            ki67_pct: Optional[float],
            age: float,
            cN_stage: str,
            regimen: str,
            has_pertuzumab: bool = False,
            has_carboplatin: bool = False,
            has_pembrolizumab: bool = False,
            genomic_risk: str = "unknown",
        ) -> dict:
            logit = self._INTERCEPT + math.log(self._SUBTYPE_OR[subtype])
            if subtype == "HR+/HER2-" and ki67_pct is not None and ki67_pct >= 20:
                logit += math.log(self._KI67_HI_OR)
            logit += (age - 50.0) * math.log(self._AGE_OR)
            if cN_stage in ("cN2", "cN3"):
                logit += math.log(self._CN23_OR)
            if has_pertuzumab and subtype in ("HR+/HER2+", "HR-/HER2+"):
                logit += math.log(self._PERTU_OR)
            if has_carboplatin and subtype == "TNBC":
                logit += math.log(self._CARBO_OR)
            if has_pembrolizumab and subtype == "TNBC":
                logit += math.log(self._PEMBRO_OR)
            if subtype == "HR+/HER2-":
                if genomic_risk == "high":
                    logit += math.log(self._GEN_HI_OR)
                elif genomic_risk == "low":
                    logit += math.log(self._GEN_LO_OR)
            prob = 1.0 / (1.0 + math.exp(-logit))
            band = "low" if prob < 0.2 else ("medium" if prob < 0.5 else "high")
            return {"probability": round(prob, 4), "risk_band": band,
                    "confidence": "medium"}


# ---------------------------------------------------------------------------
# Cohorts
# ---------------------------------------------------------------------------
COHORTS = [
    {
        "name": "ACOSOG Z1071 (Boughey 2013, JAMA)",
        "pmid": "24101169",
        "n": 525,
        "observed_pCR": 0.410,
        "case_mix": {
            "HR-/HER2+": 0.18,
            "HR+/HER2+": 0.16,
            "TNBC": 0.30,
            "HR+/HER2-": 0.36,
        },
        "median_age": 50,
        "cN_modal": "cN1",
        "regimen": "anthra_taxane",
        "has_pertuzumab": False,
        "has_carboplatin": False,
        "has_pembrolizumab": False,
        "ki67_assumption": {"HR+/HER2-": 25},
        "case_mix_notes": (
            "Pre-pertuzumab / pre-pembro era (2009-2011). 74.6% anthra+taxane."
            " Subtype proportions extrapolated from contemporaneous US distribution;"
            " HER2+ slightly enriched because trial required cN1. cN1 = 95%."
        ),
        "include_in_calibration": True,
    },
    {
        "name": "RISAS (Simons 2022, JAMA Surg)",
        "pmid": "36069889",
        "n": 212,
        "observed_pCR": 0.304,
        "case_mix": {
            "HR-/HER2+": 0.15,
            "HR+/HER2+": 0.15,
            "TNBC": 0.25,
            "HR+/HER2-": 0.45,
        },
        "median_age": 52,
        "cN_modal": "cN1",
        "regimen": "anthra_taxane",
        "has_pertuzumab": True,
        "has_carboplatin": False,
        "has_pembrolizumab": False,
        "ki67_assumption": {"HR+/HER2-": 25},
        "case_mix_notes": (
            "Dutch protocols 2017-2019. Most HER2+ received dual blockade by then."
            " Subtype mix inferred from Dutch national registry of cN+ cohorts."
            " cN1 dominant (~88%)."
        ),
        "include_in_calibration": True,
    },
    {
        "name": "RISAS sub-analysis (van Amstel 2024, BJS)",
        "pmid": "39302345",
        "n": 185,
        "observed_pCR": 0.297,
        "case_mix": {
            "HR-/HER2+": 0.14,
            "HR+/HER2+": 0.16,
            "TNBC": 0.25,
            "HR+/HER2-": 0.45,
        },
        "median_age": 52,
        "cN_modal": "cN1",
        "regimen": "anthra_taxane",
        "has_pertuzumab": True,
        "has_carboplatin": False,
        "has_pembrolizumab": False,
        "ki67_assumption": {"HR+/HER2-": 25},
        "case_mix_notes": (
            "Subset of RISAS with baseline FDG-PET/CT (n=185). Same case-mix as"
            " parent RISAS cohort assumed."
        ),
        "include_in_calibration": True,
    },
    {
        "name": "Montagna 2024 (JAMA Oncol)",
        "pmid": "38662396",
        "n": 1144,
        "observed_pCR": 1.000,       # selection bias: cohort = ypN0 only
        "case_mix": {
            "HR-/HER2+": 0.30,       # 54% HER2+ overall in the paper
            "HR+/HER2+": 0.24,
            "TNBC": 0.20,
            "HR+/HER2-": 0.26,
        },
        "median_age": 50,
        "cN_modal": "cN1",
        "regimen": "anthra_taxane",
        "has_pertuzumab": True,
        "has_carboplatin": False,
        "has_pembrolizumab": False,
        "ki67_assumption": {"HR+/HER2-": 25},
        "case_mix_notes": (
            "All 1144 patients are ypN0 by inclusion (selection bias)."
            " 619 (54%) HER2+ -> HER2+ split 30/24 by HR status."
            " EXCLUDED FROM CALIBRATION because observed rate = 1.0 by design."
        ),
        "include_in_calibration": False,
    },
    {
        "name": "KEYNOTE-522 surgical (Kuemmel 2026, Ann Surg Oncol)",
        "pmid": "41739398",
        "n": 784,
        "observed_pCR": 0.767,
        "case_mix": {"TNBC": 1.0},
        "median_age": 49,
        "cN_modal": "cN1",
        "regimen": "AC_T",
        "has_pertuzumab": False,
        "has_carboplatin": True,
        "has_pembrolizumab": True,
        "ki67_assumption": {},
        "case_mix_notes": (
            "Pembrolizumab arm. All TNBC. T1cN1-2 or T2-4 N0-2 -- modal cN1."
            " 5-7% had cN2 (assumed cN1)."
        ),
        "include_in_calibration": True,
    },
    {
        "name": "NeoSTOP arm A (Sharma 2021, Clin Cancer Res)",
        "pmid": "33208340",
        "n": 48,
        "observed_pCR": 0.540,       # breast+axilla pCR (proxy for ypN0)
        "case_mix": {"TNBC": 1.0},
        "median_age": 50,
        "cN_modal": "cN1",
        "regimen": "AC_T",
        "has_pertuzumab": False,
        "has_carboplatin": True,
        "has_pembrolizumab": False,
        "ki67_assumption": {},
        "case_mix_notes": (
            "Carboplatin + paclitaxel -> dose-dense AC. Stage I-III TNBC."
            " Includes cN0 patients (mixed cohort); observed pCR = breast+axilla pCR."
            " Median age inferred (typical TNBC trial). Caveat: observed rate"
            " is breast+axilla composite, not pure axillary pCR -> floor."
        ),
        "include_in_calibration": True,
    },
]


# ---------------------------------------------------------------------------
# Prediction logic
# ---------------------------------------------------------------------------
def _predict_for_cohort(predictor: LogisticPredictor, cohort: dict) -> dict:
    """Run weighted prediction for one cohort."""
    contributions = []
    aggregate = 0.0
    for subtype, weight in cohort["case_mix"].items():
        ki67 = cohort.get("ki67_assumption", {}).get(subtype)
        result = predictor.predict(
            subtype=subtype,
            ki67_pct=ki67,
            age=cohort["median_age"],
            cN_stage=cohort["cN_modal"],
            regimen=cohort["regimen"],
            has_pertuzumab=cohort["has_pertuzumab"],
            has_carboplatin=cohort["has_carboplatin"],
            has_pembrolizumab=cohort["has_pembrolizumab"],
            genomic_risk="unknown",
        )
        p = result["probability"]
        aggregate += weight * p
        contributions.append({
            "subtype": subtype,
            "weight": weight,
            "predicted_subtype": p,
            "weighted": weight * p,
        })
    return {"predicted": aggregate, "by_subtype": contributions}


def _verdict(diff: float) -> str:
    """Reliability flag for predicted-vs-observed difference."""
    d = abs(diff)
    if d <= 0.05:
        return "OK (excellent)"
    if d <= 0.10:
        return "OK"
    if d <= 0.20:
        return "WARN"
    return "FAIL"


# ---------------------------------------------------------------------------
# Calibration metrics
# ---------------------------------------------------------------------------
def _calibration_slope(preds: list[float], obs: list[float]) -> tuple[float, float]:
    """Return (slope, intercept) of ordinary-least-squares regression obs ~ preds."""
    n = len(preds)
    if n < 2:
        return float("nan"), float("nan")
    mx = sum(preds) / n
    my = sum(obs) / n
    num = sum((p - mx) * (o - my) for p, o in zip(preds, obs))
    den = sum((p - mx) ** 2 for p in preds)
    if den == 0:
        return float("nan"), float("nan")
    slope = num / den
    intercept = my - slope * mx
    return slope, intercept


# ---------------------------------------------------------------------------
# Driver
# ---------------------------------------------------------------------------
def run() -> dict:
    predictor = LogisticPredictor()
    rows = []
    for cohort in COHORTS:
        result = _predict_for_cohort(predictor, cohort)
        diff = result["predicted"] - cohort["observed_pCR"]
        rows.append({
            "cohort": cohort,
            "predicted": result["predicted"],
            "by_subtype": result["by_subtype"],
            "diff": diff,
            "verdict": _verdict(diff) if cohort["include_in_calibration"]
                       else "N/A (selection-biased cohort)",
        })

    calibration_rows = [r for r in rows if r["cohort"]["include_in_calibration"]]
    preds = [r["predicted"] for r in calibration_rows]
    obs = [r["cohort"]["observed_pCR"] for r in calibration_rows]
    slope, intercept = _calibration_slope(preds, obs)
    mae = sum(abs(r["diff"]) for r in calibration_rows) / max(len(calibration_rows), 1)

    summary = {
        "rows": rows,
        "n_calibration": len(calibration_rows),
        "slope": slope,
        "intercept": intercept,
        "mae": mae,
        "fallback_used": _FALLBACK,
    }
    return summary


def write_report(summary: dict, path: Path) -> None:
    lines: list[str] = []
    lines.append("# Calibration report - Nivel 1 ypN0 predictor")
    lines.append("")
    if summary["fallback_used"]:
        lines.append("> WARNING: fallback predictor used (Agent B module not available).")
        lines.append("")
    lines.append("## Predicted vs observed (per cohort)")
    lines.append("")
    lines.append("| Cohort | N | Observed | Predicted | Diff | Verdict |")
    lines.append("|--------|---|---------:|----------:|-----:|---------|")
    for row in summary["rows"]:
        c = row["cohort"]
        lines.append(
            f"| {c['name']} | {c['n']} | "
            f"{c['observed_pCR']*100:.1f}% | "
            f"{row['predicted']*100:.1f}% | "
            f"{row['diff']*100:+.1f}pp | "
            f"{row['verdict']} |"
        )
    lines.append("")
    lines.append("## Calibration metrics")
    lines.append("")
    lines.append(f"- Cohorts used in calibration: **{summary['n_calibration']}** "
                 f"(Montagna 2024 excluded -- ypN0-only cohort).")
    lines.append(f"- Mean absolute error (MAE) predicted-observed: "
                 f"**{summary['mae']*100:.1f}pp**")
    lines.append(f"- Calibration slope (regress observed on predicted, ideal=1): "
                 f"**{summary['slope']:.3f}**")
    lines.append(f"- Calibration intercept (ideal=0): "
                 f"**{summary['intercept']:.3f}**")
    lines.append("")
    lines.append("Reliability thresholds (per cohort):")
    lines.append("")
    lines.append("- |diff| <= 5pp -> OK (excellent)")
    lines.append("- |diff| <= 10pp -> OK")
    lines.append("- |diff| <= 20pp -> WARN")
    lines.append("- |diff| > 20pp -> FAIL")
    lines.append("")
    lines.append("## Per-subtype breakdown")
    lines.append("")
    for row in summary["rows"]:
        c = row["cohort"]
        lines.append(f"### {c['name']}")
        lines.append("")
        lines.append(f"- Case-mix assumptions: {c['case_mix_notes']}")
        lines.append("")
        lines.append("| Subtype | Weight | Predicted | Weighted contribution |")
        lines.append("|---------|-------:|----------:|----------------------:|")
        for s in row["by_subtype"]:
            lines.append(
                f"| {s['subtype']} | {s['weight']*100:.0f}% | "
                f"{s['predicted_subtype']*100:.1f}% | "
                f"{s['weighted']*100:.1f}pp |"
            )
        lines.append("")
    lines.append("## Cohorts where case-mix data is incomplete")
    lines.append("")
    lines.append("All non-pure-subtype cohorts (Z1071, RISAS, RISAS sub-analysis,"
                 " Montagna 2024) required assumed subtype proportions because the"
                 " publications do not report exact subtype splits. Defaults are"
                 " based on contemporaneous case-mix in similar trials.")
    lines.append("")
    lines.append("- NeoSTOP observed pCR is breast+axilla composite (not pure"
                 " axillary pCR); treat predicted-vs-observed gap as a lower bound.")
    lines.append("- ACOSOG Z1071 has 5% cN2 (modelled as cN1 here -> small upward bias).")
    path.write_text("\n".join(lines), encoding="utf-8")


def main(argv: Optional[list[str]] = None) -> int:
    summary = run()
    report_path = THIS_DIR / "calibration_report.md"
    write_report(summary, report_path)

    # Console summary
    print(f"\nCalibration report written to: {report_path}")
    print(f"  MAE = {summary['mae']*100:.1f}pp  |  slope = {summary['slope']:.3f}"
          f"  |  intercept = {summary['intercept']:.3f}")
    fails = sum(1 for r in summary["rows"]
                if r["verdict"] == "FAIL"
                and r["cohort"]["include_in_calibration"])
    warns = sum(1 for r in summary["rows"]
                if r["verdict"] == "WARN"
                and r["cohort"]["include_in_calibration"])
    print(f"  Verdicts: {fails} FAIL, {warns} WARN, "
          f"{summary['n_calibration'] - fails - warns} OK")
    return 0 if fails == 0 else 1


if __name__ == "__main__":
    raise SystemExit(main())
