import json, math, random from pathlib import Path import numpy as np SEED = 2870 rng = np.random.default_rng(SEED) def make_incidence(n=7, edges=None): if edges is None: edges = [(0,1),(1,2),(2,3),(3,4),(4,5),(5,6),(0,6),(1,5),(2,6)] B = np.zeros((n, len(edges)), dtype=np.float64) for j, (tail, head) in enumerate(edges): B[tail, j] = -1.0 B[head, j] = 1.0 return B, edges def core_sweeps(): B, edges = make_incidence() n, m = B.shape ones = np.ones(n) # Prediction 1: arbitrary parallel channels remain exactly conservative. channel_rows = [] for k in [1, 2, 4, 8, 16, 32]: Pk = rng.normal(size=(m, k)) P = Pk.sum(axis=1) residual = float(abs(ones @ (B @ P))) channel_rows.append({"K": k, "abs_total_internal_residual": residual}) # Prediction 2: if an unconstrained construction adds leakage with amplitude # lambda, the mass residual is linear in lambda. Here r is fixed, so this is # a clean quantitative slope test rather than a noisy training result. P = rng.normal(size=m) r = rng.normal(size=n) lambdas = [0.0, 0.01, 0.05, 0.1, 0.25, 0.5, 1.0, 2.0] leakage_rows = [] for lam in lambdas: q_unconstrained = B @ P + lam * r observed = float(abs(ones @ q_unconstrained)) predicted = float(lam * abs(ones @ r)) leakage_rows.append({"lambda": lam, "observed_abs_residual": observed, "predicted_abs_residual": predicted}) nonzero = [(x["lambda"], x["observed_abs_residual"]) for x in leakage_rows if x["lambda"] > 0] slope = float(sum(a*b for a,b in nonzero) / sum(a*a for a,b in nonzero)) expected_slope = float(abs(ones @ r)) # Prediction 3: with internal-only dynamics, conservative Euler updates # preserve total mass for every rollout length; constant leakage accumulates # linearly in the number of steps. dt, eps = 0.05, 0.02 x0 = rng.normal(size=n) flow = rng.normal(size=m) q_internal = B @ flow steps_list = [1, 2, 5, 10, 20, 40, 80] drift_rows = [] for steps in steps_list: xc = x0.copy() xu = x0.copy() for _ in range(steps): xc += dt * q_internal xu += dt * (q_internal + eps * np.ones(n)) cdrift = float(abs(ones @ (xc - x0))) udrift = float(abs(ones @ (xu - x0))) predicted_u = float(steps * dt * eps * n) drift_rows.append({"steps": steps, "conservative_abs_mass_drift": cdrift, "unconstrained_abs_mass_drift": udrift, "predicted_unconstrained_drift": predicted_u}) return { "incidence_shape": [n, m], "incidence_column_sum_max_abs": float(np.max(np.abs(ones @ B))), "channel_sweep": channel_rows, "leakage_sweep": leakage_rows, "leakage_fitted_slope": slope, "leakage_predicted_slope": expected_slope, "rollout_sweep": drift_rows, } def toy_model_comparison(): # Small supervised graph-flow task: true interaction is a sum of two # domain-specific edge channels. The conservative model has two scalar # channel regressors and aggregates only through B. The baseline predicts # both endpoint updates independently from the composite edge features. try: import torch torch.manual_seed(SEED) np.random.seed(SEED) device = torch.device("cuda" if torch.cuda.is_available() else "cpu") B_np, edges = make_incidence(n=6, edges=[(0,1),(1,2),(2,3),(3,4),(4,5),(0,5),(1,4)]) B = torch.tensor(B_np, dtype=torch.float32, device=device) n, m = B_np.shape N = 1200 x = torch.randn(N, n, device=device) u = torch.randn(N, m, 2, device=device) tail = torch.tensor([e[0] for e in edges], device=device) head = torch.tensor([e[1] for e in edges], device=device) d = x[:, tail] - x[:, head] target_p = 0.8*torch.tanh(d) + 0.25*u[:,:,0] + 0.35*torch.sin(d*u[:,:,1]) target_q = torch.einsum('nm,bm->bn', B, target_p) train = slice(0, 900); test = slice(900, 1200) class Conservative(torch.nn.Module): def __init__(self): super().__init__() self.c0 = torch.nn.Sequential(torch.nn.Linear(3, 16), torch.nn.Tanh(), torch.nn.Linear(16, 1)) self.c1 = torch.nn.Sequential(torch.nn.Linear(3, 16), torch.nn.Tanh(), torch.nn.Linear(16, 1)) def forward(self, xx, uu): dd = xx[:,tail] - xx[:,head] z = torch.stack((dd, uu[:,:,0], uu[:,:,1]), dim=-1) p = self.c0(z).squeeze(-1) + self.c1(z).squeeze(-1) return torch.einsum('nm,bm->bn', B, p) class Unconstrained(torch.nn.Module): def __init__(self): super().__init__() self.edge = torch.nn.Sequential(torch.nn.Linear(4, 24), torch.nn.Tanh(), torch.nn.Linear(24, 2)) def forward(self, xx, uu): dd = xx[:,tail] - xx[:,head] z = torch.stack((xx[:,tail], xx[:,head], uu[:,:,0], uu[:,:,1]), dim=-1) msg = self.edge(z) out = torch.zeros((xx.shape[0], n), device=xx.device) out.index_add_(1, tail, msg[:,:,0]) out.index_add_(1, head, msg[:,:,1]) return out models = {"conservative": Conservative().to(device), "unconstrained": Unconstrained().to(device)} results = {} for name, model in models.items(): opt = torch.optim.Adam(model.parameters(), lr=3e-3) for step in range(500): opt.zero_grad() pred = model(x[train], u[train]) loss = ((pred-target_q[train])**2).mean() loss.backward(); opt.step() with torch.no_grad(): pred = model(x[test], u[test]) err = ((pred-target_q[test])**2).mean().sqrt().item() residual = (pred.sum(dim=1)).abs().mean().item() # internal conservation is exact for conservative architecture; # baseline residual is measured independently of task error. results[name] = {"test_q_rmse": err, "mean_abs_mass_residual": residual, "parameters": sum(p.numel() for p in model.parameters())} return {"device": str(device), "results": results} except Exception as exc: return {"device": "cpu-fallback", "error": repr(exc)} def main(): out = {"seed": SEED, "core_math": core_sweeps(), "toy_model": toy_model_comparison()} Path("results.json").write_text(json.dumps(out, indent=2)) print(json.dumps(out, indent=2)) if __name__ == "__main__": main()