#!/usr/bin/env python3
"""
HCC-METASCORE v2 — Biomarker-Driven Composite Scoring Framework
Two-tier model: Standard (routine biomarkers) and Extended (advanced profiling)

Weights derived from published pooled hazard ratios via log-normalization.
Transformation functions use piecewise-linear interpolation anchored to
published clinical thresholds.

PURPOSE: Research focus and multidisciplinary discussion only.
NOT for clinical prescribing, diagnosis, or treatment decisions.

Key references embedded inline per domain.
"""

from __future__ import annotations
import random
from dataclasses import dataclass, field
from typing import Optional, List


# ─────────────────────────────────────────────────────────────────────────────
# Piecewise-linear transformation engine
# ─────────────────────────────────────────────────────────────────────────────

def piecewise_linear(x: float, anchors: list[tuple[float, float]]) -> float:
    """
    Convert a raw clinical value to a 0-100 subscale score using
    piecewise-linear interpolation between (value, score) anchor pairs.
    Clamps to [0, 100].
    """
    if x <= anchors[0][0]:
        return float(anchors[0][1])
    if x >= anchors[-1][0]:
        return float(anchors[-1][1])
    for i in range(len(anchors) - 1):
        x1, s1 = anchors[i]
        x2, s2 = anchors[i + 1]
        if x1 <= x <= x2:
            return s1 + (x - x1) * (s2 - s1) / (x2 - x1)
    return float(anchors[-1][1])


# ─────────────────────────────────────────────────────────────────────────────
# Anchor tables (published thresholds cited in paper)
# ─────────────────────────────────────────────────────────────────────────────

AFP_ANCHORS = [
    # (AFP ng/mL, score)
    # Normal ULN: 20 ng/mL
    # AFP ≥200: Mendez-Blanco et al. 2024 (independent adverse OS factor)
    # AFP ≥400: Lok et al. 2010 (MVI-associated threshold)
    # AFP ≥1000: macrovascular invasion risk
    (20.0,   0.0),
    (200.0, 50.0),
    (400.0, 75.0),
    (1000.0, 100.0),
]

DCP_ANCHORS = [
    # (DCP mAU/mL, score)
    # <40: below portal vein invasion threshold (Imamura et al. 2003)
    # ~85: median in recurrence group (Scientific Reports 2024)
    # ≥400: severe MVI marker
    # ≥1000: unresectable HCC range (Hamzah et al. 2023)
    (40.0,    5.0),
    (85.0,   55.0),
    (400.0,  80.0),
    (1000.0, 100.0),
]

NLR_ANCHORS = [
    # (NLR, score)
    # ≥3: adverse OS threshold (Mendez-Blanco et al. 2024; Ou et al. 2016)
    # >3.8: associated with MVI (Wang et al. 2021)
    # ≥5: strongly adverse per meta-analysis (Ou et al. 2016)
    (2.0,  0.0),
    (2.5, 15.0),
    (3.0, 35.0),
    (3.8, 55.0),
    (5.0, 75.0),
    (7.0, 100.0),
]

VEGF_ANCHORS = [
    # (VEGF pg/mL, score)
    # High VEGF: pooled OS HR 1.85 (Cai et al. 2015)
    (80.0,   0.0),
    (150.0, 25.0),
    (250.0, 55.0),
    (400.0, 80.0),
    (600.0, 100.0),
]

CTC_ANCHORS = [
    # (CTCs per 7.5 mL, score)
    (0,   0.0),
    (1,  15.0),
    (2,  25.0),
    (5,  55.0),
    (10, 80.0),
    (15, 100.0),
]

CTDNA_ANCHORS = [
    # (ctDNA VAF %, score)
    (0.0,   0.0),
    (0.1,  10.0),
    (0.5,  25.0),
    (2.0,  55.0),
    (5.0,  80.0),
    (10.0, 100.0),
]


# ─────────────────────────────────────────────────────────────────────────────
# Published hazard ratios and normalized weights
# ─────────────────────────────────────────────────────────────────────────────

