import numpy as np, json
src=open('cov.py').read(); exec(src.split("edges=")[0])
cue={N:np.load(f'cue_N{N}.npy') for N in [12,25,46]}
# extra CUE with more small gaps
from scipy.stats import unitary_group
rng=np.random.default_rng(5)
def sim(N,reps):
    out=[]
    for _ in range(reps):
        U=unitary_group.rvs(N,random_state=rng); th=np.sort(np.angle(np.linalg.eigvals(U)))
        z=np.exp(1j*th); D=np.abs(z[:,None]-z[None,:]); np.fill_diagonal(D,1); lp=np.log(D).sum(1)
        s=np.diff(np.r_[th,th[0]+2*np.pi])*N/(2*np.pi); out.append(np.c_[s,lp,np.roll(lp,-1)])
    return np.vstack(out)
for N in [12,25,46]: cue[N]=np.vstack([cue[N],sim(N,int(1.2e6/N))])
def cfit(s,d,smax=0.6):
    m=s<smax; s,d=s[m],d[m]; Y=d**2/s**2  # E[(x-y)^2]/s^2 (mean of x-y is 0 by symmetry)
    A=np.vstack([np.ones_like(s),s**2]).T; p=np.linalg.lstsq(A,Y,rcond=None)[0]
    rng=np.random.default_rng(1); B=[np.linalg.lstsq(A[k],Y[k],rcond=None)[0][0] for k in (rng.integers(0,len(s),len(s)) for _ in range(300))]
    return float(p[0]),float(np.std(B)),float(p[1]),int(m.sum())
out={}
for N,a in cue.items():
    s,l1,l2=a.T; c0,e,c2,n=cfit(s,l1-l2); ex=np.pi**2/15*(1-4/N**2); out[f'CUE N={N}']=dict(c0=c0,err=e,c2=c2,n=n,exact=ex); print(N,c0,e,c2,n,'exact',ex)
for k,a in G.items():
    s,z1,z2,t=a.T; c0,e,c2,n=cfit(s,np.log(z1)-np.log(z2)); L=float(np.median(np.log(t/(2*np.pi)))); out[k]=dict(c0=c0,err=e,c2=c2,n=n,L=L,exact_at_N_eq_L=np.pi**2/15*(1-4/L**2)); print(k,L,c0,e,c2,n)
json.dump(out,open('c_check.json','w'),indent=1)