1import json, math, random
 2from pathlib import Path
 3import numpy as np
 4
 5SEED = 3127
 6rng = np.random.default_rng(SEED)
 7
 8# Degree-calibrated flow: f=-a*r^m*z + (I-zz'/||z||^2)h(z).
 9def tangent_project(z, h, eps=1e-12):
10    q = np.sum(z*z, axis=-1, keepdims=True)
11    return h - z * (np.sum(z*h, axis=-1, keepdims=True) / (q + eps))
12
13def calibrated_f(z, W, a, m):
14    r = 0.5*np.sum(z*z, axis=-1, keepdims=True)
15    h = np.tanh(z @ W.T)
16    return -0.5*a*(r**m)*z + tangent_project(z, h)
17
18def simulate(z0, W, a, m, dt, steps, calibrated=True):
19    z = z0.copy()
20    norms, radial_errors = [], []
21    for _ in range(steps):
22        r = 0.5*np.sum(z*z, axis=1, keepdims=True)
23        if calibrated:
24            f = calibrated_f(z, W, a, m)
25            allowed = -a*(r**(m+1))
26        else:
27            f = np.tanh(z @ W.T)
28            allowed = np.zeros_like(r)
29        radial = np.sum(z*f, axis=1, keepdims=True)
30        radial_errors.append(float(np.max(radial-allowed)))
31        norms.append(float(np.max(np.linalg.norm(z, axis=1))))
32        z = z + dt*f
33    norms.append(float(np.max(np.linalg.norm(z, axis=1))))
34    return np.asarray(norms), np.asarray(radial_errors)
35
36def slope(x, y):
37    return float(np.polyfit(np.log(x), np.log(np.maximum(y, 1e-30)), 1)[0])
38
39def main():
40    random.seed(SEED); np.random.seed(SEED)
41    d, n, a, steps = 8, 256, 1.0, 30000
42    # Small initial states avoid numerical stiffness for m>0 while retaining a clear tail.
43    z0 = rng.normal(size=(n,d)); z0 *= rng.uniform(.6,1.4,size=(n,1)) / np.linalg.norm(z0,axis=1,keepdims=True)
44    W = rng.normal(scale=.7/math.sqrt(d), size=(d,d))
45
46    decay = {}
47    for m in (0,1,2):
48        z = z0.copy(); rs=[]
49        dt=.01
50        for _ in range(steps):
51            rs.append(float(np.mean(.5*np.sum(z*z,axis=1))))
52            z += dt*calibrated_f(z,W,a,m)
53        rs=np.asarray(rs)
54        # Fit sufficiently late, but before float underflow / Euler floor.
55        lo, hi = (100, 1500) if m == 0 else (3000, 25000)
56        if m == 0:
57            fit = float(np.polyfit(np.arange(lo,hi)*dt, np.log(rs[lo:hi]), 1)[0])
58            expected = -a
59            observed = fit
60        else:
61            fit = slope(np.arange(lo,hi)*dt, rs[lo:hi])
62            expected = -1.0/m
63            observed = fit
64        decay[str(m)] = {"observed_log_r_slope": observed, "expected": expected,
65                         "initial_r": float(rs[0]), "final_r": float(rs[-1]),
66                         "max_abs_radial_identity_error": None}
67        # Exact derivative check, independent of Euler discretization.
68        zcheck = rng.normal(size=(1000,d)); f=calibrated_f(zcheck,W,a,m)
69        rr=.5*np.sum(zcheck*zcheck,axis=1,keepdims=True)
70        err=np.max(np.abs(np.sum(zcheck*f,axis=1,keepdims=True)+a*rr**(m+1)))
71        decay[str(m)]["max_abs_radial_identity_error"] = float(err)
72
73    # Stability control: same learned field and Euler step, with/without radial calibration.
74    stability={}
75    for dt in (0.01, 0.1, 0.5, 1.0):
76        base_n, base_e = simulate(z0,W,a,1,dt,300,False)
77        idea_n, idea_e = simulate(z0,W,a,1,dt,300,True)
78        stability[str(dt)]={
79            "baseline_final_max_norm":float(base_n[-1]),
80            "idea_final_max_norm":float(idea_n[-1]),
81            "baseline_peak_norm":float(np.max(base_n)),
82            "idea_peak_norm":float(np.max(idea_n)),
83            "idea_max_continuous_radial_violation":float(np.max(idea_e)),
84            "baseline_max_radial_derivative":float(np.max(base_e))
85        }
86    # After input removal, calibrated state shrinks while unconstrained residual state persists/grows.
87    result={"seed":SEED,"dimension":d,"batch":n,"a":a,"decay":decay,"stability":stability,
88            "note":"Continuous radial identity is tested exactly; trajectory tests use explicit Euler."}
89    Path("results.json").write_text(json.dumps(result,indent=2))
90    print(json.dumps(result,indent=2))
91
92if __name__ == '__main__': main()