import math

# Standard Tier: HR from published meta-analyses (see paper Section 3)
_STANDARD_HR = {
    "mvi":  2.50,   # avg of 2.21-2.68: Xie 2025, Shirabe 2008
    "vegf": 1.85,   # Cai et al. 2015 meta-analysis
    "afp":  1.80,   # Mendez-Blanco et al. 2024
    "nlr":  1.80,   # Ou et al. 2016 meta-analysis (20,475 pts)
    "dcp":  1.78,   # proxy: Scientific Reports 2024; Imamura 2003
}

# Extended Tier: includes Standard domains + 5 additional
_EXTENDED_HR = {
    **_STANDARD_HR,
    "tme":     1.70,  # proxy: Llovet et al. 2022 (PD-L1 in HCC immunotherapy)
    "ctc_ctdna": 1.75, # proxy: Ye et al. 2019 (CTC detection)
    "genetic": 1.60,  # proxy: Schulze et al. 2015 (TP53 in HCC)
    "emt":     1.50,  # proxy: Ruiz de Galarreta et al. 2019 (CTNNB1 subtype)
    "fgfr":    1.20,  # mechanistic signal only: Ikeda 2021, Shitara 2024
}

def _log_normalize(hr_dict: dict) -> dict:
    log_hrs = {k: math.log(v) for k, v in hr_dict.items()}
    total = sum(log_hrs.values())
    return {k: round(v / total, 4) for k, v in log_hrs.items()}

STANDARD_WEIGHTS = _log_normalize(_STANDARD_HR)
EXTENDED_WEIGHTS = _log_normalize(_EXTENDED_HR)


# ─────────────────────────────────────────────────────────────────────────────
# Patient dataclass
# ─────────────────────────────────────────────────────────────────────────────

@dataclass
class HCCPatient:
    # ── Standard Tier (required for Standard scoring) ──────────────────────
    afp_ng_ml: float = 20.0          # AFP; normal ULN ≈ 20 ng/mL
    afp_l3_pct: float = 0.0          # AFP-L3 fraction %; add-on if ≥15%
    dcp_mau_ml: float = 40.0         # DCP/PIVKA-II in mAU/mL
    nlr: float = 2.5                 # Neutrophil-to-lymphocyte ratio
    vegf_pg_ml: float = 100.0        # Serum VEGF in pg/mL
    mvi_grade: str = "M0"            # "M0", "M1", "M2", "macro", or "unknown"

    # ── Extended Tier (optional; leave None if unavailable) ────────────────
    ctc_count: Optional[int] = None          # CTCs per 7.5 mL blood
    ctdna_vaf_pct: Optional[float] = None    # ctDNA variant allele fraction %
    pd_l1_tps_pct: Optional[float] = None    # PD-L1 TPS %
    til_density: Optional[str] = None        # "low", "moderate", "high"
    ctnnb1_mutation: Optional[bool] = None   # Wnt/beta-catenin mutation
    e_cadherin_loss: Optional[bool] = None   # E-cadherin downregulation
    vimentin_positive: Optional[bool] = None # Vimentin upregulation
    fgfr_amplification: Optional[bool] = None
    tp53_mutation: Optional[bool] = None
    pten_loss: Optional[bool] = None
    rb1_loss: Optional[bool] = None
    epigenetic_dysregulation: Optional[bool] = None


# ─────────────────────────────────────────────────────────────────────────────
# Domain scoring functions
# ─────────────────────────────────────────────────────────────────────────────

def score_afp(afp: float, afp_l3: float) -> tuple[float, str]:
    """
    Piecewise-linear on AFP anchors; AFP-L3 ≥15% adds +15 (capped at 100).
    Anchors: Lok et al. 2010 (Gastroenterology), Mendez-Blanco et al. 2024 (PMC).
    """
    s = piecewise_linear(afp, AFP_ANCHORS)
    detail = f"AFP={afp:.0f} ng/mL"
    if afp_l3 >= 15:
        s = min(s + 15, 100)
        detail += f", AFP-L3={afp_l3:.1f}% (≥15% add-on applied)"
    else:
        detail += f", AFP-L3={afp_l3:.1f}% (<15%)"
    return round(s, 1), detail


