1import json, math, os
  2import numpy as np
  3
  4SEED = 1402
  5RNG = np.random.default_rng(SEED)
  6
  7
  8def geom_sensitivity(delta0, rho, K):
  9    return delta0 * sum(rho ** k for k in range(K))
 10
 11
 12def gaussian_proxy(delta, sigma, delta_fail=1e-5):
 13    return delta / sigma * math.sqrt(2.0 * math.log(1.25 / delta_fail))
 14
 15
 16def stability_sweep():
 17    # For x_{k+1}=(1-beta*lambda)x_k, contraction requires |1-beta*lambda|<1.
 18    lam = 1.0
 19    betas = [0.5, 1.5, 2.0, 2.1, 2.5]
 20    rows = []
 21    for beta in betas:
 22        rho = abs(1.0-beta*lam)
 23        x = 1.0
 24        for _ in range(40):
 25            x = (1.0-beta*lam)*x
 26        rows.append({"beta": beta, "predicted_rho": rho,
 27                     "observed_final_abs": abs(x),
 28                     "predicted_contractive": rho < 1.0,
 29                     "observed_decay": abs(x) < 1.0})
 30    return rows
 31
 32def toy_predictions():
 33    # Prediction 1: a contractive recurrence has geometric sensitivity and a finite limit.
 34    delta0 = 1.0
 35    K = 100
 36    rhos = [0.2, 0.5, 0.8, 0.95, 1.05]
 37    rows = []
 38    for rho in rhos:
 39        vals = delta0 * rho ** np.arange(K)
 40        observed = float(vals.sum())
 41        predicted = float(delta0 * (1-rho**K)/(1-rho)) if rho != 1 else float(K)
 42        infinite = delta0/(1-rho) if rho < 1 else None
 43        rows.append({"rho":rho, "observed_sum":observed, "finite_geom_prediction":predicted,
 44                     "infinite_bound":infinite, "ratio_to_bound":(observed/infinite if infinite else None)})
 45
 46    # Prediction 2: 95% of infinite sensitivity mass arrives by log(.05)/log(rho).
 47    k95_rows = []
 48    for rho in [0.5, 0.8, 0.9, 0.95]:
 49        pred = math.log(0.05)/math.log(rho)
 50        vals = rho ** np.arange(1000)
 51        cumulative = np.cumsum(vals)
 52        target = 0.95/(1-rho)
 53        observed = int(np.flatnonzero(cumulative >= target)[0] + 1)
 54        k95_rows.append({"rho":rho, "predicted_k95_steps":pred, "observed_k95_steps":observed,
 55                         "relative_error":abs(observed-pred)/pred})
 56
 57    # Prediction 3: for fixed d, channel privacy proxy scales as alpha^-1/2.
 58    d = 2.0; sigma0 = 0.2; eps = 1e-6; p = 2.0; delta = 0.1
 59    alphas = np.array([0.01, 0.04, 0.16, 0.64])
 60    proxy = np.array([gaussian_proxy(delta, math.sqrt(sigma0**2+a*(d+eps)**p)) for a in alphas])
 61    scaling = proxy[0] / proxy
 62    expected = np.sqrt(alphas/alphas[0])
 63    scaling_rows = [{"alpha":float(a), "observed_proxy":float(q),
 64                     "observed_reduction_vs_first":float(r), "predicted_reduction":float(e)}
 65                    for a,q,r,e in zip(alphas,proxy,scaling,expected)]
 66    return {"contraction_sum":rows, "k95":k95_rows, "alpha_scaling":scaling_rows, "stability_boundary":stability_sweep()}
 67
 68
 69def clip(x, C):
 70    n = np.linalg.norm(x)
 71    return x if n <= C else x * (C/n)
 72
 73
 74def run_federated(mode, beta=0.25, sigma0=0.08, alpha=0.08, p=2.0, rounds=80,
 75                  clients=10, dim=8, seed=1402):
 76    rng = np.random.default_rng(seed)
 77    # Each client has a quadratic target; adjacent data changes client 0's target.
 78    targets = rng.normal(0, 1.0, size=(clients, dim))
 79    targets -= targets.mean(axis=0, keepdims=True)
 80    state = np.zeros(dim)
 81    previous_updates = np.zeros((clients, dim))
 82    losses=[]; proxy_sum=0.; noise_sum=0.; disagreement=[]
 83    C=1.5; eta=0.45; delta_fail=1e-5
 84    for k in range(rounds):
 85        updates=[]; variances=[]; ds=[]
 86        for i in range(clients):
 87            grad = state - targets[i]
 88            z = clip(-eta*grad, C)
 89            d = float(np.linalg.norm(z-previous_updates[i]))
 90            if mode == 'constant':
 91                var = sigma0**2 + alpha*(1e-6)**p
 92            elif mode == 'channel':
 93                var = sigma0**2 + alpha*(d+1e-6)**p
 94            else:
 95                var = 0.0
 96            updates.append(z + rng.normal(0, math.sqrt(var), dim))
 97            variances.append(var); ds.append(d)
 98            # The bound uses clipped per-client sensitivity; use a conservative proxy.
 99            proxy_sum += (float('inf') if var == 0.0 else gaussian_proxy(min(2*C, 0.25), math.sqrt(var), delta_fail))
100        aggregate=np.mean(updates,axis=0)
101        state=(1-beta)*state + beta*(state+aggregate)
102        previous_updates=np.array([clip(-eta*(state-targets[i]), C) for i in range(clients)])
103        losses.append(float(np.mean((state-targets)**2)))
104        noise_sum += float(np.mean(variances)); disagreement.append(float(np.mean(ds)))
105    return {"final_loss":losses[-1], "best_loss":min(losses), "privacy_proxy_sum":proxy_sum,
106            "mean_noise_variance":noise_sum/rounds, "mean_disagreement":float(np.mean(disagreement)),
107            "loss_curve":losses}
108
109
110def main():
111    math_checks = toy_predictions()
112    results={"math_checks":math_checks, "federated":{}}
113    for mode in ['none','constant','channel']:
114        results['federated'][mode]=run_federated(mode)
115    # A separate paired scalar recurrence numerically estimates rho from neighboring states.
116    rho_hat=[]; a=1.0; b=1.001; rho=0.8
117    for _ in range(20):
118        rho_hat.append(abs(rho*b-rho*a)/abs(b-a)); a,b=rho*a,rho*b
119    results['estimated_contraction']={"true_rho":rho,"median_rho_hat":float(np.median(rho_hat))}
120    with open('results.json','w') as f: json.dump(results,f,indent=2)
121    print(json.dumps(results, indent=2))
122
123if __name__ == '__main__': main()