import json, math, random, time from pathlib import Path import numpy as np SEED=1427 def seed_all(seed=SEED): random.seed(seed); np.random.seed(seed) def spectral_norm(M): return float(np.linalg.svd(M, compute_uv=False)[0]) def normalize(M, target): n=spectral_norm(M) return M*(target/n) if n else M.copy() class TwoColumnMemory: """Small MPS-inspired memory. A is the single-site transfer channel; B is a second transfer channel driven by adjacent embedded features.""" def __init__(self,p=4,chi=8,alpha=.9,beta=.9,seed=SEED): rng=np.random.default_rng(seed) self.p,self.chi=p,chi self.A=np.stack([normalize(rng.normal(size=(chi,chi)),alpha) for _ in range(p)]) self.B=np.stack([normalize(rng.normal(size=(chi,chi)),beta) for _ in range(2*p)]) def run(self,emb): h1=np.zeros(self.chi); h2=np.zeros(self.chi); prev=np.zeros(emb.shape[1]); out=[] for e in emb: At=np.tensordot(e,self.A,axes=(0,0)) z=np.concatenate([prev,e]); Bt=np.tensordot(z,self.B,axes=(0,0)) h1=At@h1+e[0]*np.ones(self.chi)/math.sqrt(self.chi) h2=Bt@h2+np.dot(prev,e)*np.ones(self.chi)/math.sqrt(self.chi) out.append(np.concatenate([e,h1,h2])); prev=e return np.asarray(out) def contraction_sweep(): """Use A=alpha*orthogonal, so the formula predicts exactly alpha**r.""" rng=np.random.default_rng(SEED); chi=12; r=20 Q,_=np.linalg.qr(rng.normal(size=(chi,chi))) rows=[] for a in [.6,.8,.95,1.0,1.1,1.25]: A=a*Q; v=rng.normal(size=chi); v/=np.linalg.norm(v) norms=[] for k in range(r+1): norms.append(np.linalg.norm(v)); v=A@v ratio=norms[-1]/norms[0] predicted=a**r # fitted per-step amplification from the whole trajectory slope=(math.log(norms[-1]+1e-30)-math.log(norms[0]))/r rows.append({'alpha':a,'r':r,'observed_ratio':ratio,'predicted_ratio':predicted, 'observed_per_step':math.exp(slope),'predicted_per_step':a}) return rows def pair_scaling_sweep(): """For a two-column perturbation, scale B by lambda. One-step response is linear.""" rng=np.random.default_rng(SEED+1); chi=10; d=6 B=normalize(rng.normal(size=(chi,chi)),.9); z=rng.normal(size=d); h=rng.normal(size=chi) base=B@h; base_norm=np.linalg.norm(base) rows=[] for lam in [0,.1,.25,.5,1.,1.5]: got=np.linalg.norm(lam*B@h) rows.append({'lambda':lam,'observed_norm':got,'predicted_norm':lam*base_norm, 'relative_error':abs(got-lam*base_norm)/(1e-12+lam*base_norm) if lam else 0.0}) return rows def boundary_sweep(): """Prediction: pair perturbations decay if gamma*lambda<1 and grow if >1.""" rows=[]; r=12; chi=10 for gamma in [.7,.9,1.0,1.1]: for lam in [.5,.9,1.0,1.1,1.4]: gain=gamma*lam # orthogonal transfer makes the finite-horizon prediction exact rng=np.random.default_rng(SEED+int(100*gamma)+int(10*lam)) Q,_=np.linalg.qr(rng.normal(size=(chi,chi))) v=rng.normal(size=chi); v/=np.linalg.norm(v) observed=np.linalg.norm(np.linalg.matrix_power(gain*Q,r)@v) predicted=gain**r rows.append({'gamma':gamma,'lambda':lam,'product':gain,'r':r, 'observed_ratio':observed,'predicted_ratio':predicted, 'predicted_regime':'decay' if gain<1 else ('neutral' if gain==1 else 'growth'), 'observed_regime':'decay' if observed<1 else ('neutral' if abs(observed-1)<1e-10 else 'growth')}) return rows def delayed_learning(seed=SEED): # Tiny nonlinear delayed-correlation task: label is XOR of bits 5 steps apart. import torch def run(device): torch.manual_seed(seed); np.random.seed(seed) T,N=18,1200; rng=np.random.default_rng(seed) x=rng.integers(0,2,size=(N,T,1)).astype('float32') y=(x[:,2,0] != x[:,7,0]).astype('float32') xt=torch.tensor(x,device=device); yt=torch.tensor(y,device=device) class GRU(torch.nn.Module): def __init__(self): super().__init__(); self.g=torch.nn.GRU(1,16,batch_first=True); self.o=torch.nn.Linear(16,1) def forward(self,x): return self.o(self.g(x)[0][:,-1,:]).squeeze(-1) class TC(torch.nn.Module): def __init__(self): super().__init__(); self.e=torch.nn.Linear(1,4); self.a=torch.nn.Parameter(torch.randn(4,8,8)*.08); self.b=torch.nn.Parameter(torch.randn(8,8,8)*.08); self.o=torch.nn.Linear(20,1) def forward(self,x): e=torch.tanh(self.e(x)); h1=torch.zeros(x.size(0),8,device=x.device); h2=h1.clone(); prev=torch.zeros_like(e[:,0]) for t in range(x.size(1)): h1=torch.bmm(torch.einsum('bp,pij->bij',e[:,t],self.a),h1.unsqueeze(-1)).squeeze(-1)+e[:,t,0:1] z=torch.cat([prev,e[:,t]],-1); h2=torch.bmm(torch.einsum('bp,pij->bij',z,self.b),h2.unsqueeze(-1)).squeeze(-1)+torch.sum(prev*e[:,t],-1,keepdim=True) prev=e[:,t] return self.o(torch.cat([e[:,-1],h1,h2],-1)).squeeze(-1) def train(m): m.to(device); opt=torch.optim.Adam(m.parameters(),lr=.01); lossfn=torch.nn.BCEWithLogitsLoss(); t0=time.time() for _ in range(80): opt.zero_grad(); loss=lossfn(m(xt),yt); loss.backward(); opt.step() with torch.no_grad(): acc=((torch.sigmoid(m(xt))>.5)==yt.bool()).float().mean().item() return {'accuracy':acc,'loss':float(loss.item()),'seconds':time.time()-t0,'device':str(device)} return {'GRU':train(GRU()),'TwoColumn':train(TC())} requested='cuda' if torch.cuda.is_available() else 'cpu' try: return run(requested) except Exception as first_error: if requested == 'cuda': try: return run('cpu') | {'fallback_reason':str(first_error)} except Exception as second_error: return {'error':str(second_error),'cuda_error':str(first_error)} return {'error':str(first_error)} def main(): seed_all(); result={'contraction':contraction_sweep(),'pair_scaling':pair_scaling_sweep(),'boundary':boundary_sweep(),'learning':delayed_learning()} Path('results.json').write_text(json.dumps(result,indent=2)) print(json.dumps(result,indent=2)) if __name__=='__main__': main()