import math, json, hashlib
from pathlib import Path
from datetime import datetime, timezone
import numpy as np, pandas as pd

OUT=Path('/mnt/data/proton_origin_resume/work')
base=pd.read_csv(OUT/'PROTON_ORIGIN_CLS_INPUT_TABLE_2026-08-16.csv')
# Bali et al. Appendix E Table 18 AWI masses (same selected 25 ensembles)
awi={
'A654':(.006501,.000083,.018108,.000090),
'H102a':(.006436,.000067,.013789,.000061),
'N101':(.003972,.000033,.018685,.000023),
'C101':(.002427,.000027,.021210,.000033),
'C102':(.002468,.000039,.023958,.000037),
'S400':(.005679,.000028,.012605,.000028),
'N401':(.003795,.000028,.016460,.000035),
'D450':(.002077,.000018,.019429,.000019),
'B452':(.005525,.000030,.022001,.000026),
'N450':(.003659,.000018,.022044,.000019),
'D451':(.002120,.000016,.021843,.000019),
'N203':(.004739,.000016,.011051,.000013),
'N200':(.003156,.000012,.014146,.000012),
'D200':(.0015521,.0000084,.0172194,.0000086),
'E250':(.0006446,.0000073,.0188850,.0000074),
'N204':(.004792,.000013,.018906,.000012),
'N201':(.003149,.000014,.018853,.000014),
'D201':(.001551,.000015,.018872,.000015),
'N302':(.0037228,.0000078,.0090864,.0000078),
'J303':(.0020511,.0000061,.0123411,.0000047),
'E300':(.0009277,.0000055,.0145418,.0000080),
'N305':(.0054860,.0000088,.0152409,.0000083),
'N304':(.0037156,.0000093,.0152961,.0000080),
'J304':(.0020500,.0000061,.0152663,.0000044),
'J501':(.0027380,.0000036,.0071717,.0000031),
}
# ALPHA 2018 Table 5: Z_M connects RGI and bare PCAC mass; excludes beta 3.34 (outside published range)
ZM={3.40:(1.9684,.0035),3.46:(1.9935,.0027),3.55:(2.0253,.0033),3.70:(2.0630,.0038),3.85:(2.0814,.0045)}
for col in ['am_l_awi','am_l_awi_err','am_s_awi','am_s_awi_err']:
    base[col]=np.nan
for i,r in base.iterrows():
    vals=awi[r.id]
    base.loc[i,['am_l_awi','am_l_awi_err','am_s_awi','am_s_awi_err']]=vals
base=base[base.beta>=3.40].copy().reset_index(drop=True)
base['ZM']=base.beta.map(lambda b:ZM[round(float(b),2)][0])
base['ZM_err']=base.beta.map(lambda b:ZM[round(float(b),2)][1])
sqrt_t=np.sqrt(base.t0_over_a2.values)
# Central RGI coordinates. Eq(2.3) mass-dependent O(a) improvement omitted centrally and stress-tested below.
base['xl_rgi']=base.am_l_awi*base.ZM*sqrt_t
base['xs_rgi']=base.am_s_awi*base.ZM*sqrt_t
base['Y_N']=base.aMN*sqrt_t
base['sigY_N']=base.aMN_err*sqrt_t
base['cut']=1/(8*base.t0_over_a2)
base.to_csv(OUT/'PROTON_ORIGIN_CLS_RGI_QUARK_INPUT_V2_2026-08-16.csv',index=False)

xl=base.xl_rgi.values; xs=base.xs_rgi.values; Y=base.Y_N.values; sig=base.sigY_N.values; cut=base.cut.values

def design(name,l,s,c):
    one=np.ones_like(l)
    if name=='L': return np.column_stack([one,l,s,c])
    if name=='Ql': return np.column_stack([one,l,s,c,l*l])
    if name=='Qmix': return np.column_stack([one,l,s,c,l*l,s*s,l*s])
    if name=='LcutM': return np.column_stack([one,l,s,c,c*l,c*s])
    if name=='QlCutM': return np.column_stack([one,l,s,c,l*l,c*l,c*s])
    raise KeyError(name)

