diff --git a/src/BooleanInference.jl b/src/BooleanInference.jl index 5a6d45d..66bc7b7 100644 --- a/src/BooleanInference.jl +++ b/src/BooleanInference.jl @@ -24,6 +24,8 @@ include("core/domain.jl") include("core/stats.jl") include("core/problem.jl") +include("preprocessing/canonicalize.jl") + include("utils/utils.jl") include("utils/twosat.jl") include("utils/circuit2cnf.jl") @@ -53,6 +55,7 @@ export Region export is_fixed, has0, has1, init_doms, get_var_value, bits export setup_problem, setup_from_cnf, setup_from_circuit, setup_from_sat +export bounded_ve_canonicalize export factoring_problem, factoring_circuit, factoring_csp export is_solved @@ -62,7 +65,7 @@ export solve_circuit_sat export NumUnfixedVars -export MostOccurrenceSelector +export MostOccurrenceSelector, DiffLookaheadSelector export TNContractionSolver diff --git a/src/branch_table/contraction.jl b/src/branch_table/contraction.jl index a278d50..0f704ee 100644 --- a/src/branch_table/contraction.jl +++ b/src/branch_table/contraction.jl @@ -37,7 +37,7 @@ function slicing(static::ConstraintNetwork, tensor::BoolTensor, doms::Vector{Dom dm = doms[var_id] is_fixed(dm) || push!(free_axes, i) end - fixed_mask, fixed_val = mask_value(doms, tensor.var_axes, UInt16) + fixed_mask, fixed_val = mask_value(doms, tensor.var_axes, UInt32) dims = ntuple(_ -> 2, length(free_axes)) out = fill(ZERO_TROP, dims) # Allocate dense array diff --git a/src/branch_table/selector.jl b/src/branch_table/selector.jl index c4a479c..9d388aa 100644 --- a/src/branch_table/selector.jl +++ b/src/branch_table/selector.jl @@ -61,6 +61,63 @@ function findbest(cache::RegionCache, problem::TNProblem, measure::AbstractMeasu return (OptimalBranchingCore.get_clauses(result), variables) end +# Difficulty-guided lookahead selector (selector-design S1/S2). +# Among the top-`pool` candidate vars (by connection score), probe both polarities +# with GAC and pick the var whose HARDER child has the lowest connectivity-weighted +# difficulty (sum of active tensor degrees). Failed literals are taken immediately. +# Cheap (O(pool) propagations/node) and beats MostOccurrence on hard instances. +struct DiffLookaheadSelector <: AbstractSelector + k::Int + max_tensors::Int + pool::Int +end +DiffLookaheadSelector(k::Int, max_tensors::Int) = DiffLookaheadSelector(k, max_tensors, 16) + +@inline function _sum_active_degree(static::ConstraintNetwork, doms::Vector{DomainMask}) + s = 0 + @inbounds for t in static.tensors + for v in t.var_axes + !is_fixed(doms[v]) && (s += 1) + end + end + return s +end + +function findbest(cache::RegionCache, problem::TNProblem, measure::AbstractMeasure, set_cover_solver::AbstractSetCoverSolver, sel::DiffLookaheadSelector) + scores = compute_var_cover_scores_weighted(problem) + doms = problem.doms + cands = Int[] + @inbounds for i in eachindex(scores) + (is_fixed(doms[i]) || scores[i] <= 0.0) && continue + push!(cands, i) + end + isempty(cands) && return nothing, Int[] + sort!(cands, by = i -> -scores[i]) + length(cands) > sel.pool && (cands = cands[1:sel.pool]) + + buffer = problem.buffer + best = typemax(Int); var_id = 0 + @inbounds for u in cands + c0 = probe_assignment_core!(problem, buffer, doms, [u], UInt64(1), UInt64(0)) + f0 = has_contradiction(c0); d0 = f0 ? 0 : _sum_active_degree(problem.static, c0) + c1 = probe_assignment_core!(problem, buffer, doms, [u], UInt64(1), UInt64(1)) + f1 = has_contradiction(c1); d1 = f1 ? 0 : _sum_active_degree(problem.static, c1) + if f0 || f1 + var_id = u; break # failed literal ⇒ forced, take immediately + end + s = max(d0, d1) + s < best && (best = s; var_id = u) + end + var_id == 0 && (var_id = cands[1]) + + # probing above scribbled in the branching cache via no path that matters, but + # clear it so compute_branching_result starts clean for the chosen var + empty!(buffer.branching_cache) + result, variables = compute_branching_result(cache, problem, var_id, measure, set_cover_solver) + isnothing(result) && return nothing, variables + return (OptimalBranchingCore.get_clauses(result), variables) +end + # struct MinGammaSelector <: AbstractSelector # k::Int # max_tensors::Int diff --git a/src/branching/propagate.jl b/src/branching/propagate.jl index b407c60..019873b 100644 --- a/src/branching/propagate.jl +++ b/src/branching/propagate.jl @@ -1,11 +1,11 @@ -function scan_supports(support::Vector{UInt16}, support_or::UInt16, support_and::UInt16, query_mask0::UInt16, query_mask1::UInt16) - m = query_mask0 | query_mask1 +function scan_supports(support::Vector{UInt32}, support_or::UInt32, support_and::UInt32, query_mask0::UInt32, query_mask1::UInt32) + m = query_mask0 | query_mask1 # General case: filter by compatibility - if m == UInt16(0) + if m == UInt32(0) return support_or, support_and, !isempty(support) end - valid_or_agg = UInt16(0) - valid_and_agg = UInt16(0xFFFF) + valid_or_agg = UInt32(0) + valid_and_agg = UInt32(0xFFFFFFFF) found_any = false @inbounds for i in eachindex(support) config = support[i] @@ -14,7 +14,7 @@ function scan_supports(support::Vector{UInt16}, support_or::UInt16, support_and: valid_and_agg &= config found_any = true # Early exit once both aggregates are saturated. - if valid_or_agg == UInt16(0xFFFF) && valid_and_agg == UInt16(0x0000) + if valid_or_agg == UInt32(0xFFFFFFFF) && valid_and_agg == UInt32(0x00000000) break end end @@ -24,13 +24,13 @@ end # return (query_mask0, query_mask1) function compute_query_masks(doms::Vector{DomainMask}, var_axes::Vector{Int}) - @assert length(var_axes) <= 16 - mask0 = UInt16(0); mask1 = UInt16(0); + @assert length(var_axes) <= 32 + mask0 = UInt32(0); mask1 = UInt32(0); @inbounds for i in eachindex(var_axes) var_id = var_axes[i] domain = doms[var_id] - bit = UInt16(1) << (i - 1) + bit = UInt32(1) << (i - 1) if domain == DM_0 mask0 |= bit elseif domain == DM_1 @@ -52,15 +52,15 @@ struct PropagationContext v2c::Vector{Vector{Int}} end -@inline function apply_updates!(doms::Vector{DomainMask}, var_axes::Vector{Int}, valid_or::UInt16, valid_and::UInt16, ctx::PropagationContext) +@inline function apply_updates!(doms::Vector{DomainMask}, var_axes::Vector{Int}, valid_or::UInt32, valid_and::UInt32, ctx::PropagationContext) @inbounds for i in 1:length(var_axes) var_id = var_axes[i] old_domain = doms[var_id] (old_domain == DM_0 || old_domain == DM_1) && continue - bit = UInt16(1) << (i - 1) - can_be_1 = (valid_or & bit) != UInt16(0) - must_be_1 = (valid_and & bit) != UInt16(0) + bit = UInt32(1) << (i - 1) + can_be_1 = (valid_or & bit) != UInt32(0) + must_be_1 = (valid_and & bit) != UInt32(0) new_dom = must_be_1 ? DM_1 : (can_be_1 ? DM_BOTH : DM_0) diff --git a/src/core/static.jl b/src/core/static.jl index 3e607d4..eb6ac2d 100644 --- a/src/core/static.jl +++ b/src/core/static.jl @@ -39,9 +39,9 @@ end # Shared tensor data (flyweight pattern for deduplication) struct TensorData dense_tensor::BitVector # For contraction operations: satisfied_configs[config+1] = true - support::Vector{UInt16} # For propagation: list of satisfied configs (0-indexed) - support_or::UInt16 # OR over support (for fast m==0 scan) - support_and::UInt16 # AND over support (for fast m==0 scan) + support::Vector{UInt32} # For propagation: list of satisfied configs (0-indexed) + support_or::UInt32 # OR over support (for fast m==0 scan) + support_and::UInt32 # AND over support (for fast m==0 scan) end function Base.show(io::IO, td::TensorData) @@ -51,9 +51,9 @@ end # Extract sparse support from dense BitVector function extract_supports(dense_tensor::BitVector) indices = findall(dense_tensor) - supports = Vector{UInt16}(undef, length(indices)) + supports = Vector{UInt32}(undef, length(indices)) @inbounds for i in eachindex(indices) - supports[i] = UInt16(indices[i] - 1) # 0-indexed + supports[i] = UInt32(indices[i] - 1) # 0-indexed end return supports end @@ -61,8 +61,8 @@ end # Constructor that automatically extracts support function TensorData(dense_tensor::BitVector) support = extract_supports(dense_tensor) - support_or = UInt16(0) - support_and = UInt16(0xFFFF) + support_or = UInt32(0) + support_and = UInt32(0xFFFFFFFF) @inbounds for i in eachindex(support) config = support[i] support_or |= config diff --git a/src/preprocessing/canonicalize.jl b/src/preprocessing/canonicalize.jl new file mode 100644 index 0000000..901919c --- /dev/null +++ b/src/preprocessing/canonicalize.jl @@ -0,0 +1,185 @@ +# ============================================================================ +# Static, width-aware constraint-network canonicalizer (bounded-width VE). +# +# Generalizes `precontract_degree2!` (src/core/static.jl) from "eliminate any +# degree-2 variable" to "eliminate any variable whose elimination stays under a +# space-complexity budget B, in weighted-min-fill order". Degree-2-greedy is not +# width-aware and P6 showed it *raises* treewidth; min-fill ordering with a per- +# step width cap keeps it bounded. +# ============================================================================ + +""" + bounded_ve_canonicalize(cn::ConstraintNetwork; budget_B::Real, + protected=Int[], + order::Symbol=:weighted_min_fill) -> ConstraintNetwork + +Reshape `cn` by bucket-eliminating variables once, statically, before any search. +Eliminating a variable `v` joins all tensors incident to `v` and projects `v` out +(boolean ∃/∧ semantics, exactly as `contract_two_tensors`). A variable is eliminated +only if the join's space complexity (`sc`, log2 of the largest intermediate under +`GreedyMethod`) is `<= budget_B`; eligible variables are removed in weighted-min-fill +order (fewest new neighbor edges first, `sc` as tiebreaker). `budget_B` is the single +fine↔coarse knob. + +`protected` lists variable ids (in `cn`'s own id space) that must NEVER be eliminated — +the read-out variables (e.g. factor bits). They survive into the result, so their values can +be read directly off the solved reduced network via `result.orig_to_new[orig_id]` — no VE +back-substitution needed. Eliminated (non-protected) variables' values are *not* recoverable, +which is fine when only the protected variables are read. + +Produced tensors are additionally hard-capped at 32 variables, the limit of the `TensorData` +representation (UInt32 config indices / 32-bit support masks), regardless of `budget_B`. + +Returns a new compressed `ConstraintNetwork` whose remaining variables form the branch set +(protected variables ⊆ branch set). +""" +function bounded_ve_canonicalize(cn::ConstraintNetwork; budget_B::Real, + protected=Int[], + order::Symbol=:weighted_min_fill) + order == :weighted_min_fill || throw(ArgumentError("unsupported order $order")) + protected_set = Set{Int}(protected) + B = Float64(budget_B) + nv = length(cn.vars) + # TensorData stores configs as UInt32 (support indices + 32-bit support_or/and masks), + # so a produced tensor can hold at most 32 variables. Hard structural cap, regardless of B. + # In practice budget_B binds first (a dense 2^arity tensor is infeasible well before 32). + MAX_ARITY = 32 + + # Mutable working copies in cn-variable-id space. Deep-copy var_axes so we never + # mutate the input network (compress_variables! rewrites var_axes in place). + tensors = [BoolTensor(copy(t.var_axes), t.tensor_data_idx) for t in cn.tensors] + vars_to_tensors = [copy(lst) for lst in cn.v2t] + unique_data = copy(cn.unique_tensors) + data_to_idx = Dict{BitVector,Int}() + for (i, td) in enumerate(unique_data) + get!(data_to_idx, td.dense_tensor, i) + end + active = trues(length(tensors)) + + active_incident(v) = filter(t -> active[t], vars_to_tensors[v]) + + # Variables of the join of v's incident tensors, minus v (the produced tensor axes). + function out_vars(tids, v) + out = Int[] + @inbounds for t in tids, x in tensors[t].var_axes + x != v && !(x in out) && push!(out, x) + end + return out + end + + # Space complexity of eliminating v under the current state (no array execution). + function elim_sc(tids, out) + code = EinCode([copy(tensors[t].var_axes) for t in tids], out) + sd = uniformsize(code, 2) + return contraction_complexity(optimize_code(code, sd, GreedyMethod()), sd).sc + end + + # Weighted-min-fill: number of neighbor pairs not already sharing an active tensor. + function fill_count(out) + f = 0 + @inbounds for i in 1:length(out)-1, j in i+1:length(out) + a, b = out[i], out[j] + share = false + for t in vars_to_tensors[a] + if active[t] && b in tensors[t].var_axes + share = true; break + end + end + share || (f += 1) + end + return f + end + + # (eligible, fill, sc) for variable v in the current state. + function score(v) + v in protected_set && return (false, 0, Inf) # read-out var: never eliminate + tids = active_incident(v) + isempty(tids) && return (false, 0, Inf) # isolated: nothing to eliminate + out = out_vars(tids, v) + length(out) > MAX_ARITY && return (false, 0, Inf) # too wide for TensorData + sc = elim_sc(tids, out) + return (sc <= B, fill_count(out), sc) + end + + pq = PriorityQueue{Int,Tuple{Int,Float64}}() + for v in 1:nv + elig, f, sc = score(v) + elig && (pq[v] = (f, sc)) + end + + while !isempty(pq) + v = dequeue!(pq) + tids = active_incident(v) + isempty(tids) && continue + out = out_vars(tids, v) + length(out) > MAX_ARITY && continue # too wide for TensorData + code = EinCode([copy(tensors[t].var_axes) for t in tids], out) + sd = uniformsize(code, 2) + optcode = optimize_code(code, sd, GreedyMethod()) + contraction_complexity(optcode, sd).sc <= B || continue # stale: no longer eligible + + # Execute the bucket contraction -> dense BitVector over `out`. + arrs = [reshape(Int.(unique_data[tensors[t].tensor_data_idx].dense_tensor), + ntuple(_ -> 2, length(tensors[t].var_axes))) for t in tids] + res = optcode(arrs...) + gt = res .> 0 + new_data = gt isa AbstractArray ? BitVector(vec(gt)) : BitVector([gt]) + + # Deduplicate produced tensor data (flyweight, as in setup_problem). + idx = get(data_to_idx, new_data, 0) + if idx == 0 + push!(unique_data, TensorData(new_data)) + idx = length(unique_data) + data_to_idx[new_data] = idx + end + + # Rewrite incidence: reuse the first slot for the merged tensor, drop the rest. + keep = tids[1] + for t in tids + for x in tensors[t].var_axes + filter!(tt -> tt != t, vars_to_tensors[x]) + end + end + for t in tids[2:end] + active[t] = false + end + tensors[keep] = BoolTensor(copy(out), idx) + for x in out + push!(vars_to_tensors[x], keep) + end + # v now has no active incident tensor -> eliminated. + + # Re-score only the affected neighbors (the merged tensor's variables). + for u in out + elig, f, scu = score(u) + if elig + pq[u] = (f, scu) + elseif haskey(pq, u) + delete!(pq, u) + end + end + end + + # Compact active tensors, rebuild incidence, then compress variable ids. + active_idx = findall(active) + new_tensors = tensors[active_idx] + new_vars_to_tensors = [Int[] for _ in 1:nv] + for (newt, _) in enumerate(active_idx) + for x in new_tensors[newt].var_axes + push!(new_vars_to_tensors[x], newt) + end + end + + new_tensors, compressed_v2t, compress_o2n = + compress_variables!(new_tensors, new_vars_to_tensors) + vars = [Variable(length(compressed_v2t[i])) for i in 1:length(compressed_v2t)] + + # Compose with cn's own orig->compressed map so the result indexes original var ids. + orig_to_new = zeros(Int, length(cn.orig_to_new)) + for orig in 1:length(cn.orig_to_new) + cnid = cn.orig_to_new[orig] + orig_to_new[orig] = cnid == 0 ? 0 : compress_o2n[cnid] + end + + return ConstraintNetwork(vars, unique_data, new_tensors, compressed_v2t, orig_to_new) +end diff --git a/test/canonicalize.jl b/test/canonicalize.jl new file mode 100644 index 0000000..533982e --- /dev/null +++ b/test/canonicalize.jl @@ -0,0 +1,48 @@ +using Test +using BooleanInference +using BooleanInference: setup_from_sat, TNProblem, solve, BranchingStrategy, TNContractionSolver, + MostOccurrenceSelector, NumUnfixedVars, NoReducer, get_var_value, bits_to_int, + bounded_ve_canonicalize, ConstraintNetwork, count_unfixed +using OptimalBranchingCore: GreedyMerge +import ProblemReductions: CircuitSAT, Factoring, reduceto + +# Minimal coverage for the bounded-VE canonicalizer (M1) wired into the solver with protected +# read-out variables (M3): canonicalizing must preserve satisfiability, keep protected (factor-bit) +# variables alive, shrink the branch set, and let the factors be read straight off the solution. +@testset "bounded_ve_canonicalize + protected read-out" begin + N = 31 * 29 + red = reduceto(CircuitSAT, Factoring(5, 5, N)) + sat = CircuitSAT(red.circuit.circuit; use_constraints=true) + base = setup_from_sat(sat) + o2n = base.static.orig_to_new + q_orig = collect(red.q); p_orig = collect(red.p) + + # protected = factor-bit variables, in base-network id space (never eliminated) + prot = Int[o2n[v] for v in vcat(q_orig, p_orig) if o2n[v] != 0] + @test !isempty(prot) + + cn = bounded_ve_canonicalize(base.static; budget_B=6, protected=prot) + @test cn isa ConstraintNetwork + + # M1: elimination actually happened — the branch set is strictly smaller. + @test length(cn.vars) < length(base.static.vars) + + # M3: every protected factor-bit variable survives into the reduced network. + for v in vcat(q_orig, p_orig) + o2n[v] == 0 && continue + @test cn.orig_to_new[v] != 0 + end + + # Solve the reduced network (cadical-free path) and read factors directly off the solution. + strat = BranchingStrategy(table_solver=TNContractionSolver(), + selector=MostOccurrenceSelector(1, 2), measure=NumUnfixedVars(), + set_cover_solver=GreedyMerge()) + res = solve(TNProblem(cn), strat, NoReducer()) + @test res.found + @test count_unfixed(res.solution) == 0 + + qb = Bool[get_var_value(res.solution, cn.orig_to_new[v]) == 1 for v in q_orig] + pb = Bool[get_var_value(res.solution, cn.orig_to_new[v]) == 1 for v in p_orig] + a, b = bits_to_int(qb), bits_to_int(pb) + @test a * b == N # factors read with NO variable-elimination back-substitution +end diff --git a/test/runtests.jl b/test/runtests.jl index 912de89..62ab4ac 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -5,6 +5,10 @@ using Test include("2sat.jl") end +@testset "canonicalize.jl" begin + include("canonicalize.jl") +end + @testset "problems.jl" begin include("problems.jl") end