#!/usr/bin/env python3
"""One self-contained evaluation page, in the style of ../HW/ehd_ai_results.html.

Reads whatever of these exist and renders a single file with figures inlined as
base64 PNG, so it can be opened or emailed with nothing else attached:
  outputs/AI_Results/bcsd_fluxonly.json          downscaling ladder D0/D1/D2
  outputs/AI_Results/exp_features_sub.json       upper-air and model family
  outputs/AI_Results/threshold{,_region}.json    out-of-sample threshold rules
  outputs/AI_Results/threshold*.w?.curve.npz     the full delta sweeps
  outputs/AI_Results/w1_bootstrap_d*.json        year-block bootstrap
"""

from __future__ import annotations

import base64
import io
import json
import os
import sys
import time
from html import escape

import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import cfs_common as C

AR = os.path.join(C.OUT, "AI_Results")
OUT = "/home/tanpv/Progs/claude/HW_CFS_raw/ehd_cfs_evaluation.html"
WK = C.WEEK_NAMES
TARGETS = {"POD": (">", 0.60), "POFD": ("<", 0.20), "FAR": ("<", 0.20),
           "CSI": (">", 0.40), "FBI": ("~", 1.0)}

CSS = """
 body{font-family:-apple-system,Segoe UI,Roboto,Helvetica,Arial,sans-serif;
      margin:0 auto;max-width:1240px;padding:28px 22px 70px;color:#1a1a1a;line-height:1.55}
 h1{font-size:23px;margin:0 0 6px} h2{font-size:17px;margin:32px 0 8px}
 h3{font-size:14.5px;margin:20px 0 6px;color:#333}
 p{font-size:13.5px;margin:8px 0} p.sub{color:#555;margin:0 0 16px}
 table{border-collapse:collapse;font-size:12.5px;margin:8px 0 18px}
 caption{text-align:left;font-weight:600;font-size:13px;padding:6px 0;color:#333}
 th,td{border:1px solid #ddd;padding:3px 9px;text-align:right}
 th{background:#f6f6f6;font-weight:600} td:first-child,th:first-child{text-align:left}
 td.good{background:#dff3e0;font-weight:600} td.bad{background:#fdecea}
 td.warn{background:#fdf4e0;font-weight:600}
 figure{margin:0 0 8px} img{width:100%;display:block}
 code{background:#f4f4f4;padding:1px 4px;border-radius:3px;font-size:12px}
 .note{background:#f7f7f9;border-left:3px solid #bbb;padding:9px 14px;margin:14px 0;
       font-size:13px}
 .warnbox{background:#fff8e6;border-left:3px solid #e0a800;padding:9px 14px;
          margin:14px 0;font-size:13px}
 .stamp{font-size:12.5px;color:#555;background:#eef4fb;border:1px solid #cfe0f0;
        border-radius:3px;padding:8px 12px;margin:10px 0 0}
 .toc{font-size:13px;margin:16px 0 0;padding:10px 14px;background:#fafafa;
      border:1px solid #e3e3e3;border-radius:3px}
 .toc a{color:#1a5f9c;text-decoration:none;margin-right:14px;white-space:nowrap}
 .toc a:hover{text-decoration:underline}
"""


def cls(key, v, judge=True):
    if not judge or key not in TARGETS or not np.isfinite(v):
        return ""
    op, ref = TARGETS[key]
    ok = v > ref if op == ">" else v < ref if op == "<" else abs(v - ref) <= 0.25
    return " class='good'" if ok else " class='bad'"


def fig(fn):
    """Render a matplotlib figure to an inline <img>."""
    buf = io.BytesIO()
    fn()
    plt.savefig(buf, format="png", dpi=130, bbox_inches="tight")
    plt.close()
    return ("<figure><img src='data:image/png;base64,%s'></figure>"
            % base64.b64encode(buf.getvalue()).decode())


def jload(p):
    return json.load(open(p)) if os.path.exists(p) else None


