"""Score Met Office (UKMO) 70% ranges from data/ukmo_forecasts.csv against HURDAT2 counts for the
forecast period (June-November unless the period says otherwise; July-November in 2007-2008).
Also central-value MSSS vs the 1995..y-1 mean for the same period."""
import csv, re, numpy as np, collections
from scipy.stats import binomtest
storms = list(csv.DictReader(open('../data/storms.csv')))
ace_fix = list(csv.DictReader(open('../data/ace_fixes.csv')))
def window(period):
    p = (period or '').lower()
    if 'whole season' in p: return ('0201' if 'feb' in p else '0101'), '1232'  # Aug updates: observed + Aug-Jan forecast (next-Jan storms: none 2021-2025)
    if p.startswith('jul'): return '0701', '1201'
    if p.startswith('aug'): return '0801', '1201'
    return '0601', '1201'
def cnt(y, k, a, b):
    key = {'ns': 't_ns', 'h': 't_h', 'mh': 't_mh'}[k]
    return sum(1 for s in storms if s[key] and int(s[key][:4]) == y and a <= s[key][4:8] < b)
ACE = collections.defaultdict(float)
for r in ace_fix: ACE[(int(r['date'][:4]), r['date'][4:8])] += float(r['ace'])
def ace(y, a, b): return sum(v for (yy, md), v in ACE.items() if yy == y and a <= md < b)
rows = list(csv.DictReader(open('../data/ukmo_forecasts.csv')))
res = collections.defaultdict(list); sq = collections.defaultdict(list)
for r in rows:
    y = int(r['year'])
    if y > 2025: continue
    a, b = window(r['period'])
    for k in ('ns', 'h', 'mh', 'ace'):
        if not r.get(k + '_lo') or not r.get(k + '_hi'): continue
        v = ace(y, a, b) if k == 'ace' else cnt(y, k, a, b)
        lo, hi = float(r[k + '_lo']), float(r[k + '_hi'])
        res[k].append((y, r['issued'], lo <= v <= hi, 'above' if v > hi else ('below' if v < lo else 'in'), v, lo, hi))
        if r.get(k + '_c'):
            base = np.mean([ace(x, a, b) if k == 'ace' else cnt(x, k, a, b) for x in range(1995, y)])
            sq[k].append(((float(r[k + '_c']) - v) ** 2, (base - v) ** 2))
for k, xs in res.items():
    hit = sum(x[2] for x in xs); n = len(xs); c = collections.Counter(x[3] for x in xs)
    msss = 1 - np.mean([e for e, _ in sq[k]]) / np.mean([b for _, b in sq[k]]) if sq[k] else float('nan')
    print(f'UKMO {k:4s} in {hit}/{n} = {hit/n:.2f} (above {c["above"]}, below {c["below"]}), binom p vs .7 = {binomtest(hit, n, .7).pvalue:.3f}; central MSSS vs 1995..y-1 mean {msss:+.2f} (n={len(sq[k])})')
    for x in xs:
        if not x[2]: print('    miss', x[0], x[1], f'obs {x[4]:.0f} vs {x[5]:.0f}-{x[6]:.0f}')

# Interval score vs climatological 15-85% intervals over 1995..y-1 for the same window
def IS(l, u, y, a=0.3): return (u - l) + (2 / a) * max(l - y, 0) + (2 / a) * max(y - u, 0)
print('\nInterval score (alpha=0.3), UKMO vs climatology 1995..y-1, same season window:')
for k in ('ns', 'h', 'mh', 'ace'):
    pairs = []
    for r in rows:
        y = int(r['year'])
        if y > 2025 or not r.get(k + '_lo'): continue
        a, b = window(r['period'])
        f = (lambda x: ace(x, a, b)) if k == 'ace' else (lambda x: cnt(x, k, a, b))
        v = f(y); hist = [f(x) for x in range(1995, y)]
        l, u = np.quantile(hist, [.15, .85])
        if k != 'ace': l, u = np.floor(l), np.ceil(u)
        pairs.append((IS(float(r[k + '_lo']), float(r[k + '_hi']), v), IS(l, u, v), float(r[k + '_hi']) - float(r[k + '_lo']), u - l, l <= v <= u))
    P = np.array(pairs)
    print(f'  {k:4s} n={len(P)}: IS UKMO {P[:,0].mean():.1f} vs clim {P[:,1].mean():.1f} (ratio {P[:,0].mean()/P[:,1].mean():.2f}); width UKMO {P[:,2].mean():.1f} vs clim {P[:,3].mean():.1f}; clim coverage {P[:,4].mean():.2f}')