1import json
  2import math
  3import random
  4from itertools import combinations
  5import numpy as np
  6from scipy.optimize import minimize_scalar
  7
  8
  9def mu(G, q):
 10    q = np.asarray(q, dtype=float)
 11    d = np.exp(q - np.max(q))
 12    return float(np.max((G @ d) / d))
 13
 14
 15def rho(G):
 16    return float(np.max(np.abs(np.linalg.eigvals(G))))
 17
 18
 19def optimize_q(G):
 20    # For this sanity check, optimize the exact nonsmooth certificate in the
 21    # only meaningful relative coordinate for two nodes.
 22    if G.shape == (2, 2):
 23        f = lambda x: mu(G, [x, -x])
 24        r = minimize_scalar(f, bounds=(-10, 10), method="bounded", options={"xatol": 1e-12})
 25        return np.array([r.x, -r.x]), float(r.fun)
 26    q = np.zeros(G.shape[0])
 27    return q, mu(G, q)
 28
 29
 30def dominant_cycle_matrix(c, n=6, reverse_fraction=0.2):
 31    G = np.zeros((n, n), dtype=float)
 32    # Dominant directed 3-cycle has product c^3. Reciprocal weak links make
 33    # the proposed symmetric pair score able to discover its node pairs.
 34    G[1, 0] = c * 4.0
 35    G[2, 1] = c
 36    G[0, 2] = c / 4.0
 37    G[0, 1] = reverse_fraction * c * 4.0
 38    G[1, 2] = reverse_fraction * c
 39    G[2, 0] = reverse_fraction * c / 4.0
 40    for i in range(3, n):
 41        G[i, i - 1] = 0.08
 42        G[i - 1, i] = 0.05
 43    return G
 44
 45
 46def edge_score(G, i, j, delta=1e-8):
 47    return math.log(G[i, j] + delta) + math.log(G[j, i] + delta)
 48
 49
 50def greedy_partition(G, max_size=3):
 51    clusters = [{i} for i in range(G.shape[0])]
 52    while True:
 53        candidates = []
 54        for a, b in combinations(range(len(clusters)), 2):
 55            if len(clusters[a] | clusters[b]) <= max_size:
 56                score = max(edge_score(G, i, j) for i in clusters[a] for j in clusters[b])
 57                candidates.append((score, a, b))
 58        if not candidates:
 59            break
 60        score, a, b = max(candidates)
 61        if score < math.log(0.08 * 0.05):
 62            break
 63        merged = clusters[a] | clusters[b]
 64        clusters = [x for k, x in enumerate(clusters) if k not in (a, b)] + [merged]
 65    return sorted([sorted(x) for x in clusters], key=lambda x: x[0])
 66
 67
 68def cycle_internal(P, cycle=(0, 1, 2)):
 69    owner = {i: k for k, S in enumerate(P) for i in S}
 70    return len({owner[i] for i in cycle}) == 1
 71
 72
 73def visible_cycle_gain(G, P, cycle=(0, 1, 2)):
 74    return 0.0 if cycle_internal(P, cycle) else float(np.prod([G[cycle[(k + 1) % 3], cycle[k]] for k in range(3)]))
 75
 76
 77def main():
 78    np.random.seed(7)
 79    random.seed(7)
 80    out = {"normalization_sweep": [], "imbalance_sweep": [], "cluster_sweep": []}
 81
 82    # Prediction 1: diagonal scaling cannot beat the cycle boundary c=1.
 83    for c in [0.25, 0.5, 0.9, 1.0, 1.1, 1.5, 2.0]:
 84        G = np.array([[0., 4*c], [c/4, 0.]])
 85        q, best = optimize_q(G)
 86        out["normalization_sweep"].append({"c": c, "rho": rho(G), "mu_q0": mu(G, [0, 0]),
 87            "mu_optimized": best, "q": q.tolist(), "predicted_optimum": c,
 88            "stable_by_certificate": best < 1.0})
 89
 90    # Prediction 2: unscaled certificate depends on imbalance, optimized one does not.
 91    for a in [1., 2., 4., 8.]:
 92        c = 0.7
 93        G = np.array([[0., a*c], [c/a, 0.]])
 94        _, best = optimize_q(G)
 95        out["imbalance_sweep"].append({"imbalance": a, "mu_unscaled": mu(G, [0, 0]),
 96            "mu_optimized": best, "predicted_optimized": c})
 97
 98    # Prediction 3: greedy gain-weighted merging internalizes the cycle when
 99    # cluster capacity is 3; the fixed random partition leaves it exposed.
100    singleton_P = [[0], [1], [2], [3], [4], [5]]
101    random_P = [[0, 2], [1], [3, 4, 5]]
102    for c in [0.5, 0.9, 1.0, 1.1, 1.5]:
103        G = dominant_cycle_matrix(c)
104        greedy = greedy_partition(G, max_size=3)
105        out["cluster_sweep"].append({"c": c, "rho": rho(G), "greedy_partition": greedy,
106            "singleton_internal": cycle_internal(singleton_P),
107            "singleton_boundary_cycle_gain": visible_cycle_gain(G, singleton_P),
108            "greedy_internal": cycle_internal(greedy),
109            "greedy_boundary_cycle_gain": visible_cycle_gain(G, greedy),
110            "random_internal": cycle_internal(random_P),
111            "random_boundary_cycle_gain": visible_cycle_gain(G, random_P),
112            "cycle_product_predicted": c**3})
113
114    norm_err = max(abs(x["mu_optimized"] - x["predicted_optimum"]) for x in out["normalization_sweep"])
115    imb_err = max(abs(x["mu_optimized"] - x["predicted_optimized"]) for x in out["imbalance_sweep"])
116    boundary_ok = all((x["c"] < 1) == x["stable_by_certificate"] for x in out["normalization_sweep"] if x["c"] != 1)
117    boundary_at_one = min(abs(x["mu_optimized"] - 1) for x in out["normalization_sweep"] if x["c"] == 1) < 1e-8
118    transition_ok = all(x["greedy_internal"] and not x["singleton_internal"] and not x["random_internal"] and
119        x["greedy_boundary_cycle_gain"] == 0.0 and
120        abs(x["singleton_boundary_cycle_gain"] - x["cycle_product_predicted"]) < 1e-10 and
121        abs(x["random_boundary_cycle_gain"] - x["cycle_product_predicted"]) < 1e-10
122        for x in out["cluster_sweep"])
123    out["predictions"] = {
124        "boundary_predicted_c": 1.0, "boundary_confirmed": boundary_ok and boundary_at_one,
125        "normalization_max_abs_error": norm_err,
126        "imbalance_max_abs_error": imb_err,
127        "cluster_transition_confirmed": transition_ok,
128        "mechanism_scope": "clustering removes the dominant channel from the boundary graph; it does not reduce full spectral radius"
129    }
130    print(json.dumps(out, indent=2, sort_keys=True))
131
132if __name__ == "__main__":
133    main()