diff --git a/Cargo.lock b/Cargo.lock index 9bfa22a..ee4c9b8 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -455,7 +455,7 @@ checksum = "d6790f58c7ff633d8771f42965289203411a5e5c68388703c06e14f24770b41e" [[package]] name = "optimal-branching" version = "0.1.0" -source = "git+https://github.com/OptimalBranching/optimal-branching.git#5a7749eca4bbd4da5ad1fc5dfaa7f3973dcc829a" +source = "git+https://github.com/OptimalBranching/optimal-branching.git#51fb29d94c46c12bb0f1b692d48c129cbeb925c7" dependencies = [ "good_lp", "num-traits", diff --git a/Cargo.toml b/Cargo.toml index 2c39941..db97513 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -7,6 +7,9 @@ license = "MIT" [dependencies] # Public OptimalBranching/optimal-branching crate (published as `optimal-branching`), # aliased via `package` to the `optimal_branching_core` name our code imports. +# Tracks main; the greedymerge size_reduction memoization (bit-identical, removes +# the redundant re-scoring our GreedyMerge path triggered) landed in main via +# OptimalBranching/optimal-branching#2. Cargo.lock pins the exact rev. optimal-branching-core = { git = "https://github.com/OptimalBranching/optimal-branching.git", package = "optimal-branching" } rustc-hash = "2" smallvec = "1" diff --git a/benchmarks/csp/README.md b/benchmarks/csp/README.md new file mode 100644 index 0000000..906f430 --- /dev/null +++ b/benchmarks/csp/README.md @@ -0,0 +1,26 @@ +# CSP benchmark generators + +Instance generators probing **where region-contraction (structural, local +inference) can beat flat CDCL** — decision families whose structure our solver +sees but a CNF-only solver does not. Each generator emits the same system two +ways where applicable: a native `.csp` (one tensor per constraint, relation seen +whole) and a DIMACS `.cnf` (the flattened encoding a CDCL solver gets). + +| Script | Family | Why it might favor structural inference | +|---|---|---| +| `gen_extcsp.py` | Bounded-treewidth extensional (random table) CSP | Random tables make unit propagation weak; bounded width keeps exact local inference cheap. | +| `gen_xor.py` | Random k-XOR-SAT | Parity is exponentially hard for resolution/CDCL (no Gaussian elimination), linear-algebra-easy. | +| `gen_tseitin.py` | Tseitin parity formulas | Resolution-hard on expanders; poly via local elimination on bounded-treewidth graphs. | +| `gen_coloring.py` | k-coloring, banded vs random graphs | Low-treewidth (banded) graphs keep regions local; tests width exploitation. | +| `gen_qcp.py` | Quasigroup completion (Latin square) | Classic CP-beats-SAT; GAC on wide exactly-one tensors. | +| `gen_x3c.py` | Exact Cover by 3-Sets | Home family for the exact-cover γ branching rules. | +| `cp_sat_solve.py` | OR-Tools CP-SAT baseline | Fair SOTA CP baseline: solves the native `.csp` relations directly (no CNF handicap). | + +Native `.csp` format: line 1 = `n_vars`; then per tensor +` ... : ...`, where each `cfg` is a bitmask over the scope +listing an allowed assignment. + +Status: exploratory. The current paper direction converged on **uniform cubing +for parallel Cube-and-Conquer on arithmetic circuits**, not on a CSP decision +win — these generators are retained as reproducible probes for revisiting the +CSP question, not as a claimed result. diff --git a/benchmarks/csp/cp_sat_solve.py b/benchmarks/csp/cp_sat_solve.py new file mode 100644 index 0000000..b0d116a --- /dev/null +++ b/benchmarks/csp/cp_sat_solve.py @@ -0,0 +1,48 @@ +#!/usr/bin/env python3 +"""Solve a native .csp (our tensor format) with OR-Tools CP-SAT — the SOTA CP +solver — using AddAllowedAssignments per tensor, i.e. the SAME relations our +solver sees, modelled natively (no CNF handicap). This is the fair SOTA baseline +for goal D: if our region contraction can't beat CP-SAT here, there is no CSP +advantage over the state of the art. + +Reads the .csp format: line 1 = n_vars; then per tensor +" ... : ...", cfg = bitmask over the scope. + +Prints: result, wall_seconds, branches (CP-SAT's search-node analogue). +Usage: cp_sat_solve.py +""" +import sys +import time + +from ortools.sat.python import cp_model + + +def main(): + path = sys.argv[1] + with open(path) as f: + lines = [l for l in f.read().splitlines() if l.strip()] + n = int(lines[0]) + model = cp_model.CpModel() + x = [model.NewBoolVar(f"x{i}") for i in range(n)] + + for line in lines[1:]: + lhs, rhs = line.split(":") + scope = [int(s) for s in lhs.split()] + allowed_cfgs = [int(s) for s in rhs.split()] + tuples = [] + for c in allowed_cfgs: + tuples.append([(c >> i) & 1 for i in range(len(scope))]) + model.AddAllowedAssignments([x[v] for v in scope], tuples) + + solver = cp_model.CpSolver() + solver.parameters.num_search_workers = 1 # single-threaded, fair vs our solver + t0 = time.perf_counter() + status = solver.Solve(model) + dt = time.perf_counter() - t0 + res = {cp_model.OPTIMAL: "SAT", cp_model.FEASIBLE: "SAT", + cp_model.INFEASIBLE: "UNSAT"}.get(status, "?") + print(f"result={res} time={dt:.3f}s branches={solver.NumBranches()} conflicts={solver.NumConflicts()}") + + +if __name__ == "__main__": + main() diff --git a/benchmarks/csp/gen_coloring.py b/benchmarks/csp/gen_coloring.py new file mode 100644 index 0000000..bc72f9d --- /dev/null +++ b/benchmarks/csp/gen_coloring.py @@ -0,0 +1,102 @@ +#!/usr/bin/env python3 +"""Generate k-coloring DIMACS CNF instances, with a graph-structure knob. + +Goal (runscribe goal D): find a decision family where our region-contraction +solver's LOCALITY exploitation beats CDCL. Graph coloring on LOW-TREEWIDTH +graphs is the mechanism bet — a region grown around a vertex stays a small +local patch, so the branching tree tracks the tree-decomposition width, while +CDCL treats the CNF as flat and cannot exploit the bounded width structurally. + +Two graph families, same n and edge budget so hardness is comparable: + banded — vertices on a line; edges only within bandwidth `w`. Treewidth <= w + (a path/band decomposition), i.e. bounded, LOCAL structure. + random — Erdos-Renyi G(n, m): the SAME number of edges placed uniformly at + random. Treewidth ~ n; no local structure. CDCL's home turf. + +Encoding (standard direct k-coloring): x[v,c] = "vertex v has color c". + - at-least-one: (x[v,0] v ... v x[v,k-1]) per vertex + - at-most-one: (~x[v,c] v ~x[v,c']) per vertex, c out.cnf + gen_coloring.py --family random --n 40 --k 4 --edges 148 --seed 1 > out.cnf +""" +import argparse +import random +import sys + + +def banded_edges(n, w, rng): + """Edges within bandwidth w on a line 0..n-1. Includes the path (i,i+1) so + the graph is connected, then fills in-band chords. Treewidth <= w.""" + edges = set() + for i in range(n): + for j in range(i + 1, min(i + w + 1, n)): + edges.add((i, j)) + return edges + + +def random_edges(n, m, rng): + """m distinct random edges (Erdos-Renyi G(n,m)).""" + edges = set() + m = min(m, n * (n - 1) // 2) + while len(edges) < m: + u = rng.randrange(n) + v = rng.randrange(n) + if u != v: + edges.add((min(u, v), max(u, v))) + return edges + + +def main(): + ap = argparse.ArgumentParser() + ap.add_argument("--family", choices=["banded", "random"], required=True) + ap.add_argument("--n", type=int, required=True, help="vertices") + ap.add_argument("--k", type=int, required=True, help="colors") + ap.add_argument("--w", type=int, default=4, help="bandwidth (banded only)") + ap.add_argument("--edges", type=int, default=0, help="edge count (random; default=match banded w)") + ap.add_argument("--seed", type=int, default=1) + ap.add_argument("--dump-graph", default="", help="also write ' ' + edge list here") + args = ap.parse_args() + + rng = random.Random(args.seed * 100003 + args.n) + if args.family == "banded": + edges = banded_edges(args.n, args.w, rng) + else: + m = args.edges or (args.n * args.w - args.w * (args.w + 1) // 2) # ~banded count + edges = random_edges(args.n, m, rng) + + n, k = args.n, args.k + + if args.dump_graph: + with open(args.dump_graph, "w") as f: + f.write(f"{n} {k}\n") + for (u, v) in sorted(edges): + f.write(f"{u} {v}\n") + + def var(v, c): # 1-based DIMACS var id + return v * k + c + 1 + + clauses = [] + for v in range(n): + clauses.append([var(v, c) for c in range(k)]) # at-least-one + for c in range(k): + for c2 in range(c + 1, k): + clauses.append([-var(v, c), -var(v, c2)]) # at-most-one + for (u, v) in sorted(edges): + for c in range(k): + clauses.append([-var(u, c), -var(v, c)]) # edge differ + + out = [f"c k-coloring family={args.family} n={n} k={k} " + f"w={args.w} edges={len(edges)} seed={args.seed}", + f"p cnf {n * k} {len(clauses)}"] + for cl in clauses: + out.append(" ".join(map(str, cl)) + " 0") + sys.stdout.write("\n".join(out) + "\n") + + +if __name__ == "__main__": + main() diff --git a/benchmarks/csp/gen_extcsp.py b/benchmarks/csp/gen_extcsp.py new file mode 100644 index 0000000..e533d1d --- /dev/null +++ b/benchmarks/csp/gen_extcsp.py @@ -0,0 +1,73 @@ +#!/usr/bin/env python3 +"""Bounded-treewidth extensional (table) CSP — the theoretical home of +structural constraint solving. n boolean vars laid on a line; each constraint's +scope is a random arity-`a` subset drawn from a sliding window of width `w`, so +the constraint graph has treewidth <= w (a path/band decomposition). Each +constraint is an explicit RANDOM table: each of the 2^a tuples is allowed +independently with prob (1 - tightness). Random tables make unit propagation +WEAK (no clause is short), while the bounded width keeps exact inference +(our region contraction / variable elimination) cheap — the regime where a +structural solver should beat CDCL, which sees only the flat encoding. + +Native .csp: each constraint = one arity-a tensor (allowed tuples). CNF: each +FORBIDDEN tuple = one clause (the direct/support encoding CDCL must use). + +Usage: gen_extcsp.py --n 60 --w 8 --a 3 --m 90 --tight 0.5 --seed 1 \ + --csp o.csp --cnf o.cnf +""" +import argparse +import random + + +def main(): + ap = argparse.ArgumentParser() + ap.add_argument("--n", type=int, required=True) + ap.add_argument("--w", type=int, default=8, help="window width (~treewidth bound)") + ap.add_argument("--a", type=int, default=3, help="constraint arity") + ap.add_argument("--m", type=int, default=0, help="constraints (default 1.5n)") + ap.add_argument("--tight", type=float, default=0.5, help="fraction of forbidden tuples") + ap.add_argument("--seed", type=int, default=1) + ap.add_argument("--csp", required=True) + ap.add_argument("--cnf", required=True) + args = ap.parse_args() + + rng = random.Random(args.seed * 100003 + args.n) + n, w, a = args.n, args.w, args.a + m = args.m or round(1.5 * n) + + cons = [] # (scope, allowed set) + for _ in range(m): + start = rng.randrange(max(1, n - w + 1)) + window = list(range(start, min(start + w, n))) + if len(window) < a: + window = list(range(n)) + scope = sorted(rng.sample(window, a)) + allowed = [c for c in range(1 << a) if rng.random() > args.tight] + if not allowed: # never emit an empty (trivially-UNSAT) table + allowed = [rng.randrange(1 << a)] + cons.append((scope, allowed)) + + with open(args.csp, "w") as f: + f.write(f"{n}\n") + for scope, allowed in cons: + f.write(" ".join(map(str, scope)) + " : " + " ".join(map(str, allowed)) + "\n") + + clauses = [] + for scope, allowed in cons: + aset = set(allowed) + for c in range(1 << a): + if c not in aset: # forbid tuple c: one clause + lits = [] + for i, v in enumerate(scope): + bit = (c >> i) & 1 + lits.append(-(v + 1) if bit else (v + 1)) + clauses.append(lits) + with open(args.cnf, "w") as f: + f.write(f"c extcsp n={n} w={w} a={a} m={m} tight={args.tight} seed={args.seed}\n") + f.write(f"p cnf {n} {len(clauses)}\n") + for cl in clauses: + f.write(" ".join(map(str, cl)) + " 0\n") + + +if __name__ == "__main__": + main() diff --git a/benchmarks/csp/gen_qcp.py b/benchmarks/csp/gen_qcp.py new file mode 100644 index 0000000..febc8b9 --- /dev/null +++ b/benchmarks/csp/gen_qcp.py @@ -0,0 +1,91 @@ +#!/usr/bin/env python3 +"""Quasigroup completion (Latin square with holes) — the classic CP-beats-SAT +family. Fill an m x m grid so each row and column is a permutation of 0..m-1, +given some prefilled cells. Generated from a full Latin square with holes +punched, so always SAT; hardest around ~60% holes. + +Native .csp uses the standard DECOMPOSED all-different (exactly-one per cell, +per row-value, per col-value) — the SAME decomposition SAT uses, because the +global all-different does not compress into a small boolean tensor. So this +tests whether our GAC on wide exactly-one tensors beats CDCL even without the +true all-different global propagation. + +Usage: gen_qcp.py --m 8 --holes 0.6 --seed 1 --csp o.csp --cnf o.cnf +""" +import argparse +import random + + +def main(): + ap = argparse.ArgumentParser() + ap.add_argument("--m", type=int, required=True) + ap.add_argument("--holes", type=float, default=0.6) + ap.add_argument("--seed", type=int, default=1) + ap.add_argument("--csp", required=True) + ap.add_argument("--cnf", required=True) + args = ap.parse_args() + + rng = random.Random(args.seed * 100003 + args.m) + m = args.m + # Full Latin square via cyclic construction + random row/col/symbol permutations. + base = [[(i + j) % m for j in range(m)] for i in range(m)] + rp = list(range(m)); rng.shuffle(rp) + cp = list(range(m)); rng.shuffle(cp) + sp = list(range(m)); rng.shuffle(sp) + L = [[sp[base[rp[i]][cp[j]]] for j in range(m)] for i in range(m)] + # Punch holes. + given = {} + for i in range(m): + for j in range(m): + if rng.random() > args.holes: + given[(i, j)] = L[i][j] + + def var(r, c, v): + return (r * m + c) * m + v # 0-based + + n_vars = m * m * m + scopes, tables = [], [] + full = list(range(1 << m)) + one_hot = [c for c in full if bin(c).count("1") == 1] + + def exactly_one(ids): + scopes.append(list(ids)) + tables.append(one_hot) + + for r in range(m): + for c in range(m): + exactly_one([var(r, c, v) for v in range(m)]) # one value per cell + for r in range(m): + for v in range(m): + exactly_one([var(r, c, v) for c in range(m)]) # value once per row + for c in range(m): + for v in range(m): + exactly_one([var(r, c, v) for r in range(m)]) # value once per col + + # CNF (1-based): exactly-one = at-least-one + at-most-one pairs; givens as units. + clauses = [] + for sc in scopes: + clauses.append([x + 1 for x in sc]) + for a in range(len(sc)): + for b in range(a + 1, len(sc)): + clauses.append([-(sc[a] + 1), -(sc[b] + 1)]) + given_units = [var(r, c, v) for (r, c), v in given.items()] + for u in given_units: + clauses.append([u + 1]) + + with open(args.csp, "w") as f: + f.write(f"{n_vars}\n") + for sc, tb in zip(scopes, tables): + f.write(" ".join(map(str, sc)) + " : " + " ".join(map(str, tb)) + "\n") + for u in given_units: # given: unit tensor fixing var to 1 + f.write(f"{u} : 1\n") + + with open(args.cnf, "w") as f: + f.write(f"c qcp m={m} holes={args.holes} given={len(given)} seed={args.seed}\n") + f.write(f"p cnf {n_vars} {len(clauses)}\n") + for cl in clauses: + f.write(" ".join(map(str, cl)) + " 0\n") + + +if __name__ == "__main__": + main() diff --git a/benchmarks/csp/gen_tseitin.py b/benchmarks/csp/gen_tseitin.py new file mode 100644 index 0000000..750c87b --- /dev/null +++ b/benchmarks/csp/gen_tseitin.py @@ -0,0 +1,119 @@ +#!/usr/bin/env python3 +"""Tseitin formulas — the canonical resolution-hard family. Each graph EDGE is a +boolean var; each VERTEX gets a parity constraint (XOR of its incident edge vars += its charge). If the total charge is odd the formula is UNSAT, and on an +EXPANDER graph every resolution/CDCL refutation is exponentially long — kissat's +worst case. On a bounded-treewidth graph (grid/path) variable elimination is +poly. This probes whether our region contraction (local elimination) beats CDCL +where CDCL is provably weak. + +Families: + regular — random d-regular graph on v vertices (expander w.h.p. for d>=3) + grid — r x c grid (treewidth ~ min(r,c); bounded if one side is small) + +Emits native .csp (one parity tensor per vertex) and DIMACS .cnf (each parity +expands to 2^(deg-1) clauses), same system both ways. + +Usage: + gen_tseitin.py --family regular --v 40 --d 3 --seed 1 --csp o.csp --cnf o.cnf + gen_tseitin.py --family grid --rows 4 --cols 20 --seed 1 --csp o.csp --cnf o.cnf +""" +import argparse +import random + + +def parity(c): + return bin(c).count("1") & 1 + + +def regular_graph(v, d, rng): + """Random d-regular graph via the pairing model (retry on failure).""" + for _ in range(200): + stubs = [x for x in range(v) for _ in range(d)] + rng.shuffle(stubs) + edges = set() + ok = True + for i in range(0, len(stubs), 2): + a, b = stubs[i], stubs[i + 1] + if a == b or (min(a, b), max(a, b)) in edges: + ok = False + break + edges.add((min(a, b), max(a, b))) + if ok: + return sorted(edges) + raise RuntimeError("failed to build regular graph") + + +def grid_graph(r, c): + def vid(i, j): + return i * c + j + edges = [] + for i in range(r): + for j in range(c): + if j + 1 < c: + edges.append((vid(i, j), vid(i, j + 1))) + if i + 1 < r: + edges.append((vid(i, j), vid(i + 1, j))) + return sorted(edges), r * c + + +def main(): + ap = argparse.ArgumentParser() + ap.add_argument("--family", choices=["regular", "grid"], required=True) + ap.add_argument("--v", type=int, default=40) + ap.add_argument("--d", type=int, default=3) + ap.add_argument("--rows", type=int, default=4) + ap.add_argument("--cols", type=int, default=20) + ap.add_argument("--seed", type=int, default=1) + ap.add_argument("--csp", required=True) + ap.add_argument("--cnf", required=True) + args = ap.parse_args() + + rng = random.Random(args.seed * 100003 + args.v + args.cols) + if args.family == "regular": + edges = regular_graph(args.v, args.d, rng) + nv = args.v + else: + edges, nv = grid_graph(args.rows, args.cols) + + # edge var id = index in `edges`. Incident edges per vertex. + inc = {v: [] for v in range(nv)} + for ei, (a, b) in enumerate(edges): + inc[a].append(ei) + inc[b].append(ei) + + # charges: all 0 except one vertex flipped to 1 -> total charge odd -> UNSAT. + charge = {v: 0 for v in range(nv)} + charge[0] = 1 + + n_edges = len(edges) + # native: one parity tensor per vertex over its incident edge vars. + csp_lines = [str(n_edges)] + cnf_clauses = [] + for v in range(nv): + scope = inc[v] + if not scope: + continue + b = charge[v] + k = len(scope) + allowed = [c for c in range(1 << k) if parity(c) == b] + csp_lines.append(" ".join(map(str, scope)) + " : " + " ".join(map(str, allowed))) + for c in range(1 << k): + if parity(c) != b: + lits = [] + for i, e in enumerate(scope): + bit = (c >> i) & 1 + lits.append(-(e + 1) if bit else (e + 1)) + cnf_clauses.append(lits) + + with open(args.csp, "w") as f: + f.write("\n".join(csp_lines) + "\n") + with open(args.cnf, "w") as f: + f.write(f"c tseitin {args.family} edges={n_edges} verts={nv} seed={args.seed}\n") + f.write(f"p cnf {n_edges} {len(cnf_clauses)}\n") + for cl in cnf_clauses: + f.write(" ".join(map(str, cl)) + " 0\n") + + +if __name__ == "__main__": + main() diff --git a/benchmarks/csp/gen_x3c.py b/benchmarks/csp/gen_x3c.py new file mode 100644 index 0000000..d43fe2b --- /dev/null +++ b/benchmarks/csp/gen_x3c.py @@ -0,0 +1,70 @@ +#!/usr/bin/env python3 +"""Exact Cover by 3-Sets (X3C) decision — the family the project's novelty (the +exact-cover gamma branching rules, issue #14 / #8-#13) is built for. Universe of +3q elements; a collection of 3-element subsets; is there a subcollection that +partitions the universe (covers each element EXACTLY once)? + +A planted exact cover guarantees SAT; `--extra` random 3-subsets are distractors +that make search hard. Constraint per element: exactly one selected subset +covers it (arity = #subsets containing it). Emits native .csp (one exactly-one +tensor per element) and DIMACS .cnf (at-least-one + at-most-one pairs). + +Usage: gen_x3c.py --q 20 --extra 40 --seed 1 --csp o.csp --cnf o.cnf +""" +import argparse +import random + + +def main(): + ap = argparse.ArgumentParser() + ap.add_argument("--q", type=int, required=True, help="universe = 3q elements") + ap.add_argument("--extra", type=int, default=0, help="distractor 3-subsets") + ap.add_argument("--seed", type=int, default=1) + ap.add_argument("--csp", required=True) + ap.add_argument("--cnf", required=True) + args = ap.parse_args() + + rng = random.Random(args.seed * 100003 + args.q) + n = 3 * args.q + elems = list(range(n)) + rng.shuffle(elems) + subsets = [tuple(sorted(elems[3 * i:3 * i + 3])) for i in range(args.q)] # planted cover + seen = set(subsets) + while len(subsets) < args.q + args.extra: + s = tuple(sorted(rng.sample(range(n), 3))) + if s not in seen: + seen.add(s) + subsets.append(s) + rng.shuffle(subsets) + + m = len(subsets) # subset vars, 0-based + covers = {e: [] for e in range(n)} # element -> subset ids containing it + for si, s in enumerate(subsets): + for e in s: + covers[e].append(si) + + def one_hot(k): + return [c for c in range(1 << k) if bin(c).count("1") == 1] + + with open(args.csp, "w") as f: + f.write(f"{m}\n") + for e in range(n): + sc = covers[e] + f.write(" ".join(map(str, sc)) + " : " + " ".join(map(str, one_hot(len(sc)))) + "\n") + + clauses = [] + for e in range(n): + sc = covers[e] + clauses.append([s + 1 for s in sc]) # at-least-one + for a in range(len(sc)): + for b in range(a + 1, len(sc)): + clauses.append([-(sc[a] + 1), -(sc[b] + 1)]) # at-most-one + with open(args.cnf, "w") as f: + f.write(f"c x3c q={args.q} subsets={m} extra={args.extra} seed={args.seed}\n") + f.write(f"p cnf {m} {len(clauses)}\n") + for cl in clauses: + f.write(" ".join(map(str, cl)) + " 0\n") + + +if __name__ == "__main__": + main() diff --git a/benchmarks/csp/gen_xor.py b/benchmarks/csp/gen_xor.py new file mode 100644 index 0000000..eb147bc --- /dev/null +++ b/benchmarks/csp/gen_xor.py @@ -0,0 +1,73 @@ +#!/usr/bin/env python3 +"""Random k-XOR-SAT: m parity equations over n boolean vars, each +x_{i0} XOR ... XOR x_{i(k-1)} = b. Emits the SAME system two ways: + + native .csp — one parity TENSOR per equation (arity k; allowed configs are + those with the right parity). Our solver sees the relation whole. + DIMACS .cnf — each equation expands to 2^(k-1) clauses. This is what a CDCL + solver (kissat) gets. + +Why this family (runscribe goal D): XOR/parity is the textbook case where +CDCL/resolution is EXPONENTIALLY weak (kissat has no Gaussian elimination), +yet the system is linear-algebra-easy. If our region contraction acts like +local elimination over the parity tensors, this is a home-turf win CDCL can't +match. + +Usage: + gen_xor.py --n 60 --ratio 0.9 --k 3 --seed 1 --csp out.csp --cnf out.cnf +""" +import argparse +import random + + +def parity(c): + return bin(c).count("1") & 1 + + +def main(): + ap = argparse.ArgumentParser() + ap.add_argument("--n", type=int, required=True) + ap.add_argument("--k", type=int, default=3) + ap.add_argument("--ratio", type=float, default=0.9, help="m/n") + ap.add_argument("--m", type=int, default=0, help="equations (overrides ratio)") + ap.add_argument("--seed", type=int, default=1) + ap.add_argument("--csp", required=True) + ap.add_argument("--cnf", required=True) + args = ap.parse_args() + + rng = random.Random(args.seed * 100003 + args.n) + n, k = args.n, args.k + m = args.m or max(1, round(args.ratio * n)) + + eqs = [] # (vars0based, b) + for _ in range(m): + vs = rng.sample(range(n), k) + b = rng.randrange(2) + eqs.append((sorted(vs), b)) + + # native: allowed configs = those whose parity == b + with open(args.csp, "w") as f: + f.write(f"{n}\n") + for (vs, b) in eqs: + allowed = [c for c in range(1 << k) if parity(c) == b] + f.write(" ".join(map(str, vs)) + " : " + " ".join(map(str, allowed)) + "\n") + + # cnf: forbid each config with parity != b (one clause each). DIMACS vars 1-based. + clauses = [] + for (vs, b) in eqs: + for c in range(1 << k): + if parity(c) != b: + lits = [] + for i, v in enumerate(vs): + bit = (c >> i) & 1 + lits.append(-(v + 1) if bit else (v + 1)) + clauses.append(lits) + with open(args.cnf, "w") as f: + f.write(f"c random {k}-XOR n={n} m={m} seed={args.seed}\n") + f.write(f"p cnf {n} {len(clauses)}\n") + for cl in clauses: + f.write(" ".join(map(str, cl)) + " 0\n") + + +if __name__ == "__main__": + main() diff --git a/src/adapter.rs b/src/adapter.rs index 88f0456..5b2afbb 100644 --- a/src/adapter.rs +++ b/src/adapter.rs @@ -41,13 +41,6 @@ struct MeasureScratch { tables: Vec, buffer: SolverBuffer, trail: Trail, - /// Per-node memoization of `apply_branch` results keyed by (clause.mask, - /// clause.val). `variables` is constant across a node's `optimal_rule` call, - /// so (mask,val) uniquely identifies a branch from the fixed node base; - /// GreedyMerge's O(rows²) merge loop re-evaluates many clauses, so caching the - /// propagated domains avoids re-propagating them. Cleared per node in - /// `with_measure_scratch`. Measured ~8% wall-clock on factoring_22x22 VE10. - cache: std::collections::HashMap<(u64, u64), Vec>, } thread_local! { @@ -74,7 +67,6 @@ pub(crate) fn with_measure_scratch( std::mem::swap(&mut s.trail, trail); s.doms.clear(); s.doms.extend_from_slice(doms); - s.cache.clear(); }); let r = f(); MEASURE_SCRATCH.with(|s| { @@ -125,35 +117,16 @@ impl BranchAndReduceProblem for RuleProblem { /// reach the same GAC fixpoint) but ~2-3x faster and allocation-free. /// Precondition (ob-core guarantee): called only single-level from the root, /// with the scratch primed by `with_measure_scratch`. + /// + /// No per-node memo here: `GreedyMerge` (the only rule solver that re-evaluates + /// the same clause) now memoizes `size_reduction` by `(mask, val)` in ob-core, + /// upstream of this call and covering the measure too, so a downstream cache + /// would never be hit. `IPSolver`/`LPSolver`/`NaiveBranch` evaluate each + /// candidate clause exactly once, so they never needed one. fn apply_branch(&self, clause: &Clause, variables: &[usize]) -> (RuleProblem, f64) { - let key = (clause.mask, clause.val); let snapshot = MEASURE_SCRATCH.with(|s| { let s = &mut *s.borrow_mut(); - // Hit: the branch was already propagated this node; return the cached - // domains. The scratch is untouched (no trail.open/restore), so it - // stays at base for the next call — correct because the miss path - // always restores and the hit path never mutates. - if let Some(cached) = s.cache.get(&key).cloned() { - // Self-verification (debug only, compiled out of release): a cache - // hit MUST equal a fresh recomputation, i.e. (mask,val) truly - // determines the result from the restored node base. This asserts - // the memoization is behavior-identical on every hit of a full solve. - #[cfg(debug_assertions)] - { - let fresh = apply_branch_fresh(&self.cn, &self.masks, s, clause, variables); - assert_eq!( - fresh, cached, - "MEMO MISMATCH for clause (mask={:#x}, val={:#x}): cached result \ - differs from fresh recomputation — the (mask,val) key does not \ - determine apply_branch's result", - clause.mask, clause.val - ); - } - return cached; - } - let snap = apply_branch_fresh(&self.cn, &self.masks, s, clause, variables); - s.cache.insert(key, snap.clone()); - snap + apply_branch_fresh(&self.cn, &self.masks, s, clause, variables) }); ( RuleProblem { diff --git a/src/measure.rs b/src/measure.rs index a13ce46..bd093f5 100644 --- a/src/measure.rs +++ b/src/measure.rs @@ -1,6 +1,6 @@ use crate::domain::DomainMask; use crate::network::ConstraintNetwork; -use crate::util::get_active_tensors; +use crate::util::active_tensors; #[derive(Clone, Copy, Debug, PartialEq, Eq)] pub enum Measure { @@ -14,10 +14,10 @@ pub enum Measure { pub fn measure_core(cn: &ConstraintNetwork, doms: &[DomainMask], m: Measure) -> f64 { match m { Measure::NumUnfixedVars => doms.iter().filter(|d| !d.is_fixed()).count() as f64, - Measure::NumUnfixedTensors => get_active_tensors(cn, doms).len() as f64, + Measure::NumUnfixedTensors => active_tensors(cn, doms).count() as f64, Measure::NumHardTensors => { let mut excess = 0usize; - for tid in get_active_tensors(cn, doms) { + for tid in active_tensors(cn, doms) { let degree = cn.tensors[tid] .var_axes .iter() diff --git a/src/region.rs b/src/region.rs index f8883ef..14ba5b0 100644 --- a/src/region.rs +++ b/src/region.rs @@ -1,6 +1,6 @@ use std::cmp::Reverse; -use rustc_hash::{FxHashMap, FxHashSet}; +use rustc_hash::FxHashMap; use crate::contract::{join_bounded, tensor_relation, Relation}; use crate::ct::TableMasks; @@ -83,13 +83,15 @@ pub fn grow_region( // borrow — no Relation is cloned just to read a length or feed a join) and // per-tensor entailment (the frontier re-examines the same tensors on // every growth step). + // Both memos take their map as an argument (rather than capturing it) so + // several helpers can share them without fighting the borrow checker. let mut rels: FxHashMap = FxHashMap::default(); + let mut ent: FxHashMap = FxHashMap::default(); let ensure = |tid: usize, rels: &mut FxHashMap| { rels.entry(tid) .or_insert_with(|| tensor_relation(cn, &cn.tensors[tid], doms)); }; - let mut ent: FxHashMap = FxHashMap::default(); - let mut entailed = |tid: usize| -> bool { + let entailed = |tid: usize, ent: &mut FxHashMap| -> bool { *ent.entry(tid) .or_insert_with(|| is_entailed(cn, tid, doms, masks)) }; @@ -102,7 +104,7 @@ pub fn grow_region( let incident: Vec = cn.v2t[focus] .iter() .copied() - .filter(|&tid| !entailed(tid)) + .filter(|&tid| !entailed(tid, &mut ent)) .collect(); for &tid in &incident { ensure(tid, &mut rels); @@ -117,33 +119,51 @@ pub fn grow_region( let mut tensors = vec![seed]; let mut acc = rels[&seed].clone(); - loop { - // Frontier: tensors incident to a region var, not yet absorbed. - let mut frontier: Vec = Vec::new(); - let mut seen: FxHashSet = FxHashSet::default(); - for &v in &acc.vars { + // Frontier maintained INCREMENTALLY across growth steps (was rebuilt from + // scratch every step — the profiled hot spot). `frontier[tid]` is the count + // of `tid`'s vars already in `acc.vars` — the ranking's shared-var key. + // Invariant: `frontier` holds exactly the non-entailed, not-yet-absorbed + // tensors incident to a region var, each with its exact shared count. + // Maintained by `absorb_vars`: every var enters `acc` exactly once (at seed + // or on one absorption), and there it bumps every incident frontier tensor — + // so the running count equals |tid.var_axes ∩ acc.vars| without rescanning. + let mut frontier: FxHashMap = FxHashMap::default(); + let absorb_vars = |vars: &[usize], + frontier: &mut FxHashMap, + rels: &mut FxHashMap, + ent: &mut FxHashMap, + tensors: &[usize]| { + for &v in vars { for &tid in &cn.v2t[v] { - if tensors.binary_search(&tid).is_err() && seen.insert(tid) && !entailed(tid) { - frontier.push(tid); + if tensors.binary_search(&tid).is_ok() || entailed(tid, ent) { + continue; } + ensure(tid, rels); + *frontier.entry(tid).or_insert(0) += 1; } } + }; + absorb_vars( + &acc.vars.clone(), + &mut frontier, + &mut rels, + &mut ent, + &tensors, + ); + + loop { if frontier.is_empty() { break; } - // Cheap ranking: most shared unfixed vars first (tight joins grow rows - // least), then smallest sliced support, then tid for determinism. - let mut scored: Vec<(Reverse, usize, usize)> = Vec::with_capacity(frontier.len()); - for &tid in &frontier { - ensure(tid, &mut rels); - let shared = cn.tensors[tid] - .var_axes - .iter() - .filter(|&&v| acc.vars.binary_search(&v).is_ok()) - .count(); - scored.push((Reverse(shared), rels[&tid].rows.len(), tid)); - } + // least), then smallest sliced support, then tid for determinism. The + // sort key ends in the unique `tid`, so it is a TOTAL order — the + // unordered hashmap iteration below produces a bit-identical ranking to + // the old scratch-rebuilt vector after sorting. + let mut scored: Vec<(Reverse, usize, usize)> = frontier + .iter() + .map(|(&tid, &shared)| (Reverse(shared), rels[&tid].rows.len(), tid)) + .collect(); scored.sort_unstable(); // Exact-join down the ranked list: pick the min-row join among the top @@ -156,7 +176,7 @@ pub fn grow_region( if best.is_some() && scanned >= JOIN_CANDIDATES { break; } - // `tid` was ensured during scoring above, so the memo has it. + // `tid` is in `frontier`, so the memo has it (ensured on insert). // Hard 64-var cap BEFORE joining: u64 rows silently corrupt past it. let rel = &rels[&tid]; let new_vars = rel @@ -183,9 +203,18 @@ pub fn grow_region( } match best { Some((tid, joined)) => { + // Vars this absorption newly adds to the region. + let new_vars: Vec = joined + .vars + .iter() + .copied() + .filter(|v| acc.vars.binary_search(v).is_err()) + .collect(); + frontier.remove(&tid); let pos = tensors.binary_search(&tid).unwrap_or_else(|e| e); tensors.insert(pos, tid); // keep `tensors` sorted acc = joined; + absorb_vars(&new_vars, &mut frontier, &mut rels, &mut ent, &tensors); } None => break, // nothing on the frontier fits: the region is done } diff --git a/src/selector.rs b/src/selector.rs index 908863e..2d8dfa3 100644 --- a/src/selector.rs +++ b/src/selector.rs @@ -10,7 +10,7 @@ use crate::network::ConstraintNetwork; use crate::problem::SolverBuffer; use crate::table::compute_branching_result; use crate::trail::Trail; -use crate::util::{get_active_tensors, is_entailed}; +use crate::util::{active_tensors, is_entailed}; /// Fill `buffer.occurrence_scores`: each active NON-ENTAILED tensor adds 1 to /// each of its unfixed variables — plain occurrence counting, no structural @@ -27,7 +27,7 @@ pub(crate) fn compute_occurrence_scores( for s in buffer.occurrence_scores.iter_mut() { *s = 0.0; } - for tid in get_active_tensors(cn, doms) { + for tid in active_tensors(cn, doms) { if is_entailed(cn, tid, doms, masks) { continue; } diff --git a/src/util.rs b/src/util.rs index 9c2bfdd..0d13f40 100644 --- a/src/util.rs +++ b/src/util.rs @@ -2,16 +2,25 @@ use crate::ct::TableMasks; use crate::domain::DomainMask; use crate::network::ConstraintNetwork; -/// Tensor ids that still have at least one unfixed variable, ascending. -/// Port of `utils.jl::get_active_tensors`. +/// Tensor ids that still have at least one unfixed variable, ascending. An +/// ACTIVE tensor is one whose scope isn't fully fixed, so it can still +/// constrain the search. Yielded lazily: hot callers (occurrence scoring, the +/// tensor measures) only iterate or count, so materializing a `Vec` per call is +/// pure waste. Port of `utils.jl::get_active_tensors`. +pub fn active_tensors<'a>( + cn: &'a ConstraintNetwork, + doms: &'a [DomainMask], +) -> impl Iterator + 'a { + cn.tensors + .iter() + .enumerate() + .filter(move |(_, t)| t.var_axes.iter().any(|&v| !doms[v].is_fixed())) + .map(|(tid, _)| tid) +} + +/// `active_tensors` collected — for the rare caller that needs an owned slice. pub fn get_active_tensors(cn: &ConstraintNetwork, doms: &[DomainMask]) -> Vec { - let mut active = Vec::with_capacity(cn.tensors.len()); - for (tid, t) in cn.tensors.iter().enumerate() { - if t.var_axes.iter().any(|&v| !doms[v].is_fixed()) { - active.push(tid); - } - } - active + active_tensors(cn, doms).collect() } /// Number of unfixed domains. Port of `utils.jl::count_unfixed`.