"""Plotting module for the Level-2 ypN0 nomogram.

Outputs (in plots/):
    forest_subtype.png   - per-subtype pooled rates with 95% CI
    calibration.png      - predicted vs observed for validation studies
    nomogram.png         - classic-style horizontal nomogram

Run after meta_regression.py has written results_pooled_subtype.csv,
results_validation.csv and data/nomogram_coefficients.csv.
"""

from __future__ import annotations

from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

HERE = Path(__file__).resolve().parent
PLOTS = HERE / "plots"
PLOTS.mkdir(exist_ok=True)
DATA = HERE.parent / "data"


def forest_subtype():
    df = pd.read_csv(HERE / "results_pooled_subtype.csv")
    df = df.sort_values("pooled_rate", ascending=True).reset_index(drop=True)
    fig, ax = plt.subplots(figsize=(9, 5))
    y = np.arange(len(df))
    ax.errorbar(df["pooled_rate"], y,
                xerr=[df["pooled_rate"] - df["ci_low"],
                      df["ci_high"] - df["pooled_rate"]],
                fmt="s", color="#2b6cb0", ecolor="#7b9fc7",
                markersize=10, capsize=4, linewidth=2)
    for i, row in df.iterrows():
        ax.text(row["ci_high"] + 0.02, i,
                f"{row['pooled_rate']*100:0.1f}% "
                f"({row['ci_low']*100:0.1f}-{row['ci_high']*100:0.1f}%)  "
                f"k={int(row['k_studies'])}, I2={row['i2_percent']:0.0f}%",
                va="center", fontsize=9)
    ax.set_yticks(y)
    ax.set_yticklabels(df["subtype"])
    ax.set_xlim(0, 1.0)
    ax.set_xlabel("Pooled ypN0 rate (random-effects DL)")
    ax.set_title("Forest plot - Pooled ypN0 rate by intrinsic subtype\n"
                 "DerSimonian-Laird random-effects on logit-transformed proportions")
    ax.grid(axis="x", alpha=0.3)
    fig.tight_layout()
    out = PLOTS / "forest_subtype.png"
    fig.savefig(out, dpi=160)
    plt.close(fig)
    return out


def calibration():
    df = pd.read_csv(HERE / "results_validation.csv")
    # Exclude selection-defined Montagna
    df_cal = df[df["observed_ypN0"] < 1.0].copy()
    fig, ax = plt.subplots(figsize=(7, 7))
    ax.plot([0, 1], [0, 1], "--", color="gray", label="Perfect calibration")
    ax.scatter(df_cal["predicted_ypN0"], df_cal["observed_ypN0"],
               s=120, c="#c53030", edgecolor="black", zorder=3)
    for _, r in df_cal.iterrows():
        ax.annotate(r["study"], (r["predicted_ypN0"], r["observed_ypN0"]),
                    xytext=(7, 7), textcoords="offset points",
                    fontsize=8.5)
    ax.set_xlim(0, 1)
    ax.set_ylim(0, 1)
    ax.set_xlabel("Predicted ypN0 (Level-2 nomogram)")
    ax.set_ylabel("Observed ypN0 in cohort")
    brier = float(np.mean((df_cal["predicted_ypN0"] -
                           df_cal["observed_ypN0"]) ** 2))
    ax.set_title("Calibration - Predicted vs observed ypN0 rates\n"
                 f"Brier score (cohort-level) = {brier:0.4f}")
    ax.grid(alpha=0.3)
    fig.tight_layout()
    out = PLOTS / "calibration.png"
    fig.savefig(out, dpi=160)
    plt.close(fig)
    return out


def _logit(p):
    return np.log(p / (1 - p))


def _inv_logit(x):
    return 1.0 / (1.0 + np.exp(-x))


