#!/usr/bin/env python
"""Stage 26 — cross-experiment comparison of yEHDE binary-forecast skill (self-contained HTML).

Consolidates the leakage-safe (LOYO-pooled, per-fold decision threshold) binary contingency
inventory of EVERY experiment in the project into one report, focused on the scores the
hybrid-design note asks for: POD, POFD, FAR, FBI (frequency bias), CSI, ETS, HSS, Peirce,
EDI, SEDI — for the primary target yEHDE (>=1 EHD in a lead block).

Sources (all already computed leakage-safe by their own verify stages):
  v1 headline     results/tables/cv_binary_scores_by_block.csv
  v2 zoo          results/tables/v2_cv_scores_by_model_block.csv        (best-by-SEDI/block)
  v2 ensemble     results/tables/v2_scores_with_ensemble.csv            (ens_stack)
  Stage-2         results/tables/v3cal_binary_by_block.csv              (calibration-recovery stack)
  v3 zoo          results/tables/v3zoo_scores_by_model_block.csv        (best-by-SEDI/block)
  §6 vote         results/tables/v6_vote_scores_by_block.csv            (m/n vote ensemble)
  v3 stack        outputs/loyo_pred_v3stack/v3stack_pooled.npz          (13-member; inventory
                  computed here at the pooled Peirce threshold, cross-checked vs
                  v3_stack_scores_by_block.csv)

Output: results/experiment_comparison.html   (no external assets; inline SVG charts).

NOTE on comparability: per-fold-thresholded experiments (v1/v2/Stage-2/v3zoo/§6) never let a
threshold see its own test year. The v3-stack row uses a single pooled Peirce threshold per
block (reproducing stage-23's own SEDI/Peirce), which is marginally optimistic vs a per-fold
threshold; this is footnoted in the report.
"""
import csv, os, sys, html
import numpy as np

ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
TAB = os.path.join(ROOT, "results", "tables")
OUT = os.path.join(ROOT, "results", "experiment_comparison.html")
BLOCKS = ["W1", "W2", "W3", "W4", "W5", "W6"]

sys.path.insert(0, os.path.join(ROOT, "src"))
from ehd import verify as vf  # noqa: E402

# canonical metric -> (label, direction: +1 higher-better / -1 lower-better, ideal)
METRICS = {
    "POD":    ("POD",   +1, 1.0),
    "POFD":   ("POFD",  -1, 0.0),
    "FAR":    ("FAR",   -1, 0.0),
    "FBI":    ("FBI",    0, 1.0),   # frequency bias, ideal = 1
    "CSI":    ("CSI",   +1, 1.0),
    "EDI":    ("EDI",   +1, 1.0),
    "SEDI":   ("SEDI",  +1, 1.0),
    "ETS":    ("ETS",   +1, 1.0),
    "HSS":    ("HSS",   +1, 1.0),
    "Peirce": ("Peirce",+1, 1.0),
}
FOCUS = ["POD", "POFD", "FAR", "FBI", "EDI", "SEDI", "CSI"]   # user's primary focus
ALLCOL = ["POD", "POFD", "FAR", "FBI", "CSI", "ETS", "HSS", "Peirce", "EDI", "SEDI"]
SRC = {  # source-column -> canonical
    "POD": "POD", "POFD": "POFD", "FAR": "FAR", "CSI": "CSI", "ETS": "ETS",
    "HSS": "HSS", "Peirce": "Peirce", "EDI": "EDI", "SEDI": "SEDI",
    "freq_bias": "FBI", "FBI": "FBI",
}


def read_csv(path):
    with open(path) as fh:
        return list(csv.DictReader(fh))


def canon(row):
    """Extract canonical metric dict from a source row (missing -> None)."""
    out = {}
    for sc, cc in SRC.items():
        if sc in row and row[sc] not in ("", None):
            try:
                out[cc] = float(row[sc])
            except ValueError:
                pass
    for m in METRICS:
        out.setdefault(m, None)
    return out


