1import json, math, random
  2from pathlib import Path
  3import numpy as np
  4
  5SEED = 2740
  6rng = np.random.default_rng(SEED)
  7
  8
  9def disagreement(x):
 10    z = x - x.mean(axis=0, keepdims=True)
 11    return float(np.mean(z*z))
 12
 13
 14def graph_laplacian(n, p, seed):
 15    r = np.random.default_rng(seed)
 16    A = (r.random((n,n)) < p).astype(float)
 17    A = np.triu(A, 1); A = A + A.T
 18    # Ensure connectedness by adding a ring.
 19    for i in range(n):
 20        A[i, (i+1)%n] = A[(i+1)%n, i] = 1.0
 21    return np.diag(A.sum(1)) - A
 22
 23
 24def toy_sweeps():
 25    # Prediction 1: exact complete-consensus contraction D'/D=(1-alpha)^2.
 26    x = rng.normal(size=(30, 4))
 27    alphas = np.linspace(0.0, 1.8, 19)
 28    complete = []
 29    for a in alphas:
 30        y = x + a*(x.mean(0, keepdims=True)-x)
 31        complete.append(disagreement(y)/disagreement(x))
 32    pred1 = (1-alphas)**2
 33    err1 = float(np.max(np.abs(np.asarray(complete)-pred1)))
 34
 35    # Prediction 2: Laplacian diffusion is stable iff alpha*lambda_max<2.
 36    L = graph_laplacian(30, .16, SEED)
 37    ev = np.linalg.eigvalsh(L)
 38    lmax, l2 = float(ev[-1]), float(ev[1])
 39    # Use the highest-frequency eigenmode; arbitrary x can hide the operator boundary.
 40    vmax = np.linalg.eigh(L)[1][:,-1:]
 41    agrid = np.linspace(0, 2.4/lmax, 97)
 42    lap_ratios = []
 43    for a in agrid:
 44        y = vmax - a*L@vmax
 45        lap_ratios.append(disagreement(y)/disagreement(vmax))
 46    observed_boundary = float(agrid[np.where(np.asarray(lap_ratios) > 1)[0][0]]) if np.any(np.asarray(lap_ratios)>1) else float('nan')
 47    predicted_boundary = 2/lmax
 48    # Prediction 3: along the slow Fiedler mode, one-step ratio is (1-alpha*l2)^2.
 49    fiedler = np.linalg.eigh(L)[1][:,1:2]
 50    f_ratios=[]; f_pred=[]
 51    for a in alphas:
 52        q=fiedler + 0.0
 53        f_ratios.append(disagreement(q-a*L@q)/disagreement(q))
 54        f_pred.append((1-a*l2)**2)
 55    err3=float(np.max(np.abs(np.asarray(f_ratios)-np.asarray(f_pred))))
 56
 57    return {
 58        'complete_consensus': {'alpha': alphas.tolist(), 'observed_ratio': complete,
 59                               'predicted_ratio_(1-alpha)^2': pred1.tolist(), 'max_abs_error': err1},
 60        'laplacian_diffusion': {'lambda2': l2, 'lambda_max': lmax,
 61                                'predicted_boundary_2/lambda_max': predicted_boundary,
 62                                'observed_first_grid_ratio_gt_1': observed_boundary,
 63                                'alpha_grid': agrid.tolist(), 'ratios': lap_ratios},
 64        'fiedler_mode': {'predicted_ratio_(1-alpha*lambda2)^2': f_pred,
 65                         'observed_ratio': f_ratios, 'max_abs_error': err3}
 66    }
 67
 68
 69def normalized_adj(A):
 70    d=A.sum(1)+1e-6
 71    return A/np.sqrt(d[:,None]*d[None,:])
 72
 73
 74def learned_outage_test():
 75    # Tiny scalar node task: target is local feature plus graph-wide mean.
 76    # Train on random edge masks; evaluate on a held-out, much sparser mask.
 77    import torch
 78    torch.manual_seed(SEED); np.random.seed(SEED)
 79    device='cuda' if torch.cuda.is_available() else 'cpu'
 80    try:
 81        n,d=24,5
 82        base=graph_laplacian(n,.22,SEED+4); baseA=np.diag(np.diag(base))-base
 83        X=torch.tensor(np.random.randn(n,d),dtype=torch.float32,device=device)
 84        target=(X[:,0:1] + .7*X.mean(0,keepdim=True).repeat(n,1)[:,0:1]).detach()
 85        class Net(torch.nn.Module):
 86            def __init__(self, consensus):
 87                super().__init__(); self.consensus=consensus
 88                self.w=torch.nn.Linear(d,12); self.out=torch.nn.Linear(12,1)
 89                self.msg=torch.nn.Linear(12,12,bias=False)
 90            def forward(self,A):
 91                H=torch.relu(self.w(X)); S=torch.tensor(normalized_adj(A),dtype=torch.float32,device=device)
 92                for _ in range(4):
 93                    m=S@H
 94                    H=H + torch.sigmoid(m.mean(1,keepdim=True))*self.msg(m)
 95                    if self.consensus: H=H + .10*(H.mean(0,keepdim=True)-H)
 96                    H=torch.relu(H)
 97                return self.out(H),H
 98        models=[Net(False),Net(True)]
 99        results={}
100        for model in models:
101            model.to(device); opt=torch.optim.Adam(model.parameters(),lr=.015)
102            for step in range(500):
103                mask=(np.random.random((n,n))<.75).astype(float); mask=np.triu(mask,1); mask=mask+mask.T
104                mask=np.maximum(mask,baseA)
105                pred,_=model(mask); loss=((pred-target)**2).mean()
106                opt.zero_grad(); loss.backward(); opt.step()
107            vals=[]; dis=[]
108            for frac in [.75,.55,.35,.20]:
109                errs=[]; ds=[]
110                for k in range(30):
111                    mask=(np.random.random((n,n))<frac).astype(float); mask=np.triu(mask,1); mask=mask+mask.T
112                    mask=np.maximum(mask, np.eye(n)*0)
113                    with torch.no_grad(): pred,H=model(mask)
114                    errs.append(float(((pred-target)**2).mean().cpu())); ds.append(disagreement(H.cpu().numpy()))
115                vals.append(float(np.mean(errs))); dis.append(float(np.mean(ds)))
116            results['consensus' if model.consensus else 'baseline']={'mse_by_edge_keep':vals,'representation_D':dis}
117        results['device']=device
118        return results
119    except Exception as e:
120        return {'error':repr(e),'device':'cpu_fallback'}
121
122
123def main():
124    out={'seed':SEED,'toy':toy_sweeps(),'learned_outage':learned_outage_test()}
125    Path('results.json').write_text(json.dumps(out,indent=2))
126    t=out['toy']; print(json.dumps({'complete_max_error':t['complete_consensus']['max_abs_error'],
127      'lap_pred_boundary':t['laplacian_diffusion']['predicted_boundary_2/lambda_max'],
128      'lap_observed_boundary':t['laplacian_diffusion']['observed_first_grid_ratio_gt_1'],
129      'fiedler_max_error':t['fiedler_mode']['max_abs_error'],
130      'learned':out['learned_outage']},indent=2))
131if __name__=='__main__': main()