1import json, math, os
  2import numpy as np
  3
  4SEED = 1195
  5rng = np.random.default_rng(SEED)
  6
  7
  8def delayed_quadratic(a, delay, steps=5000, x0=1.0):
  9    # x[t+1] = x[t] - a*x[t-delay], with zero-padded history
 10    x = np.ones(steps + delay + 1) * x0
 11    for t in range(steps):
 12        x[t+delay+1] = x[t+delay] - a*x[t]
 13    tail = np.max(np.abs(x[-500:]))
 14    return tail, np.max(np.abs(x))
 15
 16
 17def theoretical_delay_limit(d):
 18    # Exact Schur stability boundary for z[t+1]=z[t]-a z[t-d].
 19    return 2.0 * math.sin(math.pi / (4*d + 2))
 20
 21
 22def numerical_boundary(d):
 23    lo, hi = 0.0, 3.0
 24    # binary search using a long trajectory; margin avoids classifying boundary oscillations as stable
 25    for _ in range(38):
 26        mid = (lo + hi) / 2
 27        tail, peak = delayed_quadratic(mid, d)
 28        if np.isfinite(tail) and tail < 0.8 and peak < 20:
 29            lo = mid
 30        else:
 31            hi = mid
 32    return lo
 33
 34
 35def lipschitz_check():
 36    # F(x)=.5*lambda*x^2 has L=lambda, so the stated inequality should be exact.
 37    lam = 3.7
 38    xs = np.linspace(-2.0, 2.0, 101)
 39    ys = np.linspace(-2.0, 2.0, 101)
 40    ratios = []
 41    for x in xs:
 42        for y in ys:
 43            if abs(x-y) > 1e-10:
 44                ratios.append((abs(lam*x-lam*y)**2)/(lam**2*(x-y)**2))
 45    return float(np.max(ratios)), float(np.min(ratios))
 46
 47
 48def delay_error_scaling():
 49    # A slowly moving parameter trajectory x_t=v*t makes displacement exactly proportional to d.
 50    # For a quadratic, ||grad(x_t)-grad(x_{t-d})||^2 = L^2 ||x_t-x_{t-d}||^2,
 51    # hence absolute error contribution is predicted to scale as d^2.
 52    lam, v, t = 2.0, 0.013, 1000
 53    ds = np.arange(1, 9)
 54    errs = np.array([(lam*v*d)**2 for d in ds])
 55    slope = np.polyfit(np.log(ds), np.log(errs), 1)[0]
 56    return ds.tolist(), errs.tolist(), float(slope), float(errs[-1]/errs[0])
 57
 58
 59def async_stream(seed=SEED, n_parts=8, workers=8, c=4, steps=180, eta=0.16):
 60    """Small event-driven coded stream on separable logistic regression shards.
 61    Each worker has a replica assignment and random compute delays. The coordinator uses
 62    the newest available gradient per shard once every shard is at most c versions old.
 63    """
 64    rg = np.random.default_rng(seed)
 65    dim, shard_size = 10, 48
 66    true_w = rg.normal(size=dim)
 67    shards = []
 68    for i in range(n_parts):
 69        X = rg.normal(size=(shard_size, dim))
 70        y = (X @ true_w + .35*rg.normal(size=shard_size) > 0).astype(float)
 71        shards.append((X, y))
 72    # two replicas, assigned round-robin with distinct worker speeds
 73    assignments = [[i for i in range(n_parts) if i % workers == w or i % workers == (w-1) % workers]
 74                   for w in range(workers)]
 75    speed = np.exp(rg.normal(np.log(1.0), .38, size=workers))
 76    # queue entries: (completion time, worker, partition, version, gradient)
 77    queue = []
 78    table = [None] * n_parts
 79    beta_hist = [np.zeros(dim)]
 80    clock, t, updates = 0.0, 0, 0
 81    ages, wall = [], []
 82    def grad(beta, i):
 83        X, y = shards[i]
 84        z = np.clip(X @ beta, -30, 30)
 85        p = 1/(1+np.exp(-z))
 86        return X.T @ (p-y) / len(y)
 87    def schedule(w, part, version, beta):
 88        # occasional stragglers, but no unrealistic seconds are needed in normalized time
 89        delay = speed[w] * float(np.exp(rg.normal(0, .35)))
 90        if rg.random() < .07: delay *= 4.5
 91        queue.append((clock+delay, w, part, version, grad(beta, part)))
 92    # initial replicas on version 0
 93    for w in range(workers):
 94        for i in assignments[w]: schedule(w, i, 0, beta_hist[0])
 95    while t < steps and queue:
 96        queue.sort(key=lambda z: z[0])
 97        clock = queue[0][0]
 98        ready = [q for q in queue if q[0] <= clock + 1e-12]
 99        queue = [q for q in queue if q[0] > clock + 1e-12]
