1import json, math, random, time
  2from pathlib import Path
  3import numpy as np
  4
  5SEED = 1552
  6np.random.seed(SEED)
  7random.seed(SEED)
  8
  9
 10def eig_sensitivity(A, B):
 11    """d lambda/dp for A(p)=A+pB, using the stated left/right formula."""
 12    vals, vr = np.linalg.eig(A)
 13    vals_l, vl = np.linalg.eig(A.T.conj())
 14    out = []
 15    for i, lam in enumerate(vals):
 16        # left eigenvector w satisfies w^H A=lambda w^H; column eigvec of A^H
 17        j = int(np.argmin(np.abs(vals_l - np.conj(lam))))
 18        v = vr[:, i]
 19        w = vl[:, j]
 20        d = (np.conj(w) @ B @ v) / (np.conj(w) @ v)
 21        out.append(d)
 22    return vals, np.asarray(out)
 23
 24
 25def radial_slope(lam, dlam):
 26    if abs(lam) < 1e-12:
 27        return abs(dlam)
 28    return float(np.real(np.exp(-1j*np.angle(lam))*dlam))
 29
 30
 31def toy_verification():
 32    # A(p)=A0+pB has a known dominant mode and a nontrivial off-diagonal coupling.
 33    # It makes the first-order crossing prediction directly falsifiable.
 34    rows = []
 35    for margin in [0.05, 0.10, 0.20, 0.35]:
 36        for gain in [0.25, 0.50, 1.0, 1.75]:
 37            a = 1.0 - margin
 38            A = np.array([[a, 0.65], [0.0, 0.42]], dtype=float)
 39            B = np.array([[gain, -0.3], [0.15, 0.05]], dtype=float)
 40            vals, ds = eig_sensitivity(A, B)
 41            i = int(np.argmax(np.abs(vals)))
 42            lam, dlam = vals[i], ds[i]
 43            slope = radial_slope(lam, dlam)
 44            pred = (1-abs(lam))/slope if slope > 1e-10 else np.inf
 45            # exact observed crossing by a dense one-dimensional sweep
 46            ps = np.linspace(0, min(2.0, max(.02, 2.5*pred)), 20001)
 47            radii = np.array([max(abs(np.linalg.eigvals(A+p*B))) for p in ps])
 48            ix = np.where(radii >= 1.0)[0]
 49            obs = float(ps[ix[0]]) if len(ix) else float('nan')
 50            # independently check derivative with centered finite differences
 51            eps = 1e-6
 52            fd = (np.linalg.eigvals(A+eps*B)[i]-np.linalg.eigvals(A-eps*B)[i])/(2*eps)
 53            relsens = abs(dlam-fd)/max(1e-12, abs(fd))
 54            rows.append(dict(margin=margin, gain=gain, predicted_radius=float(pred),
 55                             observed_radius=obs, relative_prediction_error=(abs(obs-pred)/pred if np.isfinite(obs) else None),
 56                             sensitivity=float(abs(dlam)), finite_difference_relative_error=float(relsens)))
 57    # Three explicit predictions: formula correctness, linear gain scaling, and margin scaling.
 58    formula_err = max(r['finite_difference_relative_error'] for r in rows)
 59    valid = [r for r in rows if r['observed_radius'] == r['observed_radius']]
 60    cross_err = float(np.median([r['relative_prediction_error'] for r in valid]))
 61    # Since derivative is gain-proportional, delta_hat * gain should be constant at fixed margin.
 62    scaled = []
 63    for m in [0.05, .10, .20, .35]:
 64        x = [r['predicted_radius']*r['gain'] for r in rows if r['margin']==m]
 65        scaled.append(float(np.std(x)/np.mean(x)))
 66    gain_scaling = float(np.median(scaled))
 67    # delta_hat is proportional to margin for this affine family.
 68    margin_ratios = []
 69    for g in [.25, .5, 1., 1.75]:
 70        x = [r['predicted_radius']/r['margin'] for r in rows if r['gain']==g]
 71        margin_ratios.append(float(np.std(x)/np.mean(x)))
 72    margin_scaling = float(np.median(margin_ratios))
 73    return rows, dict(max_formula_relative_error=formula_err,
 74                      median_crossing_relative_error=cross_err,
 75                      gain_scaling_cv=gain_scaling, margin_scaling_cv=margin_scaling,
 76                      predictions={'eigen_derivative_matches_finite_difference': formula_err < 1e-5,
 77                                   'crossing_radius_is_margin_over_slope': cross_err < .03,
 78                                   'radius_scales_inverse_linearly_with_gain': gain_scaling < .02,
 79                                   'radius_scales_linearly_with_margin': margin_scaling < .02})
 80
 81
 82def train_rnn(regularized, device='cpu', steps=350):
 83    import torch
 84    torch.manual_seed(SEED + int(regularized))
 85    D, T, B = 12, 25, 64
 86    Wh = torch.nn.Parameter(torch.randn(D,D,device=device)*0.28)
 87    Wx = torch.nn.Parameter(torch.randn(D,1,device=device)*0.35)
 88    Wo = torch.nn.Parameter(torch.randn(1,D,device=device)*0.15)
 89    b = torch.nn.Parameter(torch.zeros(D,device=device))
 90    params=[Wh,Wx,Wo,b]
 91    opt=torch.optim.Adam(params, lr=.012)
 92    history=[]
 93    def jacobian(h, u, gain=1.0):
 94        # Jacobian of one state update at one representative operating point.
 95        def fn(x): return torch.tanh(gain*Wh@x + Wx@u + b)
 96        return torch.autograd.functional.jacobian(fn, h, create_graph=True)
 97    for step in range(steps):
 98        u=torch.randn(B,T,1,device=device)
 99        # Stable teacher: exponentially filtered input, trained as sequence predictor.
