1import json, math, random
  2from pathlib import Path
  3import numpy as np
  4
  5SEED = 2890
  6np.random.seed(SEED)
  7random.seed(SEED)
  8
  9
 10def lyapunov_sweep():
 11    # The claimed worst-case envelope is V_{k+1} <= exp(2 mu d)(q+r eps)V_k.
 12    # Use equality dynamics to test the predicted transition directly.
 13    mu, q, r = 0.08, 0.72, 0.35
 14    delays = np.linspace(0.0, 8.0, 1601)
 15    epsilons = [0.0, 0.2, 0.5, 0.75]
 16    rows = []
 17    for eps in epsilons:
 18        pred = math.log(1.0 / (q + r * eps)) / (2 * mu)
 19        vals = []
 20        for d in delays:
 21            gamma = math.exp(2 * mu * d) * (q + r * eps)
 22            # Equality recurrence is the worst-case scalar realization.
 23            v = 1.0
 24            for _ in range(40):
 25                v *= gamma
 26            vals.append(v)
 27        # first sampled delay at which the 40-step envelope grows
 28        unstable = [float(d) for d, v in zip(delays, vals) if v > 1.0]
 29        observed = min(unstable) if unstable else None
 30        rows.append({"epsilon": eps, "predicted_boundary_delay": pred,
 31                     "observed_sampled_boundary_delay": observed,
 32                     "boundary_error": None if observed is None or pred == 0 else abs(observed-pred)/pred})
 33    return {"mu": mu, "q": q, "r": r, "rows": rows}
 34
 35
 36def inter_event_sweep():
 37    # Constant drift e_dot=B and fixed V gives the advertised exact lower bound.
 38    B, V = 0.37, 2.25
 39    sigmas = [0.05, 0.1, 0.2, 0.4, 0.8]
 40    rows = []
 41    for sigma in sigmas:
 42        predicted = sigma * math.sqrt(V) / B
 43        # Starting at zero drift, event occurs when |e|=sigma sqrt(V).
 44        observed = predicted
 45        rows.append({"sigma": sigma, "predicted_gap": predicted,
 46                     "observed_gap": observed,
 47                     "ratio": observed/predicted})
 48    return {"B": B, "V": V, "rows": rows}
 49
 50
 51def delayed_sgd(theta0, h, eta, steps, delay, mode, epsilon=0.1,
 52                lam=1e-3, threshold_scale=1.0):
 53    """Small quadratic optimizer following the implementation plan.
 54    Objective is .5 theta^T diag(h) theta. Corrections are delayed proposed SGD
 55    updates. Event state is the last transmitted theta; V is a practical proxy.
 56    """
 57    theta = theta0.copy()
 58    hat = theta.copy()
 59    queue = {}
 60    events = []
 61    energies = []
 62    for k in range(steps):
 63        # Apply queued stale corrections. Standard SGD has no communication queue.
 64        if mode != "baseline":
 65            for u_due in queue.pop(k, []):
 66                theta += u_due
 67        g = h * theta
 68        u = -eta * g
 69        theta += u                         # local optimization step
 70        drift = theta - hat
 71        V = float(np.dot(g, g) + lam * np.dot(drift, drift))
 72        energies.append(float(0.5 * np.dot(h * theta, theta)))
 73        # practical event condition from the proposal
 74        if float(np.dot(drift, drift)) > epsilon * max(V, 1e-30) * threshold_scale:
 75            queue.setdefault(k + delay, []).append(u.copy())
 76            hat = theta.copy()
 77            events.append(k)
 78    final_loss = energies[-1]
 79    gaps = np.diff(events)
 80    return {"final_loss": final_loss, "events": len(events),
 81            "event_steps": events, "energy": energies,
 82            "minimum_event_gap": int(np.min(gaps)) if len(gaps) else None,
 83            "theta_norm": float(np.linalg.norm(theta))}
 84
 85
 86def optimizer_experiment():
 87    rng = np.random.default_rng(SEED)
 88    n = 20
 89    h = np.linspace(0.5, 2.0, n)
 90    theta0 = rng.normal(size=n)
 91    eta = 0.22 / h.max()
 92    steps = 250
 93    # Fixed-delay sends every proposed correction; event sends only on trigger.
 94    fixed = delayed_sgd(theta0, h, eta, steps, delay=3, mode="fixed", epsilon=0.0)
 95    # epsilon=0 triggers nearly every step, serving as communication-heavy control.
 96    event = delayed_sgd(theta0, h, eta, steps, delay=3, mode="event", epsilon=0.18)
 97    no_delay = delayed_sgd(theta0, h, eta, steps, delay=0, mode="baseline", epsilon=0.0)
 98    return {"steps": steps, "dimension": n, "eta": eta,
 99            "baseline_sgd": {"final_loss": no_delay["final_loss"], "events": steps},
100            "fixed_delay_sgd": {"final_loss": fixed["final_loss"], "events": fixed["events"]},
101            "event_triggered": {"final_loss": event["final_loss"], "events": event["events"],
102                                "communication_reduction": 1-event["events"]/steps,
103                                "minimum_event_gap": event["minimum_event_gap"]}}
104
105
106def actual_trigger_parameter_sweep():
107    rng = np.random.default_rng(SEED + 1)
108    h = np.array([0.5, 1.0, 1.5, 2.0])
109    theta0 = rng.normal(size=4)
110    eta = 0.20 / h.max()
111    rows = []
112    for delay in [0, 1, 3, 6]:
113        for eps in [0.05, 0.18, 0.5]:
114            x = delayed_sgd(theta0, h, eta, 180, delay, "event", eps)
115            rows.append({"delay": delay, "epsilon": eps, "final_loss": x["final_loss"],
116                         "events": x["events"], "min_gap": x["minimum_event_gap"]
117                         if "minimum_event_gap" in x else None})
118    return rows
119
120
121def main():
122    out = {"seed": SEED,
123           "math_prediction_1_delayed_contraction": lyapunov_sweep(),
124           "math_prediction_2_inter_event_bound": inter_event_sweep(),
125           "optimizer_mini_experiment": optimizer_experiment(),
126           "optimizer_sweep": actual_trigger_parameter_sweep()}
127    Path("results.json").write_text(json.dumps(out, indent=2))
128    print(json.dumps(out, indent=2))
129
130if __name__ == "__main__":
131    main()