1import numpy as np
 2
 3
 4def energy(z):
 5    n = len(z)
 6    d = z[:, None, :] - z[None, :, :]
 7    r = np.sqrt(np.sum(d*d, axis=-1))
 8    iu = np.triu_indices(n, 1)
 9    return 0.5*np.sum(z*z) - np.sum(np.log(r[iu]))/n
10
11
12def force(z, eps=1e-5):
13    # negative gradient of H: -z + repulsive Coulomb term
14    d = z[:, None, :] - z[None, :, :]
15    r2 = np.sum(d*d, axis=-1) + eps
16    np.fill_diagonal(r2, np.inf)
17    return -z + np.sum(d/r2[:, :, None], axis=1)/len(z)
18
19
20def finite_difference_check(seed=0):
21    rng=np.random.default_rng(seed)
22    z=rng.normal(size=(7,2)); v=rng.normal(size=z.shape); v/=np.linalg.norm(v)
23    vals=[]
24    for h in [1e-2, 1e-3, 3e-4, 1e-4]:
25        fd=(energy(z+h*v)-energy(z-h*v))/(2*h)
26        analytic=-np.sum(force(z,eps=0.0)*v)
27        vals.append((h, abs(fd-analytic), abs(fd), abs(analytic)))
28    return vals
29
30
31def langevin(n, steps, beta=2.0, eta=.015, seed=0, eps=1e-5, burn=0, thin=1):
32    rng=np.random.default_rng(seed)
33    a=np.linspace(0,2*np.pi,n,endpoint=False)+rng.normal(0,.03,n)
34    z=np.c_[np.cos(a),np.sin(a)]*np.sqrt(n/2.0)
35    out=[]
36    for t in range(steps):
37        z += eta*force(z,eps) + np.sqrt(2*eta/beta)*rng.normal(size=z.shape)
38        if t>=burn and (t-burn)%thin==0: out.append(z.copy())
39    return np.asarray(out)
40
41
42def acf_iat(x, maxlag=None):
43    x=np.asarray(x,float); x=x-x.mean(); n=len(x)
44    if n<10 or np.var(x)<1e-14:return np.nan
45    maxlag=min(n//3, maxlag or n//3)
46    lags=np.arange(maxlag+1)
47    # Exactly maxlag+1 entries: lags 0,...,maxlag.
48    ac=np.correlate(x,x,mode='full')[n-1:n+maxlag]/(n-lags)
49    ac/=ac[0]
50    k=1
51    while k<len(ac) and ac[k]>0:k+=1
52    return float(max(1.,1+2*np.sum(ac[1:k])))
53
54
55def bank_stats(z, seed=123):
56    rng=np.random.default_rng(seed)
57    d=np.sqrt(np.sum((z[:,None]-z[None,:])**2,axis=-1)); np.fill_diagonal(d,np.inf)
58    nn=d.min(axis=1)
59    q=rng.normal(size=(4000,2))*np.sqrt(2)
60    assign=np.argmin(np.sum((q[:,None,:]-z[None,:,:])**2,axis=-1),axis=1)
61    counts=np.bincount(assign,minlength=len(z))/len(q)
62    entropy=-np.sum(counts[counts>0]*np.log(counts[counts>0]))/np.log(len(z))
63    dead=float(np.mean(counts==0))
64    return float(np.median(nn)), float(np.min(nn)), float(entropy), dead
65
66
67def main():
68    print('FINITE_DIFFERENCE h abs_error abs_fd abs_analytic')
69    for row in finite_difference_check():print('%g %.3e %.3e %.3e'%row)
70    print('MIXING n iat_radius iat_mean_radius median_nn min_nn entropy dead')
71    for n in [8,32,128]:
72        traj=langevin(n, steps=4500, burn=1000, thin=1, seed=10+n)
73        radius=np.linalg.norm(traj,axis=2)
74        iat=np.nanmean([acf_iat(radius[:,i]) for i in range(n)])
75        mean_iat=acf_iat(radius.mean(axis=1))
76        st=bank_stats(traj[-1],seed=90+n)
77        print(n, '%.2f %.2f %.4f %.4f %.4f %.4f'%(iat,mean_iat,*st))
78    print('GAUSSIAN_CONTROL n median_nn min_nn entropy dead')
79    for n in [8,32,128]:
80        rng=np.random.default_rng(100+n); z=rng.normal(size=(n,2))*np.sqrt(n/2)
81        print(n,'%.4f %.4f %.4f %.4f'%bank_stats(z,seed=90+n))
82
83if __name__=='__main__': main()