1import json, math
 2import numpy as np
 3
 4SEED=1349
 5alpha,beta=4.0,0.2
 6
 7def drift(x,g): return alpha*x*x+g-x**3-beta*x
 8def ddrift(x): return 2*alpha*x-3*x*x-beta
 9
10def extrema(): return np.sort(np.roots([3.,-2*alpha,beta]).real)
11
12def fold_bounds():
13    e=extrema(); q=lambda x: alpha*x*x-x**3-beta*x
14    # gamma values at which a pair of fixed points coalesces
15    return float(-q(e[1])),float(-q(e[0])),e
16
17def roots(g):
18    z=np.roots([1.,-alpha,beta,-g])
19    return np.sort(z[np.abs(z.imag)<1e-7].real)
20
21def integrate(gfun,omega,dt=0.004,cycles=6,noise=0.,x0=None):
22    T=2*np.pi/omega; n=min(max(int(cycles*T/dt),4000),300000)
23    dt=cycles*T/n; x=3.9 if x0 is None else x0
24    rng=np.random.default_rng(SEED+int(1000*omega))
25    xs=[]; gs=[]
26    for k in range(n):
27        t=k*dt
28        x=max(0.,x+dt*drift(x,gfun(t))+math.sqrt(2*noise*dt)*rng.normal())
29        if k>=n//2: xs.append(x); gs.append(gfun(t+dt))
30    xs,gs=np.asarray(xs),np.asarray(gs)
31    return float(abs(np.sum(.5*(xs[1:]+xs[:-1])*(gs[1:]-gs[:-1])))),float(dt)
32
33def frequency_sweep(g0,amp):
34    # Estimate tau on the high stable branch at the center forcing.
35    rr=roots(g0); stable=rr[ddrift(rr)<0]
36    tau=float(-1/ddrift(stable[-1]))
37    omegas=np.logspace(math.log10(.03/tau),math.log10(30/tau),13)
38    rows=[]
39    for w in omegas:
40        a,dt=integrate(lambda t:g0+amp*math.sin(w*t),w)
41        rows.append({'omega':float(w),'omega_tau':float(w*tau),'area':a})
42    return tau,rows
43
44def classification():
45    rng=np.random.default_rng(SEED); N,L=320,40
46    X=np.zeros((N,L)); y=np.arange(N)%2
47    for i in range(N): X[i]=(-5.0 if y[i]==0 else -3.9)+.08*rng.normal(size=L)
48    tr=np.arange(240); te=np.arange(240,320)
49    def run(kind):
50        st=[]
51        for seq in X:
52            if kind=='cell':
53                x=3.9
54                for u in seq:
55                    for _ in range(4): x=max(0.,x+.05*drift(x,u))
56            else:
57                x=0.
58                for u in seq: x=.9*x+.1*u
59            st.append(x)
60        st=np.asarray(st); th=(st[tr[y[tr]==0]].mean()+st[tr[y[tr]==1]].mean())/2
61        return float(np.mean((st[te]>th)==y[te])),float(st[y==0].mean()),float(st[y==1].mean())
62    return run('linear'),run('cell')
63
64def main():
65    lo,hi,e=fold_bounds(); gmid=(lo+hi)/2
66    rr=roots(gmid)
67    math_check={'extrema':e.tolist(),'predicted_fold_interval':[lo,hi],
68      'midpoint_gamma':gmid,'all_roots':rr.tolist(),
69      'derivatives':[float(ddrift(x)) for x in rr],
70      'stable_root_count':int(np.sum(ddrift(rr)<0))}
71    # Sweep spans both folds, so deterministic state switching is possible.
72    tau,freq=frequency_sweep(gmid,(hi-lo)*.56)
73    wstar=1/tau
74    noises=[0.,1e-4,5e-4,1e-3,5e-3,2e-2]
75    noise=[]
76    for d in noises:
77        a,_=integrate(lambda t:gmid+(hi-lo)*.56*math.sin(wstar*t),wstar,noise=d)
78        noise.append([d,a])
79    result={'parameters':{'alpha':alpha,'beta':beta,'seed':SEED},'math_check':math_check,
80      'tau_high':tau,'frequency_sweep':freq,'noise_sweep':noise,
81      'classification':{'linear':classification()[0],'cell':classification()[1]}}
82    with open('results.json','w') as f: json.dump(result,f,indent=2)
83    print(json.dumps(result,indent=2))
84if __name__=='__main__': main()