diff --git a/.codecov.yml b/.codecov.yml new file mode 100644 index 0000000..f6d756e --- /dev/null +++ b/.codecov.yml @@ -0,0 +1,10 @@ +coverage: + status: + patch: + default: + target: 70% + threshold: 2% + project: + default: + target: auto + threshold: 10% \ No newline at end of file diff --git a/.gitignore b/.gitignore index fbce418..3dad150 100644 --- a/.gitignore +++ b/.gitignore @@ -11,6 +11,7 @@ benchmarks/Manifest.toml # IDE .vscode/ .idea/ +CLAUDE.md # OS .DS_Store @@ -19,6 +20,22 @@ Thumbs.db # Temporary files *.tmp *.temp +*.swp +*.swo statprof/ -OptimalBranching.jl/ \ No newline at end of file +OptimalBranching.jl/ +.julia/ +.claude/ + +# === Build outputs === +/benchmarks/artifacts/ +/benchmarks/third-party/**/build/ +/benchmarks/third-party/**/bin/ +/benchmarks/third-party/**/CMakeFiles/ +/benchmarks/third-party/**/Makefile +/benchmarks/third-party/**/cmake_install.cmake + +# === Logs === +*.log +*.out \ No newline at end of file diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 0000000..e61a015 --- /dev/null +++ b/.gitmodules @@ -0,0 +1,16 @@ +[submodule "benchmarks/third-party/x-sat"] + path = benchmarks/third-party/x-sat + url = https://github.com/shaowei-cai-group/X-SAT.git + ignore = all +[submodule "benchmarks/third-party/aiger"] + path = benchmarks/third-party/aiger + url = https://github.com/arminbiere/aiger.git + ignore = all +[submodule "benchmarks/third-party/abc"] + path = benchmarks/third-party/abc + url = https://github.com/berkeley-abc/abc.git + ignore = all +[submodule "benchmarks/third-party/CnC"] + path = benchmarks/third-party/CnC + url = https://github.com/marijnheule/CnC.git + ignore = all \ No newline at end of file diff --git a/Project.toml b/Project.toml index 7eebf4d..a1ef5bb 100644 --- a/Project.toml +++ b/Project.toml @@ -9,6 +9,7 @@ GenericTensorNetworks = "3521c873-ad32-4bb4-b63d-f4f178f42b49" OptimalBranchingCore = "c76e7b22-e1d2-40e8-b0f1-f659837787b8" ProblemReductions = "899c297d-f7d2-4ebf-8815-a35996def416" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" TropicalNumbers = "b3a74e9c-7526-4576-a4eb-79c0d4c32334" [compat] @@ -17,6 +18,7 @@ GenericTensorNetworks = "4" OptimalBranchingCore = "0.1" ProblemReductions = "0.3.5" SparseArrays = "1.11.0" +Statistics = "1.11.1" TropicalNumbers = "0.6.2" julia = "1" diff --git a/benchmarks/Project.toml b/benchmarks/Project.toml index 7f0e4d5..15665c7 100644 --- a/benchmarks/Project.toml +++ b/benchmarks/Project.toml @@ -8,6 +8,7 @@ BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" BooleanInference = "d12a6243-8379-43cb-9459-098119da1519" Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" Gurobi = "2e9cd046-0924-5485-92f1-d5272153d98b" +HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" JSON3 = "0f8b85d8-7281-11e9-16c2-39a750bddbf1" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" Primes = "27ebfcd6-29c5-5fa9-bf4b-fb8fc14df3ae" @@ -20,6 +21,7 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" BenchmarkTools = "1.6.0" Dates = "1.11.0" Gurobi = "1.7.6" +HiGHS = "1.19.3" JSON3 = "1.14.3" JuMP = "1.29.1" ProblemReductions = "0.3.5" diff --git a/benchmarks/examples/factoring_example.jl b/benchmarks/examples/factoring_example.jl new file mode 100644 index 0000000..400ecfa --- /dev/null +++ b/benchmarks/examples/factoring_example.jl @@ -0,0 +1,62 @@ +using Pkg +Pkg.activate(joinpath(@__DIR__, "..")) + +using BooleanInferenceBenchmarks +using Random +using Gurobi + +const env = Gurobi.Env() + +# ---------------------------------------- +# Example 1: Generate datasets +# ---------------------------------------- + +println("=" ^80) +println("Generating Factoring Datasets") +println("=" ^80) + +configs = [ + # FactoringConfig(10, 10), + # FactoringConfig(12, 12), + # FactoringConfig(14, 14), + # FactoringConfig(16, 16), + # FactoringConfig(18, 18), + # FactoringConfig(20, 20), + FactoringConfig(22, 22), + # FactoringConfig(24, 24) +] + +paths = generate_factoring_datasets(configs; per_config=5, include_solution=true, force_regenerate=false) + +# ---------------------------------------- +# Example 2: Run benchmark on a single dataset +# ---------------------------------------- + +# println("\n" * "=" ^80) +# println("Running Benchmark") +# println("=" ^80) + +# # result = benchmark_dataset(FactoringProblem, paths[1]; solver=XSATSolver(;yosys_path="/opt/homebrew/bin/yosys")) +# result = benchmark_dataset(FactoringProblem, paths[1]; solver=IPSolver(Gurobi.Optimizer, env)) + +# if !isnothing(result) +# println("\nResults:") +# println(" Dataset: $(result["dataset_path"])") +# println(" Instances tested: $(result["instances_tested"])") +# println(" Accuracy: $(round(result["accuracy_rate"]*100, digits=1))%") +# println(" Median time: $(result["median_time"]) seconds") +# end + +# ---------------------------------------- +# Example 3: Compare solvers +# ---------------------------------------- + +println("\n" * "=" ^80) +println("Comparing Solvers") +println("=" ^80) + +# results = run_solver_comparison(FactoringProblem, paths, solvers=[XSATSolver(;yosys_path="/opt/homebrew/bin/yosys"), BooleanInferenceSolver()]) +# results = run_solver_comparison(FactoringProblem, paths, solvers=[XSATSolver(;yosys_path="/opt/homebrew/bin/yosys")]) +results = run_solver_comparison(FactoringProblem, paths, solvers=[BooleanInferenceSolver()]) +print_solver_comparison_summary(results) + diff --git a/benchmarks/examples/run_benchmarks.jl b/benchmarks/examples/run_benchmarks.jl deleted file mode 100755 index 33b118d..0000000 --- a/benchmarks/examples/run_benchmarks.jl +++ /dev/null @@ -1,70 +0,0 @@ -# Comprehensive solver comparison benchmark script -# Demonstrates both automatic and manual solver comparison approaches -# Usage: julia --project=benchmark benchmark/examples/run_benchmarks.jl - -using Pkg -Pkg.activate(joinpath(@__DIR__, "..")) - -using BooleanInferenceBenchmarks - -using Gurobi -const GUROBI_ENV = Gurobi.Env() - -function main() - println("Running BooleanInference Solver Comparison Benchmarks...") - - # Show available solvers - println("\nAvailable solvers:") - list_available_solvers(FactoringProblem) - - # Configure benchmark parameters - configs = [(10,10), (12,12), (14,14)] - dataset_per_config = 10 - - println("\nComparing BooleanInference vs IP Solver...") - println("Configurations: $configs") - println("Instances per config: $dataset_per_config") - - # Option 1: Automatic comparison of all available solvers - # println("\n=== Automatic Solver Comparison (All Available Solvers) ===") - all_results = run_solver_comparison(FactoringProblem, configs; - dataset_per_config=dataset_per_config) - - # Print beautiful comparison table - print_solver_comparison_summary(all_results) - # Persist comparison results - save_solver_comparison(FactoringProblem, all_results; note="demo run") - - # # Option 2: Manual head-to-head comparison for detailed analysis - # println("\n=== Head-to-Head Comparison (Detailed Analysis) ===") - - # # Test BooleanInference solver - # println("Benchmarking BooleanInference solver...") - # bi_results = run_full_benchmark(FactoringProblem, configs; - # dataset_per_config=dataset_per_config, - # solver=BooleanInferenceSolver()) - - # # Test IP solver - # println("Benchmarking IP solver...") - # ip_results = run_full_benchmark(FactoringProblem, configs; - # dataset_per_config=dataset_per_config, - # solver=IPSolver(Gurobi.Optimizer, GUROBI_ENV)) - - # # Side-by-side comparison - # println("\nHead-to-Head Results:") - # compare_solver_results("BooleanInference", bi_results, "IP-Gurobi", ip_results) - - # println("\nAll benchmarks completed!") - # println("Datasets cached in: benchmarks/data/factoring/") - # println("Results saved in: benchmarks/data/factoring_benchmark_results.json") - - # println("\nFor quick one-liner comparisons, you can also use:") - # println(" # Compare all solvers automatically:") - # println(" run_solver_comparison(FactoringProblem, [(10,10)])") - # println(" # Or test a specific solver:") - # println(" run_full_benchmark(FactoringProblem, [(10,10)]; solver=BooleanInferenceSolver())") -end - -if abspath(PROGRAM_FILE) == @__FILE__ - main() -end diff --git a/benchmarks/src/BooleanInferenceBenchmarks.jl b/benchmarks/src/BooleanInferenceBenchmarks.jl index dff693f..18ce260 100644 --- a/benchmarks/src/BooleanInferenceBenchmarks.jl +++ b/benchmarks/src/BooleanInferenceBenchmarks.jl @@ -2,36 +2,42 @@ module BooleanInferenceBenchmarks using Random using JSON3 +using Primes +using Dates: now using BenchmarkTools using BooleanInference -using JuMP +using JuMP, HiGHS, Gurobi +using Statistics: mean, median +using SHA: bytes2hex, sha256, sha1 using ProblemReductions -using Statistics: mean, median, std, quantile -using Dates: now -using SHA: bytes2hex, sha256 - - -# Core utilities -include("utils.jl") +using ProblemReductions: Factoring, CircuitSAT, reduceto, constraints, objectives, AbstractProblem, Symbol -# Abstract types and interfaces include("abstract_types.jl") - -# Generic benchmark framework -include("generic_benchmark.jl") - -# Problem implementations -include("factoring/factoring_problem.jl") -include("factoring/factor_IP.jl") - -# Re-export main interfaces -export AbstractBenchmarkProblem, AbstractProblemConfig, AbstractSolver -export generate_instance, solve_instance, problem_id, default_configs, filename_pattern +include("utils.jl") +include("benchmark.jl") +include("formatting.jl") +include("comparison.jl") + +# CircuitIO +include("circuitIO/circuitIO.jl") +using .CircuitIO + +# Factoring problem +include("factoring/types.jl") +include("factoring/interface.jl") +include("factoring/generators.jl") +include("factoring/solver_ip.jl") +include("factoring/solvers.jl") +include("factoring/dataset.jl") +include("factoring/solver_xsat.jl") + +export AbstractBenchmarkProblem, AbstractProblemConfig, AbstractInstance, AbstractSolver +export solve_instance, verify_solution, read_instances export available_solvers, default_solver, solver_name -export generate_datasets, benchmark_problem, run_full_benchmark, benchmark_backend -export create_configs_from_tuples -export list_available_solvers, run_solver_comparison -export FactoringProblem, FactoringConfig, BooleanInferenceSolver, IPSolver -export print_solver_comparison_summary, save_solver_comparison +export benchmark_dataset, run_solver_comparison +export list_available_solvers, print_solver_comparison_summary, compare_solver_results +export FactoringProblem, FactoringConfig, FactoringInstance, generate_factoring_datasets +export CircuitSATProblem, CircuitSATConfig, CircuitSATInstance +export BooleanInferenceSolver, IPSolver, XSATSolver -end # module \ No newline at end of file +end \ No newline at end of file diff --git a/benchmarks/src/abstract_types.jl b/benchmarks/src/abstract_types.jl index 4dc94fc..31b90d7 100644 --- a/benchmarks/src/abstract_types.jl +++ b/benchmarks/src/abstract_types.jl @@ -1,98 +1,85 @@ -# Abstract types and interfaces for benchmark problems - -""" -Abstract base type for all benchmark problems. -""" abstract type AbstractBenchmarkProblem end - -""" -Abstract base type for problem configurations. -""" abstract type AbstractProblemConfig end - -""" -Abstract base type for solvers. -""" +abstract type AbstractInstance end abstract type AbstractSolver end -# Core interface methods that each problem type must implement +@kwdef struct BooleanInferenceSolver <: AbstractSolver + warmup::Bool = true +end -""" - generate_instance(problem_type::Type{<:AbstractBenchmarkProblem}, config::AbstractProblemConfig; rng, include_solution=false) +struct IPSolver <: AbstractSolver + warmup::Bool + optimizer::Any + env::Any + function IPSolver(optimizer=Gurobi.Optimizer, env=nothing) + new(true, optimizer, env) + end +end -Generate a single instance of the given problem type. -""" -function generate_instance end +struct XSATSolver <: AbstractSolver + warmup::Bool + csat_path::String + yosys_path::String + satisfiable_when_high::Bool + function XSATSolver(;csat_path=joinpath(dirname(@__DIR__), "artifacts", "bin", "csat"), yosys_path=nothing) + !isfile(csat_path) && error("File $csat_path for X-SAT solver does not exist") + # check if yosys is installed by homebrew or apt + + if isnothing(yosys_path) + yosys_path = try + yosys_path = strip(read(`which yosys`, String)) + catch + error("Yosys not found in PATH, and yosys_path is not provided") + end + else + !isfile(yosys_path) && error("File $yosys_path for Yosys does not exist") + end + new(false, csat_path, yosys_path, true) + end +end -""" - solve_instance(problem_type::Type{<:AbstractBenchmarkProblem}, instance, solver::AbstractSolver) +struct KissatSolver <: AbstractSolver + warmup::Bool + kissat_path::String + yosys_path::String + abc_path::String + satisfiable_when_high::Bool + function KissatSolver(;kissat_path=nothing, yosys_path=nothing, abc_path=nothing) + if isnothing(kissat_path) + kissat_path = try + kissat_path = strip(read(`which kissat`, String)) + catch + error("Kissat not found in PATH, and kissat_path is not provided") + end + else + !isfile(kissat_path) && error("File $kissat_path for Kissat does not exist") + end + new(true, kissat_path, yosys_path, abc_path, false) + end +end -Solve a single problem instance using the specified solver. -""" +# ---------------------------------------- +# Core Interface (must be implemented) +# ---------------------------------------- function solve_instance end +function verify_solution end +function read_instances end # Read instances from a dataset file/directory +function available_solvers end +function default_solver end -""" - solve_instance(problem_type::Type{<:AbstractBenchmarkProblem}, instance) - -Solve a single problem instance using the default solver. -""" +# Default fallback for solve_instance function solve_instance(problem_type::Type{<:AbstractBenchmarkProblem}, instance) solve_instance(problem_type, instance, default_solver(problem_type)) end -""" - problem_id(instance) - -Generate a stable ID for the instance. -""" -function problem_id end - -""" - default_configs(problem_type::Type{<:AbstractBenchmarkProblem}) - -Return default configurations for benchmarking this problem type. -""" -function default_configs end - -""" - config(problem_type::Type{<:AbstractBenchmarkProblem}, params) - -Generate a configuration for the given problem type and parameters. -""" -function config end - -""" - filename_pattern(problem_type::Type{<:AbstractBenchmarkProblem}, config::AbstractProblemConfig) - -Generate filename for dataset files. -""" -function filename_pattern end - -""" - available_solvers(problem_type::Type{<:AbstractBenchmarkProblem}) - -Return available solvers for this problem type. -""" -function available_solvers end - -""" - default_solver(problem_type::Type{<:AbstractBenchmarkProblem}) - -Return the default solver for this problem type. -""" -function default_solver end - -""" - verify_solution(problem_type::Type{<:AbstractBenchmarkProblem}, instance, result) - -Verify that the solver result is correct for the given instance. -Returns true if correct, false otherwise. -""" -function verify_solution end +function solver_name(::BooleanInferenceSolver) + return "BI" +end -""" - solver_name(solver::AbstractSolver) +function solver_name(solver::IPSolver) + return "IP-$(string(solver.optimizer)[1:end-10])" +end -Return human-readable name for the solver. -""" -function solver_name end +function solver_name(::XSATSolver) + return "X-SAT" +end \ No newline at end of file diff --git a/benchmarks/src/benchmark.jl b/benchmarks/src/benchmark.jl new file mode 100644 index 0000000..995dfec --- /dev/null +++ b/benchmarks/src/benchmark.jl @@ -0,0 +1,113 @@ +""" + benchmark_dataset(problem_type, dataset_path; solver=nothing, warmup=true) + +Run benchmark on a single dataset file/directory. + +# Arguments +- `problem_type::Type{<:AbstractBenchmarkProblem}`: The type of problem to benchmark +- `dataset_path::String`: Path to the dataset file or directory +- `solver`: The solver to use (defaults to the problem type's default solver) +- `warmup::Bool`: Whether to perform warmup/tuning (default: true). Set to false for solvers + implemented in non-Julia languages (e.g., C++), which don't benefit from Julia's JIT compilation. +""" +function benchmark_dataset(problem_type::Type{<:AbstractBenchmarkProblem}, dataset_path::String; solver=nothing) + actual_solver = isnothing(solver) ? default_solver(problem_type) : solver + warmup = actual_solver.warmup + solver_info = solver_name(actual_solver) + + @info "Benchmarking: $dataset_path" + @info "Using solver: $solver_info" + + try + if !isfile(dataset_path) && !isdir(dataset_path) + @error "Dataset not found: $dataset_path" + return nothing + end + + @info " Loading instances from: $dataset_path" + instances = read_instances(problem_type, dataset_path) + @info " Testing $(length(instances)) instances" + + if isempty(instances) + @error " No instances found in dataset" + return nothing + end + + all_times = Float64[] + all_memory = Int64[] + all_gc_times = Float64[] + successful_runs = 0 + failed_runs = 0 + + # Warmup and tune benchmark parameters if requested + b = nothing + if warmup + @info " Performing warmup and tuning..." + sample_trial_fn = () -> solve_instance(problem_type, instances[1], solver) + b = @benchmarkable $sample_trial_fn() + tune!(b) + @info " Benchmark tuned, starting verification and timing..." + else + @info " Skipping warmup (non-Julia solver), starting verification and timing..." + end + + correct_runs = 0 + incorrect_runs = 0 + + for (i, instance) in enumerate(instances) + try + result = solve_instance(problem_type, instance, solver) + is_correct = verify_solution(problem_type, instance, result) + + if !is_correct + @warn " Instance $i: Incorrect solution" + incorrect_runs += 1 + continue + end + + correct_runs += 1 + + trial_fn = () -> solve_instance(problem_type, instance, solver) + b_instance = @benchmarkable $trial_fn() + if warmup && !isnothing(b) + b_instance.params = b.params + end + timing_result = run(b_instance, samples=1) + + push!(all_times, median(timing_result).time / 1e9) + push!(all_memory, median(timing_result).memory) + push!(all_gc_times, median(timing_result).gctime / 1e9) + successful_runs += 1 + + if i % 10 == 0 || i == length(instances) + @info " Completed $i/$(length(instances)) instances ($(correct_runs) correct, $(incorrect_runs) incorrect)" + end + + catch e + @warn " Instance $i failed with error: $e" + failed_runs += 1 + end + end + + if isempty(all_times) + @warn "All instances failed" + return nothing + end + + return Dict( + "dataset_path" => dataset_path, + "instances_tested" => length(instances), + "correct_runs" => correct_runs, + "successful_runs" => successful_runs, + "accuracy_rate" => correct_runs / length(instances), + "median_time" => median(all_times), + "median_memory" => median(all_memory), + "solver" => solver_info + ) + + catch e + @error " Benchmark failed: $e" + return nothing + end +end + diff --git a/benchmarks/src/circuitIO/aig.jl b/benchmarks/src/circuitIO/aig.jl new file mode 100644 index 0000000..39f5677 --- /dev/null +++ b/benchmarks/src/circuitIO/aig.jl @@ -0,0 +1,133 @@ +# ----------------------------- +# Minimal AAG (ASCII AIGER) reader +# ----------------------------- +# Supports: header "aag M I L O A", I inputs, L latches (L must be 0 here), +# O outputs (as literals), A AND gates (lhs rhs0 rhs1). +# Literals: lit = 2*var + is_neg; even = positive; odd = negation. +# +# We deliberately keep a tiny parser here for robustness & no deps. + +struct AIG + inputs::Vector{Int} # even literals + latches::Vector{NTuple{3,Int}} # (q, d, init) -- we disallow L>0 below + outputs::Vector{Int} # any literals + ands::Vector{NTuple{3,Int}} # (lhs, rhs0, rhs1), lhs must be even +end + +# Helpers for literals +varid(l::Int) = l >>> 1 +isneg(l::Int) = (l & 1) == 1 + +function read_aag(path::AbstractString)::AIG + open(path, "r") do io + # header + header = strip(readline(io)) + startswith(header, "aig ") && error("Got binary 'aig' header. Convert to ASCII: `aigtoaig -a in.aig > out.aag`.") + @assert startswith(header, "aag ") "Unrecognized header (need 'aag ...')" + parts = split(header) + @assert length(parts) >= 6 "Header must be: aag M I L O A" + M = parse(Int, parts[2]); I = parse(Int, parts[3]); L = parse(Int, parts[4]) + O = parse(Int, parts[5]); A = parse(Int, parts[6]) + + inputs = [parse(Int, strip(readline(io))) for _ in 1:I] + latches = NTuple{3,Int}[] + for _ in 1:L + toks = split(strip(readline(io))) + q = parse(Int, toks[1]); d = parse(Int, toks[2]) + init = length(toks) >= 3 ? parse(Int, toks[3]) : 2 + push!(latches, (q, d, init)) + end + outputs = [parse(Int, strip(readline(io))) for _ in 1:O] + ands = NTuple{3,Int}[] + for _ in 1:A + toks = split(strip(readline(io))) + lhs = parse(Int, toks[1]); r0 = parse(Int, toks[2]); r1 = parse(Int, toks[3]) + push!(ands, (lhs, r0, r1)) + end + + # Basic sanity + for l in inputs + @assert iseven(l) "Input literal must be even, got $l" + end + for (lhs, _, _) in ands + @assert iseven(lhs) "AND lhs must be even, got $lhs" + end + + return AIG(inputs, latches, outputs, ands) + end +end + +# ----------------------------- +# Name mapping & literal → expression-string +# ----------------------------- +# We will generate a @circuit block as a string then `Meta.parse` + `eval`. +# This avoids touching internal Assignment types. + +""" +Turn an AIG var id to a Julia symbol name used in @circuit. +We prefix with `n_` to avoid accidental name clashes. +""" +varname(v::Int) = Symbol("n_", string(v)) + +""" +Literal -> expression snippet string for @circuit: +- 0 -> false +- 1 -> true +- even -> n_ +- odd -> ¬(n_) +""" +function lit_str(l::Int) + if l == 0 + return "false" + elseif l == 1 + return "true" + else + v = varid(l) + nm = string(varname(v)) + return isneg(l) ? "¬($nm)" : nm + end +end + +# ----------------------------- +# AIG → ProblemReductions.Circuit (combinational) +# ----------------------------- +""" +Build a ProblemReductions.Circuit from an AIG (no latches allowed). +- Each AND gate (lhs, r0, r1) becomes: n_ = ∧(expr(r0), expr(r1)) +- Each primary output becomes: out = expr(lit) +Inputs do not need explicit assignments; they are just symbols used by gates. +""" +function aig_to_circuit(aig::AIG; output_prefix::AbstractString="out") + if !isempty(aig.latches) + error("This converter handles *combinational* circuits only (L=0). Found L=$(length(aig.latches)).") + end + + # Collect all lhs var ids to ensure their symbols are bound in circuit order + # AIGER usually guarantees ANDs are topologically ordered by growing lhs. + assigns = String[] + + # AND gates + for (lhs, r0, r1) in aig.ands + lhs_nm = string(varname(varid(lhs))) + push!(assigns, "$lhs_nm = ∧($(lit_str(r0)), $(lit_str(r1)))") + end + + # Outputs: give them stable names + for (i, l) in enumerate(aig.outputs) + outnm = Symbol(output_prefix, "_", i) + push!(assigns, string(outnm, " = ", lit_str(l))) + end + + src = IOBuffer() + println(src, "ProblemReductions.@circuit begin") + for a in assigns + println(src, " ", a) + end + println(src, "end") + code_str = String(take!(src)) + + # Evaluate macro call to get a `Circuit` + # This is the public, documented entry: @circuit returns a Circuit. [oai_citation:1‡giggleliu.github.io](https://giggleliu.github.io/ProblemReductions.jl/dev/rules/spinglass_sat/) + return Base.invokelatest(eval, Meta.parse(code_str)) +end + diff --git a/benchmarks/src/circuitIO/circuitIO.jl b/benchmarks/src/circuitIO/circuitIO.jl new file mode 100644 index 0000000..56d335e --- /dev/null +++ b/benchmarks/src/circuitIO/circuitIO.jl @@ -0,0 +1,12 @@ +module CircuitIO + using ProblemReductions + using ProblemReductions: Circuit, BooleanExpr, extract_symbols!, simple_form + + include("aig.jl") + include("verilog.jl") + + export aig_to_circuit + + export write_verilog, circuit_to_verilog + +end \ No newline at end of file diff --git a/benchmarks/src/circuitIO/verilog.jl b/benchmarks/src/circuitIO/verilog.jl new file mode 100644 index 0000000..36e387f --- /dev/null +++ b/benchmarks/src/circuitIO/verilog.jl @@ -0,0 +1,292 @@ +# ---------- Verilog codegen for Circuit ---------- + +# Make a Verilog-safe identifier from a Symbol (e.g., Symbol("##var#236") -> "__var_236") + +sanitize_name(s::Symbol) = let raw = String(s) + # replace non-alnum with underscore + cleaned = replace(raw, r"[^A-Za-z0-9_]" => "_") + # if starts with digit, prefix underscore + startswith(cleaned, r"[0-9]") ? "_" * cleaned : cleaned +end + +# Natural sort key: split alpha prefix and numeric suffix so p2 < p10 +function _nat_key(s::Symbol) + str = String(s) + m = match(r"^([A-Za-z_]+)(\d+)$", str) + if m === nothing + return (str, -1) # names without numeric suffix come first among same prefix + else + return (m.captures[1], parse(Int, m.captures[2])) + end +end + +is_true_sym(s::Symbol) = s === Symbol("true") +is_false_sym(s::Symbol) = s === Symbol("false") + +function verilog_expr(ex::BooleanExpr, rename::Dict{Symbol,String}) + if ex.head == :var + s = ex.var + return is_true_sym(s) ? "1'b1" : + is_false_sym(s) ? "1'b0" : + rename[s] # Regular variable + end + + args = verilog_expr.(ex.args, Ref(rename)) + + head = ex.head + if head == :¬ + @assert length(args) == 1 "Unary NOT expects 1 argument" + return "(~" * args[1] * ")" + elseif head == :∧ + return "(" * join(args, " & ") * ")" + elseif head == :∨ + return "(" * join(args, " | ") * ")" + elseif head == :⊻ + return "(" * join(args, " ^ ") * ")" + else + error("Unsupported gate head: $head") + end +end + +# Gather symbol usage: defs (LHS outputs), uses (RHS variables), and all symbols +function collect_def_use(c::Circuit) + defs = Symbol[] + uses = Symbol[] + for ex in c.exprs + append!(defs, ex.outputs) + extract_symbols!(ex.expr, uses) # RHS variables (includes true/false which we'll drop later) + end + # drop boolean constants from uses + filter!(s -> !(is_true_sym(s) || is_false_sym(s)), uses) + return unique!(defs), unique!(uses) +end + +# Infer module I/O and internal wires. +# - inputs: used but never assigned +# - sink outputs: assigned but never used later as inputs +# - internals: assigned but also used somewhere -> wires +function infer_io(c::Circuit; top_inputs::Union{Nothing,Vector{Symbol}}=nothing, + top_outputs::Union{Nothing,Vector{Symbol}}=nothing) + defs, uses = collect_def_use(c) + defs_set = Set(defs) + uses_set = Set(uses) + + inferred_inputs = collect(setdiff(uses_set, defs_set)) + sinks = collect(setdiff(defs_set, uses_set)) # candidates for final outputs + internals = collect(intersect(defs_set, uses_set)) + + inputs = top_inputs === nothing ? sort(inferred_inputs, by=_nat_key) : top_inputs + outputs = top_outputs === nothing ? sort(sinks, by=_nat_key) : top_outputs + wires = sort(setdiff(defs, outputs), by=_nat_key) # everything assigned except chosen outputs + return inputs, outputs, wires +end + +# Build a stable rename dict for all symbols involved +function build_renames(sc::Circuit, inputs::Vector{Symbol}, outputs::Vector{Symbol}, wires::Vector{Symbol}; + use_tensor_names::Bool=false, orig_to_simplified::Union{Nothing,Dict{Int,Int}}=nothing) + all_syms = Symbol[] + append!(all_syms, inputs, outputs, wires) + # Also catch any remaining RHS-only variables (should be inputs already) + extract_symbols!(sc, all_syms) + # Remove constants + filter!(s -> !(is_true_sym(s) || is_false_sym(s)), all_syms) + unique!(all_syms) + rename = Dict{Symbol,String}() + + if use_tensor_names && !isnothing(orig_to_simplified) + # Build reverse mapping: symbol -> tensor index + symbol_to_tensor = Dict{Symbol, Int}() + for (simplified_idx, orig_idx) in orig_to_simplified + if simplified_idx <= length(sc.exprs) + expr = sc.exprs[simplified_idx] + for o in expr.outputs + symbol_to_tensor[o] = orig_idx + end + end + end + + # Rename based on tensor indices + for s in all_syms + if haskey(symbol_to_tensor, s) + # This is a gate output, name it based on tensor index + rename[s] = "w$(symbol_to_tensor[s])" + else + # This is an input or output, keep original sanitized name + rename[s] = sanitize_name(s) + end + end + else + # Original behavior: sanitize names + for s in all_syms + rename[s] = sanitize_name(s) + end + end + return rename +end + +# --- constraint helpers: detect constant RHS and split constraints --- +is_bool_const(ex::BooleanExpr) = (ex.head == :var) && (is_true_sym(ex.var) || is_false_sym(ex.var)) +const_val(ex::BooleanExpr) = ex.head == :var ? is_true_sym(ex.var) : error("not a const expr") + +""" +Return: + defs::Vector{Assignment} # keep normal assignments + constraints::Vector{Tuple{Symbol,Bool}} # (signal, must_be_value) +We treat any assignment whose RHS is a boolean constant as a constraint rather than an assignment +(to avoid collapsing the logic when exporting to Verilog/AIG). +""" +function split_defs_and_constraints(sc::Circuit) + defs = Assignment[] + cons = Tuple{Symbol,Bool}[] + for ex in sc.exprs + if is_bool_const(ex.expr) + # record constraints for *each* LHS symbol + val = const_val(ex.expr) + for o in ex.outputs + push!(cons, (o, val)) + end + # do NOT emit this as an assignment + else + push!(defs, ex) + end + end + return defs, cons +end + +function circuit_to_verilog(c::Circuit; + module_name::String="circuit", + top_inputs::Union{Nothing,Vector{Symbol}}=nothing, + top_outputs::Union{Nothing,Vector{Symbol}}=nothing, + satisfiable_when_high::Bool=true, + use_tensor_names::Bool=true) + # Ensure simplified form (introduces one-op assignments and temps) + sc = simple_form(c) + + # Split functional defs vs. constant constraints + defs, constraints = split_defs_and_constraints(sc) + + inputs, outputs, wires = infer_io(sc; top_inputs=top_inputs, top_outputs=top_outputs) + + # If there are constraints like o = 0/1, we don't override o; + # instead we add a new SAT output that encodes all constraints. + has_sat = !isempty(constraints) + if has_sat + outputs = Symbol[:sat] + def_syms, _ = collect_def_use(sc) + wires = sort(setdiff(def_syms, outputs), by=_nat_key) + end + + # Build mapping from original circuit expression index to simplified circuit expression + orig_to_simplified = nothing + if use_tensor_names + orig_to_simplified = build_original_to_simplified_mapping(c, sc) + end + + rename = build_renames(sc, inputs, outputs, wires; + use_tensor_names=use_tensor_names, + orig_to_simplified=orig_to_simplified) + if has_sat + rename[:sat] = "sat" + end + + # Header: keep a clean, deterministic port order: inputs then outputs + port_list = join([rename[s] for s in vcat(inputs, outputs)], ", ") + + lines = String[] + push!(lines, "module $(sanitize_name(Symbol(module_name))) ($port_list);") + + # Declarations + if !isempty(inputs) + push!(lines, " input " * join(getindex.(Ref(rename), inputs), ", ") * ";") + end + # outputs already include :sat if needed + if !isempty(outputs) + push!(lines, " output " * join(getindex.(Ref(rename), outputs), ", ") * ";") + end + if !isempty(wires) + push!(lines, " wire " * join(getindex.(Ref(rename), wires), ", ") * ";") + end + push!(lines, "") # blank line + + # Assignments for functional defs only + for (idx, ex) in enumerate(defs) + rhs = verilog_expr(ex.expr, rename) + for o in ex.outputs + oname = rename[o] + push!(lines, " assign $oname = $rhs;") + end + end + + # Encode constraints as a single SAT output: sat = ∧ (o == const) + if has_sat + if satisfiable_when_high + # Normal polarity: sat=1 when all constraints satisfied + terms = String[] + for (o, val) in constraints + on = rename[o] + push!(terms, val ? on : "(~" * on * ")") + end + sat_rhs = isempty(terms) ? "1'b1" : "(" * join(terms, " & ") * ")" + else + # Inverted polarity: sat=0 when all constraints satisfied + # Use OR logic: sat = (o != val_1) | (o != val_2) | ... + # This is equivalent to ~((o == val_1) & (o == val_2) & ...) + terms = String[] + for (o, val) in constraints + on = rename[o] + # Invert each constraint check + push!(terms, val ? "(~" * on * ")" : on) + end + sat_rhs = isempty(terms) ? "1'b0" : "(" * join(terms, " | ") * ")" + end + push!(lines, " assign " * rename[:sat] * " = " * sat_rhs * ";") + end + + push!(lines, "endmodule") + return join(lines, "\n") +end + +# Build a mapping from simplified expression index to original circuit expression index +# This function traces back through simple_form transformations to find the original expression +function build_original_to_simplified_mapping(c::Circuit, sc::Circuit) + # For now, we'll use a simple heuristic: + # Map each simplified expression to its index in the original circuit based on output symbol matching + mapping = Dict{Int, Int}() + + # Build a map from output symbol to original expression index + output_to_orig = Dict{Symbol, Int}() + for (i, expr) in enumerate(c.exprs) + for o in expr.outputs + output_to_orig[o] = i + end + end + + # For each simplified expression, try to find its corresponding original expression + for (j, sexpr) in enumerate(sc.exprs) + # Check if any output symbol matches an original expression + for o in sexpr.outputs + if haskey(output_to_orig, o) + mapping[j] = output_to_orig[o] + break + end + end + # If no match found, use the simplified index + if !haskey(mapping, j) + mapping[j] = j + end + end + + return mapping +end + +# Convenience writer +function write_verilog(io::IO, c::Circuit; kwargs...) + print(io, circuit_to_verilog(c; kwargs...)) +end + +function write_verilog(path::AbstractString, c::Circuit; kwargs...) + open(path, "w") do io + write_verilog(io, c; kwargs...) + end + return path +end diff --git a/benchmarks/src/comparison.jl b/benchmarks/src/comparison.jl new file mode 100644 index 0000000..dd210d7 --- /dev/null +++ b/benchmarks/src/comparison.jl @@ -0,0 +1,141 @@ +# run_full_benchmark is deprecated - implement per problem type + +function list_available_solvers(problem_type::Type{<:AbstractBenchmarkProblem}) + solvers = available_solvers(problem_type) + default = default_solver(problem_type) + default_name = solver_name(default) + + println("Available solvers for $(string(problem_type)):") + for solver in solvers + name = solver_name(solver) + is_default = name == default_name ? " (default)" : "" + println(" • $name$is_default") + end +end + +""" + run_solver_comparison(problem_type, dataset_paths; solvers=nothing) + +Compare multiple solvers on the same datasets. +""" +function run_solver_comparison(problem_type::Type{<:AbstractBenchmarkProblem}, + dataset_paths::Vector{String}; + solvers=nothing) + solver_list = isnothing(solvers) ? available_solvers(problem_type) : solvers + results = Dict() + + for solver in solver_list + solver_name_str = solver_name(solver) + @info "Running benchmark with solver: $solver_name_str" + + solver_results = [] + for path in dataset_paths + result = benchmark_dataset(problem_type, path; solver) + if !isnothing(result) + push!(solver_results, result) + end + end + + results[solver_name_str] = solver_results + end + + return results +end + +function print_solver_comparison_summary(solver_results::Dict) + if isempty(solver_results) + println("No solver results to compare!") + return + end + + println("\n" * repeat("=", 80)) + println("SOLVER COMPARISON SUMMARY") + println(repeat("=", 80)) + + # Get all successful results (non-nothing) + first_results = first(values(solver_results)) + successful_results = filter(!isnothing, first_results) + + if isempty(successful_results) + println("No successful results to compare!") + return + end + + # Get unique dataset paths + dataset_paths = unique([result["dataset_path"] for result in successful_results]) + + solver_names = collect(keys(solver_results)) + header = "| Dataset " + for solver in solver_names + header *= "| $(rpad(solver, 12)) " + end + header *= "|" + + println(header) + println(repeat("-", length(header))) + + for path in dataset_paths + # Extract a short name from path + path_parts = splitpath(path) + dataset_name = length(path_parts) > 0 ? path_parts[end] : path + dataset_str = dataset_name[1:min(length(dataset_name), 17)] + row = "| $(rpad(dataset_str, 17)) " + + for solver in solver_names + solver_result = solver_results[solver] + matching_result = findfirst(r -> !isnothing(r) && r["dataset_path"] == path, solver_result) + + if isnothing(matching_result) + time_str = "FAILED" + else + result = solver_result[matching_result] + median_time = result["median_time"] + time_str = format_time(median_time) + end + row *= "| $(rpad(time_str, 12)) " + end + row *= "|" + println(row) + end + println(repeat("=", 80)) +end + +function compare_solver_results(name1::String, results1, name2::String, results2) + println(repeat("+", 80)) + println("| $(rpad("Dataset", 15)) | $(rpad("$name1", 15)) | $(rpad("$name2", 15)) | $(rpad("Winner", 20)) |") + println(repeat("-", 80)) + + success1 = filter(!isnothing, results1) + success2 = filter(!isnothing, results2) + + for r1 in success1 + path = r1["dataset_path"] + path_parts = splitpath(path) + dataset_name = length(path_parts) > 0 ? path_parts[end] : path + dataset_str = dataset_name[1:min(length(dataset_name), 15)] + + r2 = findfirst(r -> r["dataset_path"] == path, success2) + if !isnothing(r2) + r2 = success2[r2] + + time1 = r1["median_time"] + time2 = r2["median_time"] + time1_str = format_time(time1) + time2_str = format_time(time2) + + if time1 < time2 + winner_str = "$name1 ($(round(time2/time1, digits=2))x)" + elseif time2 < time1 + winner_str = "$name2 ($(round(time1/time2, digits=2))x)" + else + winner_str = "Tie" + end + + println("| $(rpad(dataset_str, 15)) | $(rpad(time1_str, 15)) | $(rpad(time2_str, 15)) | $(rpad(winner_str, 20)) |") + else + println("| $(rpad(dataset_str, 15)) | $(rpad(format_time(r1["median_time"]), 15)) | $(rpad("FAILED", 15)) | $(rpad(name1, 20)) |") + end + end + println(repeat("+", 80)) +end + diff --git a/benchmarks/src/factoring/dataset.jl b/benchmarks/src/factoring/dataset.jl new file mode 100644 index 0000000..4ae9e15 --- /dev/null +++ b/benchmarks/src/factoring/dataset.jl @@ -0,0 +1,105 @@ +function dataset_metadata_hash(config::FactoringConfig, per_config::Int, seed::UInt64, include_solution::Bool) + parts = [ + string(config.m), + string(config.n), + string(per_config), + string(seed), + string(include_solution) + ] + content = join(parts, "|") + return bytes2hex(sha256(content))[1:16] +end + +function deterministic_seed(config::FactoringConfig) + content = "$(config.m)|$(config.n)" + hash_bytes = sha256(content) + return reinterpret(UInt64, hash_bytes[1:8])[1] +end + +function check_dataset_compatibility(path::String, config::FactoringConfig, + per_config::Int, seed::UInt64, include_solution::Bool) + if !isfile(path) + return false + end + + try + instances = read_instances(FactoringProblem, path) + isempty(instances) && return false + + if length(instances) != per_config + @info " Dataset size mismatch: expected $per_config, found $(length(instances))" + return false + end + + expected_hash = dataset_metadata_hash(config, per_config, seed, include_solution) + if has_metadata(instances[1], expected_hash) + @info " Found compatible dataset: $path" + return true + else + @info " Dataset parameter mismatch: $path" + return false + end + catch e + @warn " Error checking dataset: $e" + return false + end +end + + +function generate_factoring_datasets(configs::Vector{FactoringConfig}; + per_config::Int=100, + include_solution::Bool=false, + force_regenerate::Bool=false) + + outdir = resolve_data_dir("factoring") + isdir(outdir) || mkpath(outdir) + + paths = String[] + generated_count = 0 + reused_count = 0 + + for config in configs + filename = filename_pattern(FactoringProblem, config) + path = joinpath(outdir, filename) + + seed = deterministic_seed(config) + + if !force_regenerate && check_dataset_compatibility(path, config, per_config, seed, include_solution) + @info "Reusing existing dataset: $path" + push!(paths, path) + reused_count += 1 + continue + end + + @info "Generating new dataset: $path" + @info " Using seed: $seed, instances: $per_config" + + seeded_rng = Random.Xoshiro(seed) + metadata_hash = dataset_metadata_hash(config, per_config, seed, include_solution) + + open(path, "w") do io + for i in 1:per_config + instance = generate_instance(FactoringProblem, config; + rng=seeded_rng, + include_solution=include_solution) + + if i == 1 + instance = add_metadata(instance; + metadata_hash=metadata_hash, + generation_seed=seed, + generation_timestamp=nothing) + end + + write_instance(io, instance) + end + end + + @info "Generated dataset: $path" + push!(paths, path) + generated_count += 1 + end + + @info "Dataset generation complete: $generated_count generated, $reused_count reused" + return paths +end + diff --git a/benchmarks/src/factoring/factor_IP.jl b/benchmarks/src/factoring/factor_IP.jl deleted file mode 100644 index 14bf4e5..0000000 --- a/benchmarks/src/factoring/factor_IP.jl +++ /dev/null @@ -1,82 +0,0 @@ -""" - findmin(problem::AbstractProblem, optimizer, tag::Bool) - -Build and solve a 0-1 integer program corresponding to a CircuitSAT -`problem` using the specified `optimizer` (e.g., `SCIP.Optimizer`). - -Details: -- Creates binary variables `x` representing boolean assignments. -- For each circuit constraint, forbids assignments that violate the truth table - by adding linear inequalities (one per forbidden pattern). -- If the problem provides objectives, aggregates them and either minimizes or - maximizes depending on `tag` (true => minimize, false => maximize). If no - objectives exist, uses a dummy objective 0. - -Returns a vector of 0/1 integers representing the optimal assignment. -""" -function findmin(problem::AbstractProblem,tag::Bool; optimizer, env) - # Extract circuit constraints and sizing info - cons = constraints(problem) - nsc = ProblemReductions.num_variables(problem) - maxN = maximum([length(c.variables) for c in cons]) - combs = [ProblemReductions.combinations(2,i) for i in 1:maxN] - - objs = objectives(problem) - - # IP by JuMP - opt = isnothing(env) ? optimizer() : optimizer(env) - model = JuMP.direct_model(opt) - set_silent(model) - set_string_names_on_creation(model, false) - - # Binary assignment variables for each SAT variable - JuMP.@variable(model, 0 <= x[i = 1:nsc] <= 1, Int) - - # Add constraints that eliminate forbidden patterns in the truth tables. - for con in cons - f_vec = findall(!,con.specification) - num_vars = length(con.variables) - for f in f_vec - # Sum of matching literals must be <= num_vars - 1 (i.e., at least one differs) - JuMP.@constraint(model, sum(j-> iszero(combs[num_vars][f][j]) ? (1 - x[con.variables[j]]) : x[con.variables[j]], 1:num_vars) <= num_vars -1) - end - end - if isempty(objs) - # No objective: use a constant objective to just find a feasible point - JuMP.@objective(model, Min, 0) - else - # Aggregate objective pieces based on the boolean of the associated variable - obj_sum = sum(objs) do obj - (1-x[obj.variables[1]])*obj.specification[1] + x[obj.variables[1]]*obj.specification[2] - end - tag ? JuMP.@objective(model, Min, obj_sum) : JuMP.@objective(model, Max, obj_sum) - end - - # Solve and validate feasibility - JuMP.optimize!(model) - @assert JuMP.is_solved_and_feasible(model) "The problem is infeasible" - # Return the 0/1 assignment as integers - return round.(Int, JuMP.value.(x)) -end - -""" - factoring(m, n, N) - -Construct a `Factoring(m, n, N)` problem, reduce it to `CircuitSAT`, solve the -resulting IP, and decode the boolean assignment back into integer factors. - -Arguments: -- `m`, `n`: bit-widths of the two factors (p and q) -- `N`: the composite integer to factor - -Returns a pair `(a, b)` such that `a * b == N` on success. -""" -function factoring(m, n, N; optimizer, env) - fact3 = Factoring(m, n, N) - res3 = reduceto(CircuitSAT, fact3) - # Build CircuitSAT instance with constraints enabled - problem = CircuitSAT(res3.circuit.circuit; use_constraints=true); - # Solve for boolean assignments then reconstruct integer factors - vals = findmin(problem, true; optimizer, env) - return ProblemReductions.read_solution(fact3, [vals[res3.p]...,vals[res3.q]...]) -end \ No newline at end of file diff --git a/benchmarks/src/factoring/factoring_problem.jl b/benchmarks/src/factoring/factoring_problem.jl deleted file mode 100644 index 918ae60..0000000 --- a/benchmarks/src/factoring/factoring_problem.jl +++ /dev/null @@ -1,205 +0,0 @@ -using Random -using Primes -using JSON3 -using SHA: bytes2hex, sha1 -using BooleanInference -using JuMP -using Gurobi -using ..BooleanInferenceBenchmarks: resolve_data_dir, AbstractSolver - -# Define solver types for factoring problems -""" -BooleanInference solver - uses the original BooleanInference.solve_factoring method -""" -struct BooleanInferenceSolver <: AbstractSolver end - -""" -Integer Programming solver - uses JuMP with Gurobi optimizer -""" -struct IPSolver <: AbstractSolver - optimizer::Any - env::Any - - function IPSolver(optimizer=Gurobi.Optimizer, env=nothing) - new(optimizer, env) - end -end - -# Prime generation utilities -function random_semiprime(m::Int, n::Int; rng=Random.GLOBAL_RNG, distinct::Bool=false) - # Generate two random primes of approximately m and n bits - p = random_prime_bits(m; rng) - q = random_prime_bits(n; rng) - - if distinct && p == q - # If they're the same and we want distinct, regenerate q - while q == p - q = random_prime_bits(n; rng) - end - end - - N = p * q - return p, q, N -end - -function random_prime_bits(bits::Int; rng=Random.GLOBAL_RNG) - # Generate a random prime with approximately 'bits' bits - if bits <= 1 - return 2 - end - - # Range for n-bit numbers: 2^(n-1) to 2^n - 1 - min_val = 2^(bits-1) - max_val = 2^bits - 1 - - # Find a prime in this range - for _ in 1:1000 # Max attempts - candidate = rand(rng, min_val:max_val) - if candidate % 2 == 0 - candidate += 1 # Make it odd - end - - if isprime(candidate) - return candidate - end - end - - # Fallback: use nextprime - candidate = rand(rng, min_val:max_val) - return nextprime(candidate) -end - -""" -Factoring problem type. -""" -struct FactoringProblem <: AbstractBenchmarkProblem end - -""" -Configuration for factoring problems. -""" -struct FactoringConfig <: AbstractProblemConfig - m::Int # bits for first factor - n::Int # bits for second factor -end - -# Implement the interface for FactoringProblem - -function generate_instance(::Type{FactoringProblem}, config::FactoringConfig; - rng::AbstractRNG=Random.GLOBAL_RNG, - include_solution::Bool=false) - p, q, N = random_semiprime(config.m, config.n; rng=rng, distinct=true) - - instance = Dict( - "m" => config.m, - "n" => config.n, - "N" => string(N), - "id" => problem_id(FactoringProblem, config, N) - ) - # include the solution in the dataset if requested - if include_solution - instance["p"] = string(p) - instance["q"] = string(q) - end - - return instance -end - -# Implement solve_instance with different solvers -function solve_instance(::Type{FactoringProblem}, instance, solver::BooleanInferenceSolver) - m = instance["m"] - n = instance["n"] - N = parse(Int, instance["N"]) - return BooleanInference.solve_factoring(m, n, N) -end - -function solve_instance(::Type{FactoringProblem}, instance, solver::IPSolver) - m = instance["m"] - n = instance["n"] - N = parse(Int, instance["N"]) - return factoring(m, n, N; optimizer=solver.optimizer, env=solver.env) -end - -""" - verify_solution(::Type{FactoringProblem}, instance, result) - -Verify that the factoring result is correct. -""" -function verify_solution(::Type{FactoringProblem}, instance, result) - try - N = parse(Int, instance["N"]) - - # Handle different result formats - if result isa Tuple && length(result) == 2 - p, q = result - else - @warn "Unknown result format: $(typeof(result))" - return false - end - - # Verify the factorization - if p * q == N - return true - else - @warn "Incorrect factorization: $p × $q = $(p*q) ≠ $N" - return false - end - catch e - @warn "Error verifying solution: $e" - return false - end -end - -function run_full_benchmark(::Type{FactoringProblem}, - input_configs::Union{Vector{Tuple{Int,Int}}, Nothing}=nothing; - dataset_per_config::Int=50, - solver=nothing) - - if isnothing(input_configs) - configs = default_configs(FactoringProblem) - else - configs = [config(FactoringProblem, config_item) for config_item in input_configs] - end - - results = benchmark_backend(FactoringProblem, configs; dataset_per_config=dataset_per_config, solver=solver) - - return results -end - - -function problem_id(::Type{FactoringProblem}, config::FactoringConfig, N::Integer) - h = bytes2hex(sha1(string(config.m, "|", config.n, "|", N))) - return h[1:16] # 16-char prefix -end - -function default_configs(::Type{FactoringProblem}) - return FactoringConfig[ - FactoringConfig(10, 10), - FactoringConfig(12, 12), - FactoringConfig(14, 14), - ] -end - -function config(::Type{FactoringProblem}, params::Tuple{Int,Int}) - return FactoringConfig(params[1], params[2]) -end - -function filename_pattern(::Type{FactoringProblem}, config::FactoringConfig) - return "numbers_$(config.m)x$(config.n).jsonl" -end - -# Solver interface implementations -function available_solvers(::Type{FactoringProblem}) - return [BooleanInferenceSolver(), IPSolver()] -end - -function default_solver(::Type{FactoringProblem}) - return BooleanInferenceSolver() -end - -function solver_name(::BooleanInferenceSolver) - return "BooleanInference" -end - -function solver_name(solver::IPSolver) - return "IP-$(solver.optimizer)" -end diff --git a/benchmarks/src/factoring/generators.jl b/benchmarks/src/factoring/generators.jl new file mode 100644 index 0000000..c37d4d8 --- /dev/null +++ b/benchmarks/src/factoring/generators.jl @@ -0,0 +1,58 @@ +function random_prime_bits(bits::Int; rng=Random.GLOBAL_RNG) + if bits <= 1 + return 2 + end + min_val = 2^(bits-1) + max_val = 2^bits - 1 + for _ in 1:1000 + candidate = rand(rng, min_val:max_val) + if candidate % 2 == 0 + candidate += 1 + end + if isprime(candidate) + return candidate + end + end + candidate = rand(rng, min_val:max_val) + return nextprime(candidate) +end + +function random_semiprime(m::Int, n::Int; rng=Random.GLOBAL_RNG, distinct::Bool=false) + p = random_prime_bits(m; rng) + q = random_prime_bits(n; rng) + if distinct && p == q + while q == p + q = random_prime_bits(n; rng) + end + end + N = p * q + return p, q, N +end + +function generate_instance(::Type{FactoringProblem}, config::FactoringConfig; + rng::AbstractRNG=Random.GLOBAL_RNG, + include_solution::Bool=false) + p, q, N = random_semiprime(config.m, config.n; rng=rng, distinct=true) + id = problem_id(config, N) + + return FactoringInstance( + config.m, config.n, N, id; + p = include_solution ? p : nothing, + q = include_solution ? q : nothing + ) +end + +# Helper function to add metadata to an instance +function add_metadata(instance::FactoringInstance; + metadata_hash=nothing, + generation_seed=nothing, + generation_timestamp=nothing) + return FactoringInstance( + instance.m, instance.n, instance.N, instance.id; + p=instance.p, q=instance.q, + metadata_hash=metadata_hash, + generation_seed=generation_seed, + generation_timestamp=generation_timestamp + ) +end + diff --git a/benchmarks/src/factoring/interface.jl b/benchmarks/src/factoring/interface.jl new file mode 100644 index 0000000..9c69a85 --- /dev/null +++ b/benchmarks/src/factoring/interface.jl @@ -0,0 +1,23 @@ +# ---------------------------------------- +# Factoring Problem Interface +# ---------------------------------------- + +# Helper for generating problem IDs +function problem_id(config::FactoringConfig, N::Integer) + h = bytes2hex(sha1(string(config.m, "|", config.n, "|", N))) + return h[1:16] +end + +# Dataset filename pattern +function filename_pattern(::Type{FactoringProblem}, config::FactoringConfig) + return "numbers_$(config.m)x$(config.n).jsonl" +end + +# Solvers +function available_solvers(::Type{FactoringProblem}) + return [BooleanInferenceSolver(), IPSolver(), IPSolver(HiGHS.Optimizer), XSATSolver()] +end + +function default_solver(::Type{FactoringProblem}) + return BooleanInferenceSolver() +end \ No newline at end of file diff --git a/benchmarks/src/factoring/solver_ip.jl b/benchmarks/src/factoring/solver_ip.jl new file mode 100644 index 0000000..7617853 --- /dev/null +++ b/benchmarks/src/factoring/solver_ip.jl @@ -0,0 +1,44 @@ +function findmin_ip(problem::AbstractProblem, tag::Bool; optimizer, env) + cons = constraints(problem) + nsc = ProblemReductions.num_variables(problem) + maxN = maximum([length(c.variables) for c in cons]) + combs = [ProblemReductions.combinations(2,i) for i in 1:maxN] + objs = objectives(problem) + + opt = isnothing(env) ? optimizer() : optimizer(env) + model = JuMP.direct_model(opt) + set_silent(model) + set_string_names_on_creation(model, false) + + JuMP.@variable(model, 0 <= x[i = 1:nsc] <= 1, Int) + + for con in cons + f_vec = findall(!,con.specification) + num_vars = length(con.variables) + for f in f_vec + JuMP.@constraint(model, sum(j-> iszero(combs[num_vars][f][j]) ? (1 - x[con.variables[j]]) : x[con.variables[j]], 1:num_vars) <= num_vars -1) + end + end + + if isempty(objs) + JuMP.@objective(model, Min, 0) + else + obj_sum = sum(objs) do obj + (1-x[obj.variables[1]])*obj.specification[1] + x[obj.variables[1]]*obj.specification[2] + end + tag ? JuMP.@objective(model, Min, obj_sum) : JuMP.@objective(model, Max, obj_sum) + end + + JuMP.optimize!(model) + @assert JuMP.is_solved_and_feasible(model) "The problem is infeasible" + return round.(Int, JuMP.value.(x)) +end + +function factoring_ip(m, n, N; solver::IPSolver) + fact3 = Factoring(m, n, N) + res3 = reduceto(CircuitSAT, fact3) + problem = CircuitSAT(res3.circuit.circuit; use_constraints=true); + vals = findmin_ip(problem, true; optimizer=solver.optimizer, env=solver.env) + return ProblemReductions.read_solution(fact3, [vals[res3.p]...,vals[res3.q]...]) +end + diff --git a/benchmarks/src/factoring/solver_xsat.jl b/benchmarks/src/factoring/solver_xsat.jl new file mode 100644 index 0000000..a0a0f26 --- /dev/null +++ b/benchmarks/src/factoring/solver_xsat.jl @@ -0,0 +1,37 @@ +struct XSatResult + status::Symbol # :sat | :unsat | :unknown + model::Union{Nothing,Dict{Int,Bool}} # Variable -> value (present only when SAT) + raw::String # Raw console output for debugging/logging +end + +""" + run_xsat_and_parse(exec::AbstractString, aig_path::AbstractString) -> XSatResult + +Run XSAT (e.g., exec="/path/to/xsat") and parse the output of `-i aig_path`. +Parse the `s ...` status line and one or more `v ...` assignment lines (DIMACS style, terminated by 0). +""" +function run_xsat_and_parse(exec::AbstractString, aig_path::AbstractString)::XSatResult + raw = read(`$exec -i $aig_path`, String) + + status = + occursin(r"(?m)^s\s+SATISFIABLE\b", raw) ? :sat : + occursin(r"(?m)^s\s+UNSATISFIABLE\b", raw) ? :unsat : :unknown + + model = nothing + if status == :sat + lits = Int[] + for m in eachmatch(r"(?m)^v\s+([^\n]+)", raw) + append!(lits, parse.(Int, split(m.captures[1]))) + end + filter!(!=(0), lits) + + d = Dict{Int,Bool}() + for lit in lits + v = abs(lit) + d[v] = lit > 0 + end + model = d + end + + return XSatResult(status, model, raw) +end diff --git a/benchmarks/src/factoring/solvers.jl b/benchmarks/src/factoring/solvers.jl new file mode 100644 index 0000000..a507d6c --- /dev/null +++ b/benchmarks/src/factoring/solvers.jl @@ -0,0 +1,56 @@ +function solve_instance(::Type{FactoringProblem}, instance::FactoringInstance, solver::BooleanInferenceSolver) + return BooleanInference.solve_factoring(instance.m, instance.n, Int(instance.N)) +end + +function solve_instance(::Type{FactoringProblem}, instance::FactoringInstance, solver::IPSolver) + return factoring_ip(instance.m, instance.n, Int(instance.N); solver) +end + +function solve_instance(::Type{FactoringProblem}, instance::FactoringInstance, solver::XSATSolver) + m, n = instance.m, instance.n + fproblem = Factoring(m, n, instance.N) + circuit_sat = reduceto(CircuitSAT, fproblem) + + mktempdir() do dir + vfile = joinpath(dir, "circuit.v") + aig = joinpath(dir, "circuit.aig") + + write_verilog(vfile, circuit_sat.circuit.circuit) + + run(`$(solver.yosys_path) -q -p "read_verilog $vfile; prep -top circuit; aigmap; write_aiger -symbols $aig"`) + + res = run_xsat_and_parse(solver.csat_path, aig) + res.status != :sat && return :unsat + + model = res.model::Dict{Int,Bool} + + bits_p = [get(model, i, false) for i in 1:m] + bits_q = [get(model, i, false) for i in m+1:m+n] + + p_val = sum(Int(bit) << (i-1) for (i, bit) in enumerate(bits_p)) + q_val = sum(Int(bit) << (i-1) for (i, bit) in enumerate(bits_q)) + + return (p_val, q_val) + end +end + +function verify_solution(::Type{FactoringProblem}, instance::FactoringInstance, result) + try + if result isa Tuple + p = result[1]; q = result[2]; + else + @warn "Unknown result format: $(typeof(result))" + return false + end + if p * q == instance.N + return true + else + @warn "Incorrect factorization: $p × $q = $(p*q) ≠ $(instance.N)" + return false + end + catch e + @warn "Error verifying solution: $e" + return false + end +end + diff --git a/benchmarks/src/factoring/types.jl b/benchmarks/src/factoring/types.jl new file mode 100644 index 0000000..2771611 --- /dev/null +++ b/benchmarks/src/factoring/types.jl @@ -0,0 +1,100 @@ +struct FactoringProblem <: AbstractBenchmarkProblem end + +struct FactoringConfig <: AbstractProblemConfig + m::Int + n::Int +end + +struct FactoringInstance <: AbstractInstance + m::Int + n::Int + N::BigInt + id::String + # Optional solution fields + p::Union{BigInt, Nothing} + q::Union{BigInt, Nothing} + # Optional metadata fields + metadata_hash::Union{String, Nothing} + generation_seed::Union{UInt64, Nothing} + generation_timestamp::Union{String, Nothing} + + function FactoringInstance(m::Int, n::Int, N::Integer, id::String; + p=nothing, q=nothing, + metadata_hash=nothing, + generation_seed=nothing, + generation_timestamp=nothing) + new(m, n, BigInt(N), id, + isnothing(p) ? nothing : BigInt(p), + isnothing(q) ? nothing : BigInt(q), + metadata_hash, generation_seed, generation_timestamp) + end +end + + +# ---------------------------------------- +# Dataset I/O Implementation (JSON format) +# ---------------------------------------- + +# Write a single instance to IO (JSON format) +function write_instance(io::IO, instance::FactoringInstance) + obj = Dict{String, Any}( + "m" => instance.m, + "n" => instance.n, + "N" => string(instance.N), + "id" => instance.id + ) + + if !isnothing(instance.p) + obj["p"] = string(instance.p) + end + if !isnothing(instance.q) + obj["q"] = string(instance.q) + end + if !isnothing(instance.metadata_hash) + obj["_metadata_hash"] = instance.metadata_hash + end + if !isnothing(instance.generation_seed) + obj["_generation_seed"] = instance.generation_seed + end + if !isnothing(instance.generation_timestamp) + obj["_generation_timestamp"] = instance.generation_timestamp + end + + JSON3.write(io, obj) + write(io, '\n') +end + +# Read instances from file (JSON format) +function read_instances(::Type{FactoringProblem}, path::AbstractString) + instances = FactoringInstance[] + open(path, "r") do io + for line in eachline(io) + isempty(strip(line)) && continue + data = JSON3.read(line) + + m = data["m"] + n = data["n"] + N = parse(BigInt, data["N"]) + id = data["id"] + + p = haskey(data, "p") ? parse(BigInt, data["p"]) : nothing + q = haskey(data, "q") ? parse(BigInt, data["q"]) : nothing + metadata_hash = haskey(data, "_metadata_hash") ? data["_metadata_hash"] : nothing + generation_seed = haskey(data, "_generation_seed") ? UInt64(data["_generation_seed"]) : nothing + generation_timestamp = haskey(data, "_generation_timestamp") ? data["_generation_timestamp"] : nothing + + instance = FactoringInstance(m, n, N, id; + p=p, q=q, + metadata_hash=metadata_hash, + generation_seed=generation_seed, + generation_timestamp=generation_timestamp) + push!(instances, instance) + end + end + return instances +end + +# Check metadata +function has_metadata(instance::FactoringInstance, expected_hash::String) + return !isnothing(instance.metadata_hash) && instance.metadata_hash == expected_hash +end \ No newline at end of file diff --git a/benchmarks/src/formatting.jl b/benchmarks/src/formatting.jl new file mode 100644 index 0000000..6b08bd6 --- /dev/null +++ b/benchmarks/src/formatting.jl @@ -0,0 +1,72 @@ +function format_bytes(bytes::Number) + if bytes < 1024 + return "$(Int(bytes))B" + elseif bytes < 1024^2 + return "$(round(bytes/1024, digits=1))KB" + elseif bytes < 1024^3 + return "$(round(bytes/1024^2, digits=1))MB" + else + return "$(round(bytes/1024^3, digits=1))GB" + end +end + +function format_time(seconds::Float64) + if seconds < 1e-6 + return "$(round(seconds * 1e9, digits=1))ns" + elseif seconds < 1e-3 + return "$(round(seconds * 1e6, digits=1))us" + elseif seconds < 1 + return "$(round(seconds * 1e3, digits=1))ms" + elseif seconds < 60 + return "$(round(seconds, digits=2))s" + else + minutes = floor(seconds / 60) + secs = seconds - minutes * 60 + return "$(Int(minutes))m$(round(secs, digits=1))s" + end +end + +function print_benchmark_summary(results) + println("\n" * repeat("=", 60)) + println("BENCHMARK SUMMARY") + println(repeat("=", 60)) + + successful = filter(r -> r["status"] == "success", results) + failed = filter(r -> r["status"] == "failed", results) + + if !isempty(successful) + println("Successful benchmarks: $(length(successful))") + println(repeat("+", 90)) + println("| Config | Median Time | Memory | Instances | Correct | Accuracy | Timed |") + println(repeat("-", 90)) + + for result in successful + config = result["config"] + median_time = result["median_time"] + median_memory = result["median_memory"] + instances = get(result, "instances_tested", 0) + correct = get(result, "correct_runs", 0) + accuracy = get(result, "accuracy_rate", 0.0) + timed = get(result, "successful_runs", 0) + + config_str = "$(config.m)x$(config.n)" + time_str = format_time(median_time) + memory_str = format_bytes(median_memory) + accuracy_str = "$(round(accuracy * 100, digits=1))%" + + println("| $(rpad(config_str, 11)) | $(rpad(time_str, 11)) | $(rpad(memory_str, 11)) | $(rpad(string(instances), 9)) | $(rpad(string(correct), 7)) | $(rpad(accuracy_str, 8)) | $(rpad(string(timed), 5)) |") + end + println(repeat("+", 90)) + end + + if !isempty(failed) + println("\nFailed benchmarks: $(length(failed))") + for result in failed + config = result["config"] + reason = result["reason"] + println(" - $(config): $reason") + end + end + println(repeat("=", 60)) +end + diff --git a/benchmarks/src/generic_benchmark.jl b/benchmarks/src/generic_benchmark.jl deleted file mode 100644 index c8b0924..0000000 --- a/benchmarks/src/generic_benchmark.jl +++ /dev/null @@ -1,834 +0,0 @@ -# Generic benchmark framework using multiple dispatch - -""" - dataset_metadata_hash(problem_type, config, per_config, seed, include_solution) - -Generate a stable hash for dataset parameters to enable caching and reproducibility. -""" -function dataset_metadata_hash(problem_type::Type{<:AbstractBenchmarkProblem}, - config::AbstractProblemConfig, - per_config::Int, - seed::UInt64, - include_solution::Bool) - # Create a deterministic string representation of all parameters - parts = [ - string(problem_type), - string(config), - string(per_config), - string(seed), - string(include_solution) - ] - content = join(parts, "|") - return bytes2hex(sha256(content))[1:16] # 16-char prefix -end - -""" - deterministic_seed(problem_type, config) - -Generate a deterministic seed for reproducible dataset generation. -""" -function deterministic_seed(problem_type::Type{<:AbstractBenchmarkProblem}, - config::AbstractProblemConfig) - # Create deterministic seed from problem type and config - content = string(problem_type) * "|" * string(config) - hash_bytes = sha256(content) - # Convert first 8 bytes to UInt64 for seed - return reinterpret(UInt64, hash_bytes[1:8])[1] -end - -""" - check_dataset_compatibility(dataset_path, problem_type, config, per_config, seed, include_solution) - -Check if existing dataset is compatible with current parameters. -""" -function check_dataset_compatibility(dataset_path::String, - problem_type::Type{<:AbstractBenchmarkProblem}, - config::AbstractProblemConfig, - per_config::Int, - seed::UInt64, - include_solution::Bool) - if !isfile(dataset_path) - return false - end - - try - # Read first line to get metadata - lines = readlines(dataset_path) - if isempty(lines) - return false - end - - # Check if we have the expected number of lines - if length(lines) != per_config - @info " Dataset size mismatch: expected $per_config instances, found $(length(lines))" - return false - end - - # Try to parse first instance and check if it has metadata - first_instance = JSON3.read(lines[1]) - expected_hash = dataset_metadata_hash(problem_type, config, per_config, seed, include_solution) - - # Check if instance has metadata hash - if haskey(first_instance, "_metadata_hash") - if first_instance["_metadata_hash"] == expected_hash - @info " Found compatible dataset: $dataset_path" - return true - else - @info " Dataset parameter mismatch: $dataset_path" - return false - end - else - # Legacy dataset without metadata - consider incompatible - @info " Legacy dataset without metadata: $dataset_path" - return false - end - catch e - @warn " Error checking dataset compatibility: $e" - return false - end -end - -""" - generate_datasets(problem_type::Type{<:AbstractBenchmarkProblem}; - configs::Vector{<:AbstractProblemConfig}, - per_config=100, - include_solution=false, - force_regenerate=false) - -Generate datasets for a given problem type with reproducible, cached results. -If datasets already exist with compatible parameters, they will be reused. -Uses deterministic seeding for reproducibility - no external RNG needed. -""" -function generate_datasets(problem_type::Type{<:AbstractBenchmarkProblem}; - configs::Vector{<:AbstractProblemConfig}, - per_config::Int=100, - include_solution::Bool=false, - force_regenerate::Bool=false) - - # Get output directory based on problem type - outdir = resolve_data_dir(lowercase(string(problem_type)[1:end-7])) # Remove "Problem" suffix - isdir(outdir) || mkpath(outdir) - - paths = String[] - generated_count = 0 - reused_count = 0 - - for config in configs - filename = filename_pattern(problem_type, config) - path = joinpath(outdir, filename) - - # Generate deterministic seed for reproducibility - seed = deterministic_seed(problem_type, config) - - # Check if compatible dataset already exists - if !force_regenerate && check_dataset_compatibility(path, problem_type, config, per_config, seed, include_solution) - @info "Reusing existing dataset: $path" - push!(paths, path) - reused_count += 1 - continue - end - - @info "Generating new dataset: $path" - @info " Using deterministic seed: $seed" - - # Create seeded RNG for reproducible generation - seeded_rng = Random.Xoshiro(seed) - metadata_hash = dataset_metadata_hash(problem_type, config, per_config, seed, include_solution) - - open(path, "w") do io - for i in 1:per_config - instance = generate_instance(problem_type, config; rng=seeded_rng, include_solution=include_solution) - - # Add metadata to first instance - if i == 1 - instance = merge(instance, Dict( - "_metadata_hash" => metadata_hash, - "_generation_seed" => seed, - "_generation_timestamp" => string(now()) - )) - end - - JSON3.write(io, instance) - write(io, '\n') - end - end - - @info "Generated dataset: $path (metadata hash: $metadata_hash)" - push!(paths, path) - generated_count += 1 - end - - @info "Dataset generation complete: $(generated_count) generated, $(reused_count) reused in: $outdir" - return paths -end - -""" - benchmark_problem(problem_type::Type{<:AbstractBenchmarkProblem}; - configs=default_configs(problem_type), - solver=nothing) - -Simple and direct benchmark using pre-generated datasets. -""" -function benchmark_problem(problem_type::Type{<:AbstractBenchmarkProblem}; - configs=default_configs(problem_type), - solver=nothing) - - results = [] - - # Use specified solver or default - actual_solver = isnothing(solver) ? default_solver(problem_type) : solver - solver_info = solver_name(actual_solver) - - @info "Using solver: $solver_info" - - for config in configs - @info "Benchmarking config: $config" - - # Load instances from dataset - result = benchmark_dataset_instances(problem_type, config, actual_solver, solver_info) - push!(results, result) - end - - return results -end - -""" - benchmark_dataset_instances(problem_type, config, solver, solver_info) - -Benchmark using all instances from pre-generated dataset. -""" -function benchmark_dataset_instances(problem_type::Type{<:AbstractBenchmarkProblem}, - config::AbstractProblemConfig, - solver, - solver_info::String) - try - # Load instances from dataset file - problem_name = lowercase(string(problem_type)[1:end-7]) # Remove "Problem" suffix - filename = filename_pattern(problem_type, config) - dataset_path = joinpath(resolve_data_dir(problem_name), filename) - - if !isfile(dataset_path) - @error "Dataset file not found: $dataset_path" - return create_failed_result(config, "dataset_not_found: $dataset_path") - end - - @info " Loading instances from: $dataset_path" - instances = read_jsonl(dataset_path) - @info " Testing $(length(instances)) instances" - - if isempty(instances) - @error " No instances found in dataset file" - return create_failed_result(config, "empty_dataset") - end - - # Benchmark each instance once - all_times = Float64[] - all_memory = Int64[] - all_gc_times = Float64[] - successful_runs = 0 - failed_runs = 0 - - # Pre-tune benchmark once for all instances to save time - sample_trial_fn = () -> solve_instance(problem_type, instances[1], solver) - b = @benchmarkable $sample_trial_fn() - tune!(b) - @info " Benchmark tuned, starting verification and timing..." - - correct_runs = 0 - incorrect_runs = 0 - - for (i, instance) in enumerate(instances) - try - # First, verify solution correctness (without timing) - result = solve_instance(problem_type, instance, solver) - is_correct = verify_solution(problem_type, instance, result) - - if !is_correct - @warn " Instance $i: Incorrect solution, skipping timing" - incorrect_runs += 1 - continue - end - - correct_runs += 1 - - # Only time correct solutions - trial_fn = () -> solve_instance(problem_type, instance, solver) - b_instance = @benchmarkable $trial_fn() - # Use pre-tuned parameters - b_instance.params = b.params - timing_result = run(b_instance, samples=1) - - # Collect statistics - push!(all_times, median(timing_result).time / 1e9) # Convert to seconds - push!(all_memory, median(timing_result).memory) - push!(all_gc_times, median(timing_result).gctime / 1e9) - successful_runs += 1 - - if i % 10 == 0 || i == length(instances) - @info " Completed $i/$(length(instances)) instances ($(correct_runs) correct, $(incorrect_runs) incorrect)" - end - - catch e - @warn " Instance $i failed with error: $e" - failed_runs += 1 - end - end - - if isempty(all_times) - return create_failed_result(config, "all_instances_failed: $failed_runs errors, $incorrect_runs incorrect", - failed_runs + incorrect_runs) - end - - # Calculate statistics - return Dict( - "config" => config, - "status" => "success", - "instances_tested" => length(instances), - "successful_runs" => successful_runs, - "correct_runs" => correct_runs, - "incorrect_runs" => incorrect_runs, - "failed_runs" => failed_runs, - "accuracy_rate" => correct_runs / length(instances), - - # Timing statistics (only for correct solutions) - "mean_time" => mean(all_times), - "median_time" => median(all_times), - "min_time" => minimum(all_times), - "max_time" => maximum(all_times), - "std_time" => std(all_times), - "q25_time" => quantile(all_times, 0.25), - "q75_time" => quantile(all_times, 0.75), - - # Memory statistics - "mean_memory" => mean(all_memory), - "median_memory" => median(all_memory), - - # GC statistics - "mean_gc_time" => mean(all_gc_times), - "total_gc_time" => sum(all_gc_times), - "gc_percentage" => 100 * sum(all_gc_times) / sum(all_times) - ) - - catch e - @error " Benchmark failed: $e" - return create_failed_result(config, "benchmark_failed: $e") - end -end - - - - -"""Create result dict for failed benchmark""" -function create_failed_result(config, reason, failed_runs=0) - return Dict( - "config" => config, - "status" => "failed", - "reason" => reason, - "samples" => 0, - "failed_runs" => failed_runs - ) -end - - - -"""Format bytes in human-readable way""" -function format_bytes(bytes::Number) - if bytes < 1024 - return "$(Int(bytes))B" - elseif bytes < 1024^2 - return "$(round(bytes/1024, digits=1))KB" - elseif bytes < 1024^3 - return "$(round(bytes/1024^2, digits=1))MB" - else - return "$(round(bytes/1024^3, digits=1))GB" - end -end - -""" - benchmark_backend(problem_type::Type{<:AbstractBenchmarkProblem}, - configs::Vector{<:AbstractProblemConfig}; - solver=nothing, - dataset_per_config=50, - force_regenerate_datasets=false) - -Run complete benchmark suite using pre-generated datasets. -""" -function benchmark_backend(problem_type::Type{<:AbstractBenchmarkProblem}, - configs::Vector{<:AbstractProblemConfig}; - solver=nothing, - dataset_per_config::Int=50, - force_regenerate_datasets::Bool=false) - - problem_name = lowercase(string(problem_type)[1:end-7]) # Remove "Problem" suffix - - # Generate datasets - @info "Preparing $problem_name datasets..." - generate_datasets(problem_type; - configs=configs, - per_config=dataset_per_config, - include_solution=true, - force_regenerate=force_regenerate_datasets) - - # Run benchmarks using dataset instances - @info "Running $problem_name benchmarks..." - results = benchmark_problem(problem_type; - configs=configs, - solver=solver) - - # Save results with enhanced information - results_path = joinpath(resolve_data_dir(), "$(problem_name)_benchmark_results.json") - open(results_path, "w") do io - # Get solver info for results - actual_solver = isnothing(solver) ? default_solver(problem_type) : solver - solver_info = solver_name(actual_solver) - - enhanced_results = Dict( - "benchmark_info" => Dict( - "problem_type" => string(problem_type), - "solver" => solver_info, - "timestamp" => string(now()), - "dataset_per_config" => dataset_per_config, - "force_regenerate_datasets" => force_regenerate_datasets, - "julia_version" => string(VERSION) - ), - "results" => results - ) - JSON3.pretty(io, enhanced_results) - end - @info "Benchmark results saved to: $results_path" - - # Print summary - print_benchmark_summary(results) - - return results -end - -""" - run_full_benchmark(problem_type::Type{<:AbstractBenchmarkProblem}, - config_tuples::Vector{<:Tuple}; - solver=nothing, kwargs...) - -High-level interface: run benchmark with tuple configs (e.g., [(10,10), (12,12)]). -""" -function run_full_benchmark(problem_type::Type{<:AbstractBenchmarkProblem}, - config_tuples::Vector{<:Tuple}; - solver=nothing, kwargs...) - configs = create_configs_from_tuples(problem_type, config_tuples) - return benchmark_backend(problem_type, configs; solver=solver, kwargs...) -end - -""" - run_full_benchmark(problem_type::Type{<:AbstractBenchmarkProblem}; - solver=nothing, force_regenerate_datasets=false, kwargs...) - -High-level interface: run benchmark with optional configs (uses defaults if nothing). -""" -function run_full_benchmark(problem_type::Type{<:AbstractBenchmarkProblem}; - solver=nothing, - kwargs...) - configs = default_configs(problem_type) - return benchmark_backend(problem_type, configs; solver=solver, kwargs...) -end - -"""Create config objects from tuples""" -function create_configs_from_tuples(problem_type::Type{<:AbstractBenchmarkProblem}, tuples) - return [config(problem_type, tuple) for tuple in tuples] -end - -""" - regenerate_datasets(problem_type::Type{<:AbstractBenchmarkProblem}; - configs=nothing, per_config=100, include_solution=true) - -Force regeneration of datasets for a problem type, ignoring existing cache. -""" -function regenerate_datasets(problem_type::Type{<:AbstractBenchmarkProblem}; - configs::Union{Vector{<:AbstractProblemConfig}, Nothing}=nothing, - per_config::Int=100, - include_solution::Bool=true) - actual_configs = isnothing(configs) ? default_configs(problem_type) : configs - return generate_datasets(problem_type; - configs=actual_configs, - per_config=per_config, - include_solution=include_solution, - force_regenerate=true) -end - -"""Print a nice summary of benchmark results""" -function print_benchmark_summary(results) - println("\n" * repeat("=", 60)) - println("BENCHMARK SUMMARY") - println(repeat("=", 60)) - - successful = filter(r -> r["status"] == "success", results) - failed = filter(r -> r["status"] == "failed", results) - - if !isempty(successful) - println("Successful benchmarks: $(length(successful))") - println(repeat("+", 90)) - println("| Config | Median Time | Memory | Instances | Correct | Accuracy | Timed |") - println(repeat("-", 90)) - - for result in successful - config = result["config"] - median_time = result["median_time"] - median_memory = result["median_memory"] - instances = get(result, "instances_tested", 0) - correct = get(result, "correct_runs", 0) - accuracy = get(result, "accuracy_rate", 0.0) - timed = get(result, "successful_runs", 0) - - config_str = "$(config.m)x$(config.n)" - time_str = format_time(median_time) - memory_str = format_bytes(median_memory) - accuracy_str = "$(round(accuracy * 100, digits=1))%" - - println("| $(rpad(config_str, 11)) | $(rpad(time_str, 11)) | $(rpad(memory_str, 11)) | $(rpad(string(instances), 9)) | $(rpad(string(correct), 7)) | $(rpad(accuracy_str, 8)) | $(rpad(string(timed), 5)) |") - end - println(repeat("+", 90)) - end - - if !isempty(failed) - println("\nFailed benchmarks: $(length(failed))") - for result in failed - config = result["config"] - reason = result["reason"] - println(" - $(config): $reason") - end - end - println(repeat("=", 60)) -end - -"""Format time in human-readable way""" -function format_time(seconds::Float64) - if seconds < 1e-6 - return "$(round(seconds * 1e9, digits=1))ns" - elseif seconds < 1e-3 - return "$(round(seconds * 1e6, digits=1))us" - elseif seconds < 1 - return "$(round(seconds * 1e3, digits=1))ms" - elseif seconds < 60 - return "$(round(seconds, digits=2))s" - else - minutes = floor(seconds / 60) - secs = seconds - minutes * 60 - return "$(Int(minutes))m$(round(secs, digits=1))s" - end -end - -""" - list_available_solvers(problem_type::Type{<:AbstractBenchmarkProblem}) - -Display available solvers for a problem type. -""" -function list_available_solvers(problem_type::Type{<:AbstractBenchmarkProblem}) - solvers = available_solvers(problem_type) - default = default_solver(problem_type) - default_name = solver_name(default) - - println("Available solvers for $(string(problem_type)):") - for solver in solvers - name = solver_name(solver) - is_default = name == default_name ? " (default)" : "" - println(" • $name$is_default") - end -end - -""" - run_solver_comparison(problem_type::Type{<:AbstractBenchmarkProblem}, config_tuples; kwargs...) - -Compare all available solvers on the same problem configurations. -""" -function run_solver_comparison(problem_type::Type{<:AbstractBenchmarkProblem}, - config_tuples::Vector{<:Tuple}; kwargs...) - solvers = available_solvers(problem_type) - results = Dict() - - for solver in solvers - solver_name_str = solver_name(solver) - @info "Running benchmark with solver: $solver_name_str" - results[solver_name_str] = run_full_benchmark(problem_type, config_tuples; solver=solver, kwargs...) - end - - return results -end - -""" - save_solver_comparison(problem_type::Type{<:AbstractBenchmarkProblem}, solver_results::Dict; note="") - -Persist solver comparison results to JSON under the benchmark data directory. -Includes benchmark metadata (timestamp, Julia version) and optional note. -""" -function save_solver_comparison(problem_type::Type{<:AbstractBenchmarkProblem}, solver_results::Dict; note::AbstractString="") - if isempty(solver_results) - @warn "No solver results to save. Skipping." - return nothing - end - - # Determine problem name and output path - problem_name = lowercase(string(problem_type)[1:end-7]) # Remove "Problem" suffix - results_path = joinpath(resolve_data_dir(), "$(problem_name)_solver_comparison.json") - - # Extract common configs from the first solver - first_results = first(values(solver_results)) - configs = [result["config"] for result in first_results if haskey(result, "config")] - - # Build a normalized structure for saving - payload = Dict( - "benchmark_info" => Dict( - "problem_type" => string(problem_type), - "timestamp" => string(now()), - "julia_version" => string(VERSION), - "note" => note, - "configs" => [Dict("m"=>c.m, "n"=>c.n) for c in configs] - ), - "solvers" => Dict( - (name => solver_results[name]) for name in keys(solver_results) - ) - ) - - open(results_path, "w") do io - JSON3.pretty(io, payload) - end - @info "Solver comparison saved to: $results_path" - return results_path -end - -""" - print_solver_comparison_summary(solver_results::Dict) - -Print a comprehensive comparison table of all solver results. -""" -function print_solver_comparison_summary(solver_results::Dict) - if isempty(solver_results) - println("No solver results to compare!") - return - end - - println("\n" * repeat("=", 80)) - println("SOLVER COMPARISON SUMMARY") - println(repeat("=", 80)) - - # Get all configs from first solver (assuming all solvers tested same configs) - first_results = first(values(solver_results)) - configs = [result["config"] for result in first_results if result["status"] == "success"] - - if isempty(configs) - println("No successful results to compare!") - return - end - - # Print header - solver_names = collect(keys(solver_results)) - header = "| Config " - for solver in solver_names - header *= "| $(rpad(solver, 15)) " - end - header *= "|" - - println(repeat("+", length(header))) - println(header) - println(repeat("-", length(header))) - - # Print results for each config - for config in configs - config_str = "$(config.m)x$(config.n)" - row = "| $(rpad(config_str, 8)) " - - for solver in solver_names - solver_result = solver_results[solver] - # Find result for this config - matching_result = nothing - for result in solver_result - if result["config"] == config && result["status"] == "success" - matching_result = result - break - end - end - - if isnothing(matching_result) - time_str = "FAILED" - else - median_time = matching_result["median_time"] - time_str = format_time(median_time) - end - row *= "| $(rpad(time_str, 15)) " - end - row *= "|" - println(row) - end - println(repeat("+", length(header))) - - # Print winner analysis with accuracy information - println("\nPerformance winners by configuration (with accuracy):") - for config in configs - config_str = "$(config.m)x$(config.n)" - times = Dict{String, Float64}() - accuracies = Dict{String, Float64}() - - for solver in solver_names - solver_result = solver_results[solver] - for result in solver_result - if result["config"] == config && result["status"] == "success" - times[solver] = result["median_time"] - accuracies[solver] = get(result, "accuracy_rate", 0.0) - break - end - end - end - - if !isempty(times) - # Only consider solvers with perfect accuracy for performance comparison - perfect_solvers = filter(kv -> accuracies[kv[1]] == 1.0, times) - - if !isempty(perfect_solvers) - winner = argmin(perfect_solvers) - winner_time = perfect_solvers[winner] - println(" - $config_str: winner $winner ($(format_time(winner_time)), 100% accuracy)") - - # Show relative performance for perfect solvers - sorted_perfect = sort(collect(perfect_solvers), by=x->x[2]) - if length(sorted_perfect) > 1 - for (i, (solver, time)) in enumerate(sorted_perfect[2:end]) - speedup = time / winner_time - println(" - $solver: $(round(speedup, digits=2))x slower (100% accuracy)") - end - end - - # Show imperfect solvers separately - imperfect_solvers = filter(kv -> accuracies[kv[1]] < 1.0, times) - for (solver, time) in imperfect_solvers - accuracy_pct = round(accuracies[solver] * 100, digits=1) - println(" - $solver: $(format_time(time)) ($(accuracy_pct)% accuracy)") - end - else - println(" - $config_str: no solver achieved 100% accuracy") - for (solver, time) in sort(collect(times), by=x->x[2]) - accuracy_pct = round(accuracies[solver] * 100, digits=1) - println(" $(solver): $(format_time(time)) ($(accuracy_pct)% accuracy)") - end - end - end - end - println(repeat("=", 80)) -end - -""" - save_benchmark_results(problem_type::Type{<:AbstractBenchmarkProblem}, results; solver=nothing, dataset_per_config=nothing, force_regenerate_datasets=nothing) - -Persist single-solver benchmark results to JSON, mirroring benchmark_backend output. -""" -function save_benchmark_results(problem_type::Type{<:AbstractBenchmarkProblem}, results; - solver=nothing, - dataset_per_config=nothing, - force_regenerate_datasets=nothing) - problem_name = lowercase(string(problem_type)[1:end-7]) - results_path = joinpath(resolve_data_dir(), "$(problem_name)_benchmark_results.json") - - actual_solver = isnothing(solver) ? default_solver(problem_type) : solver - solver_info = solver_name(actual_solver) - - payload = Dict( - "benchmark_info" => Dict( - "problem_type" => string(problem_type), - "solver" => solver_info, - "timestamp" => string(now()), - "dataset_per_config" => dataset_per_config, - "force_regenerate_datasets" => force_regenerate_datasets, - "julia_version" => string(VERSION) - ), - "results" => results - ) - - open(results_path, "w") do io - JSON3.pretty(io, payload) - end - @info "Benchmark results saved to: $results_path" - return results_path -end - -""" - compare_solver_results(name1, results1, name2, results2) - -Compare two specific solver results side by side with accuracy information. -""" -function compare_solver_results(name1::String, results1, name2::String, results2) - println(repeat("+", 98)) - println("| $(rpad("Config", 8)) | $(rpad("$name1 Time", 13)) | $(rpad("$name1 Acc", 10)) | $(rpad("$name2 Time", 13)) | $(rpad("$name2 Acc", 10)) | $(rpad("Winner", 18)) |") - println(repeat("-", 98)) - - # Get successful results from both solvers - success1 = filter(r -> r["status"] == "success", results1) - success2 = filter(r -> r["status"] == "success", results2) - - # Match configs - for r1 in success1 - config = r1["config"] - config_str = "$(config.m)x$(config.n)" - - # Find matching config in results2 - r2 = findfirst(r -> r["config"] == config, success2) - if !isnothing(r2) - r2 = success2[r2] - - time1 = r1["median_time"] - time2 = r2["median_time"] - time1_str = format_time(time1) - time2_str = format_time(time2) - - acc1 = get(r1, "accuracy_rate", 0.0) - acc2 = get(r2, "accuracy_rate", 0.0) - acc1_str = "$(round(acc1 * 100, digits=1))%" - acc2_str = "$(round(acc2 * 100, digits=1))%" - - # Determine winner considering both speed and accuracy - if acc1 == 1.0 && acc2 == 1.0 - # Both perfect accuracy - compare speed - if time1 < time2 - winner_str = "$name1 ($(round(time2/time1, digits=2))x faster)" - elseif time2 < time1 - winner_str = "$name2 ($(round(time1/time2, digits=2))x faster)" - else - winner_str = "Tie" - end - elseif acc1 == 1.0 && acc2 < 1.0 - winner_str = "$name1 (100% acc)" - elseif acc1 < 1.0 && acc2 == 1.0 - winner_str = "$name2 (100% acc)" - else - # Neither has perfect accuracy - compare accuracy first, then speed - if acc1 > acc2 - winner_str = "$name1 (better acc)" - elseif acc2 > acc1 - winner_str = "$name2 (better acc)" - else - # Same accuracy - compare speed - if time1 < time2 - winner_str = "$name1 (faster)" - elseif time2 < time1 - winner_str = "$name2 (faster)" - else - winner_str = "Tie" - end - end - end - - println("| $(rpad(config_str, 8)) | $(rpad(time1_str, 13)) | $(rpad(acc1_str, 10)) | $(rpad(time2_str, 13)) | $(rpad(acc2_str, 10)) | $(rpad(winner_str, 18)) |") - else - acc1 = get(r1, "accuracy_rate", 0.0) - acc1_str = "$(round(acc1 * 100, digits=1))%" - println("| $(rpad(config_str, 8)) | $(rpad(format_time(r1["median_time"]), 13)) | $(rpad(acc1_str, 10)) | $(rpad("FAILED", 13)) | $(rpad("N/A", 10)) | $(rpad(name1, 18)) |") - end - end - - # Check for configs only in results2 - for r2 in success2 - config = r2["config"] - if !any(r -> r["config"] == config, success1) - config_str = "$(config.m)x$(config.n)" - acc2 = get(r2, "accuracy_rate", 0.0) - acc2_str = "$(round(acc2 * 100, digits=1))%" - println("| $(rpad(config_str, 8)) | $(rpad("FAILED", 13)) | $(rpad("N/A", 10)) | $(rpad(format_time(r2["median_time"]), 13)) | $(rpad(acc2_str, 10)) | $(rpad(name2, 18)) |") - end - end - println(repeat("+", 98)) -end diff --git a/benchmarks/src/utils.jl b/benchmarks/src/utils.jl index b534a98..eef0714 100644 --- a/benchmarks/src/utils.jl +++ b/benchmarks/src/utils.jl @@ -1,9 +1,3 @@ -# Generic benchmark utilities - -""" - resolve_data_dir(parts...) -Resolve data directory under benchmark package root. -""" function resolve_data_dir(parts::AbstractString...) base = normpath(joinpath(@__DIR__, "..")) dir = joinpath(base, "data", parts...) @@ -11,19 +5,11 @@ function resolve_data_dir(parts::AbstractString...) return dir end -""" - write_jsonl(io, obj) -Write one JSON object per line (JSONL format). -""" function write_jsonl(io::IO, obj) JSON3.write(io, obj) write(io, '\n') end -""" - read_jsonl(path) -Read JSONL file and return vector of parsed objects. -""" function read_jsonl(path::AbstractString) objs = Any[] open(path, "r") do io @@ -33,57 +19,4 @@ function read_jsonl(path::AbstractString) end end return objs -end - -""" - build_datasets(; outdir, generator, filename_fn, configs, per_config=100, rng) -Build separate dataset files for each config (one file per config). -""" -function build_datasets(; - outdir::AbstractString, - generator, - filename_fn, - configs, - per_config::Int=100, - rng::AbstractRNG=Random.GLOBAL_RNG) - isdir(outdir) || mkpath(outdir) - paths = String[] - for cfg in configs - path = joinpath(outdir, filename_fn(cfg)) - open(path, "w") do io - for _ in 1:per_config - obj = generator(cfg; rng=rng) - write_jsonl(io, obj) - end - end - @info "Generated dataset: $path" - push!(paths, path) - end - @info "Generated $(length(paths)) dataset files in: $outdir" - return paths -end - -""" - benchmark_configs(configs; trial_fn, samples_per_config=5) -Generic benchmark runner for different configurations. -""" -function benchmark_configs(configs; trial_fn, samples_per_config::Int=5) - results = [] - for cfg in configs - @info "Benchmarking config: $cfg" - times = Float64[] - for _ in 1:samples_per_config - t = @benchmarkable $trial_fn($cfg) - tune!(t) - run(t) - push!(times, t.times) - end - push!(results, Dict( - "config" => cfg, - "times" => times, - "mean_time" => sum(times) / length(times), - "min_time" => minimum(times) - )) - end - return results end \ No newline at end of file diff --git a/benchmarks/third-party/CnC b/benchmarks/third-party/CnC new file mode 160000 index 0000000..705b60c --- /dev/null +++ b/benchmarks/third-party/CnC @@ -0,0 +1 @@ +Subproject commit 705b60c6491ef2b61988b3ce6ac674be1b90571d diff --git a/benchmarks/third-party/abc b/benchmarks/third-party/abc new file mode 160000 index 0000000..18f6464 --- /dev/null +++ b/benchmarks/third-party/abc @@ -0,0 +1 @@ +Subproject commit 18f6464ec765c2ca14ede4f7a17f09011b0a0600 diff --git a/benchmarks/third-party/aiger b/benchmarks/third-party/aiger new file mode 160000 index 0000000..3290e7a --- /dev/null +++ b/benchmarks/third-party/aiger @@ -0,0 +1 @@ +Subproject commit 3290e7ac1ee0217ea3ca91d58e608642096b03c4 diff --git a/benchmarks/third-party/glucose b/benchmarks/third-party/glucose new file mode 160000 index 0000000..7f887ab --- /dev/null +++ b/benchmarks/third-party/glucose @@ -0,0 +1 @@ +Subproject commit 7f887abba7cf13636a5ac2d28653668a20a91b25 diff --git a/benchmarks/third-party/x-sat b/benchmarks/third-party/x-sat new file mode 160000 index 0000000..1787935 --- /dev/null +++ b/benchmarks/third-party/x-sat @@ -0,0 +1 @@ +Subproject commit 17879351ead8154ca0758026cd381c5b34ce8e77 diff --git a/examples/plot.jl b/examples/plot.jl new file mode 100644 index 0000000..744b86d --- /dev/null +++ b/examples/plot.jl @@ -0,0 +1,30 @@ +using CairoMakie + +configs = ["10x10", "12x12", "14x14", "16x16", "18x18", "20x20", "22x22"] +sizes = [20, 24, 28, 32, 36, 40, 44] + +bi = [0.0671, 0.1093, 0.1547, 0.6152, 1.03, 2.67, 20.30] +gurobi = [0.0675, 0.7289, 1.63, 7.6, NaN, NaN, NaN] +xsat = [0.0614, 0.1026, 0.1876, 0.3284, 1.26, 2.51, 11.32] + + +f = Figure(resolution = (700, 450)) +ax = Axis(f[1, 1]; + xlabel = "Semiprime length (bits)", + ylabel = "Runtime (seconds)", + title = "Runtime Comparison", + yscale = log10, +) + +lines!(ax, sizes, bi; label = "BooleanInference", linewidth = 2) +lines!(ax, sizes, gurobi; label = "IP-Gurobi", linewidth = 2) +lines!(ax, sizes, xsat; label = "X-SAT", linewidth = 2, color = :red) + +scatter!(ax, sizes, bi) +scatter!(ax, sizes, gurobi) +scatter!(ax, sizes, xsat; color = :red) + +axislegend(ax; position = :lt) + +save("runtime_comparison.png", f) +f \ No newline at end of file diff --git a/src/BooleanInference.jl b/src/BooleanInference.jl index c8e1133..b6c09ee 100644 --- a/src/BooleanInference.jl +++ b/src/BooleanInference.jl @@ -3,7 +3,7 @@ module BooleanInference using TropicalNumbers using SparseArrays using OptimalBranchingCore -using OptimalBranchingCore: AbstractProblem, select_variables, reduce_problem, _vec2int, candidate_clauses +using OptimalBranchingCore: AbstractProblem, select_variables, reduce_problem, _vec2int, candidate_clauses, Clause using OptimalBranchingCore.BitBasis using GenericTensorNetworks using GenericTensorNetworks.OMEinsum @@ -11,71 +11,86 @@ import ProblemReductions import ProblemReductions: CircuitSAT, Circuit, Factoring, reduceto, Satisfiability using DataStructures using DataStructures: PriorityQueue +using Statistics: median + +include("core/types.jl") +include("core/static.jl") +include("core/domain.jl") +include("core/stats.jl") +include("core/workspace.jl") +include("core/region.jl") +include("core/problem.jl") + +include("utils/utils.jl") + +include("algorithms/knn.jl") +include("algorithms/contraction.jl") +include("algorithms/propagate.jl") + +include("branching/measure.jl") +include("branching/selector.jl") +include("branching/branchtable.jl") +include("branching/greedymerge.jl") +include("branching/branch.jl") -# using GenericTensorNetworks: ∧, ∨, ¬ - -include("problems.jl") -include("region.jl") -include("measure.jl") -include("utils.jl") -include("knn.jl") -include("selector.jl") -include("contraction.jl") -include("branchtable.jl") -include("branch.jl") -include("propagate.jl") include("interface.jl") -# Export core types export Variable, EdgeRef, BoolTensor, TNStatic, DomainMask, TNProblem export DM_BOTH, DM_0, DM_1 -export Region, RegionContraction, RegionCacheEntry, RegionCacheState, RegionCache +export Region, RegionCacheEntry, RegionCacheState export DynamicWorkspace -# Export domain mask functions export is_fixed, has0, has1, init_doms, get_var_value -# Export problem setup functions export setup_problem, setup_from_tensor_network, setup_from_cnf, setup_from_circuit, setup_from_sat -# Export problem state functions export is_solved, cache_branch_solution!, reset_last_branch_problem!, has_last_branch_problem, last_branch_problem -# Export solving functions export solve, solve_sat_problem, solve_sat_with_assignments, solve_factoring +export solve_sat_problem_lookahead, solve_sat_problem_hybrid +export solve_sat_problem_region, solve_sat_problem_region_lookahead +export solve_sat_problem_cdcl -# Export measure types export NumUnfixedVars -# Export selector types -export LeastOccurrenceSelector, AbstractSelector +export MostOccurrenceSelector, MinGammaSelector, AbstractSelector +export RegionAwareSelector, PropagationAwareSelector +export RegionQualityMetrics, evaluate_region_quality, score_region_quality +export evaluate_propagation_power, compute_propagation_score, evaluate_propagation_for_branch -# Export table solver types export TNContractionSolver, AbstractTableSolver -# Export contraction functions export contract_region, contract_tensors, slicing, tensor_unwrapping -# Export propagation functions export propagate, get_active_tensors, build_tensor_masks export TensorMasks, PropagationBuffers +export look_ahead_propagation, look_ahead_score, LookAheadResult -# Export region management functions -export cache_region!, get_cached_region_entry, get_cached_region, cache_region_contraction! -export get_cached_region_contraction, clear_region_cache!, clear_all_region_caches! +export propagate_incremental, propagate_from_assignment, propagate_from_assignments +export propagate_adaptive, should_use_incremental + +export propagate_cdcl, propagate_after_assignment_cdcl, propagate_with_state! +export CDCLState, LearnedClause, TrailEntry +export solve_cdcl, get_solution, print_cdcl_stats +export propagate_bcp!, analyze_conflict!, make_decision!, backtrack! + +export cache_region!, get_cached_region, clear_all_region_caches! -# Export k-neighboring functions export k_neighboring, KNNWorkspace -# Export utility functions export get_unfixed_vars, count_unfixed, bits_to_int -# Export branching statistics functions export get_branching_stats, reset_branching_stats!, print_branching_stats -# Export branching table functions -export separate_fixed_free_boundary, construct_boundary_config, construct_inner_config -export extract_inner_configs, combine_configs, get_region_contraction, slice_region_contraction -export handle_no_boundary_case +export BranchingStats, DetailedStats +export record_depth!, record_branch!, record_unsat_leaf!, record_solved_leaf!, record_skipped_subproblem! +export record_propagation!, record_propagation_fixpoint!, record_domain_reduction!, record_early_unsat! +export record_branching_time!, record_contraction_time!, record_filtering_time!, record_cache_hit!, record_cache_miss! +export record_variable_selection! +export print_stats_summary +export reset! + +export extract_inner_configs, combine_configs, slice_region_contraction +export handle_no_boundary_case_unfixed end diff --git a/src/contraction.jl b/src/algorithms/contraction.jl similarity index 83% rename from src/contraction.jl rename to src/algorithms/contraction.jl index 8535772..f455dbc 100644 --- a/src/contraction.jl +++ b/src/algorithms/contraction.jl @@ -63,45 +63,53 @@ function contract_tensors(tensors::Vector{Vector{T}}, ixs::Vector{Vector{Int}}, end -function slicing(tensor::Vector{T}, doms::Vector{DomainMask}, axis_vars::Vector{Int}) where T +@inline function slicing(tensor::Vector{T}, doms::Vector{DomainMask}, axis_vars::Vector{Int}) where T # Number of axes (Boolean => size 2 per axis) k = trailing_zeros(length(tensor)) # log2 - fixed_axes = Int[]; fixed_vals = Int[]; free_axes = Int[] - - for axis in 1:k # each variable + + # Pre-allocate with size hints + fixed_axes = Int[] + fixed_vals = Int[] + free_axes = Int[] + sizehint!(fixed_axes, k) + sizehint!(fixed_vals, k) + sizehint!(free_axes, k) + + @inbounds for axis in 1:k # each variable var_id = axis_vars[axis] dom_mask = doms[var_id] if is_fixed(dom_mask) push!(fixed_axes, axis) - push!(fixed_vals, has1(dom_mask) ? 1 : 0) + push!(fixed_vals, ifelse(has1(dom_mask), 1, 0)) else push!(free_axes, axis) end end - + n_free = length(free_axes) out_len = 1 << n_free # 2^n_free out = Vector{T}(undef, out_len) - + + # Precompute fixed part of the index + fixed_idx = 0 + @inbounds for (axis, val) in zip(fixed_axes, fixed_vals) + fixed_idx |= (val << (axis-1)) + end + # Iterate over all assignments to free variables @inbounds for free_idx in 0:(out_len-1) - # Build full index - full_idx = 0 - + # Build full index starting from fixed part + full_idx = fixed_idx + # Set free variable bits for (i, axis) in enumerate(free_axes) bit = (free_idx >> (i-1)) & 0x1 # get the i-th bit of free_idx full_idx |= (bit << (axis-1)) # set free-idx in the full_idx end - - # Set fixed variable bits - for (axis, val) in zip(fixed_axes, fixed_vals) - full_idx |= (val << (axis-1)) - end - + out[free_idx+1] = tensor[full_idx+1] end - + return out end diff --git a/src/knn.jl b/src/algorithms/knn.jl similarity index 97% rename from src/knn.jl rename to src/algorithms/knn.jl index 5c9f5c2..d81e0b5 100644 --- a/src/knn.jl +++ b/src/algorithms/knn.jl @@ -57,7 +57,7 @@ function expand_one_var!( for e in tn.t2v[tensor_id] var_id = e.var # Only expand to unfixed variables - if !_check_seen_var(ws, var_id) && !is_fixed(doms[var_id]) + if !_check_seen_var(ws, var_id) _mark_var!(ws, var_id) push!(ws.collected_vars, var_id) push!(ws.next_frontier, var_id) diff --git a/src/algorithms/propagate.jl b/src/algorithms/propagate.jl new file mode 100644 index 0000000..4935eb1 --- /dev/null +++ b/src/algorithms/propagate.jl @@ -0,0 +1,270 @@ +# function get_active_tensors(static::TNStatic, doms::Vector{DomainMask}) +# active = Int[] +# sizehint!(active, length(static.tensors)) +# @inbounds for (tid, tensor) in enumerate(static.tensors) +# # Check if any variable in this tensor is unfixed +# has_unfixed = false +# for var_id in tensor.var_axes +# dm_bits = doms[var_id].bits +# # Unfixed if domain allows both 0 and 1, or is not yet determined +# has_unfixed |= dm_bits == 0x03 +# end +# has_unfixed && push!(active, tid) +# end +# return active +# end + +# Backward compatible interface: propagate all tensors when changed_vars not specified +function propagate(static::TNStatic, doms::Vector{DomainMask}, ws::Union{Nothing, DynamicWorkspace}=nothing) + # When no changed_vars specified, propagate all tensors + all_vars = collect(1:length(doms)) + return propagate(static, doms, all_vars, ws) +end + +function propagate(static::TNStatic, doms::Vector{DomainMask}, changed_vars::Vector{Int}, ws::Union{Nothing, DynamicWorkspace}=nothing) + isempty(changed_vars) && return doms + working_doms = copy(doms) + + stats = !isnothing(ws) ? ws.branch_stats : nothing + has_detailed = !isnothing(stats) && !isnothing(stats.detailed) + + # Record propagation time (nanosecond precision) + propagation_start_time = has_detailed ? time_ns() : 0 + + # Build the propagation queue with tensors affected by the changed variables + tensor_queue = Int[] + in_queue = falses(length(static.tensors)) + @inbounds for var_id in changed_vars + @inbounds for tid in static.v2t[var_id] + if !in_queue[tid] + push!(tensor_queue, tid) + in_queue[tid] = true + end + end + end + + # Reuse propagation buffers to avoid reallocating BitVectors + buffers = if !isnothing(ws) && !isnothing(ws.prop_buffers) + ws.prop_buffers + else + buf = PropagationBuffers(static) + !isnothing(ws) && (ws.prop_buffers = buf) + buf + end + + # Count initial domain assignments + initial_domain_count = has_detailed ? sum(count_ones(d.bits & 0x03) for d in doms) : 0 + + queue_pos = 1 + iteration_count = 0 + while queue_pos <= length(tensor_queue) + iteration_count += 1 + tensor_idx = tensor_queue[queue_pos] + queue_pos += 1 + in_queue[tensor_idx] = false + + tensor = static.tensors[tensor_idx] + # Use the precomputed masks directly without cache lookups + masks = static.tensor_to_masks[tensor_idx] + + # Check constraint and propagate + if !propagate_tensor!(working_doms, tensor, masks, static.v2t, + tensor_queue, in_queue, buffers) + # Contradiction detected + if has_detailed + record_early_unsat!(stats) + end + return fill(DomainMask(0x00), length(working_doms)) + end + end + + # Check whether we reached a fixed point (iteration count) + if has_detailed && iteration_count == 1 && queue_pos > length(tensor_queue) + record_propagation_fixpoint!(stats) + end + + # Record domain reduction counts and elapsed time + if has_detailed + propagation_time = (time_ns() - propagation_start_time) / 1e9 # Convert to seconds + record_propagation!(stats, propagation_time) + + final_domain_count = sum(count_ones(d.bits & 0x03) for d in working_doms) + domain_reduction = initial_domain_count - final_domain_count + if domain_reduction > 0 + record_domain_reduction!(stats, domain_reduction) + end + end + + return working_doms +end + +# Propagate constraints from a single tensor. Returns false if contradiction detected. +function propagate_tensor!(working_doms::Vector{DomainMask}, tensor::BoolTensor, masks::TensorMasks, v2t::Vector{Vector{Int}}, tensor_queue::Vector{Int}, in_queue::BitVector, buffers::PropagationBuffers) + # Step 1: Compute feasible configurations given current domains + compute_feasible_configs!(buffers.feasible, working_doms, tensor, masks) || return false + n_feasible = count(buffers.feasible) + + # Step 2: Different propagation strategies based on number of feasible configs + if n_feasible == 1 + # Only one valid config -> fix all variables to that config + return propagate_unit_constraint!(working_doms, tensor, buffers.feasible, v2t, tensor_queue, in_queue) + else + # Multiple configs -> prune unsupported values + return propagate_support_pruning!(working_doms, tensor, masks, buffers.feasible, buffers.temp, v2t, tensor_queue, in_queue) + end +end + +# Compute which tensor configurations are feasible given current variable domains. +# Returns false if no feasible configuration exists. +@inline function compute_feasible_configs!(feasible::BitVector, working_doms::Vector{DomainMask}, tensor::BoolTensor, masks::TensorMasks) + n_cfg = length(masks.sat) + + copyto!(feasible, 1, masks.sat, 1, n_cfg) + + n_words = (n_cfg + 63) >> 6 + feas_chunks = feasible.chunks + + # Clear stale bits beyond n_cfg in the final UInt64 to avoid reuse artefacts + bits_in_last_word = n_cfg & 63 # n_cfg % 64 + if bits_in_last_word != 0 + # Keep only the lowest bits_in_last_word bits in the last word + mask_val = (UInt64(1) << bits_in_last_word) - 1 + feas_chunks[n_words] &= mask_val + end + + @inbounds for (axis, var_id) in enumerate(tensor.var_axes) + dm_bits = working_doms[var_id].bits + + # Skip unconstrained variables (both 0 and 1 are allowed) + dm_bits == 0x03 && continue + # Contradiction: variable has no feasible assignment + dm_bits == 0x00 && return false + + # Pick the corresponding mask (0 or 1) + mask_chunks = dm_bits == 0x01 ? masks.axis_masks0[axis].chunks : masks.axis_masks1[axis].chunks + + # Apply the mask and ensure a feasible configuration remains + # Optimization: merge computation and checking to avoid extra memory traffic + has_nonzero = false + for i in 1:n_words + feas_chunks[i] &= mask_chunks[i] + has_nonzero |= (feas_chunks[i] != 0) + end + + # No feasible configuration left, exit early + has_nonzero || return false + end + + return true +end + +# When tensor has exactly one feasible configuration, fix all its variables to that config. +@inline function propagate_unit_constraint!(working_doms::Vector{DomainMask}, tensor::BoolTensor, feasible::BitVector, v2t::Vector{Vector{Int}}, tensor_queue::Vector{Int}, in_queue::BitVector) + first_idx = findfirst(feasible) + config = first_idx - 1 + + @inbounds for (axis, var_id) in enumerate(tensor.var_axes) + bit_val = (config >> (axis - 1)) & 1 + # Branchless: avoid conditional for better performance + required_bits = ifelse(bit_val == 1, DM_1.bits, DM_0.bits) + + old_bits = working_doms[var_id].bits + new_bits = old_bits & required_bits + + # Check for contradiction + new_bits == 0x00 && return false + + # Update domain and enqueue affected tensors + if new_bits != old_bits + working_doms[var_id] = DomainMask(new_bits) + enqueue_affected_tensors!(tensor_queue, in_queue, v2t, var_id) + end + end + return true +end + +# Prune domain values that have no support in any feasible configuration. +@inline function propagate_support_pruning!(working_doms::Vector{DomainMask}, tensor::BoolTensor, masks::TensorMasks, feasible::BitVector, temp::BitVector, v2t::Vector{Vector{Int}}, tensor_queue::Vector{Int}, in_queue::BitVector) + n_cfg = length(masks.sat) + n_words = (n_cfg + 63) >> 6 + feas_chunks = feasible.chunks + + @inbounds for (axis, var_id) in enumerate(tensor.var_axes) + dm = working_doms[var_id] + dm_bits = dm.bits + + if (dm_bits & 0x01) != 0 + has_support_0 = false + mask_chunks = masks.axis_masks0[axis].chunks + for i in 1:n_words + if (feas_chunks[i] & mask_chunks[i]) != 0 + has_support_0 = true + break + end + end + if !has_support_0 + if !update_domain!(working_doms, var_id, dm_bits, DM_1.bits, + v2t, tensor_queue, in_queue) + return false + end + end + end + + if (dm_bits & 0x02) != 0 + has_support_1 = false + mask_chunks = masks.axis_masks1[axis].chunks + for i in 1:n_words + if (feas_chunks[i] & mask_chunks[i]) != 0 + has_support_1 = true + break + end + end + if !has_support_1 + if !update_domain!(working_doms, var_id, dm_bits, DM_0.bits, v2t, tensor_queue, in_queue) + return false + end + end + end + end + + return true +end + +# Add all tensors containing the variable to the propagation queue. +@inline function enqueue_affected_tensors!( + tensor_queue::Vector{Int}, + in_queue::BitVector, + v2t::Vector{Vector{Int}}, + var_id::Int +) + @inbounds for tensor_id in v2t[var_id] + if !in_queue[tensor_id] + push!(tensor_queue, tensor_id) + in_queue[tensor_id] = true + end + end +end + +# Update variable domain by intersecting with keep_mask. Returns false if contradiction. +@inline function update_domain!( + working_doms::Vector{DomainMask}, + var_id::Int, + old_bits::UInt8, + keep_bits::UInt8, + v2t::Vector{Vector{Int}}, + tensor_queue::Vector{Int}, + in_queue::BitVector +) + new_bits = old_bits & keep_bits + + # Check for contradiction + new_bits == 0x00 && return false + + # Update if changed + if new_bits != old_bits + working_doms[var_id] = DomainMask(new_bits) + enqueue_affected_tensors!(tensor_queue, in_queue, v2t, var_id) + end + + return true +end diff --git a/src/branch.jl b/src/branch.jl deleted file mode 100644 index f9246f6..0000000 --- a/src/branch.jl +++ /dev/null @@ -1,169 +0,0 @@ -function OptimalBranchingCore.branch_and_reduce( - problem::TNProblem, - config::OptimalBranchingCore.BranchingStrategy, - reducer::OptimalBranchingCore.AbstractReducer, - result_type::Type{TR}; - show_progress::Bool=false, - tag::Vector{Tuple{Int,Int}}=Tuple{Int,Int}[] -) where TR - # Update max depth - current_depth = length(tag) - if current_depth > problem.ws.max_depth - problem.ws.max_depth = current_depth - end - - # Step 1: Check if problem is solved (all variables fixed) - if is_solved(problem) - @debug "problem is solved" - cache_branch_solution!(problem) - return one(result_type) - end - - # Step 2: Try to reduce the problem - @assert reducer isa NoReducer - # reduced_problem, reduced_value = reduce_problem(result_type, problem, reducer) - # if reduced_problem !== problem - # # Problem was reduced, recurse on reduced problem - # return branch_and_reduce( - # reduced_problem, config, reducer, result_type; - # tag=tag, show_progress=show_progress - # ) * reduced_value - # end - - # Step 3: Select variables for branching - variables = OptimalBranchingCore.select_variables(problem, config.measure, config.selector) - - # Step 4: Compute branching table - tbl = OptimalBranchingCore.branching_table(problem, config.table_solver, variables) - # Step 5: Compute optimal branching rule - result = OptimalBranchingCore.optimal_branching_rule(tbl, variables, problem, config.measure, config.set_cover_solver) - - # @show result.γ - - # Step 6: Branch and recurse - clauses = OptimalBranchingCore.get_clauses(result) - @debug "A new branch-level search starts with $(length(clauses)) clauses: $(clauses)" - - # Record branching statistics - problem.ws.total_branches += 1 - problem.ws.total_subproblems += length(clauses) - - return sum(enumerate(clauses)) do (i, branch) - show_progress && (OptimalBranchingCore.print_sequence(stdout, tag); println(stdout)) - @debug "branch=$branch, n_unfixed=$(problem.n_unfixed)" - - # Apply branch to get subproblem - subproblem, local_value = OptimalBranchingCore.apply_branch(problem, branch, variables) - - @debug "local_value=$local_value, n_unfixed=$(subproblem.n_unfixed)" - - # If branch led to contradiction (UNSAT), skip this branch - if local_value == 0 || subproblem.n_unfixed == 0 && any(dm -> dm.bits == 0x00, subproblem.doms) - @debug "Returning zero: local_value=$local_value, n_unfixed=$(subproblem.n_unfixed), has_contradiction=$(any(dm -> dm.bits == 0x00, subproblem.doms))" - return zero(result_type) - end - - # Recursively solve subproblem - # Always update tag for depth tracking, only print when show_progress=true - new_tag = [tag..., (i, length(clauses))] - sub_result = OptimalBranchingCore.branch_and_reduce(subproblem, config, reducer, result_type;tag=new_tag, show_progress=show_progress) - - # Combine results - sub_result * result_type(local_value) - end -end - -function OptimalBranchingCore.optimal_branching_rule( - tbl::OptimalBranchingCore.BranchingTable, - variables::Vector{T}, - problem::TNProblem, - measure::OptimalBranchingCore.AbstractMeasure, - solver::OptimalBranchingCore.AbstractSetCoverSolver -) where T - candidates = OptimalBranchingCore.bit_clauses(tbl) - return OptimalBranchingCore.greedymerge(candidates, problem, variables, measure) -end - - -function OptimalBranchingCore.apply_branch( - problem::TNProblem, - clause::OptimalBranchingCore.Clause{INT}, - variables::Vector{Int} -) where {INT<:Integer} - new_doms = copy(problem.doms) - # Apply clause: fix variables according to mask and values - n_fixed = 0 - for i in 1:length(variables) - if OptimalBranchingCore.readbit(clause.mask, i) == 1 - # This variable is fixed by the clause - var_id = variables[i] - bit_val = OptimalBranchingCore.readbit(clause.val, i) - new_val = (bit_val == 1) ? DM_1 : DM_0 - - if !is_fixed(problem.doms[var_id]) - new_doms[var_id] = new_val - n_fixed += 1 - end - end - end - - @debug "apply_branch: Fixed $n_fixed variables" - - # Apply propagation (unit propagation) - propagated_doms = propagate(problem.static, new_doms) - - # Check for contradiction (all domains set to 0x00) - if any(dm -> dm.bits == 0x00, propagated_doms) - # UNSAT: contradiction detected during propagation - @debug "apply_branch: Contradiction detected during propagation" - return (TNProblem(problem.static, fill(DomainMask(0x00), length(propagated_doms)), 0, problem.ws), 0) - end - - # Count unfixed variables - new_n_unfixed = count_unfixed(propagated_doms) - - @debug "apply_branch: n_unfixed: $(problem.n_unfixed) -> $new_n_unfixed" - - # Safety check: problem must have gotten smaller OR we fixed at least one variable - if new_n_unfixed == problem.n_unfixed && n_fixed == 0 - @debug "apply_branch: No progress made (n_unfixed same and n_fixed=0)" - return (TNProblem(problem.static, fill(DomainMask(0x00), length(propagated_doms)), 0, problem.ws), 0) - end - - # Create new problem with updated domains - new_problem = TNProblem(problem.static, propagated_doms, new_n_unfixed, problem.ws) - clear_region_cache!(problem) - return (new_problem, 1) # local_value = 1 (no scoring for now) -end - -function OptimalBranchingCore.reduce_problem(::Type{T}, problem::TNProblem, ::OptimalBranchingCore.NoReducer) where T - return (problem, one(T)) -end - -# function reduce_problem(::Type{T}, problem::TNProblem, reducer::UnitPropagationReducer) where T -# propagated = propagate(problem.static, problem.doms) - -# doms = problem.doms -# changed = false -# n_unfixed::Int = 0 - -# @inbounds for i in eachindex(doms) -# dm_new = propagated[i] -# bits = dm_new.bits - -# if bits == 0x00 -# clear_region_cache!(problem) -# return (TNProblem(problem.static, propagated, 0, problem.ws), zero(T)) -# end - -# changed |= bits != doms[i].bits -# n_unfixed += is_fixed(dm_new) ? 0 : 1 -# end - -# if !changed -# return (problem, one(T)) -# end - -# clear_region_cache!(problem) -# return (TNProblem(problem.static, propagated, n_unfixed, problem.ws), one(T)) -# end diff --git a/src/branching/branch.jl b/src/branching/branch.jl new file mode 100644 index 0000000..83df429 --- /dev/null +++ b/src/branching/branch.jl @@ -0,0 +1,349 @@ +function OptimalBranchingCore.branch_and_reduce( + problem::TNProblem, + config::OptimalBranchingCore.BranchingStrategy, + reducer::OptimalBranchingCore.AbstractReducer, + result_type::Type{TR}; + show_progress::Bool=false, + tag::Vector{Tuple{Int,Int}}=Tuple{Int,Int}[] +) where TR + try + stats = problem.ws.branch_stats + current_depth = length(tag) + + # Step 1: Check if problem is solved (all variables fixed) + if is_solved(problem) + record_solved_leaf!(stats, current_depth, problem.ws.current_path) + @debug "problem is solved at depth $current_depth" + cache_branch_solution!(problem) + return one(result_type) + end + + record_depth!(stats, current_depth) + @debug "======= Decision Level: $(current_depth + 1) =======" + # println("======= Decision Level: $(current_depth + 1) =======") + # Step 2: Try to reduce the problem + @assert reducer isa NoReducer + + # record branching time (only if detailed stats are enabled) + has_detailed = !isnothing(stats.detailed) + branching_start_time = has_detailed ? time_ns() : 0 + + # Step 3: Select variables for branching + variable = OptimalBranchingCore.select_variables(problem, config.measure, config.selector) + + # Record variable selection statistics and push to current path (only if tracking) + record_variable_selection!(stats, variable, problem.n_unfixed, current_depth) + needs_path_tracking(stats) && push!(problem.ws.current_path, variable) + @debug "Variable $variable selected at depth $current_depth" + + # Step 4: Compute branching table + tbl, variables = OptimalBranchingCore.branching_table(problem, config.table_solver, variable) + + # Check if table is empty (UNSAT - no valid configurations) + if isempty(tbl.table) + @debug "Empty branching table - UNSAT" + has_detailed && record_branching_time!(stats, (time_ns() - branching_start_time) / 1e9) + record_unsat_leaf!(stats, current_depth) + # Pop current variable before returning (restore path for parent) + needs_path_tracking(stats) && !isempty(problem.ws.current_path) && pop!(problem.ws.current_path) + return zero(result_type) + end + + # Step 5: Compute optimal branching rule + result = OptimalBranchingCore.optimal_branching_rule(tbl, variables, problem, config.measure, config.set_cover_solver) + + # record branching time (only if detailed stats are enabled) + has_detailed && record_branching_time!(stats, (time_ns() - branching_start_time) / 1e9) + + # @show tbl + # Step 6: Branch and recurse + clauses = OptimalBranchingCore.get_clauses(result) + @debug "A new branch-level search starts with $(length(clauses)) clauses: $(clauses)" + # @show clauses + # Record branching statistics + record_branch!(stats, length(clauses), current_depth) + # @show clauses + + # Use explicit loop instead of sum() to avoid closure allocation overhead + accum = zero(result_type) + + @inbounds for (i, branch) in enumerate(clauses) + show_progress && (OptimalBranchingCore.print_sequence(stdout, tag); println(stdout)) + @debug "branch=$branch, n_unfixed=$(problem.n_unfixed)" + + # Apply branch to get subproblem + subproblem, local_value, changed_vars = OptimalBranchingCore.apply_branch(problem, branch, variables) + + @debug "local_value=$local_value, n_unfixed=$(subproblem.n_unfixed)" + + # If branch led to contradiction (UNSAT), skip this branch + if local_value == 0 || subproblem.n_unfixed == 0 && any(dm -> dm.bits == 0x00, subproblem.doms) + @debug "Skipping branch: local_value=$local_value, n_unfixed=$(subproblem.n_unfixed), has_contradiction=$(any(dm -> dm.bits == 0x00, subproblem.doms))" + record_skipped_subproblem!(stats) + record_unsat_leaf!(stats, current_depth + 1) + continue # Skip to next branch instead of creating zero value + end + + # Recursively solve subproblem + # Use mutable tag buffer: push before recursion, pop after + # This avoids allocating new arrays on each recursive call + push!(tag, (i, length(clauses))) + # Save path length before recursion to restore after (only if tracking) + path_length_before = needs_path_tracking(stats) ? length(problem.ws.current_path) : 0 + sub_result = OptimalBranchingCore.branch_and_reduce(subproblem, config, reducer, result_type; tag=tag, show_progress=show_progress) + # Restore path to length before recursion (pop variables added during recursion) + if needs_path_tracking(stats) + while length(problem.ws.current_path) > path_length_before + pop!(problem.ws.current_path) + end + end + pop!(tag) + + # Combine results + accum = accum + sub_result * result_type(local_value) + + # Early exit: if we found a solution (accum > 0), return immediately + # This ensures we only find one solution path + if accum > zero(result_type) + @debug "Found solution, early exit at depth $current_depth" + needs_path_tracking(stats) && !isempty(problem.ws.current_path) && pop!(problem.ws.current_path) + return accum + end + end + + # Pop current variable before returning (restore path for parent) + needs_path_tracking(stats) && !isempty(problem.ws.current_path) && pop!(problem.ws.current_path) + return accum + finally + clear_branch_cache!(problem.ws, Base.objectid(problem.doms)) + end +end + +@inline function lookup_branch_cache(problem::TNProblem, + clause::OptimalBranchingCore.Clause{INT}, + variables::Vector{Int} +) where {INT} + ws = problem.ws + doms_id = Base.objectid(problem.doms) + inner = get(ws.branch_cache, doms_id, nothing) + inner === nothing && return nothing + key = (Base.objectid(variables), clause_key(clause)) + return get(inner, key, nothing) +end + +@inline function store_branch_cache!( + problem::TNProblem, + clause::OptimalBranchingCore.Clause{INT}, + variables::Vector{Int}, + doms::Vector{DomainMask}, + n_unfixed::Int, + local_value::Int +) where {INT} + ws = problem.ws + doms_id = Base.objectid(problem.doms) + inner = get!(ws.branch_cache, doms_id) do + Dict{Tuple{UInt, Any}, Any}() + end + key = (Base.objectid(variables), clause_key(clause)) + inner[key] = (doms=doms, n_unfixed=n_unfixed, local_value=local_value) + return nothing +end + +function OptimalBranchingCore.optimal_branching_rule(tbl::OptimalBranchingCore.BranchingTable, variables::Vector{T}, problem::TNProblem, measure::OptimalBranchingCore.AbstractMeasure, solver::OptimalBranchingCore.GreedyMerge) where T + candidates = OptimalBranchingCore.bit_clauses(tbl) + return OptimalBranchingCore.greedymerge(candidates, problem, variables, measure) +end + +function OptimalBranchingCore.optimal_branching_rule( + tbl::OptimalBranchingCore.BranchingTable{INT}, + variables::Vector{T}, + problem::TNProblem, + measure::OptimalBranchingCore.AbstractMeasure, + solver::OptimalBranchingCore.AbstractSetCoverSolver +) where {INT<:Integer, T} + candidates = OptimalBranchingCore.candidate_clauses(tbl) + valid_clauses = Vector{OptimalBranchingCore.Clause{INT}}() + reductions = Float64[] + + for clause in candidates + reduction = Float64(OptimalBranchingCore.size_reduction(problem, measure, clause, variables)) + if isfinite(reduction) && reduction > 0 + push!(valid_clauses, clause) + push!(reductions, reduction) + end + end + + if isempty(valid_clauses) + empty_clauses = Vector{OptimalBranchingCore.Clause{INT}}() + return OptimalBranchingCore.OptimalBranchingResult( + OptimalBranchingCore.DNF(empty_clauses), + Float64[], + Inf, + ) + end + + covered_mask = BitVector(undef, length(tbl.table)) + @inbounds for (idx, group) in pairs(tbl.table) + covered_mask[idx] = any(valid_clauses) do clause + any(config -> OptimalBranchingCore.covered_by(config, clause), group) + end + end + + if !all(covered_mask) + if !any(covered_mask) + empty_clauses = Vector{OptimalBranchingCore.Clause{INT}}() + return OptimalBranchingCore.OptimalBranchingResult( + OptimalBranchingCore.DNF(empty_clauses), + Float64[], + Inf, + ) + end + filtered_groups = tbl.table[findall(covered_mask)] + tbl = OptimalBranchingCore.BranchingTable{INT}(tbl.bit_length, filtered_groups) + end + + return OptimalBranchingCore.minimize_γ(tbl, valid_clauses, reductions, solver) +end + +function OptimalBranchingCore.optimal_branching_rule( + tbl::OptimalBranchingCore.BranchingTable{INT}, + variables::Vector{T}, + problem::TNProblem, + measure::OptimalBranchingCore.AbstractMeasure, + solver::OptimalBranchingCore.NaiveBranch +) where {INT<:Integer, T} + return invoke( + OptimalBranchingCore.optimal_branching_rule, + Tuple{ + OptimalBranchingCore.BranchingTable, + Vector, + OptimalBranchingCore.AbstractProblem, + OptimalBranchingCore.AbstractMeasure, + OptimalBranchingCore.NaiveBranch, + }, + tbl, + variables, + problem, + measure, + solver, + ) +end + +# Apply a clause to domain masks, fixing variables according to the clause's mask and values. Returns the number of variables that were fixed. +function apply_clause_to_doms!(new_doms::Vector{DomainMask}, clause::OptimalBranchingCore.Clause{INT}, variables::Vector{Int}, original_doms::Vector{DomainMask}, changed_flags::BitVector, changed_indices::Vector{Int}) where {INT<:Integer} + n_fixed = 0 + + @inbounds for i in 1:length(variables) + mask_bit = OptimalBranchingCore.readbit(clause.mask, i) + if mask_bit == 1 + # This variable is fixed by the clause + var_id = variables[i] + bit_val = OptimalBranchingCore.readbit(clause.val, i) + new_val = ifelse(bit_val == 1, DM_1, DM_0) + + old_bits = original_doms[var_id].bits + # Branchless: check if not fixed (bits == 0x03) + if old_bits == 0x03 + new_doms[var_id] = new_val + n_fixed += 1 + changed_flags[var_id] = true + push!(changed_indices, var_id) + end + end + end + + return n_fixed +end + +function OptimalBranchingCore.apply_branch( + problem::TNProblem, + clause::OptimalBranchingCore.Clause{INT}, + variables::Vector{Int} +) where {INT<:Integer} + if (cached = lookup_branch_cache(problem, clause, variables)) !== nothing + subproblem = TNProblem(problem.static, cached.doms, cached.n_unfixed, problem.ws) + return (subproblem, cached.local_value, Int[]) + end + + new_doms = copy(problem.doms) + + changed_flags = problem.ws.changed_vars_flags # BitVector: tracking variable changes + changed_indices = problem.ws.changed_vars_indices + empty!(changed_indices) + + # Apply clause: fix variables according to mask and values + n_fixed = apply_clause_to_doms!(new_doms, clause, variables, problem.doms, changed_flags, changed_indices) + + # @debug "apply_branch: Clause $(clause) Fixed $n_fixed variables" + + propagated_doms = propagate(problem.static, new_doms, changed_indices, problem.ws) + + @inbounds for i in eachindex(propagated_doms) + if propagated_doms[i].bits != problem.doms[i].bits && !changed_flags[i] + changed_flags[i] = true + push!(changed_indices, i) + end + end + + if has_contradiction(propagated_doms) + # UNSAT: contradiction detected during propagation + @debug "apply_branch: Clause $(clause) Contradiction detected during propagation" + doms_zero = fill(DomainMask(0x00), length(propagated_doms)) + store_branch_cache!(problem, clause, variables, doms_zero, problem.n_unfixed, 0) + return (TNProblem(problem.static, doms_zero, problem.n_unfixed, problem.ws), 0, Int[]) + end + + # Count unfixed variables + new_n_unfixed = count_unfixed(propagated_doms) + + @debug "apply_branch: Clause $(clause) n_unfixed: $(problem.n_unfixed) -> $new_n_unfixed" + + # Safety check: problem must have gotten smaller OR we fixed at least one variable + if new_n_unfixed == problem.n_unfixed && n_fixed == 0 + @debug "apply_branch: No progress made (n_unfixed same and n_fixed=0)" + doms_zero = fill(DomainMask(0x00), length(propagated_doms)) + store_branch_cache!(problem, clause, variables, doms_zero, 0, 0) + return (TNProblem(problem.static, doms_zero, 0, problem.ws), 0, Int[]) + end + + new_problem = TNProblem(problem.static, propagated_doms, new_n_unfixed, problem.ws) + + @inbounds for idx in changed_indices + changed_flags[idx] = false + end + + store_branch_cache!(problem, clause, variables, propagated_doms, new_n_unfixed, 1) + return (new_problem, 1, Int[]) # local_value = 1 (no scoring for now) +end + +function OptimalBranchingCore.reduce_problem(::Type{T}, problem::TNProblem, ::OptimalBranchingCore.NoReducer) where T + return (problem, one(T)) +end + +# function reduce_problem(::Type{T}, problem::TNProblem, reducer::UnitPropagationReducer) where T +# propagated = propagate(problem.static, problem.doms) + +# doms = problem.doms +# changed = false +# n_unfixed::Int = 0 + +# @inbounds for i in eachindex(doms) +# dm_new = propagated[i] +# bits = dm_new.bits + +# if bits == 0x00 +# clear_region_cache!(problem) +# return (TNProblem(problem.static, propagated, 0, problem.ws), zero(T)) +# end + +# changed |= bits != doms[i].bits +# n_unfixed += is_fixed(dm_new) ? 0 : 1 +# end + +# if !changed +# return (problem, one(T)) +# end + +# clear_region_cache!(problem) +# return (TNProblem(problem.static, propagated, n_unfixed, problem.ws), one(T)) +# end diff --git a/src/branching/branchtable.jl b/src/branching/branchtable.jl new file mode 100644 index 0000000..eb6b04a --- /dev/null +++ b/src/branching/branchtable.jl @@ -0,0 +1,246 @@ +struct TNContractionSolver <: AbstractTableSolver + k::Int + max_tensors::Int +end +TNContractionSolver() = TNContractionSolver(1, 2) + +@inline function bit_extract(config::UInt64, positions::Vector{Int})::UInt64 + out::UInt64 = 0 + @inbounds for (newpos, oldpos) in enumerate(positions) + out |= ((config >> (oldpos-1)) & 0x1) << (newpos-1) + end + return out +end + +@inline bit_match(config::UInt64, mask::UInt64, pattern::UInt64) = (config & mask) == pattern + +function _build_axismap(output_vars::Vector{Int}) + isempty(output_vars) && return Int[] + max_var = maximum(output_vars) + axismap = zeros(Int, max_var) + @inbounds for (idx, var_id) in enumerate(output_vars) + axismap[var_id] = idx + end + return axismap +end + +function extract_inner_configs(contracted::AbstractArray{T, N}, n_inner::Int) where {T, N} + @assert n_inner == N + configs = Vector{UInt64}() + one_tropical = one(Tropical{Float64}) + + @inbounds for lin in eachindex(contracted) + if contracted[lin] == one_tropical + linear_idx = LinearIndices(contracted)[lin] - 1 + config_bits = UInt64(linear_idx) + push!(configs, config_bits) + end + end + return configs +end + +function combine_configs(boundary_bits::UInt64, inner_configs::Vector{UInt64}, n_boundary::Int) + n_configs = length(inner_configs) + full_configs = Vector{UInt64}(undef, n_configs) + + @inbounds for i in 1:n_configs + full_configs[i] = boundary_bits | (inner_configs[i] << n_boundary) + end + + return full_configs +end + +function slice_region_contraction( + tensor::AbstractArray{Tropical{Float64}}, + assignments::Vector{Tuple{Int,Bool}}, + axismap::Vector{Int}, +) + isempty(assignments) && return tensor + nd = ndims(tensor) + nd == 0 && return tensor + + indices = Any[Colon() for _ in 1:nd] + + @inbounds for (var_id, value) in assignments + axis = axismap[var_id] + @assert axis > 0 "Boundary variable $var_id not found in contraction axes" + indices[axis] = value ? 2 : 1 + end + + return @view tensor[indices...] +end + +function handle_no_boundary_case_unfixed( + region::Region, + contracted::AbstractArray{Tropical{Float64}}, + inner_output_vars::Vector{Int}, +) + n_inner = length(region.inner_vars) + free_inner_configs = extract_inner_configs(contracted, length(inner_output_vars)) + + if !isempty(free_inner_configs) + return BranchingTable(n_inner, [free_inner_configs]) + else + return BranchingTable(0, [UInt64[]]) + end +end + +function create_region(problem::TNProblem, variable::Int, solver::TNContractionSolver) + # Compute k-neighboring region using all-unfixed domains + # This ensures the region is consistent across different branches + all_unfixed_doms = fill(DM_BOTH, length(problem.doms)) + k_neighboring(problem.static, all_unfixed_doms, variable; max_tensors = solver.max_tensors, k = solver.k) +end + +# Filter cached BranchingTable: filter rows by fixed values, extract columns for unfixed vars +function filter_branching_table(region::Region, table::BranchingTable, problem::TNProblem) + stats = problem.ws.branch_stats + has_detailed = !isnothing(stats) && !isnothing(stats.detailed) + filtering_start_time = has_detailed ? time_ns() : 0 + + var_ids = vcat(region.boundary_vars, region.inner_vars) + n_vars = length(var_ids) + + fixed_positions = Tuple{Int, Bool}[] + unfixed_positions = Int[] + + @inbounds for (i, var_id) in enumerate(var_ids) + if is_fixed(problem.doms[var_id]) + required_value = has1(problem.doms[var_id]) + push!(fixed_positions, (i, required_value)) + else + push!(unfixed_positions, i) + end + end + + if isempty(fixed_positions) + if has_detailed + filtering_time = (time_ns() - filtering_start_time) / 1e9 + record_filtering_time!(stats, filtering_time) + end + return table, var_ids + end + + n_unfixed = length(unfixed_positions) + + mask::UInt64 = 0x0 + pattern::UInt64 = 0x0 + @inbounds for (pos, val) in fixed_positions + pos_mask = UInt64(1) << (pos - 1) + mask |= pos_mask + if val + pattern |= pos_mask + end + end + + filtered_table = similar(table.table, 0) + sizehint!(filtered_table, length(table.table)) + + @inbounds for config_group in table.table + filtered_group = similar(config_group, 0) + sizehint!(filtered_group, length(config_group)) + + for config_bits in config_group + if bit_match(UInt64(config_bits), mask, pattern) + new_config = bit_extract(UInt64(config_bits), unfixed_positions) + push!(filtered_group, typeof(config_bits)(new_config)) + end + end + + if !isempty(filtered_group) + push!(filtered_table, filtered_group) + end + end + + if has_detailed + filtering_time = (time_ns() - filtering_start_time) / 1e9 + record_filtering_time!(stats, filtering_time) + end + + isempty(filtered_table) && return BranchingTable(0, eltype(table.table)[]), Int[] + + unfixed_var_ids = [var_ids[i] for i in unfixed_positions] + + return BranchingTable(n_unfixed, filtered_table), unfixed_var_ids +end + +function OptimalBranchingCore.branching_table(problem::TNProblem, solver::TNContractionSolver, variable::Int) + stats = problem.ws.branch_stats + + cached_region, cached_table = get_cached_region(variable) + if !isnothing(cached_region) && !isnothing(cached_table) + record_cache_hit!(stats) + filtered_table, unfixed_vars = filter_branching_table(cached_region, cached_table, problem) + return filtered_table, unfixed_vars + end + + record_cache_miss!(stats) + region = create_region(problem, variable, solver) + n_boundary = length(region.boundary_vars) + n_inner = length(region.inner_vars) + n_total = n_boundary + n_inner + + # Contract with all-unfixed doms for consistent caching + all_unfixed_doms = fill(DM_BOTH, length(problem.doms)) + + stats = problem.ws.branch_stats + contraction_start_time = time_ns() + contracted_tensor, output_vars = contract_region(problem.static, region, all_unfixed_doms) + contraction_time = (time_ns() - contraction_start_time) / 1e9 + record_contraction_time!(stats, contraction_time) + axismap = _build_axismap(output_vars) + + inner_output_vars = Int[] + @inbounds for var_id in output_vars + if var_id in region.inner_vars + push!(inner_output_vars, var_id) + end + end + + if n_boundary == 0 + table = handle_no_boundary_case_unfixed(region, contracted_tensor, inner_output_vars) + variables = vcat(region.boundary_vars, region.inner_vars) + cache_region!(region, table) + filtered_table, unfixed_vars = filter_branching_table(region, table, problem) + return filtered_table, unfixed_vars + end + + free_boundary_vars = region.boundary_vars + n_free_boundary = length(free_boundary_vars) + + valid_config_groups = Vector{Vector{UInt64}}() + assignments = Tuple{Int,Bool}[] + resize!(assignments, n_free_boundary) + + for free_config in 0:(2^n_free_boundary - 1) + for (j, var_id) in enumerate(free_boundary_vars) + bit = (free_config >> (j-1)) & 0x1 + assignments[j] = (var_id, bit == 1) + end + + contracted_slice = slice_region_contraction(contracted_tensor, assignments, axismap) + free_inner_configs = extract_inner_configs(contracted_slice, length(inner_output_vars)) + + isempty(free_inner_configs) && continue + + boundary_bits = UInt64(free_config) + full_configs = combine_configs(boundary_bits, free_inner_configs, n_boundary) + + push!(valid_config_groups, full_configs) + end + + if isempty(valid_config_groups) + table = BranchingTable(0, [UInt64[]]) + cache_region!(region, table) + return table, Int[] + end + + table = BranchingTable(n_total, valid_config_groups) + cache_region!(region, table) + + filtered_table, unfixed_vars = filter_branching_table(region, table, problem) + return filtered_table, unfixed_vars +end + +# Constructor for MinGammaSelector (defined here after TNContractionSolver is available) +MinGammaSelector() = MinGammaSelector(TNContractionSolver(2,5), OptimalBranchingCore.GreedyMerge()) diff --git a/src/branching/greedymerge.jl b/src/branching/greedymerge.jl new file mode 100644 index 0000000..793510d --- /dev/null +++ b/src/branching/greedymerge.jl @@ -0,0 +1,222 @@ +# Calculate size reduction after applying a clause to the problem +function OptimalBranchingCore.size_reduction(p::TNProblem, m::OptimalBranchingCore.AbstractMeasure, cl::OptimalBranchingCore.Clause{INT}, variables::Vector{Int}) where {INT<:Integer} + subproblem, _, _ = OptimalBranchingCore.apply_branch(p, cl, variables) + if any(dm -> dm.bits == 0x00, subproblem.doms) + return -Inf + end + reduction = OptimalBranchingCore.measure(p, m) - OptimalBranchingCore.measure(subproblem, m) + return reduction +end + +# Generate a key for a clause based on its fields +@inline clause_key(cl::Clause{INT}) where {INT} = (getfield(cl, 1), getfield(cl, 2)) + +# Get size reduction from cache or compute and cache it +function cached_size_reduction!( + cache::Dict{Tuple{INT,INT},Float64}, + problem::OptimalBranchingCore.AbstractProblem, + m::OptimalBranchingCore.AbstractMeasure, + variables::Vector{Int}, + clause::Clause{INT} +) where {INT} + key = clause_key(clause) + return get!(cache, key) do + Float64(OptimalBranchingCore.size_reduction(problem, m, clause, variables)) + end +end + +# Select the best representative clause from each row based on size reduction +function select_representatives!( + cls::Vector{Vector{Clause{INT}}}, + cache::Dict{Tuple{INT,INT},Float64}, + problem::OptimalBranchingCore.AbstractProblem, + variables::Vector{Int}, + m::OptimalBranchingCore.AbstractMeasure +) where {INT} + size_reductions = Vector{Float64}(undef, length(cls)) + @inbounds for (idx, row) in pairs(cls) + best_val = -Inf + best_pos = 0 + @inbounds for k in eachindex(row) + reduction = cached_size_reduction!(cache, problem, m, variables, row[k]) + if reduction > best_val && reduction > 0 && isfinite(reduction) + best_val = reduction + best_pos = k + end + end + size_reductions[idx] = best_val + if best_pos > 1 + row[1], row[best_pos] = row[best_pos], row[1] + end + end + return size_reductions +end + +# Remove rows with invalid (non-finite or non-positive) size reductions +function drop_invalid_rows!(cls::Vector{Vector{Clause{INT}}}, size_reductions::Vector{Float64}) where {INT} + keep = isfinite.(size_reductions) .& (size_reductions .> 0) + all(keep) && return cls, size_reductions + return cls[keep], size_reductions[keep] +end + +# Remove duplicate singleton clauses from the clause set +function deduplicate_singletons!(cls::Vector{Vector{Clause{INT}}}, size_reductions::Vector{Float64}) where {INT} + seen = Dict{Tuple{INT,INT},Int}() + mask = trues(length(cls)) + @inbounds for i in 1:length(cls) + if length(cls[i]) == 1 + key = clause_key(cls[i][1]) + if haskey(seen, key) + mask[i] = false + else + seen[key] = i + end + end + end + if all(mask) + return cls, size_reductions + end + return cls[mask], size_reductions[mask] +end + +# Find the best clause to merge from two clause vectors +function best_merge_clause( + cache::Dict{Tuple{INT,INT},Float64}, + problem::OptimalBranchingCore.AbstractProblem, + m::OptimalBranchingCore.AbstractMeasure, + variables::Vector{Int}, + left::Vector{Clause{INT}}, + right::Vector{Clause{INT}} +) where {INT} + best_clause = nothing + best_reduction = -Inf + nvars = length(variables) + @inbounds for cl_left in left + @inbounds for cl_right in right + merged = OptimalBranchingCore.gather2(nvars, cl_left, cl_right) + iszero(merged.mask) && continue + reduction = cached_size_reduction!(cache, problem, m, variables, merged) + if reduction > best_reduction && reduction > 0 && isfinite(reduction) + best_clause = merged + best_reduction = reduction + end + end + end + return best_clause, best_reduction +end + +# Enqueue a merge pair if it's beneficial (finite positive reduction with negative energy change) +@inline function maybe_enqueue!(queue::PriorityQueue{NTuple{2,Int},Float64}, idx_a::Int, idx_b::Int, reduction::Float64, weights::Vector{Float64}, γ::Float64) + if !isfinite(reduction) || reduction <= 0 + return + end + dE = γ^(-reduction) - weights[idx_a] - weights[idx_b] + dE <= -1e-12 && enqueue!(queue, (idx_a, idx_b), dE) +end + +# Enqueue all beneficial merge pairs from the clause set +function enqueue_beneficial_merges!( + queue::PriorityQueue{NTuple{2,Int},Float64}, + cls::Vector{Vector{Clause{INT}}}, + γ::Float64, + weights::Vector{Float64}, + cache::Dict{Tuple{INT,INT},Float64}, + problem::OptimalBranchingCore.AbstractProblem, + variables::Vector{Int}, + m::OptimalBranchingCore.AbstractMeasure +) where {INT} + nc = length(cls) + @inbounds for i in 1:nc-1 + row_i = cls[i] + for j in i+1:nc + clause, reduction = best_merge_clause(cache, problem, m, variables, row_i, cls[j]) + clause === nothing && continue + maybe_enqueue!(queue, i, j, reduction, weights, γ) + end + end +end + +# Remove queue entries related to a specific row +function purge_queue_entries!(queue::PriorityQueue{NTuple{2,Int},Float64}, mask::BitVector, rowid::Int, nc::Int) + mask[rowid] = false + @inbounds for l in 1:nc + if mask[l] + a, b = minmax(rowid, l) + haskey(queue, (a, b)) && delete!(queue, (a, b)) + end + end +end + +# Process the merge queue and perform beneficial merges +function process_merge_queue!( + queue::PriorityQueue{NTuple{2,Int},Float64}, + cls::Vector{Vector{Clause{INT}}}, size_reductions::Vector{Float64}, + weights::Vector{Float64}, mask::BitVector, γ::Float64, + cache::Dict{Tuple{INT,INT},Float64}, + problem::OptimalBranchingCore.AbstractProblem, + variables::Vector{Int}, + m::OptimalBranchingCore.AbstractMeasure +) where {INT} + nc = length(cls) + any_merge = false + while !isempty(queue) + (i, j) = dequeue!(queue) + if !mask[i] || !mask[j] + continue + end + clause, new_reduction = best_merge_clause(cache, problem, m, variables, cls[i], cls[j]) + clause === nothing && continue + + purge_queue_entries!(queue, mask, i, nc) + purge_queue_entries!(queue, mask, j, nc) + + mask[i] = true + cls[i] = [clause] + size_reductions[i] = new_reduction + weights[i] = γ^(-new_reduction) + any_merge = true + + @inbounds for l in 1:nc + if i !== l && mask[l] + a, b = minmax(i, l) + _, reduction = best_merge_clause(cache, problem, m, variables, cls[a], cls[b]) + maybe_enqueue!(queue, a, b, reduction, weights, γ) + end + end + end + return any_merge +end + +# Greedy merge algorithm to find optimal branching clauses +function OptimalBranchingCore.greedymerge(cls::Vector{Vector{Clause{INT}}}, problem::TNProblem, variables::Vector{T}, m::OptimalBranchingCore.AbstractMeasure) where {INT<:Integer, T} + cache = Dict{Tuple{INT,INT},Float64}() + working = copy(cls) + size_reductions = select_representatives!(working, cache, problem, variables, m) + working, size_reductions = drop_invalid_rows!(working, size_reductions) + if isempty(working) + @debug "greedymerge: no valid clauses remain after filtering; returning empty result" + return OptimalBranchingResult(DNF(Clause{INT}[]), Float64[], Inf) + end + + while true + working, size_reductions = deduplicate_singletons!(working, size_reductions) + nc = length(working) + mask = trues(nc) + γ = OptimalBranchingCore.complexity_bv(size_reductions) + weights = map(s -> γ^(-s), size_reductions) + queue = PriorityQueue{NTuple{2, Int}, Float64}() + + enqueue_beneficial_merges!(queue, working, γ, weights, cache, problem, variables, m) + if isempty(queue) + return OptimalBranchingResult(DNF(first.(working)), size_reductions, γ) + end + + any_merge = process_merge_queue!(queue, working, size_reductions, weights, mask, γ, cache, problem, variables, m) + if !any_merge + return OptimalBranchingResult(DNF(first.(working)), size_reductions, γ) + end + + working = working[mask] + size_reductions = size_reductions[mask] + end +end diff --git a/src/measure.jl b/src/branching/measure.jl similarity index 100% rename from src/measure.jl rename to src/branching/measure.jl diff --git a/src/branching/selector.jl b/src/branching/selector.jl new file mode 100644 index 0000000..84378b2 --- /dev/null +++ b/src/branching/selector.jl @@ -0,0 +1,87 @@ +struct MostOccurrenceSelector <: AbstractSelector end + +function OptimalBranchingCore.select_variables( + problem::TNProblem, + measure::AbstractMeasure, + selector::MostOccurrenceSelector +) + unfixed_vars = get_unfixed_vars(problem.doms) + isempty(unfixed_vars) && return Int[] + least_show_var_idx = argmax(length(problem.static.v2t[u]) for u in unfixed_vars) + + return unfixed_vars[least_show_var_idx] +end + +""" + MinGammaSelector <: AbstractSelector + +Selector that chooses the variable with minimum branching gamma. +For each unfixed variable, it computes the branching gamma on the corresponding +region and selects the variable with the smallest gamma value. + +This selector does not consider performance optimization and will iterate through +all unfixed variables to find the best one. +""" +struct MinGammaSelector <: AbstractSelector + table_solver::AbstractTableSolver + set_cover_solver::OptimalBranchingCore.AbstractSetCoverSolver +end + +# Constructor will be defined after TNContractionSolver is available + +function OptimalBranchingCore.select_variables( + problem::TNProblem, + measure::AbstractMeasure, + selector::MinGammaSelector +) + unfixed_vars = get_unfixed_vars(problem.doms) + isempty(unfixed_vars) && return Int[] + + best_var = unfixed_vars[1] + best_gamma = Inf + + @debug "Selecting variables with MinGammaSelector, n_unfixed=$(length(unfixed_vars))" + # println("Selecting variables with MinGammaSelector, n_unfixed=$(length(unfixed_vars))") + # Iterate through all unfixed variables + for var in unfixed_vars + # Compute branching table for this variable + tbl, variables = OptimalBranchingCore.branching_table( + problem, + selector.table_solver, + var + ) + # @show tbl + # Skip if table is empty (UNSAT) + if isempty(tbl.table) + continue + end + + # Compute optimal branching rule for this variable + result = OptimalBranchingCore.optimal_branching_rule( + tbl, + variables, + problem, + measure, + selector.set_cover_solver + ) + + @debug "Optimal branching rule for variable $var: $(OptimalBranchingCore.get_clauses(result))" + # println("Optimal branching rule for variable $var: $(OptimalBranchingCore.get_clauses(result))") + + # Get the gamma value + gamma = result.γ + @debug "Gamma for variable $var: $gamma" + # println("Gamma for variable $var: $gamma") + + # Update best variable if this gamma is smaller + if gamma < best_gamma + @debug "Updating best variable to $var with gamma $gamma" + # println("Updating best variable to $var with gamma $gamma") + best_gamma = gamma + best_var = var + end + end + + return best_var +end + diff --git a/src/branchtable.jl b/src/branchtable.jl deleted file mode 100644 index 2f8fa0a..0000000 --- a/src/branchtable.jl +++ /dev/null @@ -1,223 +0,0 @@ -using OptimalBranchingCore: BranchingTable - -struct TNContractionSolver <: AbstractTableSolver end - -function _build_axismap(output_vars::Vector{Int}) - isempty(output_vars) && return Int[] - max_var = maximum(output_vars) - axismap = zeros(Int, max_var) - @inbounds for (idx, var_id) in enumerate(output_vars) - axismap[var_id] = idx - end - return axismap -end - -function separate_fixed_free_boundary(region::Region, doms::Vector{DomainMask}) - fixed_vars = Int[]; fixed_vals = Int[]; free_vars = Int[]; free_indices = Int[] - - for (i, var_id) in enumerate(region.boundary_vars) - if is_fixed(doms[var_id]) - push!(fixed_vars, var_id) - push!(fixed_vals, has1(doms[var_id]) ? 1 : 0) - else - push!(free_vars, var_id) - push!(free_indices, i) - end - end - return fixed_vars, fixed_vals, free_vars, free_indices -end - -function construct_boundary_config(region::Region, doms::Vector{DomainMask}, free_pos_of_boundary::Vector{Int}, free_config::Int) - n_boundary = length(region.boundary_vars) - boundary_config = Vector{Bool}(undef, n_boundary) - - @inbounds for i in 1:n_boundary - var_id = region.boundary_vars[i] - if is_fixed(doms[var_id]) - boundary_config[i] = has1(doms[var_id]) - else - j = free_pos_of_boundary[i] - bit = (free_config >> (j-1)) & 0x1 - boundary_config[i] = (bit == 1) - end - end - - return boundary_config -end - -function construct_inner_config(region::Region, doms::Vector{DomainMask}, inner_posmap::Vector{Int}, free_inner_configs::Vector{Vector{Bool}}) - n_inner = length(region.inner_vars) - n_config = length(free_inner_configs) - inner_configs = Vector{Vector{Bool}}(undef, n_config) - - @inbounds for config_idx in 1:n_config - inner_config = Vector{Bool}(undef, n_inner) - for i in 1:n_inner - var_id = region.inner_vars[i] - if is_fixed(doms[var_id]) - inner_config[i] = has1(doms[var_id]) - else - j = inner_posmap[var_id] - inner_config[i] = free_inner_configs[config_idx][j] - end - end - inner_configs[config_idx] = inner_config - end - - return inner_configs -end - -function extract_inner_configs(contracted::AbstractArray{T, N}, n_inner::Int) where {T, N} - @assert n_inner == N - configs = Vector{Vector{Bool}}() - one_tropical = one(Tropical{Float64}) - - @inbounds for lin in eachindex(contracted) - if contracted[lin] == one_tropical - config = Vector{Bool}(undef, n_inner) - linear_idx = LinearIndices(contracted)[lin] - 1 - for i in 1:n_inner - config[i] = ((linear_idx >> (i-1)) & 0x1) == 1 - end - push!(configs, config) - end - end - return configs -end - -function combine_configs(boundary_config::Vector{Bool}, inner_configs::Vector{Vector{Bool}}) - n_configs = length(inner_configs) - n_boundary = length(boundary_config) - n_inner = isempty(inner_configs) ? 0 : length(inner_configs[1]) - n_total = n_boundary + n_inner - - full_configs = Vector{Vector{Bool}}(undef, n_configs) - - @inbounds for i in 1:n_configs - full_config = Vector{Bool}(undef, n_total) - copyto!(full_config, 1, boundary_config, 1, n_boundary) - copyto!(full_config, n_boundary + 1, inner_configs[i], 1, n_inner) - full_configs[i] = full_config - end - - return full_configs -end - -function get_region_contraction(problem::TNProblem, region::Region) - cached = get_cached_region_contraction(problem; region_id=region.id) - if isnothing(cached) - @debug "contract_region_cache_miss" focus=region.id - contracted, output_vars = contract_region(problem.static, region, problem.doms) - cache_region_contraction!(problem, contracted, output_vars; region_id=region.id) - return contracted, output_vars - else - @debug "contract_region_cache_hit" focus=region.id - return cached - end -end - -function slice_region_contraction( - tensor::AbstractArray{Tropical{Float64}}, - assignments::Vector{Tuple{Int,Bool}}, - axismap::Vector{Int}, -) - isempty(assignments) && return tensor - nd = ndims(tensor) - nd == 0 && return tensor - - indices = Any[Colon() for _ in 1:nd] - - @inbounds for (var_id, value) in assignments - axis = axismap[var_id] - @assert axis > 0 "Boundary variable $var_id not found in contraction axes" - indices[axis] = value ? 2 : 1 - end - - return @view tensor[indices...] -end - -function handle_no_boundary_case( - problem::TNProblem, - region::Region, - contracted::AbstractArray{Tropical{Float64}}, - inner_output_vars::Vector{Int}, -) - n_inner = length(region.inner_vars) - - free_inner_configs = extract_inner_configs(contracted, length(inner_output_vars)) - - if !isempty(free_inner_configs) - inner_posmap = _build_axismap(inner_output_vars) - inner_configs = construct_inner_config( - region, problem.doms, inner_posmap, free_inner_configs - ) - return BranchingTable(n_inner, [inner_configs]) - else - return BranchingTable(0, [Int[]]) - end -end - -function OptimalBranchingCore.branching_table( - problem::TNProblem, - solver::TNContractionSolver, - variables::Vector{T} -) where T - region = get_cached_region(problem) - isnothing(region) && error("No cached region found. Make sure `select_variables` is called first.") - - n_boundary = length(region.boundary_vars) - n_inner = length(region.inner_vars) - n_total = n_boundary + n_inner - @assert n_total == length(variables) - - contracted_tensor, output_vars = get_region_contraction(problem, region) - axismap = _build_axismap(output_vars) - - inner_output_vars = Int[] - @inbounds for var_id in output_vars - if var_id in region.inner_vars - push!(inner_output_vars, var_id) - end - end - - n_boundary == 0 && return handle_no_boundary_case(problem, region, contracted_tensor, inner_output_vars) - - _, _, free_boundary_vars, free_boundary_indices = separate_fixed_free_boundary(region, problem.doms) - - n_free_boundary = length(free_boundary_vars) - - free_pos_of_boundary = zeros(Int, n_boundary) - @inbounds for (j, boundary_idx) in enumerate(free_boundary_indices) - free_pos_of_boundary[boundary_idx] = j - end - - inner_posmap = _build_axismap(inner_output_vars) - - valid_config_groups = Vector{Vector{Vector{Bool}}}() - assignments = Tuple{Int,Bool}[] - resize!(assignments, n_free_boundary) - - for free_config in 0:(2^n_free_boundary - 1) - for (j, var_id) in enumerate(free_boundary_vars) - bit = (free_config >> (j-1)) & 0x1 - assignments[j] = (var_id, bit == 1) - end - - contracted_slice = slice_region_contraction(contracted_tensor, assignments, axismap) - - free_inner_configs = extract_inner_configs(contracted_slice, length(inner_output_vars)) - - isempty(free_inner_configs) && continue - - boundary_config = construct_boundary_config(region, problem.doms, free_pos_of_boundary, free_config) - - inner_configs = construct_inner_config(region, problem.doms, inner_posmap, free_inner_configs) - - full_configs = combine_configs(boundary_config, inner_configs) - - push!(valid_config_groups, full_configs) - end - - isempty(valid_config_groups) && return BranchingTable(0, [Int[]]) - return BranchingTable(n_total, valid_config_groups) -end diff --git a/src/core/domain.jl b/src/core/domain.jl new file mode 100644 index 0000000..9903926 --- /dev/null +++ b/src/core/domain.jl @@ -0,0 +1,21 @@ +struct DomainMask + bits::UInt8 +end + +const DM_BOTH = DomainMask(0x03) +const DM_0 = DomainMask(0x01) +const DM_1 = DomainMask(0x02) + +@inline is_fixed(dm::DomainMask) = (dm.bits == 0x01) || (dm.bits == 0x02) +@inline has0(dm::DomainMask)::Bool = (dm.bits & 0x01) != 0 +@inline has1(dm::DomainMask)::Bool = (dm.bits & 0x02) != 0 +function get_var_value(dms::Vector{DomainMask}, var_id::Int) + dm = dms[var_id] + dm.bits == 0x01 && return 0 + dm.bits == 0x02 && return 1 + return -1 # not fixed +end + +init_doms(static::TNStatic) = fill(DM_BOTH, length(static.vars)) + +@inline has_contradiction(doms::Vector{DomainMask}) = any(dm -> dm.bits == 0x00, doms) \ No newline at end of file diff --git a/src/core/problem.jl b/src/core/problem.jl new file mode 100644 index 0000000..e5429a2 --- /dev/null +++ b/src/core/problem.jl @@ -0,0 +1,54 @@ +struct TNProblem <: AbstractProblem + static::TNStatic + doms::Vector{DomainMask} + n_unfixed::Int + ws::DynamicWorkspace +end + +function TNProblem(static::TNStatic; verbose::Bool = false)::TNProblem + doms = init_doms(static) + ws = DynamicWorkspace(length(static.vars), verbose) + doms = propagate(static, doms, collect(1:length(static.vars)), ws) + n_unfixed = count_unfixed(doms) + return TNProblem(static, doms, n_unfixed, ws) +end + +function Base.show(io::IO, problem::TNProblem) + print(io, "TNProblem(unfixed=$(problem.n_unfixed)/$(length(problem.static.vars)))") +end + +get_var_value(problem::TNProblem, var_id::Int) = get_var_value(problem.doms, var_id) +get_var_value(problem::TNProblem, var_ids::Vector{Int}) = Bool[get_var_value(problem.doms, var_id) for var_id in var_ids] + +is_solved(problem::TNProblem) = problem.n_unfixed == 0 + +function cache_branch_solution!(problem::TNProblem) + ws = problem.ws + copyto!(ws.cached_doms, problem.doms) + ws.has_cached_solution = true + return nothing +end + +function reset_last_branch_problem!(problem::TNProblem) + problem.ws.has_cached_solution = false + return nothing +end + +@inline has_last_branch_problem(problem::TNProblem) = problem.ws.has_cached_solution + +function last_branch_problem(problem::TNProblem) + has_last_branch_problem(problem) || return nothing + doms = copy(problem.ws.cached_doms) + fixed = count(x -> is_fixed(x), doms) + @assert fixed == length(doms) + return TNProblem(problem.static, doms, 0, problem.ws) +end + +function get_branching_stats(problem::TNProblem) + return copy(problem.ws.branch_stats) +end + +function reset_branching_stats!(problem::TNProblem) + reset!(problem.ws.branch_stats) + return nothing +end diff --git a/src/core/region.jl b/src/core/region.jl new file mode 100644 index 0000000..bc3e78f --- /dev/null +++ b/src/core/region.jl @@ -0,0 +1,41 @@ +struct Region + id::Int + tensors::Vector{Int} + inner_vars::Vector{Int} + boundary_vars::Vector{Int} +end + +function Base.show(io::IO, region::Region) + print(io, "Region(focus=$(region.id), tensors=$(length(region.tensors)), inner_vars=$(region.inner_vars), boundary_vars=$(region.boundary_vars))") +end + +struct RegionCacheEntry + region::Region + table::BranchingTable +end + +mutable struct RegionCacheState + entries::Dict{Int, RegionCacheEntry} +end + +RegionCacheState() = RegionCacheState(Dict{Int, RegionCacheEntry}()) + +# global region cache +const REGION_CACHE = RegionCacheState() + +function cache_region!(region::Region, table::BranchingTable) + REGION_CACHE.entries[region.id] = RegionCacheEntry(region, table) + return nothing +end + +function get_cached_region(region_id::Int) + entry = get(REGION_CACHE.entries, region_id, nothing) + (isnothing(entry) || isnothing(entry.table)) && return nothing, nothing + return entry.region, entry.table +end + +function clear_all_region_caches!() + empty!(REGION_CACHE.entries) + return nothing +end + diff --git a/src/core/static.jl b/src/core/static.jl new file mode 100644 index 0000000..db7c839 --- /dev/null +++ b/src/core/static.jl @@ -0,0 +1,100 @@ +struct TNStatic + vars::Vector{Variable} + tensors::Vector{BoolTensor} + v2t::Vector{Vector{Int}} + t2v::Vector{Vector{EdgeRef}} + axis_of_t::Dict{Tuple{Int,Int},Int} + # Precomputed tensor masks: maps tensor.tensor to TensorMasks + # Reuse masks for identical tensor shapes to avoid recomputation + precomputed_masks::Dict{Vector{Tropical{Float64}}, TensorMasks} + # Map each tensor instance to its precomputed masks (indexing optimization) + tensor_to_masks::Vector{TensorMasks} +end + +function Base.show(io::IO, tn::TNStatic) + print(io, "TNStatic(vars=$(length(tn.vars)), tensors=$(length(tn.tensors)))") +end + +function Base.show(io::IO, ::MIME"text/plain", tn::TNStatic) + println(io, "TNStatic:") + println(io, " variables: $(length(tn.vars))") + println(io, " tensors: $(length(tn.tensors))") + println(io, " variable-tensor incidence: $(length(tn.v2t))") + println(io, " tensor-variable incidence: $(length(tn.t2v))") + println(io, " axis mappings: $(length(tn.axis_of_t))") +end + +function setup_problem(var_num::Int, + tensors_to_vars::Vector{Vector{Int}}, + tensor_data::Vector{Vector{Tropical{Float64}}}) + F = length(tensors_to_vars) + tensors = Vector{BoolTensor}(undef, F) + vars_to_tensors = [Int[] for _ in 1:var_num] + tensors_to_edges = [EdgeRef[] for _ in 1:F] + for i in 1:F + var_axes = tensors_to_vars[i] + @assert length(tensor_data[i]) == 1 << length(var_axes) + tensors[i] = BoolTensor(i, var_axes, tensor_data[i]) + for (j, v) in enumerate(var_axes) + push!(vars_to_tensors[v], i) + push!(tensors_to_edges[i], EdgeRef(v, j)) + end + end + + axis_of_t = Dict{Tuple{Int,Int},Int}() + for (fid, tensor) in enumerate(tensors) + for (axis, vid) in enumerate(tensor.var_axes) + axis_of_t[(fid, vid)] = axis + end + end + + vars = Vector{Variable}(undef, var_num) + for i in 1:var_num + vars[i] = Variable(i, 2, length(vars_to_tensors[i])) + end + + # Precompute masks for every unique tensor + # Strategy: use hash values to quickly find identical tensors and avoid recomputation + tensor_to_masks = Vector{TensorMasks}(undef, F) + + # Dictionary keyed by hash for fast lookup + hash_to_masks = Dict{UInt, TensorMasks}() + hash_to_tensor = Dict{UInt, Vector{Tropical{Float64}}}() + + for (i, tensor) in enumerate(tensors) + # Compute the hash of the tensor contents + tensor_hash = hash(tensor.tensor) + + if haskey(hash_to_masks, tensor_hash) + # Hash collision check: ensure the contents truly match + if hash_to_tensor[tensor_hash] == tensor.tensor + # Reuse the previously computed masks + tensor_to_masks[i] = hash_to_masks[tensor_hash] + else + error("Hash collision: $tensor_hash") + end + else + # First time we see this hash, compute the masks + masks = build_tensor_masks(tensor) + hash_to_masks[tensor_hash] = masks + hash_to_tensor[tensor_hash] = tensor.tensor + tensor_to_masks[i] = masks + end + end + + # Build the precomputed_masks dictionary for statistics (from unique tensors) + precomputed_masks = Dict{Vector{Tropical{Float64}}, TensorMasks}( + hash_to_tensor[h] => hash_to_masks[h] for h in keys(hash_to_masks) + ) + + return TNStatic(vars, tensors, vars_to_tensors, tensors_to_edges, axis_of_t, + precomputed_masks, tensor_to_masks) +end + +function setup_from_tensor_network(tn)::TNStatic + t2v = getixsv(tn.code) + tensors = GenericTensorNetworks.generate_tensors(Tropical(1.0), tn) + vec_tensors = [vec(t) for t in tensors] + new_tensors = [replace(t, Tropical(1.0) => zero(Tropical{Float64})) for t in vec_tensors] + return setup_problem(length(tn.problem.symbols), t2v, new_tensors) +end diff --git a/src/core/stats.jl b/src/core/stats.jl new file mode 100644 index 0000000..a2b0761 --- /dev/null +++ b/src/core/stats.jl @@ -0,0 +1,369 @@ +# Detailed statistics (optional, only allocated when needed) +mutable struct DetailedStats + # Depth distribution (node count per depth) + depth_distribution::Vector{Int} + + # Branching factor history (actual child count per branch) + branching_factors::Vector{Int} + + # Propagation efficiency + propagation_calls::Int + domain_reductions::Int + early_unsat_detections::Int + + # Time breakdown (seconds) + time_propagation::Float64 + time_branching::Float64 + time_contraction::Float64 + time_filtering::Float64 + + # Cache efficiency + cache_hits::Int + cache_misses::Int + + # Branch decision quality + variable_selection_counts::Dict{Int, Int} + remaining_vars_at_branch::Vector{Int} + depth_at_selection::Vector{Int} # Depth when each variable was selected + variable_selection_sequence::Vector{Int} # Variable selection order (by pick sequence) + successful_paths::Vector{Vector{Int}} # Successful paths from root to solution (variable sequence) +end + +DetailedStats() = DetailedStats( + Int[], + Int[], + 0, 0, 0, + 0.0, 0.0, 0.0, 0.0, + 0, 0, + Dict{Int, Int}(), + Int[], + Int[], + Int[], + Vector{Int}[] +) + +mutable struct BranchingStats + total_branches::Int + total_subproblems::Int + max_depth::Int + solved_leaves::Int + + detailed::Union{Nothing, DetailedStats} +end + +BranchingStats(verbose::Bool = false) = BranchingStats( + 0, 0, 0, 0, + verbose ? DetailedStats() : nothing +) + +function reset!(stats::BranchingStats) + stats.total_branches = 0 + stats.total_subproblems = 0 + stats.max_depth = 0 + stats.solved_leaves = 0 + + if !isnothing(stats.detailed) + d = stats.detailed + empty!(d.depth_distribution) + empty!(d.branching_factors) + d.propagation_calls = 0 + d.domain_reductions = 0 + d.early_unsat_detections = 0 + d.time_propagation = 0.0 + d.time_branching = 0.0 + d.time_contraction = 0.0 + d.time_filtering = 0.0 + d.cache_hits = 0 + d.cache_misses = 0 + empty!(d.variable_selection_counts) + empty!(d.remaining_vars_at_branch) + empty!(d.depth_at_selection) + empty!(d.variable_selection_sequence) + empty!(d.successful_paths) + end + + return stats +end + +function Base.copy(stats::BranchingStats) + detailed_copy = if stats.detailed !== nothing + d = stats.detailed + DetailedStats( + copy(d.depth_distribution), + copy(d.branching_factors), + d.propagation_calls, + d.domain_reductions, + d.early_unsat_detections, + d.time_propagation, + d.time_branching, + d.time_contraction, + d.time_filtering, + d.cache_hits, + d.cache_misses, + copy(d.variable_selection_counts), + copy(d.remaining_vars_at_branch), + copy(d.depth_at_selection), + copy(d.variable_selection_sequence), + [copy(path) for path in d.successful_paths] + ) + else + nothing + end + + return BranchingStats( + stats.total_branches, + stats.total_subproblems, + stats.max_depth, + stats.solved_leaves, + detailed_copy + ) +end + +function Base.getproperty(stats::BranchingStats, name::Symbol) + if name === :avg_branching_factor + return stats.total_branches > 0 ? + stats.total_subproblems / stats.total_branches : + 0.0 + else + return getfield(stats, name) + end +end + +function Base.propertynames(::BranchingStats; private::Bool=false) + return (fieldnames(BranchingStats)..., :avg_branching_factor) +end + +@inline function record_depth!(stats::BranchingStats, depth::Int) + depth > stats.max_depth && (stats.max_depth = depth) + return nothing +end + +@inline function record_branch!(stats::BranchingStats, subproblem_count::Int, depth::Int) + stats.total_branches += 1 + stats.total_subproblems += subproblem_count + record_depth!(stats, depth) + + # Record detailed statistics + if stats.detailed !== nothing + d = stats.detailed + push!(d.branching_factors, subproblem_count) + + # Update the depth distribution + while length(d.depth_distribution) <= depth + push!(d.depth_distribution, 0) + end + d.depth_distribution[depth + 1] += 1 + end + + return nothing +end + +@inline function record_unsat_leaf!(stats::BranchingStats, depth::Int) + record_depth!(stats, depth) + return nothing +end + +@inline function record_solved_leaf!(stats::BranchingStats, depth::Int, current_path::Vector{Int} = Int[]) + stats.solved_leaves += 1 + record_depth!(stats, depth) + if stats.detailed !== nothing && !isempty(current_path) + push!(stats.detailed.successful_paths, copy(current_path)) + end + return nothing +end + +@inline function record_skipped_subproblem!(stats::BranchingStats) + return nothing +end + +# === Detailed statistics recording functions === + +@inline function record_propagation!(stats::BranchingStats, time_elapsed::Float64 = 0.0) + if stats.detailed !== nothing + stats.detailed.propagation_calls += 1 + stats.detailed.time_propagation += time_elapsed + end + return nothing +end + +@inline function record_domain_reduction!(stats::BranchingStats, reduction_count::Int) + if stats.detailed !== nothing + stats.detailed.domain_reductions += reduction_count + end + return nothing +end + +@inline function record_early_unsat!(stats::BranchingStats) + if stats.detailed !== nothing + stats.detailed.early_unsat_detections += 1 + end + return nothing +end + +@inline function record_branching_time!(stats::BranchingStats, time_elapsed::Float64) + if stats.detailed !== nothing + stats.detailed.time_branching += time_elapsed + end + return nothing +end + +@inline function record_contraction_time!(stats::BranchingStats, time_elapsed::Float64) + if stats.detailed !== nothing + stats.detailed.time_contraction += time_elapsed + end + return nothing +end + +@inline function record_filtering_time!(stats::BranchingStats, time_elapsed::Float64) + if stats.detailed !== nothing + stats.detailed.time_filtering += time_elapsed + end + return nothing +end + +@inline function record_cache_hit!(stats::BranchingStats) + if stats.detailed !== nothing + stats.detailed.cache_hits += 1 + end + return nothing +end + +@inline function record_cache_miss!(stats::BranchingStats) + if stats.detailed !== nothing + stats.detailed.cache_misses += 1 + end + return nothing +end + +@inline function record_variable_selection!(stats::BranchingStats, var::Int, remaining_vars::Int, depth::Int) + if stats.detailed !== nothing + d = stats.detailed + d.variable_selection_counts[var] = get(d.variable_selection_counts, var, 0) + 1 + push!(d.remaining_vars_at_branch, remaining_vars) + push!(d.depth_at_selection, depth) + push!(d.variable_selection_sequence, var) + end + return nothing +end + +# Check if we need to track paths (only when detailed stats are enabled) +@inline needs_path_tracking(stats::BranchingStats) = !isnothing(stats.detailed) + +# === Statistics analysis and visualization === +function print_stats_summary(stats::BranchingStats; io::IO = stdout) + println(io, "=== Branching Statistics ===") + println(io, "Branch decisions: ", stats.total_branches) + println(io, "Child nodes: ", stats.total_subproblems) + println(io, "Max depth: ", stats.max_depth + 1) + println(io, "Solved leaves: ", stats.solved_leaves) + + if stats.detailed !== nothing + println(io, "\n=== Detailed Statistics ===") + d = stats.detailed + + # Propagation statistics + println(io, "\nPropagation:") + println(io, " Calls: ", d.propagation_calls) + println(io, " Domain reductions: ", d.domain_reductions) + println(io, " Early UNSAT detections: ", d.early_unsat_detections) + if d.propagation_calls > 0 + println(io, " Avg domains reduced per call: ", + round(d.domain_reductions / d.propagation_calls, digits=2)) + end + + # Timing statistics + total_time = d.time_propagation + d.time_branching + d.time_contraction + d.time_filtering + println(io, "\nTime breakdown:") + if total_time > 0 + println(io, " Propagation: ", round(d.time_propagation, digits=5), "s (", + round(100 * d.time_propagation / total_time, digits=5), "%)") + println(io, " Branching: ", round(d.time_branching, digits=5), "s (", + round(100 * d.time_branching / total_time, digits=5), "%)") + println(io, " Contraction: ", round(d.time_contraction, digits=5), "s (", + round(100 * d.time_contraction / total_time, digits=5), "%)") + println(io, " Filtering: ", round(d.time_filtering, digits=5), "s (", + round(100 * d.time_filtering / total_time, digits=5), "%)") + println(io, " Total: ", round(total_time, digits=5), "s") + else + println(io, " (Time statistics not recorded)") + end + + # Cache statistics + cache_total = d.cache_hits + d.cache_misses + if cache_total > 0 + println(io, "\nCache efficiency:") + println(io, " Hits: ", d.cache_hits) + println(io, " Misses: ", d.cache_misses) + println(io, " Hit rate: ", + round(100 * d.cache_hits / cache_total, digits=1), "%") + end + + # Branching factor distribution + if !isempty(d.branching_factors) + println(io, "\nBranching factor distribution:") + println(io, " Min: ", minimum(d.branching_factors)) + println(io, " Max: ", maximum(d.branching_factors)) + println(io, " Mean: ", round(sum(d.branching_factors) / length(d.branching_factors), digits=2)) + println(io, " Median: ", median(d.branching_factors)) + end + + # Depth distribution + if !isempty(d.depth_distribution) + println(io, "\nDepth distribution:") + # Compute cumulative node counts per depth to align with branching_factors + cumulative_nodes = 0 + for (depth_idx, node_count) in enumerate(d.depth_distribution) + if node_count > 0 + depth = depth_idx # Display depth starting from 1 + # Branching factor indices at this depth: cumulative_nodes+1 to cumulative_nodes+node_count + start_idx = cumulative_nodes + 1 + end_idx = cumulative_nodes + node_count + + if start_idx <= length(d.branching_factors) + end_idx = min(end_idx, length(d.branching_factors)) + depth_branching_factors = d.branching_factors[start_idx:end_idx] + avg_bf = round(sum(depth_branching_factors) / length(depth_branching_factors), digits=2) + total_subproblems = sum(depth_branching_factors) + println(io, " Depth ", depth, ": ", node_count, " nodes, ", + total_subproblems, " subproblems, avg branching factor: ", avg_bf) + else + println(io, " Depth ", depth, ": ", node_count, " nodes") + end + cumulative_nodes += node_count + end + end + end + + # Variable selection statistics + if !isempty(d.variable_selection_counts) + println(io, "\nVariable selection (top 10):") + sorted_vars = sort(collect(d.variable_selection_counts), by=x->x[2], rev=true) + # Compute the average remaining variable count and depth per variable + if !isempty(d.variable_selection_sequence) && + length(d.variable_selection_sequence) == length(d.remaining_vars_at_branch) && + length(d.variable_selection_sequence) == length(d.depth_at_selection) + for (var, count) in sorted_vars[1:min(10, length(sorted_vars))] + # Locate every position where this variable was selected + var_indices = [i for (i, v) in enumerate(d.variable_selection_sequence) if v == var] + avg_remaining = round(sum(d.remaining_vars_at_branch[i] for i in var_indices) / length(var_indices), digits=1) + avg_depth = round(sum(d.depth_at_selection[i] for i in var_indices) / length(var_indices), digits=1) + println(io, " Var ", var, ": ", count, " times, ", + "avg remaining vars: ", avg_remaining, ", ", + "avg depth: ", avg_depth + 1) + end + else + # Fallback display: show selection counts only + for (var, count) in sorted_vars[1:min(10, length(sorted_vars))] + println(io, " Var ", var, ": ", count, " times") + end + end + end + + if !isempty(d.successful_paths) + println(io, "\nSuccessful paths:") + for (path_idx, path) in enumerate(d.successful_paths) + println(io, " Path ", path_idx, ": ", join([string(v) for v in path], " -> ")) + end + end + end +end diff --git a/src/core/types.jl b/src/core/types.jl new file mode 100644 index 0000000..a4c8783 --- /dev/null +++ b/src/core/types.jl @@ -0,0 +1,63 @@ +struct Variable + id::Int + dom_size::Int + deg::Int +end + +function Base.show(io::IO, v::Variable) + print(io, "Variable($(v.id), dom_size=$(v.dom_size), deg=$(v.deg))") +end + +struct EdgeRef + var::Int + axis::Int +end + +function Base.show(io::IO, e::EdgeRef) + print(io, "EdgeRef($(e.var), axis=$(e.axis))") +end + +struct BoolTensor + id::Int + var_axes::Vector{Int} + tensor::Vector{Tropical{Float64}} +end + +function Base.show(io::IO, f::BoolTensor) + print(io, "BoolTensor($(f.id), vars=[$(join(f.var_axes, ", "))], size=$(length(f.tensor)))") +end + +struct TensorMasks + sat::BitVector + axis_masks0::Vector{BitVector} + axis_masks1::Vector{BitVector} +end + +# Build support masks for a tensor by enumerating all configurations. +function build_tensor_masks(tensor::BoolTensor) + nvars = length(tensor.var_axes) + n_cfg = 1 << nvars + + sat = falses(n_cfg) + axis_masks0 = [falses(n_cfg) for _ in 1:nvars] + axis_masks1 = [falses(n_cfg) for _ in 1:nvars] + + @inbounds for cfg in 0:(n_cfg-1) + # Check if this configuration satisfies the constraint + if tensor.tensor[cfg+1] == Tropical(0.0) # one(Tropical{Float64}) + sat[cfg+1] = true + end + + # Record which value each variable takes in this config + @inbounds for axis in 1:nvars + bit = (cfg >> (axis - 1)) & 1 + if bit == 0 + axis_masks0[axis][cfg+1] = true + else + axis_masks1[axis][cfg+1] = true + end + end + end + return TensorMasks(sat, axis_masks0, axis_masks1) +end + diff --git a/src/core/workspace.jl b/src/core/workspace.jl new file mode 100644 index 0000000..56a5aeb --- /dev/null +++ b/src/core/workspace.jl @@ -0,0 +1,49 @@ +mutable struct PropagationBuffers + feasible::BitVector + temp::BitVector + max_configs::Int +end + +function PropagationBuffers(static::TNStatic) + max_nvars = maximum(length(t.var_axes) for t in static.tensors; init=0) + max_configs = max_nvars > 0 ? (1 << max_nvars) : 1 + return PropagationBuffers(falses(max_configs), falses(max_configs), max_configs) +end + +mutable struct DynamicWorkspace + # Cache the full solution from the last branch for quick restoration + cached_doms::Vector{DomainMask} + has_cached_solution::Bool + # Branching statistics + branch_stats::BranchingStats + var_values::PriorityQueue{Int, Float64} + # O(1) membership test to avoid O(n) scans + changed_vars_flags::BitVector + changed_vars_indices::Vector{Int} + # Temporary BitVector cache for propagation to avoid reallocations + prop_buffers::Union{Nothing, PropagationBuffers} + # Cache of branch applications to avoid recomputing apply_branch + branch_cache::Dict{UInt, Dict{Tuple{UInt, Any}, Any}} + # Current search path (tracked only when verbose=true) + current_path::Vector{Int} +end + +DynamicWorkspace(var_num::Int, verbose::Bool = false) = DynamicWorkspace( + Vector{DomainMask}(undef, var_num), + false, + BranchingStats(verbose), + PriorityQueue{Int, Float64}(), + falses(var_num), + Int[], + nothing, + Dict{UInt, Dict{Tuple{UInt, Any}, Any}}(), + Int[] +) + +@inline function clear_branch_cache!(ws::DynamicWorkspace, doms_id::UInt) + inner = pop!(ws.branch_cache, doms_id, nothing) + if !isnothing(inner) + empty!(inner) + end + return nothing +end diff --git a/src/interface.jl b/src/interface.jl index 423c6de..b1dbf6d 100644 --- a/src/interface.jl +++ b/src/interface.jl @@ -1,36 +1,48 @@ -function setup_from_cnf(cnf::CNF) - return setup_from_sat(Satisfiability(cnf; use_constraints=true)) +# Performance Note: +# For maximum performance on large problems, set verbose=false and show_stats=false. +# This disables all detailed statistics collection (timing, path tracking, etc.) +# and only collects basic counters. + +function setup_from_cnf(cnf::CNF; verbose::Bool = false) + return setup_from_sat(Satisfiability(cnf; use_constraints=true); verbose=verbose) end -function setup_from_circuit(cir::Circuit) - return setup_from_sat(CircuitSAT(cir; use_constraints=true)) +function setup_from_circuit(cir::Circuit; verbose::Bool = false) + return setup_from_sat(CircuitSAT(cir; use_constraints=true); verbose=verbose) end -function setup_from_sat(sat::ConstraintSatisfactionProblem) +function setup_from_sat(sat::ConstraintSatisfactionProblem; verbose::Bool = false) tn = GenericTensorNetwork(sat) static = setup_from_tensor_network(tn) - TNProblem(static) + TNProblem(static; verbose=verbose) end -function solve(problem::TNProblem, bsconfig::BranchingStrategy, reducer::AbstractReducer) +function solve(problem::TNProblem, bsconfig::BranchingStrategy, reducer::AbstractReducer; show_stats::Bool=false) reset_branching_stats!(problem) # Reset stats before solving - depth = branch_and_reduce(problem, bsconfig, reducer, Tropical{Float64}; show_progress=false) + depth = OptimalBranchingCore.branch_and_reduce(problem, bsconfig, reducer, Tropical{Float64}; show_progress=false) res = last_branch_problem(problem) stats = get_branching_stats(problem) + clear_all_region_caches!() + if show_stats + print_stats_summary(stats) + end return (res, depth, stats) end function solve_sat_problem( sat::ConstraintSatisfactionProblem; bsconfig::BranchingStrategy=BranchingStrategy( - table_solver=TNContractionSolver(), - selector=LeastOccurrenceSelector(1, 2), + table_solver=TNContractionSolver(1,2), + selector=MostOccurrenceSelector(), measure=NumUnfixedVars() ), - reducer::AbstractReducer=NoReducer() + reducer::AbstractReducer=NoReducer(), + verbose::Bool = false, + show_stats::Bool=false ) - tn_problem = setup_from_sat(sat) - result, depth, stats = solve(tn_problem, bsconfig, reducer) + verbose = verbose || show_stats + tn_problem = setup_from_sat(sat; verbose=verbose) + result, depth, stats = solve(tn_problem, bsconfig, reducer; show_stats=show_stats) satisfiable = !isnothing(result) return satisfiable, result, depth, stats end @@ -38,13 +50,15 @@ end function solve_sat_with_assignments( sat::ConstraintSatisfactionProblem; bsconfig::BranchingStrategy=BranchingStrategy( - table_solver=TNContractionSolver(), - selector=LeastOccurrenceSelector(1, 2), + table_solver=TNContractionSolver(1,2), + selector=MostOccurrenceSelector(), measure=NumUnfixedVars() - ), - reducer::AbstractReducer=NoReducer() + ), + reducer::AbstractReducer=NoReducer(), + verbose::Bool = false, + show_stats::Bool=false ) - satisfiable, result, depth, stats = solve_sat_problem(sat; bsconfig, reducer) + satisfiable, result, depth, stats = solve_sat_problem(sat; bsconfig, reducer, verbose=verbose, show_stats=show_stats) if satisfiable && !isnothing(result) # Convert TNProblem result to variable assignments assignments = Dict{Symbol, Int}() @@ -59,16 +73,41 @@ end function solve_factoring( n::Int, m::Int, N::Int; - bsconfig::BranchingStrategy=BranchingStrategy(table_solver=TNContractionSolver(), selector=LeastOccurrenceSelector(1, 2), measure=NumUnfixedVars()), - reducer::AbstractReducer=NoReducer() + bsconfig::BranchingStrategy=BranchingStrategy( + table_solver=TNContractionSolver(1,2), + # selector=MinGammaSelector(TNContractionSolver(1,2), GreedyMerge()), + selector=MostOccurrenceSelector(), + measure=NumUnfixedVars(), + set_cover_solver=GreedyMerge() + ), + reducer::AbstractReducer=NoReducer(), + verbose::Bool=false, + show_stats::Bool=false ) + verbose = verbose || show_stats fproblem = Factoring(n, m, N) circuit_sat = reduceto(CircuitSAT, fproblem) problem = CircuitSAT(circuit_sat.circuit.circuit; use_constraints=true) - tn_problem = setup_from_sat(problem) - res, _, stats = solve(tn_problem, bsconfig, reducer) + tn_problem = setup_from_sat(problem; verbose=verbose) + if verbose + unique_tensors = length(tn_problem.static.precomputed_masks) + total_tensors = length(tn_problem.static.tensors) + println("Unique tensor types: $unique_tensors (out of $total_tensors total)") + end + res, _, stats = solve(tn_problem, bsconfig, reducer; show_stats=show_stats) isnothing(res) && return nothing, nothing, stats a = get_var_value(res, circuit_sat.q) b = get_var_value(res, circuit_sat.p) return bits_to_int(a), bits_to_int(b), stats end + + +# function solve_circuit_sat( +# circuit::CircuitSAT; +# bsconfig::BranchingStrategy=BranchingStrategy(table_solver=TNContractionSolver(), selector=LeastOccurrenceSelector(1, 2), measure=NumUnfixedVars()), +# reducer::AbstractReducer=NoReducer() +# ) +# tn_problem = setup_from_circuit(circuit.circuit) +# res, _, stats = solve(tn_problem, bsconfig, reducer) +# return res, stats +# end \ No newline at end of file diff --git a/src/problems.jl b/src/problems.jl deleted file mode 100644 index 8bf45e7..0000000 --- a/src/problems.jl +++ /dev/null @@ -1,198 +0,0 @@ -struct Variable - id::Int - dom_size::Int - deg::Int -end -function Base.show(io::IO, v::Variable) - print(io, "Variable($(v.id), dom_size=$(v.dom_size), deg=$(v.deg))") -end - -struct EdgeRef - var::Int # variable id - axis::Int # inside a factor -end -function Base.show(io::IO, e::EdgeRef) - print(io, "EdgeRef($(e.var), axis=$(e.axis))") -end - -struct BoolTensor - id::Int - var_axes::Vector{Int} - tensor::Vector{Tropical{Float64}} -end - -function Base.show(io::IO, f::BoolTensor) - print(io, "BoolTensor($(f.id), vars=[$(join(f.var_axes, ", "))], size=$(length(f.tensor)))") -end - -struct TNStatic - vars::Vector{Variable} - tensors::Vector{BoolTensor} - v2t::Vector{Vector{Int}} # var id -> tensor ids - t2v::Vector{Vector{EdgeRef}} - axis_of_t::Dict{Tuple{Int,Int},Int} -end - -function Base.show(io::IO, tn::TNStatic) - print(io, "TNStatic(vars=$(length(tn.vars)), tensors=$(length(tn.tensors)))") -end - -function Base.show(io::IO, ::MIME"text/plain", tn::TNStatic) - println(io, "TNStatic:") - println(io, " variables: $(length(tn.vars))") - println(io, " tensors: $(length(tn.tensors))") - println(io, " variable-tensor incidence: $(length(tn.v2t))") - println(io, " tensor-variable incidence: $(length(tn.t2v))") - println(io, " axis mappings: $(length(tn.axis_of_t))") -end - -function setup_problem(var_num::Int, - tensors_to_vars::Vector{Vector{Int}}, - tensor_data::Vector{Vector{Tropical{Float64}}}) - F = length(tensors_to_vars) - tensors = Vector{BoolTensor}(undef, F) - vars_to_tensors = [Int[] for _ in 1:var_num] - tensors_to_edges = [EdgeRef[] for _ in 1:F] - for i in 1:F - var_axes = tensors_to_vars[i] - @assert length(tensor_data[i]) == 1 << length(var_axes) "Boolean tensor length mismatch" - tensors[i] = BoolTensor(i, var_axes, tensor_data[i]) - for (j, v) in enumerate(var_axes) - push!(vars_to_tensors[v], i) - push!(tensors_to_edges[i], EdgeRef(v, j)) - end - end - - axis_of_t = Dict{Tuple{Int,Int},Int}() - for (fid, tensor) in enumerate(tensors) - for (axis, vid) in enumerate(tensor.var_axes) - axis_of_t[(fid, vid)] = axis - end - end - - vars = Vector{Variable}(undef, var_num) - for i in 1:var_num - vars[i] = Variable(i, 2, length(vars_to_tensors[i])) - end - - return TNStatic(vars, tensors, vars_to_tensors, tensors_to_edges, axis_of_t) -end - -function setup_from_tensor_network(tn)::TNStatic - t2v = getixsv(tn.code) - tensors = GenericTensorNetworks.generate_tensors(Tropical(1.0), tn) - vec_tensors = [vec(t) for t in tensors] - new_tensors = [replace(t, Tropical(1.0) => zero(Tropical{Float64})) for t in vec_tensors] - return setup_problem(length(tn.problem.symbols), t2v, new_tensors) -end - -struct DomainMask - bits::UInt8 -end -const DM_BOTH = DomainMask(0x03) -const DM_0 = DomainMask(0x01) -const DM_1 = DomainMask(0x02) - -@inline is_fixed(dm::DomainMask) = (dm.bits == 0x01) | (dm.bits == 0x02) -@inline has0(dm::DomainMask)::Bool = (dm.bits & 0x01) != 0 -@inline has1(dm::DomainMask)::Bool = (dm.bits & 0x02) != 0 - -function init_doms(static::TNStatic) - return fill(DM_BOTH, length(static.vars)) -end - -@inline get_var_value(dms::Vector{DomainMask}, var_id::Int) = has1(dms[var_id]) ? 1 : 0 - -mutable struct DynamicWorkspace - cached_doms::Vector{DomainMask} - has_cached_solution::Bool - # Branching statistics - total_branches::Int - total_subproblems::Int - max_depth::Int -end -DynamicWorkspace(var_num::Int) = DynamicWorkspace(Vector{DomainMask}(undef, var_num), false, 0, 0, 0) - -struct Region - id::Int - tensors::Vector{Int} - inner_vars::Vector{Int} - boundary_vars::Vector{Int} -end -function Base.show(io::IO, region::Region) - print(io, "Region(focus=$(region.id), tensors=$(length(region.tensors)), inner=$(length(region.inner_vars)), boundary=$(length(region.boundary_vars)))") -end - -struct TNProblem <: AbstractProblem - static::TNStatic - doms::Vector{DomainMask} - n_unfixed::Int - ws::DynamicWorkspace -end -function TNProblem(static::TNStatic)::TNProblem - doms = init_doms(static) - doms = propagate(static, doms) - n_unfixed = count_unfixed(doms) - ws = DynamicWorkspace(length(static.vars)) - return TNProblem(static, doms, n_unfixed, ws) -end -function Base.show(io::IO, problem::TNProblem) - print(io, "TNProblem(unfixed=$(problem.n_unfixed)/$(length(problem.static.vars)))") -end - -get_var_value(problem::TNProblem, var_id::Int) = get_var_value(problem.doms, var_id) -get_var_value(problem::TNProblem, var_ids::Vector{Int}) = Bool[get_var_value(problem.doms, var_id) for var_id in var_ids] - -@inline function is_solved(problem::TNProblem)::Bool - return problem.n_unfixed == 0 -end - -function cache_branch_solution!(problem::TNProblem) - ws = problem.ws - copyto!(ws.cached_doms, problem.doms) - ws.has_cached_solution = true - return nothing -end - -function reset_last_branch_problem!(problem::TNProblem) - problem.ws.has_cached_solution = false - return nothing -end - -@inline has_last_branch_problem(problem::TNProblem) = problem.ws.has_cached_solution - -function last_branch_problem(problem::TNProblem) - has_last_branch_problem(problem) || return nothing - doms = copy(problem.ws.cached_doms) - fixed = count(x -> is_fixed(x), doms) - @assert fixed == length(doms) - return TNProblem(problem.static, doms, 0, problem.ws) -end - -# Branching statistics functions -function get_branching_stats(problem::TNProblem) - return ( - total_branches = problem.ws.total_branches, - total_subproblems = problem.ws.total_subproblems, - max_depth = problem.ws.max_depth, - avg_branching_factor = problem.ws.total_branches > 0 ? - problem.ws.total_subproblems / problem.ws.total_branches : 0.0 - ) -end - -function reset_branching_stats!(problem::TNProblem) - problem.ws.total_branches = 0 - problem.ws.total_subproblems = 0 - problem.ws.max_depth = 0 - return nothing -end - -function print_branching_stats(io::IO, problem::TNProblem) - stats = get_branching_stats(problem) - println(io, "Branching Statistics:") - println(io, " Total branch points: $(stats.total_branches)") - println(io, " Total subproblems: $(stats.total_subproblems)") - println(io, " Maximum depth: $(stats.max_depth)") - println(io, " Average branching factor: $(round(stats.avg_branching_factor, digits=3))") -end -print_branching_stats(problem::TNProblem) = print_branching_stats(stdout, problem) diff --git a/src/propagate.jl b/src/propagate.jl deleted file mode 100644 index c0d27fc..0000000 --- a/src/propagate.jl +++ /dev/null @@ -1,330 +0,0 @@ -""" - TensorMasks - -Precomputed bitmasks for efficient propagation: -- `sat`: which configurations satisfy the tensor constraint -- `axis_masks0[i]`: which configurations have variable i = 0 -- `axis_masks1[i]`: which configurations have variable i = 1 -""" -struct TensorMasks - sat::BitVector - axis_masks0::Vector{BitVector} - axis_masks1::Vector{BitVector} -end - -struct PropagationBuffers - feasible::BitVector - allowed_axis::BitVector - temp::BitVector -end - -function PropagationBuffers(static::TNStatic) - max_nvars = maximum(length(t.var_axes) for t in static.tensors; init=0) - max_configs = max_nvars > 0 ? (1 << max_nvars) : 1 - return PropagationBuffers( - falses(max_configs), - falses(max_configs), - falses(max_configs) - ) -end - - -function get_active_tensors(static::TNStatic, doms::Vector{DomainMask}) - active = Int[] - sizehint!(active, length(static.tensors)) # Pre-allocate to avoid resizing - - @inbounds for (tid, tensor) in enumerate(static.tensors) - # Check if any variable in this tensor is unfixed - has_unfixed = false - for var_id in tensor.var_axes - dm_bits = doms[var_id].bits - # Unfixed if domain allows both 0 and 1, or is not yet determined - if dm_bits == 0x03 # DM_BOTH - has_unfixed = true - break # Early exit: found one unfixed variable - end - end - - has_unfixed && push!(active, tid) - end - return active -end - -function propagate(static::TNStatic, doms::Vector{DomainMask}) - working_doms = copy(doms) - masks_cache = _masks_cache!(static) - - # Initialize propagation queue with only active tensors (containing unfixed variables) - tensor_queue = get_active_tensors(static, doms) - in_queue = falses(length(static.tensors)) - @inbounds for tid in tensor_queue - in_queue[tid] = true - end - queue_pos = 1 - - # Pre-allocate buffers to reduce allocations during propagation - buffers = PropagationBuffers(static) - - while queue_pos <= length(tensor_queue) - tensor_idx = tensor_queue[queue_pos] - queue_pos += 1 - in_queue[tensor_idx] = false - - tensor = static.tensors[tensor_idx] - masks = get_mask!(masks_cache, tensor_idx, tensor) - - # Check constraint and propagate - if !propagate_tensor!(working_doms, tensor, masks, static.v2t, - tensor_queue, in_queue, buffers) - # Contradiction detected - return fill(DomainMask(0x00), length(working_doms)) - end - end - - return working_doms -end - - -# Propagate constraints from a single tensor. Returns false if contradiction detected. -function propagate_tensor!( - working_doms::Vector{DomainMask}, - tensor::BoolTensor, - masks::TensorMasks, - v2t::Vector{Vector{Int}}, - tensor_queue::Vector{Int}, - in_queue::BitVector, - buffers::PropagationBuffers -) - # Step 1: Compute feasible configurations given current domains - if !compute_feasible_configs!(buffers.feasible, buffers.allowed_axis, working_doms, tensor, masks) - return false # No feasible config -> contradiction - end - - n_feasible = count(buffers.feasible) - - # Step 2: Different propagation strategies based on number of feasible configs - if n_feasible == 1 - # Only one valid config -> fix all variables to that config - return propagate_unit_constraint!(working_doms, tensor, buffers.feasible, v2t, tensor_queue, in_queue) - else - # Multiple configs -> prune unsupported values - return propagate_support_pruning!(working_doms, tensor, masks, - buffers.feasible, buffers.temp, - v2t, tensor_queue, in_queue) - end -end - - -# Compute which tensor configurations are feasible given current variable domains. -# Returns false if no feasible configuration exists. -@inline function compute_feasible_configs!(feasible::BitVector, allowed_axis::BitVector, working_doms::Vector{DomainMask}, tensor::BoolTensor, masks::TensorMasks) - n_cfg = length(masks.sat) - - # Resize buffers to match this tensor's configuration space - resize!(feasible, n_cfg) - resize!(allowed_axis, n_cfg) - - # Start with satisfying configurations - copyto!(feasible, masks.sat) - - # Filter by each variable's domain - @inbounds for (axis, var_id) in enumerate(tensor.var_axes) - dm = working_doms[var_id] - dm_bits = dm.bits - - # Skip filtering if domain is unrestricted (both 0 and 1 allowed) - if dm_bits == 0x03 # DM_BOTH - continue - end - - # Directly intersect with the appropriate mask(s) - if dm_bits == 0x01 # DM_0: only 0 allowed - feasible .&= masks.axis_masks0[axis] - elseif dm_bits == 0x02 # DM_1: only 1 allowed - feasible .&= masks.axis_masks1[axis] - else # dm_bits == 0x00: contradiction, but shouldn't happen - return false - end - - # Early exit if no feasible config remains - any(feasible) || return false - end - - return true -end - -# When tensor has exactly one feasible configuration, fix all its variables to that config. -@inline function propagate_unit_constraint!( - working_doms::Vector{DomainMask}, - tensor::BoolTensor, - feasible::BitVector, - v2t::Vector{Vector{Int}}, - tensor_queue::Vector{Int}, - in_queue::BitVector -) - first_idx = findfirst(feasible) - config = first_idx - 1 - - @inbounds for (axis, var_id) in enumerate(tensor.var_axes) - bit_val = (config >> (axis - 1)) & 1 - # Branchless: avoid conditional for better performance - required_bits = ifelse(bit_val == 1, DM_1.bits, DM_0.bits) - - old_bits = working_doms[var_id].bits - new_bits = old_bits & required_bits - - # Check for contradiction - new_bits == 0x00 && return false - - # Update domain and enqueue affected tensors - if new_bits != old_bits - working_doms[var_id] = DomainMask(new_bits) - enqueue_affected_tensors!(tensor_queue, in_queue, v2t, var_id) - end - end - - return true -end - -# Prune domain values that have no support in any feasible configuration. -@inline function propagate_support_pruning!( - working_doms::Vector{DomainMask}, - tensor::BoolTensor, - masks::TensorMasks, - feasible::BitVector, - temp::BitVector, - v2t::Vector{Vector{Int}}, - tensor_queue::Vector{Int}, - in_queue::BitVector -) - n_cfg = length(feasible) - resize!(temp, n_cfg) - - @inbounds for (axis, var_id) in enumerate(tensor.var_axes) - dm = working_doms[var_id] - dm_bits = dm.bits - - # Check support for value 0 (avoid temp array allocation) - if (dm_bits & 0x01) != 0 # has0(dm) - # Use copyto! + .&= to avoid RHS temporary array - copyto!(temp, feasible) - temp .&= masks.axis_masks0[axis] - if !any(temp) # No support for 0 - if !update_domain!(working_doms, var_id, dm_bits, DM_1.bits, - v2t, tensor_queue, in_queue) - return false - end - end - end - - # Check support for value 1 - if (dm_bits & 0x02) != 0 # has1(dm) - copyto!(temp, feasible) - temp .&= masks.axis_masks1[axis] - if !any(temp) # No support for 1 - if !update_domain!(working_doms, var_id, dm_bits, DM_0.bits, - v2t, tensor_queue, in_queue) - return false - end - end - end - end - - return true -end - -# Add all tensors containing the variable to the propagation queue. -@inline function enqueue_affected_tensors!( - tensor_queue::Vector{Int}, - in_queue::BitVector, - v2t::Vector{Vector{Int}}, - var_id::Int -) - @inbounds for tensor_id in v2t[var_id] - if !in_queue[tensor_id] - push!(tensor_queue, tensor_id) - in_queue[tensor_id] = true - end - end -end - -# Update variable domain by intersecting with keep_mask. Returns false if contradiction. -@inline function update_domain!( - working_doms::Vector{DomainMask}, - var_id::Int, - old_bits::UInt8, - keep_bits::UInt8, - v2t::Vector{Vector{Int}}, - tensor_queue::Vector{Int}, - in_queue::BitVector -) - new_bits = old_bits & keep_bits - - # Check for contradiction - new_bits == 0x00 && return false - - # Update if changed - if new_bits != old_bits - working_doms[var_id] = DomainMask(new_bits) - enqueue_affected_tensors!(tensor_queue, in_queue, v2t, var_id) - end - - return true -end - -# Build support masks for a tensor by enumerating all configurations. -function build_tensor_masks(tensor::BoolTensor) - nvars = length(tensor.var_axes) - n_cfg = 1 << nvars - - sat = falses(n_cfg) - axis_masks0 = [falses(n_cfg) for _ in 1:nvars] - axis_masks1 = [falses(n_cfg) for _ in 1:nvars] - - @inbounds for cfg in 0:(n_cfg-1) - # Check if this configuration satisfies the constraint - if tensor.tensor[cfg+1] == Tropical(0.0) # one(Tropical{Float64}) - sat[cfg+1] = true - end - - # Record which value each variable takes in this config - @inbounds for axis in 1:nvars - bit = (cfg >> (axis - 1)) & 1 - if bit == 0 - axis_masks0[axis][cfg+1] = true - else - axis_masks1[axis][cfg+1] = true - end - end - end - return TensorMasks(sat, axis_masks0, axis_masks1) -end - - -@inline function get_mask!(cache::Vector{Union{Nothing, TensorMasks}}, idx::Int, tensor::BoolTensor) - m = cache[idx] - isnothing(m) && (m = build_tensor_masks(tensor); cache[idx] = m) - return m -end - - -const _MASKS_CACHE = IdDict{UInt, Vector{Union{Nothing, TensorMasks}}}() - -function _masks_cache!(static::TNStatic) - key = objectid(static) - vec = get!(_MASKS_CACHE, key) do - v = Vector{Union{Nothing, TensorMasks}}(undef, length(static.tensors)) - fill!(v, nothing) - v - end - - # Resize cache if problem structure changed - if length(vec) != length(static.tensors) - v = Vector{Union{Nothing, TensorMasks}}(undef, length(static.tensors)) - fill!(v, nothing) - _MASKS_CACHE[key] = v - return v - end - - return vec -end diff --git a/src/region.jl b/src/region.jl deleted file mode 100644 index acdcbbe..0000000 --- a/src/region.jl +++ /dev/null @@ -1,80 +0,0 @@ -struct RegionContraction # tensor and vars that are contracted to form the region - tensor::AbstractArray{Tropical{Float64}} - vars::Vector{Int} -end - -mutable struct RegionCacheEntry # one single region and its contraction - region::Region - contraction::Union{Nothing, RegionContraction} -end -RegionCacheEntry(region::Region) = RegionCacheEntry(region, nothing) - -mutable struct RegionCacheState # a collection of region entries - entries::Dict{Int, RegionCacheEntry} # region id -> region entry - last_region_id::Union{Nothing, Int} -end -RegionCacheState() = RegionCacheState(Dict{Int, RegionCacheEntry}(), nothing) - -struct RegionCache # a collection of region cache states - data::IdDict{UInt, RegionCacheState} -end - -const REGION_CACHE = RegionCache(IdDict{UInt, RegionCacheState}()) - -function cache_region!(problem::TNProblem, region::Region) - # get or create the region cache state for the problem - state = get!(REGION_CACHE.data, objectid(problem)) do - RegionCacheState() - end - state.entries[region.id] = RegionCacheEntry(region) - state.last_region_id = region.id - return nothing -end - -function get_cached_region_entry(problem::TNProblem, region_id::Int) - state = get(REGION_CACHE.data, objectid(problem), nothing) - isnothing(state) && return nothing - get(state.entries, region_id, nothing) -end -function get_cached_region_entry(problem::TNProblem) - # get the last region id from the state - state = get(REGION_CACHE.data, objectid(problem), nothing) - (isnothing(state) || isnothing(state.last_region_id)) && return nothing - get(state.entries, state.last_region_id, nothing) -end - -function get_cached_region(problem::TNProblem) - entry = get_cached_region_entry(problem) - isnothing(entry) && return nothing - entry.region -end - -function cache_region_contraction!(problem::TNProblem, tensor::AbstractArray{Tropical{Float64}}, vars::Vector{Int}; region_id::Union{Nothing, Int}=nothing) - entry = isnothing(region_id) ? get_cached_region_entry(problem) : get_cached_region_entry(problem, region_id) - isnothing(entry) && error("No cached region entry available for contraction caching.") - entry.contraction = RegionContraction(tensor, vars) - return entry.contraction -end - -function get_cached_region_contraction( - problem::TNProblem; - region_id::Union{Nothing,Int}=nothing, -) - entry = isnothing(region_id) ? get_cached_region_entry(problem) : - get_cached_region_entry(problem, region_id) - isnothing(entry) && return nothing - isnothing(entry.contraction) && return nothing - contraction = entry.contraction - return contraction.tensor, contraction.vars -end - -function clear_region_cache!(problem::TNProblem) - delete!(REGION_CACHE.data, objectid(problem)) - return nothing -end - -function clear_all_region_caches!() - empty!(REGION_CACHE.data) - return nothing -end - diff --git a/src/selector.jl b/src/selector.jl deleted file mode 100644 index ada49c2..0000000 --- a/src/selector.jl +++ /dev/null @@ -1,32 +0,0 @@ -struct LeastOccurrenceSelector <: AbstractSelector - k::Int # number of hops - max_tensors::Int -end - -LeastOccurrenceSelector(k::Int) = LeastOccurrenceSelector(k, 2) - -function OptimalBranchingCore.select_variables( - problem::TNProblem, - measure::AbstractMeasure, - selector::LeastOccurrenceSelector -) - unfixed_vars = get_unfixed_vars(problem.doms) - isempty(unfixed_vars) && return Int[] - least_show_var_idx = argmin(length(problem.static.v2t[u]) for u in unfixed_vars) - - # Compute k-neighboring region (only expanding to unfixed variables) - # The returned region will only contain unfixed variables - region = k_neighboring( - problem.static, - problem.doms, - unfixed_vars[least_show_var_idx]; - max_tensors = selector.max_tensors, - k = selector.k - ) - cache_region!(problem, region) - - # Return all variables in the region (boundary first, then inner) - # These are guaranteed to be unfixed since k_neighboring filters them - return vcat(region.boundary_vars, region.inner_vars) -end - diff --git a/src/utils.jl b/src/utils/utils.jl similarity index 100% rename from src/utils.jl rename to src/utils/utils.jl diff --git a/test/branch.jl b/test/branch.jl index a8dc47f..728998e 100644 --- a/test/branch.jl +++ b/test/branch.jl @@ -1,6 +1,6 @@ using Test using BooleanInference -using BooleanInference: TNProblem, TNContractionSolver, LeastOccurrenceSelector, NumUnfixedVars, setup_from_tensor_network, setup_problem, select_variables, get_cached_region, last_branch_problem, reset_last_branch_problem!, has_last_branch_problem, get_var_value, bits_to_int +using BooleanInference: TNProblem, TNContractionSolver, MostOccurrenceSelector, NumUnfixedVars, setup_from_tensor_network, setup_problem, select_variables, get_cached_region, last_branch_problem, reset_last_branch_problem!, has_last_branch_problem, get_var_value, bits_to_int using BooleanInference: BranchingStrategy, NoReducer using OptimalBranchingCore: Clause using OptimalBranchingCore: branching_table, branch_and_reduce, apply_branch @@ -12,8 +12,9 @@ using TropicalNumbers: Tropical # ENV["JULIA_DEBUG"] = "BooleanInference" @testset "branch" begin - fproblem = Factoring(12, 12, 10395529) - # fproblem = Factoring(2, 2, 4) + # fproblem = Factoring(16, 16, 3363471157) + fproblem = Factoring(10, 10, 559619) + circuit_sat = reduceto(CircuitSAT, fproblem) problem = CircuitSAT(circuit_sat.circuit.circuit; use_constraints=true) @@ -21,8 +22,9 @@ using TropicalNumbers: Tropical tn_static = setup_from_tensor_network(tn) tn_problem = TNProblem(tn_static) @test !has_last_branch_problem(tn_problem) - br_strategy = BranchingStrategy(table_solver = TNContractionSolver(), selector = LeastOccurrenceSelector(1, 2), measure = NumUnfixedVars()) - @time res = branch_and_reduce(tn_problem, br_strategy, NoReducer(), Tropical{Float64}; show_progress=false) + br_strategy = BranchingStrategy(table_solver = TNContractionSolver(1,2), selector = MostOccurrenceSelector(), measure = NumUnfixedVars()) + res = branch_and_reduce(tn_problem, br_strategy, NoReducer(), Tropical{Float64}; show_progress=false) + @show res if res != Tropical(-Inf) @test has_last_branch_problem(tn_problem) res = last_branch_problem(tn_problem) @@ -34,6 +36,7 @@ using TropicalNumbers: Tropical # @show a, b @show bits_to_int(a), bits_to_int(b) end + clear_all_region_caches!() end @testset "apply_branch" begin diff --git a/test/contraction.jl b/test/contraction.jl index 211c091..72cc971 100644 --- a/test/contraction.jl +++ b/test/contraction.jl @@ -1,17 +1,15 @@ using Test using BooleanInference -using BooleanInference: setup_from_tensor_network, TNProblem, DynamicWorkspace, setup_problem, select_variables, get_cached_region, LeastOccurrenceSelector, NumUnfixedVars +using BooleanInference: setup_from_tensor_network, TNProblem, DynamicWorkspace, setup_problem, select_variables, get_cached_region, MostOccurrenceSelector, NumUnfixedVars using BooleanInference: Region, slicing, tensor_unwrapping, DomainMask using BooleanInference: DM_BOTH, DM_0, DM_1, has0, has1, is_fixed using BooleanInference: contract_tensors, contract_region, TNContractionSolver -using BooleanInference: separate_fixed_free_boundary, construct_boundary_config, construct_inner_config -using BooleanInference: extract_inner_configs, combine_configs, handle_no_boundary_case -using BooleanInference: get_cached_region_contraction, clear_region_cache! +using BooleanInference: extract_inner_configs, combine_configs, handle_no_boundary_case_unfixed +using BooleanInference: clear_all_region_caches! using OptimalBranchingCore: branching_table using TropicalNumbers: Tropical using ProblemReductions: Factoring, reduceto, CircuitSAT using GenericTensorNetworks -using BooleanInference: separate_fixed_free_boundary, construct_boundary_config, construct_inner_config, extract_inner_configs, combine_configs @testset "Region constructor" begin region = Region(1, @@ -218,197 +216,212 @@ end tn = GenericTensorNetwork(problem) tn_static = setup_from_tensor_network(tn) tn_problem = TNProblem(tn_static) - select_variables(tn_problem, NumUnfixedVars(), LeastOccurrenceSelector(2, 5)) - region = get_cached_region(tn_problem) - @show region - @test region != nothing - @test isnothing(get_cached_region_contraction(tn_problem)) + variable = select_variables(tn_problem, NumUnfixedVars(), MostOccurrenceSelector()) + @test variable isa Int - contracted, _ = contract_region(tn_static, region, tn_problem.doms) - @test contracted != nothing - @test length(size(contracted)) == length(region.boundary_vars) + length(region.inner_vars) + # Create region for the selected variable + solver = TNContractionSolver(1, 2) + all_unfixed_doms = fill(DM_BOTH, length(tn_problem.doms)) + region = k_neighboring(tn_static, all_unfixed_doms, variable; max_tensors=solver.max_tensors, k=solver.k) + @test region !== nothing - table = branching_table(tn_problem, TNContractionSolver(), vcat([v for v in region.boundary_vars], [v for v in region.inner_vars])) + contracted, output_vars = contract_region(tn_static, region, tn_problem.doms) + @test contracted !== nothing + @test length(output_vars) <= length(region.boundary_vars) + length(region.inner_vars) + + table = branching_table(tn_problem, TNContractionSolver(), region.id) @show vcat([v for v in region.boundary_vars], [v for v in region.inner_vars]) @show table - - cached_contraction = get_cached_region_contraction(tn_problem) - @test !isnothing(cached_contraction) - - clear_region_cache!(tn_problem) - @test isnothing(get_cached_region(tn_problem)) - @test isnothing(get_cached_region_contraction(tn_problem)) end -@testset "separate_fixed_free_boundary" begin - region = Region(1, - [1, 2], - [5, 6], - [3, 4]) - - doms = fill(DM_BOTH, 10) - fixed, fixed_vals, free, free_indices = separate_fixed_free_boundary(region, doms) - @test isempty(fixed) - @test isempty(fixed_vals) - @test length(free) == 2 - @test free == Int[3, 4] - @test free_indices == [1, 2] - - doms[3] = DM_0 - doms[4] = DM_1 - fixed, fixed_vals, free, free_indices = separate_fixed_free_boundary(region, doms) - @test length(fixed) == 2 - @test fixed == Int[3, 4] - @test fixed_vals == Int[0, 1] - @test isempty(free) - @test isempty(free_indices) - - doms[3] = DM_0 - doms[4] = DM_BOTH - fixed, fixed_vals, free, free_indices = separate_fixed_free_boundary(region, doms) - @test length(fixed) == 1 - @test fixed == Int[3] - @test fixed_vals == Int[0] - @test length(free) == 1 - @test free == Int[4] - @test free_indices == [2] -end +# These test sets are commented out because the functions they test +# (separate_fixed_free_boundary, construct_boundary_config, construct_inner_config) +# are no longer in the codebase - they appear to have been refactored. +# The functionality is now handled internally by other functions. +# @testset "separate_fixed_free_boundary" begin +# region = Region(1, +# [1, 2], +# [5, 6], +# [3, 4]) +# +# doms = fill(DM_BOTH, 10) +# fixed, fixed_vals, free, free_indices = separate_fixed_free_boundary(region, doms) +# @test isempty(fixed) +# @test isempty(fixed_vals) +# @test length(free) == 2 +# @test free == Int[3, 4] +# @test free_indices == [1, 2] +# +# doms[3] = DM_0 +# doms[4] = DM_1 +# fixed, fixed_vals, free, free_indices = separate_fixed_free_boundary(region, doms) +# @test length(fixed) == 2 +# @test fixed == Int[3, 4] +# @test fixed_vals == Int[0, 1] +# @test isempty(free) +# @test isempty(free_indices) +# +# doms[3] = DM_0 +# doms[4] = DM_BOTH +# fixed, fixed_vals, free, free_indices = separate_fixed_free_boundary(region, doms) +# @test length(fixed) == 1 +# @test fixed == Int[3] +# @test fixed_vals == Int[0] +# @test length(free) == 1 +# @test free == Int[4] +# @test free_indices == [2] +# end -@testset "construct_boundary_config" begin - region = Region(1, - [1], - [5], - [3, 4]) - - doms = fill(DM_BOTH, 10) - free_boundary_indices = [1, 2] - # Build free_pos_of_boundary: maps boundary position to free position - n_boundary = length(region.boundary_vars) - free_pos_of_boundary = zeros(Int, n_boundary) - for (j, boundary_idx) in enumerate(free_boundary_indices) - free_pos_of_boundary[boundary_idx] = j - end - - config = construct_boundary_config(region, doms, free_pos_of_boundary, 0) - @test config == [false, false] - - config = construct_boundary_config(region, doms, free_pos_of_boundary, 1) - @test config == [true, false] - - config = construct_boundary_config(region, doms, free_pos_of_boundary, 2) - @test config == [false, true] - - config = construct_boundary_config(region, doms, free_pos_of_boundary, 3) - @test config == [true, true] - - doms[3] = DM_1 - free_boundary_indices = [2] - # Rebuild free_pos_of_boundary with new free_boundary_indices - free_pos_of_boundary = zeros(Int, n_boundary) - for (j, boundary_idx) in enumerate(free_boundary_indices) - free_pos_of_boundary[boundary_idx] = j - end - - config = construct_boundary_config(region, doms, free_pos_of_boundary, 0) - @test config == [true, false] - - config = construct_boundary_config(region, doms, free_pos_of_boundary, 1) - @test config == [true, true] -end +# @testset "construct_boundary_config" begin +# region = Region(1, +# [1], +# [5], +# [3, 4]) +# +# doms = fill(DM_BOTH, 10) +# free_boundary_indices = [1, 2] +# # Build free_pos_of_boundary: maps boundary position to free position +# n_boundary = length(region.boundary_vars) +# free_pos_of_boundary = zeros(Int, n_boundary) +# for (j, boundary_idx) in enumerate(free_boundary_indices) +# free_pos_of_boundary[boundary_idx] = j +# end +# +# config = construct_boundary_config(region, doms, free_pos_of_boundary, 0) +# @test config == [false, false] +# +# config = construct_boundary_config(region, doms, free_pos_of_boundary, 1) +# @test config == [true, false] +# +# config = construct_boundary_config(region, doms, free_pos_of_boundary, 2) +# @test config == [false, true] +# +# config = construct_boundary_config(region, doms, free_pos_of_boundary, 3) +# @test config == [true, true] +# +# doms[3] = DM_1 +# free_boundary_indices = [2] +# # Rebuild free_pos_of_boundary with new free_boundary_indices +# free_pos_of_boundary = zeros(Int, n_boundary) +# for (j, boundary_idx) in enumerate(free_boundary_indices) +# free_pos_of_boundary[boundary_idx] = j +# end +# +# config = construct_boundary_config(region, doms, free_pos_of_boundary, 0) +# @test config == [true, false] +# +# config = construct_boundary_config(region, doms, free_pos_of_boundary, 1) +# @test config == [true, true] +# end @testset "extract_inner_configs" begin T0 = zero(Tropical{Float64}) T1 = one(Tropical{Float64}) + # extract_inner_configs now returns Vector{UInt64}, not Vector{Bool} contracted_1d = reshape([T1, T0], 2) configs = extract_inner_configs(contracted_1d, 1) @test length(configs) == 1 - @test configs[1] == [false] + @test configs[1] == UInt64(0) # 0 in binary = [false] contracted_2d = reshape([T1, T0, T1, T1], 2, 2) configs = extract_inner_configs(contracted_2d, 2) @test length(configs) == 3 - @test [false, false] in configs - @test [false, true] in configs - @test [true, true] in configs - @test !([true, false] in configs) + # Julia uses column-major order, so: + # Linear index 1 -> [1,1] -> T1 -> config 0 (0b00 = [false, false]) + # Linear index 2 -> [2,1] -> T0 -> config 1 (0b01 = [false, true]) - NOT satisfied, so not in configs + # Linear index 3 -> [1,2] -> T1 -> config 2 (0b10 = [true, false]) + # Linear index 4 -> [2,2] -> T1 -> config 3 (0b11 = [true, true]) + @test UInt64(0) in configs # [false, false] + @test !(UInt64(1) in configs) # [false, true] - NOT satisfied + @test UInt64(2) in configs # [true, false] + @test UInt64(3) in configs # [true, true] contracted_3d = reshape([T1, T1, T1, T1, T1, T1, T1, T1], 2, 2, 2) configs = extract_inner_configs(contracted_3d, 3) @test length(configs) == 8 end -@testset "construct_inner_config" begin - region = Region(1, - [1], - [5, 6, 7], - [3]) - - doms = fill(DM_BOTH, 10) - inner_var_ids = Int[5, 6, 7] - # Build inner_posmap: maps var_id to position in free inner vars - inner_posmap = zeros(Int, maximum(inner_var_ids)) - for (j, var_id) in enumerate(inner_var_ids) - inner_posmap[var_id] = j - end - free_inner_configs = [[false, false, false], [true, false, true]] - - inner_configs = construct_inner_config(region, doms, inner_posmap, free_inner_configs) - @test length(inner_configs) == 2 - @test inner_configs[1] == [false, false, false] - @test inner_configs[2] == [true, false, true] - - doms[5] = DM_1 - inner_var_ids = Int[6, 7] - # Rebuild inner_posmap with new inner_var_ids - inner_posmap = zeros(Int, maximum(inner_var_ids)) - for (j, var_id) in enumerate(inner_var_ids) - inner_posmap[var_id] = j - end - free_inner_configs = [[false, false], [false, true]] - - inner_configs = construct_inner_config(region, doms, inner_posmap, free_inner_configs) - @test length(inner_configs) == 2 - @test inner_configs[1] == [true, false, false] - @test inner_configs[2] == [true, false, true] - - doms[6] = DM_0 - inner_var_ids = Int[7] - # Rebuild inner_posmap with new inner_var_ids - inner_posmap = zeros(Int, maximum(inner_var_ids)) - for (j, var_id) in enumerate(inner_var_ids) - inner_posmap[var_id] = j - end - free_inner_configs = [[false], [true]] - - inner_configs = construct_inner_config(region, doms, inner_posmap, free_inner_configs) - @test length(inner_configs) == 2 - @test inner_configs[1] == [true, false, false] - @test inner_configs[2] == [true, false, true] -end +# construct_inner_config is no longer in the codebase - it appears to have been refactored. +# @testset "construct_inner_config" begin +# region = Region(1, +# [1], +# [5, 6, 7], +# [3]) +# +# doms = fill(DM_BOTH, 10) +# inner_var_ids = Int[5, 6, 7] +# # Build inner_posmap: maps var_id to position in free inner vars +# inner_posmap = zeros(Int, maximum(inner_var_ids)) +# for (j, var_id) in enumerate(inner_var_ids) +# inner_posmap[var_id] = j +# end +# free_inner_configs = [[false, false, false], [true, false, true]] +# +# inner_configs = construct_inner_config(region, doms, inner_posmap, free_inner_configs) +# @test length(inner_configs) == 2 +# @test inner_configs[1] == [false, false, false] +# @test inner_configs[2] == [true, false, true] +# +# doms[5] = DM_1 +# inner_var_ids = Int[6, 7] +# # Rebuild inner_posmap with new inner_var_ids +# inner_posmap = zeros(Int, maximum(inner_var_ids)) +# for (j, var_id) in enumerate(inner_var_ids) +# inner_posmap[var_id] = j +# end +# free_inner_configs = [[false, false], [false, true]] +# +# inner_configs = construct_inner_config(region, doms, inner_posmap, free_inner_configs) +# @test length(inner_configs) == 2 +# @test inner_configs[1] == [true, false, false] +# @test inner_configs[2] == [true, false, true] +# +# doms[6] = DM_0 +# inner_var_ids = Int[7] +# # Rebuild inner_posmap with new inner_var_ids +# inner_posmap = zeros(Int, maximum(inner_var_ids)) +# for (j, var_id) in enumerate(inner_var_ids) +# inner_posmap[var_id] = j +# end +# free_inner_configs = [[false], [true]] +# +# inner_configs = construct_inner_config(region, doms, inner_posmap, free_inner_configs) +# @test length(inner_configs) == 2 +# @test inner_configs[1] == [true, false, false] +# @test inner_configs[2] == [true, false, true] +# end @testset "combine_configs" begin - boundary_config = [true, false] - inner_configs = [[false, true], [true, false]] - - full_configs = combine_configs(boundary_config, inner_configs) + # combine_configs now takes (boundary_bits::UInt64, inner_configs::Vector{UInt64}, n_boundary::Int) + # and returns Vector{UInt64} + # Note: bit representation uses LSB first, so 0b01 = [true, false], 0b10 = [false, true] + boundary_bits = UInt64(0b01) # [true, false] = 1 in binary (bit0=1, bit1=0) + inner_configs = [UInt64(0b01), UInt64(0b10)] # [true, false] and [false, true] + n_boundary = 2 + + full_configs = combine_configs(boundary_bits, inner_configs, n_boundary) @test length(full_configs) == 2 - @test full_configs[1] == [true, false, false, true] - @test full_configs[2] == [true, false, true, false] + # boundary=01 (1), inner=01 (1) -> full = 01 | (01 << 2) = 1 | 4 = 5 (0b101 = [true, false, true, false]) + # boundary=01 (1), inner=10 (2) -> full = 01 | (10 << 2) = 1 | 8 = 9 (0b1001 = [true, false, false, true]) + @test UInt64(5) in full_configs # [true, false, true, false] + @test UInt64(9) in full_configs # [true, false, false, true] - boundary_config = [true] - inner_configs = [[false]] + boundary_bits = UInt64(0b1) # [true] + inner_configs = [UInt64(0b0)] # [false] + n_boundary = 1 - full_configs = combine_configs(boundary_config, inner_configs) + full_configs = combine_configs(boundary_bits, inner_configs, n_boundary) @test length(full_configs) == 1 - @test full_configs[1] == [true, false] + @test full_configs[1] == UInt64(1) # [true, false] - boundary_config = Bool[] - inner_configs = [[true, true], [false, false]] + boundary_bits = UInt64(0b0) # [] + inner_configs = [UInt64(0b11), UInt64(0b00)] # [true, true] and [false, false] + n_boundary = 0 - full_configs = combine_configs(boundary_config, inner_configs) + full_configs = combine_configs(boundary_bits, inner_configs, n_boundary) @test length(full_configs) == 2 - @test full_configs[1] == [true, true] - @test full_configs[2] == [false, false] + @test UInt64(3) in full_configs # [true, true] + @test UInt64(0) in full_configs # [false, false] end diff --git a/test/greedymerge.jl b/test/greedymerge.jl new file mode 100644 index 0000000..7d3f853 --- /dev/null +++ b/test/greedymerge.jl @@ -0,0 +1,225 @@ +using Test +using BooleanInference +using BooleanInference: TNProblem, NumUnfixedVars, setup_problem +using OptimalBranchingCore +using OptimalBranchingCore: Clause, greedymerge, bit_clauses, BranchingTable +using TropicalNumbers: Tropical +using DataStructures: PriorityQueue + +# Helper function to create a simple test problem +function create_test_problem() + dummy_tensors_to_vars = [[1, 2], [2, 3]] + dummy_tensor_data = [ + fill(Tropical(0.0), 4), + fill(Tropical(0.0), 4) + ] + static = BooleanInference.setup_problem(3, dummy_tensors_to_vars, dummy_tensor_data) + return TNProblem(static) +end + + +@testset "size_reduction" begin + problem = create_test_problem() + variables = [1, 2] + measure = NumUnfixedVars() + + # Test with a valid clause + clause = Clause(0b11, 0b10) # variable 1 = 1 + reduction = OptimalBranchingCore.size_reduction(problem, measure, clause, variables) + @test isfinite(reduction) + + # Test with an invalid clause that causes UNSAT + clause_unsat = Clause(0b11, 0b00) # This might cause UNSAT depending on problem + reduction_unsat = OptimalBranchingCore.size_reduction(problem, measure, clause_unsat, variables) + # Either finite or -Inf depending on whether it causes UNSAT + @test isfinite(reduction_unsat) || reduction_unsat == -Inf +end + +@testset "clause_key" begin + # This is an internal function, but we can test it indirectly + clause1 = Clause(0b11, 0b10) + clause2 = Clause(0b11, 0b10) + clause3 = Clause(0b11, 0b01) + + # Keys should be equal for identical clauses + key1 = (getfield(clause1, 1), getfield(clause1, 2)) + key2 = (getfield(clause2, 1), getfield(clause2, 2)) + @test key1 == key2 + + # Keys should differ for different clauses + key3 = (getfield(clause3, 1), getfield(clause3, 2)) + @test key1 != key3 +end + +@testset "cached_size_reduction!" begin + problem = create_test_problem() + variables = [1, 2] + measure = NumUnfixedVars() + clause = Clause(0b11, 0b10) + + cache = Dict{Tuple{UInt64, UInt64}, Float64}() + + # First call should compute and cache + reduction1 = OptimalBranchingCore.size_reduction(problem, measure, clause, variables) + + # Second call with same clause should use cache (if implemented) + reduction2 = OptimalBranchingCore.size_reduction(problem, measure, clause, variables) + @test reduction1 == reduction2 +end + +@testset "greedymerge - drop invalid rows" begin + problem = create_test_problem() + variables = [1, 2] + measure = NumUnfixedVars() + + # Test that clauses with invalid reductions are filtered out + # This tests drop_invalid_rows! indirectly + clauses = [ + [Clause(0b11, 0b10)], # Valid + [Clause(0b11, 0b00)] # Might be invalid + ] + result = OptimalBranchingCore.greedymerge(clauses, problem, variables, measure) + @test !isnothing(result) + result_clauses = OptimalBranchingCore.get_clauses(result) + # Result should only contain valid clauses + @test length(result_clauses) <= length(clauses) +end + +@testset "greedymerge - deduplicate singletons" begin + problem = create_test_problem() + variables = [1, 2] + measure = NumUnfixedVars() + + # Test deduplication of duplicate singleton clauses + clause1 = Clause(0b11, 0b10) + clauses = [ + [clause1], + [clause1], # Duplicate + [Clause(0b11, 0b01)] + ] + result = OptimalBranchingCore.greedymerge(clauses, problem, variables, measure) + @test !isnothing(result) + result_clauses = OptimalBranchingCore.get_clauses(result) + # Should deduplicate, so result clauses should be <= input + @test length(result_clauses) <= length(clauses) +end + +@testset "greedymerge - basic functionality" begin + problem = create_test_problem() + variables = [1, 2] + measure = NumUnfixedVars() + + # Create a simple branching table + # Use a simple case with valid clauses + clauses = [[Clause(0b11, 0b10)], [Clause(0b11, 0b01)]] + + result = OptimalBranchingCore.greedymerge(clauses, problem, variables, measure) + + @test !isnothing(result) + @test isfinite(result.γ) || result.γ == Inf + result_clauses = OptimalBranchingCore.get_clauses(result) + @test !isempty(result_clauses) || isempty(clauses) +end + +@testset "greedymerge - empty clauses" begin + problem = create_test_problem() + variables = [1, 2] + measure = NumUnfixedVars() + + # Test with empty clause set + empty_clauses = Vector{Vector{Clause{UInt64}}}() + result = OptimalBranchingCore.greedymerge(empty_clauses, problem, variables, measure) + @test result.γ == Inf + result_clauses = OptimalBranchingCore.get_clauses(result) + @test isempty(result_clauses) +end + +@testset "greedymerge - invalid clauses" begin + problem = create_test_problem() + variables = [1, 2] + measure = NumUnfixedVars() + + # Test with clauses that might cause UNSAT or invalid reductions + invalid_clauses = [[Clause(0b11, 0b00)], [Clause(0b11, 0b11)]] + result = OptimalBranchingCore.greedymerge(invalid_clauses, problem, variables, measure) + # Result should handle invalid clauses gracefully + @test !isnothing(result) +end + +@testset "greedymerge - with branching table" begin + problem = create_test_problem() + variables = [1, 2] + measure = NumUnfixedVars() + + # Create a branching table and test the full flow + # This tests the integration with optimal_branching_rule + table = BranchingTable(2, [[0x0, 0x1], [0x2, 0x3]]) + candidates = OptimalBranchingCore.bit_clauses(table) + + result = OptimalBranchingCore.greedymerge(candidates, problem, variables, measure) + @test !isnothing(result) + @test isfinite(result.γ) || result.γ == Inf + result_clauses = OptimalBranchingCore.get_clauses(result) + @test !isempty(result_clauses) || isempty(candidates) +end + +@testset "greedymerge - multiple clauses per row" begin + problem = create_test_problem() + variables = [1, 2, 3] + measure = NumUnfixedVars() + + # Test with multiple clauses in a row (should select best representative) + # This tests select_representatives! + clauses = [ + [Clause(0b111, 0b100), Clause(0b111, 0b010)], # Multiple options + [Clause(0b111, 0b001)] + ] + + result = OptimalBranchingCore.greedymerge(clauses, problem, variables, measure) + @test !isnothing(result) + @test isfinite(result.γ) || result.γ == Inf + result_clauses = OptimalBranchingCore.get_clauses(result) + @test !isempty(result_clauses) || isempty(clauses) +end + +@testset "greedymerge - merge queue processing" begin + problem = create_test_problem() + variables = [1, 2] + measure = NumUnfixedVars() + + # Test that merging queue is processed correctly + # This tests enqueue_beneficial_merges!, process_merge_queue!, etc. + clauses = [ + [Clause(0b11, 0b10)], + [Clause(0b11, 0b01)], + [Clause(0b11, 0b00)] + ] + + result = OptimalBranchingCore.greedymerge(clauses, problem, variables, measure) + @test !isnothing(result) + @test isfinite(result.γ) || result.γ == Inf + result_clauses = OptimalBranchingCore.get_clauses(result) + @test length(result_clauses) <= length(clauses) +end + +@testset "greedymerge - merge behavior" begin + problem = create_test_problem() + variables = [1, 2] + measure = NumUnfixedVars() + + # Test that merging actually happens when beneficial + # Create clauses that might benefit from merging + clauses = [ + [Clause(0b11, 0b10)], + [Clause(0b11, 0b01)] + ] + + result = OptimalBranchingCore.greedymerge(clauses, problem, variables, measure) + @test !isnothing(result) + # Result should have valid structure + @test isfinite(result.γ) || result.γ == Inf + result_clauses = OptimalBranchingCore.get_clauses(result) + @test length(result_clauses) <= length(clauses) +end + + diff --git a/test/interface.jl b/test/interface.jl index 1ba4f79..c81979d 100644 --- a/test/interface.jl +++ b/test/interface.jl @@ -87,8 +87,8 @@ end # Solve and check stats are recorded result, depth, stats = BooleanInference.solve(tn_problem, - BranchingStrategy(table_solver=TNContractionSolver(), - selector=LeastOccurrenceSelector(1, 2), + BranchingStrategy(table_solver=TNContractionSolver(1, 2), + selector=MostOccurrenceSelector(), measure=NumUnfixedVars()), NoReducer()) @@ -100,7 +100,7 @@ end # Print stats for debugging println("\nBranching Statistics:") - print_branching_stats(tn_problem) + print_stats_summary(stats) # Test reset functionality reset_branching_stats!(tn_problem) diff --git a/test/runtests.jl b/test/runtests.jl index ca1ea47..8991a84 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -29,6 +29,10 @@ end include("selector.jl") end -@testset "tablesolver.jl" begin - include("tablesolver.jl") -end \ No newline at end of file +@testset "stats.jl" begin + include("stats.jl") +end + +@testset "greedymerge.jl" begin + include("greedymerge.jl") +end diff --git a/test/selector.jl b/test/selector.jl index 8acd6bf..b23c9f1 100644 --- a/test/selector.jl +++ b/test/selector.jl @@ -1,92 +1,75 @@ using Test using BooleanInference -using BooleanInference: LeastOccurrenceSelector -using ProblemReductions: Factoring, reduceto, CircuitSAT -using GenericTensorNetworks -using BooleanInference: setup_from_tensor_network, TNProblem, DynamicWorkspace, NumUnfixedVars -using BooleanInference: get_cached_region, clear_region_cache!, clear_all_region_caches!, setup_from_cnf, k_neighboring -using OptimalBranchingCore: select_variables -using BooleanInference.GenericTensorNetworks: ∧, ∨, ¬ -using BooleanInference.OptimalBranchingCore: select_variables, apply_branch, Clause +using BooleanInference: TNProblem, MostOccurrenceSelector, MinGammaSelector, TNContractionSolver, NumUnfixedVars, setup_problem +using OptimalBranchingCore +using TropicalNumbers -function generate_example_problem(n::Int = 12) - fproblem = Factoring(n, n, 6) - res = reduceto(CircuitSAT, fproblem) - problem = CircuitSAT(res.circuit.circuit; use_constraints=true) - return problem +const dummy_tensors_to_vars = [[1, 2], [2, 3], [2]] +const dummy_tensor_data = [ + fill(Tropical(0.0), 4), + fill(Tropical(0.0), 4), + fill(Tropical(0.0), 2) +] + +function build_dummy_problem() + static = BooleanInference.setup_problem(3, dummy_tensors_to_vars, dummy_tensor_data) + return TNProblem(static) end -@testset "LeastOccurrenceSelector" begin - selector = LeastOccurrenceSelector(1) - @test selector.k == 1 - # Default max_tensors is 2 - @test selector.max_tensors == 2 - selector = LeastOccurrenceSelector(2, 20) - @test selector.k == 2 - @test selector.max_tensors == 20 +@testset "MostOccurrenceSelector" begin + problem = build_dummy_problem() + selector = MostOccurrenceSelector() + chosen = OptimalBranchingCore.select_variables(problem, NumUnfixedVars(), selector) + + # Variable 2 appears in all 3 tensors, so it should be selected + # Variable 1 appears in tensor 1, variable 3 appears in tensor 2 + @test chosen == 2 end -@testset "select_variables" begin - problem = generate_example_problem() - tn = GenericTensorNetwork(problem) - tn_static = setup_from_tensor_network(tn) - tn_problem = TNProblem(tn_static) - selector = LeastOccurrenceSelector(1) - select_variables(tn_problem, NumUnfixedVars(), selector) +# Mock implementations for MinGammaSelector testing +struct MockTableSolver <: OptimalBranchingCore.AbstractTableSolver end +struct MockSetCoverSolver <: OptimalBranchingCore.AbstractSetCoverSolver end - region = get_cached_region(tn_problem) - first_var = region.id - should_involved_tensors = tn_static.v2t[first_var] - @test length(region.tensors) == length(should_involved_tensors) - @test all(t in region.tensors for t in should_involved_tensors) - @test length(region.inner_vars) == 1 - @test region.inner_vars[1] == first_var - # After propagation, some boundary variables may be fixed - # Check that all boundary variables in the region are unfixed and neighboring the focus var - for v in region.boundary_vars - @test !BooleanInference.is_fixed(tn_problem.doms[v]) +# Gamma values for testing: var 2 has lower gamma, var 3 has higher gamma +const gamma_values = Dict(2 => 0.25, 3 => 0.9) + +function OptimalBranchingCore.branching_table( + ::TNProblem, + ::MockTableSolver, + var::Int +)::Tuple{OptimalBranchingCore.BranchingTable, Vector{Int}} + # Return empty table for var 1 (UNSAT), valid tables for others + if var == 1 + return OptimalBranchingCore.BranchingTable(0, Vector{UInt64}[]), Int[] + else + # Return a simple valid branching table + return OptimalBranchingCore.BranchingTable(1, [[0x0, 0x1]]), [var] end - # Check that focus variable is properly in inner vars - @test first_var in region.inner_vars end -@testset "clear_region_cache!" begin - problem = generate_example_problem() - tn = GenericTensorNetwork(problem) - tn_static = setup_from_tensor_network(tn) - tn_problem = TNProblem(tn_static) - @test get_cached_region(tn_problem) == nothing - select_variables(tn_problem, NumUnfixedVars(), LeastOccurrenceSelector(1)) - @test get_cached_region(tn_problem) != nothing - - problem2 = generate_example_problem() - tn2 = GenericTensorNetwork(problem2) - tn_static2 = setup_from_tensor_network(tn2) - tn_problem2 = TNProblem(tn_static2) - @test get_cached_region(tn_problem2) == nothing - select_variables(tn_problem2, NumUnfixedVars(), LeastOccurrenceSelector(1)) - @test get_cached_region(tn_problem2) != nothing - - clear_region_cache!(tn_problem) - @test get_cached_region(tn_problem) == nothing - clear_all_region_caches!() - @test get_cached_region(tn_problem) == nothing - @test get_cached_region(tn_problem2) == nothing +function OptimalBranchingCore.optimal_branching_rule( + ::OptimalBranchingCore.BranchingTable, + variables::Vector{Int}, + ::TNProblem, + ::OptimalBranchingCore.AbstractMeasure, + ::MockSetCoverSolver +) + var = first(variables) + gamma = get(gamma_values, var, Inf) + # Return a result with the gamma value + # The actual structure depends on OptimalBranchingCore, but we need γ field + return (γ=gamma, clauses=[[var]]) end -@testset "KNeighborSelector" begin - @bools a b c d e f g - cnf = ∧(∨(a, b, ¬d, ¬e), ∨(¬a, d, e, ¬f), ∨(f, g), ∨(¬b, c)) - problem = setup_from_cnf(cnf) - # After initial propagation, some variables may be fixed - # Use a larger max_tensors to include all tensors - vars = select_variables(problem, NumUnfixedVars(), LeastOccurrenceSelector(4, 10)) - # After propagation, only unfixed variables are returned - unfixed_vars = BooleanInference.get_unfixed_vars(problem.doms) - @test vars == unfixed_vars +OptimalBranchingCore.get_clauses(result) = result.clauses - region = get_cached_region(problem) - # Check that the region contains tensors related to unfixed variables - @test length(region.tensors) > 0 - @test region.boundary_vars == [] -end \ No newline at end of file +@testset "MinGammaSelector" begin + problem = build_dummy_problem() + selector = MinGammaSelector(MockTableSolver(), MockSetCoverSolver()) + chosen = OptimalBranchingCore.select_variables(problem, NumUnfixedVars(), selector) + + # Variable 2 has gamma 0.25, variable 3 has gamma 0.9 + # Variable 1 should be skipped (empty table) + # So variable 2 should be selected + @test chosen == 2 +end diff --git a/test/stats.jl b/test/stats.jl new file mode 100644 index 0000000..409b71b --- /dev/null +++ b/test/stats.jl @@ -0,0 +1,92 @@ +using Test +using BooleanInference +using BooleanInference: NumUnfixedVars + +@testset "BranchingStats detailed metrics" begin + stats = BranchingStats(true) + @test BooleanInference.needs_path_tracking(stats) + @test stats.avg_branching_factor == 0.0 + + record_branch!(stats, 3, 2) + record_branch!(stats, 2, 1) + + @test stats.total_branches == 2 + @test stats.total_subproblems == 5 + @test stats.max_depth == 2 + @test stats.avg_branching_factor ≈ 2.5 + + d = stats.detailed + @test d.branching_factors == [3, 2] + @test d.depth_distribution[3] == 1 + @test d.depth_distribution[2] == 1 + + record_propagation!(stats, 0.25) + record_propagation!(stats, 0.5) + @test d.propagation_calls == 2 + @test d.time_propagation ≈ 0.75 + + record_domain_reduction!(stats, 7) + record_early_unsat!(stats) + @test d.domain_reductions == 7 + @test d.early_unsat_detections == 1 + + record_branching_time!(stats, 0.1) + record_contraction_time!(stats, 0.2) + record_filtering_time!(stats, 0.3) + @test d.time_branching ≈ 0.1 + @test d.time_contraction ≈ 0.2 + @test d.time_filtering ≈ 0.3 + + record_cache_hit!(stats) + record_cache_hit!(stats) + record_cache_miss!(stats) + @test d.cache_hits == 2 + @test d.cache_misses == 1 + + record_variable_selection!(stats, 5, 10, 2) + record_variable_selection!(stats, 5, 6, 1) + record_variable_selection!(stats, 7, 4, 0) + @test d.variable_selection_counts[5] == 2 + @test d.variable_selection_counts[7] == 1 + @test d.remaining_vars_at_branch == [10, 6, 4] + @test d.depth_at_selection == [2, 1, 0] + @test d.variable_selection_sequence == [5, 5, 7] + + record_solved_leaf!(stats, 3, [1, 2, 3]) + @test !isempty(d.successful_paths) + @test d.successful_paths[1] == [1, 2, 3] + + record_unsat_leaf!(stats, 4) + @test stats.max_depth == 4 + + record_skipped_subproblem!(stats) + + buf = IOBuffer() + print_stats_summary(stats; io=buf) + summary_output = String(take!(buf)) + @test occursin("Branch decisions: 2", summary_output) + @test occursin("Successful paths", summary_output) + + copied = copy(stats) + @test copied !== stats + @test copied.detailed !== d + @test copied.detailed.successful_paths[1] == d.successful_paths[1] + copied.detailed.successful_paths[1][1] = 99 + @test d.successful_paths[1][1] == 1 + + reset!(stats) + @test stats.total_branches == 0 + @test stats.total_subproblems == 0 + @test stats.max_depth == 0 + @test stats.solved_leaves == 0 + @test stats.avg_branching_factor == 0.0 + @test isempty(d.branching_factors) + @test isempty(d.depth_distribution) + @test d.propagation_calls == 0 + @test d.time_propagation == 0.0 + @test isempty(d.variable_selection_counts) + @test isempty(d.successful_paths) + + plain_stats = BranchingStats() + @test !BooleanInference.needs_path_tracking(plain_stats) +end diff --git a/test/tablesolver.jl b/test/tablesolver.jl deleted file mode 100644 index fb01e12..0000000 --- a/test/tablesolver.jl +++ /dev/null @@ -1,27 +0,0 @@ -using BooleanInference -using Test -using BooleanInference.GenericTensorNetworks -using BooleanInference.GenericTensorNetworks: ∧, ∨, ¬ -using BooleanInference.OptimalBranchingCore: select_variables -using BooleanInference.OptimalBranchingCore.BitBasis -# using BooleanInference:SubBIP - -# @testset "branching_table" begin -# @bools a b c d e f g -# cnf = ∧(∨(a, b, ¬d, ¬e), ∨(¬a, d, e, ¬f), ∨(f, g), ∨(¬b, c)) -# bip,syms = convert_cnf_to_bip(cnf) -# bs = initialize_branching_status(bip) - -# subbip = SubBIP(bip,bs,[1,2,3,4,5]) -# tbl = BooleanInference.branching_table(bip, bs, BooleanInference.TNContractionSolver(), subbip) - -# @test subbip.outside_vs_ind == [2,5] -# test_tag = true -# for vec in tbl.table -# for cl in vec -# test_tag = (readbit(cl,2) == readbit(vec[1],2)) && test_tag -# test_tag = (readbit(cl,5) == readbit(vec[1],5)) && test_tag -# end -# end -# @test test_tag -# end \ No newline at end of file