Pullback random-attractor monitor / pullback_monitor.py
Failed on benchmark
1"""Pullback random-attractor monitor for a linear noisy damped oscillator.
2
3The shared additive noise cancels between pullback replicas, so their diameter
4measures contraction of the stochastic one-step map rather than noise variance.
5"""
6import json
7from pathlib import Path
8import numpy as np
9
10
11def oscillator_matrix(h, omega, damping, method):
12 """Return the 2x2 deterministic state transition matrix."""
13 # x' = v, v' = -omega^2*x - damping*v + noise
14 if method == "euler":
15 return np.array([[1.0, h], [-h * omega**2, 1.0 - h * damping]])
16 if method == "semi_implicit":
17 # Backward Euler only on velocity damping/force, then symplectic x update.
18 den = 1.0 + h * damping
19 v_x = -h * omega**2 / den
20 v_v = 1.0 / den
21 return np.array([[1.0 + h * v_x, h * v_v], [v_x, v_v]])
22 raise ValueError(method)
23
24
25def predicted_rho(M):
26 return float(np.max(np.abs(np.linalg.eigvals(M))))
27
28
29def monitor(h, method, omega, damping, K=32, steps=128, seed=0, sigma=0.15):
30 rng = np.random.default_rng(seed)
31 # Fixed bounded ball: random directions with radius in [0.8, 1.2].
32 z = rng.normal(size=(K, 2))
33 z /= np.linalg.norm(z, axis=1, keepdims=True)
34 z *= rng.uniform(0.8, 1.2, size=(K, 1))
35 M = oscillator_matrix(h, omega, damping, method)
36 diam = [pairwise_rms(z)]
37 # One identical noise vector is replayed for every replica at every step.
38 noises = rng.normal(size=(steps, 2)) * sigma * np.sqrt(h)
39 for noise in noises:
40 z = z @ M.T + noise[None, :]
41 diam.append(pairwise_rms(z))
42 return np.asarray(diam), M
43
44
45def pairwise_rms(z):
46 # RMS of all pair distances, a stable diameter proxy.
47 d = z[:, None, :] - z[None, :, :]
48 return float(np.sqrt(np.mean(np.sum(d * d, axis=-1))))
49
50
51def fit_slope(diam, h, start=8, end=None):
52 end = len(diam) if end is None else end
53 y = np.log(np.maximum(diam[start:end], 1e-30))
54 x = np.arange(start, end) * h
55 return float(np.polyfit(x, y, 1)[0])
56
57
58def independent_growth_rate(M, h, steps=300, seed=991):
59 rng = np.random.default_rng(seed)
60 v = rng.normal(size=2); v /= np.linalg.norm(v)
61 logs = []
62 for _ in range(steps):
63 v = M @ v
64 n = np.linalg.norm(v)
65 logs.append(np.log(n))
66 v /= n
67 # This is the accumulated Lyapunov sum for repeatedly renormalized vectors.
68 return float(np.mean(np.asarray(logs)[50:]) / h)
69
70
71def stability_sweep(omega=1.0, damping=0.2):
72 # Prediction: Euler is stable iff rho(M)<1, with boundary h= damping/omega^2.
73 hs = [0.10, 0.18, 0.20, 0.22, 0.40]
74 rows = []
75 for method in ("euler", "semi_implicit"):
76 for h in hs:
77 d, M = monitor(h, method, omega, damping, steps=128, seed=10)
78 rho = predicted_rho(M)
79 slope = fit_slope(d, h, start=10, end=100)
80 rows.append({"method": method, "h": h, "predicted_rho": rho,
81 "measured_slope": slope, "D0": d[0], "D128": d[-1],
82 "contracting_predicted": bool(rho < 1.0),
83 "contracting_observed": bool(d[-1] < d[0])})
84 return rows
85
86
87def rate_sweep(omega=1.0, h=0.05):
88 # Prediction: as damping varies, lambda_h=log(rho(M))/h and measured log D
89 # slope track it; near the boundary contraction time is ~1/|lambda_h|.
90 rows = []
91 for damping in [0.05, 0.10, 0.20, 0.40, 0.80]:
92 d, M = monitor(h, "euler", omega, damping, steps=500, seed=20)
93 rho = predicted_rho(M)
94 theory = np.log(rho) / h
95 measured = fit_slope(d, h, start=12, end=150)
96 independent = independent_growth_rate(M, h)
97 target = d[0] * 0.1
98 hit = np.where(d <= target)[0]
99 hit_time = float(hit[0] * h) if len(hit) else None
100 predicted_time = float(np.log(10.0) / (-theory)) if theory < 0 else None
101 rows.append({"damping": damping, "predicted_lambda": theory,
102 "measured_pullback_slope": measured,
103 "independent_vector_slope": independent,
104 "predicted_time_to_0.1": predicted_time,
105 "observed_time_to_0.1": hit_time})
106 return rows
107
108
109def noise_control(h=0.05, damping=0.4, omega=1.0, K=32, steps=600):
110 """Compare the monitor's shared-noise replay with independent noise.
111 Prediction: shared noise contracts; independent noise reaches a noise floor.
112 """
113 rng = np.random.default_rng(123)
114 z = rng.normal(size=(K, 2)); z /= np.linalg.norm(z, axis=1, keepdims=True)
115 z *= rng.uniform(0.8, 1.2, size=(K, 1))
116 M = oscillator_matrix(h, omega, damping, "euler")
117 shared = z.copy(); independent = z.copy()
118 shared_d, independent_d = [pairwise_rms(shared)], [pairwise_rms(independent)]
119 for _ in range(steps):
120 eps = rng.normal(size=(K, 2)) * 0.15 * np.sqrt(h)
121 shared = shared @ M.T + eps[0][None, :]
122 independent = independent @ M.T + eps
123 shared_d.append(pairwise_rms(shared)); independent_d.append(pairwise_rms(independent))
124 return {"shared_initial": shared_d[0], "shared_final": shared_d[-1],
125 "independent_initial": independent_d[0], "independent_final": independent_d[-1],
126 "shared_slope": fit_slope(np.asarray(shared_d), h, 10, 90),
127 "independent_slope": fit_slope(np.asarray(independent_d), h, 10, 90)}
128
129def main():
130 out = {"stability_boundary_prediction": "Euler rho=1 at h=damping/omega^2=0.20; semi-implicit should remain stable here",
131 "rate_scaling_prediction": "pullback log-diameter slope and independent vector growth equal log(rho)/h; time scales as 1/abs(lambda)",
132 "noise_replay_prediction": "shared replay contracts while independent noise produces a nonzero diameter floor",
133 "noise_control": noise_control(),
134 "stability_sweep": stability_sweep(), "rate_sweep": rate_sweep()}
135 Path("results.json").write_text(json.dumps(out, indent=2))
136 print(json.dumps(out, indent=2))
137
138if __name__ == "__main__":
139 main()