def fit(name,mask=None,l_override=None,s_override=None):
    if mask is None: mask=np.ones(len(Y),bool)
    ll=xl if l_override is None else l_override
    ssx=xs if s_override is None else s_override
    X=design(name,ll[mask],ssx[mask],cut[mask]); yy=Y[mask]; sy=sig[mask]
    W=1/(sy*sy); H=(X.T*W)@X
    cov0=np.linalg.pinv(H); coef=cov0@((X.T*W)@yy)
    res=yy-X@coef; chi2=float(np.sum((res/sy)**2)); dof=len(yy)-X.shape[1]
    cov=cov0*max(chi2/dof,1.0)
    k=X.shape[1]; n=len(yy); aicc=chi2+2*k+2*k*(k+1)/(n-k-1)
    return coef,cov,chi2,dof,aicc

names=['L','Ql','Qmix','LcutM','QlCutM']
fits={m:fit(m) for m in names}
a=np.array([fits[m][4] for m in names]); w=np.exp(-.5*(a-a.min())); w/=w.sum()

# Independent upstream Claim-A quark inputs: ALPHA 2020 Nf=3 RGI masses.
Mud=.00470; Mud_sig=math.sqrt(.00015**2+.00012**2)
Ms=.1270; Ms_sig=math.sqrt(.0031**2+.0032**2)
Lam=.3397960719428473; Lam_sig=.0042
sTL=.820/3.277 # sqrt(t0)*Lambda3, same frozen scale bridge
# Note Lambda cancels algebraically if using RGI masses in GeV and sqrt(t0)=sTL/Lambda, but we retain the frozen chain explicitly.
sqrt_t0=sTL/Lam
xl_phys=sqrt_t0*Mud
xs_phys=sqrt_t0*Ms

def eval_model(name,coef,l=xl_phys,s=xs_phys,st=sqrt_t0):
    X=design(name,np.array([l]),np.array([s]),np.array([0.0]))[0]
    return float((X@coef)/st)

recs={}
for i,m in enumerate(names):
    co,cv,ch,dof,ai=fits[m]; X=design(m,np.array([xl_phys]),np.array([xs_phys]),np.array([0.0]))[0]
    recs[m]={'chi2':ch,'dof':dof,'chi2_dof':ch/dof,'AICc':ai,'weight':float(w[i]),'M_GeV':eval_model(m,co),'fit_sigma_GeV':float(math.sqrt(max(0,X@cv@X))/sqrt_t0),'coef':co.tolist()}
central=sum(w[i]*recs[m]['M_GeV'] for i,m in enumerate(names))

# Robustness preferred AICc model, leave one beta and ensemble.
best=names[int(np.argmin(a))]
def pred_mask(mask):
    co,*_=fit(best,mask); return eval_model(best,co)
loo=[]
for j in range(len(base)):
    mask=np.ones(len(base),bool); mask[j]=False; loo.append((base.id.iloc[j],pred_mask(mask)))
lob=[]
for b in sorted(base.beta.unique()):
    mask=(base.beta.values!=b); lob.append((float(b),int(mask.sum()),pred_mask(mask)))

# Stress test omitted mass-dependent O(a) improvement: independently perturb l/s coordinates by beta-dependent <=1% envelopes.
rng=np.random.default_rng(20260816); N=120000
# central parameter uncertainty: quark RGI masses + scale bridge anchors + Lambda; model/fit covariance
Mud_samp=rng.normal(Mud,Mud_sig,N); Ms_samp=rng.normal(Ms,Ms_sig,N)
r0L=rng.normal(.820,.028,N); r0s=rng.normal(3.277,.039,N); Lam_samp=rng.normal(Lam,Lam_sig,N)
st=(r0L/r0s)/Lam_samp
lp=st*Mud_samp; sp=st*Ms_samp
choices=rng.choice(len(names),N,p=w); pred=np.empty(N)
for i,m in enumerate(names):
    mask=choices==i; n=mask.sum(); co,cv,*_=fits[m]
    cf=rng.multivariate_normal(co,cv,size=n)
    X=design(m,lp[mask],sp[mask],np.zeros(n)); pred[mask]=np.sum(X*cf,axis=1)/st[mask]