def block_of(row):
    b = row.get("block", "").strip()
    return b if b in BLOCKS else None


# ---------------------------------------------------------------- gather records
# records[exp_key] = {"label":.., "note":.., "rows": {block: {"model":.., metrics..}},
#                     "all": [ (model, block, metricdict) ... ]}
experiments = []


def best_by_sedi(rows, model_field):
    """From many-model rows, keep best-SEDI model per block."""
    per = {}
    for r in rows:
        b = block_of(r)
        if not b:
            continue
        c = canon(r)
        s = c.get("SEDI")
        if s is None:
            continue
        if b not in per or s > per[b][1]:
            per[b] = (r.get(model_field, "?"), s, c, r.get("base_rate"))
    return {b: {"model": v[0], "metrics": v[2], "base_rate": v[3]} for b, v in per.items()}


def full_rows(rows, model_field):
    out = []
    for r in rows:
        b = block_of(r)
        if not b:
            continue
        out.append((r.get(model_field, "?"), b, canon(r), r.get("base_rate")))
    return out


def add(label, note, primary, all_rows):
    experiments.append({"label": label, "note": note, "rows": primary, "all": all_rows})


# --- v1 headline (single model / block) ---
p = os.path.join(TAB, "cv_binary_scores_by_block.csv")
if os.path.exists(p):
    rr = read_csv(p)
    prim = {block_of(r): {"model": "hgb_recal", "metrics": canon(r),
                          "base_rate": r.get("base_rate")}
            for r in rr if block_of(r)}
    add("v1 headline", "hgb + isotonic recalibration; forecast-Tx statistics only (Stage 1).",
        prim, [("hgb_recal", b, prim[b]["metrics"], prim[b]["base_rate"]) for b in prim])

# --- v2 zoo (best-by-SEDI across family x cutoff) ---
p = os.path.join(TAB, "v2_cv_scores_by_model_block.csv")
if os.path.exists(p):
    rr = read_csv(p)
    for r in rr:
        r["_m"] = f"{r.get('family','?')}·c{r.get('cutoff','?')}"
    add("v2 zoo (best)", "Severity predictand + expanded zoo + cascade; best family×cutoff per block.",
        best_by_sedi(rr, "_m"), full_rows(rr, "_m"))

# --- v2 ensemble (ens_stack) ---
p = os.path.join(TAB, "v2_scores_with_ensemble.csv")
if os.path.exists(p):
    rr = read_csv(p)
    prim = {}
    for r in rr:
        b = block_of(r)
        if b and r.get("model") == "ens_stack":
            prim[b] = {"model": "ens_stack", "metrics": canon(r), "base_rate": r.get("base_rate")}
    add("v2 ens_stack", "v2 stacked ensemble (9 members incl. CNN/ANN).",
        prim, full_rows(rr, "model"))

# --- Stage-2 calibration-recovery stack ---
p = os.path.join(TAB, "v3cal_binary_by_block.csv")
if os.path.exists(p):
    rr = read_csv(p)
    prim = {}
    for r in rr:
        b = block_of(r)
        if b and r.get("model") == "stack":
            prim[b] = {"model": "stack", "metrics": canon(r), "base_rate": r.get("base_rate")}
    add("Stage-2 stack", "Lead-dependent extended predictors (ENSO/MJO/Z500/persistence) + OOF calibration stack.",
        prim, full_rows(rr, "model"))

# --- v3 zoo (best-by-SEDI) ---
p = os.path.join(TAB, "v3zoo_scores_by_model_block.csv")
if os.path.exists(p):
    rr = read_csv(p)
    add("v3 zoo (best single)", "Stage-3 hybrid zoo on lead-dependent predictors; best single family per block.",
        best_by_sedi(rr, "family"), full_rows(rr, "family"))

# --- v4 full-field true-spatial ResU-Net (single member) ---
p = os.path.join(TAB, "v4ff_scores_by_model_block.csv")
if os.path.exists(p):
    rr = read_csv(p)
    add("v4 full-field ResU-Net", "Stage-4 TRUE spatial DL (55 gridded channels, Z500/T850/winds as fields); null result — no lift over box-index.",
        best_by_sedi(rr, "family"), full_rows(rr, "family"))

