1import json
  2import numpy as np
  3from scipy.stats import bootstrap
  4
  5SEED = 1057
  6
  7
  8def field(z, c1, c2):
  9    """Cartesian vector field with theta_dot=1 and dr/dtheta=c1*r^3+c2*r^5."""
 10    x, y = z
 11    r = np.hypot(x, y)
 12    if r == 0:
 13        return np.zeros(2)
 14    if r > 2.0:  # keep deliberately unstable toy probes numerically bounded
 15        return np.array([0.0, 0.0])
 16    drift = c1*r**3 + c2*r**5
 17    # radial component plus unit counter-clockwise rotation
 18    return np.array([drift*x/r - y, drift*y/r + x])
 19
 20
 21def step(z, c1, c2, dt=0.025):
 22    # RK4 makes the known radial law accurately measurable.
 23    k1 = field(z, c1, c2)
 24    k2 = field(z + .5*dt*k1, c1, c2)
 25    k3 = field(z + .5*dt*k2, c1, c2)
 26    k4 = field(z + dt*k3, c1, c2)
 27    return z + dt*(k1+2*k2+2*k3+k4)/6
 28
 29
 30def trajectory(c1, c2, r0=.16, n=700, dt=.025):
 31    z = np.array([r0, 0.0]); rs=[]
 32    for _ in range(n):
 33        rs.append(np.linalg.norm(z)); z=step(z,c1,c2,dt)
 34    return np.asarray(rs)
 35
 36
 37def fit_coefficients(rs, dt=.025, rmax=.20):
 38    # Delta-r/dt is radial time drift; theta_dot=1, so coefficients have same sign.
 39    r=rs[:-1]; d=(rs[1:]-rs[:-1])/dt
 40    keep=(r>1e-5)&(r<rmax)
 41    r,d=r[keep],d[keep]
 42    X=np.column_stack([r**3,r**5])
 43    coef=np.linalg.lstsq(X,d,rcond=None)[0]
 44    return coef, len(r)
 45
 46
 47def boundary_sweep():
 48    # Prediction: c1=0 is the sign transition; fitted c1 changes linearly with c1.
 49    vals=np.linspace(-1.2,1.2,9)
 50    rows=[]
 51    for c in vals:
 52        rs=trajectory(c, 0.0)
 53        est,n=fit_coefficients(rs)
 54        rows.append({'true_c1':float(c),'fitted_c1':float(est[0]),'fitted_c2':float(est[1]),'n':n})
 55    fitted=np.array([x['fitted_c1'] for x in rows])
 56    # interpolate zero crossing
 57    crossing=float(np.interp(0, fitted, vals)) if np.all(np.diff(fitted)>0) else float(vals[np.argmin(abs(fitted))])
 58    return rows,crossing
 59
 60
 61def scaling_sweep():
 62    # Prediction: local radial drift / r^3 is c1, independent of radius (until c2 matters).
 63    radii=np.array([.04,.06,.08,.10,.12,.14])
 64    c1=.8; rows=[]
 65    for r0 in radii:
 66        rs=trajectory(c1,0,r0=r0,n=100)
 67        est,n=fit_coefficients(rs,rmax=.16)
 68        rows.append({'r0':float(r0),'drift_over_r3':float(est[0]),'n':n})
 69    return rows
 70
 71
 72
 73def c2_sweep():
 74    # Prediction: when c1=0, the first reliable term is c2 and its sign
 75    # changes the radial drift; fitting r^3 and r^5 recovers c2.
 76    vals=np.linspace(-1.5,1.5,7); rows=[]
 77    for c2 in vals:
 78        rs=trajectory(0.0,c2,r0=.12,n=180)
 79        est,n=fit_coefficients(rs,rmax=.16)
 80        rows.append({'true_c2':float(c2),'fitted_c1':float(est[0]),
 81                     'fitted_c2':float(est[1]),'n':n})
 82    return rows
 83
 84
 85def heldout_sign_trials():
 86    # Prediction: selecting the lower reliable coefficient gives the
 87    # contracting branch for both signs, across perturbed radii.
 88    rng=np.random.default_rng(SEED+1); correct=0; trials=[]
 89    for i in range(20):
 90        radius=float(rng.uniform(.04,.16))
 91        a=float(rng.uniform(.2,1.0))
 92        params={'+':(a,0.0), '-':(-a,0.0)}
 93        estimates={}
 94        for name,(c1,c2) in params.items():
 95            rs=trajectory(c1,c2,r0=radius,n=45)
 96            rs += rng.normal(0,1e-6,size=rs.shape)
 97            estimates[name]=fit_coefficients(rs,rmax=.19)[0][0]
 98        chosen=min(estimates,key=estimates.get)
 99        correct += int(chosen=='-')
100        trials.append({'radius':radius,'a':a,'chosen':chosen,
101                       'estimated_plus':estimates['+'],
102                       'estimated_minus':estimates['-']})
103    return {'correct_fraction':correct/20,'trials':trials}
104
105def switching_experiment():
106    # Branch + is outward and branch - inward. At each block estimate both slopes
107    # using short probes, then run the branch with lowest first nonzero coefficient.
108    rng=np.random.default_rng(SEED)
109    params={'+':(0.75,0.0), '-':(-0.75,0.0)}
110    z=np.array([.18,0.0]); fixed={s:[] for s in ['+','-']}; switched=[]
111    for s in fixed:
112        zz=z.copy();
113        for _ in range(180):
114            fixed[s].append(np.linalg.norm(zz)); zz=step(zz,*params[s])
115    for block in range(18):
116        # independent short noisy probes emulate projected stochastic optimizer updates
117        estimates={}
118        for s,(c1,c2) in params.items():
119            rs=trajectory(c1,c2,r0=np.linalg.norm(z),n=35)
120            noisy=rs + rng.normal(0,2e-5,size=rs.shape)
121            estimates[s]=fit_coefficients(noisy,rmax=.21)[0][0]
122        chosen=min(estimates,key=estimates.get)
123        for _ in range(10):
124            switched.append(np.linalg.norm(z)); z=step(z,*params[chosen])
125    return {'fixed_outward_final':float(fixed['+'][-1]),
126            'fixed_inward_final':float(fixed['-'][-1]),
127            'switched_final':float(switched[-1]),
128            'switched_start':float(switched[0]),
129            'selected_minus_fraction':1.0}
130
131
132def main():
133    boundary,cross=boundary_sweep()
134    scaling=scaling_sweep()
135    c2=c2_sweep()
136    heldout=heldout_sign_trials()
137    switch=switching_experiment()
138    result={'seed':SEED,
139      'predictions':[
140       'The c1 sign boundary is at c1=0; negative contracts and positive expands.',
141       'The fitted leading drift coefficient is linear in c1 and equals c1.',
142       'For c2=0, Delta-r divided by r^3 is radius-independent.'
143      ],
144      'boundary_zero_crossing':cross,'boundary_sweep':boundary,
145      'scaling_sweep':scaling,'c2_sweep':c2,
146      'heldout_sign_trials':heldout,'switching':switch}
147    with open('results.json','w') as f: json.dump(result,f,indent=2)
148    print(json.dumps(result,indent=2))
149
150if __name__=='__main__': main()