1import math
 2import numpy as np
 3
 4SEED = 2517
 5A, B, KAPPA = 0.8, 3.0, 0.35
 6XMIN = 0.20
 7
 8def phi(x): return np.sin(2*np.pi*x)
 9def phip(x): return 2*np.pi*np.cos(2*np.pi*x)
10def W(x, n):
11    x=np.asarray(x,float); y=np.zeros_like(x)
12    for j in range(n+1): y += A**j*phi(B**j*x)
13    return y
14def Wprime(x,n):
15    x=np.asarray(x,float); y=np.zeros_like(x)
16    for j in range(n+1): y += (A*B)**j*phip(B**j*x)
17    return y
18def dq_direct(x,n): return (W(x,n)-A*W(B*x,n))/x
19def dq_closed(x,n): return (phi(x)-A**(n+1)*phi(B**(n+1)*x))/x
20def V(x): return .5*(x-1.5)**2
21def Vp(x): return x-1.5
22def U(x,n): return V(x)+KAPPA*W(x,n)
23def force(x,n,kind):
24    return Vp(x) + KAPPA*(Wprime(x,n) if kind=='classical' else dq_closed(x,n))
25
26def proposal(x,p,n,kind,eps,L):
27    x0,p0=x,p
28    for _ in range(L):
29        p -= .5*eps*force(x,n,kind); x += eps*p
30        if x <= XMIN or not np.isfinite(x): return x0,p0,False
31        p -= .5*eps*force(x,n,kind)
32    dH=U(x,n)+.5*p*p-U(x0,n)-.5*p0*p0
33    ok=math.log(np.random.random()) < min(0.,-dH)
34    return (x,p,True) if ok else (x0,p0,False)
35
36def chain(n,kind,eps=.01,L=4,steps=5000,burn=1000):
37    x=1.5; xs=[]; accepted=0; valid=0
38    for t in range(steps):
39        p=np.random.normal(); xn,pn,ok=proposal(x,p,n,kind,eps,L)
40        valid += int(ok)
41        if ok: accepted += int(xn != x or pn != p); x=xn
42        if t>=burn: xs.append(x)
43    z=np.asarray(xs); zc=z-z.mean(); ac1=np.dot(zc[:-1],zc[1:])/max(1e-12,np.dot(zc,zc))
44    ess=len(z)*(1-ac1)/(1+ac1)
45    return float(accepted/steps),float(valid/steps),float(z.mean()),float(ess)
46
47def main():
48    np.random.seed(SEED)
49    grid=np.linspace(XMIN,4.0,4000)
50    print('CORE_MATH')
51    for n in [0,2,4,8,12,16]:
52        err=float(np.max(np.abs(dq_direct(grid,n)-dq_closed(grid,n))))
53        md=float(np.max(np.abs(dq_closed(grid,n))))
54        mp=float(np.max(np.abs(Wprime(grid,n))))
55        print(f'N={n:2d} closure_maxerr={err:.3e} quotient_max={md:.3f} derivative_max={mp:.3f}')
56    print('SCALING_PREDICTION')
57    # Predicted quotient envelope <= 2/xmin; derivative leading envelope grows (ab)^N.
58    for n in [2,4,6,8,10,12,14,16]:
59        q=float(np.percentile(np.abs(dq_closed(grid,n)),99))
60        d=float(np.percentile(np.abs(Wprime(grid,n)),99))
61        print(f'N={n:2d} q99={q:.3f} classical_dq99={d:.3f} ratio={d/max(q,1e-9):.2f}')
62    print('SAMPLER eps=.01 L=4, exact target MH correction')
63    for n in [4,8,12,16]:
64        np.random.seed(SEED+n)
65        c=chain(n,'classical'); np.random.seed(SEED+n)
66        q=chain(n,'quotient')
67        print(f'N={n:2d} classical acc={c[0]:.3f} mean={c[2]:.3f} ESS={c[3]:.1f} quotient acc={q[0]:.3f} mean={q[2]:.3f} ESS={q[3]:.1f}')
68    print('STEP_SWEEP_N16')
69    for eps in [.002,.005,.01,.02]:
70        np.random.seed(SEED+int(eps*10000)); c=chain(16,'classical',eps=eps,steps=2500,burn=500)
71        np.random.seed(SEED+int(eps*10000)); q=chain(16,'quotient',eps=eps,steps=2500,burn=500)
72        print(f'eps={eps:.3f} classical_acc={c[0]:.3f} quotient_acc={q[0]:.3f}')
73
74if __name__=='__main__': main()