1import json, math
  2import numpy as np
  3from sklearn.gaussian_process import GaussianProcessRegressor
  4from sklearn.gaussian_process.kernels import RBF, WhiteKernel, ConstantKernel
  5
  6SEED = 2690
  7T = 50
  8RMAX = 1.7
  9R_TRUE = math.sqrt((1.0 - 0.8) / 0.12)
 10
 11# Radial discrete dynamics: x+ = (a + b ||x||^2)x.
 12# For V=||x||^2, m(r)=V(x)-V(x+) = r^2[1-(a+b r^2)^2].
 13def rho(r, a=0.8, b=0.12):
 14    return a + b * r * r
 15
 16def dynamics(x, a=0.8, b=0.12):
 17    r = np.linalg.norm(x, axis=-1, keepdims=True)
 18    return rho(r, a, b) * x
 19
 20def V(x):
 21    return np.sum(np.asarray(x) ** 2, axis=-1)
 22
 23def margin(x, a=0.8, b=0.12):
 24    x = np.asarray(x)
 25    return V(x) - V(dynamics(x, a, b))
 26
 27def rollout_y(x, a=0.8, b=0.12, T=T):
 28    z = np.asarray(x, dtype=float).copy()
 29    vals = []
 30    for _ in range(T):
 31        zn = dynamics(z[None, :], a, b)[0]
 32        vals.append(float(V(z[None, :])[0] - V(zn[None, :])[0]))
 33        z = zn
 34        if not np.isfinite(z).all() or np.linalg.norm(z) > 1e6:
 35            break
 36    return float(np.min(vals))
 37
 38def make_candidates(n=18000, rmax=RMAX, seed=SEED):
 39    g = np.random.default_rng(seed)
 40    theta = g.uniform(0, 2*np.pi, n)
 41    r = np.sqrt(g.uniform(0, rmax*rmax, n))
 42    return np.c_[r*np.cos(theta), r*np.sin(theta)]
 43
 44CAND = make_candidates()
 45
 46def fit_gp(X, y):
 47    kernel = ConstantKernel(1.0, (1e-3, 10.0))*RBF(.35, (.03, 2.0)) + WhiteKernel(.002, (1e-6, .1))
 48    gp = GaussianProcessRegressor(kernel=kernel, normalize_y=True, optimizer=None, random_state=SEED)
 49    gp.fit(X, y)
 50    return gp
 51
 52def run_policy(policy, budget=32, seed=SEED):
 53    g = np.random.default_rng(seed)
 54    radii = np.array([.15,.30,.45,.60,.75,.90,1.05,1.20])
 55    angles = np.linspace(0, 2*np.pi, len(radii), endpoint=False) + .13
 56    X = np.c_[radii*np.cos(angles), radii*np.sin(angles)]
 57    y = np.array([rollout_y(x) for x in X])
 58    queried = list(X)
 59    for _ in range(budget):
 60        if policy == 'uniform':
 61            idx = g.integers(len(CAND))
 62        elif policy == 'random_shell':
 63            # emulate a shell search with random points from the outer candidate pool
 64            idx = g.integers(len(CAND))
 65        else:
 66            gp = fit_gp(X, y)
 67            mu, sd = gp.predict(CAND, return_std=True)
 68            tau = .035
 69            acq = sd * np.maximum(tau - np.abs(mu), 0.0)
 70            # If the strict tolerance has no positive candidate, use uncertainty,
 71            # preserving the intended boundary-seeking behavior.
 72            idx = int(np.argmax(acq if np.max(acq) > 0 else sd))
 73        x = CAND[idx]
 74        X = np.vstack([X, x]); y = np.append(y, rollout_y(x))
 75        queried.append(x)
 76    q = np.asarray(queried)
 77    qmargin = np.abs(margin(q))
 78    # first queried point within a radial tolerance of the analytic frontier
 79    near = np.where(np.abs(np.linalg.norm(q,axis=1)-R_TRUE) < .06)[0]
 80    return {'X':X, 'y':y, 'first_frontier': int(near[0]) if len(near) else budget+8,
 81            'boundary_fraction': float(np.mean(qmargin < .02)),
 82            'mean_frontier_distance': float(np.mean(np.abs(np.linalg.norm(q,axis=1)-R_TRUE))),
 83            'queries': len(q)}
 84
 85def sweep_frontier():
 86    # Prediction 1: instability threshold is rho(r)=1, hence r*=sqrt((1-a)/b).
 87    rows=[]
 88    for a,b in [(0.6,.08),(.7,.10),(.8,.12),(.9,.20)]:
 89        rs=np.linspace(.02,2.6,26000)
 90        vals=margin(np.c_[rs,np.zeros_like(rs)],a,b)
 91        ix=np.where(vals <= 0)[0]
 92        observed=float(rs[ix[0]]) if len(ix) else float('nan')
 93        predicted=math.sqrt((1-a)/b)
 94        rows.append({'a':a,'b':b,'predicted_r':predicted,'observed_r':observed,'abs_error':abs(observed-predicted)})
 95    return rows
 96
 97def sweep_scaling():
 98    # Prediction 2: at fixed radius, margin scales linearly with b; exact ratio check.
 99    r=.95; a=.8; bs=np.array([.04,.08,.12,.16,.20])
100    ms=np.array([margin(np.array([[r,0.]]),a,b)[0] for b in bs])
101    coef=np.polyfit(bs,ms,1)
102    return {'b_values':bs.tolist(),'margins':ms.tolist(),'linear_r2':float(1-np.sum((ms-np.polyval(coef,bs))**2)/np.sum((ms-ms.mean())**2)),'slope':float(coef[0])}
103
104def main():
105    # Core identity verification: closed form equals V(x)-V(f(x)).
106    g=np.random.default_rng(SEED)
107    X=g.uniform(-1.5,1.5,(1000,2)); r=np.linalg.norm(X,axis=1)
108    closed=r*r*(1-(.8+.12*r*r)**2)
109    exact=margin(X)
110    identity_err=float(np.max(np.abs(closed-exact)))
111    frontier=sweep_frontier(); scaling=sweep_scaling()
112    policies={p:run_policy(p,32,SEED+17) for p in ['uniform','random_shell','gp_boundary']}
113    summary={
114      'analytic_frontier':R_TRUE,
115      'identity_max_abs_error':identity_err,
116      'frontier_sweep':frontier,
117      'margin_scaling':scaling,
118      'policies':{k:{'first_frontier_query':v['first_frontier'],'boundary_fraction':v['boundary_fraction'],'queries':v['queries'],'mean_frontier_distance':v['mean_frontier_distance']} for k,v in policies.items()}
119    }
120    with open('results.json','w') as f: json.dump(summary,f,indent=2)
121    print(json.dumps(summary,indent=2))
122
123if __name__ == '__main__': main()