# --- §6 vote ensemble ---
p = os.path.join(TAB, "v6_vote_scores_by_block.csv")
if os.path.exists(p):
    rr = read_csv(p)
    prim = {block_of(r): {"model": "vote", "metrics": canon(r), "base_rate": r.get("base_rate")}
            for r in rr if block_of(r)}
    add("§6 vote ensemble", "23-member m/n vote, per-fold optimal P0 (EDI/ETS/HSS not in source table).",
        prim, [(prim[b]["model"], b, prim[b]["metrics"], prim[b]["base_rate"]) for b in prim])

# --- v3 stack (13-member) computed from pooled npz ---
STACK_KEY = None
p = os.path.join(ROOT, "outputs", "loyo_pred_v3stack", "v3stack_pooled.npz")
if os.path.exists(p):
    z = np.load(p, allow_pickle=True)
    y = z["y"].astype(float)
    bid = z["block_id"].astype(int)
    pstk = z["p_stack"].astype(float)
    prim = {}
    allr = []
    for k, b in enumerate(BLOCKS, start=1):
        m = (bid == k) & np.isfinite(y) & np.isfinite(pstk)  # finite labels only, as stage-25
        if not m.any():
            continue
        yy, pp = y[m].astype(int), pstk[m]
        thr, _ = vf.best_threshold(yy, pp, "Peirce")
        bs = vf.binary_scores(yy, pp, thr)
        met = {"POD": bs["POD"], "POFD": bs["POFD"], "FAR": bs["FAR"], "CSI": bs["CSI"],
               "ETS": bs["ETS"], "HSS": bs["HSS"], "Peirce": bs["Peirce"], "EDI": bs["EDI"],
               "SEDI": bs["SEDI"], "FBI": bs["freq_bias"]}
        br = f"{yy.mean():.4f}"
        prim[b] = {"model": "ens_stack", "metrics": met, "base_rate": br}
        allr.append(("ens_stack", b, met, br))
    add("v3 stack (13-member)",
        "Stage-3 error-diverse logistic stack over 13 members (incl. spatial ResU-Net/CNN); "
        "pooled Peirce threshold.", prim, allr)
    STACK_KEY = len(experiments) - 1

    # sanity vs existing v3_stack_scores_by_block.csv (SEDI): expect tiny delta because that
    # table scored NaN-label cells as negatives whereas we filter to finite labels (as stage-25).
    chk = os.path.join(TAB, "v3_stack_scores_by_block.csv")
    if os.path.exists(chk):
        ref = {(r["block"], r["combo"]): r for r in read_csv(chk)}
        msg = []
        for b in BLOCKS:
            r = ref.get((b, "ens_stack"))
            if r and b in prim:
                d = abs(float(r["SEDI"]) - prim[b]["metrics"]["SEDI"])
                msg.append(f"{b} ΔSEDI={d:.4f}")
        print("[26] v3-stack SEDI vs stage-23 table (finite-label filtering):", "  ".join(msg))

# ---------------------------------------------------------------- HTML helpers
def fmt(v):
    return "—" if v is None else f"{v:.3f}"


def cell_color(metric, v):
    """Green (good) -> red (bad) background for a metric value."""
    if v is None:
        return "#f2f2f2"
    _, direction, ideal = METRICS[metric]
    if direction == 0:            # FBI: distance from 1
        q = max(0.0, 1.0 - min(abs(v - ideal), 1.0))
    else:
        lo, hi = (0.0, 1.0)
        q = (v - lo) / (hi - lo)
        q = max(0.0, min(1.0, q))
        if direction < 0:
            q = 1.0 - q
    # q in [0,1], 1=good. interpolate red(0)->yellow(.5)->green(1)
    if q < 0.5:
        r, g = 220, int(80 + 300 * q)
    else:
        r, g = int(220 - 300 * (q - 0.5)), 200
    r = max(0, min(230, r)); g = max(0, min(210, g))
    return f"rgb({r},{g},90)"


