1import json
  2from pathlib import Path
  3import numpy as np
  4
  5SEED = 1168
  6rng = np.random.default_rng(SEED)
  7
  8
  9def phi(u, a, b, c, order="12"):
 10    """Path energy for F1=-a*u1+b*u2 and F2=-a*u2+c*u1."""
 11    u = np.asarray(u)
 12    cross = c if order == "12" else b
 13    return -0.5 * a * np.sum(u * u, axis=-1) + cross * u[..., 0] * u[..., 1]
 14
 15
 16def field(u, a, b, c):
 17    u = np.asarray(u)
 18    return np.stack((-a*u[..., 0] + b*u[..., 1],
 19                     -a*u[..., 1] + c*u[..., 0]), axis=-1)
 20
 21
 22def grid_distribution(a, b, c, gamma, order, n=501, lo=-1.0, hi=1.0):
 23    x = np.linspace(lo, hi, n)
 24    X, Y = np.meshgrid(x, x, indexing="ij")
 25    U = np.stack((X, Y), axis=-1)
 26    logits = phi(U, a, b, c, order) / gamma
 27    logits -= logits.max()
 28    weights = np.exp(logits)
 29    weights /= weights.sum()
 30    return x, X, Y, weights
 31
 32
 33def moments(a, b, c, gamma, order, n=501):
 34    _, X, Y, W = grid_distribution(a, b, c, gamma, order, n=n)
 35    return {
 36        "cross": float(np.sum(X * Y * W)),
 37        "var1": float(np.sum(X * X * W)),
 38        "var2": float(np.sum(Y * Y * W)),
 39    }
 40
 41
 42def sample_joint(a, b, c, gamma, order, count=50000, n=301):
 43    _, X, Y, W = grid_distribution(a, b, c, gamma, order, n=n)
 44    indices = rng.choice(W.size, size=count, p=W.ravel())
 45    return np.stack((X.ravel()[indices], Y.ravel()[indices]), axis=1)
 46
 47
 48def sample_independent(a, gamma, count=50000, n=10001):
 49    # Independent Gibbs baseline from each self-quadratic q_i=-a*u_i^2/2.
 50    x = np.linspace(-1., 1., n)
 51    p = np.exp(-0.5*a*x*x/gamma)
 52    p /= p.sum()
 53    i = rng.choice(n, size=(count, 2), p=p)
 54    return x[i]
 55
 56
 57def main():
 58    a = 2.0
 59    records = {}
 60
 61    # Prediction 1: Delta Phi=(c-b)u1u2, whose uniform-box RMS is |c-b|/3.
 62    disagreement = []
 63    grid = np.linspace(-1., 1., 401)
 64    X, Y = np.meshgrid(grid, grid, indexing="ij")
 65    U = np.stack((X, Y), axis=-1)
 66    for d in [0., .2, .5, 1.0, 1.5]:
 67        b, c = .4 - d/2, .4 + d/2
 68        observed = float(np.sqrt(np.mean((phi(U,a,b,c,"12") - phi(U,a,b,c,"21"))**2)))
 69        predicted = abs(c-b)/3.
 70        disagreement.append({"asymmetry_abs_c_minus_b": d,
 71                             "predicted_rms": predicted,
 72                             "observed_rms": observed})
 73    records["path_disagreement_scaling"] = disagreement
 74
 75    # Prediction 2: if b=c, the field is conservative and the two paths coincide.
 76    conservative = []
 77    for k in [-1., -.2, 0., .4, 1.0]:
 78        d = float(np.max(np.abs(phi(U,a,k,k,"12") - phi(U,a,k,k,"21"))))
 79        conservative.append({"b_equals_c": k, "max_path_difference": d})
 80    records["conservative_field_check"] = conservative
 81
 82    # Prediction 3: on an unbounded quadratic Gibbs model, Cov(u1,u2)=gamma*k/(a^2-k^2).
 83    # Compare this scaling against exact bounded-box quadrature for a stable k.
 84    b = 0.0
 85    c = 0.5
 86    temperature = []
 87    for gamma in [.20, .35, .50, .75, 1.0, 1.5]:
 88        observed = moments(a,b,c,gamma,"12",n=501)["cross"]
 89        predicted = gamma*c/(a*a-c*c)
 90        temperature.append({"gamma": gamma, "predicted_unbounded_cov": predicted,
 91                           "observed_bounded_cov": observed,
 92                           "ratio_observed_to_predicted": observed/predicted})
 93    records["temperature_covariance_scaling"] = temperature
 94
 95    # Direct numerical gradient verification of the path construction.
 96    u = np.array([.31, -.47])
 97    eps = 1e-5
 98    grad_num = []
 99    for j in range(2):
100        plus, minus = u.copy(), u.copy()
101        plus[j] += eps; minus[j] -= eps
102        grad_num.append((phi(plus,a,.1,.7,"12") - phi(minus,a,.1,.7,"12"))/(2*eps))
103    expected = np.array([-a*u[0] + .7*u[1], -a*u[1] + .7*u[0]])
104    records["gradient_check"] = {"numerical": grad_num, "expected_path_gradient": expected.tolist(),
105                                  "max_abs_error": float(np.max(np.abs(np.array(grad_num)-expected)))}
106
107    # Secondary coordination comparison: coherent path Gibbs versus independent Gibbs.
108    joint = sample_joint(a, .0, .8, .5, "12")
109    independent = sample_independent(a, .5, count=len(joint))
110    records["coordination_comparison"] = {
111        "path_joint_mean_product": float(np.mean(joint[:,0]*joint[:,1])),
112        "independent_mean_product": float(np.mean(independent[:,0]*independent[:,1])),
113        "path_joint_abs_product": float(np.mean(np.abs(joint[:,0]*joint[:,1]))),
114        "independent_abs_product": float(np.mean(np.abs(independent[:,0]*independent[:,1])))
115    }
116
117    out = Path("results.json")
118    out.write_text(json.dumps(records, indent=2))
119    print(json.dumps(records, indent=2))
120
121
122if __name__ == "__main__":
123    main()