1import json, math, random
 2import numpy as np
 3import mpmath as mp
 4
 5SEED=1234
 6np.random.seed(SEED); random.seed(SEED)
 7
 8def kappa(m,beta):
 9    if m*beta <= 1: raise ValueError('requires m beta > 1')
10    return beta*beta/(8*(m*beta-1)*(2*m+1))
11
12def V(a):
13    a=np.asarray(a); return sum((a[i]-a[j])**2 for i in range(len(a)) for j in range(i+1,len(a)))
14
15def rho_det_mp(x):
16    mp.mp.dps=80
17    x=[mp.mpf(str(v)) for v in x]; m=len(x)
18    A=mp.matrix(m)
19    for i in range(m):
20        for j in range(m):
21            d=x[i]-x[j]
22            A[i,j]=1/(2*mp.pi) if d==0 else mp.sin(d/2)/(mp.pi*d)
23    return mp.det(A)
24
25def normalized_ratio_mp(e,a,beta=2):
26    e=mp.mpf(str(e)); aa=[mp.mpf(str(v)) for v in a]; m=len(aa)
27    lead=abs(e)**(beta*m*(m-1)/2)
28    for i in range(m):
29        for j in range(i+1,m): lead*=abs(aa[i]-aa[j])**beta
30    return rho_det_mp([e*v for v in aa])/lead
31
32def prior_score(x,beta=2,delta=.02,include_correction=True):
33    x=np.asarray(x,float); m=len(x); kap=kappa(m,beta) if include_correction else 0.
34    d=x[:,None]-x[None,:]; den=d*d+delta*delta
35    rep=np.divide(d,den,out=np.zeros_like(d),where=den>0)
36    np.fill_diagonal(rep,0); np.fill_diagonal(d,0)
37    return beta*np.sum(rep,axis=1)-2*kap*np.sum(d,axis=1)
38
39def langevin(kind, n=8, beta=2, steps=2800, burn=800, dt=.0015, delta=.03):
40    x=np.sort(np.random.randn(n)); samples=[]
41    for it in range(steps):
42        if kind=='baseline': g=-x
43        elif kind=='repulsion': g=-x+prior_score(x,beta,delta,False)
44        else: g=-x+prior_score(x,beta,delta,True)
45        x=x+dt*g+math.sqrt(2*dt)*np.random.randn(n)
46        if it>=burn: samples.append(np.sort(x.copy()))
47    return np.asarray(samples)
48
49def metrics(samples):
50    gaps=np.diff(samples,axis=1).ravel(); gaps=gaps/(np.mean(gaps)+1e-12)
51    return {'collision_rate_gap_lt_0.10':float(np.mean(gaps<.10)), 'gap_std':float(np.std(gaps)), 'min_gap_p05':float(np.quantile(gaps,.05))}
52
53def theorem_check(a):
54    m=len(a); kap=kappa(m,2); es=np.array([.004,.006,.009,.013,.018,.025,.035,.05])
55    rs=np.array([float(normalized_ratio_mp(e,a)) for e in es])
56    # R(e)=C + B e^2 + D e^4; theorem predicts -B/(C V)=kappa.
57    z=es**2; A=np.column_stack([np.ones(len(z)),z,z*z])
58    C,B,_=np.linalg.lstsq(A,rs,rcond=None)[0]
59    est=-B/(C*V(a))
60    return {'m':m,'theory_kappa':kap,'estimated_kappa':float(est),'relative_error':float(abs(est-kap)/kap), 'leading_constant':float(C)}
61
62def main():
63    checks=[theorem_check([0.,1.]),theorem_check([0.,1.,2.])]
64    x=np.array([-.71,-.11,.63]); b=1.7; h=1e-6
65    def logp(y): return b*sum(math.log(abs(y[i]-y[j])) for i in range(3) for j in range(i+1,3))-kappa(3,b)*V(y)
66    fd=np.array([(logp(x+np.eye(3)[i]*h)-logp(x-np.eye(3)[i]*h))/(2*h) for i in range(3)])
67    analytic=prior_score(x,b,0.0,True)
68    score={'max_abs_error':float(np.max(abs(fd-analytic))),'finite_difference':fd.tolist(),'analytic':analytic.tolist()}
69    sam={k:metrics(langevin(k)) for k in ['baseline','repulsion','full']}
70    result={'seed':SEED,'theorem_checks':checks,'score_check':score,'sampler':sam}
71    print(json.dumps(result,indent=2)); open('results.json','w').write(json.dumps(result,indent=2))
72
73if __name__=='__main__': main()