1import json, math, time
  2from pathlib import Path
  3import numpy as np
  4
  5SEED=7
  6rng=np.random.default_rng(SEED)
  7
  8# ---------- Core phase controller ----------
  9def control(theta,k):
 10    z=np.exp(1j*theta); r=z.mean()
 11    # formula in the prompt
 12    return (2*k/len(theta))*np.imag(z*np.conj(r))
 13
 14def grad(theta):
 15    return -control(theta,1.0)
 16
 17def phase_run(theta0,k=1.,dt=0.002,T=4.,delta=None,omega=None):
 18    th=theta0.copy(); n=int(T/dt); V=[]; us=[]; events=[]; held=np.zeros(len(th)); last=0
 19    if omega is None: omega=np.zeros(len(th))
 20    for q in range(n+1):
 21        ue=control(th,k)
 22        if delta is None:
 23            held=ue
 24        elif q==0 or np.linalg.norm(ue-held)>=delta:
 25            held=ue.copy(); events.append(q*dt)
 26        z=np.exp(1j*th); V.append(abs(z.mean())**2); us.append(ue.copy())
 27        if q<n: th += dt*(held+omega)
 28    return np.asarray(V),np.asarray(us),np.asarray(events),th
 29
 30def math_checks():
 31    N=24; theta=rng.normal(0,.22,N) # synchronized perturbation, clear descent
 32    dt=.0002; T=.25
 33    rows=[]
 34    for k in [.25,.5,1.,2.]:
 35        V,u,ev,_=phase_run(theta,k,dt,T,None)
 36        # finite difference at initial point compared with exact -k||grad||2
 37        observed=(V[1]-V[0])/dt
 38        predicted=-k*np.sum(grad(theta)**2)
 39        rows.append({'k':k,'initial_Vdot_observed':float(observed),'predicted':float(predicted),'ratio':float(observed/predicted)})
 40    # event tolerance: error and event rate versus delta
 41    evrows=[]
 42    theta2=rng.normal(0,.3,N); omega=np.linspace(-1.0,1.0,N)
 43    # finite-difference Lipschitz estimate for the exact control under uncontrolled drift
 44    _,ud,_,_=phase_run(theta2,1.,.001,10.,None,omega)
 45    M=float(np.percentile(np.linalg.norm(np.diff(ud,axis=0)/.001,axis=1),99))
 46    for d in [.005,.02,.08,.2]:
 47        V,u,ev,_=phase_run(theta2,1.,.001,10.,d,omega)
 48        dw=np.diff(ev)
 49        mindwell=float(dw.min()) if len(dw) else None
 50        # held/exact discrepancy reconstructed from held intervals, max just-before event
 51        # use trigger sampling; report event count and final V
 52        evrows.append({'delta':d,'events':int(len(ev)),'rate':float(len(ev)/2),'final_V':float(V[-1]),'min_dwell':mindwell,'predicted_dwell_lower_bound':float(d/M)})
 53    return rows,evrows
 54
 55# ---------- tiny neural experiment ----------
 56def train_model(mode, seed=11, epochs=18, force_cpu=False):
 57    import torch
 58    torch.manual_seed(seed); np.random.seed(seed)
 59    device='cpu' if force_cpu else ('cuda' if torch.cuda.is_available() else 'cpu')
 60    try:
 61        dev=torch.device(device)
 62        # sequence label is sign of a noisy temporal sum; common input encourages synchrony
 63        ntr,nte,L=640,256,24
 64        g=torch.Generator().manual_seed(seed)
 65        X=torch.randn(ntr+nte,L,1,device=dev)
 66        labels=(X.sum(1)[:,0]>0).long()
 67        # fixed split, same data for all modes
 68        tr=slice(0,ntr); te=slice(ntr,None)
 69        N=16; H=2*N; dt=.15; k=1.2; delta=.035
 70        W=torch.randn(H,H,device=dev)*.22
 71        U=torch.randn(1,H,device=dev)*.35
 72        b=torch.zeros(H,device=dev)
 73        W.requires_grad_(); U.requires_grad_(); b.requires_grad_()
 74        out=torch.randn(H,2,device=dev)*.15; out.requires_grad_()
 75        opt=torch.optim.Adam([W,U,b,out],lr=.018)
 76        event_total=0; steps_total=0; amp_drift=[]
 77        for ep in range(epochs):
 78            perm=torch.randperm(ntr,device=dev)
 79            for st in range(0,ntr,64):
 80                ids=perm[st:st+64]; xb=X[ids]; yb=labels[ids]; B=xb.shape[0]
 81                h=torch.zeros(B,H,device=dev); held=torch.zeros(B,N,device=dev)
 82                for t in range(L):
 83                    base=torch.tanh(h@W + xb[:,t]@U + b)
 84                    pairs=base.view(B,N,2)
 85                    z=pairs[...,0]+1j*pairs[...,1]; z=z/(torch.abs(z)+1e-6); r=z.mean(1,keepdim=True)
 86                    ue=(2*k/N)*torch.imag(z*torch.conj(r))
 87                    if mode=='continuous': held=ue
 88                    elif mode=='event':
 89                        trigger=torch.linalg.vector_norm(ue-held,dim=1)>=delta
 90                        held=torch.where(trigger[:,None],ue,held); event_total+=int(trigger.sum())
 91                    # exact tangent rotation, preserving each pair norm
 92                    ang=dt*held; x,yy=pairs[...,0],pairs[...,1]
 93                    rot=torch.stack([x*torch.cos(ang)-yy*torch.sin(ang),x*torch.sin(ang)+yy*torch.cos(ang)],-1)
 94                    h=rot.reshape(B,H); steps_total+=B
 95                logits=h@out; loss=torch.nn.functional.cross_entropy(logits,yb)
 96                opt.zero_grad(); loss.backward(); torch.nn.utils.clip_grad_norm_([W,U,b,out],1.); opt.step()
 97        with torch.no_grad():
 98            h=torch.zeros(nte,H,device=dev); held=torch.zeros(nte,N,device=dev); norm0=None
 99            for t in range(L):