100        for _, w, i, v, g in ready:
101            if table[i] is None or v >= table[i][0]: table[i] = (v, g)
102        if all(x is not None and 0 <= t-x[0] < c for x in table):
103            used = [x[0] for x in table]
104            gmean = np.mean([x[1] for x in table], axis=0)
105            beta = beta_hist[-1] - eta*gmean
106            beta_hist.append(beta); ages.extend([t-v for v in used]); wall.append(clock)
107            t += 1; updates += 1
108            # each worker begins one new partition evaluation on the new snapshot
109            for w in range(workers):
110                schedule(w, (t+w) % n_parts, t, beta)
111    # objective on fresh data
112    Xv = rg.normal(size=(800, dim)); yv=(Xv@true_w+.35*rg.normal(size=800)>0).astype(float)
113    z=np.clip(Xv@beta_hist[-1],-30,30)
114    loss=float(np.mean(np.logaddexp(0,z)-yv*z))
115    return dict(loss=loss, updates=updates, wall=float(wall[-1]) if wall else float('inf'),
116                mean_age=float(np.mean(ages)) if ages else float('inf'), max_age=int(max(ages)) if ages else -1)
117
118
119def sync_stream(seed=SEED, steps=180, eta=.16):
120    # Same stochastic problem and delay distribution, but one round waits for all workers.
121    # Approximate synchronous wall time by max of eight independent worker times per round.
122    rg=np.random.default_rng(seed); dim=10; n=8; m=48
123    tw=rg.normal(size=dim); shards=[]
124    for i in range(n):
125        X=rg.normal(size=(m,dim)); y=(X@tw+.35*rg.normal(size=m)>0).astype(float); shards.append((X,y))
126    beta=np.zeros(dim); total=0.
127    def g(i):
128        X,y=shards[i]; z=np.clip(X@beta,-30,30); p=1/(1+np.exp(-z)); return X.T@(p-y)/m
129    for t in range(steps):
130        gs=[g(i) for i in range(n)]; beta-=eta*np.mean(gs,axis=0)
131        delays=rg.lognormal(0,.38,size=8); delays[rg.random(8)<.07]*=4.5; total+=float(np.max(delays))
132    X=rg.normal(size=(800,dim)); y=(X@tw+.35*rg.normal(size=800)>0).astype(float); z=np.clip(X@beta,-30,30)
133    return dict(loss=float(np.mean(np.logaddexp(0,z)-y*z)), updates=steps, wall=total, mean_age=0., max_age=0)
134
135
136def main():
137    stability=[]
138    for d in [0,1,2,3,5,8]:
139        pred=theoretical_delay_limit(d); obs=numerical_boundary(d)
140        stability.append(dict(delay=d, predicted=pred, observed=obs, relative_error=abs(obs-pred)/pred))
141    lip=lipschitz_check(); ds, errs, slope, growth=delay_error_scaling()
142    base=sync_stream(); idea={str(c):async_stream(c=c) for c in [2,4,8]}
143    out={'seed':SEED, 'stability_boundary':stability,
144         'lipschitz_bound_max_ratio':lip[0], 'lipschitz_bound_min_ratio':lip[1],
145         'delay_error_delays':ds, 'delay_error_values':errs, 'delay_error_loglog_slope':slope,
146         'delay_error_growth_d1_to_d8':growth, 'baseline_sync':base, 'pipelined':idea}
147    with open('results.json','w') as f: json.dump(out,f,indent=2)
148    print(json.dumps(out,indent=2))
149
150if __name__=='__main__': main()