"""Compare each published Cycle 25 forecast with the simple baseline available at the START of its publication year
(climatology for <= 2020; early-rise 13-mo smoothed SN regression, cycles 10-24, issued by Dec of the previous year).
Baseline distributions: Gaussian with sigma = half the 68% interval. Reads scored.csv; prints a table."""
import csv, os, numpy as np
from score_karak import crps_gauss, OBS
HERE = os.path.dirname(os.path.abspath(__file__))
S = list(csv.DictReader(open(os.path.join(HERE, 'scored.csv'), encoding='utf-8')))
L = [r for r in csv.DictReader(open(os.path.join(HERE, '..', 'baselines', 'leadtime_curve.csv')))
     if r['fit'] == '10-24' and r['predictor'] == 'ms' and r['log'] == 'False']
def base(year):
    if year <= 2020:
        return 'climatology 1-24', 178.7, (236.8 - 115.7) / 2
    cand = [r for r in L if r['issue'] <= f'{year - 1}-12']
    r = cand[-1]
    return f"early-rise k={r['k']} ({r['issue']})", float(r['point']), (float(r['hi68']) - float(r['lo68'])) / 2
print(f'{"pub year":9s} {"baseline":32s} {"b.point":>7s} {"b.sig":>6s} {"b.CRPS":>6s} {"n":>3s} {"pub med CRPS":>12s} {"frac pub better":>15s}')
tot = []
for lo, hi in [(2006, 2018), (2019, 2019), (2020, 2020), (2021, 2021), (2022, 2022), (2023, 2023), (2024, 2024)]:
    P = [s for s in S if lo <= int(s['year']) <= hi]
    name, mu, sig = base(hi)
    bc = crps_gauss(mu, sig, OBS)
    c = np.array([float(s['crps']) for s in P])
    tot += list(c < bc)
    print(f'{lo}-{hi:<4d} {name:32s} {mu:7.1f} {sig:6.1f} {bc:6.1f} {len(P):3d} {np.median(c):12.1f} {np.mean(c < bc):15.0%}')
print(f'overall fraction of published forecasts beating the time-matched baseline: {np.mean(tot):.0%} (n={len(tot)})')