"""Case-specific mathematical checks for Cohler & Soon (2026).
All experiments are SYNTHETIC, not atmospheric-data replication.
Uses total carbon atom amounts in arbitrary carbon-amount units.
"""
import csv, hashlib, json, math, platform, sys
from pathlib import Path
from datetime import datetime, timezone

ROOT=Path(__file__).resolve().parent
STD=0.0111802
def r(d): return STD*(1+d/1000)
def q(d): return r(d)/(1+r(d))
def delta(x): return 1000*(x/(1-x)/STD-1)
def ols(xs,ys):
    n=len(xs); xm=sum(xs)/n; ym=sum(ys)/n
    xx=sum((x-xm)**2 for x in xs); yy=sum((y-ym)**2 for y in ys)
    xy=sum((x-xm)*(y-ym) for x,y in zip(xs,ys))
    b=xy/xx; a=ym-b*xm
    return dict(intercept=a,slope=b,pearson_r=xy/math.sqrt(xx*yy),pearson_r2=xy*xy/(xx*yy))
def writecsv(name,fields,rows):
    with (ROOT/name).open('w',newline='',encoding='utf-8') as f:
        w=csv.writer(f); w.writerow(fields); w.writerows(rows)

result=dict(kind="synthetic mathematical verification; not empirical replication",
            utc=datetime.now(timezone.utc).isoformat(),python=sys.version,platform=platform.platform(),
            units="total-carbon atom amount units; time in years; isotope delta in per mil",
            script_sha256=hashlib.sha256(Path(__file__).read_bytes()).hexdigest())
# Exact mixing weights use abundance q, not delta.
conditional=[]
for dn in [-13.2,-13,-8,0,1.6]:
    f_lin=(-13.2-dn)/(-28-dn)
    f_exact=(q(-13.2)-q(dn))/(q(-28)-q(dn))
    conditional.append(dict(delta_n=dn,linear_fraction=f_lin,exact_total_amount_fraction=f_exact))
result["conditional_mixing"]=conditional
result["delta_i_interval_minus13point4_to_minus13point0"]={
    "at_delta_n_minus8":[(-13.0+8)/(-28+8),(-13.4+8)/(-28+8)],
    "at_delta_n_minus13":[0,(-13.4+13)/(-28+13)],
    "at_delta_n_zero":[13.0/28,13.4/28]}
# Constant exact qI creates a rational delta-vs-1/N curve with tiny ordinary-range bias.
n0=600.; d0=-6.5; di=-13.2
nvals=[600+i*3 for i in range(101)]
dvals=[delta(q(di)+(q(d0)-q(di))*n0/n) for n in nvals]
fit=ols([1/n for n in nvals],dvals)
chord=ols([1/nvals[0],1/nvals[-1]],[dvals[0],dvals[-1]])
predicted_chord_bias=-STD/1000/(1+r(di))*(dvals[0]-di)*(dvals[-1]-di)
result["exact_constant_net_signature_curve"]=dict(fit=fit,intercept_error=fit["intercept"]-di,
  min_delta=min(dvals),max_delta=max(dvals),chord_error=chord["intercept"]-di,
  analytical_chord_error=predicted_chord_bias,
  maximum_deviation_from_paper_curve=max(abs(d-(di+(d0-di)*n0/n)) for n,d in zip(nvals,dvals)))
assert abs((chord["intercept"]-di)-predicted_chord_bias)<1e-9
writecsv("constant-net-signature.csv",["N","delta","inverse_N"],zip(nvals,dvals,[1/n for n in nvals]))
# Exact signed-source/sink example: rise 5; fossil emission 10; natural net sink -5.
E=10.; I=95.; O=100.; F=E+I-O; df=-28.; din=-20.
qout=(E*q(df)+I*q(din)-F*q(di))/O
dout=delta(qout)
natural_qnet=(I*q(din)-O*qout)/(I-O)
qnet=(E*q(df)+I*q(din)-O*qout)/F
result["signed_flux_counterexample"]=dict(fossil_input=E,natural_input=I,natural_output=O,
  fossil_delta=df,natural_input_delta=din,natural_output_delta=dout,
  total_growth=F,natural_net_flux=I-O,net_delta=delta(qnet),
  natural_net_equivalent_delta=delta(natural_qnet),
  emission_over_net_growth=E/F,
  counterfactual_growth_if_fossil_removed_and_other_fluxes_fixed=I-O,
  isotope13_balance_residual=E*q(df)+I*q(din)-O*qout-F*q(di),
  isotope12_balance_residual=E*(1-q(df))+I*(1-q(din))-O*(1-qout)-F*(1-q(di)))
