1import json
  2import math
  3import numpy as np
  4
  5SEED = 2069
  6rng = np.random.default_rng(SEED)
  7
  8
  9def companion(a, alpha, delay):
 10    """Companion matrix for x[t+1] = a*x[t] + alpha*x[t-delay]."""
 11    n = delay + 1
 12    C = np.zeros((n, n), dtype=complex if np.iscomplexobj(alpha) else float)
 13    C[0, 0] = a
 14    C[0, -1] = alpha
 15    C[1:, :-1] = np.eye(delay)
 16    return C
 17
 18
 19def rho(a, alpha, delay):
 20    return float(np.max(np.abs(np.linalg.eigvals(companion(a, alpha, delay)))))
 21
 22
 23def empirical_exponent(a, alpha, delay, steps=5000):
 24    C = companion(a, alpha, delay)
 25    v = rng.normal(size=delay + 1) + 1j * rng.normal(size=delay + 1)
 26    v /= np.linalg.norm(v)
 27    logs = []
 28    for _ in range(steps):
 29        v = C @ v
 30        n = np.linalg.norm(v)
 31        logs.append(math.log(max(n, 1e-300)))
 32        v /= max(n, 1e-300)
 33    # discard transient; renormalization makes this a Lyapunov estimate
 34    return float(np.mean(logs[500:]))
 35
 36
 37def nontrivial_eigenvalues(W):
 38    vals = np.linalg.eigvals(W)
 39    # For row-stochastic W, remove the eigenvalue nearest one.
 40    j = int(np.argmin(np.abs(vals - 1)))
 41    return np.delete(vals, j)
 42
 43
 44def disagreement_rollout(W, a, sigma, delay, steps=180):
 45    n = W.shape[0]
 46    # History is a transverse random perturbation around a nonzero synchronized value.
 47    hist = np.zeros((delay + 1, n))
 48    hist[:] = 1.0 + 0.03 * rng.normal(size=(delay + 1, n))
 49    hist -= hist.mean(axis=1, keepdims=True)
 50    ds = []
 51    for t in range(steps):
 52        x = hist[-1]
 53        delayed = hist[0]
 54        y = a * x + sigma * W @ delayed
 55        ds.append(float(np.mean((x - x.mean()) ** 2)))
 56        hist = np.concatenate([hist[1:], y[None]], axis=0)
 57    tail = np.array(ds[20:])
 58    slope = float(np.polyfit(np.arange(len(tail)), np.log(tail + 1e-30), 1)[0])
 59    return ds, slope
 60
 61
 62def make_symmetric_ring(n):
 63    W = np.zeros((n, n))
 64    for i in range(n):
 65        W[i, (i - 1) % n] = 0.5
 66        W[i, (i + 1) % n] = 0.5
 67    return W
 68
 69
 70def make_directed_heterogeneous(n, trials=3000):
 71    # Directed, nonreciprocal, heterogeneous indegrees; retain row normalization and W1=1.
 72    best, best_score = None, float('inf')
 73    for _ in range(trials):
 74        W = rng.dirichlet(0.22 * np.ones(n), size=n)
 75        vals = nontrivial_eigenvalues(W)
 76        score = float(np.max(np.abs(vals)))
 77        if score < best_score:
 78            best, best_score = W, score
 79    return best
 80
 81
 82def main():
 83    a = 0.72
 84    # Prediction 1: exact transition at rho=1 for every delay.
 85    boundary_rows = []
 86    for d in [0, 1, 2, 4, 7]:
 87        # bisection on positive real alpha
 88        lo, hi = 0.0, 2.0
 89        for _ in range(55):
 90            mid = (lo + hi) / 2
 91            if rho(a, mid, d) < 1:
 92                lo = mid
 93            else:
 94                hi = mid
 95        alpha_star = (lo + hi) / 2
 96        boundary_rows.append({"delay": d, "predicted_alpha_rho1": alpha_star,
 97                              "rho_at_boundary": rho(a, alpha_star, d)})
 98
 99    # Prediction 2: measured finite-rollout exponent equals log spectral radius.
100    exponent_rows = []
101    for d in [1, 2, 4]:
102        for alpha in [0.05, 0.15, 0.30, 0.50]:
103            r = rho(a, alpha, d)
104            e = empirical_exponent(a, alpha, d)
105            exponent_rows.append({"delay": d, "alpha": alpha, "log_rho": math.log(r),
106                                  "measured_exponent": e, "abs_error": abs(e-math.log(r))})
107
108    # Prediction 3: increasing delay shrinks the stable positive-alpha interval.
109    delay_rows = []
110    for d in range(0, 9):
111        lo, hi = 0.0, 2.0
112        for _ in range(50):
113            mid = (lo + hi) / 2
114            if rho(a, mid, d) < 1: lo = mid
115            else: hi = mid
116        delay_rows.append({"delay": d, "stable_alpha_max": (lo+hi)/2})
117
118    n = 8
119    W_base = make_symmetric_ring(n)
120    W_idea = make_directed_heterogeneous(n)
121    vals_base = nontrivial_eigenvalues(W_base)
122    vals_idea = nontrivial_eigenvalues(W_idea)
123    sigma = 0.42
124    delay = 2
125    ds_b, slope_b = disagreement_rollout(W_base, a, sigma, delay)
126    ds_i, slope_i = disagreement_rollout(W_idea, a, sigma, delay)
127    network = {
128        "a": a, "sigma": sigma, "delay": delay,
129        "baseline_max_transverse_abs_eigenvalue": float(np.max(np.abs(vals_base))),
130        "idea_max_transverse_abs_eigenvalue": float(np.max(np.abs(vals_idea))),
131        "baseline_predicted_worst_log_rho": max(math.log(rho(a, sigma*v, delay)) for v in vals_base),
132        "idea_predicted_worst_log_rho": max(math.log(rho(a, sigma*v, delay)) for v in vals_idea),
133        "baseline_disagreement_log_slope": slope_b,
134        "idea_disagreement_log_slope": slope_i,
135        "baseline_final_disagreement": ds_b[-1],
136        "idea_final_disagreement": ds_i[-1],
137        "baseline_W": W_base.tolist(), "idea_W": W_idea.tolist()
138    }
139    result = {"seed": SEED, "predictions": {
140        "boundary_rho_equals_one": boundary_rows,
141        "exponent_log_rho": exponent_rows,
142        "delay_stability_boundary": delay_rows}, "network_comparison": network}
143    with open("results.json", "w") as f: json.dump(result, f, indent=2)
144    print(json.dumps({"boundary": boundary_rows, "max_exponent_error": max(x["abs_error"] for x in exponent_rows),
145                      "delay_boundary": delay_rows, "network": {k:v for k,v in network.items() if not k.endswith("_W")}}, indent=2))
146
147if __name__ == "__main__": main()