KISS/DRY cleanup: hot-path micro-opts, dedup consolidation, CSP bench generators

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"

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
+`<v0> <v1> ... : <cfg0> <cfg1> ...`, 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.

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
+"<v0> <v1> ... : <cfg0> <cfg1> ...", cfg = bitmask over the scope.
+
+Prints: result, wall_seconds, branches (CP-SAT's search-node analogue).
+Usage: cp_sat_solve.py <problem.csp>
+"""
+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()

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<c'
+  - edge-differ:   (~x[u,c] v ~x[v,c])                  per edge, per color
+Same encoding for both families and for our solver + kissat (one CNF, two
+solvers) — no reduction routing, maximum reuse.
+
+Usage:
+  gen_coloring.py --family banded --n 40 --k 4 --w 4 --seed 1 > 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 '<n> <k>' + 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()

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()

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()