"""Compare NOAA 70% ranges with climatological 70% intervals available at issue time.
Baselines (for outlook year y):
  clim30   : empirical 15th-85th percentile of y-30..y-1 full-year counts
  era1995  : same over 1995..y-1 (the 'active era' NOAA itself cites)
  todate   (August only): count already formed by Aug 8 of year y + 15th-85th pct of the
             Aug 8-Dec 31 counts over y-30..y-1 (i.e. 'what has happened + climatology of the rest')
  todate95 (August only): same, rest-of-season climatology over 1995..y-1
Scores: coverage, mean width, and the interval score IS_0.3 (Gneiting & Raftery 2007)."""
import csv, numpy as np, collections
storms = list(csv.DictReader(open('../data/storms.csv')))
obs = {int(r['year']): r for r in csv.DictReader(open('../data/season_counts.csv'))}
ace = collections.defaultdict(float); ace_before = collections.defaultdict(float)
for r in csv.DictReader(open('../data/ace_fixes.csv')):
    y = int(r['date'][:4]); ace[y] += float(r['ace'])
    if r['date'][4:] < '0808': ace_before[y] += float(r['ace'])
def count(y, k, before=None, after=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   # storms counted in the year they first reach the category
        md = t[4:8]
        if before and md >= before: continue
        if after and md < after: continue
        n += 1
    return n
# note: HURDAT sid year == year of first named fix for all modern storms; check
def q(vals, lo=0.15, hi=0.85):
    v = np.array(vals, float); return np.floor(np.quantile(v, lo)), np.ceil(np.quantile(v, hi))
def IS(l, u, y, a=0.3): return (u - l) + (2 / a) * max(l - y, 0) + (2 / a) * max(y - u, 0)
noaa = list(csv.DictReader(open('../data/noaa_outlooks.csv')))
res = collections.defaultdict(list)
for r in noaa:
    y = int(r['year']); m = r['month']
    for k in ('ns', 'h', 'mh'):
        v = int(obs[y][k])
        cands = {'NOAA': (int(r[k + '_lo']), int(r[k + '_hi'])),
                 'clim30': q([int(obs[x][k]) for x in range(y - 30, y)]),
                 'era1995': q([int(obs[x][k]) for x in range(1995, y)])}
        if m == 'Aug':
            done = count(y, k, before='0808')
            rest = [count(x, k, after='0808') for x in range(y - 30, y)]
            l, u = q(rest); cands['todate'] = (done + l, done + u)
            l, u = q([count(x, k, after='0808') for x in range(1995, y)]); cands['todate95'] = (done + l, done + u)
        for name, (l, u) in cands.items():
            res[(name, m, k, y >= 2008)].append((l <= v <= u, u - l, IS(l, u, v)))
    # ACE as raw units
    v = ace[y]; med = 87.5 if y <= 2010 else (92.4 if y <= 2020 else 96.7)
    if r['ace_lo_pct']:
        cands = {'NOAA': (float(r['ace_lo_pct']) * med / 100, float(r['ace_hi_pct']) * med / 100),
                 'clim30': tuple(np.quantile([ace[x] for x in range(y - 30, y)], [.15, .85])),
                 'era1995': tuple(np.quantile([ace[x] for x in range(1995, y)], [.15, .85]))}
        if m == 'Aug':
            rest = [ace[x] - ace_before[x] for x in range(y - 30, y)]
            l, u = np.quantile(rest, [.15, .85]); cands['todate'] = (ace_before[y] + l, ace_before[y] + u)
            l, u = np.quantile([ace[x] - ace_before[x] for x in range(1995, y)], [.15, .85]); cands['todate95'] = (ace_before[y] + l, ace_before[y] + u)
        for name, (l, u) in cands.items():
            res[(name, m, 'ace', y >= 2008)].append((l <= v <= u, u - l, IS(l, u, v)))
out = []
for era in (True, False):
    print('\n=== ' + ('2008-2025 (probability stated)' if era else '2001-2007'))
    for m in ('May', 'Aug'):
        for k in ('ns', 'h', 'mh', 'ace'):
            line = f'{m} {k:4s}'
            for name in ('NOAA', 'clim30', 'era1995', 'todate', 'todate95'):
                xs = res.get((name, m, k, era))
                if not xs: continue
                cov = np.mean([x[0] for x in xs]); wid = np.mean([x[1] for x in xs]); isc = np.mean([x[2] for x in xs])
                line += f' | {name}: cov {cov:.2f} w {wid:5.1f} IS {isc:6.1f}'
                out.append(dict(era='2008-2025' if era else '2001-2007', month=m, var=k, method=name, n=len(xs), coverage=round(cov, 3), mean_width=round(wid, 2), interval_score=round(isc, 2)))
            print(line)
with open('../data/interval_baselines.csv', 'w', newline='') as fh:
    w = csv.DictWriter(fh, out[0].keys()); w.writeheader(); w.writerows(out)