"""Reanalyse public Scripps station series; no author data substituted or inferred.
C: ppm CO2; delta: per mil on Scripps isotope scale. Exact abundance uses paper Rstd.
Current downloaded vintage, not byte-for-byte reproduction of the 2023 retrieval.
OLS uncertainties condition on x and omit calibration/measurement errors.
"""
import csv, datetime, hashlib, io, json, platform, sys
from pathlib import Path
import numpy as np
ROOT=Path(__file__).resolve().parents[1]; SRC=ROOT/"sources"; OUT=ROOT/"analysis"/"scripps"
OUT.mkdir(exist_ok=True)
RSTD=0.0111802
def read(name):
    text=(SRC/name).read_text(encoding="utf-8-sig")
    lines=[s for s in text.splitlines() if len(s)>4 and s[:4].isdigit()]
    rows=np.array([[float(v) for v in r] for r in csv.reader(lines)],float)
    assert rows.shape[1]==10
    keys=[(int(x[0]),int(x[1])) for x in rows]
    assert len(keys)==len(set(keys))
    return {k:r for k,r in zip(keys,rows)}
def fit(c,d,lag):
    # Center/scale x for numerical conditioning, transform intercept covariance.
    x=1000/c; mu=x.mean(); z=x-mu
    X=np.column_stack([np.ones(len(x)),z])
    bread=np.linalg.inv(X.T@X); beta=bread@X.T@d
    residual=d-X@beta; n=len(x); g=np.array([1.,-mu])
    b=float(g@beta); slope=float(beta[1])
    iid=bread*(residual@residual/(n-2))
    xu=X*residual[:,None]; meat=xu.T@xu
    for k in range(1,min(lag,n-1)+1):
        cross=xu[k:].T@xu[:-k]; meat+=(1-k/(lag+1))*(cross+cross.T)
    hac=bread@meat@bread*n/(n-2)
    r2=1-float(residual@residual)/float(((d-d.mean())**2).sum())
    se=float(np.sqrt(g@iid@g)); hse=float(np.sqrt(max(0,g@hac@g)))
    return {"n":n,"intercept_permil":b,"slope_permil_per_1000_over_ppm":slope,
            "r2":r2,"iid_se":se,"hac_se":hse,"hac_lag":lag,
            "hac_normal_95_interval":[b-1.96*hse,b+1.96*hse],
            "residual_lag1_correlation":float(np.corrcoef(residual[1:],residual[:-1])[0,1])}
def exact_abundance_intercept(c,d):
    r=RSTD*(1+d/1000); q=r/(1+r)
    a=fit(c,q,0)["intercept_permil"]
    return 1000*((a/(1-a))/RSTD-1)
