"""Month-block bootstrap 95% intervals (200 resamples) for the scale_test.py odds ratios (M-class, lead 1)."""
import numpy as np, pandas as pd
from sklearn.linear_model import LogisticRegression
import os; R = os.environ.get('FLARE_ROOT', '/work/flare-forecasts')
rs = pd.read_csv(f'{R}/data/rsga.csv', parse_dates=['issue_date'])[['issue_date', 'f107_90d']].rename(columns={'issue_date': 'issue'})
a = pd.read_csv(f'{R}/data/flare_pairs.csv', parse_dates=['issue']); b = pd.read_csv(f'{R}/data/flare_pairs_oldscale.csv', parse_dates=['issue'])
d = a[(a.event == 'M') & (a.lead == 1)][['issue', 'p', 'y']].merge(b[(b.event == 'M') & (b.lead == 1)][['issue', 'y']].rename(columns={'y': 'y_old'}), on='issue').merge(rs, on='issue').dropna().reset_index(drop=True)
d['period'] = pd.cut(d.issue, pd.to_datetime(['1996-01-01', '2009-01-01', '2019-12-09', '2025-01-01', '2027-01-01']), labels=['C23', 'C24', 'C25a', 'C25b'], right=False)
Xd = pd.get_dummies(d.period).drop(columns='C24').astype(float); Xd['logF'] = np.log(d.f107_90d); X = Xd.values
def fit(idx):
    o = {}
    for nm, y in (('oper', d.y.values), ('old', d.y_old.values)):
        c = LogisticRegression(C=1e4, max_iter=5000).fit(X[idx], y[idx]).coef_[0]; o[f'{nm}_C25a'], o[f'{nm}_C25b'] = np.exp(c[1]), np.exp(c[2])
    XX = np.vstack([X[idx], X[idx]]); yy = np.r_[np.ones(len(idx)), np.zeros(len(idx))]; ww = np.r_[d.p.values[idx], 1 - d.p.values[idx]]
    c = LogisticRegression(C=1e4, max_iter=5000).fit(XX, yy, sample_weight=ww).coef_[0]; o['swpc_C25a'], o['swpc_C25b'] = np.exp(c[1]), np.exp(c[2])
    return o
mo = d.issue.dt.to_period('M'); groups = {m: np.where(mo == m)[0] for m in mo.unique()}; keys = list(groups); rng = np.random.default_rng(1)
pt = fit(np.arange(len(d))); bs = pd.DataFrame([fit(np.concatenate([groups[keys[i]] for i in rng.integers(0, len(keys), len(keys))])) for _ in range(200)])
out = pd.DataFrame({'point': pt, 'lo': bs.quantile(0.025), 'hi': bs.quantile(0.975)}); out.to_csv(f'{R}/results/scale_ci.csv'); print(out.round(2))
