1import json
  2import numpy as np
  3
  4
  5def make_adapter(C, b, M, tol=1e-12):
  6    n = C.shape[1]
  7    Minv = np.linalg.inv(M)
  8    G = C @ Minv @ C.T
  9    up = Minv @ C.T @ np.linalg.solve(G, b)
 10    # Euclidean nullspace, then weighted orthonormalize via Gram eigendecomposition.
 11    U, s, Vh = np.linalg.svd(C, full_matrices=True)
 12    rank = int(np.sum(s > tol * (s[0] if len(s) else 1.0)))
 13    N = Vh[rank:].T
 14    gram = N.T @ M @ N
 15    ew, Q = np.linalg.eigh(gram)
 16    keep = ew > tol * max(1.0, ew.max())
 17    V = N @ Q[:, keep] @ np.diag(1.0 / np.sqrt(ew[keep]))
 18    return up, V, rank, s
 19
 20
 21def mass_norm(x, M):
 22    return float(np.sqrt(max(0.0, x @ M @ x)))
 23
 24
 25def run():
 26    rng = np.random.default_rng(7)
 27    n, m = 12, 3
 28    # A nontrivial SPD mass matrix and full-row-rank sampled trace operator.
 29    A = rng.normal(size=(n, n))
 30    M = A.T @ A + 0.5 * np.eye(n)
 31    C = np.zeros((m, n))
 32    C[0, 0] = 1; C[0, 1] = .4
 33    C[1, 4] = 1; C[1, 7] = -.25
 34    C[2, 9] = 1; C[2, 10] = .3; C[2, 11] = -.2
 35    b = np.array([0.7, -0.2, 0.35])
 36    up, V, rank, s = make_adapter(C, b, M)
 37    # Core algebra verification.
 38    algebra = {
 39        "constraint_up": float(np.linalg.norm(C @ up - b)),
 40        "constraint_V": float(np.linalg.norm(C @ V)),
 41        "weighted_orthogonality": float(np.linalg.norm(V.T @ M @ V - np.eye(V.shape[1]))),
 42        "rank": rank, "latent_dim": int(V.shape[1]),
 43    }
 44
 45    # Frozen ambient mechanism: deliberately has normal and tangent components.
 46    # The tangent part is a skew operator, hence reduced dynamics conserve M-energy.
 47    R = rng.normal(size=(n, n)); K = R - R.T
 48    # F(u)=M^{-1}K(u-up)+ normal forcing; tangent projection removes normal forcing.
 49    Minv = np.linalg.inv(M)
 50    normal = rng.normal(size=n)
 51    def F(u):
 52        return Minv @ (K @ (u - up)) + 2.0 * normal
 53    z0 = rng.normal(size=V.shape[1])
 54    u0 = up + V @ z0
 55    dt, steps = 0.035, 600
 56
 57    def reduced_rollout():
 58        z = z0.copy(); us=[]
 59        for _ in range(steps):
 60            u = up + V @ z
 61            z = z + dt * (V.T @ M @ F(u))
 62            us.append(up + V @ z)
 63        return np.asarray(us)
 64
 65    # Same ambient Euler update, then exact affine projection in the M metric.
 66    # Projection is the standard control and uses the same F and timestep.
 67    def project(u):
 68        return u + Minv @ C.T @ np.linalg.solve(C @ Minv @ C.T, b - C @ u)
 69    def projected_rollout():
 70        u=u0.copy(); us=[]
 71        for _ in range(steps):
 72            u = project(u + dt * F(u)); us.append(u.copy())
 73        return np.asarray(us)
 74
 75    # Soft penalty control: ambient update with a finite penalty and no projection.
 76    # The penalty is applied in M^{-1} C^T coordinates.
 77    def penalty_rollout(lam=30.0):
 78        u=u0.copy(); us=[]
 79        for _ in range(steps):
 80            u = u + dt * (F(u) + lam * (Minv @ C.T @ (b - C @ u)))
 81            us.append(u.copy())
 82        return np.asarray(us)
 83
 84    red, proj, pen = reduced_rollout(), projected_rollout(), penalty_rollout()
 85    def stats(arr):
 86        residual=np.linalg.norm((arr @ C.T) - b[None,:], axis=1)
 87        e0=mass_norm(u0,M)
 88        drift=np.abs(np.array([mass_norm(x,M) for x in arr])-e0)/e0
 89        return {"max_constraint_residual":float(residual.max()),
 90                "final_constraint_residual":float(residual[-1]),
 91                "max_relative_mass_norm_drift":float(drift.max()),
 92                "final_relative_mass_norm_drift":float(drift[-1])}
 93    # With a complete weighted nullspace basis, reduced Euler equals the
 94    # corresponding M-metric projection of one ambient Euler step.
 95    u_test = u0 + dt * F(u0)
 96    z_test = z0 + dt * (V.T @ M @ F(u0))
 97    one_step_difference = np.linalg.norm((up + V @ z_test) - project(u_test))
 98    penalty_sweep = {str(lam): stats(penalty_rollout(lam))
 99                     for lam in (1.0, 10.0, 30.0, 100.0)}
100    result={"seed":7,"dt":dt,"steps":steps,"algebra":algebra,
101            "one_step_reduced_vs_projection_difference":float(one_step_difference),
102            "reduced_coordinate":stats(red),"projection_control":stats(proj),
103            "soft_penalty_lambda_30":stats(pen),"soft_penalty_sweep":penalty_sweep}
104    # A deliberately unconstrained version demonstrates why the adapter matters.
105    u=u0.copy(); raw=[]
106    for _ in range(steps):
107        u=u+dt*F(u); raw.append(u.copy())
108    result["unconstrained_ambient"] = stats(np.asarray(raw))
109    print(json.dumps(result, indent=2))
110
111if __name__ == '__main__': run()