import json, math, random from pathlib import Path import numpy as np from scipy.stats import norm from sklearn.linear_model import Lasso import torch import torch.nn as nn SEED=1639 np.random.seed(SEED); random.seed(SEED); torch.manual_seed(SEED) try: device=torch.device('cuda' if torch.cuda.is_available() else 'cpu') # small tensors only; fall back if CUDA initialization fails torch.zeros(1, device=device) except Exception: device=torch.device('cpu') # True oscillator has a known linear base and an unobserved, slowly varying forcing. # The observed acceleration residual is predominantly -0.8*x^3 + 0.2*u. def true_acc(x,v,u,h): return -x - .8*x**3 + .2*u + .18*h def make_data(ntraj=80, steps=40, dt=.04, noise=.015, seed=SEED): rng=np.random.default_rng(seed); rows=[]; test=[] for q in range(ntraj): x,v,h=rng.uniform(-1.0,1.0),rng.uniform(-.7,.7),rng.uniform(-.4,.4) for k in range(steps): u=np.sin(.07*k+q*.31) a=true_acc(x,v,u,h) # observed derivative target, with measurement/derivative noise rows.append([x,v,u,a + rng.normal(0,noise)]) x=x+dt*v; v=v+dt*a; h=.995*h + .025*np.sin(.11*k+q) # independent trajectories for rollout evaluation for q in range(40): x,v,h=rng.uniform(-1.0,1.0),rng.uniform(-.7,.7),rng.uniform(-.4,.4) traj=[] for k in range(100): u=np.sin(.07*k+(q+100)*.31); traj.append((x,v,u)) a=true_acc(x,v,u,h); x=x+dt*v; v=v+dt*a; h=.995*h+.025*np.sin(.11*k+q+100) test.append(np.asarray(traj)) return np.asarray(rows), test def features(z): x,v,u=z[:,0],z[:,1],z[:,2] return np.column_stack([x,x**3,u,x*u,v,v**3]) class Residual(nn.Module): def __init__(self): super().__init__(); self.net=nn.Sequential(nn.Linear(3,32),nn.Tanh(),nn.Linear(32,32),nn.Tanh(),nn.Linear(32,1)) def forward(self,x): return self.net(x).squeeze(-1) def fit_neural(X, y, epochs=350): m=Residual().to(device); opt=torch.optim.Adam(m.parameters(),lr=.008) tx=torch.tensor(X,dtype=torch.float32,device=device); ty=torch.tensor(y,dtype=torch.float32,device=device) for _ in range(epochs): opt.zero_grad(); loss=((m(tx)-ty)**2).mean(); loss.backward(); opt.step() return m def sparse_project(P, target, alpha=.002, tau=2.5): # elastic-net-like L1 projection; standardize only for optimization, then unscale. mu=P.mean(0); sd=P.std(0)+1e-8; Q=(P-mu)/sd reg=Lasso(alpha=alpha,fit_intercept=True,max_iter=10000).fit(Q,target) b=reg.coef_/sd # OLS refit on active terms, with empirical covariance and t-statistics. active=np.flatnonzero(np.abs(b)>1e-5) selected=[]; bref=np.zeros(P.shape[1]); t=np.zeros(P.shape[1]) if len(active): A=np.column_stack([np.ones(len(P)),P[:,active]]) coef=np.linalg.lstsq(A,target,rcond=None)[0] resid=target-A@coef; dof=max(1,len(P)-A.shape[1]); s2=(resid@resid)/dof cov=s2*np.linalg.pinv(A.T@A); se=np.sqrt(np.maximum(np.diag(cov),1e-14)) for i,j in enumerate(active): bref[j]=coef[i+1]; t[j]=coef[i+1]/se[i+1] if abs(t[j])>tau: selected.append(j) return bref, selected, t def acceleration_symbolic(z,b): return features(np.asarray(z).reshape(1,3))[0]@b def rollout(model_kind,b=None,m=None,dt=.04,seed=999): rng=np.random.default_rng(seed); errs=[]; one=[] for q in range(40): x,v,h=rng.uniform(-1,1),rng.uniform(-.7,.7),rng.uniform(-.4,.4); xp,vp=x,v for k in range(100): u=np.sin(.07*k+(q+100)*.31); a=true_acc(x,v,u,h) if model_kind=='neural': with torch.no_grad(): ar=float(m(torch.tensor([[x,v,u]],dtype=torch.float32,device=device)).cpu()) else: ar=float(features(np.array([[x,v,u]]))[0]@b) # both models share known base; idea replaces neural correction by sparse field ap=-x+ar; at=a x=x+dt*v; v=v+dt*ap; h=.995*h+.025*np.sin(.11*k+q+100) xp=xp+dt*v # only used as a benign numerical reference one.append((ap-at)**2); errs.append((ap-at)**2) return float(np.sqrt(np.mean(one))), float(np.sqrt(np.mean(errs))) def mechanism_sweep(): # OLS theory predicts t(signal) proportional to sqrt(N); null t is N(0,1), # so a two-sided threshold tau has acceptance 2*Phi(-tau). rng=np.random.default_rng(44); out=[] beta=.8; sigma=.25; tau=2.5 for n in [40,80,160,320,640]: ts=[]; null_accept=[] for rep in range(180): x=rng.normal(size=n); X=np.column_stack([x,x**3,x**2]) y=beta*x**3+rng.normal(0,sigma,n) c=np.linalg.lstsq(np.column_stack([np.ones(n),X]),y,rcond=None)[0] e=y-np.column_stack([np.ones(n),X])@c; s2=e@e/(n-4); cov=s2*np.linalg.pinv(np.column_stack([np.ones(n),X]).T@np.column_stack([np.ones(n),X])) ts.append(c[2]/math.sqrt(cov[2,2])); null_accept.append(abs(c[3]/math.sqrt(cov[3,3]))>tau) out.append({'n':n,'signal_t_mean':float(np.mean(ts)),'sqrt_n_scaled':float(np.mean(ts)/math.sqrt(n)),'null_accept_rate':float(np.mean(null_accept))}) ratio=out[-1]['signal_t_mean']/out[0]['signal_t_mean'] slope=np.polyfit(np.log([x['n'] for x in out]),np.log([x['signal_t_mean'] for x in out]),1)[0] return {'predicted_signal_t_ratio_640_over_40':math.sqrt(640/40),'observed_signal_t_ratio_640_over_40':float(ratio),'predicted_t_scaling_exponent':0.5,'observed_t_scaling_exponent':float(slope),'predicted_null_accept_rate_tau_2.5':2*norm.sf(2.5),'observed_null_accept_rate_mean':float(np.mean([x['null_accept_rate'] for x in out])),'observed':out} def main(): data,test=make_data(); X=data[:,:3]; y=data[:,3] # Train on derivative residual relative to known base F_obs acceleration=-x. residual_target=y + X[:,0] m=fit_neural(X,residual_target) with torch.no_grad(): pred=m(torch.tensor(X,dtype=torch.float32,device=device)).cpu().numpy() P=features(X); b,sel,t=sparse_project(P,pred) # Refit accepted coefficients directly to noisy physical residual for fair final field. if sel: A=np.column_stack([np.ones(len(P)),P[:,sel]]); co=np.linalg.lstsq(A,residual_target,rcond=None)[0]; bf=np.zeros(6); bf[sel]=co[1:]; b=bf # Sparse regression directly on the measured residual is the non-neural baseline. b0,sel0,t0=sparse_project(P,residual_target) bs=np.zeros(6) if sel0: A0=np.column_stack([np.ones(len(P)),P[:,sel0]]) co0=np.linalg.lstsq(A0,residual_target,rcond=None)[0] bs[sel0]=co0[1:] nr=rollout('neural',m=m); sr=rollout('symbolic',b=b); sparse_r=rollout('symbolic',b=bs) result={'device':str(device),'n_train':len(X),'selected_terms':[ ['x','x3','u','xu','v','v3'][j] for j in sel], 'coefficients':b.tolist(),'t_statistics':t.tolist(), 'neural_train_rmse':float(np.sqrt(np.mean((pred-residual_target)**2))), 'rollout_rmse_acceleration':{'neural':nr[1],'sparse_only':sparse_r[1],'residual_to_symbolic':sr[1]}, 'sparse_only_terms':[['x','x3','u','xu','v','v3'][j] for j in sel0], 'mechanism_sweep':mechanism_sweep()} Path('results.json').write_text(json.dumps(result,indent=2)) print(json.dumps(result,indent=2)) if __name__=='__main__': main()