def score_dcp(dcp: float) -> tuple[float, str]:
    """
    Piecewise-linear on DCP anchors.
    Anchors: Imamura et al. 2003; Scientific Reports 2024; Hamzah et al. 2023.
    Note: values below 40 mAU/mL return a floor score of 5 (not 0) since
    DCP may be mildly elevated even in early HCC.
    """
    if dcp < 40:
        return 5.0, f"DCP={dcp:.0f} mAU/mL [below threshold <40]"
    s = piecewise_linear(dcp, DCP_ANCHORS)
    return round(s, 1), f"DCP={dcp:.0f} mAU/mL"


def score_mvi(grade: str) -> tuple[float, str]:
    """
    Categorical: M0=0, M1=50, M2=80, macro=100.
    Ribero et al. 2024 (5-yr OS data); Xie et al. 2025 (HR 2.68).
    Returns (score, detail, is_missing).
    """
    mapping = {
        "M0":    (0.0,   "M0: no MVI — 5-yr OS ~60.7% [Ribero et al. 2024]"),
        "M1":    (50.0,  "M1: mild MVI — 5-yr OS ~57.4%"),
        "M2":    (80.0,  "M2: severe MVI — 5-yr OS ~29.7% [Ribero et al. 2024]"),
        "macro": (100.0, "Macrovascular invasion — BCLC stage C [EASL 2018]"),
    }
    if grade.lower() == "unknown":
        return (None, "MVI status unknown — domain excluded from composite", True)
    result = mapping.get(grade.upper(), mapping.get(grade.lower()))
    if result is None:
        return (None, f"Unrecognised MVI grade '{grade}' — domain excluded", True)
    score, detail = result
    return (score, detail, False)


def score_nlr(nlr: float) -> tuple[float, str]:
    """
    Piecewise-linear on NLR anchors.
    NLR ≥3 threshold: Mendez-Blanco et al. 2024; pooled HR 1.80: Ou et al. 2016.
    NLR >3.8 and MVI: Wang et al. 2021.
    """
    s = piecewise_linear(nlr, NLR_ANCHORS)
    return round(s, 1), f"NLR={nlr:.1f}"


def score_vegf(vegf: float) -> tuple[float, str]:
    """
    Piecewise-linear on VEGF anchors.
    High VEGF OS HR 1.85 (Cai et al. 2015, meta-analysis of sorafenib-treated HCC).
    """
    s = piecewise_linear(vegf, VEGF_ANCHORS)
    return round(s, 1), f"VEGF={vegf:.0f} pg/mL"


def score_ctc_ctdna(ctc: Optional[int], vaf: Optional[float]) -> tuple[float, str]:
    """
    Combined CTC and ctDNA VAF signal.
    CTC detection: Ye et al. 2019 (Molecular Cancer) review.
    """
    if ctc is None and vaf is None:
        return (None, "CTC/ctDNA not tested — domain excluded", True)
    s = 0.0
    parts = []
    if ctc is not None:
        cs = piecewise_linear(float(ctc), CTC_ANCHORS)
        s += cs * 0.5
        parts.append(f"CTCs={ctc}/7.5mL")
    if vaf is not None:
        vs = piecewise_linear(vaf, CTDNA_ANCHORS)
        s += vs * 0.5
        parts.append(f"ctDNA VAF={vaf:.1f}%")
    if ctc is not None and vaf is None:
        s = s * 2  # scale to full range if only one component
    elif vaf is not None and ctc is None:
        s = s * 2
    return round(min(s, 100), 1), ", ".join(parts), False


