1import json, math, random
  2from pathlib import Path
  3import numpy as np
  4
  5SEED = 1379
  6np.random.seed(SEED)
  7random.seed(SEED)
  8
  9
 10def simulate_front(D=0.25, r=1.0, K=1.0, kappa=0.0, nx=1501, nt=3000, dx=1.0,
 11                   dt=0.2, seed=0, impulse_x=150, impulse_amplitude=None, noise=True):
 12    """Euler-Maruyama Fisher-KPP field with reflecting boundaries and clipping."""
 13    if D > 0:
 14        dt = min(dt, 0.8 * dx * dx / (2 * D))
 15    rng = np.random.default_rng(seed)
 16    u = np.zeros(nx, dtype=np.float64)
 17    u[impulse_x] = K if impulse_amplitude is None else impulse_amplitude
 18    history = [u.copy()]
 19    for _ in range(nt):
 20        lap = np.zeros_like(u)
 21        if D > 0:
 22            lap[1:-1] = (u[2:] - 2*u[1:-1] + u[:-2]) / dx**2
 23            lap[0] = 2*(u[1]-u[0]) / dx**2
 24            lap[-1] = 2*(u[-2]-u[-1]) / dx**2
 25        react = r * u * (1.0 - u / K)
 26        stochastic = 0.0
 27        if noise and kappa > 0 and r > 0:
 28            # Formula in the prompt, with spatial white-noise scaling 1/dx.
 29            var = np.maximum(0.0, 2*kappa*r*dt/dx * u*(1-u/K))
 30            stochastic = np.sqrt(var) * rng.standard_normal(nx)
 31        u = np.clip(u + dt*(D*lap + react) + stochastic, 0.0, K)
 32        history.append(u.copy())
 33    return np.asarray(history), dt
 34
 35
 36def front_position(profile, threshold=0.1, dx=1.0):
 37    ids = np.where(profile >= threshold)[0]
 38    if len(ids) == 0: return 0.0
 39    j = ids[-1]
 40    if j >= len(profile)-1: return j*dx
 41    if j == 0: return 0.0
 42    a,b=profile[j],profile[j+1]
 43    return (j + ((threshold-a)/(b-a) if b != a else 0.0))*dx
 44
 45
 46def estimate_speed(hist, dt, threshold=0.1, start_frac=.30, end_frac=.75):
 47    ts=np.arange(len(hist))*dt
 48    xs=np.array([front_position(p,threshold) for p in hist])
 49    lo,hi=int(len(ts)*start_frac),int(len(ts)*end_frac)
 50    return float(np.polyfit(ts[lo:hi],xs[lo:hi],1)[0]),ts,xs
 51
 52
 53def math_checks():
 54    # Prediction 1: v* = 2 sqrt(D r), tested after a long transient.
 55    speed_rows=[]
 56    for D in (0.10,0.25,0.50):
 57        h,dt=simulate_front(D=D,r=1.0,nt=3000)
 58        v,_,_=estimate_speed(h,dt)
 59        pred=2*math.sqrt(D)
 60        speed_rows.append({'D':D,'predicted':pred,'observed':v,'relative_error':abs(v-pred)/pred})
 61    # Prediction 2: linearized zero-state growth rate is r. Use a small perturbation.
 62    stability=[]
 63    for r in (-0.5,0.0,0.5):
 64        h,dt=simulate_front(D=0.0,r=r,K=1.0,nt=30,impulse_amplitude=1e-3)
 65        amp=h[:,150]
 66        slope=float(np.polyfit(np.arange(1,20)*dt,np.log(np.maximum(amp[1:20],1e-30)),1)[0])
 67        stability.append({'r':r,'predicted_log_slope':r,'observed_log_slope':slope,'absolute_error':abs(slope-r)})
 68    # Prediction 3: increasing kappa broadens the stochastic front.
 69    widths=[]; means=[]
 70    for kap in (0.0,0.15,0.5):
 71        vals=[]; early=[]
 72        for seed in range(24):
 73            h,dt=simulate_front(D=.25,r=.8,kappa=kap,nt=120,seed=seed,nx=401,impulse_x=40)
 74            p=h[-1]; right=np.where(p>=.1)[0]; core=np.where(p>=.9)[0]
 75            vals.append((right[-1]-core[0]) if len(right) and len(core) else 0)
 76            early.append(float(h[15].sum()))
 77        widths.append({'kappa':kap,'width_mean':float(np.mean(vals)),'width_std':float(np.std(vals))})
 78        means.append({'kappa':kap,'early_total_mean':float(np.mean(early))})
 79    return {'speed_scaling':speed_rows,'stability_boundary':stability,'noise_broadening':widths,'noise_early_mean':means}
 80
 81
 82def task_experiment(device='cpu'):
 83    import torch
 84    torch.manual_seed(SEED); np.random.seed(SEED)
 85    n,length,d,layers=96,16,24,6
 86    x=torch.randint(0,2,(n,length),dtype=torch.float32,device=device); y=(x.sum(1)%2).long()
 87    class Net(torch.nn.Module):
 88        def __init__(self,gated):
 89            super().__init__(); self.gated=gated
 90            self.inp=torch.nn.Linear(1,d)
 91            self.blocks=torch.nn.ModuleList([torch.nn.Sequential(torch.nn.Linear(d,d),torch.nn.Tanh(),torch.nn.Linear(d,d)) for _ in range(layers)])
 92            self.out=torch.nn.Linear(d,2)
 93            if gated: self.rho=torch.nn.Parameter(torch.tensor(-.5))
 94        def forward(self,x):
 95            h=self.inp(x.unsqueeze(-1))
 96            if self.gated:
 97                u=torch.clamp(x*.5+.1,0,1).mean(0).detach().clone()
 98                D=.12; r=torch.nn.functional.softplus(self.rho); dt=.4
 99                for b in self.blocks:
100                    lap=torch.zeros_like(u); lap[1:-1]=u[2:]-2*u[1:-1]+u[:-2]; lap[0]=u[1]-u[0]; lap[-1]=u[-2]-u[-1]
101                    u=torch.clamp(u+dt*(D*lap+r*u*(1-u)),0,1)
102                    h=h+(u/(1+u))[None,:,None]*b(h)
103            else:
104                for b in self.blocks: h=h+b(h)
105            return self.out(h.mean(1))
106    out={}
107    for gated in (False,True):
108        torch.manual_seed(SEED); net=Net(gated).to(device); opt=torch.optim.Adam(net.parameters(),lr=.01); losses=[]
109        for _ in range(180):
110            opt.zero_grad(); loss=torch.nn.functional.cross_entropy(net(x),y); loss.backward(); opt.step(); losses.append(float(loss))
111        out['gated' if gated else 'baseline']={'final_loss':losses[-1],'best_loss':min(losses),'params':sum(p.numel() for p in net.parameters())}
112    return out
113
114
115def main():
116    result={'math':math_checks()}
117    try:
118        import torch
119        device='cuda' if torch.cuda.is_available() else 'cpu'
120        try: result['task']=task_experiment(device)
121        except Exception as e: result['task_error']=str(e); result['task']=task_experiment('cpu')
122        result['device']=device
123    except Exception as e: result['task_error']=str(e)
124    Path('results.json').write_text(json.dumps(result,indent=2)); print(json.dumps(result,indent=2))
125
126if __name__=='__main__': main()