1import json, math
 2import numpy as np
 3from run_experiment import make_graph, collision_count, pairwise_c
 4
 5SEED = 421
 6A, _ = make_graph()
 7n = len(A)
 8incoming = A.T.astype(np.uint8)
 9rng = np.random.default_rng(SEED + 20)
10c = pairwise_c(A)
11
12def fast_greedy(max_s):
13    remaining = np.arange(n)
14    chosen = []
15    Z = np.zeros((n, 0), dtype=np.uint8)
16    for _ in range(min(max_s, n)):
17        if collision_count(Z) == 0:
18            break
19        _, inv = np.unique(Z, axis=0, return_inverse=True)
20        groups = [np.flatnonzero(inv == g) for g in range(int(inv.max()) + 1)]
21        cand = incoming[:, remaining]
22        scores = np.zeros(len(remaining), dtype=np.int64)
23        for inds in groups:
24            if len(inds) > 1:
25                a = cand[inds].sum(axis=0, dtype=np.int64)
26                scores += a * (len(inds) - a)
27        j = int(np.argmax(scores))
28        if scores[j] <= 0:
29            break
30        q = int(remaining[j])
31        chosen.append(q)
32        remaining = np.delete(remaining, j)
33        Z = np.column_stack((Z, incoming[:, q]))
34    return chosen
35
36rows = []
37for s in [8, 16, 24, 32, 40, 49, 64, 96, 128, 160]:
38    trials = 100 if s < 100 else 50
39    vals = []
40    for _ in range(trials):
41        S = rng.choice(n, s, replace=False)
42        vals.append(collision_count(incoming[:, S]))
43    Sg = fast_greedy(s)
44    bound = math.comb(n, 2) * (math.comb(n-c, s) / math.comb(n, s) if n-c >= s else 0.0)
45    rows.append({'s': s, 'uniform_mean_collisions': float(np.mean(vals)),
46                 'uniform_zero_fraction': float(np.mean(np.asarray(vals) == 0)),
47                 'greedy_landmarks_used': len(Sg),
48                 'greedy_collisions': collision_count(incoming[:, Sg]),
49                 'bound_using_c': bound})
50out = {'n': n, 'c_exact': c,
51       'theoretical_s_scale': math.ceil(2*n*math.log(n)/c), 'sweep': rows}
52with open('sweep_results.json', 'w') as f:
53    json.dump(out, f, indent=2)
54print(json.dumps(out, indent=2))