#!/usr/bin/env python3
"""Year-block bootstrap for the chosen W1 operating point (EXP1 regional scale).

The threshold comes from the training years and is held fixed; only the
validation sample is resampled. Blocks are whole YEARS, because inits a few days
apart share weather and resampling them independently would understate the
spread by a lot. Three validation years is a coarse block set -- that is a real
limit of the design, not something the bootstrap can repair -- so the per-year
numbers are printed alongside.
"""
from __future__ import annotations
import argparse, json, os, sys
import numpy as np
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import cfs_common as C
from threshold import contingency, regional, in_box

FLD = os.path.join(C.OUT, "AI_Results", "bcsd_fluxonly.json")
NB = 5000
rng = np.random.default_rng(20260903)

_ap = argparse.ArgumentParser()
_ap.add_argument("--delta", type=float, default=0.40)
_a = _ap.parse_args()
w, DELTA = 0, _a.delta                   # W1, threshold fitted on training years
z = np.load("%s.w%d.fields.npz" % (FLD, w))
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_all = np.zeros_like(E)
ok = np.isfinite(score)
fc_all[ok] = score[ok] >= DELTA
FR, OR = regional(fc_all, reg), regional(E, reg)

print("W1 EXP1, delta = %+.2f chot tren nam huan luyen\n" % DELTA)
print("%-8s %6s %6s %6s %6s %6s %6s" % ("nam", "n", "POD", "POFD", "FAR", "FBI", "CSI"))
for y in C.VALID_YEARS:
    m = valid & (yrs == y)
    r = contingency(FR[m], OR[m])
    print("%-8s %6d %6.3f %6.3f %6.3f %6.3f %6.3f"
          % (y, m.sum(), r["POD"], r["POFD"], r["FAR"], r["FBI"], r["CSI"]))
full = contingency(FR[valid], OR[valid])
print("%-8s %6d %6.3f %6.3f %6.3f %6.3f %6.3f"
      % ("gop", valid.sum(), full["POD"], full["POFD"], full["FAR"],
         full["FBI"], full["CSI"]))

yr_idx = [np.where(valid & (yrs == y))[0] for y in C.VALID_YEARS]
keys = ("POD", "POFD", "FAR", "FBI", "CSI", "EDI", "SEDI")
draws = {k: [] for k in keys}
nbox = 0
for _ in range(NB):
    take = np.concatenate([yr_idx[i] for i in rng.integers(0, len(yr_idx), len(yr_idx))])
    r = contingency(FR[take], OR[take])
    for k in keys:
        draws[k].append(r[k])
    nbox += in_box(r)

print("\nbootstrap khoi NAM, %d lan (khoang tin cay 90%%)" % NB)
print("%-6s %8s %8s %8s   %s" % ("", "5%", "trung vi", "95%", "chi tieu"))
tgt = {"POD": "> 0.60", "POFD": "< 0.20", "FAR": "< 0.20", "FBI": "~ 1",
       "CSI": "> 0.40", "EDI": "-", "SEDI": "-"}
res = {}
for k in keys:
    v = np.array(draws[k]); q = np.percentile(v, [5, 50, 95])
    res[k] = dict(lo=q[0], med=q[1], hi=q[2])
    print("%-6s %8.3f %8.3f %8.3f   %s" % (k, q[0], q[1], q[2], tgt[k]))
print("\nty le lan bootstrap dat TRON hop chi tieu: %.1f%%" % (100 * nbox / NB))
res["box_rate"] = nbox / NB
res["delta"] = DELTA
json.dump(res, open(os.path.join(C.OUT, "AI_Results", "w1_bootstrap_d%.2f.json" % DELTA), "w"),
          indent=1)