def esc(s):
    return html.escape(str(s))


# ---------------------------------------------------------------- SVG grouped bar chart
PALETTE = ["#1f77b4", "#ff7f0e", "#2ca02c", "#d62728", "#9467bd", "#8c564b", "#17becf"]


def svg_chart(metric):
    label, direction, ideal = METRICS[metric]
    exps = experiments
    n_exp = len(exps)
    W, H = 760, 300
    padL, padR, padT, padB = 46, 12, 34, 40
    plotW, plotH = W - padL - padR, H - padT - padB
    # scale: use [0, max] for +metrics/FBI, [0, max] for -metrics too (raw values)
    vals = [e["rows"].get(b, {}).get("metrics", {}).get(metric)
            for e in exps for b in BLOCKS]
    vals = [v for v in vals if v is not None]
    vmax = max(vals + [ideal if direction == 0 else 0.0]) if vals else 1.0
    vmax = max(vmax, 0.05) * 1.12
    gw = plotW / len(BLOCKS)              # group width per block
    bw = gw * 0.8 / max(n_exp, 1)         # bar width
    s = [f'<svg viewBox="0 0 {W} {H}" width="100%" style="max-width:820px">']
    s.append(f'<text x="{padL}" y="18" font-size="13" font-weight="bold">{label} by lead block '
             f'({"lower=better" if direction<0 else "ideal=1" if direction==0 else "higher=better"})</text>')
    # y gridlines
    for gy in range(5):
        val = vmax * gy / 4
        yy = padT + plotH - plotH * gy / 4
        s.append(f'<line x1="{padL}" y1="{yy:.1f}" x2="{W-padR}" y2="{yy:.1f}" stroke="#eee"/>')
        s.append(f'<text x="{padL-4}" y="{yy+3:.1f}" font-size="9" text-anchor="end" fill="#666">{val:.2f}</text>')
    if direction == 0:  # FBI ideal=1 line
        yy = padT + plotH - plotH * min(ideal, vmax) / vmax
        s.append(f'<line x1="{padL}" y1="{yy:.1f}" x2="{W-padR}" y2="{yy:.1f}" stroke="#999" stroke-dasharray="4,3"/>')
    for bi, b in enumerate(BLOCKS):
        gx = padL + bi * gw
        s.append(f'<text x="{gx+gw/2:.1f}" y="{H-padB+16}" font-size="11" text-anchor="middle">{b}</text>')
        for ei, e in enumerate(exps):
            v = e["rows"].get(b, {}).get("metrics", {}).get(metric)
            if v is None:
                continue
            bh = plotH * min(v, vmax) / vmax
            x = gx + gw * 0.1 + ei * bw
            yy = padT + plotH - bh
            s.append(f'<rect x="{x:.1f}" y="{yy:.1f}" width="{bw*0.9:.1f}" height="{bh:.1f}" '
                     f'fill="{PALETTE[ei % len(PALETTE)]}"><title>{esc(e["label"])} {b} {label}={v:.3f}</title></rect>')
    s.append('</svg>')
    return "".join(s)


def legend():
    s = ['<div class="legend">']
    for ei, e in enumerate(experiments):
        s.append(f'<span class="lg"><i style="background:{PALETTE[ei%len(PALETTE)]}"></i>{esc(e["label"])}</span>')
    s.append('</div>')
    return "".join(s)


# ---------------------------------------------------------------- champion table
def champion_rows():
    """Per block: which experiment wins each focus metric, + the winning value."""
    out = {}
    for b in BLOCKS:
        out[b] = {}
        for m in FOCUS:
            best = None
            for e in experiments:
                v = e["rows"].get(b, {}).get("metrics", {}).get(m)
                if v is None:
                    continue
                _, direction, ideal = METRICS[m]
                key = -abs(v - ideal) if direction == 0 else (v if direction > 0 else -v)
                if best is None or key > best[0]:
                    best = (key, v, e["label"])
            out[b][m] = best
    return out