assert abs(delta(qnet)-di)<1e-9
# Provenance reading of f: identical non-fractionating removal by stock share.
# Fossil-tag stock share p=.05, removal 100, input 10 => net fossil-tag growth 5.
# Natural-tag input/removal both 95; their isotope exchange need not be zero.
p=.05; atmosphere_delta=-8.5
natural_stock_q=(q(atmosphere_delta)-p*q(df))/(1-p)
natural_input_q=natural_stock_q+(5/95)*(q(di)-q(df))
tag_net13=10*q(df)-5*q(df)
natural_net13=95*(natural_input_q-natural_stock_q)
result['origin_tag_counterexample']=dict(atmospheric_fossil_tag_share=p,
    atmosphere_delta=atmosphere_delta,fossil_input=10,natural_input=95,total_output=100,
    fossil_output=5,natural_output=95,fossil_net_growth=5,natural_net_growth=0,
    total_growth=5,fossil_net_fraction=1,
    fossil_tag_and_emission_delta=df,natural_stock_and_output_delta=delta(natural_stock_q),
    natural_input_delta=delta(natural_input_q),
    natural_net_isotope13_flux=natural_net13,
    observed_net_delta=delta((tag_net13+natural_net13)/5),
    statement='Identical nonfractionating per-atom removal for both provenance tags. Natural net carbon is zero while natural net isotope flux is nonzero; a natural residual delta cannot be assigned.')
assert abs(delta((tag_net13+natural_net13)/5)-di)<1e-9
# Same intercept while fossil fraction increases within the paper's own admissible range.
constancy=[]
for j in range(27):
    f=j/100
    dn_linear=(-13.2+28*f)/(1-f)
    dn_exact=delta((q(di)-f*q(df))/(1-f))
    constancy.append([f,dn_linear,dn_exact])
writecsv("constant-intercept-varying-fraction.csv",["f","natural_delta_linear","natural_delta_exact"],constancy)
result["constant_intercept_does_not_fix_f"]=dict(start=constancy[0],end=constancy[-1],
  statement="N signature changes within approximately [-13.2,-8]; intercept fixed; f need not remain zero")
# Arbitrary good regression does not imply instantaneous net signature constancy.
ns=[300+i*0.1 for i in range(1201)]
eps=.005; omega=2*math.pi/5
ys=[-13.2+2010/n+eps*math.sin(omega*(n-300)) for n in ns]
dy=[-2010/n**2+eps*omega*math.cos(omega*(n-300)) for n in ns]
qIs=[q(y)+n*STD/1000/(1+r(y))**2*z for n,y,z in zip(ns,ys,dy)]
local=[delta(z) for z in qIs]
result["near_linear_not_constant_net_signature"]=dict(regression=ols([1/n for n in ns],ys),
  delta_perturbation_amplitude=eps,net_delta_min=min(local),net_delta_max=max(local))
writecsv("near-linear-regression-varying-net-signature.csv",["N","delta","exact_instantaneous_net_delta"],zip(ns,ys,local))
# A regression intercept is not generally the endpoint carbon-flux-weighted mean.
N=[1.,2.,3.]
Q=[q(-7),(q(-7)+q(-10))/2,(q(-7)+q(-10)+q(-30))/3]
D=list(map(delta,Q))
result["variable_source_OLS_not_flux_mean"]=dict(N=N,delta=D,
  regression=ols([1/n for n in N],D),exact_integrated_net_delta=delta((q(-10)+q(-30))/2))
# Same transfer rates for both isotopes, return flux retained.
k=1/3.8; ell=1/38.; ainf=ell/(k+ell)
two_box=[]
for t in [0,1,3.8,10,30,60,100]:
    a=ainf+(1-ainf)*math.exp(-(k+ell)*t)
    b=1-a
    two_box.append([t,a,b,math.exp(-k*t)])
writecsv("two-box-tag-return.csv",["year","atmospheric_tag","other_box_tag","one_way_exponential"],two_box)
result["same_operator_two_box_return"]=dict(k=k,return_rate=ell,mean_initial_residence=1/k,
 eigenvalues=[0,-(k+ell)],long_run_atmospheric_fraction=ainf,rows=two_box,
 statement="No isotope-specific rates. Eq15 omits ell*B_f; local residence does not determine system perturbation response.")
# Distinct slower modes can also exist without isotope-dependent transport.
s=.001
trace=-(k+ell+s); det=k*s; disc=trace*trace-4*det
e1=(trace+math.sqrt(disc))/2; e2=(trace-math.sqrt(disc))/2
result["two_box_with_slow_deep_loss"]=dict(deep_loss_rate=s,eigenvalues=[e1,e2],timescales=[-1/e1,-1/e2])
# Eq15 linear ramp check (not the unavailable author convolution).
T=175; alpha=10.4/T
result["eq15_linear_ramp"]=[]
for tau in [3.8,50,60]:
    nf=alpha*tau*(T-tau)+alpha*tau*tau*math.exp(-T/tau)
    growth=10.4-nf/tau
    result["eq15_linear_ramp"].append(dict(tau=tau,Nf=nf,net_fossil=growth,
       f_if_total_growth5=growth/5,predicted_delta_under_paper_mixture=-13-15*growth/5))
# Monotonicity alone is not a correlation tautology.
xs=[i/1000 for i in range(1001)]
result["two_monotone_curves"]=ols(xs,[x**10 for x in xs])
# Literal mass vs atomic ratio issue using mass numbers, not precision isotope masses.
result["literal_mass_number_ratio_delta_for_atom_delta_minus8point5"]=1000*((13/12)*(1-8.5/1000)-1)
(ROOT/"results.json").write_text(json.dumps(result,indent=2),encoding='utf-8')
print(json.dumps(result,indent=2))