"""Score NOAA May/August Atlantic outlooks 2001-2025 against HURDAT2 (full-year counts).
Outputs: ../data/noaa_scored.csv and a printed summary."""
import csv, math, collections
from scipy.stats import binomtest
obs = {int(r['year']): r for r in csv.DictReader(open('../data/season_counts.csv'))}
MED = lambda y: 87.5 if y <= 2010 else (92.4 if y <= 2020 else 96.7)  # ACE median each outlook refers to
rows = []
for r in csv.DictReader(open('../data/noaa_outlooks.csv')):
    y = int(r['year']); o = obs[y]; d = dict(year=y, month=r['month'], stated=r['range_prob_stated'])
    for k in ('ns', 'h', 'mh'):
        lo, hi, v = int(r[k + '_lo']), int(r[k + '_hi']), int(o[k])
        d[k + '_obs'] = v; d[k + '_lo'] = lo; d[k + '_hi'] = hi
        d[k + '_res'] = 'in' if lo <= v <= hi else ('above' if v > hi else 'below')
    if r['ace_lo_pct']:
        pct = float(o['ace']) / MED(y) * 100
        d['ace_pct_obs'] = round(pct); d['ace_lo'] = r['ace_lo_pct']; d['ace_hi'] = r['ace_hi_pct']
        d['ace_res'] = 'in' if float(r['ace_lo_pct']) <= pct <= float(r['ace_hi_pct']) else ('above' if pct > float(r['ace_hi_pct']) else 'below')
    else:
        d['ace_pct_obs'] = d['ace_lo'] = d['ace_hi'] = d['ace_res'] = ''
    rows.append(d)
with open('../data/noaa_scored.csv', 'w', newline='') as fh:
    w = csv.DictWriter(fh, rows[0].keys()); w.writeheader(); w.writerows(rows)
def summ(sel, label):
    print(f'\n## {label} (n outlooks = {len(sel)})')
    allin = 0
    for k in ('ns', 'h', 'mh', 'ace'):
        rs = [d[k + '_res'] for d in sel if d[k + '_res']]
        c = collections.Counter(rs); n = len(rs); hit = c['in']
        p = binomtest(hit, n, 0.7).pvalue if n else float('nan')
        print(f'  {k:4s} in {hit:2d}/{n:2d} = {hit/n:.2f}   above {c["above"]:2d}  below {c["below"]:2d}   binom p vs 0.70 = {p:.3f}')
    a3 = sum(all(d[k + '_res'] == 'in' for k in ('ns', 'h', 'mh')) for d in sel)
    print(f'  all three counts in range: {a3}/{len(sel)}')
summ(rows, 'All NOAA outlooks 2001-2025')
summ([d for d in rows if d['year'] >= 2008], 'Stated-probability era, 2008-2025')
summ([d for d in rows if d['year'] >= 2008 and d['month'] == 'May'], 'May, 2008-2025')
summ([d for d in rows if d['year'] >= 2008 and d['month'] == 'Aug'], 'August, 2008-2025')
summ([d for d in rows if d['year'] < 2008], '2001-2007 (no stated probability)')
print('\nPer-year misses (2008+):')
for d in rows:
    if d['year'] >= 2008:
        m = [f"{k}:{d[k+'_obs']} vs {d[k+'_lo']}-{d[k+'_hi']}" for k in ('ns', 'h', 'mh') if d[k + '_res'] != 'in']
        if d['ace_res'] and d['ace_res'] != 'in': m.append(f"ace:{d['ace_pct_obs']}% vs {d['ace_lo']}-{d['ace_hi']}")
        print(' ', d['year'], d['month'], '; '.join(m) or '-')