def score_tme(pd_l1: Optional[float], til: Optional[str],
              nlr_score: float, ctnnb1: Optional[bool]) -> tuple:
    """
    PD-L1 TPS ≥10% adds 50; ≥1% adds 20.
    TIL-high adds 30; TIL-moderate adds 15.
    NLR score >65 (NLR >5) penalises by -20 (systemic inflammation).
    CTNNB1 mutation noted as caveat (Ruiz de Galarreta et al. 2019).
    Llovet et al. 2022 (Nat Rev Clin Oncol).
    """
    if pd_l1 is None and til is None:
        return (None, "TME not assessed — domain excluded", True)
    s = 0.0
    parts = []
    if pd_l1 is not None:
        if pd_l1 >= 10:
            s += 50
        elif pd_l1 >= 1:
            s += 20
        parts.append(f"PD-L1 TPS={pd_l1:.0f}%")
    if til is not None:
        if til.lower() == "high":
            s += 30
        elif til.lower() == "moderate":
            s += 15
        parts.append(f"TIL={til}")
    if nlr_score > 65:
        s = max(s - 20, 0)
        parts.append("NLR penalty applied (NLR >5)")
    if ctnnb1:
        parts.append("CTNNB1 mut [note: associated with reduced immune infiltration]")
    return round(min(s, 100), 1), ", ".join(parts), False


def score_emt(e_cad: Optional[bool], vim: Optional[bool],
              ctnnb1: Optional[bool]) -> tuple:
    """
    E-cadherin loss: +35; vimentin: +35; CTNNB1 mutation: +25.
    Schulze et al. 2015 (Nat Genet); Ruiz de Galarreta et al. 2019 (Cancer Discov).
    """
    if e_cad is None and vim is None and ctnnb1 is None:
        return (None, "EMT markers not assessed — domain excluded", True)
    s = 0.0
    parts = []
    if e_cad:
        s += 35; parts.append("E-cad loss")
    if vim:
        s += 35; parts.append("vimentin+")
    if ctnnb1:
        s += 25; parts.append("CTNNB1 mut")
    if not parts:
        parts.append("No EMT markers detected")
    return round(min(s, 100), 1), ", ".join(parts), False


def score_fgfr(amp: Optional[bool]) -> tuple:
    """
    FGFR amplification/dysregulation present → 85.
    Absent → 0.
    Ikeda et al. 2021 (REFLECT biomarker analysis); Shitara et al. 2024.
    Note: lowest weight (0.03) reflecting mechanistic signal only, not validated OS HR.
    """
    if amp is None:
        return (None, "FGFR status not tested — domain excluded", True)
    if amp:
        return 85.0, "FGFR amplification/dysregulation detected", False
    return 0.0, "No FGFR amplification", False


def score_genetic(tp53: Optional[bool], pten: Optional[bool],
                  rb1: Optional[bool]) -> tuple:
    """
    TP53 mutation: +40; PTEN loss: +35; RB1 loss: +25.
    Schulze et al. 2015 (Nat Genet); TCGA-LIHC (Cancer Cell 2017).
    """
    if tp53 is None and pten is None and rb1 is None:
        return (None, "Genetic drivers not tested — domain excluded", True)
    s = 0.0
    parts = []
    if tp53:
        s += 40; parts.append("TP53 mut")
    if pten:
        s += 35; parts.append("PTEN loss")
    if rb1:
        s += 25; parts.append("RB1 loss")
    if not parts:
        parts.append("No high-risk driver mutations detected in tested genes")
    return round(min(s, 100), 1), ", ".join(parts), False


# ─────────────────────────────────────────────────────────────────────────────
# Composite computation
# ─────────────────────────────────────────────────────────────────────────────

def compute_composite(scores: dict, weights: dict, present_domains: set) -> float:
    """
    Compute composite score over present (non-missing) domains,
    re-normalizing weights to those domains only.
    """
    present_weights = {k: v for k, v in weights.items() if k in present_domains}
    total_weight = sum(present_weights.values())
    if total_weight == 0:
        return 0.0
    composite = sum(scores[k] * present_weights[k] / total_weight
                    for k in present_domains)
    return round(min(composite, 100), 1)


