"""Paired test of NOAA vs baseline interval scores (2008-2025), resampling seasons.
Scores May and August separately; pools NS/H/MH (not ACE) by normalising each variable's IS by the baseline's mean IS."""
import csv, numpy as np, importlib.util, sys, collections
spec = importlib.util.spec_from_file_location('b', 'baseline_intervals.py')
import io, contextlib
with contextlib.redirect_stdout(io.StringIO()):
    b = importlib.util.module_from_spec(spec); spec.loader.exec_module(b)
# rebuild per-year IS
obs, q, IS, count, ace, ace_before = b.obs, b.q, b.IS, b.count, b.ace, b.ace_before
recs = collections.defaultdict(dict)
for r in b.noaa:
    y = int(r['year']); m = r['month']
    if y < 2008: continue
    for k in ('ns', 'h', 'mh'):
        v = int(obs[y][k])
        recs[(y, m, k)]['NOAA'] = IS(int(r[k + '_lo']), int(r[k + '_hi']), v)
        recs[(y, m, k)]['era1995'] = IS(*q([int(obs[x][k]) for x in range(1995, y)]), v)
        recs[(y, m, k)]['clim30'] = IS(*q([int(obs[x][k]) for x in range(y - 30, y)]), v)
        if m == 'Aug':
            done = count(y, k, before='0808'); l, u = q([count(x, k, after='0808') for x in range(y - 30, y)])
            recs[(y, m, k)]['todate'] = IS(done + l, done + u, v)
            l, u = q([count(x, k, after='0808') for x in range(1995, y)])
            recs[(y, m, k)]['todate95'] = IS(done + l, done + u, v)
rng = np.random.default_rng(0); years = sorted({k[0] for k in recs})
for base in ('era1995', 'clim30', 'todate', 'todate95'):
    for m in ('May', 'Aug'):
        keys = [k for k in recs if k[1] == m and base in recs[k]]
        if not keys: continue
        # relative IS ratio: sum NOAA / sum base, pooled over variables after normalising per variable
        norm = {v: np.mean([recs[k][base] for k in keys if k[2] == v]) for v in ('ns', 'h', 'mh')}
        def stat(ys):
            a = sum(recs[(y, m, v)]['NOAA'] / norm[v] for y in ys for v in ('ns', 'h', 'mh'))
            c = sum(recs[(y, m, v)][base] / norm[v] for y in ys for v in ('ns', 'h', 'mh'))
            return a / c
        boots = [stat(rng.choice(years, len(years))) for _ in range(4000)]
        wins = sum(recs[k]['NOAA'] < recs[k][base] for k in keys); ties = sum(recs[k]['NOAA'] == recs[k][base] for k in keys)
        print(f'{m} NOAA/{base}: IS ratio {stat(years):.2f}  95% CI [{np.quantile(boots, .025):.2f}, {np.quantile(boots, .975):.2f}]  NOAA better in {wins}/{len(keys)} (ties {ties})')