1"""Directed-path synchronization coupling: math checks and a tiny experiment.
  2
  3Convention: C has C[i,j]=a_ij for j != i and C[i,i]=-sum_j a_ij,
  4so xdot = f(x) - eta*C_pos?  We instead construct the positive
  5row Laplacian B with B[i,i]=sum a and B[i,j]=-a, and use
  6xdot = f(x) - eta B x.  Every non-root node receives its parent.
  7"""
  8import json
  9import numpy as np
 10
 11SEED = 2470
 12rng = np.random.default_rng(SEED)
 13
 14
 15def directed_star(n, weight=1.0):
 16    # n-1 directed parent paths (root 0 -> each child), represented as
 17    # receiver rows: child receives the root state.
 18    B = np.zeros((n, n))
 19    for i in range(1, n):
 20        B[i, i] = weight
 21        B[i, 0] = -weight
 22    return B
 23
 24
 25def dense_symmetric(n, weight=1.0):
 26    # Complete balanced diffusive graph, same row-sum normalization.
 27    B = np.full((n, n), -weight / (n - 1))
 28    np.fill_diagonal(B, weight)
 29    return B
 30
 31
 32def chain(n, weights):
 33    # Imbalanced directed path: each node receives only its predecessor.
 34    B = np.zeros((n, n))
 35    for i, w in enumerate(weights, start=1):
 36        B[i, i] = w
 37        B[i, i - 1] = -w
 38    return B
 39
 40
 41def graph_gain(B, tol=1e-9):
 42    """Conservative spectral gain: smallest positive real Laplacian mode."""
 43    ev = np.linalg.eigvals(B)
 44    positive = [z.real for z in ev if z.real > tol and abs(z.imag) < 1e-7]
 45    if not positive:
 46        # For general directed matrices, use the minimum positive real part.
 47        positive = [z.real for z in ev if z.real > tol]
 48    return float(min(positive))
 49
 50
 51def simulate(B, lam, eta, dt=0.002, T=5.0, x0=None):
 52    n = B.shape[0]
 53    steps = int(T / dt)
 54    x = np.array(x0 if x0 is not None else rng.normal(size=n), dtype=float)
 55    # For directed graphs use the left-null weighted consensus coordinate.
 56    left = np.linalg.svd(B.T)[2][-1]
 57    left = left / np.sum(left)
 58    V = []
 59    for _ in range(steps + 1):
 60        mean = np.sum(left * x)
 61        V.append(0.5 * np.sum((x - mean) ** 2))
 62        x = x + dt * (lam * x - eta * (B @ x))
 63    return np.asarray(V)
 64
 65
 66def slope(V, dt, tail_fraction=0.6):
 67    k = np.arange(len(V))
 68    first = int(len(V) * (1-tail_fraction))
 69    y = np.log(np.maximum(V, 1e-300))
 70    return float(np.polyfit(k[first:] * dt, y[first:], 1)[0])
 71
 72
 73def threshold_sweep(B, lam, gamma, etas):
 74    # Stable means final V is below its initial V, with a substantial margin.
 75    x0 = np.linspace(-1.0, 1.0, B.shape[0])
 76    rows = []
 77    for eta in etas:
 78        V = simulate(B, lam, eta, x0=x0)
 79        rows.append({"eta": float(eta), "final_over_initial": float(V[-1]/V[0]),
 80                     "slope": slope(V, .002)})
 81    return rows
 82
 83
 84def main():
 85    n = 8
 86    lam = 0.8                 # exact Lipschitz constant of f(x)=lam*x
 87    star = directed_star(n)
 88    dense = dense_symmetric(n)
 89    imbalanced = chain(n, [0.35, 0.55, 0.75, 0.95, 1.15, 1.35, 1.55])
 90    graphs = {"directed_n_minus_1_star": star, "dense_symmetric": dense,
 91              "severely_imbalanced_path": imbalanced}
 92    gains = {name: graph_gain(B) for name, B in graphs.items()}
 93    # Prediction 1: boundary eta_c = L/gamma.
 94    boundary = {name: lam/g for name, g in gains.items()}
 95    etas = np.array([0.3, 0.6, 0.79, 0.81, 1.0, 1.4, 2.0])
 96    sweeps = {name: threshold_sweep(B, lam, gains[name], etas)
 97              for name, B in graphs.items()}
 98
 99    # Prediction 2: rate d log(V)/dt approximately -2(eta*gamma-L).
100    rate_rows = []
101    for eta in [1.0, 1.4, 2.0]:
102        V = simulate(star, lam, eta, x0=np.linspace(-1, 1, n))
103        observed = slope(V, .002)
104        predicted = -2*(eta*gains["directed_n_minus_1_star"] - lam)
105        rate_rows.append({"eta": eta, "predicted": predicted, "observed": observed,
106                          "relative_error": abs(observed-predicted)/max(abs(predicted),1e-9)})
107
108    # Prediction 3: threshold scales linearly with L (gamma fixed).
109    scaling = []
110    for L in [0.2, 0.5, 0.8, 1.2]:
111        g = gains["directed_n_minus_1_star"]
112        candidate = L/g
113        # Estimate the threshold from the first negative fitted log-V slope.
114        candidates = np.linspace(max(.02, candidate*.5), candidate*1.5, 61)
115        found = None
116        for eta in candidates:
117            V = simulate(star, L, eta, x0=np.linspace(-1, 1, n))
118            if slope(V, .002) < 0:
119                found = float(eta); break
120        scaling.append({"L": L, "predicted_eta_c": candidate,
121                        "observed_negative_slope_eta": found,
122                        "relative_boundary_error": (abs(found-candidate)/candidate
123                            if found is not None else None)})
124
125    # Standard baselines at the same lambda and eta.
126    x0 = np.linspace(-1, 1, n)
127    baseline = simulate(np.zeros((n,n)), lam, 1.0, x0=x0)
128    idea = simulate(star, lam, 1.0, x0=x0)
129    dense_run = simulate(dense, lam, 1.0, x0=x0)
130    out = {"seed": SEED, "L": lam, "graph_gains": gains,
131           "predicted_thresholds": boundary, "threshold_sweep": sweeps,
132           "rate_test": rate_rows, "L_scaling": scaling,
133           "comparison_eta_1": {
134               "uncoupled_final_over_initial": float(baseline[-1]/baseline[0]),
135               "directed_path_final_over_initial": float(idea[-1]/idea[0]),
136               "dense_symmetric_final_over_initial": float(dense_run[-1]/dense_run[0]),
137               "uncoupled_slope": slope(baseline,.002),
138               "directed_path_slope": slope(idea,.002),
139               "dense_symmetric_slope": slope(dense_run,.002)}}
140    print(json.dumps(out, indent=2))
141
142if __name__ == "__main__":
143    main()