def pss_formulas(scores: dict, tier: str) -> dict:
    """Compute Pathway Signal Scores for each systemic agent."""
    mvi_s = scores.get("mvi", 0) or 0
    afp_s = scores.get("afp", 0) or 0
    vegf_s = scores.get("vegf", 0) or 0
    fgfr_s = scores.get("fgfr", 0) or 0
    tme_s = scores.get("tme", 0) or 0
    ctc_s = scores.get("ctc_ctdna", 0) or 0
    nlr_s = scores.get("nlr", 0) or 0

    if tier == "extended":
        pss_sor = (vegf_s*0.50) + (afp_s*0.30) + (mvi_s*0.20)
        pss_len = (vegf_s*0.40) + (fgfr_s*0.30) + (afp_s*0.20) + (mvi_s*0.10)
        pss_atz = (tme_s*0.55) + (vegf_s*0.30) + (ctc_s*0.15)
    else:
        # Standard tier: no TME, FGFR, ctDNA
        nlr_inv = max(0, 100 - nlr_s)  # lower NLR = better immune milieu
        pss_sor = (vegf_s*0.50) + (afp_s*0.30) + (mvi_s*0.20)
        pss_len = (vegf_s*0.60) + (afp_s*0.25) + (mvi_s*0.15)
        pss_atz = (vegf_s*0.55) + (nlr_inv*0.45)

    def level(x):
        if x >= 60: return "VERY HIGH"
        if x >= 40: return "HIGH"
        if x >= 20: return "MODERATE"
        return "LOW"

    return {
        "Sorafenib":  (round(pss_sor, 1), level(pss_sor)),
        "Lenvatinib": (round(pss_len, 1), level(pss_len)),
        "Atezo/Bev":  (round(pss_atz, 1), level(pss_atz)),
    }