def year_table(delta):
    """Per-year validation scores, computed rather than transcribed: a hand-typed
    table silently goes stale the first time the pipeline is re-run."""
    from threshold import contingency, regional
    f = os.path.join(AR, "bcsd_fluxonly.json.w0.fields.npz")
    if not os.path.exists(f):
        return ""
    z = np.load(f)
    score = z["d1"] - z["thr"]
    E, yrs, reg = z["E"].astype(bool), z["yrs"], z["reg"]
    valid = np.isin(yrs, C.VALID_YEARS) & np.isfinite(score).any(1)
    fc = np.zeros_like(E)
    ok = np.isfinite(score)
    fc[ok] = score[ok] >= delta
    FR, OR = regional(fc, reg), regional(E, reg)
    ks = ("POD", "POFD", "FAR", "FBI", "CSI")
    out = ["<table><caption>Từng năm kiểm định tại δ = %+.2f</caption>"
           "<tr><th>Năm</th><th>n init</th>%s</tr>"
           % (delta, "".join("<th>%s</th>" % k for k in ks))]
    for y in list(C.VALID_YEARS) + ["gộp"]:
        m = valid if y == "gộp" else (valid & (yrs == y))
        r = contingency(FR[m], OR[m])
        b = ("<b>%s</b>" % y) if y == "gộp" else str(y)
        out.append("<tr><td>%s</td><td>%d</td>%s</tr>"
                   % (b, m.sum(),
                      "".join("<td%s>%.3f</td>" % (cls(k, r[k]), r[k]) for k in ks)))
    out.append("</table>")
    return "".join(out)


# ------------------------------------------------------------------ sections
def sec_ladder(H):
    d = jload(os.path.join(AR, "bcsd_fluxonly.json"))
    if not d:
        return
    H.append("<h2 id='ladder'>1. Thang downscaling — D0 thô / D1 BCSD / D2 hồi quy ML</h2>")
    H.append("<p class='sub'>Toàn bộ 13 năm flux, 3086 init, kiểm định 2001/2006/2009. "
             "Cắt ngưỡng tại chính ngưỡng EHD địa phương (chưa dò δ).</p>")
    keys = ["ME", "MAE", "RMSE", "CORR", "exc_ratio", "POD", "FAR", "CSI", "SEDI"]
    H.append("<table><caption>Chấm hai lần: theo giá trị (°C) và theo sự kiện</caption>")
    H.append("<tr><th>Tuần</th><th>Bậc</th>" +
             "".join("<th>%s</th>" % k for k in keys) + "</tr>")
    for w in WK:
        for tag in ("D0", "D1", "D2"):
            r = d[w][tag]
            cells = "".join("<td%s>%.3f</td>" % (cls(k, r[k]), r[k]) for k in keys)
            H.append("<tr><td>%s</td><td>%s</td>%s</tr>" % (w, tag, cells))
    H.append("</table>")
    H.append("<div class='note'><b>Điều quan trọng nhất ở bảng này là cột "
             "<code>exc_ratio</code></b> (tỷ số tần suất vượt ngưỡng mô hình/quan trắc). "
             "D2 có MAE ngang hoặc tốt hơn D1 nhưng exc_ratio sụp còn 0.55–0.78 ở mọi "
             "tuần: hồi quy co trường về khí hậu học, trường bị co rồi cắt ngưỡng Tx90p "
             "thì gần như không sinh sự kiện, nên CSI thua D1 khắp nơi. "
             "RMSE hoàn toàn không nhìn thấy điều này. "
             "<b>Kết luận: chọn BCSD (D1). Muốn dùng ML thì phải là bộ phân loại "
             "P(EHD), không phải hồi quy giá trị rồi cắt ngưỡng.</b></div>")

    def f():
        x = np.arange(len(WK))
        fg, ax = plt.subplots(1, 2, figsize=(11, 3.4))
        for tag, c in (("D0", "#999"), ("D1", "#1a5f9c"), ("D2", "#c0504d")):
            ax[0].plot(x, [d[w][tag]["MAE"] for w in WK], "o-", color=c, label=tag)
            ax[1].plot(x, [d[w][tag]["exc_ratio"] for w in WK], "o-", color=c, label=tag)
        ax[0].set_ylabel("MAE (°C)"); ax[1].set_ylabel("exc_ratio")
        ax[1].axhline(1.0, ls="--", lw=1, color="k")
        for a in ax:
            a.set_xticks(x); a.set_xticklabels(WK); a.grid(alpha=.3); a.legend(fontsize=8)
        ax[0].set_title("Sai số giá trị", fontsize=10)
        ax[1].set_title("Tần suất sự kiện — D2 sụp ở đây", fontsize=10)
    H.append(fig(f))


