import numpy as np, json, matplotlib; matplotlib.use('Agg'); import matplotlib.pyplot as plt
exec(open('cexact_s.py').read().split("out={}")[0].replace("GAPFREE=len(sys.argv)>1","GAPFREE=True"))
src=open('cov.py').read(); exec(src.split("edges=")[0])
edges=np.array([0,.1,.2,.3,.4,.5,.6,.7,.8,.9,1.0]); res={}
sets={f'CUE N={N}':(N,np.load(f'cue_N{N}.npy')) for N in [12,25,46]}
for k,a in G.items():
    s,z1,z2,t=a.T; L=float(np.median(np.log(t/(2*np.pi)))); sets[k]=(int(round(L)),np.c_[s,np.log(z1),np.log(z2)])
fig,ax=plt.subplots(1,2,figsize=(13,5))
for i,(k,(N,a)) in enumerate(sets.items()):
    s,x,y=a[:,0],a[:,1],a[:,2]; Y=(x-y)**2/s**2; rows=[]
    for lo,hi in zip(edges[:-1],edges[1:]):
        m=(s>=lo)&(s<hi); n=m.sum()
        if n<20: continue
        sm=float(s[m].mean()); ex=np.mean([c_of_s(N,q)[0] for q in np.quantile(s[m],[.1,.3,.5,.7,.9])])
        rows.append((sm,float(Y[m].mean()),float(Y[m].std()/np.sqrt(n)),ex,int(n)))
    r=np.array(rows); chi=float((((r[:,1]-r[:,3])/r[:,2])**2).sum()); res[k]=dict(N_eff=N,chi2=chi,dof=len(r),bins=r.tolist()); print(k,'N',N,'chi2',round(chi,1),'/',len(r),[(round(q[0],2),round(q[1],3),round(q[2],3),round(q[3],3)) for q in r])
    j=0 if 'CUE' in k else 1; c=f'C{i}'
    ax[j].errorbar(r[:,0],r[:,1],r[:,2],fmt='o',color=c,label=f'{k} (data)'); ax[j].plot(r[:,0],r[:,3],'-',color=c,label=f'exact CUE N={N}')
for j in range(2): ax[j].axhline(np.pi**2/15,ls=':',c='k',label='π²/15'); ax[j].set_xlabel('s'); ax[j].set_ylabel("E[(x−y)²|s]/s²"); ax[j].legend(fontsize=7)
ax[0].set_title('CUE simulations vs exact Toeplitz formula'); ax[1].set_title('zeta zeros vs exact CUE (N = L)')
plt.tight_layout(); plt.savefig('c_check.png',dpi=130); json.dump(res,open('c_compare.json','w'),indent=1)