1import json
  2import random
  3from pathlib import Path
  4import numpy as np
  5
  6
  7def simulate(theta0, theta1, t_end, dt, kappa, theta_c, slope=1.0, mode="baseline"):
  8    n = int(round(t_end / dt)) + 1
  9    ts = np.arange(n) * dt
 10    desired = theta0 + (theta1 - theta0) * np.minimum(ts / t_end, 1.0)
 11    r = (theta1 - theta0) / t_end
 12    if mode == "baseline":
 13        command = desired.copy()
 14    elif mode == "minus":
 15        command = desired - r / kappa
 16    elif mode == "plus":
 17        command = desired + r / kappa
 18    else:
 19        raise ValueError(mode)
 20    eff = np.empty(n)
 21    eff[0] = theta0
 22    for i in range(1, n):
 23        eff[i] = eff[i - 1] + dt * kappa * (command[i - 1] - eff[i - 1])
 24    lam_eff = slope * (eff - theta_c)
 25    lam_cmd = slope * (command - theta_c)
 26
 27    def first_cross(x):
 28        ix = np.flatnonzero(x >= 0)
 29        if len(ix) == 0:
 30            return float("nan")
 31        i = int(ix[0])
 32        if i == 0:
 33            return float(ts[0])
 34        # Linear interpolation makes timing less dependent on dt.
 35        return float(ts[i - 1] + dt * (-x[i - 1]) / (x[i] - x[i - 1]))
 36
 37    target_cross = (theta_c - theta0) / r
 38    return {
 39        "target_cross": target_cross,
 40        "command_cross": first_cross(lam_cmd),
 41        "effective_cross": first_cross(lam_eff),
 42        "delay": first_cross(lam_eff) - target_cross,
 43        "ts": ts, "desired": desired, "eff": eff,
 44    }
 45
 46
 47def linear_fit(x, y):
 48    x, y = np.asarray(x), np.asarray(y)
 49    coef = np.polyfit(x, y, 1)
 50    pred = np.polyval(coef, x)
 51    ssr = np.sum((y - pred) ** 2)
 52    sst = np.sum((y - y.mean()) ** 2)
 53    return float(coef[0]), float(coef[1]), float(1 - ssr / sst)
 54
 55
 56def math_verification():
 57    # Prediction 1: after ramp transients, theta_eff = theta_desired-r/kappa.
 58    kappa = 4.0
 59    t_end, dt = 20.0, 0.002
 60    out = simulate(0, 2, t_end, dt, kappa, 99, mode="baseline")
 61    r = 2 / t_end
 62    tail_error = np.mean((out["desired"][-1000:] - out["eff"][-1000:]))
 63    predicted_offset = r / kappa
 64
 65    # Prediction 2: threshold delay is (approximately) r/(kappa*slope).
 66    kappas = np.array([1., 2., 4., 8., 16., 32.])
 67    delays = np.array([simulate(0, 2, 20, .002, k, 1, mode="baseline")["delay"] for k in kappas])
 68    slope, intercept, r2 = linear_fit(1 / kappas, delays)
 69    predicted_slope = 1.0  # crossing delay=(r/kappa)/(|d theta_c/dt|)=1/kappa
 70
 71    # Prediction 3: inverse feed-forward reduces delay; printed proposal sign
 72    # is tested separately because for an increasing ramp it has the wrong sign.
 73    base = np.array([simulate(0, 2, 20, .002, k, 1, mode="baseline")["delay"] for k in kappas])
 74    minus = np.array([simulate(0, 2, 20, .002, k, 1, mode="minus")["delay"] for k in kappas])
 75    plus = np.array([simulate(0, 2, 20, .002, k, 1, mode="plus")["delay"] for k in kappas])
 76    return {
 77        "offset_observed": float(tail_error), "offset_predicted": predicted_offset,
 78        "delay_fit_slope": slope, "delay_fit_intercept": intercept,
 79        "delay_predicted_slope": predicted_slope, "delay_r2": r2,
 80        "kappas": kappas.tolist(), "baseline_delays": base.tolist(),
 81        "minus_delays": minus.tolist(), "plus_delays": plus.tolist(),
 82        "minus_reduction_fraction": float(1 - np.nanmean(np.abs(minus)) / np.nanmean(np.abs(base))),
 83        "plus_reduction_fraction": float(1 - np.nanmean(np.abs(plus)) / np.nanmean(np.abs(base))),
 84    }
 85
 86
 87def train_toy(seed=7, kappa=4.0, mode="baseline", steps=500):
 88    # 2-D binary classification; learning rate itself is EMA-filtered.
 89    rng = np.random.default_rng(seed)
 90    x0 = rng.normal([-0.7, -0.7], .8, (96, 2))
 91    x1 = rng.normal([0.7, 0.7], .8, (96, 2))
 92    X = np.vstack([x0, x1]); y = np.r_[np.zeros(96), np.ones(96)]
 93    w = np.zeros(2); b = 0.0; lr_eff = .01; lr0, lr1 = .01, .08
 94    losses = []
 95    dt = 1.0
 96    for step in range(steps):
 97        frac = min(step / 150, 1.0)
 98        desired = lr0 + (lr1 - lr0) * frac
 99        rate = (lr1 - lr0) / 150 if step < 150 else 0.0
100        if mode == "baseline": command = desired
101        elif mode == "minus": command = max(.0001, desired - rate / kappa)
102        elif mode == "plus": command = desired + rate / kappa
103        alpha = min(0.9, 0.1 * kappa)
104        lr_eff += alpha * (command - lr_eff)
105        z = X @ w + b
106        p = 1 / (1 + np.exp(-np.clip(z, -30, 30)))
107        loss = -np.mean(y*np.log(p+1e-8) + (1-y)*np.log(1-p+1e-8))
108        losses.append(float(loss))
109        gw = X.T @ (p-y) / len(y); gb = float(np.mean(p-y))
110        w -= lr_eff * gw; b -= lr_eff * gb
111        if not np.isfinite(loss) or np.linalg.norm(w) > 1e6:
112            return {"final_loss": float("inf"), "best_loss": float(np.min(losses)), "exploded": True}
113    return {"final_loss": losses[-1], "best_loss": float(np.min(losses)), "exploded": False}
114
115
116def main():
117    random.seed(7); np.random.seed(7)
118    math = math_verification()
119    training = {m: train_toy(mode=m, steps=300) for m in ["baseline", "minus", "plus"]}
120    result = {"math_verification": math, "training": training}
121    Path("results.json").write_text(json.dumps(result, indent=2))
122    print(json.dumps(result, indent=2))
123
124
125if __name__ == "__main__":
126    main()