def nomogram():
    """Classic Harrell-style nomogram.

    Each predictor gets a horizontal axis whose tick positions are scaled so
    that 1 unit of "points" equals a fixed amount of beta (log-odds). The
    maximum-effect predictor is anchored at 100 points.
    """
    coef = pd.read_csv(DATA / "nomogram_coefficients.csv")
    coef = coef.set_index("term")

    # Build axis spec: (label, list of (value_label, beta))
    # Subtype as a 4-level axis
    axes_spec = []
    axes_spec.append((
        "Subtipo molecular",
        [
            ("HR+/HER2- (ref)", 0.0),
            ("TNBC",  float(coef.loc["subtype_TNBC", "beta"])),
            ("HR+/HER2+", float(coef.loc["subtype_HR+/HER2+", "beta"])),
            ("HR-/HER2+", float(coef.loc["subtype_HR-/HER2+", "beta"])),
        ],
    ))
    axes_spec.append((
        "Ki67 (solo HR+/HER2-)",
        [
            ("<20%", 0.0),
            ("20-30%", float(coef.loc["ki67_20_30", "beta"])),
            (">=30%",  float(coef.loc["ki67_gt30", "beta"])),
        ],
    ))
    axes_spec.append((
        "Riesgo genomico\n(HR+/HER2-)",
        [
            ("Bajo (ref)", 0.0),
            ("Alto (RS>25 / MP-H)", float(coef.loc["genomic_high", "beta"])),
        ],
    ))
    axes_spec.append((
        "Edad (anos)",
        [
            ("70", 2.0 * float(coef.loc["age_per_decade", "beta"])),
            ("60", 1.0 * float(coef.loc["age_per_decade", "beta"])),
            ("50 (ref)", 0.0),
            ("40", -1.0 * float(coef.loc["age_per_decade", "beta"])),
            ("30", -2.0 * float(coef.loc["age_per_decade", "beta"])),
        ],
    ))
    axes_spec.append((
        "Estadio clinico cN",
        [
            ("cN1 (ref)", 0.0),
            ("cN2-3", float(coef.loc["cN_stage_cN23", "beta"])),
        ],
    ))
    axes_spec.append((
        "Pembrolizumab\n(solo TNBC)",
        [("No (ref)", 0.0), ("Si", float(coef.loc["regimen_pembro", "beta"]))],
    ))
    axes_spec.append((
        "Pertuzumab\n(solo HER2+)",
        [("No (ref)", 0.0), ("Si", float(coef.loc["regimen_pertuzumab", "beta"]))],
    ))
    axes_spec.append((
        "Carboplatino\n(solo TNBC)",
        [("No (ref)", 0.0), ("Si", float(coef.loc["regimen_carbo", "beta"]))],
    ))

    # Determine the scaling: 100 points = largest absolute beta swing
    max_swing = max(max(b for _, b in ax[1]) - min(b for _, b in ax[1])
                    for ax in axes_spec)
    pts_per_beta = 100.0 / max_swing

    intercept = float(coef.loc["intercept", "beta"])

    # Total points axis: range of possible points
    sum_min = sum(min(b for _, b in ax[1]) for ax in axes_spec)
    sum_max = sum(max(b for _, b in ax[1]) for ax in axes_spec)
    total_min = sum_min * pts_per_beta
    total_max = sum_max * pts_per_beta

    # Probability axis: from intercept + sum_min to intercept + sum_max
    p_min_logit = intercept + sum_min
    p_max_logit = intercept + sum_max

    # Draw the nomogram
    n_axes = len(axes_spec)
    fig, ax = plt.subplots(figsize=(13, 1.2 + 0.7 * (n_axes + 2)))
    ax.set_xlim(-5, 110)
    ax.set_ylim(0, n_axes + 2.5)
    ax.invert_yaxis()

    # Points reference axis at top
    pts_y = 0.4
    ax.hlines(pts_y, 0, 100, color="black", lw=1.2)
    for pt in range(0, 101, 10):
        ax.vlines(pt, pts_y - 0.08, pts_y + 0.08, color="black", lw=1)
        ax.text(pt, pts_y - 0.18, f"{pt}", ha="center", va="bottom", fontsize=8)
    ax.text(-2, pts_y, "Puntos", ha="right", va="center",
            fontweight="bold", fontsize=10)

    # Each predictor axis
    for i, (label, levels) in enumerate(axes_spec):
        y = i + 1.5
        betas = [b for _, b in levels]
        beta_min, beta_max = min(betas), max(betas)
        # Map betas to x in 0..100 relative to local axis (each axis starts at 0)
        # We'll center each axis at its own scale starting from 0 points
        # Place each level: x = (beta - beta_min) * pts_per_beta
        ax.hlines(y, 0, (beta_max - beta_min) * pts_per_beta,
                  color="#2b6cb0", lw=1.2)
        ax.text(-2, y, label, ha="right", va="center", fontsize=9)
        for value_label, b in levels:
            x = (b - beta_min) * pts_per_beta
            ax.vlines(x, y - 0.12, y + 0.12, color="#2b6cb0", lw=1)
            ax.text(x, y + 0.22, value_label, ha="center", va="top",
                    fontsize=7.5, color="#1a365d")

    # Total points axis
    tot_y = n_axes + 1.5
    ax.hlines(tot_y, 0, 100, color="black", lw=1.2)
    n_ticks = 11
    for j in range(n_ticks):
        x = j * 10
        pts_value = total_min + (total_max - total_min) * j / (n_ticks - 1)
        ax.vlines(x, tot_y - 0.08, tot_y + 0.08, color="black", lw=1)
        ax.text(x, tot_y - 0.18, f"{pts_value:0.0f}", ha="center",
                va="bottom", fontsize=8)
    ax.text(-2, tot_y, "Puntos totales", ha="right", va="center",
            fontweight="bold", fontsize=10)

    # Probability axis (probabilities at common cut-offs)
    prob_y = tot_y + 0.8
    target_probs = [0.05, 0.10, 0.20, 0.30, 0.50, 0.70, 0.85, 0.95]
    # Compute x position for each probability
    x_per_total = 100.0 / (total_max - total_min) if total_max > total_min else 1
    ax.hlines(prob_y, 0, 100, color="#c53030", lw=1.4)
    for p in target_probs:
        target_logit = np.log(p / (1 - p))
        # total points such that intercept + (points / pts_per_beta) ... wait
        # Convert: P -> logit -> "betaSum" needed = logit - intercept;
        # but our total axis spans betaSum from sum_min to sum_max.
        # So total_pts = (betaSum - sum_min) * pts_per_beta
        betaSum = target_logit - intercept
        total_pts = (betaSum - sum_min) * pts_per_beta
        if 0 <= total_pts <= 100:
            ax.vlines(total_pts, prob_y - 0.10, prob_y + 0.10,
                      color="#c53030", lw=1.4)
            ax.text(total_pts, prob_y + 0.25, f"{p:0.2f}", ha="center",
                    va="top", fontsize=8.5, color="#c53030",
                    fontweight="bold")
    ax.text(-2, prob_y, "P(ypN0)", ha="right", va="center",
            fontweight="bold", color="#c53030", fontsize=10)

    ax.set_axis_off()
    fig.suptitle("Nomograma Nivel-2 - Probabilidad de ypN0 tras NAST en cN+ M0\n"
                 "Coeficientes agrupados (random-effects DL) - "
                 "ref: HR+/HER2-, Ki67<20%, edad 50, cN1, sin pembro/pertuzumab/carbo",
                 fontsize=11, y=0.99)
    fig.tight_layout(rect=[0, 0, 1, 0.96])
    out = PLOTS / "nomogram.png"
    fig.savefig(out, dpi=160)
    plt.close(fig)
    return out


if __name__ == "__main__":
    f1 = forest_subtype()
    f2 = calibration()
    f3 = nomogram()
    print(f"Forest plot:    {f1}")
    print(f"Calibration:    {f2}")
    print(f"Nomogram:       {f3}")