def build_html():
    css = """
    body{font-family:-apple-system,Segoe UI,Roboto,Helvetica,Arial,sans-serif;margin:0;color:#1a1a1a;background:#fafafa}
    .wrap{max-width:1180px;margin:0 auto;padding:22px 26px}
    h1{font-size:23px;margin:0 0 4px} h2{font-size:18px;margin:30px 0 8px;border-bottom:2px solid #ddd;padding-bottom:4px}
    h3{font-size:14px;margin:18px 0 6px;color:#333}
    p.sub{color:#555;margin:2px 0 14px;font-size:13px}
    table{border-collapse:collapse;font-size:12px;margin:6px 0 14px;width:100%}
    th,td{border:1px solid #d8d8d8;padding:4px 7px;text-align:center}
    th{background:#f0f0f0;position:sticky;top:0} td.l,th.l{text-align:left}
    td.metric{color:#111;font-variant-numeric:tabular-nums}
    .note{color:#666;font-size:11px;font-style:italic;margin:4px 0 18px}
    details{margin:6px 0} summary{cursor:pointer;font-weight:600;font-size:13px;padding:4px 0}
    .legend{margin:4px 0 10px;font-size:12px} .lg{display:inline-block;margin-right:14px}
    .lg i{display:inline-block;width:12px;height:12px;border-radius:2px;margin-right:5px;vertical-align:middle}
    .grid{display:flex;flex-wrap:wrap;gap:10px} .grid>div{flex:1 1 380px}
    .win{font-size:10px;color:#0a7} .best{font-weight:700}
    code{background:#eee;padding:1px 4px;border-radius:3px;font-size:11px}
    """
    H = ['<!doctype html><html><head><meta charset="utf-8">',
         '<title>yEHDE binary-forecast skill — experiment comparison</title>',
         f'<style>{css}</style></head><body><div class="wrap">']
    H.append('<h1>Cross-experiment comparison — yEHDE binary-forecast skill</h1>')
    H.append('<p class="sub">Primary target <b>yEHDE</b> (≥1 extreme-heat day in a lead block). '
             'All scores are leave-one-year-out pooled (2003–2022), out-of-sample, at each '
             'experiment\'s leakage-safe per-region×block decision threshold. Focus metrics: '
             'POD, POFD, FAR, FBI (frequency bias), CSI, ETS, HSS, Peirce, EDI, SEDI.</p>')

    # legend
    H.append('<h2>1. Experiments compared</h2>')
    H.append('<table><tr><th class="l">Experiment</th><th class="l">Representative model shown</th></tr>')
    for e in experiments:
        H.append(f'<tr><td class="l"><b>{esc(e["label"])}</b></td><td class="l">{esc(e["note"])}</td></tr>')
    H.append('</table>')

    # champion table
    H.append('<h2>2. Best experiment per lead block, per focus metric</h2>')
    H.append('<p class="sub">Winner among all experiments for each metric×block (FBI scored by closeness to 1).</p>')
    ch = champion_rows()
    H.append('<table><tr><th class="l">Metric</th>' + "".join(f'<th>{b}</th>' for b in BLOCKS) + '</tr>')
    for m in FOCUS:
        H.append(f'<tr><td class="l"><b>{METRICS[m][0]}</b></td>')
        for b in BLOCKS:
            best = ch[b].get(m)
            if best is None:
                H.append('<td>—</td>')
            else:
                _, v, lab = best
                H.append(f'<td class="metric" style="background:{cell_color(m, v)}">'
                         f'<span class="best">{v:.3f}</span><br><span class="win">{esc(lab)}</span></td>')
        H.append('</tr>')
    H.append('</table>')

    # charts
    H.append('<h2>3. Visual comparison (representative model per experiment)</h2>')
    H.append(legend())
    H.append('<div class="grid">')
    for m in FOCUS:
        H.append('<div>' + svg_chart(m) + '</div>')
    H.append('</div>')

    # side-by-side focus tables per block
    H.append('<h2>4. Focus-metric tables, all experiments × block</h2>')
    for b in BLOCKS:
        H.append(f'<h3>{b}</h3>')
        H.append('<table><tr><th class="l">Experiment</th><th class="l">model</th><th>base rate</th>'
                 + "".join(f'<th>{METRICS[m][0]}</th>' for m in ALLCOL) + '</tr>')
        for e in experiments:
            row = e["rows"].get(b)
            if not row:
                continue
            br = row.get("base_rate")
            try:
                brs = f"{float(br):.3f}"
            except (TypeError, ValueError):
                brs = "—"
            H.append(f'<tr><td class="l">{esc(e["label"])}</td><td class="l"><code>{esc(row["model"])}</code></td><td>{brs}</td>')
            for m in ALLCOL:
                v = row["metrics"].get(m)
                H.append(f'<td class="metric" style="background:{cell_color(m, v)}">{fmt(v)}</td>')
            H.append('</tr>')
        H.append('</table>')

    # full drilldown
    H.append('<h2>5. Full drill-down — every model in every experiment</h2>')
    for e in experiments:
        H.append(f'<details><summary>{esc(e["label"])} — {len(e["all"])} rows</summary>')
        H.append('<table><tr><th class="l">model</th><th>block</th><th>base rate</th>'
                 + "".join(f'<th>{METRICS[m][0]}</th>' for m in ALLCOL) + '</tr>')
        for model, b, met, br in sorted(e["all"], key=lambda t: (t[1], t[0])):
            try:
                brs = f"{float(br):.3f}"
            except (TypeError, ValueError):
                brs = "—"
            H.append(f'<tr><td class="l"><code>{esc(model)}</code></td><td>{b}</td><td>{brs}</td>')
            for m in ALLCOL:
                v = met.get(m)
                H.append(f'<td class="metric" style="background:{cell_color(m, v)}">{fmt(v)}</td>')
            H.append('</tr>')
        H.append('</table></details>')

    H.append('<h2>6. Notes on comparability &amp; provenance</h2>')
    H.append('<ul style="font-size:12.5px;color:#333;line-height:1.5">'
             '<li>All numbers are <b>out-of-sample</b> (20-fold leave-one-year-out, pooled 2003–2022), '
             'honouring the project leakage rule: no climatology, calibration, or decision threshold '
             'ever sees its own held-out year.</li>'
             '<li>Per-fold-thresholded experiments (v1, v2, Stage-2, v3 zoo, §6 vote) accumulate '
             'contingency counts using <b>each fold\'s own</b> per-region×block threshold. '
             'The <b>v3 stack</b> row uses a single <b>pooled Peirce</b> threshold per block '
             '(reproducing stage-23\'s reported SEDI/Peirce); it is marginally optimistic versus a '
             'per-fold threshold and should be read as such.</li>'
             '<li>The §6 vote source table does not carry EDI/ETS/HSS; those cells show “—”.</li>'
             '<li>Sources: <code>cv_binary_scores_by_block.csv</code>, '
             '<code>v2_cv_scores_by_model_block.csv</code>, <code>v2_scores_with_ensemble.csv</code>, '
             '<code>v3cal_binary_by_block.csv</code>, <code>v3zoo_scores_by_model_block.csv</code>, '
             '<code>v6_vote_scores_by_block.csv</code>, '
             '<code>outputs/loyo_pred_v3stack/v3stack_pooled.npz</code>. '
             'Regenerate with <code>python scripts/26_compare_experiments_html.py</code>.</li>'
             '</ul>')
    H.append('</div></body></html>')
    return "".join(H)


if __name__ == "__main__":
    if not experiments:
        sys.exit("[26] no experiment tables found under results/tables/")
    with open(OUT, "w") as fh:
        fh.write(build_html())
    print(f"[26] wrote {os.path.relpath(OUT, ROOT)}  ({len(experiments)} experiments)")