1import json
 2import numpy as np
 3from scipy.optimize import brentq
 4
 5
 6def transition(g, taps):
 7    taps=np.asarray(taps,float); T=np.zeros((1+len(taps),1+len(taps)))
 8    T[0,0]=1.; T[0,1:]=g*taps; T[1:,:-1]=np.eye(len(taps)); return T
 9
10
11def margin(g,taps,sigma=1,window=None):
12    if not np.isfinite(g): return float('nan')
13    taps=np.asarray(taps,float); window=len(taps)+1 if window is None else window
14    R=np.linalg.matrix_power(transition(g,taps),window)-sigma*np.eye(len(taps)+1)
15    return float(np.linalg.svd(R,compute_uv=False)[-1])
16
17
18def determinant(g,taps,sigma):
19    T=transition(g,taps); n=T.shape[0]
20    return float(np.linalg.det(np.linalg.matrix_power(T,n)-sigma*np.eye(n)))
21
22
23def roots_in_range(taps,sigma,lo=-8.,hi=8.):
24    xs=np.linspace(lo,hi,32001); ys=np.array([determinant(x,taps,sigma) for x in xs]); roots=[]
25    for i in range(len(xs)-1):
26        if not (np.isfinite(ys[i]) and np.isfinite(ys[i+1])): continue
27        if ys[i]==0: roots.append(xs[i])
28        elif ys[i]*ys[i+1]<0:
29            roots.append(brentq(lambda z:determinant(z,taps,sigma),xs[i],xs[i+1]))
30    return sorted(set(round(float(r),12) for r in roots))
31
32
33def first_nontrivial_root(taps,sigma):
34    roots=[r for r in roots_in_range(taps,sigma) if abs(r)>1e-4]
35    return (min(roots,key=abs),'determinant_root') if roots else (float('nan'),'none')
36
37
38def grid_minimum(taps,sigma,lo=-8.,hi=8.):
39    xs=np.linspace(lo,hi,32001); ys=np.array([margin(x,taps,sigma) for x in xs])
40    i=int(np.nanargmin(ys)); return float(xs[i]),float(ys[i])
41
42
43def numerical_check():
44    base=np.array([.7,-.25,.12]); rp,_=first_nontrivial_root(base,1); rm,km=first_nontrivial_root(base,-1)
45    rows=[]
46    for a in [.5,1.,2.]:
47        r,method=first_nontrivial_root(a*base,1)
48        rows.append({'tap_scale':a,'predicted_signed_gain':rp/a,'observed_signed_gain':r,'relative_error':None if not np.isfinite(r) else abs(r-rp/a)/(abs(rp/a)+1e-12),'method':method})
49    am,amin=grid_minimum(base,-1)
50    return {'predictions':['Fixed tap ratios imply roots scale as 1/a because Phi(g,a*taps)=Phi(a*g,taps).','A true periodic/antiperiodic bifurcation has det(Phi-sigma I)=0 and margin s_min=0.','The sigma sign selects distinct modes; absence of a root is itself a negative result.'],'base_periodic_gain':rp,'base_antiperiodic_gain':None if not np.isfinite(rm) else rm,'periodic_margin_at_root':margin(rp,base,1),'antiperiodic_root_margin':None if not np.isfinite(rm) else margin(rm,base,-1),'antiperiodic_grid_min_gain':am,'antiperiodic_grid_min_margin':amin,'scaling_rows':rows,'periodic_roots':roots_in_range(base,1),'antiperiodic_roots':roots_in_range(base,-1)}
51
52
53def train_comparison(seed=7):
54    import torch
55    torch.manual_seed(seed); np.random.seed(seed); device='cuda' if torch.cuda.is_available() else 'cpu'
56    try:
57        n,L,delay=96,24,3; X=torch.randn(n,L,1,device=device); Y=torch.zeros_like(X); Y[:,delay:]=X[:,:-delay]
58        class Delayed(torch.nn.Module):
59            def __init__(self):
60                super().__init__(); self.inp=torch.nn.Linear(1,8); self.out=torch.nn.Linear(8,1); self.g=torch.nn.Parameter(torch.tensor(.65)); self.taps=torch.nn.Parameter(torch.tensor([.55,-.25,.10]))
61            def forward(self,x):
62                h=[torch.zeros(x.shape[0],8,device=x.device) for _ in range(4)]; ys=[]
63                for t in range(L):
64                    cur=torch.tanh(self.inp(x[:,t])+self.g*(self.taps[0]*h[-1]+self.taps[1]*h[-2]+self.taps[2]*h[-3])); h.append(cur); h=h[-4:]; ys.append(self.out(cur))
65                return torch.stack(ys,1)
66            def penalty(self):
67                T=torch.zeros((4,4),device=device); T[0,0]=1.; T[0,1:]=self.g*self.taps; T[1:,:-1]=torch.eye(3,device=device); P=torch.matrix_power(T,4)
68                return sum(torch.relu(.08-torch.linalg.svdvals(P-s*torch.eye(4,device=device))[-1])**2 for s in (1.,-1.))
69        results={}
70        for name,monitored in [('baseline',False),('idea',True)]:
71            model=Delayed().to(device); opt=torch.optim.Adam(model.parameters(),lr=.025); obj=0.
72            for _ in range(180):
73                task=torch.mean((model(X)-Y)**2); loss=task+(.8*model.penalty() if monitored else 0.); opt.zero_grad(); loss.backward(); torch.nn.utils.clip_grad_norm_(model.parameters(),2.); opt.step(); obj=float(loss.detach().cpu())
74            with torch.no_grad(): task=float(torch.mean((model(X)-Y)**2).cpu()); g=float(model.g.cpu()); t=model.taps.detach().cpu().numpy()
75            results[name]={'task_mse':task,'final_objective':obj,'gain':g,'taps':t.tolist(),'periodic_margin':margin(g,t,1),'antiperiodic_margin':margin(g,t,-1)}
76        return {'device':device,'results':results}
77    except Exception as e: return {'device':'cpu-fallback','error':repr(e)}
78
79if __name__=='__main__':
80    out={'math_check':numerical_check(),'training':train_comparison()}
81    with open('results.json','w') as f: json.dump(out,f,indent=2)
82    print(json.dumps(out,indent=2))