# Stress beta-level 1% multiplicative shifts to both quark coords, refit best each draw (few draws enough)
stress=[]
for _ in range(2000):
    factors={b:rng.normal(1,0.01) for b in base.beta.unique()}
    fac=np.array([factors[b] for b in base.beta])
    co,*_=fit(best,l_override=xl*fac,s_override=xs*fac)
    stress.append(eval_model(best,co))

q=np.quantile(pred,[.025,.16,.5,.84,.975])
record={
 'classification':'PROVENANCE_HARDENING_V2__CLS_RGI_QUARK_COORDINATES__DIRECT_ALPHA_S_TO_LAMBDA3',
 'timestamp_utc':datetime.now(timezone.utc).isoformat(),
 'n_ensembles':len(base),'betas':sorted(map(float,base.beta.unique())),
 'excluded':'A654 beta=3.34 excluded because published ALPHA Z_M Table 5 starts at beta=3.40; no extrapolation used.',
 'renormalization':{'ZM':{str(k):{'value':v[0],'sigma':v[1]} for k,v in ZM.items()},'mass_dependent_improvement':'central omitted; paper states numerically small; 1% beta-correlated coordinate stress test carried separately'},
 'physical_quark_inputs':{'Mud_RGI_GeV':Mud,'Mud_sigma_GeV':Mud_sig,'Ms_RGI_GeV':Ms,'Ms_sigma_GeV':Ms_sig,'source':'ALPHA 2020 Nf=3 RGI masses; independent upstream inputs under Claim A'},
 'scale_bridge':{'Lambda3_GeV':Lam,'Lambda3_sigma':Lam_sig,'Lambda3_method':'four-loop alpha_s running + three-loop heavy-threshold matching from frozen alpha_s(MZ)=0.1184; sigma=4.2 MeV conservative RG/matching sensitivity, target-free','sqrt_t0_Lambda3':sTL,'sqrt_t0_GeVinv_central':sqrt_t0,'xl_phys':xl_phys,'xs_phys':xs_phys},
 'models':recs,'best_AICc_model':best,'model_avg_central_GeV':float(central),
 'MC':{'mean_GeV':float(pred.mean()),'sigma_GeV':float(pred.std()),'median_GeV':float(q[2]),'68':[float(q[1]),float(q[3])],'95':[float(q[0]),float(q[4])]},
 'robustness':{'leave_one_ensemble_range_GeV':[float(min(v for _,v in loo)),float(max(v for _,v in loo))], 'leave_one_beta':lob, 'one_percent_mass_improvement_stress_mean_GeV':float(np.mean(stress)),'one_percent_mass_improvement_stress_sigma_GeV':float(np.std(stress)),'one_percent_mass_improvement_stress_range_GeV':[float(np.min(stress)),float(np.max(stress))]},
 'holdout_used':False,
 'note':'This is a post-hoc provenance hardening cross-check after historical holdout exposure, not a pristine blind prediction. No proton/neutron value enters the fit, model selection, quark physical point, or scale bridge.'
}
out=OUT/'PROTON_ORIGIN_CLS_RGI_QUARK_HARDENING_V2_2026-08-16.json'; out.write_text(json.dumps(record,indent=2,sort_keys=True)+'\n')
sha=hashlib.sha256(out.read_bytes()).hexdigest(); (OUT/'PROTON_ORIGIN_CLS_RGI_QUARK_HARDENING_2026-08-16.sha256').write_text(f'{sha}  {out.name}\n')
print(json.dumps(record,indent=2)[:18000]); print('SHA',sha)
