1import json, math
 2from pathlib import Path
 3import numpy as np
 4from scipy.optimize import least_squares
 5SEED=2295
 6
 7def mixture(p,w,mu,sig):
 8    z=p[:,None,:]-np.asarray(mu)[None,:,:]
 9    q=np.sum(z*z,axis=2); s=np.asarray(sig)
10    return np.sum(np.asarray(w)*(2*np.pi*s*s)**-1*np.exp(-q/(2*s*s)),axis=1)
11
12def make_grid(n=7,extent=4.):
13    a=np.linspace(-extent,extent,n); x,y=np.meshgrid(a,a,indexing='ij')
14    return np.c_[x.ravel(),y.ravel()],(a[1]-a[0])**2
15
16def projected_collision(p,gamma):
17    n=len(p); A=np.c_[np.ones(n),np.sum(p*p,axis=1)]
18    P=A@np.linalg.inv(A.T@A)@A.T
19    return gamma*(np.eye(n)-P)
20
21def response(L,p,k,omega):
22    n=len(p); M=L+1j*(k*p[:,0]-omega)*np.eye(n)
23    return np.mean(np.linalg.solve(M,np.ones(n)))
24
25def fit_test():
26    p,_=make_grid(11,4.5)
27    tw=np.array([.58,.27,.15]); tm=np.array([[.4,-.25],[-1,.8],[1.2,1.] ]); ts=np.array([.55,.8,.38])
28    y=mixture(p,tw,tm,ts)
29    def unpack(z):
30        q=z[:3]; w=np.exp(q-q.max()); w/=w.sum(); mu=z[3:9].reshape(3,2); sig=np.log1p(np.exp(z[9:12]))+1e-4
31        return w,mu,sig
32    def fun(z):
33        w,m,s=unpack(z); return mixture(p,w,m,s)-y
34    z=np.zeros(12); z[3:9]=np.array([[.2,0],[-.5,.4],[.7,.5]]).ravel(); z[9:]=np.log(np.expm1(.7))
35    sol=least_squares(fun,z,max_nfev=60)
36    pred=mixture(p,*unpack(sol.x));
37    c,_=make_grid(7,4.5); yc=mixture(c,tw,tm,ts)
38    ii=np.clip(np.round((p[:,0]+4.5)/(9/6)).astype(int),0,6); jj=np.clip(np.round((p[:,1]+4.5)/(9/6)).astype(int),0,6)
39    gp=yc[ii*7+jj]
40    return {'mixture_rmse':float(np.sqrt(np.mean((pred-y)**2))), 'grid_rmse':float(np.sqrt(np.mean((gp-y)**2))), 'mixture_min':float(pred.min()),'grid_min':float(gp.min())}
41
42def main():
43    p,dp=make_grid(); w=np.array([.2,.5,.3]); mu=np.array([[.4,.2],[-.8,.5],[1.1,-.7]]); sig=np.array([.45,.7,.35])
44    f=mixture(p,w,mu,sig)
45    eig=[]
46    for g in [.25,.5,1.,2.]:
47        ev=np.linalg.eigvalsh(projected_collision(p,g)); nz=ev[ev>1e-8]
48        eig.append({'gamma':g,'zero_modes':int((ev<1e-8).sum()),'median_nonzero':float(np.median(nz)),'relative_error':float(abs(np.median(nz)-g)/g)})
49    L=projected_collision(p,1.); ratios=np.array([.25,.5,1.,2.,4.]); widths=[]; peak_freq=[]
50    ws=np.linspace(-3,3,25)
51    for r in ratios:
52        amp=np.array([abs(response(L,p,r,om)) for om in ws]); peak=amp.max(); ii=np.where(amp>=peak/2)[0]
53        widths.append(float(ws[ii[-1]]-ws[ii[0]])); peak_freq.append(float(ws[amp.argmax()]))
54    out={'seed':SEED,'quadrature_points':len(p),'normalization_integral_truncated':float(f.sum()*dp),'normalization_error_due_to_box':float(abs(f.sum()*dp-1)),'min_density':float(f.min()),'eigen_sweep':eig,'crossover':{'predicted_ratio':1.,'ratios':ratios.tolist(),'fwhm':widths,'peak_frequency':peak_freq},'fit':fit_test()}
55    Path('results.json').write_text(json.dumps(out,indent=2)); print(json.dumps(out,indent=2))
56if __name__=='__main__': main()