def sec_features(H):
    d = jload(os.path.join(AR, "exp_features_sub.json"))
    if not d:
        return
    H.append("<h2 id='feat'>2. Biến cao không và mô hình phi tuyến mua được gì</h2>")
    H.append("<div class='warnbox'><b>SƠ BỘ.</b> Tập pgb hiện chỉ có 1999–2004 trên đĩa "
             "local, nên năm kiểm định rút còn 2001 tháng 1–8 (phủ trọn mùa nóng, thiếu "
             "tháng 9–12 tức phần “ấm bất thường mùa đông”). Sẽ chạy lại trên đủ 13 năm "
             "khi việc trích pgb qua NAS hoàn tất.</div>")
    H.append("<table><caption>AUC trên năm kiểm định. L = hồi quy logistic, "
             "G = LightGBM; flux = 16 biến bề mặt, both = thêm 27 biến cao không</caption>")
    H.append("<tr><th>Tuần</th><th>L_flux</th><th>L_both</th><th>G_flux</th>"
             "<th>G_both</th><th>AUC cần</th><th>Khoảng cách</th>"
             "<th>cao không mua</th><th>phi tuyến mua</th></tr>")
    for w in WK:
        r = d[w]
        best = max(r["L_flux"], r["L_both"], r["G_flux"], r["G_both"])
        gap = best - r["AUC_req"]
        H.append("<tr><td>%s</td>%s<td>%.4f</td><td class='bad'>%+.4f</td>"
                 "<td>%+.4f</td><td>%+.4f</td></tr>"
                 % (w, "".join("<td>%.4f</td>" % r[k] for k in
                               ("L_flux", "L_both", "G_flux", "G_both")),
                    r["AUC_req"], gap, r["L_both"] - r["L_flux"],
                    r["G_flux"] - r["L_flux"]))
    H.append("</table>")
    H.append("<div class='note'>Tốt nhất ở W1 là 0.787 so với 0.858 cần có — còn thiếu "
             "0.070. Biến cao không mua được +0.016 ở W1 và ~0 ở nơi khác; phi tuyến âm "
             "ở W1–W2. W4–W5 có vẻ được +0.03…+0.05 nhưng đó là dao động quanh mức ngẫu "
             "nhiên (AUC 0.51–0.56), không có ý nghĩa nghiệp vụ. "
             "<b>Đây là giới hạn THÔNG TIN, không phải giới hạn phương pháp.</b></div>")


def _curve(scale, w):
    p = os.path.join(AR, "threshold%s.w%d.curve.npz"
                     % ("_region" if scale == "region" else "", w))
    return np.load(p) if os.path.exists(p) else None


def sec_threshold(H):
    H.append("<h2 id='thr'>3. Chốt ngưỡng ngoài mẫu</h2>")
    H.append("<p class='sub'>Điểm số = <code>D1 − ngưỡng EHD địa phương</code>; phát báo "
             "khi điểm số ≥ δ. δ là tham số tự do duy nhất, <b>chốt trên năm huấn luyện "
             "rồi áp nguyên si sang năm kiểm định</b>. Quy mô vùng theo EXP1: một vùng có "
             "EHDE = 1 khi hơn 20% số ô của nó có EHD = 1.</p>")
    for scale, name in (("region", "Quy mô VÙNG (EXP1)"), ("point", "Quy mô ĐIỂM")):
        d = jload(os.path.join(AR, "threshold%s.json"
                               % ("_region" if scale == "region" else "")))
        if not d:
            continue
        H.append("<h3>%s</h3>" % name)
        H.append("<table><caption>Bốn quy tắc chốt δ. Hàng VA là kết quả "
                 "ngoài mẫu — đó là con số duy nhất đáng tin.</caption>")
        H.append("<tr><th>Tuần</th><th>Quy tắc</th><th>δ</th><th>Tập</th>"
                 "<th>POD</th><th>POFD</th><th>FAR</th><th>FBI</th><th>CSI</th>"
                 "<th>EDI</th><th>SEDI</th><th>nền</th></tr>")
        for w in WK:
            for rule in ("csi", "fbi", "far", "box"):
                r = d[w].get(rule)
                if r is None:
                    H.append("<tr><td>%s</td><td>%s</td><td colspan='10' "
                             "style='text-align:left'>không có ngưỡng nào đạt hộp chỉ "
                             "tiêu trên năm huấn luyện</td></tr>" % (w, rule))
                    continue
                for tag, sub in (("train", "huấn luyện"), ("valid", "KIỂM ĐỊNH")):
                    m = r[tag]
                    judge = tag == "valid"
                    cells = "".join("<td%s>%.3f</td>" % (cls(k, m[k], judge), m[k])
                                    for k in ("POD", "POFD", "FAR", "FBI", "CSI",
                                              "EDI", "SEDI"))
                    H.append("<tr><td>%s</td><td>%s</td><td>%+.2f</td><td>%s</td>%s"
                             "<td>%.4f</td></tr>"
                             % (w, rule, r["delta"], sub, cells, m["base"]))
        H.append("</table>")

    def f():
        fg, ax = plt.subplots(1, 2, figsize=(11, 3.8))
        for k, (scale, ttl) in enumerate((("region", "Vùng (EXP1)"), ("point", "Điểm"))):
            for w in range(6):
                z = _curve(scale, w)
                if z is None:
                    continue
                ax[k].plot(z["va_POFD"], z["va_POD"], "-", lw=1.6, label=WK[w])
            ax[k].add_patch(plt.Rectangle((0, .6), .2, .4, fc="#dff3e0", ec="#5a5",
                                          lw=1, zorder=0))
            ax[k].plot([0, 1], [0, 1], "k--", lw=.8)
            ax[k].set_xlabel("POFD"); ax[k].set_ylabel("POD")
            ax[k].set_title("%s — ROC trên năm kiểm định" % ttl, fontsize=10)
            ax[k].grid(alpha=.3); ax[k].legend(fontsize=8, loc="lower right")
            ax[k].set_xlim(0, 1); ax[k].set_ylim(0, 1)
    H.append(fig(f))
    H.append("<p class='sub'>Ô xanh là vùng đạt POD &gt; 0.6 và POFD &lt; 0.2. "
             "Chỉ đường W1 quy mô vùng đi qua nó.</p>")


