"""Score the Cycle 25 forecasts in Karak 2026 Table 1 against the outcome.

Outcome: SILSO SN v2 13-month smoothed max 160.9 (2024-10).
Inputs: karak_table1.csv (parse_karak.py); optional corrections.csv (row, value, err_plus, err_minus, err_meaning, sn_version,
        issued) from primary-source verification. Baseline intervals from ../baselines/results_c25.csv.
Outputs: scored.csv, summary printed to stdout.
"""
import csv, os, re, sys
import numpy as np
from scipy.stats import norm, spearmanr

HERE = os.path.dirname(os.path.abspath(__file__))
OBS = 160.9
V1_TO_V2 = float(os.environ.get("V1_TO_V2", 1.43))   # Karak/Nandy scaling for pre-2016 forecasts; SILSO factor is 1/0.6 = 1.67


def method_class(m):
    m = m.lower()
    if re.search(r'machine learning|\bml\b|neural|deep|lstm|long short-term|informer|neuro-fuzzy|gaussian process|learning', m):
        return 'ML / neural'
    if re.search(r'geomagnetic|aa-index|ap-index', m):
        return 'geomagnetic'
    if re.search(r'polar|dipole|sft|dynamo|faculae|helicity|rush-to-the-poles|extended solar cycle|interplanetary|heliospheric|'
                 r'flux cancelled|physics based|data assimilation|supergranular', m):
        return 'polar field / dynamo / SFT'
    if re.search(r'rise rate|increase rate|rate of increase|growth rate|slope of rising|ascending inflection|curve fitting|'
                 r'fitting the solar cycle|curve-fitting|logistic|hathaway function|superposed-epoch|bimodal', m):
        return 'rising-phase fit'
    return 'other statistical'


def crps_gauss(mu, sig, y):
    z = (y - mu) / sig
    return sig * (z * (2 * norm.cdf(z) - 1) + 2 * norm.pdf(z) - 1 / np.sqrt(np.pi))


def load(corr_path=None, scale_rule='karak'):
    rows = list(csv.DictReader(open(os.path.join(HERE, 'karak_table1.csv'), encoding='utf-8')))
    corr = {}
    if corr_path and os.path.exists(corr_path):
        corr = {r['row']: r for r in csv.DictReader(open(corr_path, encoding='utf-8'))}
    seen, out = set(), []
    for r in rows:
        if not r['value']:
            continue
        key = (r['citation'], r['value_str'])
        if key in seen:          # exact duplicate line in Karak's table
            continue
        seen.add(key)
        c = corr.get(r['row'], {})
        val = float(c.get('value') or r['value'])
        ep = c.get('err_plus', r['err_plus']); em = c.get('err_minus', r['err_minus'])
        ver = c.get('sn_version') or ('v1' if int(r['year']) <= 2015 else 'v2')
        s = V1_TO_V2 if ver == 'v1' else 1.0
        d = dict(row=int(r['row']), citation=r['citation'], year=int(r['year']), method=r['method'],
                 cls=method_class(r['method']), value_v2=val * s, sn_version=ver,
                 err_meaning=c.get('err_meaning', ''), corrected=bool(c))
        d['err_plus'] = float(ep) * s if ep else np.nan
        d['err_minus'] = float(em) * s if em else np.nan
        out.append(d)
    return out


