import json from pathlib import Path import numpy as np SEED = 2046 rng = np.random.default_rng(SEED) def fit_koopman(X, Y): return Y @ np.linalg.pinv(X) def rollout_risk(K, z, H, Q): zz = np.asarray(z, dtype=float).copy() vals = [] for _ in range(H): zz = K @ zz vals.append(float(zz @ Q @ zz)) return max(vals), np.asarray(vals) def spectral_decay_check(): # Prediction 1: for normal K=rI, perturbations decay exactly as r^j. rows = [] v = np.array([1.0, -0.7]) for r in (.50, .80, .95): K = r * np.eye(2) j = np.arange(1, 11) observed = np.array([np.linalg.norm(np.linalg.matrix_power(K, k) @ v) / np.linalg.norm(v) for k in j]) predicted = r ** j rows.append({'spectral_radius': r, 'observed_log_slope': float(np.polyfit(j, np.log(observed), 1)[0]), 'predicted_log_slope': float(np.log(r)), 'max_relative_error': float(np.max(np.abs(observed-predicted)/predicted))}) return rows def make_system(n=6000): # Bounded, burst-driven latent trajectory. Pulses create rising-risk episodes; # clipping prevents irrelevant numerical divergence. A = np.array([[0.96, 0.05], [0.0, 0.88]]) x = np.zeros((n + 1, 2)) x[0] = [.02, .0] forcing = np.zeros((n, 2)) for start in range(180, n, 530): length = 30 forcing[start:start+length, 0] += np.linspace(.16, .0, length) forcing[start:start+length, 1] += .035 forcing += rng.normal(0, .0025, forcing.shape) for t in range(n): x[t+1] = np.clip(A @ x[t] + forcing[t], -2.0, 2.0) return A, x def preview_and_gate(states, K, H, tau, Q): N = len(states) - H - 1 gate = np.zeros(N, dtype=bool) preview = np.zeros(N) # Actual risk is a property of the current state. Preview forecasts risk H steps ahead. actual = np.einsum('ij,jk,ik->i', states[:N], Q, states[:N]) for t in range(N): preview[t], _ = rollout_risk(K, states[t], H, Q) gate[t] = preview[t] > tau return gate, preview, actual def analytic_preview_lead_check(): # Prediction 2: for a scalar growing mode x(t)=r^t x0, preview # activation precedes current-risk threshold crossing by approximately H. rows = [] for r in (1.02, 1.05, 1.10): x = np.array([0.01 * r**t for t in range(500)]) tau = 0.25 crossing = int(np.where(x*x > tau)[0][0]) for H in (2, 5, 8): activation = int(np.where((x * r**H)**2 > tau)[0][0]) rows.append({'growth_factor': r, 'horizon': H, 'predicted_lead': H, 'observed_lead': crossing-activation, 'absolute_error': abs((crossing-activation)-H)}) return rows def horizon_lead_check(states, K, Q, tau): rows = [] for H in (2, 5, 8, 12): g, _, actual = preview_and_gate(states, K, H, tau, Q) high = actual > tau crossings = np.where(high & ~np.r_[False, high[:-1]])[0] activations = np.where(g & ~np.r_[False, g[:-1]])[0] leads = [] for c in crossings: prior = activations[activations <= c] if len(prior): leads.append(c - prior[-1]) rows.append({'horizon': H, 'predicted_max_lead': H, 'observed_median_lead': float(np.median(leads)) if leads else None, 'observed_mean_lead': float(np.mean(leads)) if leads else None, 'events': len(leads)}) return rows def threshold_sweep(states, K, Q, H=8): # Prediction 2: increasing tau monotonically decreases refinement calls. _, p, actual0 = preview_and_gate(states, K, H, 0.0, Q) out = [] for q in (.50, .70, .82, .90): # Threshold is calibrated on measured current risk, while p is a future preview. tau = float(np.quantile(actual0, q)) g, _, actual = preview_and_gate(states, K, H, tau, Q) high = actual > tau out.append({'quantile': q, 'tau': tau, 'gate_rate': float(g.mean()), 'recall': float(np.sum(g & high) / max(np.sum(high), 1)), 'precision': float(np.sum(g & high) / max(np.sum(g), 1))}) return out def adaptive_toy(states, K, tau, H=8): Q = np.eye(2) gates, preview, actual = preview_and_gate(states, K, H, tau, Q) true = states[:len(gates)] # Cheap base degrades in high-risk regions; correction is expensive but accurate. severity = np.maximum(0.0, np.linalg.norm(true, axis=1) - np.sqrt(tau)) base = true + rng.normal(0, .005, true.shape) + severity[:, None] * np.array([.15, .05]) refine = true + rng.normal(0, .004, true.shape) rho = .2 gated = np.empty_like(true) previous = base[0] for i, on in enumerate(gates): if on: gated[i] = refine[i] else: gated[i] = rho * previous + (1-rho) * base[i] previous = gated[i] base_mse = float(np.mean((base-true)**2)) always_mse = float(np.mean((refine-true)**2)) gated_mse = float(np.mean((gated-true)**2)) high = actual > tau tp = np.sum(gates & high) return {'koopman_matrix': K.tolist(), 'spectral_radius_fit': float(max(abs(np.linalg.eigvals(K)))), 'horizon': H, 'tau': tau, 'gate_rate': float(gates.mean()), 'high_risk_rate': float(high.mean()), 'precision': float(tp/max(np.sum(gates),1)), 'recall': float(tp/max(np.sum(high),1)), 'base_mse': base_mse, 'always_refine_mse': always_mse, 'gated_mse': gated_mse, 'relative_mse_increase_vs_always': float(gated_mse/always_mse-1), 'refinement_calls_saved_vs_always': float(1-gates.mean()), 'rho_interpolation': rho, 'mean_preview_risk': float(preview.mean())} def main(): A, states = make_system() split = 3000 K = fit_koopman(states[:split].T, states[1:split+1].T) Q = np.eye(2) _, _, actual_train = preview_and_gate(states[split-1:split+1500], K, 8, 0.0, Q) tau = float(np.quantile(actual_train, .82)) result = {'seed': SEED, 'true_system_matrix': A.tolist(), 'decay_sweep': spectral_decay_check(), 'analytic_preview_lead_sweep': analytic_preview_lead_check(), 'horizon_lead_sweep': horizon_lead_check(states[split:], K, Q, tau), 'threshold_sweep': threshold_sweep(states[split:], K, Q), 'adaptive_experiment': adaptive_toy(states[split:], K, tau)} Path('results.json').write_text(json.dumps(result, indent=2)) print(json.dumps(result, indent=2)) if __name__ == '__main__': main()