100                pairs=torch.tanh(h@W+X[te,t]@U+b).view(nte,N,2)
101                if norm0 is None: norm0=torch.linalg.vector_norm(pairs,dim=-1)
102                z=pairs[...,0]+1j*pairs[...,1]; z=z/(torch.abs(z)+1e-6); r=z.mean(1,keepdim=True); ue=(2*k/N)*torch.imag(z*torch.conj(r))
103                if mode=='continuous': held=ue
104                elif mode=='event': held=torch.where((torch.linalg.vector_norm(ue-held,dim=1)>=delta)[:,None],ue,held)
105                ang=dt*held; x,yy=pairs[...,0],pairs[...,1]
106                h=torch.stack([x*torch.cos(ang)-yy*torch.sin(ang),x*torch.sin(ang)+yy*torch.cos(ang)],-1).reshape(nte,H)
107            acc=(h@out).argmax(1).eq(labels[te]).float().mean().item()
108            # compare base pair norms before/after rotation is exact up to float error
109            normerr=float((torch.linalg.vector_norm(h.view(nte,N,2),dim=-1)-torch.linalg.vector_norm(pairs,dim=-1)).abs().mean())
110        return {'accuracy':acc,'events_per_sequence':event_total/max(1,ntr*epochs),'rotation_norm_error':normerr,'device':str(dev)}
111    except Exception as e:
112        if device=='cuda':
113            return train_model(mode, seed, epochs, force_cpu=True)
114        return {'error':repr(e)}
115
116def main():
117    mathrows,eventrows=math_checks()
118    results={}
119    for mode in ['none','continuous','event']:
120        results[mode]=train_model(mode)
121    out={'seed':SEED,'math_gain_sweep':mathrows,'event_delta_sweep':eventrows,'neural_results':results}
122    Path('results.json').write_text(json.dumps(out,indent=2))
123    print(json.dumps(out,indent=2))
124if __name__=='__main__': main()