import json, math from pathlib import Path import numpy as np SEED = 1487 rng = np.random.default_rng(SEED) def orthogonal(n, seed): q, _ = np.linalg.qr(np.random.default_rng(seed).normal(size=(n, n))) return q def boundary_sweep(): # F(x)=gamma*A*x, ||A||_2=Lambda. Prediction: decay iff gamma*Lambda<1. n, steps, lam = 8, 40, 1.2 Q = orthogonal(n, 1) x0 = np.ones(n) / np.sqrt(n) rows = [] for gamma in [0.60, 0.75, 0.80, 0.833333, 0.90, 1.00]: F = gamma * lam * Q x = x0.copy() norms = [] for _ in range(steps): x = F @ x norms.append(np.linalg.norm(x)) observed = np.mean(np.array(norms[-5:]) / np.array(norms[-6:-1])) predicted = gamma * lam rows.append({"gamma": gamma, "predicted_rho": predicted, "observed_ratio": float(observed), "classification_pred": "decay" if predicted < 1 else "grow_or_neutral", "classification_obs": "decay" if observed < 0.999 else "grow_or_neutral"}) return rows def decay_sweep(): # F=rho*Q is an isometry scaled by rho. Prediction log(norm_k/norm_0)/k=log(rho). n, steps = 10, 30 Q = orthogonal(n, 2) x0 = np.ones(n) / np.sqrt(n) rows = [] for rho in [0.50, 0.70, 0.90, 0.97]: x = x0.copy() for _ in range(steps): x = rho * Q @ x measured = math.log(np.linalg.norm(x) / np.linalg.norm(x0)) / steps rows.append({"rho": rho, "predicted_log_slope": math.log(rho), "measured_log_slope": float(measured), "abs_error": float(abs(measured - math.log(rho)))}) return rows def idm_acc(v, gap, dv): # Safe numerical IDM implementation; dv=v-v_leader. amax, v0, delta, s0, T, b = 1.2, 1.4, 4.0, 0.5, 1.1, 1.5 gap = max(float(gap), 0.15) sstar = s0 + v*T + v*dv/(2*np.sqrt(amax*b)) return amax * (1 - (max(v, 0)/v0)**delta - (sstar/gap)**2) def make_data(T=700, n=8): # Ring traffic-like system: position and speed, only every other position observed. dt = 0.08 pos = np.zeros((T, n)); vel = np.zeros((T, n)) pos[0] = np.linspace(0, 8, n, endpoint=False) vel[0] = 0.8 + 0.05*rng.normal(size=n) for t in range(T-1): for i in range(n): lead = (i + 1) % n gap = (pos[t, lead] - pos[t, i]) % 8.0 dv = vel[t, i] - vel[t, lead] acc = idm_acc(vel[t,i], gap, dv) vel[t+1,i] = np.clip(vel[t,i] + dt*acc + 0.008*rng.normal(), -0.2, 2.0) pos[t+1,i] = (pos[t,i] + dt*vel[t+1,i]) % 8.0 # normalize physical channels for stable reservoir driving y = np.concatenate([pos/8.0, vel/2.0], axis=1) obs_idx = np.arange(0, 2*n, 2) # positions of half the cars only u = y[:, obs_idx] return u, y def graph_matrix(n, scale=0.88): # Mechanistic directed predecessor coupling, normalized to prescribed norm. M = np.zeros((2*n, 2*n)) for i in range(n): lead = (i+1) % n # position error and relative velocity drive each vehicle's two channels M[n+i, i] = 0.20 M[n+i, lead] = -0.20 M[n+i, n+i] = 0.72 M[n+i, n+lead] = 0.12 M[i, i] = 0.95 M[i, n+i] = 0.12 # scale spectral norm, maintaining mechanism orientation return M * (scale / np.linalg.svd(M, compute_uv=False)[0]) def reservoir_features(u, y, kind, width=32, leak=0.75): # x update is a contractive graph/mechanistic map or vanilla random ESN. T, inp = u.shape outdim = y.shape[1] if kind == 'mechanistic': base = graph_matrix(outdim//2, 0.88) # Lift physical state to width with fixed nonlinear features; graph state is 2n. P = rng.normal(0, 0.35, (width, outdim)) A = rng.normal(0, 0.08, (width, width)) # Make A contractive and retain graph signal through P/base/P pseudo-lift. A = A * (0.22 / np.linalg.svd(A, compute_uv=False)[0]) C = rng.normal(0, 0.25, (width, inp)) mech = P @ base else: A = rng.normal(0, 1/np.sqrt(width), (width, width)) A *= 0.92 / np.linalg.svd(A, compute_uv=False)[0] C = rng.normal(0, 0.35, (width, inp)) mech = rng.normal(0, 0.18, (width, outdim)) x = np.zeros(width) Z = [] for t in range(T): # undersensed input plus mechanistic relative-state forcing (available as latent model prior) # Both models receive only undersensed channels; hidden channels are unavailable. y_mask = np.zeros_like(y[t]) y_mask[::2] = y[t, ::2] forcing = mech @ y_mask x = (1-leak)*x + leak*np.tanh(A @ x + C @ u[t] + forcing) Z.append(np.r_[x, u[t], 1.0]) return np.asarray(Z) def ridge_fit_predict(Ztr, Ytr, Zte, reg=1e-4): G = Ztr.T @ Ztr + reg*np.eye(Ztr.shape[1]) W = np.linalg.solve(G, Ztr.T @ Ytr) return Zte @ W def prediction_experiment(): global rng rng = np.random.default_rng(SEED) u, y = make_data() split = 450 # fit one-step state prediction from current reservoir features to next full state results = {} for kind in ['random', 'mechanistic']: Z = reservoir_features(u[:-1], y[:-1], kind) pred = ridge_fit_predict(Z[:split-1], y[1:split], Z[split-1:]) mse = np.mean((pred - y[split:])**2) results[kind] = float(mse) return results def main(): boundary = boundary_sweep() decay = decay_sweep() pred_ok_boundary = all((r['observed_ratio'] < 1) == (r['predicted_rho'] < 1) for r in boundary) pred_ok_decay = max(r['abs_error'] for r in decay) < 1e-10 comparison = prediction_experiment() report = { 'seed': SEED, 'math_predictions': { 'boundary': 'decay iff gamma*Lambda < 1; observed ratio should equal gamma*Lambda', 'geometric_decay': 'log norm slope should equal log(rho)', 'boundary_sweep': boundary, 'decay_sweep': decay, 'boundary_confirmed': pred_ok_boundary, 'decay_confirmed': pred_ok_decay }, 'masked_traffic_readout_mse': comparison, 'mechanistic_relative_improvement': float((comparison['random']-comparison['mechanistic'])/comparison['random']) } Path('results.json').write_text(json.dumps(report, indent=2)) print(json.dumps(report, indent=2)) if __name__ == '__main__': main()