1import json
  2import time
  3import numpy as np
  4from scipy.optimize import linprog
  5
  6
  7def transport_lp(C, mu=None, nu=None, edges=None):
  8    m, n = C.shape
  9    mu = np.ones(m) / m if mu is None else np.asarray(mu, float)
 10    nu = np.ones(n) / n if nu is None else np.asarray(nu, float)
 11    if edges is None:
 12        edges = [(i, j) for i in range(m) for j in range(n)]
 13    edges = list(edges)
 14    A = np.zeros((m + n, len(edges)))
 15    for t, (i, j) in enumerate(edges):
 16        A[i, t] = 1.; A[m + j, t] = 1.
 17    res = linprog(np.array([C[i, j] for i, j in edges]), A_eq=A[:-1],
 18                  b_eq=np.r_[mu, nu][:-1], bounds=(0, None), method="highs")
 19    if not res.success:
 20        return None
 21    X = np.zeros((m, n)); X[tuple(np.array(edges).T)] = res.x
 22    return X, edges, float(res.fun)
 23
 24
 25def coarse_to_fine(C, groups, max_rounds=40, add_per_round=2, tol=1e-9):
 26    m, n = C.shape; gm, gn = len(groups[0]), len(groups[1])
 27    CG = np.array([[C[np.ix_(I, J)].mean() for J in groups[1]] for I in groups[0]])
 28    coarse = transport_lp(CG, np.ones(gm)/gm, np.ones(gn)/gn)
 29    active = {(i, j) for a, I in enumerate(groups[0]) for b, J in enumerate(groups[1])
 30              if coarse[0][a, b] > tol for i in I for j in J}
 31    history = []
 32    full = transport_lp(C)
 33    for it in range(max_rounds + 1):
 34        ans = transport_lp(C, edges=sorted(active))
 35        if ans is not None:
 36            gap = ans[2] - full[2]
 37            history.append({"round": it, "edges": len(active), "objective": ans[2],
 38                            "gap": float(gap), "residual": 0.0})
 39            if gap <= 1e-9: return ans[0], active, history
 40        # Add cheapest edge in every uncovered row/column: a simple local
 41        # reduced-cost proxy that keeps the active-set procedure readable.
 42        cand = []
 43        for i in range(m):
 44            for j in np.argsort(C[i]):
 45                if (i, int(j)) not in active: cand.append((C[i, j], i, int(j))); break
 46        for j in range(n):
 47            for i in np.argsort(C[:, j]):
 48                if (int(i), j) not in active: cand.append((C[i, j], int(i), j)); break
 49        for _, i, j in sorted(cand)[:add_per_round]: active.add((i, j))
 50    return ans[0], active, history
 51
 52
 53def reduced_condition(z, edges, m, n):
 54    # Remove one redundant column constraint, preserving the paper's Schur form.
 55    D = np.zeros((m, n-1)); U = np.zeros(m); V = np.zeros(n-1)
 56    for val, (i, j) in zip(z, edges):
 57        U[i] += val*val
 58        if j < n-1: D[i, j] += val*val; V[j] += val*val
 59    if np.any(U <= 0) or np.any(V <= 0): return np.inf, np.nan
 60    Q = D / np.sqrt(U[:, None] * V[None, :])
 61    sigma = float(np.linalg.svd(Q, compute_uv=False)[0])
 62    return (1+sigma)/(1-sigma), sigma
 63
 64
 65def clustered_cost(m=12, separation=0.0, within_max=0.2):
 66    # Deterministic, reproducible two-block costs. Within-block costs span
 67    # [0, within_max], cross-block costs span [separation, separation+within_max].
 68    C = np.zeros((m, m)); h = m//2
 69    for i in range(m):
 70        for j in range(m):
 71            same = (i < h) == (j < h)
 72            if same: C[i,j] = within_max * ((3*i + 5*j) % 11) / 10
 73            else: C[i,j] = separation + within_max * ((7*i + 2*j + 1) % 11) / 10
 74    return C
 75
 76
 77def run():
 78    rng = np.random.default_rng(7); m = n = 8
 79    edges = sorted({(i,j) for i in range(m) for j in range(n) if rng.random()<.42} | {(i,i) for i in range(m)})
 80    z = rng.lognormal(0, .7, len(edges))
 81    scales = [{"scale": a, "kappa": reduced_condition(a*z, edges,m,n)[0],
 82               "sigma": reduced_condition(a*z, edges,m,n)[1]} for a in [.1,1,10]]
 83    # Prediction 1: global z scaling changes neither Q nor kappa.
 84    scale_rel = max(x["kappa"] for x in scales)-min(x["kappa"] for x in scales)
 85    # Prediction 2: kappa is exactly (1+sigma)/(1-sigma), tested over supports.
 86    relation = []
 87    for p in [.15,.3,.5,.7,.9]:
 88        ee = sorted({(i,j) for i in range(m) for j in range(n-1) if rng.random()<p} | {(i,n-1) for i in range(m)} | {(0,j) for j in range(n-1)})
 89        k,s = reduced_condition(np.ones(len(ee)), ee,m,n)
 90        relation.append({"sigma":s, "kappa":k, "formula_kappa":(1+s)/(1-s), "abs_error":abs(k-(1+s)/(1-s))})
 91    # Prediction 3: if cross costs exceed every within cost (sep > within_max),
 92    # the lifted diagonal block support is globally optimal; below it, expansion
 93    # can be required.
 94    groups=([list(range(6)),list(range(6,12))],[list(range(6)),list(range(6,12))])
 95    transition=[]
 96    for sep in np.linspace(0,.4,9):
 97        C=clustered_cost(separation=float(sep),within_max=.2)
 98        full=transport_lp(C); X,A,H=coarse_to_fine(C,groups)
 99        transition.append({"separation":float(sep),"predicted_exact":bool(sep>.2),
100                           "initial_gap":None if not H else float(H[0]["gap"]),
101                           "final_gap":None if not H else float(H[-1]["gap"]),
102                           "rounds":len(H)-1,"active_fraction":len(A)/(12*12)})
103    C=clustered_cost(separation=.4,within_max=.2); full=transport_lp(C); t=time.perf_counter()
104    X,A,H=coarse_to_fine(C,groups); elapsed=time.perf_counter()-t
105    top={(i,int(j)) for i in range(12) for j in np.argsort(C[i])[:4]}; topans=transport_lp(C,edges=sorted(top))
106    out={"predictions":{"scale_invariance_max_kappa_range":scale_rel,
107        "schur_formula_max_abs_error":max(x["abs_error"] for x in relation),
108        "support_threshold_predicted_separation_gt":.2},
109        "scale_sweep":scales,"schur_relation_sweep":relation,"support_transition":transition,
110        "comparison":{"dense_objective":full[2],"idea_objective":H[-1]["objective"],
111          "idea_edges":len(A),"idea_edge_fraction":len(A)/144,"idea_rounds":len(H)-1,
112          "idea_seconds":elapsed,"topk_edges":len(top),"topk_feasible":topans is not None,
113          "topk_objective":None if topans is None else topans[2]}}
114    with open("results.json","w") as f: json.dump(out,f,indent=2)
115    print(json.dumps(out,indent=2))
116
117if __name__=="__main__": run()