1import json, math, random
  2import numpy as np
  3
  4# ISS-Constrained Modular Recurrent Network: toy verification + tiny task comparison.
  5# The analytic toy recurrence is x_{k+1}=q*x_k+c, q=1+dt*(lam-gamma).
  6
  7def set_seed(seed=7):
  8    random.seed(seed); np.random.seed(seed)
  9
 10
 11def toy_verification():
 12    dt, lam = 0.1, 2.0
 13    # Prediction 1: stability boundary is gamma=lam (for positive Euler q),
 14    # with divergence for gamma < lam and contraction for gamma > lam.
 15    gammas = np.array([0.0, 1.0, 1.9, 2.0, 2.1, 3.0, 10.0])
 16    boundary_rows=[]
 17    for g in gammas:
 18        q=1+dt*(lam-g)
 19        x=1.0
 20        for _ in range(80): x=q*x
 21        boundary_rows.append({'gamma':float(g),'q_pred':float(q),'abs_q':float(abs(q)),
 22                              'final_abs':float(abs(x)),'stable_pred':bool(abs(q)<1)})
 23    # Prediction 2: log perturbation slope equals log |q| in stable regime.
 24    slope_rows=[]
 25    for g in [2.1, 3.0, 5.0]:
 26        q=1+dt*(lam-g); vals=[]; d=1.0
 27        for _ in range(45):
 28            vals.append(abs(d)); d=q*d
 29        slope=np.polyfit(np.arange(8,45), np.log(np.maximum(vals[8:],1e-300)), 1)[0]
 30        slope_rows.append({'gamma':g,'pred_log_abs_q':float(np.log(abs(q))),
 31                           'observed_slope':float(slope),'relative_error':float(abs(slope-math.log(abs(q)))/abs(math.log(abs(q))))})
 32    # Prediction 3: constant forcing reaches c/(gamma-lambda), and geometric bound
 33    # using measured contraction a=|q| is c*dt/(1-a), exactly equal here.
 34    force_rows=[]
 35    c=0.7
 36    for g in [2.1, 3.0, 5.0]:
 37        q=1+dt*(lam-g); x=0.0
 38        for _ in range(500): x=q*x+dt*c
 39        measured=abs(x); exact=c/(g-lam); bound=dt*c/(1-abs(q))
 40        force_rows.append({'gamma':g,'observed_limit':float(measured),'predicted_limit':float(exact),
 41                           'geometric_bound':float(bound),'relative_error':float(abs(measured-exact)/exact)})
 42    # Perceptual claim: z_{k+1}=alpha*z_k has decay slope log(alpha).
 43    alpha=0.93; z=1.; zs=[]
 44    for _ in range(80): zs.append(abs(z)); z*=alpha
 45    z_slope=np.polyfit(np.arange(10,80),np.log(np.maximum(zs[10:],1e-300)),1)[0]
 46    return {'boundary_sweep':boundary_rows,'decay_sweep':slope_rows,
 47            'iss_forcing_sweep':force_rows,
 48            'perception':{'alpha':alpha,'pred_log_alpha':math.log(alpha),'observed_slope':float(z_slope)}}
 49
 50
 51def task_comparison():
 52    # Small sequence regression: predict normalized sum of inputs from the final state.
 53    try:
 54        import torch
 55        import torch.nn as nn
 56        torch.set_num_threads(4)
 57        torch.manual_seed(7); np.random.seed(7)
 58        device=torch.device('cpu')  # deterministic and avoids shared cuDNN allocation failures
 59        try:
 60            # Probe CUDA and fall back on any allocation/runtime issue.
 61            if device.type=='cuda': torch.zeros(1,device=device)
 62        except Exception: device=torch.device('cpu')
 63        T,N,B,H,Z=40,1,64,24,12
 64        class Modular(nn.Module):
 65            def __init__(self):
 66                super().__init__(); self.h=H; self.z=Z; self.alpha=.92; self.dt=.1; self.gamma=1.5
 67                self.pz=nn.Linear(Z,Z,bias=False); self.pu=nn.Linear(N,Z)
 68                self.fx=nn.Linear(H,H,bias=False); self.fz=nn.Linear(Z,H); self.fu=nn.Linear(N,H); self.head=nn.Linear(H,1)
 69                nn.init.orthogonal_(self.pz.weight); self.pz.weight.data.mul_(self.alpha)
 70            def forward(self,u):
 71                b=u.shape[1]; z=torch.zeros(b,Z,device=u.device); x=torch.zeros(b,H,device=u.device)
 72                for k in range(u.shape[0]):
 73                    z=torch.tanh(self.pz(z)+self.pu(u[k]))
 74                    f=torch.tanh(self.fx(x)+self.fz(z)+self.fu(u[k]))
 75                    x=x+self.dt*(f-self.gamma*x)
 76                return self.head(x).squeeze(-1)
 77        class Vanilla(nn.Module):
 78            def __init__(self):
 79                super().__init__(); self.r=nn.GRU(N,H); self.head=nn.Linear(H,1)
 80            def forward(self,u): return self.head(self.r(u)[0][-1]).squeeze(-1)
 81        def run(model):
 82            model.to(device); opt=torch.optim.Adam(model.parameters(),lr=3e-3); losses=[]
 83            for step in range(90):
 84                u=torch.randn(T,B,N,device=device); y=u.sum(0).squeeze(-1)/math.sqrt(T)
 85                opt.zero_grad(); pred=model(u); loss=((pred-y)**2).mean(); loss.backward(); opt.step(); losses.append(float(loss.detach().cpu()))
 86            with torch.no_grad():
 87                u=torch.randn(T,256,N,device=device); y=u.sum(0).squeeze(-1)/math.sqrt(T); test=float(((model(u)-y)**2).mean().cpu())
 88            return {'final_train_mse':losses[-1],'test_mse':test,'parameters':sum(p.numel() for p in model.parameters())}
 89        try:
 90            return {'device':str(device),'vanilla_gru':run(Vanilla()),'iss_modular':run(Modular())}
 91        except Exception as first_error:
 92            # Shared CUDA environments can fail during cuDNN workspace allocation;
 93            # rerun from fresh CPU modules, preserving the experiment definition.
 94            if device.type == 'cuda':
 95                device=torch.device('cpu')
 96                return {'device':str(device),'cuda_error':repr(first_error),
 97                        'vanilla_gru':run(Vanilla()),'iss_modular':run(Modular())}
 98            raise
 99    except Exception as e:
100        return {'error':repr(e)}
101
102if __name__=='__main__':
103    set_seed(7)
104    out={'toy':toy_verification(),'task':task_comparison()}
105    with open('results.json','w') as f: json.dump(out,f,indent=2)
106    print(json.dumps(out,indent=2))