"""Parse Pesnell 2012 (Sol. Phys. 281, 507; NTRS 20120008362) Table 2 (75 Cycle 24 predictions, SN v1, '1-sigma error')
and score interval coverage against the Cycle 24 outcome: SN v1 13-month smoothed max 81.9 (2014-04)
[v2 116.4 / 1.42; Pesnell 2016 used annual 2014 ~79]. Output: c24_table2.csv + stdout summary."""
import re, csv, numpy as np, os, sys
HERE = os.path.dirname(os.path.abspath(__file__))
sys.path[:0] = [HERE, os.path.join(HERE, '..', 'audit')]   # score_karak.py lives next to this file or in ../audit
from score_karak import crps_gauss
t = open(os.path.join(HERE, 'pesnell2012.txt'), encoding='utf-8').read().splitlines()
i0 = next(i for i, l in enumerate(t) if 'Predictions of Solar Cycle 24' in l and 'able 2' in l)
rows = []
pat = re.compile(r'^(\d+(?:\. ?\d+)?)\s*±\s*(\d+(?:\. ?\d+)?)\s+(\S+)\s+([A-Z—]|—)\s+(.*)$')
pat0 = re.compile(r'^(\d+(?:\.\d+)?)\s+(\S+)\s+([A-Z—])\s+(.*)$')
for l in t[i0:i0 + 200]:
    if l.startswith('Notes:'):
        break
    m = pat.match(l)
    if m:
        rows.append(dict(value=float(m[1].replace(' ', '')), err=float(m[2].replace(' ', '')), timing=m[3], cat=m[4], text=m[5]))
        continue
    m = pat0.match(l)
    if m and float(m[1]) > 20:
        rows.append(dict(value=float(m[1]), err='', timing=m[2], cat=m[3], text=m[4]))
with open(os.path.join(HERE, 'c24_table2.csv'), 'w', newline='') as f:
    w = csv.DictWriter(f, fieldnames=list(rows[0])); w.writeheader(); w.writerows(rows)
v = np.array([r['value'] for r in rows])
E = [r for r in rows if r['err'] != '']
for OBS in (81.9, 79.3):
    z = np.array([(OBS - r['value']) / r['err'] for r in E])
    print(f'C24 outcome {OBS}: {len(rows)} numeric (of 75), mean {v.mean():.1f}, above outcome {np.mean(v > OBS):.0%}; '
          f'{len(E)} with err: cover ±1 {np.mean(abs(z) <= 1):.0%}, ±2 {np.mean(abs(z) <= 2):.0%}, ±3 {np.mean(abs(z) <= 3):.0%}; '
          f'median z {np.median(z):.2f}, median |z| {np.median(abs(z)):.2f}, x for 68%: {np.quantile(abs(z), .68):.2f}; '
          f'median err {np.median([r["err"] for r in E]):.1f} ({np.median([r["err"]/r["value"] for r in E]):.2f} rel)')
    crps = np.array([crps_gauss(r['value'], r['err'], OBS) if r['err'] != '' else abs(OBS - r['value']) for r in rows])
    clim = crps_gauss(115, 40, OBS)   # Pesnell's climatology 115 ± 40 (v1)
    print(f'   mean CRPS {crps.mean():.1f}; climatology 115±40 CRPS {clim:.1f}; frac published worse than climatology {np.mean(crps > clim):.0%}')
from scipy.stats import spearmanr
print('Spearman(err, |error|):', spearmanr([r['err'] for r in E], [abs(81.9 - r['value']) for r in E]))