1import json, math
  2from pathlib import Path
  3import numpy as np
  4
  5SEED = 1522
  6
  7
  8def cycle_laplacian(n):
  9    L = 2*np.eye(n)
 10    for i in range(n):
 11        L[i, (i-1) % n] = -1
 12        L[i, (i+1) % n] = -1
 13    return L
 14
 15
 16def run_rollout(L, alpha, gamma, sigma, T, d=3, burn=0, reps=1, seed=0, initial=None):
 17    """Linear contractive recurrent GNN: H[t+1]=alpha(I-gamma L)H[t]+sigma Xi[t]."""
 18    local = np.random.default_rng(seed)
 19    n = L.shape[0]
 20    P = alpha * (np.eye(n) - gamma * L)
 21    vals = []
 22    for _ in range(reps):
 23        H = np.zeros((n, d)) if initial is None else initial.copy()
 24        es = []
 25        for t in range(T):
 26            H = P @ H + sigma * local.standard_normal((n, d))
 27            if t >= burn:
 28                es.append(0.5 * np.trace(H.T @ L @ H))
 29        vals.append(np.mean(es))
 30    return float(np.mean(vals)), float(np.std(vals, ddof=1) / math.sqrt(reps)) if reps > 1 else 0.0
 31
 32
 33def stationary_energy_formula(eigs, alpha, gamma, sigma, d):
 34    q = alpha * (1 - gamma * eigs)
 35    return float(0.5 * d * np.sum(eigs[1:] * sigma**2 / (1-q[1:]**2)))
 36
 37
 38def main():
 39    alpha, gamma, d = 0.82, 0.20, 3
 40    n = 12
 41    L = cycle_laplacian(n)
 42    eigs = np.linalg.eigvalsh(L)
 43    rho = np.max(np.abs(alpha * (1 - gamma * eigs)))
 44    assert rho < 1, rho
 45
 46    # Prediction 1: stationary energy is proportional to sigma^2, and agrees
 47    # with the exact sum over graph Fourier modes.
 48    sigmas = np.array([0.01, 0.02, 0.04, 0.08])
 49    rows_sigma = []
 50    for s in sigmas:
 51        observed, se = run_rollout(L, alpha, gamma, s, T=3500, burn=1000, d=d, reps=80, seed=100+int(s*10000))
 52        predicted = stationary_energy_formula(eigs, alpha, gamma, s, d)
 53        rows_sigma.append({'sigma': float(s), 'sigma2_lambda2': float(s*s*eigs[1]), 'observed': observed, 'predicted': predicted, 'relative_error': abs(observed-predicted)/predicted, 'se': se})
 54    x = sigmas**2
 55    y = np.array([r['observed'] for r in rows_sigma])
 56    slope = float(np.dot(x, y) / np.dot(x, x))
 57    slope_r2 = float(1 - np.sum((y-slope*x)**2) / np.sum((y-y.mean())**2))
 58
 59    # Prediction 2: a pure lambda_2 mode decays as q_2^(2t), giving the
 60    # half-life log(1/2)/log(q_2^2).
 61    v2 = np.linalg.eigh(L)[1][:, 1]
 62    initial = np.outer(v2, np.ones(d))
 63    q2 = alpha * (1 - gamma * eigs[1])
 64    pred_half = math.log(0.5) / math.log(q2*q2)
 65    P = alpha * (np.eye(n) - gamma * L)
 66    H = initial.copy(); energies = []
 67    for t in range(80):
 68        energies.append(0.5*np.trace(H.T@L@H)); H = P@H
 69    target = energies[0]/2
 70    obs_half = next((i for i,e in enumerate(energies) if e <= target), 80)
 71
 72    # Prediction 3: for one fixed graph topology, uniformly scaling edge
 73    # weights increases lambda_2 and increases the stationary energy. This
 74    # avoids confounding gap with graph size/mode count.
 75    rows_gap = []
 76    sigma = 0.04
 77    for scale in [0.25, 0.5, 0.75, 1.0, 1.25, 1.5]:
 78        LL = scale * L
 79        ee = np.linalg.eigvalsh(LL)
 80        observed, se = run_rollout(LL, alpha, gamma, sigma, T=4000, burn=1200, d=d, reps=50, seed=500+int(scale*100))
 81        predicted = stationary_energy_formula(ee, alpha, gamma, sigma, d)
 82        rows_gap.append({'scale': scale, 'lambda2': float(ee[1]), 'observed': observed, 'predicted': predicted, 'se': se})
 83    gap_monotonic = all(rows_gap[i+1]['observed'] > rows_gap[i]['observed'] for i in range(len(rows_gap)-1))
 84
 85    # Baseline comparison: same initial nonconstant state, with and without
 86    # persistent noise, at long horizon.
 87    baseline_long, _ = run_rollout(L, alpha, gamma, 0.0, T=100, burn=99, d=d, reps=1, seed=7, initial=initial)
 88    idea_energy, idea_se = run_rollout(L, alpha, gamma, 0.04, T=3500, burn=1000, d=d, reps=80, seed=777, initial=initial)
 89
 90    result = {
 91        'setup': {'n': n, 'd': d, 'alpha': alpha, 'gamma': gamma, 'contraction_rho': float(rho), 'lambda2': float(eigs[1]), 'seed': SEED},
 92        'prediction_1_sigma_squared_scaling': {'rows': rows_sigma, 'fit_slope_energy_per_sigma2': slope, 'R2_through_origin': slope_r2},
 93        'prediction_2_deterministic_decay': {'q_lambda2': float(q2), 'predicted_half_life_steps': pred_half, 'observed_first_half_step': int(obs_half), 'initial_energy': float(energies[0]), 'energy_at_20': float(energies[20])},
 94        'prediction_3_fixed_topology_gap_sweep': {'rows': rows_gap, 'strictly_monotone_observed': gap_monotonic},
 95        'baseline_vs_persistent_noise': {'deterministic_energy_at_100': baseline_long, 'persistent_noise_stationary_energy': idea_energy, 'persistent_noise_se': idea_se},
 96    }
 97    Path('results.json').write_text(json.dumps(result, indent=2))
 98    print(json.dumps(result, indent=2))
 99
100if __name__ == '__main__':
101    main()