def sec_operating(H):
    z = _curve("region", 0)
    b30 = jload(os.path.join(AR, "w1_bootstrap.json"))
    b40 = jload(os.path.join(AR, "w1_bootstrap_d0.40.json"))
    if z is None:
        return
    H.append("<h2 id='op'>4. Điểm làm việc đề nghị cho W1 theo vùng</h2>")

    def f():
        d = z["delta"]
        m = (d >= -0.5) & (d <= 1.5)
        fg, ax = plt.subplots(figsize=(7.5, 4))
        for k, c in (("POD", "#1a5f9c"), ("POFD", "#c0504d"),
                     ("FAR", "#e08214"), ("CSI", "#4a7c59")):
            ax.plot(d[m], z["tr_" + k][m], "--", lw=1, color=c, alpha=.55)
            ax.plot(d[m], z["va_" + k][m], "-", lw=1.9, color=c, label=k)
        ax.axhline(.6, color="#1a5f9c", lw=.7, ls=":")
        ax.axhline(.2, color="#c0504d", lw=.7, ls=":")
        ax.axvline(.40, color="k", lw=1.2)
        ax.annotate("δ = +0.40", (.40, .93), fontsize=9, ha="left",
                    xytext=(.46, .93))
        ax.set_xlabel("δ (°C trên ngưỡng EHD)"); ax.set_ylabel("giá trị")
        ax.set_title("W1 vùng — nét liền: kiểm định, nét đứt: huấn luyện", fontsize=10)
        ax.grid(alpha=.3); ax.legend(fontsize=9, ncol=4)
    H.append(fig(f))

    if b40:
        H.append("<table><caption>Bootstrap khối NĂM, 5000 lần, ngưỡng giữ cố định "
                 "(khoảng tin cậy 90%)</caption>")
        H.append("<tr><th>Chỉ số</th><th>δ=+0.30 (quy tắc box)</th>"
                 "<th>δ=+0.40 (có biên an toàn)</th><th>Chỉ tiêu</th></tr>")
        tgt = {"POD": "&gt; 0.60", "POFD": "&lt; 0.20", "FAR": "&lt; 0.20",
               "FBI": "~ 1", "CSI": "&gt; 0.40", "EDI": "—", "SEDI": "—"}
        for k in ("POD", "POFD", "FAR", "FBI", "CSI", "EDI", "SEDI"):
            a, b = b30[k], b40[k]
            H.append("<tr><td>%s</td><td>%.3f [%.3f–%.3f]</td>"
                     "<td%s>%.3f [%.3f–%.3f]</td><td>%s</td></tr>"
                     % (k, a["med"], a["lo"], a["hi"], cls(k, b["med"]),
                        b["med"], b["lo"], b["hi"], tgt[k]))
        H.append("<tr><td>đạt <b>trọn</b> hộp chỉ tiêu</td><td>%.1f%% số lần</td>"
                 "<td><b>%.1f%% số lần</b></td><td>—</td></tr>"
                 % (100 * b30["box_rate"], 100 * b40["box_rate"]))
        H.append("</table>")

    H.append(year_table(0.40))
    H.append("<div class='warnbox'><b>Hai điều phải cân nhắc trước khi dùng.</b><br>"
             "<b>(1) “FBI ≈ 1” và “POFD &lt; 0.2” mâu thuẫn nhau trên đường cong này.</b> "
             "Tại FBI = 1 (δ = 0) thì POFD là 0.328, gấp 1.8 lần chỉ tiêu. Ngược lại điểm "
             "đạt chỉ tiêu có FBI 0.75, tức phát báo thiếu 25%. Tôi ưu tiên POFD vì nó là "
             "con số định lượng trong bộ chỉ tiêu, còn “FBI ~ 1” nêu định tính — "
             "<b>nhưng đây là lựa chọn của người dùng, không phải của tôi.</b><br>"
             "<b>(2) 2006 chi phối đuôi trên của POFD.</b> Năm này chỉ có 88 init (năm "
             "flux thưa) và POFD 0.330; một mình nó đẩy bách phân vị 95 của POFD lên "
             "0.238. Ba năm kiểm định là quá ít để bootstrap khối năm cho kết luận chắc — "
             "đó là giới hạn của cách chia mẫu, bootstrap không sửa được.</div>")


