import numpy as np, json, matplotlib; matplotlib.use('Agg'); import matplotlib.pyplot as plt
r=json.load(open('results.json')); sl=np.load('s_low.npy'); sh=np.load('s_high.npy')
w=lambda s:32/np.pi**2*s**2*np.exp(-4*s**2/np.pi)
fig,ax=plt.subplots(1,2,figsize=(13,5)); x=np.linspace(0,3.5,400)
ax[0].hist(sl,bins=200,range=(0,3.5),density=True,alpha=.6,label='first 2,001,052 zeros (Odlyzko)')
ax[0].hist(sh,bins=70,range=(0,3.5),density=True,histtype='step',lw=1.5,label='3x10^4 zeros near 1e12/1e21/1e22')
ax[0].plot(x,w(x),'k',label='GUE Wigner surmise'); ax[0].set_xlabel('normalized gap s'); ax[0].legend(); ax[0].set_title('Gap distribution')
b=r['fits']['low']['0.2']; x2=np.linspace(0,.3,200)
ax[1].hist(sl[sl<.3],bins=30,range=(0,.3),density=False,weights=np.full((sl<.3).sum(),1/(len(sl)*0.01)),alpha=.6,label='data (low)')
ax[1].plot(x2,w(x2),'k',label='Wigner surmise'); ax[1].plot(x2,np.pi**2/3*x2**2,'r--',label='GUE (pi^2/3)s^2')
ax[1].set_title(f'small-s zoom; MLE beta={b["beta"]:.3f}±{b["err"]:.3f} (s<0.2)'); ax[1].legend(); ax[1].set_xlabel('s')
plt.tight_layout(); plt.savefig('gap_hist_gue.png',dpi=130)
a=np.load('zp_pairs.npy'); s,z1,z2,t=a.T; L=np.log(t/(2*np.pi))
zn=np.sqrt(z1*z2)/L**2  # normalized slope (geometric mean of both ends)
m=s<0.3; A=np.vstack([np.ones(m.sum()),np.log(s[m])]).T; c,res_,_,_=np.linalg.lstsq(A,np.log(zn[m]),rcond=None)
e=np.log(zn[m])-A@c; cov=np.linalg.inv(A.T@A)*e.var(); 
mall=np.ones(len(s),bool); Aa=np.vstack([np.ones(len(s)),np.log(s)]).T; ca=np.linalg.lstsq(Aa,np.log(zn),rcond=None)[0]
rho=np.corrcoef(np.log(s[m]),np.log(zn[m]))[0,1]
plt.figure(figsize=(7,6)); plt.loglog(s,zn,'.',ms=2,alpha=.4); xx=np.logspace(-2.3,np.log10(.3)); plt.loglog(xx,np.exp(c[0])*xx**c[1],'r',label=f'fit s<0.3: slope {c[1]:.3f}±{np.sqrt(cov[1,1]):.3f}')
plt.xlabel('normalized gap s'); plt.ylabel("sqrt(|Z'(g1)||Z'(g2)|)/log^2(t/2pi)"); plt.legend(); plt.title("|Z'| at zero pair vs gap"); plt.savefig('zprime_vs_gap.png',dpi=130)
out=dict(slope=c[1],slope_err=np.sqrt(cov[1,1]),intercept=c[0],n=int(m.sum()),rho=rho,slope_allsample=ca[1],resid_sd=e.std())
# smallest pair
i=np.argmin(s); out['min_pair']=dict(s=s[i],z1=z1[i],z2=z2[i])
print(out); json.dump(out,open('zprime_fit.json','w'),default=float,indent=1)