"""Scoped empirical reconstruction, not execution of unavailable author analysis.py."""
import csv,hashlib,json,math,sys,platform,datetime,io
from pathlib import Path
import numpy as np,openpyxl
ROOT=Path(__file__).resolve().parents[1]; S=ROOT/"sources"; O=ROOT/"analysis"/"extended";O.mkdir(exist_ok=True)
book=S/"GCB2025-verified-v1.0.xlsx"
assert hashlib.sha256(book.read_bytes()).hexdigest()=="a928cf06c57576b66761d1fec8224c9a41a781eb63442afc03ba902858dd64a9"
w=openpyxl.load_workbook(book,read_only=True,data_only=True)
def numeric(sheet):
    return {int(r[0]):list(r[:8]) for r in w[sheet].iter_rows(values_only=True) if isinstance(r[0],(int,float)) and 1700<r[0]<2100}
budget=numeric("Global Carbon Budget"); hist=numeric("Historical Budget")
assert all(y in hist for y in range(1850,2025))
years=np.arange(1850,2025); emissions=np.array([hist[y][1] for y in years])
annual_co2=np.loadtxt(S/"noaa-mlo-annual.txt")
bmw=np.array([[int(r[1]),int(r[2]),float(r[3])] for line in (S/"noaa-bmw-month.txt").read_text().splitlines() if (r:=line.split()) and r[0]=="BMW" and float(r[3])>-99])
text=(S/"law2018d13c-co2.txt").read_text(); lines=[x for x in text.splitlines() if x and not x.startswith("#")]
rows=list(csv.DictReader(lines,delimiter="\t"))
firn=np.array([[float(r["age_CO2"]),float(r["d13CO2"])] for r in rows if r["age_ice"]=="firn" and 1963<=float(r["age_CO2"])<=1991])
fyr=np.array([[y,firn[firn[:,0]==y,1].mean()] for y in sorted(set(firn[:,0]))])
byr=np.array([[y,bmw[bmw[:,0]==y,2].mean(),len(bmw[bmw[:,0]==y])] for y in sorted(set(bmw[:,0])) if 1992<=y<=2014])
obs=np.vstack([fyr,byr[:,:2]])
def corr(a,b):
    r=float(np.corrcoef(a,b)[0,1]);return {"r":r,"r2":r*r}
def resid(t,y,degree):
    z=(t-t.mean())/max(1,np.ptp(t))
    return y-np.polynomial.polynomial.polyval(z,np.polynomial.polynomial.polyfit(z,y,degree))
def conv(tau,mode="constant",e=emissions):
    # value at end of each calendar year; constant input is exact within-year integration.
    decay=math.exp(-1/tau); gain=tau*(-math.expm1(-1/tau)) if mode=="constant" else 1.
    out=[]; value=0.
    for emission in e:value=decay*value+gain*emission;out.append(value)
    return np.array(out)
def model(gtcppm,mode):
    c=280+conv(60,mode)/gtcppm
    d=-6.5+(-28+6.5)*conv(8,mode)/gtcppm/c
    return c,d
runs=[]
for factor in [2.124,2.13,2.12361]:
 for mode in ["constant","year_end_impulse"]:
    c,d=model(factor,mode)
    mask=(annual_co2[:,0]>=1959)&(annual_co2[:,0]<=2023);data=annual_co2[mask]
    t=data[:,0]; o=data[:,1]; m=c[(t-1850).astype(int)]
    di=d[(obs[:,0]-1850).astype(int)]
    growtho=np.diff(o);growthm=np.diff(m)
    stats={"gtc_per_ppm":factor,"within_year_input":mode,"co2_2023":float(c[2023-1850]),"co2_raw":corr(m,o),
      "co2_cubic_residual":corr(resid(t,m,3),resid(t,o,3)),
      "co2_growth_raw":corr(growthm,growtho),
      "co2_growth_linear_residual":corr(resid(t[1:],growthm,1),resid(t[1:],growtho,1)),
      "growth_residual_sd_ratio":float(np.std(resid(t[1:],growtho,1),ddof=1)/np.std(resid(t[1:],growthm,1),ddof=1)),
      "isotope_raw":corr(di,obs[:,1]),"isotope_cubic_residual":corr(resid(obs[:,0],di,3),resid(obs[:,0],obs[:,1],3)),
      "isotope_n":len(obs),"firn_unique_years":len(fyr),"bmw_years":len(byr),
      "isotope_residual_sensitivities":[{"degree":degree,"bmw_offset_permil":offset,**corr(resid(obs[:,0],di,degree),resid(obs[:,0],obs[:,1]+np.where(obs[:,0]>=1992,offset,0),degree))} for degree in [1,2,3,4] for offset in [-.1,-.05,0,.05,.1]]}
    runs.append(stats)
# Convolutions report 2024 observed forcing; 2025 appended10.4 is stated assumption, not downloaded observation.
stocks=[]
for tau in [3.8,17.2,50,60]:
    vals=conv(tau);v25=conv(tau,e=np.r_[emissions,10.4])[-1]
    stocks.append({"tau_years":tau,"stock_end2024_GtC":float(vals[-1]),"stock_end2025_GtC_assumed_E2025_10point4":float(v25),"fraction_of_paper_907GtC":float(v25/907)})