def sec_concl(H):
    H.append("<h2 id='kl'>5. Kết luận</h2>")
    H.append("<table><caption>Trạng thái so với bộ chỉ tiêu</caption>"
             "<tr><th>Sản phẩm</th><th>Kết luận</th></tr>"
             "<tr><td>W1, quy mô vùng (EXP1)</td><td class='good'>ĐẠT tại δ = +0.40: "
             "POD 0.653, POFD 0.185, FAR 0.133, CSI 0.593 — ngoài mẫu</td></tr>"
             "<tr><td>W1, quy mô điểm</td><td class='bad'>KHÔNG đạt ở bất kỳ ngưỡng nào</td></tr>"
             "<tr><td>W2–W6, cả hai quy mô</td><td class='bad'>KHÔNG đạt; không tồn tại "
             "δ nào cho FAR &lt; 0.2 ngay cả trên năm huấn luyện</td></tr></table>")
    H.append("<div class='note'><b>Đề nghị.</b> W1 phát báo nhị phân theo vùng ở "
             "δ = +0.40. W2–W6 và toàn bộ sản phẩm theo điểm phát dưới dạng "
             "<b>trường xác suất đã hiệu chuẩn</b>, kèm ghi rõ trần tín hiệu — vì khoảng "
             "cách AUC 0.25–0.36 là giới hạn thông tin của CFS chứ không phải chuyện dò "
             "ngưỡng hay đổi lớp mô hình.<br><br>"
             "<b>Hướng còn lại chưa thử và có cơ sở:</b> dùng bộ <b>phân loại P(EHD)</b> "
             "rồi mới gộp vùng, thay vì gộp vùng từ trường tất định D1 — vì mục 1 đã chỉ "
             "ra hồi quy giá trị là sai lớp mô hình cho bài toán sự kiện.</div>")


def main():
    H = ["<!doctype html><html lang='vi'><meta charset='utf-8'>",
         "<title>Đánh giá dự báo EHD từ CFSv2</title><style>%s</style>" % CSS,
         "<h1>Dự báo S2S ngày nóng cực đoan từ CFSv2 — kết quả đánh giá</h1>",
         "<p class='sub'>Downscaling BCSD về lưới 0.25° và chốt ngưỡng ngoài mẫu, "
         "hạn 0–45 ngày.</p>",
         "<div class='stamp'>Tập mẫu <b>3086 init</b> (flux ∩ pgb), 13 năm, 541 ô lưới "
         "đất, 7 vùng. Kiểm định <b>2001 / 2006 / 2009</b> (726 init), huấn luyện 10 năm "
         "còn lại (2360 init). Mọi khí hậu học, quantile map và ngưỡng đều chỉ khớp trên "
         "năm huấn luyện. Tuần: W1 = hạn 0–10, W2 = 11–17, W3 = 18–24, W4 = 25–31, "
         "W5 = 32–38, W6 = 39–45. Dựng lúc %s.</div>"
         % time.strftime("%Y-%m-%d %H:%M"),
         "<div class='toc'><a href='#ladder'>1. Thang downscaling</a>"
         "<a href='#feat'>2. Biến cao không &amp; phi tuyến</a>"
         "<a href='#thr'>3. Chốt ngưỡng ngoài mẫu</a>"
         "<a href='#op'>4. Điểm làm việc W1</a>"
         "<a href='#kl'>5. Kết luận</a></div>"]
    sec_ladder(H)
    sec_features(H)
    sec_threshold(H)
    sec_operating(H)
    sec_concl(H)
    H.append("</html>")
    open(OUT, "w").write("\n".join(H))
    print("viet %s (%.1f MB)" % (OUT, os.path.getsize(OUT) / 1e6))


if __name__ == "__main__":
    main()