1import json
  2from pathlib import Path
  3import numpy as np
  4
  5SEED = 2046
  6rng = np.random.default_rng(SEED)
  7
  8
  9def fit_koopman(X, Y):
 10    return Y @ np.linalg.pinv(X)
 11
 12
 13def rollout_risk(K, z, H, Q):
 14    zz = np.asarray(z, dtype=float).copy()
 15    vals = []
 16    for _ in range(H):
 17        zz = K @ zz
 18        vals.append(float(zz @ Q @ zz))
 19    return max(vals), np.asarray(vals)
 20
 21
 22def spectral_decay_check():
 23    # Prediction 1: for normal K=rI, perturbations decay exactly as r^j.
 24    rows = []
 25    v = np.array([1.0, -0.7])
 26    for r in (.50, .80, .95):
 27        K = r * np.eye(2)
 28        j = np.arange(1, 11)
 29        observed = np.array([np.linalg.norm(np.linalg.matrix_power(K, k) @ v) /
 30                             np.linalg.norm(v) for k in j])
 31        predicted = r ** j
 32        rows.append({'spectral_radius': r,
 33                     'observed_log_slope': float(np.polyfit(j, np.log(observed), 1)[0]),
 34                     'predicted_log_slope': float(np.log(r)),
 35                     'max_relative_error': float(np.max(np.abs(observed-predicted)/predicted))})
 36    return rows
 37
 38
 39def make_system(n=6000):
 40    # Bounded, burst-driven latent trajectory. Pulses create rising-risk episodes;
 41    # clipping prevents irrelevant numerical divergence.
 42    A = np.array([[0.96, 0.05], [0.0, 0.88]])
 43    x = np.zeros((n + 1, 2))
 44    x[0] = [.02, .0]
 45    forcing = np.zeros((n, 2))
 46    for start in range(180, n, 530):
 47        length = 30
 48        forcing[start:start+length, 0] += np.linspace(.16, .0, length)
 49        forcing[start:start+length, 1] += .035
 50    forcing += rng.normal(0, .0025, forcing.shape)
 51    for t in range(n):
 52        x[t+1] = np.clip(A @ x[t] + forcing[t], -2.0, 2.0)
 53    return A, x
 54
 55
 56def preview_and_gate(states, K, H, tau, Q):
 57    N = len(states) - H - 1
 58    gate = np.zeros(N, dtype=bool)
 59    preview = np.zeros(N)
 60    # Actual risk is a property of the current state. Preview forecasts risk H steps ahead.
 61    actual = np.einsum('ij,jk,ik->i', states[:N], Q, states[:N])
 62    for t in range(N):
 63        preview[t], _ = rollout_risk(K, states[t], H, Q)
 64        gate[t] = preview[t] > tau
 65    return gate, preview, actual
 66
 67
 68def analytic_preview_lead_check():
 69    # Prediction 2: for a scalar growing mode x(t)=r^t x0, preview
 70    # activation precedes current-risk threshold crossing by approximately H.
 71    rows = []
 72    for r in (1.02, 1.05, 1.10):
 73        x = np.array([0.01 * r**t for t in range(500)])
 74        tau = 0.25
 75        crossing = int(np.where(x*x > tau)[0][0])
 76        for H in (2, 5, 8):
 77            activation = int(np.where((x * r**H)**2 > tau)[0][0])
 78            rows.append({'growth_factor': r, 'horizon': H,
 79                         'predicted_lead': H, 'observed_lead': crossing-activation,
 80                         'absolute_error': abs((crossing-activation)-H)})
 81    return rows
 82
 83
 84def horizon_lead_check(states, K, Q, tau):
 85    rows = []
 86    for H in (2, 5, 8, 12):
 87        g, _, actual = preview_and_gate(states, K, H, tau, Q)
 88        high = actual > tau
 89        crossings = np.where(high & ~np.r_[False, high[:-1]])[0]
 90        activations = np.where(g & ~np.r_[False, g[:-1]])[0]
 91        leads = []
 92        for c in crossings:
 93            prior = activations[activations <= c]
 94            if len(prior):
 95                leads.append(c - prior[-1])
 96        rows.append({'horizon': H, 'predicted_max_lead': H,
 97                     'observed_median_lead': float(np.median(leads)) if leads else None,
 98                     'observed_mean_lead': float(np.mean(leads)) if leads else None,
 99                     'events': len(leads)})
