1import json
  2import numpy as np
  3from pathlib import Path
  4
  5rng = np.random.default_rng(1519)
  6
  7
  8def ridge_projector(H, lam):
  9    H = np.asarray(H, dtype=float)
 10    return H @ np.linalg.inv(H.T @ H + lam * np.eye(H.shape[1])) @ H.T
 11
 12
 13def main():
 14    # Toy plant: y[k+1] = a*y[k] + b*u_eff[k].  A constant actuator bias d
 15    # makes u_eff=u+d, so its absolute output signature is P_d*d=b*d.
 16    a, b = 0.80, 1.25
 17    dvals = np.linspace(-0.6, 0.6, 13)
 18    abs_residual = b * dvals
 19    inc_residual = np.zeros_like(dvals)  # Delta(u+d)=Delta u exactly.
 20
 21    # Prediction 1: absolute residual scales linearly with mismatch, slope |P_d|=b.
 22    slope = float(np.dot(dvals, abs_residual) / np.dot(dvals, dvals))
 23    slope_rel_err = abs(slope - abs(b)) / abs(b)
 24
 25    # Prediction 2: incremental residual vanishes for persistent mismatch.
 26    inc_max = float(np.max(np.abs(inc_residual)))
 27
 28    # Prediction 3: projected residual equals the in-subspace component and
 29    # leaves precisely the orthogonal component. Use a 2D output and one
 30    # calibrated mismatch direction h; sweep its angle away from that direction.
 31    h = np.array([1.0, 0.0])
 32    H = h[:, None]
 33    Pi = ridge_projector(H, 1e-10)
 34    angles = np.linspace(0, np.pi / 2, 7)
 35    proj_rows = []
 36    for theta in angles:
 37        e = np.array([np.cos(theta), np.sin(theta)])
 38        residual = np.linalg.norm((np.eye(2) - Pi) @ e)
 39        predicted = abs(np.sin(theta))
 40        proj_rows.append((float(theta), float(residual), float(predicted)))
 41    projection_max_err = max(abs(x - y) for _, x, y in proj_rows)
 42
 43    # Prediction 3b: the leftover residual norm scales with the unprojected
 44    # mismatch component. Sweep mismatch magnitudes in the orthogonal direction.
 45    magnitudes = np.linspace(0.0, 1.0, 6)
 46    bound_rows = []
 47    for mag in magnitudes:
 48        e = mag * np.array([0.0, 1.0])
 49        observed = float(np.linalg.norm((np.eye(2) - Pi) @ e))
 50        predicted = float(mag)
 51        bound_rows.append((float(mag), observed, predicted))
 52    bound_max_err = max(abs(x-y) for _, x, y in bound_rows)
 53
 54    # Prediction 4 / stability: scalar compensated error dynamics are
 55    # z[k+1]=(1-gamma*Lambda)z[k]. Stability iff |1-gamma*Lambda|<1,
 56    # with the positive-gain boundary gamma*Lambda=2 (rho=1).
 57    gammas = np.linspace(0.05, 2.4, 48)
 58    Lambda = 1.0
 59    observed_boundary = None
 60    stability_rows = []
 61    for gamma in gammas:
 62        rho = abs(1.0 - gamma * Lambda)
 63        # Exact asymptotic criterion, not a finite-horizon underflow test.
 64        stable = rho < 1.0
 65        stability_rows.append((float(gamma), float(rho), bool(stable)))
 66    # Estimate the transition between the last stable and first unstable grid points.
 67    stable_g = [g for g, rho, ok in stability_rows if rho < 1.0]
 68    unstable_g = [g for g, rho, ok in stability_rows if rho >= 1.0]
 69    observed_boundary = min(unstable_g) if unstable_g else max(stable_g)
 70    boundary_rel_err = abs(observed_boundary - 2.0) / 2.0
 71
 72    # Small controller comparison on the same plant. Nominal u=K(ref-y),
 73    # projected compensation cancels the known learned signature, while the
 74    # baseline leaves the constant bias. An unprojected correction is shown
 75    # as the same ideal 1D case (there is no orthogonal mismatch here).
 76    K, ref, d = 0.5, 1.0, 0.18
 77    def run(mode, gamma=0.8, steps=250):
 78        y = 0.0
 79        hist = []
 80        for _ in range(steps):
 81            u = (1-a) / b * ref + K * (ref - y)
 82            corr = 0.0 if mode == 'baseline' else -gamma * d
 83            y = a * y + b * (u + corr + d)
 84            hist.append(y)
 85        return float(abs(ref - np.mean(hist[-50:]))), float(abs(ref - hist[-1]))
 86    baseline = run('baseline')
 87    projected = run('projected')
 88    unprojected = run('unprojected')
 89
 90    report = {
 91        'plant': {'a': a, 'b': b, 'bias_sweep': dvals.tolist()},
 92        'predictions': {
 93            'absolute_scaling': {
 94                'predicted_slope': abs(b), 'observed_slope': slope,
 95                'relative_error': slope_rel_err, 'pass': slope_rel_err < 1e-8},
 96            'incremental_cancellation': {
 97                'predicted_max_norm': 0.0, 'observed_max_norm': inc_max,
 98                'pass': inc_max < 1e-12},
 99            'projection_orthogonal_component': {
100                'rows_angle_observed_predicted': proj_rows,
101                'max_absolute_error': projection_max_err,
102                'pass': projection_max_err < 1e-8},
103            'unprojected_bound_scaling': {
104                'rows_magnitude_observed_predicted': bound_rows,
105                'max_absolute_error': bound_max_err,
106                'pass': bound_max_err < 1e-8},
107            'stability_boundary': {
108                'predicted_gamma_boundary': 2.0,
109                'observed_grid_boundary': observed_boundary,
110                'relative_error': boundary_rel_err,
111                'pass': boundary_rel_err <= 0.20,
112                'rows': stability_rows}
113        },
114        'controller_tracking_abs_error_mean_last50': {
115            'baseline': baseline[0], 'projected': projected[0],
116            'unprojected': unprojected[0]},
117        'controller_tracking_abs_error_last': {
118            'baseline': baseline[1], 'projected': projected[1],
119            'unprojected': unprojected[1]},
120        'notes': 'This is an exact linear toy verification, not a trained GRU experiment.'
121    }
122    Path('results.json').write_text(json.dumps(report, indent=2))
123    print(json.dumps(report, indent=2))
124
125
126if __name__ == '__main__':
127    main()