"""Bounded checks of retained analyses; does not rerun acquisitions or first reviews."""
import csv, hashlib, json, math, platform, sys
from pathlib import Path
import numpy as np

ROOT = Path(__file__).resolve().parents[1]
S = ROOT / 'sources'
O = ROOT / 'analysis' / 'final-checks'
O.mkdir(exist_ok=True)
checks = {}
def sha(p): return hashlib.sha256(p.read_bytes()).hexdigest()
sr = json.loads((ROOT/'analysis/scripps/results.json').read_text())
er = json.loads((ROOT/'analysis/extended/results.json').read_text())
for label, result, script, field in [
    ('scripps',sr,'scripps_reanalysis.py','data_hashes'),
    ('extended',er,'extended_reanalysis.py','source_hashes')]:
    checks[label+'_script_identity'] = sha(ROOT/'analysis'/script) == result['script_sha256']
    checks[label+'_source_identities'] = all(sha(S/name)==h for name,h in result[field].items())
fits=[]
for station in ['ptb','ljo','mlo','spo']:
    a=np.loadtxt(ROOT/'analysis/scripps'/f'{station}-baseline-annual.csv',delimiter=',',skiprows=1)
    x=1000/a[:,1]; y=a[:,2]
    # A separate centered covariance formula, compared to saved normal-equation result.
    slope=np.dot(x-x.mean(),y-y.mean())/np.dot(x-x.mean(),x-x.mean())
    intercept=float(y.mean()-slope*x.mean())
    saved=next(r for r in sr['fits'] if r['station']==station and r['mode']=='seasonally_adjusted' and r['window']==[1978,2022] and r['aggregation'].startswith('annual'))
    fits.append({'station':station,'intercept':intercept,'difference':intercept-saved['intercept_permil'],
                 'calendar_gaps':np.diff(a[:,0]).astype(int).tolist()})
checks['baseline_regressions_alternate_formula']=all(abs(r['difference'])<1e-10 for r in fits)
def law_rows(name):
    lines=[x for x in (S/name).read_text().splitlines() if x and not x.startswith('#')]
    return list(csv.DictReader(lines,delimiter='\t'))
dr=law_rows('law2018d13c-co2.txt'); cr=law_rows('law2018co2.txt')
matched=[]
for d in dr:
    candidates=[c for c in cr if c['SampleID']==d['SampleID']]
    if len(candidates)==1 and float(candidates[0]['CO2ppm'])>0:
        c=candidates[0]
        matched.append({'id':d['SampleID'],'isotope_gas_year':float(d['age_CO2']),
                        'co2_gas_year':float(c['age_CO2']),
                        'difference':float(d['age_CO2'])-float(c['age_CO2'])})
age_mismatches=[r for r in matched if r['difference']!=0]
firn=[r for r in dr if r['age_ice']=='firn' and 1963<=float(r['age_CO2'])<=1991]
checks['firn_baseline_years_are_integers']=all(float(r['age_CO2']).is_integer() for r in firn)

# Assess fixed predictions separately from squared Pearson correlation.
forcing=np.loadtxt(ROOT/'analysis/extended/gcb-extracted.csv',delimiter=',',skiprows=1)
years=forcing[:,0]; emissions=forcing[:,1]
def conv(tau):
    value=0.; out=[]
    for e in emissions:
        value=value*math.exp(-1/tau)+e*tau*(-math.expm1(-1/tau)); out.append(value)
    return np.array(out)
c=280+conv(60)/2.124; d=-6.5-21.5*conv(8)/2.124/c
noaa=np.loadtxt(S/'noaa-mlo-annual.txt'); noaa=noaa[(noaa[:,0]>=1959)&(noaa[:,0]<=2023)]
obsiso=np.loadtxt(ROOT/'analysis/extended/combined-isotope-observations.csv',delimiter=',',skiprows=1)
def residual(t,y,degree):
    z=(t-t.mean())/np.ptp(t)
    return y-np.polynomial.polynomial.polyval(z,np.polynomial.polynomial.polyfit(z,y,degree))
def skill(pred,obs):
    err=pred-obs
    return {'pearson_r':float(np.corrcoef(pred,obs)[0,1]),
            'prediction_R2':float(1-np.sum(err**2)/np.sum((obs-obs.mean())**2)),
            'bias':float(err.mean()),'rmse':float(np.sqrt(np.mean(err**2)))}
t=noaa[:,0]; obs=noaa[:,1]; pred=c[(t-1850).astype(int)]
ti=obsiso[:,0]; oi=obsiso[:,1]; pi=np.interp(ti,years,d)
skills={'co2_raw_ppm':skill(pred,obs),
        'co2_cubic_residual_ppm':skill(residual(t,pred,3),residual(t,obs,3)),
        'growth_linear_residual_ppm_per_year':skill(residual(t[1:],np.diff(pred),1),residual(t[1:],np.diff(obs),1)),
        'isotope_raw_permil':skill(pi,oi),
        'isotope_cubic_residual_permil':skill(residual(ti,pi,3),residual(ti,oi,3))}
checks['co2_raw_correlation_matches']=abs(skills['co2_raw_ppm']['pearson_r']**2-er['model_runs'][0]['co2_raw']['r2'])<1e-12
checks['isotope_correlation_matches']=abs(skills['isotope_cubic_residual_permil']['pearson_r']-er['model_runs'][0]['isotope_cubic_residual']['r'])<1e-12
sens=er['model_runs'][0]['isotope_residual_sensitivities']
result={'status':'COMPLETED_BOUNDED_CHECKS','environment':{'python':sys.version,'numpy':np.__version__,'platform':platform.platform()},
        'script_sha256':sha(Path(__file__)), 'checks':checks,'all_checks_pass':all(checks.values()),
        'independent_baseline_formula':fits,'law_dome_matching':{'n':len(matched),'gas_age_mismatches':age_mismatches,
        'note':'Matching sample identifiers is a diagnostic. It does not account for age distributions, calibration, smoothing, repeated physical samples, or original 2006 data vintage.'},
        'fixed_prediction_skill':skills,'isotope_degree_offset_sensitivities':sens,
        'interpretation':'Prediction skill is evaluated for the stated reconstruction, not unavailable author code; bias and RMSE have units named in each key. Offset sensitivities are hypothetical stress tests. Monthly/year-end alignment and regional splice remain limitations.'}
(O/'results.json').write_text(json.dumps(result,indent=2,allow_nan=False),encoding='utf-8')
print(json.dumps(result,indent=2))
assert result['all_checks_pass']