"""Paired ERA5 sensitivity at the fitted d and d=0.45; other fitted quantities fixed."""
import io,json,sys,zipfile
from pathlib import Path
import numpy as np,pandas as pd
from scipy import signal,stats
root=Path(sys.argv[1]);sys.path.insert(0,str(root/'code'))
import fr,tmax
P=json.loads((root/'data/parameters.json').read_text());R=json.loads((root/'results/R2.json').read_text())
with zipfile.ZipFile(root/'data/foster-rahmstorf/AdjustData.zip') as z: obs=pd.read_csv(io.BytesIO(z.read('ObservedData.csv')))
ok=np.isfinite(obs.ERA5);tt=obs.t.to_numpy()[ok];x=obs.ERA5.to_numpy()[ok]
A=fr.Adjuster(tt,np.nan_to_num(obs.volc.to_numpy(),nan=0)[ok],obs.NINO34.to_numpy()[ok],obs.ssn.to_numpy()[ok],fit=tuple(P['fit_window']),nback=P['nback'])
adj,k,coef=A.adjust(x);post=tt>=1970;m=post & (tt<2025)
b,a=np.polyfit(tt[m]-1970,adj[m],1);smooth=fr.lowess_r(tt,adj,20/(tt.max()-tt.min()))
g0=a+b*(tt-1970);first=np.argmax(post);g0[:first]=smooth[:first]+g0[first]-smooth[first]
exo=sum(coef[j+1]*A.L[j][k[j]] for j in range(3))
prm=R['noise_models']['long_memory']['ERA5']['noise'];test=tmax.TmaxTest(55)
observed=float(test.stats(fr.annual(tt,adj,1970,2024)[1][None])[0][0]);N=2000;burn=1200
out={'dataset':'ERA5','sims_per_d':N,'seed':2026101401,'observed_tmax':observed,'fixed_ar':prm['ar'],'fixed_sd':prm['sd'],'burn_in_months':burn,'results':[]}
freq=2*np.pi*np.arange(1,26)/660;regressor=np.log(4*np.sin(freq/2)**2)
out['gph_asymptotic_regression_se_d']=float(np.sqrt(np.pi**2/6 / np.sum((regressor-regressor.mean())**2)))
for d in [prm['d'],.45]:
    rng=np.random.default_rng(out['seed']);w=np.ones(len(tt)+burn)
    for j in range(1,len(w)):w[j]=w[j-1]*(j-1+d)/j
    exceed=0;statlist=[]
    for i in range(N):
        u=signal.lfilter([1.],[1.,-prm['ar']],prm['sd']*rng.standard_normal(len(w)))
        e=signal.fftconvolve(u,w)[:len(w)][burn:]
        z=A.adjust(g0+exo+e)[0];T=float(test.stats(fr.annual(tt,z,1970,2024)[1][None])[0][0])
        statlist.append(T);exceed+=T>=observed
    out['results'].append({'d':d,'exceedances':int(exceed),'p_plus_one':(exceed+1)/(N+1),'ci95_tail':[float(stats.beta.ppf(.025,exceed,N-exceed+1)),float(stats.beta.ppf(.975,exceed+1,N-exceed))],'tmax_statistics':statlist})
print(json.dumps(out,indent=2))