def block_bootstrap(c,d):
    rng=np.random.default_rng(20261105); n=len(c); vals=[]
    for i in range(2000):
        starts=rng.integers(0,n,size=(n+4)//5)
        idx=np.concatenate([(s+np.arange(5))%n for s in starts])[:n]
        if np.ptp(c[idx])>0: vals.append(fit(c[idx],d[idx],0)["intercept_permil"])
    return list(np.quantile(vals,[.025,.975]))
results=[]; pooled={}
for station in ["ptb","ljo","mlo","spo"]:
    cfile=("monthly_flask_co2_" if station=="mlo" else "monthly_merge_co2_")+station+".csv"
    dfile="monthly_flask_c13_"+station+".csv"
    cs,ds=read(cfile),read(dfile)
    for mode,col in [("raw",4),("seasonally_adjusted",5),("filled_raw",8),("filled_adjusted",9)]:
        paired=[]
        for key in sorted(set(cs)&set(ds)):
            c,d=cs[key][col],ds[key][col]
            if np.isfinite(c) and np.isfinite(d) and c>250 and -20<d<0:
                paired.append([key[0],key[1],c,d])
        arr=np.array(paired)
        with (OUT/(station+"-"+mode+"-paired.csv")).open("w",newline="") as f:
            w=csv.writer(f); w.writerow(["year","month","co2_ppm","delta13c_permil"]);w.writerows(paired)
        for lo,hi in [(1978,2022),(1978,2023),(1978,2024),(1985,2000),(2001,2022)]:
            a=arr[(arr[:,0]>=lo)&(arr[:,0]<=hi)]
            if len(a)<24:continue
            monthly=fit(a[:,2],a[:,3],12)
            monthly.update(station=station,mode=mode,aggregation="monthly",window=[lo,hi],first=a[0,:2].tolist(),last=a[-1,:2].tolist())
            results.append(monthly)
            annual=[]
            for year in sorted(set(a[:,0])):
                rows=a[a[:,0]==year]
                if len(rows)>=9: annual.append([year,rows[:,2].mean(),rows[:,3].mean(),len(rows)])
            y=np.array(annual)
            if len(y)<8:continue
            stats=fit(y[:,1],y[:,2],3)
            stats.update(station=station,mode=mode,aggregation="annual_at_least_9_paired_months",window=[lo,hi],first_year=y[0,0],last_year=y[-1,0],
                         abundance_fit_intercept_permil=exact_abundance_intercept(y[:,1],y[:,2]))
            if mode=="seasonally_adjusted" and hi==2022 and lo==1978:
                stats["paired_block_bootstrap_5year_95_interval"]=block_bootstrap(y[:,1],y[:,2])
                with (OUT/(station+"-baseline-annual.csv")).open("w",newline="") as f:
                    w=csv.writer(f);w.writerow(["year","co2_ppm","delta13c_permil","paired_months"]);w.writerows(annual)
                pooled[station]=y
            results.append(stats)
# Deliberately label pooled standard errors as naive: stations share global signals.
p=np.concatenate([a[:,1:3] for a in pooled.values()])
pool=fit(p[:,0],p[:,1],0);pool["warning"]="Naive unequal-station pooled annual OLS; no independent-station inference or valid pooled HAC."
j={"analysis":"Scripps current-vintage sensitivity, not original author code replay","timestamp_utc":datetime.datetime.now(datetime.timezone.utc).isoformat(),
"environment":{"python":sys.version,"executable":sys.executable,"numpy":np.__version__,"platform":platform.platform()},
"script_sha256":hashlib.sha256(Path(__file__).read_bytes()).hexdigest(),
"data_hashes":{p.name:hashlib.sha256(p.read_bytes()).hexdigest() for p in sorted(SRC.glob("monthly*.csv"))},
"methods":{"missing":"Only finite co2>250 ppm and -20<delta<0; no imputation in raw/adjusted modes. Filled modes use publisher supplied cols9/10.",
"annual":"Means over SAME paired months; minimum9/months peryear.","regression":"OLS delta vs1000/CO2, centered x; intercept is x=0.",
"uncertainty":"IID and Newey-West Bartlett lag12 monthly/lag3 annual; normal95 descriptive intervals. Circular paired 5-year block bootstrap2000replicates for baseline annual; seed20261105. These omit x-error/calibration, do not establish source-identification error, and are not a trend stationarity test.",
"windows":"Fixed1978-2022 baseline, extended2023/2024, split1985-2000/2001-2022; use available paired observations only.",
"abundance":"q=R/(1+R), R=.0111802*(1+delta/1000). Fit q vs1/C then transform fitted intercept back to delta. This checks small-ratio sensitivity, not physical identification."},
"fits":results,"naive_pool":pool}
(OUT/"results.json").write_text(json.dumps(j,indent=2,allow_nan=False),encoding="utf-8")
summary=[r for r in results if r["mode"]=="seasonally_adjusted" and r["aggregation"].startswith("annual") and r["window"]==[1978,2022]]
for r in summary: print(json.dumps(r))
print("Naive pool",pool)
print("SUCCESS: saved",len(results),"fits; each paired series and baseline annual data retained")