1import json
  2import time
  3import numpy as np
  4
  5
  6def energy(S, mu, v):
  7    return float(-0.5 * v @ S @ v - mu @ v)
  8
  9
 10def reduce_blocks(S, mu, v_H, H, Q):
 11    """Exact conditional objective represented by the active block and field."""
 12    H = np.asarray(H, dtype=int); Q = np.asarray(Q, dtype=int)
 13    return S[np.ix_(Q, Q)], mu[Q] + S[np.ix_(Q, H)] @ v_H
 14
 15
 16def certified_candidates(S, mu, v, H, Q, eps=1e-12):
 17    """Return active indices certified for every assignment of the other active spins."""
 18    H = np.asarray(H, dtype=int); Q = np.asarray(Q, dtype=int)
 19    out = []
 20    for i in Q:
 21        fixed = mu[i] + (S[i, H] @ v[H] if len(H) else 0.0)
 22        unresolved_bound = np.abs(S[i, Q]).sum() - abs(S[i, i])
 23        if abs(fixed) > unresolved_bound + eps and np.sign(v[i]) == np.sign(fixed):
 24            out.append(int(i))
 25    return out
 26
 27
 28def verify_math(seed=7):
 29    rng = np.random.default_rng(seed)
 30    n = 7
 31    A = rng.normal(size=(n, n)); S = (A + A.T) / 2; np.fill_diagonal(S, 0)
 32    mu = rng.normal(size=n); v = rng.choice([-1, 1], size=n)
 33    # Exact conditional identity for several arbitrary H/Q assignments.
 34    H = np.array([0, 3, 5]); Q = np.array([1, 2, 4, 6])
 35    block, meff = reduce_blocks(S, mu, v[H], H, Q)
 36    C = energy(S, mu, v) - (-0.5 * v[Q] @ block @ v[Q] - meff @ v[Q])
 37    identity_max_err = 0.0
 38    for bits in range(1 << len(Q)):
 39        q = np.array([1 if bits >> k & 1 else -1 for k in range(len(Q))])
 40        full = np.zeros(n, dtype=int); full[H] = v[H]; full[Q] = q
 41        reduced = C - 0.5 * q @ block @ q - meff @ q
 42        identity_max_err = max(identity_max_err, abs(energy(S, mu, full) - reduced))
 43    # Construct a guaranteed certified case and exhaustively check all Q assignments.
 44    m = 6; B = rng.normal(size=(m, m)); T = (B + B.T) / 2; np.fill_diagonal(T, 0)
 45    b = np.zeros(m); b[0] = 5.0; b[1] = -4.5
 46    state = np.ones(m, dtype=int); state[1] = -1
 47    H2 = np.array([], dtype=int); Q2 = np.arange(m)
 48    cand = certified_candidates(T, b, state, H2, Q2)
 49    violations = 0
 50    for i in cand:
 51        for bits in range(1 << len(Q2)):
 52            q = np.array([1 if bits >> k & 1 else -1 for k in range(len(Q2))])
 53            if state[i] != np.sign(b[i] + T[i] @ q):
 54                violations += 1
 55    return {"conditional_identity_max_abs_error": identity_max_err,
 56            "certified_indices": cand, "certification_violations": violations}
 57
 58
 59def run_solver(S, mu, v0, steps, reduce=False, eps=1e-10):
 60    n = len(mu); v = v0.copy(); frozen = np.zeros(n, dtype=bool)
 61    active = np.arange(n); active_bias = mu.copy(); work = 0
 62    cert_viol = 0; active_trace = []
 63    t0 = time.perf_counter()
 64    for _ in range(steps):
 65        if reduce:
 66            # Certification uses the exact current frozen field and remaining tail.
 67            H = np.flatnonzero(frozen); Q = np.flatnonzero(~frozen)
 68            cand = certified_candidates(S, mu, v, H, Q, eps)
 69            for i in cand:
 70                frozen[i] = True
 71            active = np.flatnonzero(~frozen)
 72            if len(active):
 73                # This is the same induced field as incremental maintenance.
 74                active_bias = mu[active] + S[np.ix_(active, np.flatnonzero(frozen))] @ v[frozen]
 75        if len(active):
 76            fields = active_bias + S[np.ix_(active, active)] @ v[active]
 77            v[active] = np.where(fields >= 0, 1, -1)
 78            work += len(active) * len(active)
 79        active_trace.append(int(len(active)))
 80        # Count post-update violations for frozen spins against the actual full state.
 81        if reduce and frozen.any():
 82            h = mu + S @ v
 83            cert_viol += int(np.sum(frozen & (v * h <= 0)))
 84    return {"energy": energy(S, mu, v), "state": v, "seconds": time.perf_counter()-t0,
 85            "work_matvec_entries": int(work), "active_trace": active_trace,
 86            "frozen": int(frozen.sum()), "certification_violations": cert_viol}
 87
 88
 89def benchmark(seed=11, n=1200, steps=35):
 90    rng = np.random.default_rng(seed)
 91    A = rng.normal(size=(n, n)).astype(np.float64)
 92    S = (A + A.T) * (0.5 / np.sqrt(n)); np.fill_diagonal(S, 0)
 93    # Strong polarized coordinates make the claimed late-stage tail reduction visible,
 94    # while the remaining coordinates retain a nontrivial dense Ising problem.
 95    mu = rng.normal(0, 0.10, n); strong = int(0.70*n)
 96    mu[:strong] += 30.0 * rng.choice([-1, 1], strong)
 97    v0 = rng.choice([-1, 1], n)
 98    base = run_solver(S, mu, v0, steps, reduce=False)
 99    idea = run_solver(S, mu, v0, steps, reduce=True)
100    return {"n": n, "steps": steps, "strong_field_fraction": strong/n,
101            "baseline": {k:v for k,v in base.items() if k != "state"},
102            "idea": {k:v for k,v in idea.items() if k != "state"},
103            "final_energy_difference_idea_minus_baseline": idea["energy"]-base["energy"],
104            "work_ratio_idea_over_baseline": idea["work_matvec_entries"]/base["work_matvec_entries"],
105            "final_state_agreement": float(np.mean(idea["state"] == base["state"]))}
106
107
108if __name__ == "__main__":
109    result = {"math": verify_math(), "benchmark": benchmark()}
110    print(json.dumps(result, indent=2))