100    return rows
101
102
103def threshold_sweep(states, K, Q, H=8):
104    # Prediction 2: increasing tau monotonically decreases refinement calls.
105    _, p, actual0 = preview_and_gate(states, K, H, 0.0, Q)
106    out = []
107    for q in (.50, .70, .82, .90):
108        # Threshold is calibrated on measured current risk, while p is a future preview.
109        tau = float(np.quantile(actual0, q))
110        g, _, actual = preview_and_gate(states, K, H, tau, Q)
111        high = actual > tau
112        out.append({'quantile': q, 'tau': tau, 'gate_rate': float(g.mean()),
113                    'recall': float(np.sum(g & high) / max(np.sum(high), 1)),
114                    'precision': float(np.sum(g & high) / max(np.sum(g), 1))})
115    return out
116
117
118def adaptive_toy(states, K, tau, H=8):
119    Q = np.eye(2)
120    gates, preview, actual = preview_and_gate(states, K, H, tau, Q)
121    true = states[:len(gates)]
122    # Cheap base degrades in high-risk regions; correction is expensive but accurate.
123    severity = np.maximum(0.0, np.linalg.norm(true, axis=1) - np.sqrt(tau))
124    base = true + rng.normal(0, .005, true.shape) + severity[:, None] * np.array([.15, .05])
125    refine = true + rng.normal(0, .004, true.shape)
126    rho = .2
127    gated = np.empty_like(true)
128    previous = base[0]
129    for i, on in enumerate(gates):
130        if on:
131            gated[i] = refine[i]
132        else:
133            gated[i] = rho * previous + (1-rho) * base[i]
134        previous = gated[i]
135    base_mse = float(np.mean((base-true)**2))
136    always_mse = float(np.mean((refine-true)**2))
137    gated_mse = float(np.mean((gated-true)**2))
138    high = actual > tau
139    tp = np.sum(gates & high)
140    return {'koopman_matrix': K.tolist(), 'spectral_radius_fit': float(max(abs(np.linalg.eigvals(K)))),
141            'horizon': H, 'tau': tau, 'gate_rate': float(gates.mean()),
142            'high_risk_rate': float(high.mean()), 'precision': float(tp/max(np.sum(gates),1)),
143            'recall': float(tp/max(np.sum(high),1)), 'base_mse': base_mse,
144            'always_refine_mse': always_mse, 'gated_mse': gated_mse,
145            'relative_mse_increase_vs_always': float(gated_mse/always_mse-1),
146            'refinement_calls_saved_vs_always': float(1-gates.mean()), 'rho_interpolation': rho,
147            'mean_preview_risk': float(preview.mean())}
148
149
150def main():
151    A, states = make_system()
152    split = 3000
153    K = fit_koopman(states[:split].T, states[1:split+1].T)
154    Q = np.eye(2)
155    _, _, actual_train = preview_and_gate(states[split-1:split+1500], K, 8, 0.0, Q)
156    tau = float(np.quantile(actual_train, .82))
157    result = {'seed': SEED, 'true_system_matrix': A.tolist(),
158              'decay_sweep': spectral_decay_check(),
159              'analytic_preview_lead_sweep': analytic_preview_lead_check(),
160              'horizon_lead_sweep': horizon_lead_check(states[split:], K, Q, tau),
161              'threshold_sweep': threshold_sweep(states[split:], K, Q),
162              'adaptive_experiment': adaptive_toy(states[split:], K, tau)}
163    Path('results.json').write_text(json.dumps(result, indent=2))
164    print(json.dumps(result, indent=2))
165
166
167if __name__ == '__main__':
168    main()