Dynamics-Matched Contractive Reservoir / experiment.py
Failed on benchmark
1import json, math
2from pathlib import Path
3import numpy as np
4
5SEED = 1487
6rng = np.random.default_rng(SEED)
7
8
9def orthogonal(n, seed):
10 q, _ = np.linalg.qr(np.random.default_rng(seed).normal(size=(n, n)))
11 return q
12
13
14def boundary_sweep():
15 # F(x)=gamma*A*x, ||A||_2=Lambda. Prediction: decay iff gamma*Lambda<1.
16 n, steps, lam = 8, 40, 1.2
17 Q = orthogonal(n, 1)
18 x0 = np.ones(n) / np.sqrt(n)
19 rows = []
20 for gamma in [0.60, 0.75, 0.80, 0.833333, 0.90, 1.00]:
21 F = gamma * lam * Q
22 x = x0.copy()
23 norms = []
24 for _ in range(steps):
25 x = F @ x
26 norms.append(np.linalg.norm(x))
27 observed = np.mean(np.array(norms[-5:]) / np.array(norms[-6:-1]))
28 predicted = gamma * lam
29 rows.append({"gamma": gamma, "predicted_rho": predicted,
30 "observed_ratio": float(observed),
31 "classification_pred": "decay" if predicted < 1 else "grow_or_neutral",
32 "classification_obs": "decay" if observed < 0.999 else "grow_or_neutral"})
33 return rows
34
35
36def decay_sweep():
37 # F=rho*Q is an isometry scaled by rho. Prediction log(norm_k/norm_0)/k=log(rho).
38 n, steps = 10, 30
39 Q = orthogonal(n, 2)
40 x0 = np.ones(n) / np.sqrt(n)
41 rows = []
42 for rho in [0.50, 0.70, 0.90, 0.97]:
43 x = x0.copy()
44 for _ in range(steps):
45 x = rho * Q @ x
46 measured = math.log(np.linalg.norm(x) / np.linalg.norm(x0)) / steps
47 rows.append({"rho": rho, "predicted_log_slope": math.log(rho),
48 "measured_log_slope": float(measured),
49 "abs_error": float(abs(measured - math.log(rho)))})
50 return rows
51
52
53def idm_acc(v, gap, dv):
54 # Safe numerical IDM implementation; dv=v-v_leader.
55 amax, v0, delta, s0, T, b = 1.2, 1.4, 4.0, 0.5, 1.1, 1.5
56 gap = max(float(gap), 0.15)
57 sstar = s0 + v*T + v*dv/(2*np.sqrt(amax*b))
58 return amax * (1 - (max(v, 0)/v0)**delta - (sstar/gap)**2)
59
60
61def make_data(T=700, n=8):
62 # Ring traffic-like system: position and speed, only every other position observed.
63 dt = 0.08
64 pos = np.zeros((T, n)); vel = np.zeros((T, n))
65 pos[0] = np.linspace(0, 8, n, endpoint=False)
66 vel[0] = 0.8 + 0.05*rng.normal(size=n)
67 for t in range(T-1):
68 for i in range(n):
69 lead = (i + 1) % n
70 gap = (pos[t, lead] - pos[t, i]) % 8.0
71 dv = vel[t, i] - vel[t, lead]
72 acc = idm_acc(vel[t,i], gap, dv)
73 vel[t+1,i] = np.clip(vel[t,i] + dt*acc + 0.008*rng.normal(), -0.2, 2.0)
74 pos[t+1,i] = (pos[t,i] + dt*vel[t+1,i]) % 8.0
75 # normalize physical channels for stable reservoir driving
76 y = np.concatenate([pos/8.0, vel/2.0], axis=1)
77 obs_idx = np.arange(0, 2*n, 2) # positions of half the cars only
78 u = y[:, obs_idx]
79 return u, y
80
81
82def graph_matrix(n, scale=0.88):
83 # Mechanistic directed predecessor coupling, normalized to prescribed norm.
84 M = np.zeros((2*n, 2*n))
85 for i in range(n):
86 lead = (i+1) % n
87 # position error and relative velocity drive each vehicle's two channels
88 M[n+i, i] = 0.20
89 M[n+i, lead] = -0.20
90 M[n+i, n+i] = 0.72
91 M[n+i, n+lead] = 0.12
92 M[i, i] = 0.95
93 M[i, n+i] = 0.12
94 # scale spectral norm, maintaining mechanism orientation
95 return M * (scale / np.linalg.svd(M, compute_uv=False)[0])
96
97
98def reservoir_features(u, y, kind, width=32, leak=0.75):
99 # x update is a contractive graph/mechanistic map or vanilla random ESN.
100 T, inp = u.shape
101 outdim = y.shape[1]
102 if kind == 'mechanistic':
103 base = graph_matrix(outdim//2, 0.88)
104 # Lift physical state to width with fixed nonlinear features; graph state is 2n.
105 P = rng.normal(0, 0.35, (width, outdim))
106 A = rng.normal(0, 0.08, (width, width))
107 # Make A contractive and retain graph signal through P/base/P pseudo-lift.
108 A = A * (0.22 / np.linalg.svd(A, compute_uv=False)[0])
109 C = rng.normal(0, 0.25, (width, inp))
110 mech = P @ base
111 else:
112 A = rng.normal(0, 1/np.sqrt(width), (width, width))
113 A *= 0.92 / np.linalg.svd(A, compute_uv=False)[0]
114 C = rng.normal(0, 0.35, (width, inp))
115 mech = rng.normal(0, 0.18, (width, outdim))
116 x = np.zeros(width)
117 Z = []
118 for t in range(T):
119 # undersensed input plus mechanistic relative-state forcing (available as latent model prior)
120 # Both models receive only undersensed channels; hidden channels are unavailable.
121 y_mask = np.zeros_like(y[t])
122 y_mask[::2] = y[t, ::2]
123 forcing = mech @ y_mask
124 x = (1-leak)*x + leak*np.tanh(A @ x + C @ u[t] + forcing)
125 Z.append(np.r_[x, u[t], 1.0])
126 return np.asarray(Z)
127
128
129def ridge_fit_predict(Ztr, Ytr, Zte, reg=1e-4):
130 G = Ztr.T @ Ztr + reg*np.eye(Ztr.shape[1])
131 W = np.linalg.solve(G, Ztr.T @ Ytr)
132 return Zte @ W
133
134
135def prediction_experiment():
136 global rng
137 rng = np.random.default_rng(SEED)
138 u, y = make_data()
139 split = 450
140 # fit one-step state prediction from current reservoir features to next full state
141 results = {}
142 for kind in ['random', 'mechanistic']:
143 Z = reservoir_features(u[:-1], y[:-1], kind)
144 pred = ridge_fit_predict(Z[:split-1], y[1:split], Z[split-1:])
145 mse = np.mean((pred - y[split:])**2)
146 results[kind] = float(mse)
147 return results
148
149
150def main():
151 boundary = boundary_sweep()
152 decay = decay_sweep()
153 pred_ok_boundary = all((r['observed_ratio'] < 1) == (r['predicted_rho'] < 1) for r in boundary)
154 pred_ok_decay = max(r['abs_error'] for r in decay) < 1e-10
155 comparison = prediction_experiment()
156 report = {
157 'seed': SEED,
158 'math_predictions': {
159 'boundary': 'decay iff gamma*Lambda < 1; observed ratio should equal gamma*Lambda',
160 'geometric_decay': 'log norm slope should equal log(rho)',
161 'boundary_sweep': boundary,
162 'decay_sweep': decay,
163 'boundary_confirmed': pred_ok_boundary,
164 'decay_confirmed': pred_ok_decay
165 },
166 'masked_traffic_readout_mse': comparison,
167 'mechanistic_relative_improvement': float((comparison['random']-comparison['mechanistic'])/comparison['random'])
168 }
169 Path('results.json').write_text(json.dumps(report, indent=2))
170 print(json.dumps(report, indent=2))
171
172if __name__ == '__main__':
173 main()