"""Mean-squared-error skill of point forecasts (range midpoints for range forecasts) by issuer and
issue window, 2001-2025, against (a) trailing 1995..y-1 mean, (b) trailing 30-yr mean, and for
issues from July on (c) count-to-issue-date + trailing-30-yr mean of the remainder.
MSSS = 1 - MSE_f / MSE_ref over the same seasons. Issuer names normalised; UKMO rows that are
June-November only are scored against June-November counts."""
import csv, re, numpy as np, collections, datetime
obs = {int(r['year']): r for r in csv.DictReader(open('../data/season_counts.csv'))}
storms = list(csv.DictReader(open('../data/storms.csv')))
def cnt(y, k, after=None, before=None):
    key = {'ns': 't_ns', 'h': 't_h', 'mh': 't_mh'}[k]; n = 0
    for s in storms:
        t = s[key]
        if not t or int(t[:4]) != y: continue
        if after and t[4:8] < after: continue
        if before and t[4:8] >= before: continue
        n += 1
    return n
JUNENOV = ('TSR', 'UKMO')  # issuers whose forecasts cover June-November only (primary-source check, verify_results.jsonl)
NORM = {'TWC/WSI': 'TWC', 'WSI/TWC': 'TWC', 'WSI': 'TWC', 'UKMET': 'UKMO', 'UKMO*': 'UKMO',
        'FSU COAPS': 'FSU', 'INSMET': 'InsMet'}
def window(issued, y):
    d = datetime.date.fromisoformat(issued)
    if d.year < y or d.month <= 2: return '1 Dec-Feb'
    return {3: '2 Mar-Apr', 4: '2 Mar-Apr', 5: '3 May-Jun', 6: '3 May-Jun', 7: '4 Jul', 8: '5 Aug'}.get(d.month, '6 Sep+')
F = collections.defaultdict(dict)  # (issuer, window) -> year -> row (latest in window)
# NOAA: use the curated table from CPC pages (Wikipedia lacks Aug 2003/2004); issue dates from Wikipedia where present
wdates = {(r['year'], r['issued'][5:7]): r['issued'] for r in csv.DictReader(open('../data/forecasts.csv')) if r['source'].strip() == 'NOAA'}
curated = []
for r in csv.DictReader(open('../data/noaa_outlooks.csv')):
    mm = '05' if r['month'] == 'May' else '08'
    d = dict(year=r['year'], source='NOAA', issued=wdates.get((r['year'], mm), f"{r['year']}-{mm}-{'20' if mm == '05' else '08'}"), note='')
    for k in ('ns', 'h', 'mh'): d[k + '_pt'] = (int(r[k + '_lo']) + int(r[k + '_hi'])) / 2
    curated.append(d)
for r in [x for x in csv.DictReader(open('../data/forecasts.csv')) if x['source'].strip() != 'NOAA'] + curated:
    y = int(r['year'])
    if not 2001 <= y <= 2025: continue
    src = NORM.get(r['source'].strip(), r['source'].strip().replace('*', ''))
    w = window(r['issued'], y)
    prev = F[(src, w)].get(y)
    if prev is None or r['issued'] >= prev['issued']: F[(src, w)][y] = r
out = []
print(f'{"issuer":7s} {"window":10s} var  n   MSSS_vs_era1995  MSSS_vs_clim30  MSSS_vs_todate  bias')
for (src, w), yr in sorted(F.items(), key=lambda kv: (kv[0][1], kv[0][0])):
    for k in ('ns', 'h', 'mh'):
        e, e95, e30, etd, bias = [], [], [], [], []
        for y, r in yr.items():
            if r[k + '_pt'] == '': continue
            junenov = 'june-nov' in r['note'] or src in JUNENOV
            if junenov: v = cnt(y, k, after='0601', before='1201'); clim = lambda x: cnt(x, k, after='0601', before='1201')
            else: v = int(obs[y][k]); clim = lambda x: int(obs[x][k])
            f = float(r[k + '_pt'])
            e.append((f - v) ** 2); bias.append(f - v)
            e95.append((np.mean([clim(x) for x in range(1995, y)]) - v) ** 2)
            e30.append((np.mean([clim(x) for x in range(y - 30, y)]) - v) ** 2)
            if w >= '4' and not junenov:
                md = r['issued'][5:7] + r['issued'][8:10]
                etd.append((cnt(y, k, before=md) + np.mean([cnt(x, k, after=md) for x in range(y - 30, y)]) - v) ** 2)
        n = len(e)
        if n < 8: continue
        m = np.mean(e); s95 = 1 - m / np.mean(e95); s30 = 1 - m / np.mean(e30); std = (1 - m / np.mean(etd)) if etd else float('nan')
        out.append(dict(issuer=src, window=w[2:], var=k, n=n, msss_era1995=round(s95, 2), msss_clim30=round(s30, 2), msss_todate=round(std, 2), mean_bias=round(np.mean(bias), 2), rmse=round(m ** .5, 2)))
        print(f'{src:7s} {w[2:]:10s} {k:3s} {n:2d}   {s95:+.2f}            {s30:+.2f}           {std:+.2f}         {np.mean(bias):+.2f}')
with open('../data/point_skill.csv', 'w', newline='') as fh:
    wr = csv.DictWriter(fh, out[0].keys()); wr.writeheader(); wr.writerows(out)