def analyse(F, label=''):
    v = np.array([f['value_v2'] for f in F]); yr = np.array([f['year'] for f in F])
    err = OBS - v
    print(f'\n=== {label}: {len(F)} numeric forecasts (deduplicated)')
    print(f'mean {v.mean():.1f}  median {np.median(v):.1f}  sd {v.std(ddof=1):.1f}; below outcome: {np.mean(v < OBS):.0%}; '
          f'MAE {np.mean(abs(err)):.1f}')
    E = [f for f in F if not np.isnan(f['err_plus'])]
    sig = np.array([f['err_plus'] if OBS > f['value_v2'] else f['err_minus'] for f in E])
    mu = np.array([f['value_v2'] for f in E])
    z = (OBS - mu) / sig
    print(f'{len(E)} with stated error. Coverage of ±1 err: {np.mean(abs(z) <= 1):.0%}; ±2 err: {np.mean(abs(z) <= 2):.0%}; '
          f'±3 err: {np.mean(abs(z) <= 3):.0%}. Median z {np.median(z):.2f}; median |z| {np.median(abs(z)):.2f}')
    # inflation factor so that 68% / 95% of |z| fall inside
    print(f'error-bar inflation needed for 68% coverage: x{np.quantile(abs(z), 0.68):.2f}; for 95%: x{np.quantile(abs(z), 0.95):.2f}')
    rho, p = spearmanr(sig, abs(OBS - mu))
    print(f'Spearman(stated error, |actual error|) = {rho:.2f} (p={p:.2f})')
    # PIT under Gaussian reading
    pit = norm.cdf(z)
    h, _ = np.histogram(pit, bins=10, range=(0, 1))
    print('PIT histogram (10 bins, Gaussian reading of ±err as 1 sigma):', h.tolist())
    crps = crps_gauss(mu, sig, OBS)
    print(f'mean CRPS (Gaussian, sigma = stated err): {crps.mean():.1f}; median {np.median(crps):.1f}')
    for f, zz, cc in zip(E, z, crps):
        f['z'] = zz; f['crps'] = cc
    for f in F:
        f.setdefault('z', np.nan); f.setdefault('crps', abs(OBS - f['value_v2']))
    return E


def by_group(F, key, groups=None):
    gs = groups or sorted(set(f[key] for f in F))
    print(f'\n-- by {key}')
    print(f'{"group":28s} {"n":>3s} {"med pred":>8s} {"MAE":>6s} {"n_err":>5s} {"cov1":>5s} {"cov2":>5s} {"med|z|":>6s}')
    for g in gs:
        S = [f for f in F if f[key] == g] if not callable(g) else [f for f in F if g(f)]
        if not S:
            continue
        v = np.array([f['value_v2'] for f in S]); Z = np.array([f['z'] for f in S if not np.isnan(f['z'])])
        name = g if not callable(g) else g.__doc__
        c1 = f'{np.mean(abs(Z) <= 1):.0%}' if len(Z) else ''
        c2 = f'{np.mean(abs(Z) <= 2):.0%}' if len(Z) else ''
        mz = f'{np.median(abs(Z)):.2f}' if len(Z) else ''
        print(f'{name:28s} {len(S):3d} {np.median(v):8.1f} {np.mean(abs(OBS - v)):6.1f} {len(Z):5d} {c1:>5s} {c2:>5s} {mz:>6s}')


if __name__ == '__main__':
    corr = os.path.join(HERE, 'corrections.csv')
    F = load(corr if '--corrected' in sys.argv else None)
    analyse(F, 'Karak Table 1' + (' + corrections' if '--corrected' in sys.argv else ' as printed (pre-2016 x1.43)'))
    by_group(F, 'cls')
    def y1(f): '2006-2018'; return f['year'] <= 2018
    def y2(f): '2019'; return f['year'] == 2019
    def y3(f): '2020-2021'; return 2020 <= f['year'] <= 2021
    def y4(f): '2022'; return f['year'] == 2022
    def y5(f): '2023-2024'; return f['year'] >= 2023
    for g in (y1, y2, y3, y4, y5):
        g.__doc__ = g.__code__.co_consts[0]
    by_group(F, 'year', [y1, y2, y3, y4, y5])
    with open(os.path.join(HERE, 'scored.csv'), 'w', newline='', encoding='utf-8') as fh:
        keys = ['row', 'citation', 'year', 'method', 'cls', 'sn_version', 'value_v2', 'err_plus', 'err_minus', 'err_meaning',
                'corrected', 'z', 'crps']
        w = csv.DictWriter(fh, fieldnames=keys, extrasaction='ignore'); w.writeheader()
        for f in F:
            w.writerow({k: (round(f[k], 3) if isinstance(f[k], float) else f[k]) for k in keys})