1import json
  2import itertools
  3import numpy as np
  4
  5
  6def make_problem(n, rho, seed=0):
  7    rng = np.random.default_rng(seed)
  8    d = rng.uniform(-1.0, 1.0, n)
  9    Q = np.zeros((n, n))
 10    for i in range(n):
 11        for j in range(i + 1, n):
 12            Q[i, j] = Q[j, i] = 0.35 * rho ** (j-i)
 13    return d, Q
 14
 15
 16def objective(mask, d, Q):
 17    z = np.asarray(mask, dtype=float)
 18    return float(d @ z + 0.5 * z @ Q @ z)
 19
 20
 21def exact_optimum(d, Q, k):
 22    n = len(d); best = float('inf'); best_mask = None
 23    for comb in itertools.combinations(range(n), k):
 24        z = np.zeros(n, dtype=np.int8); z[list(comb)] = 1
 25        val = objective(z, d, Q)
 26        if val < best:
 27            best, best_mask = val, z
 28    return best, best_mask
 29
 30
 31def boundary_dp(d, Q, k, eta):
 32    # Approximate decision diagram. The state retains selected count and the
 33    # last eta decisions; interactions outside that boundary are truncated.
 34    # Each surviving state also stores a feasible backtracked mask.
 35    n = len(d); width = min(eta, n)
 36    states = {(0, 0): (0.0, np.zeros(0, dtype=np.int8))}
 37    counts = [1]
 38    for i in range(n):
 39        nxt = {}
 40        for (used, hist), (cost, prefix) in states.items():
 41            for bit in (0, 1):
 42                nu = used + bit
 43                if nu > k: continue
 44                add = d[i] * bit
 45                for dist in range(1, min(width, i) + 1):
 46                    if bit and ((hist >> (dist-1)) & 1):
 47                        add += Q[i-dist, i]
 48                nh = ((hist << 1) | bit) & ((1 << width) - 1) if width else 0
 49                key = (nu, nh); val = cost + add
 50                if key not in nxt or val < nxt[key][0]:
 51                    nxt[key] = (val, np.append(prefix, bit))
 52        states = nxt; counts.append(len(states))
 53    candidates = [(v[0], v[1]) for key, v in states.items() if key[0] == k]
 54    _, mask = min(candidates, key=lambda x: x[0])
 55    return objective(mask, d, Q), counts, max(counts), mask
 56
 57
 58def magnitude(d, Q, k):
 59    z = np.zeros(len(d), dtype=np.int8)
 60    z[np.argsort(d)[:k]] = 1
 61    return objective(z, d, Q)
 62
 63
 64def path_boundary(n, eta):
 65    return min(n, 2 * eta)
 66
 67
 68def grid_boundary(side, eta):
 69    n = side * side; best = 0
 70    for cut in range(n + 1):
 71        near = set()
 72        for a in range(cut):
 73            ra, ca = divmod(a, side)
 74            for b in range(cut, n):
 75                rb, cb = divmod(b, side)
 76                if abs(ra-rb) + abs(ca-cb) <= eta:
 77                    near.update((a, b))
 78        best = max(best, len(near))
 79    return best
 80
 81
 82def main():
 83    n, k, seed = 20, 10, 7
 84    d0, Q0 = make_problem(n, 0.8, seed)
 85    exact, exact_mask = exact_optimum(d0, Q0, k)
 86    eta_sweep = []
 87    for eta in range(0, 11):
 88        approx, counts, peak, mask = boundary_dp(d0, Q0, k, eta)
 89        eta_sweep.append({'eta': eta, 'full_objective_gap': approx-exact,
 90                          'peak_states': peak, 'predicted_boundary': path_boundary(n, eta),
 91                          'theorem_upper_bound': (n+1) * 2**path_boundary(n, eta)})
 92
 93    rho_sweep = []
 94    for rho in [0.0, 0.2, 0.4, 0.6, 0.8, 0.9]:
 95        d, Q = make_problem(n, rho, seed)
 96        ex, _ = exact_optimum(d, Q, k)
 97        ap, _, peak, _ = boundary_dp(d, Q, k, 3)
 98        rho_sweep.append({'rho': rho, 'eta': 3, 'full_objective_gap': ap-ex,
 99                          'peak_states': peak, 'omitted_tail_scale': rho**4})
100
101    size_sweep = []
102    for nn in [10, 14, 18, 22, 26, 30, 40, 50]:
103        dd, QQ = make_problem(nn, 0.8, seed); kk = nn // 2
104        _, _, peak, _ = boundary_dp(dd, QQ, kk, 3)
105        size_sweep.append({'n': nn, 'eta': 3, 'peak_states': peak,
106                           'states_per_n': peak / nn})
107
108    graph = [{'eta': eta, 'path_boundary': path_boundary(n, eta),
109              'grid8_boundary': grid_boundary(8, eta), 'dense_boundary': n,
110              'path_state_bound': (n+1)*2**path_boundary(n, eta),
111              'grid_state_bound': (n+1)*2**grid_boundary(8, eta)}
112             for eta in [1, 2, 3, 4]]
113    baseline = magnitude(d0, Q0, k)
114    out = {'config': {'n': n, 'k': k, 'seed': seed},
115           'exact_objective': exact, 'magnitude_objective': baseline,
116           'magnitude_gap': baseline-exact, 'eta_sweep': eta_sweep,
117           'rho_sweep': rho_sweep, 'size_sweep': size_sweep,
118           'graph_boundary_sweep': graph}
119    with open('results.json', 'w') as f: json.dump(out, f, indent=2)
120    print(json.dumps(out, indent=2))
121
122if __name__ == '__main__': main()