import json, math, random from pathlib import Path import numpy as np SEED = 2218 np.random.seed(SEED); random.seed(SEED) try: import torch torch.manual_seed(SEED) if torch.cuda.is_available(): torch.cuda.manual_seed_all(SEED) except Exception: torch = None def weighted_scores(J, groups, noise=1.0): """Return tr(J_g'WJ_g) and Frobenius leave-one-group residuals.""" A = np.asarray(J, float) / noise traces=[]; residuals=[] blocks=[A[:, g] if isinstance(g, slice) else A[:, list(g)] for g in groups] for i, block in enumerate(blocks): traces.append(float(np.sum(block * block))) others=[blocks[j] for j in range(len(blocks)) if j != i] if not others: residuals.append(float(np.sum(block * block))); continue other=np.concatenate(others, axis=1) q,_=np.linalg.qr(other, mode='reduced') proj=q @ (q.T @ block) residuals.append(float(np.sum((block-proj)**2))) return np.array(traces), np.array(residuals) def mechanism_check(): # Three quantitative predictions directly implied by I=J'WJ and projection. rng=np.random.default_rng(SEED) base=rng.normal(size=(30,3)); base[:,1]=base[:,0] # exact redundant pair groups=[slice(0,1),slice(1,2),slice(2,3)] # Prediction 1: scaling sensitivity by alpha scales information and residual by alpha^2. alphas=np.array([.25,.5,1.,2.,4.]) tr=[]; rr=[] for a in alphas: t,r=weighted_scores(a*base, groups) tr.append(t[2]); rr.append(r[2]) tr=np.array(tr); rr=np.array(rr) slope_I=np.polyfit(np.log(alphas), np.log(tr), 1)[0] slope_R=np.polyfit(np.log(alphas), np.log(rr), 1)[0] # Prediction 2: exact collinearity has zero incremental rank; residual rises as sin^2(angle). angles=np.array([0., .1, .2, .4, .8, 1.2]) residual_angle=[] x=base[:,0] # Use one fixed orthogonal direction so only the prescribed angle changes. z=rng.normal(size=x.size); z-=z@x/(x@x)*x; z/=np.linalg.norm(z) for ang in angles: J=np.column_stack([x, math.cos(ang)*x+math.sin(ang)*z, base[:,2]]) residual_angle.append(weighted_scores(J, groups)[1][1]) residual_angle=np.array(residual_angle) ratio=np.divide(residual_angle[1:], np.sin(angles[1:])**2, out=np.zeros(5), where=np.sin(angles[1:])!=0) cv=float(np.std(ratio)/np.mean(ratio)) # Prediction 3: inverse noise std weighting makes information scale as 1/sigma^2. sigmas=np.array([.5,1.,2.,4.]); infos=[] for s in sigmas: infos.append(weighted_scores(base, groups, noise=s)[0][2]) noise_slope=np.polyfit(np.log(sigmas), np.log(infos), 1)[0] return { 'alpha_values':alphas.tolist(), 'information_values':tr.tolist(), 'residual_values':rr.tolist(), 'predicted_alpha_exponent':2.0, 'observed_information_exponent':float(slope_I), 'observed_residual_exponent':float(slope_R), 'angles_rad':angles.tolist(), 'angle_residuals':residual_angle.tolist(), 'predicted_angle_law':'r = C sin(angle)^2', 'angle_law_ratio_cv':cv, 'noise_sigmas':sigmas.tolist(), 'noise_information':infos, 'predicted_noise_exponent':-2.0, 'observed_noise_exponent':float(noise_slope), 'pass': bool(abs(slope_I-2)<.03 and abs(slope_R-2)<.03 and cv<.03 and abs(noise_slope+2)<.03) } class UnitODE(torch.nn.Module): def __init__(self, width): super().__init__(); self.width=width self.w=torch.nn.Parameter(torch.randn(width)*.5) self.b=torch.nn.Parameter(torch.randn(width)*.1) self.a=torch.nn.Parameter(torch.randn(width)*.5/width**.5) def rhs(self,x): return (self.a*torch.tanh(self.w*x+self.b)).sum(-1,keepdim=True) def forward(self,x0,dt,steps): x=x0 out=[x] for _ in range(steps): k1=self.rhs(x); k2=self.rhs(x+dt*k1/2); k3=self.rhs(x+dt*k2/2); k4=self.rhs(x+dt*k3) x=x+dt*(k1+2*k2+2*k3+k4)/6; out.append(x) return torch.cat(out,dim=1) def train(width, x0, target, epochs=350, keep=None, device='cpu'): model=UnitODE(width).to(device) opt=torch.optim.Adam(model.parameters(),lr=.025) dt=float(2.5/(target.shape[1]-1)) x0=x0.to(device); target=target.to(device) for _ in range(epochs): opt.zero_grad(); pred=model(x0,dt,target.shape[1]-1) loss=((pred-target)**2).mean(); loss.backward(); opt.step() return model def unit_sensitivity(model, x0, times, device): # Observed-output sensitivity columns using finite differences of actual trajectories. # This is equivalent to local S propagated through the fixed-step ODE integrator. model.eval(); x0=x0.to(next(model.parameters()).device); dt=float(2.5/(len(times)-1)); base=model(x0,dt,len(times)-1).detach().cpu().numpy().ravel() cols=[] with torch.no_grad(): for j in range(model.width): vals=[] for p in [model.w,model.b,model.a]: old=float(p[j]); eps=1e-4*max(1.,abs(old)); p[j]=old+eps vals.append((model(x0,dt,len(times)-1).detach().cpu().numpy().ravel()-base)/eps) p[j]=old cols.append(np.stack(vals,axis=1)) return np.stack(cols,axis=1), base def prune_and_score(kind, width=16, retain=8, device='cpu'): t=np.linspace(0,2.5,31); x0=torch.tensor([[.7]],dtype=torch.float32); true=torch.zeros(1,len(t)); true[0,0]=.7 # stable nonlinear target generated by known scalar dynamics, integrated with RK4 externally x=.7 for k in range(30): f=lambda q: -0.8*q+0.35*q*q k1=f(x); k2=f(x+(t[1]-t[0])*k1/2); k3=f(x+(t[1]-t[0])*k2/2); k4=f(x+(t[1]-t[0])*k3) x += (t[1]-t[0])*(k1+2*k2+2*k3+k4)/6; true[0,k+1]=x model=train(width,x0,true,device=device) J,_=unit_sensitivity(model,x0,t,device) # J shape observations, units, 3 params; group score sums parameter columns groups=[slice(3*j,3*j+3) for j in range(width)] Jflat=J.reshape(J.shape[0],-1) traces,res=weighted_scores(Jflat,groups) if kind=='sensitivity': chosen=np.argsort(traces+res)[:width-retain] elif kind=='random': chosen=np.random.default_rng(SEED).choice(width, width-retain, replace=False) elif kind=='magnitude': chosen=np.argsort(np.array([abs(model.a[j].item())+abs(model.w[j].item())+abs(model.b[j].item()) for j in range(width)]))[:width-retain] else: chosen=np.arange(width-retain) mask=np.ones(width,dtype=bool); mask[chosen]=False # Evaluate the selected compact subnetwork, then retrain it at the same budget. compact=UnitODE(retain).to(device) with torch.no_grad(): keep=np.flatnonzero(mask) compact.w.copy_(model.w[keep]); compact.b.copy_(model.b[keep]); compact.a.copy_(model.a[keep]) opt=torch.optim.Adam(compact.parameters(),lr=.025) xd=x0.to(device); yd=true.to(device) for _ in range(350): opt.zero_grad(); pred0=compact(xd,2.5/30,30); loss=((pred0-yd)**2).mean(); loss.backward(); opt.step() with torch.no_grad(): pred=compact(xd,2.5/30,30).cpu().numpy()[0] mse=float(np.mean((pred-true.numpy()[0])**2)); longerr=float(abs(pred[-1]-true.numpy()[0,-1])) return {'retained_units':retain,'parameters':3*retain,'mse':mse,'terminal_error':longerr, 'removed_mean_information':float(np.mean(traces[chosen])),'removed_mean_residual':float(np.mean(res[chosen]))} def main(): device='cuda' if torch is not None and torch.cuda.is_available() else 'cpu' try: mech=mechanism_check(); exp={k:prune_and_score(k,device=device) for k in ['sensitivity','magnitude','random']} except Exception as e: device='cpu'; mech=mechanism_check(); exp={k:prune_and_score(k,device=device) for k in ['sensitivity','magnitude','random']} exp['error_fallback']=repr(e) out={'seed':SEED,'device':device,'mechanism':mech,'pruning_experiment':exp} Path('results.json').write_text(json.dumps(out,indent=2)) print(json.dumps(out,indent=2)) if __name__=='__main__': main()