"""Score RSGA 1-3 day point forecasts of F10.7 (Penticton) and Ap against GFZ Kp_ap_Ap_SN_F107 file.
Baselines known at issue (2200Z on D): persistence (value on D as reported in the RSGA: f107_obs, ap_est),
27-day recurrence (GFZ value at D+lead-27), and running 81-day mean (F10.7) / 365-day median (Ap) ending D."""
import numpy as np, pandas as pd
import os; R = os.environ.get('FLARE_ROOT', '/work/flare-forecasts')
g = pd.read_csv(f'{R}/raw/Kp_ap_Ap_SN_F107_since_1932.txt', comment='#', sep=r'\s+', header=None)
g = pd.DataFrame({'date': pd.to_datetime(dict(year=g[0], month=g[1], day=g[2])), 'Ap': g[23], 'F107': g[25]}).set_index('date')
g.loc[g.F107 <= 0, 'F107'] = np.nan; g.loc[g.Ap < 0, 'Ap'] = np.nan
g['F81'] = g.F107.rolling(81, min_periods=60).mean(); g['Apmed'] = g.Ap.rolling(365, min_periods=300).median()
rs = pd.read_csv(f'{R}/data/rsga.csv', parse_dates=['issue_date', 'day1'])
rows = []
for _, r in rs.iterrows():
    D = r.issue_date
    if pd.isna(r.day1) or (r.day1 - D).days != 1 or D not in g.index: continue
    for lead in (1, 2, 3):
        T = D + pd.Timedelta(days=lead); R27 = T - pd.Timedelta(days=27)
        if T not in g.index: continue
        rows.append(dict(issue=D, lead=lead, var='F107', fc=r[f'f107_p{lead}'], obs=g.at[T, 'F107'], pers=r.f107_obs, rec27=g.F107.get(R27), mean=g.at[D, 'F81']))
        rows.append(dict(issue=D, lead=lead, var='Ap', fc=r[f'ap_p{lead}'], obs=g.at[T, 'Ap'], pers=r.ap_est, rec27=g.Ap.get(R27), mean=g.at[D, 'Apmed']))
p = pd.DataFrame(rows).dropna(); p['year'] = p.issue.dt.year
p.to_csv(f'{R}/data/index_pairs.csv', index=False)
out = []
for (v, lead), h in p.groupby(['var', 'lead']):
    d = dict(var=v, lead=lead, n=len(h), MAE_swpc=(h.fc - h.obs).abs().mean(), RMSE_swpc=np.sqrt(((h.fc - h.obs) ** 2).mean()), bias=(h.fc - h.obs).mean())
    for b in ('pers', 'rec27', 'mean'):
        d[f'MAE_{b}'] = (h[b] - h.obs).abs().mean(); d[f'skill_vs_{b}'] = 1 - d['MAE_swpc'] / d[f'MAE_{b}']
        d[f'MSEskill_vs_{b}'] = 1 - ((h.fc - h.obs) ** 2).mean() / ((h[b] - h.obs) ** 2).mean()
    out.append(d)
s = pd.DataFrame(out); s.to_csv(f'{R}/results/index_summary.csv', index=False)
pd.set_option('display.width', 250); pd.set_option('display.max_columns', 30)
print(s.round(3).to_string(index=False))
yr = p[p.lead == 1].groupby(['var', 'year']).apply(lambda h: pd.Series(dict(n=len(h), MAE=(h.fc - h.obs).abs().mean(),
     skill_pers=1 - (h.fc - h.obs).abs().mean() / (h.pers - h.obs).abs().mean(), skill_rec27=1 - (h.fc - h.obs).abs().mean() / (h.rec27 - h.obs).abs().mean())), include_groups=False).reset_index()
yr.to_csv(f'{R}/results/index_by_year.csv', index=False); print(yr.round(3).to_string(index=False))