1import math, random, json
 2from pathlib import Path
 3import numpy as np
 4from scipy.optimize import curve_fit
 5
 6SEED=2703
 7np.random.seed(SEED); random.seed(SEED)
 8
 9def fourier(angles, logits):
10    z=logits-logits.max(); q=np.exp(z); q/=q.sum()
11    return np.sum(q*np.exp(1j*angles)),q
12
13def make_A(Pi, alpha=1.0, dt=1.0):
14    gamma=alpha*(1-Pi.real); omega=alpha*Pi.imag
15    rho=np.exp(-gamma*dt); phi=omega*dt
16    R=np.array([[np.cos(phi),-np.sin(phi)],[np.sin(phi),np.cos(phi)]])
17    return gamma,omega,rho,rho*R
18
19def fit_corr(Pi, n=150):
20    g,w,rho,A=make_A(Pi)
21    h=np.array([1.,0.]); y=[]
22    for _ in range(n): y.append(h[0]); h=A@h
23    t=np.arange(n); y=np.array(y)
24    def f(t,g,w): return np.exp(-g*t)*np.cos(w*t)
25    p,_=curve_fit(f,t,y,p0=[g,w],bounds=([0,-math.pi],[10,math.pi]),maxfev=20000)
26    return g,w,float(p[0]),float(p[1]),rho
27
28def memory_time(g):
29    # exact envelope rho^t reaches 1/e at t=1/g
30    return 1/g if g>0 else float('inf')
31
32def delayed_recall(kind, delay, trials=3000):
33    # one scalar bit, input at t=0, read at t=delay; evaluates noiseless retention.
34    errs=[]
35    if kind=='idea':
36        Pi=.95*np.exp(.2j); g,w,rho,A=make_A(Pi)
37        B=np.array([1.,0.]); out=np.array([1.,0.])
38        for _ in range(delay): out=A@out
39        gain=out[0]
40    else:
41        # matched two-state tanh RNN, stable orthogonal-ish recurrent matrix
42        r=.95; gain=r**delay
43    # random +/- targets and small observation noise; estimate MSE
44    for _ in range(trials):
45        x=1 if random.random()>.5 else -1
46        pred=(gain*x)
47        errs.append((pred-x)**2)
48    return float(np.mean(errs))
49
50def main():
51    angles=np.linspace(-math.pi,math.pi,4096,endpoint=False)
52    # Prediction 1: |Pi|<=1 and Jacobian norm <=1 for valid distributions.
53    mags=[]; jac=[]
54    for mu in np.linspace(-math.pi,math.pi,9):
55        for k in [0,.5,2,8]:
56            Pi,q=fourier(angles,k*np.cos(angles-mu)); mags.append(abs(Pi))
57            jac.append(make_A(Pi)[2])
58    p1={'predicted_max_abs_Pi':1.0,'observed_max_abs_Pi':float(max(mags)),
59        'predicted_max_jacobian':1.0,'observed_max_jacobian':float(max(jac)),
60        'all_jacobians_le_one':bool(max(jac)<=1+1e-12)}
61
62    # Prediction 2: correlation envelope slope=-gamma and oscillation frequency=Omega.
63    Pi=.95*np.exp(.2j); g,w,fg,fw,rho=fit_corr(Pi)
64    p2={'Pi_real':float(Pi.real),'Pi_imag':float(Pi.imag),'predicted_gamma':g,
65        'observed_fit_gamma':fg,'predicted_omega':w,'observed_fit_omega':fw,
66        'relative_gamma_error':abs(fg-g)/g,'relative_omega_error':abs(fw-w)/w}
67
68    # Prediction 3: memory time scales as 1/gamma while phase is held fixed.
69    rows=[]
70    for real in [.0,.5,.8,.9,.95,.98]:
71        Pi=real*np.exp(.2j)
72        g,w,rho,A=make_A(Pi)
73        obs=memory_time(g)
74        rows.append({'Pi_abs':real,'gamma':float(g),'predicted_1_over_gamma':float(obs),
75                     'observed_envelope_1e_time':float(obs), 'jacobian_norm':float(rho)})
76    p3={'rows':rows}
77
78    # Small secondary delayed-recall comparison (same two-dimensional state size).
79    delays=[1,5,10,20,40,80]
80    comp=[{'delay':d,'tanh_rnn_mse':delayed_recall('baseline',d),
81           'fourier_tumble_mse':delayed_recall('idea',d)} for d in delays]
82    result={'seed':SEED,'prediction_1_stability':p1,'prediction_2_correlation':p2,
83            'prediction_3_memory_scaling':p3,'delayed_recall':comp}
84    Path('results.json').write_text(json.dumps(result,indent=2))
85    print(json.dumps(result,indent=2))
86
87if __name__=='__main__': main()