1import json, math, random, time
  2from collections import defaultdict, deque
  3import numpy as np
  4
  5SEED = 2954
  6rng = np.random.default_rng(SEED)
  7
  8
  9def build_graph(n, kind, density=0.35):
 10    # Bipartite vertices are x=0..n-1 and y=n..2n-1.
 11    edges = []
 12    if kind == "chain":
 13        # A connected sparse graph with many local cycles.
 14        for x in range(n):
 15            edges.append((x, x))
 16            if x + 1 < n: edges.append((x, x + 1))
 17            if x + 2 < n: edges.append((x, x + 2))
 18    elif kind == "grid":
 19        # Dense enough connected random support, with reproducible seed.
 20        for x in range(n):
 21            for y in range(n):
 22                if (x == y or (x + 1) % n == y or (x + 2) % n == y or rng.random() < density):
 23                    edges.append((x, y))
 24    else:
 25        for x in range(n):
 26            for y in range(n):
 27                if rng.random() < density or (x == y): edges.append((x, y))
 28    # remove duplicates and isolate-free components are handled by forest
 29    return list(dict.fromkeys(edges))
 30
 31
 32def fundamental_basis(n, edges):
 33    """Return cycle coefficient vectors and cycle rank for a bipartite graph.
 34    Each vector gives signed edge coefficients for the oriented edge x -> y.
 35    """
 36    m = len(edges)
 37    adj = defaultdict(list)
 38    for i, (x, y) in enumerate(edges):
 39        u, v = x, n + y
 40        adj[u].append((v, i, 1))       # traversing x -> y agrees with edge orientation
 41        adj[v].append((u, i, -1))      # traversing y -> x reverses it
 42    parent = {}
 43    parent_edge = {}
 44    parent_sign = {}
 45    component = {}
 46    tree = set()
 47    for root in list(range(2*n)):
 48        if root in component or not adj[root]: continue
 49        component[root] = root
 50        stack = [root]
 51        while stack:
 52            u = stack.pop()
 53            for v, ei, sign in adj[u]:
 54                if v not in component:
 55                    component[v] = root
 56                    parent[v] = u
 57                    parent_edge[v] = ei
 58                    parent_sign[v] = sign
 59                    tree.add(ei)
 60                    stack.append(v)
 61    non_tree = [i for i in range(m) if i not in tree]
 62    cycles = []
 63    # Return the tree path start -> end with signs relative to x -> y.
 64    def tree_path(start, end):
 65        anc = {}
 66        u = start
 67        while True:
 68            anc[u] = True
 69            if u not in parent: break
 70            u = parent[u]
 71        start_up = []
 72        u = start
 73        end_up = []
 74        u = end
 75        while u not in anc:
 76            end_up.append((parent_edge[u], parent_sign[u]))
 77            u = parent[u]
 78        lca = u
 79        start_up = []
 80        u = start
 81        while u != lca:
 82            start_up.append((parent_edge[u], -parent_sign[u]))
 83            u = parent[u]
 84        end_down = [(ei, sign) for ei, sign in reversed(end_up)]
 85        return start_up + end_down
 86
 87    cycles = []
 88    for ei in non_tree:
 89        x, y = edges[ei]
 90        coeff = np.zeros(m, dtype=float)
 91        coeff[ei] = 1.0  # x -> y
 92        for ej, sign in tree_path(n + y, x):
 93            coeff[ej] += sign
 94        cycles.append(coeff)
 95    active_vertices = len(component)
 96    comps = len(set(component.values()))
 97    rank = m - active_vertices + comps
 98    return np.asarray(cycles), rank
 99
100
101def exhaustive_simple_cycles(n, edges, cap=5000):
102    """Enumerate oriented simple-cycle coefficient vectors, deduplicated."""
103    adj = defaultdict(list)
104    for i, (x, y) in enumerate(edges):
105        u, v = x, n + y
106        adj[u].append((v, i, 1))
107        adj[v].append((u, i, -1))
108    found = {}
109    def dfs(start, u, visited, coeff):
110        if len(found) >= cap:
111            return
112        for v, ei, sign in adj[u]:
113            if v == start and len(visited) >= 4:
114                key = tuple(sorted(np.flatnonzero(coeff).tolist()))
115                if key not in found:
116                    c = coeff.copy()
117                    c[ei] += sign
118                    found[key] = c
119            elif v not in visited and v >= start:
120                coeff[ei] += sign
121                dfs(start, v, visited | {v}, coeff)
122                coeff[ei] -= sign
123    for start in range(2*n):
124        dfs(start, start, {start}, np.zeros(len(edges)))
125    return list(found.values())
126
127
128def main():
129    np.set_printoptions(precision=6, suppress=True)
130    report = {"seed": SEED, "graphs": []}
131    for kind,n in [("chain",12),("grid",9),("random",12)]:
132        edges=build_graph(n,kind)
133        basis,rank=fundamental_basis(n,edges)
134        # Compatible field a = node potential on x minus node potential on y.
135        ux=rng.normal(size=n); vy=rng.normal(size=n)
136        a=np.array([ux[x]-vy[y] for x,y in edges])
137        assert np.max(np.abs(basis @ a)) < 1e-10, "invalid fundamental-cycle orientation"
138        compatible_max=float(np.max(np.abs(basis @ a))) if len(basis) else 0.0
139        # Perturbed field: basis residual is the exact independent constraint vector.
140        noise=rng.normal(scale=.1,size=len(edges)); ap=a+noise
141        bres=basis@ap
142        # Rank and cycle-space reconstruction: random cycle vectors are row-space combinations.
143        if len(basis):
144            coeff=rng.normal(size=(min(20,len(basis)*2),len(basis)))
145            held=coeff@bres
146            held95=float(np.percentile(np.abs(held),95))
147            basis95=float(np.percentile(np.abs(bres),95))
148        else: held95=basis95=0.0
149        exhaustive=exhaustive_simple_cycles(n,edges)
150        exhaustive_compat_max = (max((abs(float(c @ a)) for c in exhaustive), default=0.0))
151        exhaustive_perturbed_95 = (float(np.percentile([abs(float(c @ ap)) for c in exhaustive], 95))
152                                   if exhaustive else 0.0)
153        # Cost proxy is number of residual scalar sums, plus actual timing.
154        reps=10
155        t0=time.perf_counter()
156        for _ in range(reps): _=basis@ap
157        tb=time.perf_counter()-t0
158        t0=time.perf_counter()
159        # Exhaustive cycle residuals represented by summing edge values; count only.
160        for _ in range(reps):
161            _=[float(c @ ap) for c in exhaustive]
162        te=time.perf_counter()-t0
163        report["graphs"].append({"kind":kind,"n":n,"edges":len(edges),"cycle_rank":rank,
164          "basis_constraints":len(basis),"simple_cycles":len(exhaustive),
165          "compatible_max_basis_residual":compatible_max,"exhaustive_compatible_max":exhaustive_compat_max,
166          "perturbed_basis_95":basis95,"exhaustive_perturbed_95":exhaustive_perturbed_95,
167          "heldout_linear_combination_95":held95,"basis_time_sec":tb,"exhaustive_time_sec":te,
168          "speedup":(te/tb if tb else None)})
169    print(json.dumps(report, indent=2))
170
171if __name__ == "__main__": main()