1import json
  2import math
  3import numpy as np
  4from scipy.optimize import brentq
  5
  6SEED = 2029
  7np.random.seed(SEED)
  8
  9
 10def ring_adj(n):
 11    A = np.zeros((n, n))
 12    for i in range(n):
 13        A[i, (i - 1) % n] = 1.0
 14        A[i, (i + 1) % n] = 1.0
 15    return A
 16
 17
 18def complete_adj(n):
 19    A = np.ones((n, n)) - np.eye(n)
 20    return A
 21
 22
 23def delayed_linear(A, eps, k, delay_steps, h, steps, x0=None, targets=None):
 24    n = len(A)
 25    d = A.sum(1)
 26    x = np.ones(n) if x0 is None else np.asarray(x0, float).copy()
 27    targets = np.zeros(n) if targets is None else np.asarray(targets, float)
 28    # Constant prehistory is the standard well-posed history for this test.
 29    hist = [x.copy() for _ in range(delay_steps + 1)]
 30    xs = np.empty((steps + 1, n)); xs[0] = x
 31    for t in range(steps):
 32        stale = hist[0]
 33        v = -eps * (x - targets) + k * (A @ stale - d * x)
 34        x = x + h * v
 35        hist.pop(0); hist.append(x.copy())
 36        xs[t + 1] = x
 37    return xs
 38
 39
 40def rate_from_trace(trace, h, lo, hi):
 41    y = np.abs(trace[lo:hi])
 42    t = np.arange(lo, hi) * h
 43    slope = np.polyfit(t, np.log(np.maximum(y, 1e-14)), 1)[0]
 44    return float(-slope)
 45
 46
 47def exact_slow_rate(eps, kd, tau):
 48    if tau == 0: return eps
 49    # q is the positive decay rate of q=eps+kd*(1-exp(-q*tau)).
 50    # For x(t)=exp(-q t), q=eps+kd*(1-exp(q*tau)); the
 51    # principal real decay root lies between 0 and eps.
 52    f = lambda q: q - eps - kd * (1.0 - math.exp(q * tau))
 53    return brentq(f, 0.0, eps)
 54
 55
 56def main():
 57    h = 0.001
 58    n = 4
 59    A = ring_adj(n); degree = 2.0
 60    eps = 0.2
 61    # Prediction 1: slow collective rate scales as 1/(1+k*d*tau).
 62    speed_rows = []
 63    for k in [0.25, 0.5, 1.0, 2.0]:
 64        for tau in [0.0, 0.25, 0.5, 1.0]:
 65            delay_steps = round(tau / h)
 66            tr = delayed_linear(A, eps, k, delay_steps, h, 12000)
 67            # discard transient, and use mean (consensus) mode
 68            q = rate_from_trace(tr.mean(1), h, 4000, 11000)
 69            approx = eps / (1 + k * degree * tau)
 70            exact = exact_slow_rate(eps, k * degree, tau)
 71            speed_rows.append(dict(k=k, tau=tau, observed=q, approx=approx, exact=exact,
 72                                   approx_rel_error=abs(q-approx)/approx))
 73
 74    # Prediction 2: the relevant dimensionless synchronization ratio is
 75    # rho=eps*||J||/(k*lambda_2); disagreement should rise around rho ~ 1.
 76    # Heterogeneous quadratic minima create a persistent disagreement signal.
 77    targets = np.array([-1., 1., -1., 1.])
 78    lam2 = 2.0  # ring of four
 79    sync_rows = []
 80    for k in [0.05, 0.1, 0.2, 0.4, 0.8]:
 81        tr = delayed_linear(A, eps, k, 0, h, 10000, x0=np.zeros(n), targets=targets)
 82        tail = tr[-1000:]
 83        disagreement = np.mean(np.sum((tail - tail.mean(1, keepdims=True))**2, axis=1))
 84        rho = eps / (k * lam2)
 85        sync_rows.append(dict(k=k, rho=rho, disagreement=disagreement))
 86
 87    # Prediction 3: with no local field, disagreement decay is set by k*lambda_2.
 88    gap_rows = []
 89    for name, B, gap in [('ring', A, 2.0), ('complete', complete_adj(n), 4.0)]:
 90        for k in [0.1, 0.2, 0.4]:
 91            x0 = np.array([1., -1., 0.5, 0.2])
 92            tr = delayed_linear(B, 0.0, k, 0, h, 5000, x0=x0)
 93            dis = np.sqrt(np.mean((tr - tr.mean(1, keepdims=True))**2, axis=1))
 94            q = rate_from_trace(dis, h, 500, 4000)
 95            gap_rows.append(dict(graph=name, k=k, lambda2=gap, observed_rate=q,
 96                                 predicted_rate=k*gap))
 97
 98    # Small optimizer-facing comparison: four workers minimize local quadratics.
 99    # Independent SGD is the standard local baseline; delayed diffusion is the idea.
100    opt = {}
101    for name, delay in [('local', None), ('delayed_consensus', 5)]:
102        x = np.zeros(n); hist = [x.copy()] * (delay + 1 if delay is not None else 1)
103        for _ in range(3000):
104            if delay is None:
105                x += 0.01 * (-eps * (x - targets))  # gradient descent on 0.5*(x-target)^2
106            else:
107                stale = hist[0]
108                x += 0.01 * (-0.2 * (x-targets) + 0.8 * (A@stale-degree*x))
109                hist.pop(0); hist.append(x.copy())
110        obj = float(np.mean(0.5*(x-targets)**2))
111        dis = float(np.mean((x-x.mean())**2))
112        opt[name] = {'final_mean_local_objective': obj, 'final_disagreement': dis}
113
114    # Aggregate checks deliberately distinguish the first-order approximation
115    # from the exact characteristic-root prediction.
116    mean_approx_err = float(np.mean([r['approx_rel_error'] for r in speed_rows if r['tau'] > 0]))
117    mean_exact_err = float(np.mean([abs(r['observed']-r['exact'])/r['exact'] for r in speed_rows if r['tau'] > 0]))
118    # A useful mechanism check: rate decreases monotonically with k*tau.
119    monotone = all(speed_rows[i]['observed'] >= speed_rows[i+1]['observed'] - 1e-4
120                   for i in range(len(speed_rows)-1) if speed_rows[i]['k'] == speed_rows[i+1]['k'])
121    out = {'seed': SEED, 'step_h': h, 'speed_scaling': speed_rows,
122           'synchronization_ratio_sweep': sync_rows, 'spectral_gap_sweep': gap_rows,
123           'optimizer_comparison': opt,
124           'summary': {'mean_first_order_relative_error': mean_approx_err,
125                       'mean_exact_root_relative_error': mean_exact_err,
126                       'speed_monotone_in_delay': monotone,
127                       'rho_one_disagreement_at': min(sync_rows, key=lambda r: abs(r['rho']-1))}}
128    with open('results.json', 'w') as f: json.dump(out, f, indent=2)
129    print(json.dumps(out['summary'], indent=2))
130
131if __name__ == '__main__': main()