def compute_hcc_score(patient: HCCPatient, n_sims: int = 5000, seed: int = 42):
    """Main scoring function. Returns composite score, CI, PSP, and domain details."""
    # ── Determine tier ─────────────────────────────────────────────────────
    has_extended = any([
        patient.ctc_count is not None,
        patient.ctdna_vaf_pct is not None,
        patient.pd_l1_tps_pct is not None,
        patient.fgfr_amplification is not None,
        patient.tp53_mutation is not None,
    ])
    tier = "extended" if has_extended else "standard"
    weights = EXTENDED_WEIGHTS if has_extended else STANDARD_WEIGHTS

    # ── Score each domain ───────────────────────────────────────────────────
    afp_s, afp_d = score_afp(patient.afp_ng_ml, patient.afp_l3_pct)
    dcp_s, dcp_d = score_dcp(patient.dcp_mau_ml)
    mvi_result    = score_mvi(patient.mvi_grade)
    nlr_s, nlr_d  = score_nlr(patient.nlr)
    vegf_s, vegf_d = score_vegf(patient.vegf_pg_ml)

    mvi_s    = mvi_result[0]
    mvi_d    = mvi_result[1]
    mvi_miss = mvi_result[2] if len(mvi_result) > 2 else False

    scores = {"afp": afp_s, "dcp": dcp_s, "nlr": nlr_s, "vegf": vegf_s}
    details = {"afp": afp_d, "dcp": dcp_d, "nlr": nlr_d, "vegf": vegf_d}
    if not mvi_miss:
        scores["mvi"] = mvi_s
        details["mvi"] = mvi_d

    if has_extended:
        ctc_result = score_ctc_ctdna(patient.ctc_count, patient.ctdna_vaf_pct)
        tme_result = score_tme(patient.pd_l1_tps_pct, patient.til_density,
                               nlr_s, patient.ctnnb1_mutation)
        emt_result = score_emt(patient.e_cadherin_loss, patient.vimentin_positive,
                               patient.ctnnb1_mutation)
        fgfr_result = score_fgfr(patient.fgfr_amplification)
        gen_result = score_genetic(patient.tp53_mutation, patient.pten_loss,
                                   patient.rb1_loss)

        for key, result in [("ctc_ctdna", ctc_result), ("tme", tme_result),
                             ("emt", emt_result), ("fgfr", fgfr_result),
                             ("genetic", gen_result)]:
            s, d = result[0], result[1]
            missing = result[2] if len(result) > 2 else False
            if not missing:
                scores[key] = s
                details[key] = d

    present_domains = set(scores.keys())
    composite = compute_composite(scores, weights, present_domains)

    # ── Monte Carlo CI (continuous inputs only) ─────────────────────────────
    rng = random.Random(seed)
    sims = []
    for _ in range(n_sims):
        def perturb(v, cv=0.12):
            return max(0.0, v * (1 + rng.gauss(0, cv)))
        noisy_afp = perturb(patient.afp_ng_ml)
        noisy_afp_l3 = min(100, perturb(patient.afp_l3_pct))
        noisy_dcp = perturb(patient.dcp_mau_ml)
        noisy_nlr = perturb(patient.nlr, 0.10)
        noisy_vegf = perturb(patient.vegf_pg_ml)
        noisy_scores = {
            "afp":  score_afp(noisy_afp, noisy_afp_l3)[0],
            "dcp":  score_dcp(noisy_dcp)[0],
            "nlr":  score_nlr(noisy_nlr)[0],
            "vegf": score_vegf(noisy_vegf)[0],
        }
        if not mvi_miss:
            noisy_scores["mvi"] = mvi_s  # categorical, not perturbed
        if has_extended:
            noisy_ctc_vaf = perturb(patient.ctdna_vaf_pct or 0, 0.15)
            noisy_ctc_result = score_ctc_ctdna(patient.ctc_count, noisy_ctc_vaf)
            if len(noisy_ctc_result) > 2 and not noisy_ctc_result[2]:
                noisy_scores["ctc_ctdna"] = noisy_ctc_result[0]
            for k in ["tme", "emt", "fgfr", "genetic"]:
                if k in scores:
                    noisy_scores[k] = scores[k]  # binary inputs, not perturbed
        sims.append(compute_composite(noisy_scores, weights,
                                      set(noisy_scores.keys())))
    sims.sort()
    ci_lower = round(sims[int(0.025 * n_sims)], 1)
    ci_upper = round(sims[int(0.975 * n_sims)], 1)

    # ── Category ────────────────────────────────────────────────────────────
    if composite < 20:   category = "LOW"
    elif composite < 40: category = "MODERATE"
    elif composite < 60: category = "HIGH"
    else:                category = "VERY HIGH"

    psp = pss_formulas(scores, tier)

    # ── Interpretive notes ──────────────────────────────────────────────────
    notes = []
    if mvi_miss:
        notes.append("MVI status unavailable — weight re-normalized across remaining domains. Score may underestimate risk if MVI is present.")
    if scores.get("fgfr", 0) >= 70 :
        notes.append("FGFR amplification/dysregulation detected — mechanistically differentiates lenvatinib from sorafenib [Shitara et al. 2024; Ikeda et al. 2021].")
    if scores.get("tme", 0) >= 50 and scores.get("vegf", 0) >= 40:
        notes.append("Elevated immune and angiogenic signals co-occur — profile may be relevant to combined anti-PD-L1/anti-VEGF mechanistic exploration.")
    if patient.ctnnb1_mutation and scores.get("tme", 0) < 30:
        notes.append("CTNNB1 mutation + low immune signal: consistent with Wnt-activated, non-inflamed TME phenotype [Ruiz de Galarreta et al. 2019]. This does not exclude immunotherapy consideration but is noted as a research signal.")
    if scores.get("nlr", 0) < 20 and not has_extended:
        notes.append("Low NLR proxy suggests favourable systemic immune milieu — Extended Tier profiling (PD-L1, TIL density) recommended to refine Atezo/Bev signal.")
    if tier == "standard":
        notes.append("Standard Tier only. Atezo/Bev PSP is estimated from NLR proxy. Extended Tier data (PD-L1, TIL, FGFR) would substantially improve pathway signal resolution.")

    return {
        "composite_score": composite,
        "ci_lower": ci_lower,
        "ci_upper": ci_upper,
        "category": category,
        "tier": tier,
        "domains_scored": list(present_domains),
        "domain_details": [
            {"domain": k, "raw_score": scores[k],
             "weight": weights.get(k, "re-normalized"),
             "detail": details[k]}
            for k in present_domains if k in details
        ],
        "pathway_signal_profile": psp,
        "interpretive_notes": notes,
    }


