1import json
  2import numpy as np
  3
  4SEED = 2104
  5rng = np.random.default_rng(SEED)
  6
  7
  8def ring_matrix(n):
  9    A = np.zeros((n, n))
 10    for i in range(n):
 11        A[i, i] = 0.5
 12        A[i, (i - 1) % n] += 0.25
 13        A[i, (i + 1) % n] += 0.25
 14    return A
 15
 16
 17def info_update(J, h, H, y, R):
 18    Ri = np.linalg.inv(R)
 19    return J + H.T @ Ri @ H, h + H.T @ Ri @ y
 20
 21
 22def gramian_sweep():
 23    # Two agents each see one coordinate. q controls how often agent 2 reports.
 24    # With unit noise and identity dynamics, lambda_min(G)=min(1,q).
 25    rows = []
 26    for q in [0.0, 0.05, 0.1, 0.25, 0.5, 1.0, 2.0]:
 27        H1 = np.array([[1., 0.]])
 28        H2 = np.array([[0., 1.]])
 29        G = H1.T @ H1 + q * (H2.T @ H2)
 30        observed = np.linalg.eigvalsh(G).min()
 31        predicted = min(1.0, q)
 32        rows.append((q, observed, predicted, abs(observed-predicted)))
 33    return rows
 34
 35
 36def noise_scaling_sweep():
 37    # Collective information from complementary coordinate sensors.
 38    # For sigma^2 << prior variance, P_ii is predicted to scale as sigma^2.
 39    rows = []
 40    prior = np.eye(2) * 1e6
 41    for sigma in [0.1, 0.2, 0.5, 1.0, 2.0, 4.0]:
 42        J = np.linalg.inv(prior)
 43        J += np.diag([1/sigma**2, 1/sigma**2])
 44        P = np.linalg.inv(J)
 45        observed = np.linalg.eigvalsh(P).min()
 46        predicted = sigma**2
 47        rows.append((sigma, observed, predicted, observed/predicted))
 48    return rows
 49
 50
 51def diffusion_sweep():
 52    # Natural parameters initialized at different nodes; consensus error obeys
 53    # ||A^t x - mean(x)|| <= rho^t ||x-mean(x)|| for symmetric A.
 54    A = ring_matrix(4)
 55    eig = np.linalg.eigvalsh(A)
 56    rho = np.max(np.abs(eig[np.abs(eig-1) > 1e-9]))
 57    x = np.array([2., -1., 3., 0.])
 58    mean = x.mean()
 59    e0 = np.linalg.norm(x - mean)
 60    rows = []
 61    for t in [1, 2, 3, 5, 8, 12, 20]:
 62        xt = np.linalg.matrix_power(A, t) @ x
 63        ratio = np.linalg.norm(xt - mean) / e0
 64        predicted = rho**t
 65        rows.append((t, ratio, predicted, ratio/predicted if predicted else np.nan))
 66    return rho, rows
 67
 68
 69def fusion_trial(n_trials=400):
 70    # Two asynchronous agents observe complementary coordinates. Agent 1 reports
 71    # every step, agent 2 every other step. Each communication uses a ring mix.
 72    # Compare natural-parameter fusion with arithmetic mean of local posterior means.
 73    A = ring_matrix(2)
 74    H = [np.array([[1., 0.]]), np.array([[0., 1.]])]
 75    sigma = 0.7
 76    R = np.array([[sigma**2]])
 77    priorJ = np.eye(2) * 0.05
 78    errs_info, errs_mean = [], []
 79    for _ in range(n_trials):
 80        z = rng.normal(size=2)
 81        Js = [priorJ.copy(), priorJ.copy()]
 82        hs = [np.zeros(2), np.zeros(2)]
 83        ms = [np.zeros(2), np.zeros(2)]
 84        for k in range(12):
 85            for l in range(2):
 86                if l == 0 or k % 2 == 1:
 87                    y = H[l] @ z + rng.normal(scale=sigma, size=(1,))
 88                    Js[l], hs[l] = info_update(Js[l], hs[l], H[l], y, R)
 89                    ms[l] = np.linalg.solve(Js[l], hs[l])
 90            # diffuse natural parameters; arithmetic baseline diffuses means only
 91            Jnew = [sum(A[l,j] * Js[j] for j in range(2)) for l in range(2)]
 92            hnew = [sum(A[l,j] * hs[j] for j in range(2)) for l in range(2)]
 93            Js, hs = Jnew, hnew
 94            info_m = np.linalg.solve(Js[0], hs[0])
 95            mean_m = sum(A[0,j] * ms[j] for j in range(2))
 96            errs_info.append(np.mean((info_m-z)**2))
 97            errs_mean.append(np.mean((mean_m-z)**2))
 98    return float(np.mean(errs_info)), float(np.mean(errs_mean))
 99
100
101def main():
102    gram = gramian_sweep()
103    noise = noise_scaling_sweep()
104    rho, diffusion = diffusion_sweep()
105    info_err, mean_err = fusion_trial()
106    out = {
107        "seed": SEED,
108        "predictions": {
109            "gramian": "lambda_min=min(1,q), with positivity iff q>0",
110            "noise_scaling": "small-posterior covariance eigenvalue approximately sigma^2",
111            "diffusion": "consensus disagreement ratio is bounded by rho^t"
112        },
113        "gramian": [{"q":q,"observed":o,"predicted":p,"abs_error":e} for q,o,p,e in gram],
114        "noise_scaling": [{"sigma":s,"observed":o,"predicted":p,"ratio":r} for s,o,p,r in noise],
115        "diffusion": {"rho":rho,"rows":[{"rounds":t,"observed_ratio":o,"predicted_bound":p,"ratio_to_bound":r} for t,o,p,r in diffusion]},
116        "fusion_mse": {"information_diffusion":info_err,"arithmetic_mean":mean_err},
117    }
118    with open("results.json", "w") as f:
119        json.dump(out, f, indent=2)
120    print(json.dumps(out, indent=2))
121
122if __name__ == "__main__":
123    main()