1import json
  2import math
  3import numpy as np
  4
  5
  6def rhs(x, lam=0.8, q=0.7):
  7    return -lam * x + q * x * x
  8
  9
 10def exact_lift(x, d):
 11    return np.array([x ** k for k in range(1, d + 1)], dtype=float)
 12
 13
 14def K_matrix(d, lam=0.8, q=0.7):
 15    # For z_k=x^k, d/dt z_k=-k*lam*z_k + k*q*z_{k+1}.
 16    K = np.zeros((d, d))
 17    for k in range(1, d + 1):
 18        K[k - 1, k - 1] = -k * lam
 19        if k < d:
 20            K[k - 1, k] = k * q
 21    return K
 22
 23
 24def residual(x, d, lam=0.8, q=0.7):
 25    z = exact_lift(x, d)
 26    dz = np.array([k * x ** (k - 1) * rhs(x, lam, q) for k in range(1, d + 1)])
 27    return dz - K_matrix(d, lam, q) @ z
 28
 29
 30def scaling_sweep():
 31    # At fixed shape x=a*x0, ||r_d|| must scale as a^(d+1).
 32    amps = np.geomspace(0.15, 1.2, 9)
 33    rows = []
 34    for d in (2, 3, 4):
 35        vals = np.array([np.linalg.norm(residual(a, d)) for a in amps])
 36        slope = np.polyfit(np.log(amps), np.log(vals), 1)[0]
 37        rows.append({'d': d, 'predicted_exponent': d + 1, 'observed_exponent': float(slope),
 38                     'max_relative_fit_error': float(np.max(np.abs(vals / (vals[0] * (amps / amps[0]) ** (d + 1)) - 1)))})
 39    # q enters linearly, including the q=0 vanishing prediction.
 40    qs = np.geomspace(0.1, 1.6, 8)
 41    vals = np.array([np.linalg.norm(residual(0.55, 3, q=q)) for q in qs])
 42    q_slope = np.polyfit(np.log(qs), np.log(vals), 1)[0]
 43    return {'amplitude': rows, 'q_scaling': {'predicted_exponent': 1, 'observed_exponent': float(q_slope)},
 44            'zero_q_residual_norm': float(np.linalg.norm(residual(0.55, 3, q=0.0)))}
 45
 46
 47def rho(x, d, lam=0.8, q=0.7, eps=1e-10):
 48    z = exact_lift(x, d)
 49    return np.linalg.norm(residual(x, d, lam, q)) / (np.linalg.norm(K_matrix(d, lam, q) @ z) + eps)
 50
 51
 52def euler_truth(x, dt, lam, q):
 53    return x + dt * rhs(x, lam, q)
 54
 55
 56def lifted_step(z, d, dt, lam, q):
 57    return z + dt * (K_matrix(d, lam, q) @ z)
 58
 59
 60def rollout(x0, d, steps, dt, lam, q, gate=False, threshold=0.0, dwell=2):
 61    x = float(x0)
 62    active = 2 if gate else d
 63    z = exact_lift(x, active)
 64    over = 0
 65    pred = []
 66    truth = []
 67    active_counts = []
 68    for t in range(steps):
 69        truth.append(x)
 70        # Observation-based certificate; in a real deployment x_next is measured.
 71        x_next = euler_truth(x, dt, lam, q)
 72        if gate and active == 2:
 73            r = rho(x, 2, lam, q)
 74            over = over + 1 if r > threshold else 0
 75            if over >= dwell:
 76                active = 3
 77                z = exact_lift(x, 3)
 78        if active == 2:
 79            z = lifted_step(z, 2, dt, lam, q)
 80        else:
 81            z = lifted_step(z, 3, dt, lam, q)
 82        pred.append(float(z[0]))
 83        active_counts.append(active)
 84        x = x_next
 85    return np.array(pred), np.array(truth), np.array(active_counts)
 86
 87
 88def experiment():
 89    lam, q, dt, steps, x0 = 0.8, 0.7, 0.04, 50, 0.35
 90    # Calibrate up threshold as requested from a low-amplitude training window.
 91    train_x = np.linspace(0.05, 0.35, 100)
 92    train_rhos = np.array([rho(x, 2, lam, q) for x in train_x])
 93    tau = float(np.quantile(train_rhos, 0.95))
 94    results = {}
 95    for name, d, gate in [('fixed_d2', 2, False), ('fixed_d3', 3, False), ('gated', 2, True)]:
 96        pred, truth, active = rollout(x0, d, steps, dt, lam, q, gate, tau, 2)
 97        results[name] = {'one_step_like_final_error': float(abs(pred[0] - truth[1])),
 98                         'rollout_rmse': float(np.sqrt(np.mean((pred - truth) ** 2))),
 99                         'final_abs_error': float(abs(pred[-1] - truth[-1])),
100                         'mean_active_degree': float(np.mean(active)),
101                         'active_feature_count_mean': float(np.mean(active))}
102    results['settings'] = {'tau_up_95pct': tau, 'initial_rho': rho(x0, 2, lam, q),
103                           'final_rho_d2': rho(0.12, 2, lam, q), 'steps': steps, 'dt': dt}
104    # Falsifiable gate transition: predict activation from the certificate sequence,
105    # then compare with the controller's actual first activation time.
106    transition = []
107    for x_init in np.linspace(0.08, 0.55, 10):
108        xs = []
109        x_tmp = float(x_init)
110        for _ in range(50):
111            xs.append(x_tmp)
112            x_tmp = euler_truth(x_tmp, dt, lam, q)
113        cert = np.array([rho(xx, 2, lam, q) for xx in xs])
114        pred_t = None
115        for t in range(len(cert) - 1):
116            if cert[t] > tau and cert[t + 1] > tau:
117                pred_t = t
118                break
119        _, _, act = rollout(x_init, 2, 50, dt, lam, q, True, tau, 2)
120        obs = np.where(act >= 3)[0]
121        obs_t = int(obs[0]) if len(obs) else None
122        transition.append({'x0': float(x_init), 'predicted_activation_step': pred_t,
123                           'observed_activation_step': obs_t,
124                           'rho_x0': float(cert[0])})
125    results['gate_transition_sweep'] = transition
126    return results
127
128
129def finite_difference_check():
130    lam, q, dt = 0.8, 0.7, 1e-6
131    rows = []
132    for d in (2, 3, 4):
133        x = 0.41
134        fd = (exact_lift(euler_truth(x, dt, lam, q), d) - exact_lift(x, d)) / dt
135        analytic = K_matrix(d, lam, q) @ exact_lift(x, d) + residual(x, d, lam, q)
136        rows.append({'d': d, 'relative_error': float(np.linalg.norm(fd - analytic) / np.linalg.norm(analytic))})
137    return rows
138
139
140if __name__ == '__main__':
141    out = {'math_verification': scaling_sweep(), 'finite_difference_check': finite_difference_check(), 'mini_experiment': experiment()}
142    print(json.dumps(out, indent=2))
143
144    out = {'math_verification': scaling_sweep(), 'mini_experiment': experiment()}
145    print(json.dumps(out, indent=2))