100        y=torch.zeros(B,T,1,device=device)
101        for t in range(1,T): y[:,t]=.82*y[:,t-1]+.18*u[:,t]
102        h=torch.zeros(B,D,device=device); loss=0.
103        for t in range(T):
104            h=torch.tanh(h@Wh.T + u[:,t]@Wx.T + b)
105            loss=loss+torch.mean((h@Wo.T-y[:,t])**2)
106        loss=loss/T
107        reg=torch.tensor(0.,device=device)
108        if regularized and step % 8 == 0:
109            # The perturbation is recurrent gain p. Penalize predicted/actual local
110            # spectral radius at p=0 and p=+0.08, a differentiable proxy for margin.
111            h0=h[0].detach().requires_grad_(True); u0=u[0,-1].detach()
112            A=jacobian(h0,u0,1.0)
113            Ap=jacobian(h0,u0,1.08)
114            rho=torch.max(torch.abs(torch.linalg.eigvals(A)))
115            rhop=torch.max(torch.abs(torch.linalg.eigvals(Ap)))
116            # margin target and sensitivity-aware forward perturbation penalty
117            reg=0.8*torch.relu(rho-.88)**2 + 1.5*torch.relu(rhop-.93)**2
118        total=loss+reg
119        opt.zero_grad(); total.backward(); torch.nn.utils.clip_grad_norm_(params, 1.0); opt.step()
120        if step % 50 == 0: history.append(float(loss.detach().cpu()))
121    return [x.detach().cpu() for x in params], history
122
123
124def evaluate(params, device='cpu'):
125    import torch
126    Wh,Wx,Wo,b=params
127    Wh,Wx,Wo,b=[x.to(device) for x in params]
128    D=Wh.shape[0]
129    def run(gain, noise, T=180):
130        h=torch.zeros(1,D,device=device); norms=[]; outs=[]
131        u=torch.randn(1,T,1,device=device)*noise
132        for t in range(T):
133            h=torch.tanh(gain*(h@Wh.T)+u[:,t]@Wx.T+b)
134            norms.append(float(torch.linalg.vector_norm(h).item())); outs.append(float((h@Wo.T).item()))
135        return max(norms), float(np.mean(np.square(outs[-40:])))
136    # local eigenvalue and first-order prediction for gain perturbation, via numpy fd
137    h=torch.zeros(D,device=device,requires_grad=True); u=torch.zeros(1,device=device)
138    def fn(x,g): return torch.tanh(g*(Wh@x)+Wx@u+b)
139    A=torch.autograd.functional.jacobian(lambda x:fn(x,1.),h).detach().cpu().numpy()
140    Ap=torch.autograd.functional.jacobian(lambda x:fn(x,1.001),h).detach().cpu().numpy()
141    vals=np.linalg.eigvals(A); rho=max(abs(vals)); dr=(max(abs(np.linalg.eigvals(Ap)))-rho)/.001
142    pred=(1-rho)/max(1e-9,dr)
143    return dict(rho=float(rho), predicted_radius=float(pred), gain1_noise1=run(1.,1.), gain1p_noise1=run(1.08,1.), gain1_noise2=run(1.,2.))
144
145
146def main():
147    rows,toy=toy_verification()
148    result={'seed':SEED,'toy_summary':toy,'toy_rows':rows}
149    try:
150        import torch
151        device='cuda' if torch.cuda.is_available() else 'cpu'
152        try:
153            base,hb=train_rnn(False,device); idea,hi=train_rnn(True,device)
154        except Exception as e:
155            device='cpu'; base,hb=train_rnn(False,device); idea,hi=train_rnn(True,device)
156        result.update({'device':device,'baseline':evaluate(base,device),'idea':evaluate(idea,device),
157                       'training_loss_samples':{'baseline':hb,'idea':hi}})
158    except Exception as e:
159        result['rnn_error']=repr(e)
160    Path('results.json').write_text(json.dumps(result,indent=2))
161    print(json.dumps(result,indent=2))
162
163if __name__=='__main__': main()