"""Bootstrap CIs (resampling seasons) for MSSS of each issuer x window, pooled over NS/H/MH by
summing MSE normalised per variable, against several naive baselines:
  era1995 (mean 1995..y-1), trail15 (mean y-15..y-1), trail10, clim30 (y-30..y-1);
  for July/August issues also todate30 / todate95: storms formed before the issue date + mean
  count after that date over y-30..y-1 / 1995..y-1 (same season window as the issuer).
Writes ../data/skill_ci.csv."""
import csv, numpy as np, importlib.util, io, contextlib, collections
spec = importlib.util.spec_from_file_location('p', 'point_skill.py')
with contextlib.redirect_stdout(io.StringIO()):
    p = importlib.util.module_from_spec(spec); spec.loader.exec_module(p)
obs, F = p.obs, p.F
BASE = {'era1995': lambda y: range(1995, y), 'trail15': lambda y: range(y - 15, y),
        'trail10': lambda y: range(y - 10, y), 'clim30': lambda y: range(y - 30, y)}
rng = np.random.default_rng(1); out = []
for (src, w), yr in sorted(F.items(), key=lambda kv: (kv[0][1], kv[0][0])):
    if src == 'UKMO': continue  # scored separately with its 70% ranges (score_ukmo.py)
    ys = sorted(y for y, r in yr.items() if all(r[k + '_pt'] != '' for k in ('ns', 'h', 'mh')))
    if len(ys) < 10: continue
    jn = src in p.JUNENOV
    V = (lambda x, k: p.cnt(x, k, after='0601', before='1201')) if jn else (lambda x, k: int(obs[x][k]))
    E = {y: {k: (float(yr[y][k + '_pt']) - V(y, k)) ** 2 for k in ('ns', 'h', 'mh')} for y in ys}
    line = f'{src:5s} {w[2:]:8s} n={len(ys):2d}'
    rec = dict(issuer=src, window=w[2:], n=len(ys))
    for b, rg in BASE.items():
        R = {y: {k: (np.mean([V(x, k) for x in rg(y)]) - V(y, k)) ** 2 for k in ('ns', 'h', 'mh')} for y in ys}
        def msss(sample):
            num = sum(sum(E[y][k] for y in sample) / max(sum(R[y][k] for y in sample), 1e-9) for k in ('ns', 'h', 'mh'))
            return 1 - num / 3
        est = msss(ys); bs = [msss(list(rng.choice(ys, len(ys)))) for _ in range(2000)]
        lo, hi = np.quantile(bs, [.05, .95])
        rec[b] = round(est, 2); rec[b + '_lo90'] = round(lo, 2); rec[b + '_hi90'] = round(hi, 2)
        line += f' | {b} {est:+.2f} [{lo:+.2f},{hi:+.2f}]'
    if w[0] in '45':  # July/August issues: storms formed before the issue date + mean of the rest of season
        a0, b0 = ('0601', '1201') if jn else ('0101', '1232')
        for b, rg in (('todate30', lambda y: range(y - 30, y)), ('todate95', lambda y: range(1995, y))):
            def RB(y, k):
                md = yr[y]['issued'][5:7] + yr[y]['issued'][8:10]
                done = p.cnt(y, k, after=a0, before=md)
                rest = np.mean([p.cnt(x, k, after=md, before=b0) for x in rg(y)])
                return (done + rest - V(y, k)) ** 2
            R = {y: {k: RB(y, k) for k in ('ns', 'h', 'mh')} for y in ys}
            def msss(sample):
                num = sum(sum(E[y][k] for y in sample) / max(sum(R[y][k] for y in sample), 1e-9) for k in ('ns', 'h', 'mh'))
                return 1 - num / 3
            est = msss(ys); bs = [msss(list(rng.choice(ys, len(ys)))) for _ in range(2000)]
            lo, hi = np.quantile(bs, [.05, .95])
            rec[b] = round(est, 2); rec[b + '_lo90'] = round(lo, 2); rec[b + '_hi90'] = round(hi, 2)
            line += f' | {b} {est:+.2f} [{lo:+.2f},{hi:+.2f}]'
    out.append(rec); print(line)
with open('../data/skill_ci.csv', 'w', newline='') as fh:
    keys = list(dict.fromkeys(k for r in out for k in r)); wr = csv.DictWriter(fh, keys); wr.writeheader(); wr.writerows(out)