import io,json,sys,zipfile
from pathlib import Path
import numpy as np,pandas as pd
from scipy.stats import beta
root=Path(sys.argv[1]);out=Path(sys.argv[2]);sys.path.insert(0,str(root/'code'))
import fr,tmax
P=json.loads((root/'data/parameters.json').read_text())
with zipfile.ZipFile(root/'data/foster-rahmstorf/AdjustData.zip') as z: obs=pd.read_csv(io.BytesIO(z.read('ObservedData.csv')))
pub=pd.read_csv(root/'data/foster-rahmstorf/GMST_adjusted.csv',comment='#')
rows=[];checks={}
for name in P['datasets']:
    ok=np.isfinite(obs[name]);tt=obs.t.to_numpy()[ok];raw=obs[name].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=A.adjust(raw)[0]
    published=pub.set_index(np.round(pub.time,4))[name+'.adj'].reindex(np.round(tt,4)).to_numpy()
    delta=published-adj
    checks[name]={'raw_max_abs_difference':float(np.nanmax(abs(delta))), 'offset':float(np.nanmean(delta)), 'centered_max_abs_difference':float(np.nanmax(abs(delta-np.nanmean(delta))))}
    for label,y in [('raw',raw),('adjusted',adj)]:
        years,ann=fr.annual(tt,y,*P['test_years']);_,_,Tk=tmax.TmaxTest(len(ann)).stats(ann[None])
        for year,v in zip(years,ann):rows.append({'dataset':name,'version':label,'year':int(year),'value':float(v)})
        checks[name][label+'_gls_tk']={str(int(k)):float(v) for k,v in zip(tmax.ks(len(ann)),Tk[0])}
pd.DataFrame(rows).to_csv(out/'annual.csv',index=False)
(out/'offset-and-gls.json').write_text(json.dumps(checks,indent=2))
draws=pd.read_csv(root/'results/null_draws.csv');observed=json.loads((root/'results/R1.json').read_text())
intervals=[]
for (model,name),g in draws.groupby(['noise_model','dataset']):
    for field in ['tmax','zq_monthly','zq_annual']:
        count=int((g[field].abs()>=abs(observed[name]['adjusted'][field])).sum());n=len(g)
        ci=[float(beta.ppf(.025,count,n-count+1)) if count else 0.,float(beta.ppf(.975,count+1,n-count)) if count<n else 1.]
        intervals.append({'model':model,'dataset':name,'statistic':field,'exceedances':count,'n':n,'plus_one_p':(count+1)/(n+1),'tail_probability_binomial_ci95':ci})
(out/'monte-carlo-intervals.json').write_text(json.dumps(intervals,indent=2))