1import json
  2import numpy as np
  3from pathlib import Path
  4
  5
  6def stats(k, p):
  7    k = np.asarray(k, float); p = np.asarray(p, float)
  8    mean_k = p @ k
  9    return mean_k, p @ (k*k), (p @ (k*k))/mean_k
 10
 11
 12def order_parameter(m, p):
 13    mbar = p @ m
 14    return float(p @ ((m - mbar)**2)), float(mbar)
 15
 16
 17def jacobian(m, k, p, gain, reg=0.0, mode="degree"):
 18    # gain = beta*J; this is the exact Jacobian of the supplied ODE.
 19    mean_k = p @ k
 20    z = gain * k * ((p * k) @ m) / mean_k
 21    d = 1.0 - np.tanh(z)**2
 22    A = -np.eye(len(k)) + gain * d[:, None] * (k[:, None] * (p*k)[None, :] / mean_k)
 23    if reg:
 24        if mode == "degree":
 25            A -= 2.0 * reg * (np.eye(len(k)) - np.ones((len(k), 1)) @ p[None, :])
 26        elif mode == "global":
 27            A -= 2.0 * reg * np.eye(len(k))
 28    return A
 29
 30
 31def integrate(k, p, gain, reg=0.0, mode="none", seed=0, steps=6000, dt=.05):
 32    rng = np.random.default_rng(seed)
 33    # Small degree-correlated perturbation selects one of the pitchfork branches.
 34    m = 1e-3 * (k - p @ k) + 1e-5 * rng.normal(size=len(k))
 35    mean_k = p @ k
 36    for _ in range(steps):
 37        u = (p*k) @ m / mean_k
 38        dm = -m + np.tanh(gain*k*u)
 39        if reg == 0 or mode == "none":
 40            pass
 41        elif mode == "degree":
 42            dm -= 2.0*reg*(m - p @ m)
 43        elif mode == "global":
 44            dm -= 2.0*reg*m
 45        m += dt*dm
 46        if not np.all(np.isfinite(m)):
 47            raise FloatingPointError("nonfinite state")
 48    S, mb = order_parameter(m, p)
 49    return m, S, mb
 50
 51
 52def restricted_eigenvalue(A, p):
 53    # Diagnostic specified in the idea: project onto p-weighted zero-mean vectors.
 54    n = len(p)
 55    P = np.eye(n) - np.ones((n, 1)) @ p[None, :]
 56    # Nullspace basis of p^T using SVD; projected operator is represented there.
 57    _, _, vh = np.linalg.svd(p[None, :])
 58    Q = vh[1:].T
 59    vals = np.linalg.eigvals(Q.T @ P @ A @ Q)
 60    return float(np.max(np.real(vals)))
 61
 62
 63def run():
 64    rows = []
 65    # Equal mean degree, increasing heterogeneity: p=(1/2,1/2), k=(3-d,3+d).
 66    gains = np.linspace(.08, .70, 32)
 67    for d in [0.0, 1.0, 2.0, 4.0]:
 68        k = np.array([3.0-d, 3.0+d]); p = np.array([.5, .5])
 69        mk, mk2, ratio = stats(k,p)
 70        predicted = 1.0/ratio
 71        for g in gains:
 72            m,S,mb = integrate(k,p,g,seed=123)
 73            A = jacobian(np.zeros(2),k,p,g)
 74            rows.append(dict(kind="sweep", heterogeneity=d, gain=float(g), S=S,
 75                             mean=float(mb), lambda_full=float(np.max(np.linalg.eigvals(A).real)),
 76                             lambda_sep=restricted_eigenvalue(A,p), predicted_gain=predicted,
 77                             k2_over_k=ratio))
 78    # Regularizer comparison at a supercritical gain.
 79    k=np.array([1.,5.]); p=np.array([.5,.5]); g=.55
 80    for mode, reg in [("none",0.),("global",.20),("degree",.20),("degree",.50)]:
 81        m,S,mb=integrate(k,p,g,reg=reg,mode=mode,seed=123)
 82        rows.append(dict(kind="regularizer", mode=mode, reg=reg, gain=g, S=S, mean=mb,
 83                         norm=float(np.sqrt(p@(m*m)))))
 84    out=Path("results.json")
 85    out.write_text(json.dumps(rows, indent=2))
 86    # Compact summary used by the report and easy reproduction.
 87    summary=[]
 88    for d in [0.,1.,2.,4.]:
 89        sub=[r for r in rows if r.get("kind")=="sweep" and r["heterogeneity"]==d]
 90        pred=sub[0]["predicted_gain"]
 91        # first sampled point whose S exceeds a robust numerical threshold
 92        active=[r["gain"] for r in sub if r["S"]>1e-8]
 93        summary.append({"heterogeneity":d,"predicted":pred,"observed_grid_onset":min(active) if active else None,
 94                        "max_S":max(r["S"] for r in sub)})
 95    print(json.dumps({"threshold_summary":summary,
 96                      "regularizers":[r for r in rows if r.get("kind")=="regularizer"]}, indent=2))
 97
 98
 99if __name__ == "__main__":
100    run()