bern=.2173*np.cumsum(emissions)
for a,tau in [(0.224,394.4),(.2824,36.54),(.2763,4.304)]:bern+=a*conv(tau)
# paper2012-2024 trend, GCB column B+C-G.
ts=np.arange(2012,2025);net=np.array([budget[y][1]+budget[y][2]-budget[y][6] for y in ts])
X=np.c_[np.ones(len(ts)),ts-ts.mean()];b=np.linalg.lstsq(X,net,rcond=None)[0];res=net-X@b
cov=np.linalg.inv(X.T@X)*(res@res/(len(ts)-2));se=math.sqrt(cov[1,1]);tstat=b[1]/se;df=len(ts)-2
xx=np.linspace(0,abs(tstat),10001);density=math.gamma((df+1)/2)/(math.sqrt(df*math.pi)*math.gamma(df/2))*(1+xx**2/df)**(-(df+1)/2)
pv=1-2*float(np.trapezoid(density,xx))
# Law Dome exact sample-ID CO2 matches, retaining records, independent of K2006 archived data.
clines=[x for x in (S/"law2018co2.txt").read_text().splitlines() if x and not x.startswith("#")]
cr=list(csv.DictReader(clines,delimiter="\t")); print("LawCO2fields",cr[0].keys())
matched=[]
for r in rows:
    candidates=[c for c in cr if c.get("SampleID")==r["SampleID"]]
    if len(candidates)==1:
        c=candidates[0]; key="CO2ppm"
        if key in c and c[key] and float(c[key])>0:matched.append([float(r["age_CO2"]),float(c[key]),float(r["d13CO2"]),r["SampleID"]])
lawfits=[]
for lo,hi in [(1500,1993),(1800,1993),(1500,1800)]:
 a=np.array([[r[0],r[1],r[2]] for r in matched if lo<=r[0]<=hi])
 if len(a)>=4:
    coef=np.polyfit(1/a[:,1],a[:,2],1)
    lawfits.append({"window":[lo,hi],"n":len(a),"intercept":float(coef[1]),**corr(1/a[:,1],a[:,2]),"note":"Rubino2019 matched sample IDs, not Kohler2006 reproduction; no age-kernel uncertainty propagation."})
for name,data,header in [("gcb-extracted.csv",[[int(y),float(e)] for y,e in zip(years,emissions)],["year","fossil_gtc_per_year"]),("combined-isotope-observations.csv",obs.tolist(),["year","delta13c_permil"]),("lawdome-sample-matches.csv",matched,["gas_year","co2_ppm","delta13c_permil","sample_id"])]:
 with (O/name).open("w",newline="") as f:ww=csv.writer(f);ww.writerow(header);ww.writerows(data)
isotope_sampling=[]
c,d=model(2.124,"constant")
for firn_mode,fdata in [("unique_gas_year_means",fyr),("individual_firn_samples",firn)]:
 for bmw_mode,bdata in [("annual_means",byr[:,:2]),("monthly_values",np.array([[r[0]+(r[1]-.5)/12,r[2]] for r in bmw if 1992<=r[0]<=2014]))]:
    paired=np.vstack([fdata,bdata]);paired=paired[np.argsort(paired[:,0])]
    prediction=np.interp(paired[:,0],years,d)
    isotope_sampling.append({"firn_sampling":firn_mode,"bmw_sampling":bmw_mode,"n":len(paired),"raw":corr(prediction,paired[:,1]),"cubic_residual":corr(resid(paired[:,0],prediction,3),resid(paired[:,0],paired[:,1],3))})
j={"status":"INDEPENDENT_RECONSTRUCTION_NOT_AUTHOR_PACKAGE_REPLAY","generated_utc":datetime.datetime.now(datetime.timezone.utc).isoformat(),
"environment":{"python":sys.version,"numpy":np.__version__,"openpyxl":openpyxl.__version__,"platform":platform.platform()},
"source_hashes":{p.name:hashlib.sha256(p.read_bytes()).hexdigest() for p in [book,S/"noaa-mlo-annual.txt",S/"noaa-bmw-month.txt",S/"law2018d13c-co2.txt",S/"law2018co2.txt"]},
"script_sha256":hashlib.sha256(Path(__file__).read_bytes()).hexdigest(),
"method":"CO2=280+conv60/Econversion; isotope=-6.5+(-28+6.5)*conv8/(Econversion*CO2), from W-H governing Eq57 (PDF), retaining signature factor absent in printed Eq74. Annual exact constant-input integration plus endpoint-impulse sensitivity. Start1850 zero perturbation. Forcing GCB2025 finalv1.0 historical fossil including cement and excluding carbonation subtraction. NOAA1959-2023, Rubino firn1963-1991 averaged within gas-year, BMW annual1992-2014 all available valid months. Both series separately detrended at paired times. No calibration harmonization/firn age smoothing added. These choices are explicit reconstructions; submitted analysis.py/raw-vintage choices unavailable.",
"data_caveats":{"bmw_month_counts":byr.tolist(),"firn_raw_n":len(firn),"gcb_end_year":int(years[-1]),"2025":"not present in verified historical series;10.4GtC/year appended only as paper-specified scenario"},
"model_runs":runs,"isotope_sampling_sensitivity":isotope_sampling,"convolutions":stocks,"bern_end2024_GtC":float(bern[-1]),
"emissions_trend":{"years":[2012,2024],"n":len(ts),"definition":"fossil + land use - cement carbonation","slope_GtC_per_year2":float(b[1]),"iid_se":se,"t":float(tstat),"df":df,"two_sided_student_p_numeric":pv,"warning":"Failure to reject slope0 does not establish equivalence/constancy; annual errors correlated and inventory uncertainties not propagated."},
"law_dome_current_vintage":lawfits}
(O/"results.json").write_text(json.dumps(j,indent=2,allow_nan=False),encoding="utf-8")
print(json.dumps(runs[0],indent=2));print("Convolutions",stocks);print("Trend",j["emissions_trend"]);print("LawDome",lawfits)