"""Separate nonblind challenge checks. Synthetic cases are not empirical attribution."""
from pathlib import Path
from fractions import Fraction as F
import csv, hashlib, json, math, sys
import numpy as np

R=Path(__file__).resolve().parents[1]; O=R/'analysis/challenge'; O.mkdir(exist_ok=True)
s=F('0.0111802')
def q(delta):
    r=s*(1+F(str(delta))/1000)
    return r/(1+r)
def delta(qv): return float(1000*(qv/(1-qv)/s-1))
qf=q(-28); qi=q(-13.2); qa=q(-8.5)
qn=(qa-F('0.05')*qf)/F('0.95')
qn_in=qn+F(5,95)*(qi-qf)
exact_residual=10*qf+95*qn_in-5*qf-95*qn-5*qi
assert exact_residual==0

# Finite interval: N0=1000,F0=50, natural stock=950, F'=G=5.
# Total outflow100, each tag removed in proportion to stock.
# Fossil emission adjusts to5+100*p. Natural input=100*(1-p).
# Exact qn evolution supplies constant net signature qi at every time.
rows=[]
for t in range(21):
    N=F(1000+5*t); NF=F(50+5*t); NN=F(950); p=NF/N
    qnt=qn+F(5*t,950)*(qi-qf)
    fossil_out=100*p; natural_out=100*(1-p)
    fossil_in=5+fossil_out; natural_in=natural_out
    qnit=qnt+5*(qi-qf)/natural_in
    H=fossil_in*qf+natural_in*qnit-fossil_out*qf-natural_out*qnt
    assert H==5*qi
    qat=(NF*qf+NN*qnt)/N
    assert qat==(1000*qa+5*t*qi)/N
    assert 0<qnt<1 and 0<qnit<1
    rows.append({'t':t,'N':float(N),'f_net':1,'atmospheric_delta':delta(qat),
                 'natural_stock_delta':delta(qnt),'natural_input_delta':delta(qnit),
                 'fossil_input':float(fossil_in),'natural_input':float(natural_in)})

# Nonzero positive natural net carbon avoids relying on the undefined endpoint.
J=F(5,9); G=5+J
qnin_positive=(95*qn+G*qi-5*qf)/(95+J)
qnet_natural=(G*qi-5*qf)/J
assert 10*qf+(95+J)*qnin_positive-5*qf-95*qn==G*qi
positive={'f_net':float(5/G),'natural_net_amount':float(J),
          'gross_natural_input_delta':delta(qnin_positive),
          'effective_natural_net_delta':delta(qnet_natural)}

# Independent direct regression on original Scripps files, baseline filled-adjusted monthly.
def read_scripps(path):
    out={}
    for row in csv.reader(path.read_text(encoding='utf-8-sig').splitlines()):
        if row and row[0].strip().isdigit():
            a=list(map(float,row)); out[(int(a[0]),int(a[1]))]=a
    return out
saved=json.loads((R/'analysis/scripps/results.json').read_text())
fits=[]
for st in ['ptb','ljo','mlo','spo']:
    cf=('monthly_flask_co2_' if st=='mlo' else 'monthly_merge_co2_')+st+'.csv'
    C=read_scripps(R/'sources'/cf); D=read_scripps(R/'sources'/('monthly_flask_c13_'+st+'.csv'))
    a=np.array([[C[k][9],D[k][9]] for k in sorted(C.keys()&D.keys()) if 1978<=k[0]<=2022 and math.isfinite(C[k][9]) and math.isfinite(D[k][9]) and C[k][9]>250 and -20<D[k][9]<0])
    coef=np.polyfit(1/a[:,0],a[:,1],1)
    old=next(x for x in saved['fits'] if x['station']==st and x['mode']=='filled_adjusted' and x['aggregation']=='monthly' and x['window']==[1978,2022])
    assert len(a)==old['n'] and abs(coef[1]-old['intercept_permil'])<1e-9
    fits.append({'station':st,'n':len(a),'intercept':float(coef[1]),'difference':float(coef[1]-old['intercept_permil'])})

# Independent closed-form ramp and response calculation.
ramp60=F('10.4')/175*60*(175-60)+F('10.4')/175*60**2*F(str(math.exp(-175/60)))
k=1/3.8; ell=1/38
response=ell/(k+ell)+k/(k+ell)*math.exp(-(k+ell)*30)
result={'kind':'Nonblind independent challenge checks; synthetic balances plus raw-source regression check',
        'python':sys.version,'numpy':np.__version__,
        'snapshot_exact_fraction_balance_residual':str(exact_residual),
        'finite_interval_exact_balances_pass':True,'finite_interval_rows':rows,
        'positive_natural_net_example':positive,'direct_scripps_fits':fits,
        'ramp_tau60':float(ramp60),'two_box_response_30yr':response,
        'input_hashes':{str(p.relative_to(R)):hashlib.sha256(p.read_bytes()).hexdigest() for p in
          [R/'reviews/MATHEMATICAL-REVIEW.txt',R/'reviews/MATHEMATICAL-ADDENDUM-WH.txt',R/'reviews/CARBON-SOURCE-REVIEW.txt',R/'sources/koutsoyiannis-2024.pdf',R/'sources/kohler-2006.pdf',R/'sources/SCC-Vol.6.4-04_Cohler.pdf']},
        'script_sha256':hashlib.sha256(Path(__file__).read_bytes()).hexdigest()}
(O/'results.json').write_text(json.dumps(result,indent=2),encoding='utf-8')
print(json.dumps(result,indent=2))