"""Replication: NOAA Eastern Pacific May outlook ranges 2008-2025 vs HURDAT2 (EP-numbered storms only,
full year). 2008-2018 ranges from CPC pages via Wayback (raw/cpc/), 2019-2025 from Wikipedia (noaa_epac.csv).
Climatological 70% interval baseline: 15-85th pct of 1995..y-1 (and y-30..y-1)."""
import csv, numpy as np
from scipy.stats import binomtest
lines = open('hurdat2-nepac-1949-2025-092926.txt').read().splitlines(); i = 0
C = {}
while i < len(lines):
    h = [x.strip() for x in lines[i].split(',')]; n = int(h[2]); sid = h[0]
    fx = [[x.strip() for x in l.split(',')] for l in lines[i+1:i+1+n]]; i += 1 + n
    if not sid.startswith('EP'): continue
    y = int(sid[4:8]); st = [f[3] for f in fx]; w = [int(f[6]) for f in fx]
    c = C.setdefault(y, [0, 0, 0])
    c[0] += any(s in ('TS', 'SS', 'HU') for s in st); c[1] += any(s == 'HU' for s in st)
    c[2] += any(s == 'HU' and v >= 96 for s, v in zip(st, w))
F = {2008: '11-16 5-8 1-3', 2009: '13-18 6-10 2-5', 2010: '9-15 4-8 1-3', 2011: '9-15 5-8 1-3', 2012: '12-18 5-9 2-5',
     2013: '11-16 5-8 1-4', 2014: '14-20 7-11 3-6', 2015: '15-22 7-12 5-8', 2016: '13-20 6-11 3-6', 2017: '14-20 6-11 3-7',
     2018: '14-20 7-12 3-7'}
for r in csv.DictReader(open('noaa_epac.csv')):
    if int(r['year']) >= 2019: F[int(r['year'])] = f"{r['ns']} {r['h']} {r['mh']}"
def IS(l, u, v, a=.3): return (u - l) + (2 / a) * max(l - v, 0) + (2 / a) * max(v - u, 0)
print('year  NOAA ranges        observed')
res = {k: [] for k in range(3)}; isn = {k: [] for k in range(3)}; isc = {k: [] for k in range(3)}
for y in sorted(F):
    rg = [tuple(map(int, s.split('-'))) for s in F[y].split()]
    print(y, F[y].ljust(18), C[y])
    for k in range(3):
        lo, hi = rg[k]; v = C[y][k]; res[k].append('in' if lo <= v <= hi else ('above' if v > hi else 'below'))
        hist = [C[x][k] for x in range(1995, y)]; l, u = np.floor(np.quantile(hist, .15)), np.ceil(np.quantile(hist, .85))
        isn[k].append(IS(lo, hi, v)); isc[k].append(IS(l, u, v))
for k, name in enumerate(('named storms', 'hurricanes', 'majors')):
    r = res[k]; hit = r.count('in')
    print(f'{name:12s} in {hit}/{len(r)} = {hit/len(r):.2f} above {r.count("above")} below {r.count("below")}  binom p vs .7 = {binomtest(hit, len(r), .7).pvalue:.3f}  IS NOAA/clim1995 = {np.mean(isn[k])/np.mean(isc[k]):.2f}')
print('all three in:', sum(all(res[k][j] == 'in' for k in range(3)) for j in range(len(F))), '/', len(F))
# season-resampled 95% CI for the pooled interval-score ratio (NS, H, MH normalised by baseline mean)
rng = np.random.default_rng(0); n = len(F)
norm = [np.mean(isc[k]) for k in range(3)]
def ratio(ix): return sum(np.sum(np.array(isn[k])[ix]) / norm[k] for k in range(3)) / sum(np.sum(np.array(isc[k])[ix]) / norm[k] for k in range(3))
bs = [ratio(rng.integers(0, n, n)) for _ in range(4000)]
print(f'pooled IS ratio NOAA/clim1995 = {ratio(np.arange(n)):.2f}  95% CI [{np.quantile(bs, .025):.2f}, {np.quantile(bs, .975):.2f}]')