# ─────────────────────────────────────────────────────────────────────────────
# Output printer
# ─────────────────────────────────────────────────────────────────────────────

def print_result(result: dict, label: str):
    SEP = "=" * 72
    print(f"\n{SEP}\n{label}\n{SEP}")
    print(f"Tier: {result['tier'].upper()}")
    print(f"Composite score: {result['composite_score']}/100 [{result['category']}]")
    print(f"95% CI: [{result['ci_lower']}, {result['ci_upper']}]  "
          f"(reflects continuous input measurement variability only)")
    print(f"\nDomains scored ({len(result['domains_scored'])}): "
          f"{', '.join(result['domains_scored'])}")
    print("\nDomain breakdown:")
    for d in result["domain_details"]:
        print(f"  {d['domain']:15s}  score={d['raw_score']:5.1f}  | {d['detail']}")
    print("\nPathway Signal Profile:")
    for agent, (score, level) in result["pathway_signal_profile"].items():
        print(f"  {agent:12s}: {score:.1f} [{level}]")
    if result["interpretive_notes"]:
        print("\nInterpretive notes (hypothesis-generating only):")
        for note in result["interpretive_notes"]:
            print(f"  * {note}")
    print(f"\n{'─'*72}")
    print("PURPOSE: Research focus and multidisciplinary discussion only.")
    print("NOT for clinical prescribing, diagnosis, or treatment decisions.")
    print(f"{'─'*72}")


# ─────────────────────────────────────────────────────────────────────────────
# Demo: three scenarios from published cohort summary statistics
# ─────────────────────────────────────────────────────────────────────────────

def demo():
    scenarios = [
        (
            "Scenario 1 — Standard Tier (low-burden profile)\n"
            "  Source: Favourable subgroup of Mendez-Blanco et al. 2024 sorafenib cohort\n"
            "  (0-1 adverse factors; median OS 17.4 months)",
            HCCPatient(
                afp_ng_ml=45, afp_l3_pct=6, dcp_mau_ml=28,
                mvi_grade="M1", nlr=2.2, vegf_pg_ml=140,
            ),
        ),
        (
            "Scenario 2 — Standard Tier (high-burden, angiogenic-dominant)\n"
            "  Source: Worst prognostic subgroup (3 adverse factors) of Mendez-Blanco\n"
            "  et al. 2024; values consistent with Hamzah et al. 2023 unresectable HCC cohort",
            HCCPatient(
                afp_ng_ml=1240, afp_l3_pct=34, dcp_mau_ml=520,
                mvi_grade="macro", nlr=5.8, vegf_pg_ml=420,
            ),
        ),
        (
            "Scenario 3 — Extended Tier (immune-inflamed + FGFR signal)\n"
            "  Source: Profile consistent with Atezo/Bev arm responder subgroup\n"
            "  characteristics from Finn et al. IMbrave150 2020 (NEJM)",
            HCCPatient(
                afp_ng_ml=310, afp_l3_pct=9, dcp_mau_ml=88,
                mvi_grade="M1", nlr=2.3, vegf_pg_ml=195,
                ctc_count=2, ctdna_vaf_pct=1.8,
                pd_l1_tps_pct=12, til_density="high",
                ctnnb1_mutation=False, e_cadherin_loss=True,
                vimentin_positive=False, fgfr_amplification=True,
                tp53_mutation=False, pten_loss=False, rb1_loss=False,
                epigenetic_dysregulation=False,
            ),
        ),
    ]
    for label, patient in scenarios:
        result = compute_hcc_score(patient)
        print_result(result, label)


if __name__ == "__main__":
    demo()