import json, math, random from pathlib import Path import numpy as np SEED = 1305 np.random.seed(SEED); random.seed(SEED) # Overdamped bistable diffusion: dx=b(x)dt+sqrt(eps)*G dW, # U=(x^2-1)^2/4, b=-U', G=sqrt(2), hence a=2. def U(x): return 0.25*(x*x-1.0)**2 def drift(x): return x-x**3 def exact_V(x): return U(x) # V(-1)=0 def exact_grad_V(x): return x**3-x def hj_check(): x=np.linspace(-1.5,1.5,2001) b=drift(x); a=2.0; vp=exact_grad_V(x) r=b*vp+0.5*a*vp*vp # The nonzero numerical error is only floating point cancellation. return {"max_abs_residual":float(np.max(np.abs(r))), "barrier_predicted":float(exact_V(0)-exact_V(-1)), "gradient_relation_max_error":float(np.max(np.abs(vp+2*b/a)))} def pinn_fit(): # Fit a scalar neural quasipotential to HJ residual plus anchor and # positivity. Drift/covariance are deliberately estimated from bins. import torch torch.manual_seed(SEED) device="cuda" if torch.cuda.is_available() else "cpu" try: dtype=torch.float32 centers=np.linspace(-1.45,1.45,81) # Synthetic increment samples emulate measured optimizer increments. rng=np.random.default_rng(SEED); dt=.01; eps=.08; n=50000 x=rng.uniform(-1.45,1.45,n) dx=drift(x)*dt+np.sqrt(eps*2*dt)*rng.normal(size=n) bid=np.digitize(x,centers)-1 bh=[]; ah=[]; cc=[] for i,c in enumerate(centers): q=bid==i if q.sum()>50: bh.append(dx[q].mean()/dt); ah.append(dx[q].var()/dt/eps); cc.append(c) # Estimated a is 2, as it should be; interpolation in training grid. z=torch.tensor(np.array(cc)[:,None],device=device,dtype=torch.float32) bt=torch.tensor(np.array(bh),device=device,dtype=torch.float32); at=torch.tensor(np.array(ah),device=device,dtype=torch.float32) net=torch.nn.Sequential(torch.nn.Linear(1,32),torch.nn.Tanh(),torch.nn.Linear(32,32),torch.nn.Tanh(),torch.nn.Linear(32,1)).to(device) opt=torch.optim.Adam(net.parameters(),lr=.01) for _ in range(1800): zz=z.clone().requires_grad_(True); v=net(zz).squeeze(1) g=torch.autograd.grad(v.sum(),zz,create_graph=True)[0].squeeze(1) res=bt*g+.5*at*g*g pos=torch.relu(-v).square().mean() # Anchor at the left attractor; positivity is meaningful for this reference. loss=res.square().mean()+8*v[torch.argmin(torch.abs(z[:,0]+1))].square()+.2*pos opt.zero_grad(); loss.backward(); opt.step() with torch.no_grad(): pred=net(torch.tensor([[0.]],device=device,dtype=torch.float32)).item()-net(torch.tensor([[-1.]],device=device,dtype=torch.float32)).item() # Evaluate an independent collocation residual using autodiff. zz=torch.tensor(np.linspace(-1.4,1.4,401)[:,None],device=device,dtype=torch.float32,requires_grad=True) vv=net(zz); gg=torch.autograd.grad(vv.sum(),zz)[0].squeeze(1) rr=drift(zz.detach().cpu().numpy().ravel())*gg.detach().cpu().numpy()+gg.detach().cpu().numpy()**2 return {"device":device,"estimated_covariance_mean":float(np.mean(ah)), "pinn_barrier":float(pred),"pinn_independent_rms_HJ":float(np.sqrt(np.mean(rr**2))), "exact_barrier":.25} except Exception as e: return {"device":"cpu-fallback","error":str(e)} def transition_sweep(): # Vectorized Euler simulation, with reset after every left-to-right crossing. # Rate = crossings / time spent in left basin; this gives an empirical # renewal escape rate and tests log(rate)~ -DeltaV/eps. rng=np.random.default_rng(SEED+1); dt=.005; ntraj=500; steps=30000 epses=np.array([.06,.075,.09,.11,.14,.18]) rates=[]; counts=[]; exposures=[] for eps in epses: x=np.full(ntraj,-1.0); count=0; exposure=0.0 for _ in range(steps): left=x<0 exposure += left.sum()*dt noise=np.sqrt(eps*2*dt)*rng.normal(size=ntraj) x += drift(x)*dt+noise hit=left & (x>=0) count += int(hit.sum()) x[hit]=-1.0 # keep all trajectories in the renewal experiment in the left basin x[x>0]=-1.0 rates.append(count/max(exposure,1e-12)); counts.append(count); exposures.append(exposure) rates=np.array(rates); inv=1/epses slope,intercept=np.polyfit(inv,np.log(rates),1) # fit only rates with observed events (all should have events here) # Second parameter sweep: scale a=2*c. The same HJ equation predicts # DeltaV(c)=0.25/c, so the fitted slope versus 1/eps is -0.25/c. scale_rows=[] for c in np.array([0.5, 1.0, 2.0]): ee=np.array([0.08,0.11,0.15,0.20,0.28]) xx=np.full(300,-1.0); ct=0; ex=0.0 for _ in range(20000): left=xx<0; ex += left.sum()*dt xx += drift(xx)*dt+np.sqrt(epses[0]*0 + ee[0]*2*c*dt)*rng.normal(size=300) # Use a per-scale fixed epsilon list over independent blocks below. hit=left & (xx>=0); ct += int(hit.sum()); xx[hit]=-1.; xx[xx>0]=-1. # Re-simulate each epsilon; this keeps the fit transparent and seeded. rr=[] for e in ee: xx=np.full(300,-1.0); hits=0; exposure=0.0 for _ in range(16000): left=xx<0; exposure += left.sum()*dt xx += drift(xx)*dt+np.sqrt(e*2*c*dt)*rng.normal(size=300) hit=left & (xx>=0); hits += int(hit.sum()); xx[hit]=-1.; xx[xx>0]=-1. rr.append(hits/max(exposure,1e-12)) sl=np.polyfit(1/ee,np.log(np.maximum(rr,1e-15)),1)[0] scale_rows.append({"covariance_scale":float(c),"fitted_slope":float(sl), "predicted_slope":float(-.25/c), "relative_error":float(abs(sl+.25/c)/(.25/c)),"rates":rr}) return {"eps":epses.tolist(),"rates":rates.tolist(),"counts":counts, "slope_log_rate_vs_inverse_eps":float(slope),"predicted_slope":-.25, "relative_slope_error":float(abs(slope+.25)/.25),"intercept":float(intercept), "covariance_scale_sweep":scale_rows} def main(): out={"seed":SEED,"math_check":hj_check(),"pinn":pinn_fit(),"transitions":transition_sweep()} Path("results.json").write_text(json.dumps(out,indent=2)) print(json.dumps(out,indent=2)) if __name__=="__main__": main()