1import json, math, os, random
  2import numpy as np
  3
  4SEED=1729
  5np.random.seed(SEED); random.seed(SEED)
  6
  7# ---------- Markov cocycle utilities ----------
  8def stationary(P):
  9    w,v=np.linalg.eig(P.T)
 10    i=np.argmin(abs(w-1)); p=np.real(v[:,i]); p=p/p.sum()
 11    return p
 12
 13def lyap(A,P,n=30000,seed=7):
 14    rng=np.random.default_rng(seed); k=len(A); d=A[0].shape[0]
 15    x=np.eye(d); smax=0.; smin=0.; z=rng.integers(k)
 16    p=stationary(P)
 17    for t in range(n):
 18        z=rng.choice(k,p=P[z]) if t else z
 19        x=A[z]@x
 20        if (t+1)%8==0:
 21            q,r=np.linalg.qr(x); diag=np.maximum(abs(np.diag(r)),1e-300)
 22            smax += math.log(diag.max()); smin += math.log(diag.min()); x=q
 23    return np.array([smax,smin])/(n)
 24
 25def spectral(M): return np.linalg.svd(M,compute_uv=False)[0]
 26
 27# Conformal zero-gap cocycle: scalar rotations have identical exponents exactly.
 28def rotation(th): return np.array([[math.cos(th),-math.sin(th)],[math.sin(th),math.cos(th)]])
 29P=np.array([[.82,.18],[.27,.73]])
 30r=np.array([.985,1.015]); th=np.array([.31,-.47])
 31A=np.array([r[i]*rotation(th[i]) for i in range(2)])
 32base=lyap(A,P)
 33# perturb mode 0 radially; d is exactly matrix spectral distance for this family
 34sweep=[]
 35for d in np.logspace(-5,-1,9):
 36    At=A.copy(); At[0]=(r[0]+d)*rotation(th[0])
 37    e=lyap(At,P, n=18000, seed=11)
 38    delta=abs(e[0]-base[0])
 39    sweep.append({'d':float(d),'delta_lambda':float(delta),'gap':float(e[0]-e[1]),
 40                  'linear_ratio':float(delta/d), 'invlog_half':float(1/math.sqrt(abs(math.log(d)))),
 41                  'invlog_one':float(1/abs(math.log(d)))})
 42# Estimate empirical constants from finite perturbations, then test held-out sizes.
 43C1=max(x['delta_lambda']*abs(math.log(x['d']))**.5 for x in sweep)
 44Cconf=max(x['delta_lambda']*abs(math.log(x['d'])) for x in sweep)
 45# directional random perturbations for a more meaningful conservative calibration
 46cal=[]
 47for j in range(12):
 48    d=10**np.random.uniform(-5,-2)
 49    H=np.random.randn(2,2); H=H/spectral(H)*d
 50    At=A.copy(); At[0]=A[0]+H
 51    e=lyap(At,P,n=12000,seed=j+30)
 52    cal.append(abs(e[0]-base[0])*abs(math.log(d)))
 53Cconf=max(Cconf, float(np.quantile(cal,.9)))
 54# Verify three predictions: zero gap, delta scales linearly locally, inverse-log envelope.
 55small=sweep[:5]
 56logd=np.log([x['d'] for x in small]); logdel=np.log(np.maximum([x['delta_lambda'] for x in small],1e-14))
 57slope=float(np.polyfit(logd,logdel,1)[0])
 58envelope=[x['delta_lambda'] <= Cconf/abs(math.log(x['d'])) for x in sweep]
 59# Trust-region acceptance sweep: proposal is accepted iff observed change is under modulus.
 60accept=[bool(x['delta_lambda'] <= Cconf/abs(math.log(x['d']))) for x in sweep]
 61
 62# ---------- Small RNN optimizer sanity experiment ----------
 63try:
 64 import torch
 65 torch.manual_seed(SEED); torch.set_num_threads(4)
 66 device='cuda' if torch.cuda.is_available() else 'cpu'
 67except Exception:
 68 torch=None; device='cpu'
 69
 70def train(control, seed):
 71 if torch is None: return {'loss':None,'acc':None,'backtracks':None,'device':'numpy'}
 72 try:
 73  torch.manual_seed(seed); np.random.seed(seed)
 74  dev=torch.device(device)
 75  n=512; T=18; h=24
 76  X=torch.randn(n,T,4,device=dev); y=(X[:,:T//2].sum((1,2))>0).long()
 77  class Net(torch.nn.Module):
 78   def __init__(self):
 79    super().__init__(); self.W=torch.nn.Parameter(torch.randn(h,h,device=dev)*.16)
 80    self.U=torch.nn.Parameter(torch.randn(h,4,device=dev)*.18); self.v=torch.nn.Parameter(torch.randn(h,device=dev)*.1); self.b=torch.nn.Parameter(torch.zeros(2,device=dev))
 81   def forward(self,x):
 82    z=torch.zeros(x.size(0),h,device=dev)
 83    for t in range(T): z=torch.tanh(z@self.W.T+x[:,t]@self.U.T)
 84    return torch.stack([z@self.v+self.b[0],-z@self.v+self.b[1]],1)
 85  net=Net(); opt=torch.optim.AdamW(net.parameters(),lr=.035,weight_decay=1e-4)
 86  oldW=net.W.detach().cpu().numpy().copy(); bt=0; accepted=0
 87  # Conservative C translated from exact toy calibration; stability budget is zero log spectral growth.
 88  for step in range(90):
 89   opt.zero_grad(); loss=torch.nn.functional.cross_entropy(net(X),y); loss.backward()
 90   torch.nn.utils.clip_grad_norm_(net.parameters(),5.0)
 91   saved={k:v.detach().clone() for k,v in net.state_dict().items()}
 92   opt.step()
 93   Wn=net.W.detach().cpu().numpy(); d=float(spectral(Wn-oldW))
 94   # For the recurrent matrix, use a two-mode proxy with identical matrices.
 95   oldexp=math.log(max(spectral(oldW),1e-8)); newexp=math.log(max(spectral(Wn),1e-8)); dl=abs(newexp-oldexp)
 96   bound=float(Cconf/max(abs(math.log(max(d,1e-8))),1e-8))
 97   bad=control and (dl>bound or newexp>.08)
 98   if bad:
 99    net.load_state_dict(saved); # undo and halve this optimizer's effective step by interpolation
100    for p in net.parameters():
101     if p.grad is not None: p.data.add_(p.grad,alpha=-.0175) # half lr SGD fallback
102    bt+=1
103   else: accepted+=1
104   oldW=net.W.detach().cpu().numpy().copy()
105  with torch.no_grad(): pred=net(X).argmax(1); acc=float((pred==y).float().mean().cpu())
106  return {'loss':float(loss.detach().cpu()),'acc':acc,'backtracks':bt,'accepted':accepted,'device':str(dev)}
107 except Exception as e:
108  return {'error':str(e),'device':'fallback'}
109
110runs=[]
111for s in [0,1,2]:
112 runs.append({'seed':s,'baseline':train(False,s),'idea':train(True,s)})
113
114def mean(key):
115 vals=[r[key].get('loss') for r in runs if r[key].get('loss') is not None]
116 return float(np.mean(vals)) if vals else None
117
118def meanacc(key):
119 vals=[r[key].get('acc') for r in runs if r[key].get('acc') is not None]
120 return float(np.mean(vals)) if vals else None
121out={'seed':SEED,'base_exponents':base.tolist(),'sweep':sweep,'predictions':{
122 'P1_zero_gap_max_abs_gap':float(max(abs(x['gap']) for x in sweep)),
123 'P2_local_loglog_slope_expected_1':slope,
124 'P3_inverse_log_envelope_fraction':float(np.mean(envelope)),
125 'P3_C_conformal':float(Cconf),'acceptance_by_d':accept},
126 'rnn_runs':runs,'rnn_summary':{'baseline_loss':mean('baseline'),'idea_loss':mean('idea'),'baseline_acc':meanacc('baseline'),'idea_acc':meanacc('idea')}}
127with open('results.json','w') as f: json.dump(out,f,indent=2)
128print(json.dumps(out,indent=2))