problems.py

"""
problems.py -- Max-Cut / Ising instances, certified references, and fair baselines.

Everything the reproduction scripts need to build instances, get an honest reference, and
compare the runtime against classical methods under a wall-clock deadline.
"""
import math, random

# Wall-clock cost model (ms). The premise being reproduced: a configuration measurement is
# ONE parallel acquisition, constant-cost in n (the runtime's round cost), while digital
# sequential search pays real per-flip cost that scales with node degree. On a CPU both run
# on your processor; we charge each its modelled cost to expose the crossover the hardware
# implies. These are the SAME constants used across the study.
T_MEAS = 0.05    # ms per parallel configuration measurement (constant in n)
T_OP = 0.001     # ms per digital scalar op


# ----------------------------- instance generators ----------------------------- #
def graph_random(n, seed, density=0.5):
    rng = random.Random(seed); E = {}
    for i in range(n):
        for j in range(i + 1, n):
            if rng.random() < density:
                E[(i, j)] = 0.1 + 0.9 * rng.random()
    return E


def graph_spinglass(n, seed, density=0.12):
    """Frustrated signed-weight graph (Ising spin glass) -- the genuinely hard case."""
    rng = random.Random(seed); E = {}
    for i in range(n):
        for j in range(i + 1, n):
            if rng.random() < density:
                E[(i, j)] = 2 * rng.random() - 1
    return E


def n_nodes(E):
    m = 0
    for (i, j) in E:
        m = max(m, i, j)
    return m + 1


def adjacency(E, n):
    adj = {i: [] for i in range(n)}
    for (i, j), w in E.items():
        adj[i].append((j, w)); adj[j].append((i, w))
    return adj


# ----------------------------- cut value ----------------------------- #
def cut_value(assignment, E):
    c = 0.0
    for (i, j), w in E.items():
        if assignment[i] != assignment[j]:
            c += w
    return c


def cut_delta(assignment, adj, i):
    d = 0.0
    for j, w in adj[i]:
        d += (-w) if assignment[i] != assignment[j] else w
    return d


# ----------------------------- references ----------------------------- #
def offline_reference(E, n, seed, restarts=6, sweeps=60):
    """Strong anytime-optimal reference: multi-restart local search + SA polish, NO clock.
    The practical 'achievable optimum' used to normalize solution quality."""
    adj = adjacency(E, n); best = -1e18
    for r in range(restarts):
        rng = random.Random(seed * 131 + r)
        a = [1 if rng.random() < 0.5 else 0 for _ in range(n)]
        c = cut_value(a, E)
        improved = True
        while improved:
            improved = False
            for i in range(n):
                d = cut_delta(a, adj, i)
                if d > 1e-12:
                    a[i] ^= 1; c += d; improved = True
        for s in range(sweeps):
            T = max(0.01, 1.0 * (1 - s / sweeps))
            for i in range(n):
                d = cut_delta(a, adj, i)
                if d > 0 or rng.random() < math.exp(d / max(T, 1e-6)):
                    a[i] ^= 1; c += d
        improved = True
        while improved:
            improved = False
            for i in range(n):
                d = cut_delta(a, adj, i)
                if d > 1e-12:
                    a[i] ^= 1; c += d; improved = True
        best = max(best, c)
    return best


def random_stats(E, n, seed, tries=64):
    """Mean and std of random-cut value -- the 'did no useful work' baseline."""
    rng = random.Random(seed * 999 + 1); vals = []
    for _ in range(tries):
        a = [1 if rng.random() < 0.5 else 0 for _ in range(n)]
        vals.append(cut_value(a, E))
    mu = sum(vals) / len(vals)
    var = sum((v - mu) ** 2 for v in vals) / len(vals)
    return mu, math.sqrt(var)


# ----------------------------- baselines under a deadline ----------------------------- #
def blind_parallel(E, n, deadline_ms, seed, batch=8):
    """Blind parallel sampling: same measurement cost as the runtime, no steering."""
    rng = random.Random(seed * 7 + 1); best = -1e18; t = 0.0
    while t < deadline_ms:
        for _ in range(batch):
            a = [1 if rng.random() < 0.5 else 0 for _ in range(n)]
            best = max(best, cut_value(a, E))
        t += batch * T_MEAS
    return best


def digital_timed(E, adj, n, deadline_ms, seed, kind="greedy"):
    """Digital sequential search charged O(degree) per flip. Returns (best_cut, flips)."""
    rng = random.Random(seed * 13 + 5)
    a = [1 if rng.random() < 0.5 else 0 for _ in range(n)]
    c = cut_value(a, E); best = c; t = 0.0; flips = 0
    tabu = {}; it = 0
    while t < deadline_ms:
        it += 1
        i = rng.randrange(n); deg = len(adj[i])
        d = cut_delta(a, adj, i)
        t += (1 + deg) * T_OP; flips += 1
        if kind == "greedy":
            acc = d > 0
        elif kind == "sa":
            acc = d > 0 or rng.random() < math.exp(min(20, d))
        elif kind == "tabu":
            acc = (d > 0) or (tabu.get(i, 0) <= it and d > -abs(c) * 0.001)
            if acc:
                tabu[i] = it + max(3, n // 20)
        else:
            acc = d > 0
        if acc:
            a[i] ^= 1; c += d; best = max(best, c)
    return best, flips