"""Score NOAA's above/near/below-normal probabilities (2001-2025) with the ranked probability
score (RPS) against baselines. Outcome classes:
  def2021: ACE > 126.1 above, < 73.0 below (CPC Background page, 1951-2020 basis)
  def2015: above = ACE > 111 and >=2 of (NS>=13, H>=7, MH>=3); below = ACE < 66 or (NS<=9, H<=4, MH<=1)
  at_issue: ACE-only thresholds in force at issue: 2001-2010 75%/117% of 87.5; 2011-2020 66/111;
            2021+ 73.0/126.1 (older count conditions ignored)
ACE from current HURDAT2 (tropical + subtropical)."""
import csv, numpy as np
obs = {int(r['year']): r for r in csv.DictReader(open('../data/season_counts.csv'))}
def cls(y, d):
    a = float(obs[y]['ace']); ns, h, mh = (int(obs[y][k]) for k in ('ns', 'h', 'mh'))
    if d == 'def2021': return 0 if a > 126.1 else (2 if a < 73.0 else 1)
    if d == 'at_issue':  # ACE-only thresholds of the definition in force when the outlook was issued
        lo, hi = (0.75 * 87.5, 1.17 * 87.5) if y <= 2010 else ((66.0, 111.0) if y <= 2020 else (73.0, 126.1))
        return 0 if a > hi else (2 if a < lo else 1)
    if a > 111 and (ns >= 13) + (h >= 7) + (mh >= 3) >= 2: return 0
    if a < 66 or (ns <= 9 and h <= 4 and mh <= 1): return 2
    return 1
def rps(p, o):
    P = np.cumsum(p); O = np.cumsum(np.eye(3)[o]); return float(np.sum((P - O) ** 2)) / 2
noaa = list(csv.DictReader(open('../data/noaa_outlooks.csv')))
for d in ('def2021', 'def2015', 'at_issue'):
    print(f'\n== outcome definition {d}')
    for era, lo, hi in (('2001-2025', 2001, 2025), ('2008-2025', 2008, 2025), ('2001-2007', 2001, 2007)):
        for m in ('May', 'Aug'):
            S = {'NOAA': [], 'equal': [], 'clim1951': [], 'clim1995': []}; top_hit = []
            for r in noaa:
                y = int(r['year'])
                if r['month'] != m or not lo <= y <= hi: continue
                o = cls(y, d); p = np.array([float(r['p_above']), float(r['p_near']), float(r['p_below'])]) / 100
                S['NOAA'].append(rps(p, o)); S['equal'].append(rps(np.ones(3) / 3, o))
                for name, start in (('clim1951', 1951), ('clim1995', 1995)):
                    c = np.bincount([cls(x, d) for x in range(start, y)], minlength=3) + 0.5
                    S[name].append(rps(c / c.sum(), o))
                top_hit.append(int(np.argmax(p) == o))
            mean = {k: np.mean(v) for k, v in S.items()}
            line = f'  {era} {m} n={len(S["NOAA"])}: RPS NOAA {mean["NOAA"]:.3f}'
            for b in ('equal', 'clim1951', 'clim1995'):
                line += f' | {b} {mean[b]:.3f} (RPSS {1 - mean["NOAA"] / mean[b]:+.2f})'
            print(line + f' | most-likely class right {sum(top_hit)}/{len(top_hit)}')
# reliability of p(above), by period and outcome definition
for d in ('def2021', 'at_issue'):
    for lo_y, hi_y in ((2001, 2025), (2008, 2025)):
        print(f'\nReliability of p(above), {d}, outlooks {lo_y}-{hi_y}:')
        ps, os_ = [], []
        for r in noaa:
            y = int(r['year'])
            if not lo_y <= y <= hi_y: continue
            ps.append(float(r['p_above']) / 100); os_.append(int(cls(y, d) == 0))
        ps, os_ = np.array(ps), np.array(os_)
        for a, b in ((0, .3), (.3, .55), (.55, .75), (.75, 1.01)):
            m = (ps >= a) & (ps < b)
            if m.any(): print(f'  p in [{a:.2f},{b:.2f}): n={m.sum():2d} mean p={ps[m].mean():.2f} above-normal {os_[m].sum()}/{m.sum()}')
print('class counts 2001-2025 def2021:', np.bincount([cls(y, 'def2021') for y in range(2001, 2026)], minlength=3), '(above, near, below)')