1"""Minimal numerical test of a power-preserving formation GNN.
  2
  3Run with:
  4  /home/maxwelhelp/main/bin/python3 power_formation_gnn.py
  5"""
  6import json
  7import numpy as np
  8
  9
 10def formation_matrix(n, d, edges):
 11    B = np.zeros((len(edges), n))
 12    for k, (i, j) in enumerate(edges):
 13        B[k, i] = 1.0
 14        B[k, j] = -1.0
 15    return np.kron(B, np.eye(d))
 16
 17
 18def port_hamiltonian_matrices(F, damping_node=0.15, damping_edge=0.15):
 19    nd, me = F.shape
 20    J = np.block([[np.zeros((nd, nd)), -F.T], [F, np.zeros((me, me))]])
 21    R = np.diag(np.r_[np.full(nd, damping_node), np.full(me, damping_edge)])
 22    return J, R
 23
 24
 25def energy(x, mass):
 26    return 0.5 * np.sum(x * x / mass)
 27
 28
 29def verify_math(rng):
 30    F = formation_matrix(4, 2, [(0, 1), (1, 2), (2, 3), (0, 3)])
 31    J, R = port_hamiltonian_matrices(F)
 32    mass = np.ones(J.shape[0])
 33    x = rng.normal(size=J.shape[0])
 34    grad = x / mass
 35    interconnection_power = float(grad @ J @ grad)
 36    dissipation_power = float(grad @ R @ grad)
 37    dx = (J - R) @ grad
 38    analytic_dh = float(grad @ dx)
 39    eps = 1e-7
 40    finite_difference = (energy(x + eps * dx, mass) - energy(x, mass)) / eps
 41    return {
 42        "skew_error": float(np.max(np.abs(J + J.T))),
 43        "power_cancellation_abs": abs(interconnection_power),
 44        "dissipation_power": dissipation_power,
 45        "analytic_dH": analytic_dh,
 46        "finite_difference_dH": float(finite_difference),
 47        "derivative_match_abs": abs(analytic_dh - finite_difference),
 48    }
 49
 50
 51def simulate_ph(J, R, x0, dt, steps):
 52    x = x0.copy()
 53    hist = [energy(x, np.ones_like(x))]
 54    for _ in range(steps):
 55        x = x + dt * ((J - R) @ x)
 56        hist.append(energy(x, np.ones_like(x)))
 57    return np.asarray(hist), x
 58
 59
 60def stability_sweep(rng):
 61    F = formation_matrix(4, 2, [(0, 1), (1, 2), (2, 3), (0, 3)])
 62    J, R = port_hamiltonian_matrices(F, 0.2, 0.2)
 63    x0 = rng.normal(size=J.shape[0])
 64    stated_bound = 2.0 / np.linalg.eigvalsh(R).max()
 65    # Actual linear explicit-Euler condition is rho(I + dt(J-R)) <= 1.
 66    A = J - R
 67    eig = np.linalg.eigvals(A)
 68    actual_limit = np.inf
 69    for lam in eig:
 70        if lam.real < 0:
 71            candidate = -2.0 * lam.real / (abs(lam) ** 2)
 72            actual_limit = min(actual_limit, candidate)
 73    dts = [0.001, 0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1.0, 10.0]
 74    records = []
 75    for dt in dts:
 76        h, _ = simulate_ph(J, R, x0, dt, 100)
 77        records.append({
 78            "dt": dt,
 79            "spectral_radius": float(max(abs(1.0 + dt * eig))),
 80            "final_energy": float(h[-1]),
 81            "max_energy": float(h.max()),
 82        })
 83    return {
 84        "stated_damping_only_bound": float(stated_bound),
 85        "actual_linear_euler_limit": float(actual_limit),
 86        "records": records,
 87    }
 88
 89
 90def run_depth_experiment(rng):
 91    n, d = 12, 2
 92    edges = [(i, i + 1) for i in range(n - 1)] + [(0, n - 1)]
 93    F = formation_matrix(n, d, edges)
 94    J, R = port_hamiltonian_matrices(F, damping_node=0.12, damping_edge=0.12)
 95    x0 = rng.normal(size=J.shape[0])
 96    node0 = x0[: n * d]
 97    L = F.T @ F
 98    A = np.eye(n * d) - L / max(np.linalg.eigvalsh(L).max(), 1e-12)
 99    unconstrained = node0.copy()
100    ph = x0.copy()
101    rows = []
102    for k in range(51):
103        rows.append({
104            "step": k,
105            "unconstrained_node_norm": float(np.linalg.norm(unconstrained)),
106            "ph_total_energy": float(energy(ph, np.ones_like(ph))),
107            "ph_node_norm": float(np.linalg.norm(ph[: n * d])),
108        })
109        unconstrained = 1.18 * (A @ unconstrained)
110        ph = ph + 0.08 * ((J - R) @ ph)
111    return rows
112
113
114def main():
115    rng = np.random.default_rng(3165)
116    math = verify_math(rng)
117    stability = stability_sweep(rng)
118    depth = run_depth_experiment(rng)
119    out = {
120        "math": math,
121        "stability": stability,
122        "depth_summary": {
123            "unconstrained_norm_step0": depth[0]["unconstrained_node_norm"],
124            "unconstrained_norm_step50": depth[-1]["unconstrained_node_norm"],
125            "ph_energy_step0": depth[0]["ph_total_energy"],
126            "ph_energy_step50": depth[-1]["ph_total_energy"],
127            "ph_max_energy": max(x["ph_total_energy"] for x in depth),
128            "ph_node_norm_step0": depth[0]["ph_node_norm"],
129            "ph_node_norm_step50": depth[-1]["ph_node_norm"],
130        },
131    }
132    with open("results.json", "w") as f:
133        json.dump(out, f, indent=2)
134    print(json.dumps(out, indent=2))
135
136
137if __name__ == "__main__":
138    main()