diff --git a/.gitignore b/.gitignore index 25397db..fbce418 100644 --- a/.gitignore +++ b/.gitignore @@ -1,10 +1,24 @@ +# Julia *.jl.*.cov *.jl.cov *.jl.mem /Manifest.toml -/docs/Manifest.toml -/docs/build/ + +# Benchmark data +benchmarks/data/ +benchmarks/Manifest.toml + +# IDE .vscode/ -backup/ -datas/ -.DS_Store \ No newline at end of file +.idea/ + +# OS +.DS_Store +Thumbs.db + +# Temporary files +*.tmp +*.temp + +statprof/ +OptimalBranching.jl/ \ No newline at end of file diff --git a/Project.toml b/Project.toml index 9cd203f..7eebf4d 100644 --- a/Project.toml +++ b/Project.toml @@ -1,29 +1,24 @@ name = "BooleanInference" uuid = "d12a6243-8379-43cb-9459-098119da1519" -authors = ["nzy1997"] version = "1.0.0-DEV" +authors = ["nzy1997"] [deps] +DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" GenericTensorNetworks = "3521c873-ad32-4bb4-b63d-f4f178f42b49" -HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" -JuMP = "4076af6c-e467-56ae-b986-b466b2749572" -KaHyPar = "2a6221f6-aa48-11e9-3542-2d9e0ef01880" OptimalBranchingCore = "c76e7b22-e1d2-40e8-b0f1-f659837787b8" -Primes = "27ebfcd6-29c5-5fa9-bf4b-fb8fc14df3ae" ProblemReductions = "899c297d-f7d2-4ebf-8815-a35996def416" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" TropicalNumbers = "b3a74e9c-7526-4576-a4eb-79c0d4c32334" [compat] -GenericTensorNetworks = "3" -HiGHS = "1.12.2" -JuMP = "1.23.5" -KaHyPar = "0.3.1" +DataStructures = "0.18.22" +GenericTensorNetworks = "4" OptimalBranchingCore = "0.1" -Primes = "0.5.6" +ProblemReductions = "0.3.5" SparseArrays = "1.11.0" TropicalNumbers = "0.6.2" -julia = "1.11" +julia = "1" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/README.md b/README.md index 823d3cb..a0498a5 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,114 @@ -# BooleanInference +# BooleanInference.jl [![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://nzy1997.github.io/BooleanInference.jl/stable/) [![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://nzy1997.github.io/BooleanInference.jl/dev/) [![Build Status](https://github.com/nzy1997/BooleanInference.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/nzy1997/BooleanInference.jl/actions/workflows/CI.yml?query=branch%3Amain) [![Coverage](https://codecov.io/gh/nzy1997/BooleanInference.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/nzy1997/BooleanInference.jl) + +A high-performance Julia package for solving Boolean satisfiability problems using tensor network contraction and optimal branching strategies. + +## Features + +- **Tensor Network Representation**: Efficiently represents Boolean satisfiability problems as tensor networks +- **Optimal Branching**: Uses advanced branching strategies to minimize search space +- **Multiple Problem Types**: Supports CNF, circuit, and factoring problems +- **High Performance**: Optimized for speed with efficient propagation and contraction algorithms +- **Flexible Interface**: Easy-to-use API for various constraint satisfaction problems + +## Installation + +```julia +using Pkg +Pkg.add("BooleanInference") +``` + +## Quick Start + +### Solving SAT Problems + +```julia +using BooleanInference +using GenericTensorNetworks: ∧, ∨, ¬ + +# Define a CNF formula +@bools a b c d e f g +cnf = ∧(∨(a, b, ¬d, ¬e), ∨(¬a, d, e, ¬f), ∨(f, g), ∨(¬b, c), ∨(¬a)) + +# Solve and get assignments +sat = Satisfiability(cnf; use_constraints=true) +satisfiable, assignments, depth = solve_sat_with_assignments(sat) + +println("Satisfiable: ", satisfiable) +println("Assignments: ", assignments) +``` + +### Solving Factoring Problems + +```julia +# Factor a semiprime number +a, b = solve_factoring(5, 5, 31*29) +println("Factors: $a × $b = $(a*b)") +``` + +### Circuit Problems + +```julia +# Solve circuit satisfiability +circuit = @circuit begin + c = x ∧ y +end +push!(circuit.exprs, Assignment([:c], BooleanExpr(true))) + +tnproblem = setup_from_circuit(circuit) +result, depth = solve(tnproblem, BranchingStrategy(), NoReducer()) +``` + +## Core Components + +### Problem Types +- `TNProblem`: Main problem representation +- `TNStatic`: Static problem structure +- `DomainMask`: Variable domain representation + +### Solvers +- `TNContractionSolver`: Tensor network contraction-based solver +- `LeastOccurrenceSelector`: Variable selection strategy +- `NumUnfixedVars`: Measurement strategy + +### Key Functions +- `solve()`: Main solving function +- `setup_from_cnf()`: Setup from CNF formulas +- `setup_from_circuit()`: Setup from circuit descriptions +- `solve_factoring()`: Solve integer factoring problems + +## Advanced Usage + +### Custom Branching Strategy + +```julia +using OptimalBranchingCore: BranchingStrategy + +# Configure custom solver +bsconfig = BranchingStrategy( + table_solver=TNContractionSolver(), + selector=LeastOccurrenceSelector(2, 10), + measure=NumUnfixedVars() +) + +# Solve with custom configuration +result, depth = solve(problem, bsconfig, NoReducer()) +``` + +### Benchmarking + +The package includes comprehensive benchmarking tools: + +```julia +using BooleanInferenceBenchmarks + +# Compare different solvers +configs = [(10,10), (12,12), (14,14)] +results = run_solver_comparison(FactoringProblem, configs) +print_solver_comparison_summary(results) +``` + diff --git a/benchmarks/Project.toml b/benchmarks/Project.toml new file mode 100644 index 0000000..7f0e4d5 --- /dev/null +++ b/benchmarks/Project.toml @@ -0,0 +1,27 @@ +name = "BooleanInferenceBenchmarks" +uuid = "3611a03f-6421-4259-a011-a603192356a0" +version = "0.1.0" +authors = ["Xiwei Pan "] + +[deps] +BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +BooleanInference = "d12a6243-8379-43cb-9459-098119da1519" +Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" +Gurobi = "2e9cd046-0924-5485-92f1-d5272153d98b" +JSON3 = "0f8b85d8-7281-11e9-16c2-39a750bddbf1" +JuMP = "4076af6c-e467-56ae-b986-b466b2749572" +Primes = "27ebfcd6-29c5-5fa9-bf4b-fb8fc14df3ae" +ProblemReductions = "899c297d-f7d2-4ebf-8815-a35996def416" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce" +Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" + +[compat] +BenchmarkTools = "1.6.0" +Dates = "1.11.0" +Gurobi = "1.7.6" +JSON3 = "1.14.3" +JuMP = "1.29.1" +ProblemReductions = "0.3.5" +Random = "1.11.0" +Statistics = "1.11.1" diff --git a/benchmarks/README.md b/benchmarks/README.md new file mode 100644 index 0000000..5063fd1 --- /dev/null +++ b/benchmarks/README.md @@ -0,0 +1,106 @@ +# BooleanInference Benchmarks + +This directory contains benchmarking code for BooleanInference.jl, organized as a separate Julia package with a modern **multiple dispatch architecture** to keep benchmark dependencies isolated from the main package. + +## Architecture + +The benchmark system uses Julia's **multiple dispatch** with abstract types to provide a clean, extensible interface for different problem types. + +## Structure + +```text +benchmark/ +├── Project.toml # Benchmark package dependencies +├── src/ +│ ├── BooleanInferenceBenchmarks.jl # Main benchmark module +│ ├── abstract_types.jl # Abstract type definitions and interfaces +│ ├── generic_benchmark.jl # Generic benchmark framework +│ ├── factoring_problem.jl # Factoring problem implementation +│ └── utils.jl # Generic utilities +├── scripts/ +│ ├── run_benchmarks.jl # Standalone benchmark runner +│ └── example_usage.jl # Usage examples +├── data/ # Generated datasets (gitignored) +└── README.md # This file +``` + +## Usage + +### Quick Start + +```bash +# Run example usage +julia --project=benchmark benchmark/scripts/example_usage.jl +``` + +### Programmatic Usage + +```julia +using Pkg; Pkg.activate("benchmark") +using BooleanInferenceBenchmarks + +# Create problem configurations +configs = [FactoringConfig(10, 10), FactoringConfig(12, 12)] + +# Generate datasets using multiple dispatch +generate_datasets(FactoringProblem; configs=configs, per_config=100) + +# Run benchmarks using multiple dispatch +results = benchmark_problem(FactoringProblem; configs=configs, samples_per_config=5) + +# Run complete benchmark suite +full_results = run_full_benchmark(FactoringProblem) +``` + +## Adding New Problem Types + +The multiple dispatch architecture makes adding new problem types extremely simple: + +```julia +# 1. Define problem and config types +struct YourProblem <: AbstractBenchmarkProblem end +struct YourConfig <: AbstractProblemConfig + param1::Int + param2::String +end + +# 2. Implement the 5 required interface methods +function generate_instance(::Type{YourProblem}, config::YourConfig; rng, include_solution=false) + # Generate problem instance +end + +function solve_instance(::Type{YourProblem}, instance) + # Solve the instance +end + +function problem_id(::Type{YourProblem}, config::YourConfig, data) + # Generate unique ID +end + +function default_configs(::Type{YourProblem}) + # Return default configurations +end + +function filename_pattern(::Type{YourProblem}, config::YourConfig) + # Generate filename pattern +end + +# 3. That's it! Use the same generic functions: +generate_datasets(YourProblem) +benchmark_problem(YourProblem) +run_full_benchmark(YourProblem) +``` + +## Key Advantages + +- **DRY Principle**: Write benchmark logic once, use for all problem types +- **Type Safety**: Julia's type system catches errors at compile time +- **Extensibility**: Adding new problems requires minimal code +- **Consistency**: All problem types use the same interface +- **Performance**: Multiple dispatch enables efficient, optimized code + +## Data Management + +- Datasets are generated in `benchmark/data/` +- Add `benchmark/data/` to `.gitignore` to avoid committing large files +- Use JSONL format for datasets (one JSON object per line) diff --git a/benchmarks/examples/run_benchmarks.jl b/benchmarks/examples/run_benchmarks.jl new file mode 100755 index 0000000..33b118d --- /dev/null +++ b/benchmarks/examples/run_benchmarks.jl @@ -0,0 +1,70 @@ +# 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 new file mode 100644 index 0000000..dff693f --- /dev/null +++ b/benchmarks/src/BooleanInferenceBenchmarks.jl @@ -0,0 +1,37 @@ +module BooleanInferenceBenchmarks + +using Random +using JSON3 +using BenchmarkTools +using BooleanInference +using JuMP +using ProblemReductions +using Statistics: mean, median, std, quantile +using Dates: now +using SHA: bytes2hex, sha256 + + +# Core utilities +include("utils.jl") + +# 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 +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 + +end # module \ No newline at end of file diff --git a/benchmarks/src/abstract_types.jl b/benchmarks/src/abstract_types.jl new file mode 100644 index 0000000..4dc94fc --- /dev/null +++ b/benchmarks/src/abstract_types.jl @@ -0,0 +1,98 @@ +# 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 AbstractSolver end + +# Core interface methods that each problem type must implement + +""" + generate_instance(problem_type::Type{<:AbstractBenchmarkProblem}, config::AbstractProblemConfig; rng, include_solution=false) + +Generate a single instance of the given problem type. +""" +function generate_instance end + +""" + solve_instance(problem_type::Type{<:AbstractBenchmarkProblem}, instance, solver::AbstractSolver) + +Solve a single problem instance using the specified solver. +""" +function solve_instance end + +""" + solve_instance(problem_type::Type{<:AbstractBenchmarkProblem}, instance) + +Solve a single problem instance using the default solver. +""" +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 + +""" + solver_name(solver::AbstractSolver) + +Return human-readable name for the solver. +""" +function solver_name end diff --git a/benchmarks/src/factoring/factor_IP.jl b/benchmarks/src/factoring/factor_IP.jl new file mode 100644 index 0000000..14bf4e5 --- /dev/null +++ b/benchmarks/src/factoring/factor_IP.jl @@ -0,0 +1,82 @@ +""" + 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 new file mode 100644 index 0000000..918ae60 --- /dev/null +++ b/benchmarks/src/factoring/factoring_problem.jl @@ -0,0 +1,205 @@ +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/generic_benchmark.jl b/benchmarks/src/generic_benchmark.jl new file mode 100644 index 0000000..c8b0924 --- /dev/null +++ b/benchmarks/src/generic_benchmark.jl @@ -0,0 +1,834 @@ +# 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 new file mode 100644 index 0000000..b534a98 --- /dev/null +++ b/benchmarks/src/utils.jl @@ -0,0 +1,89 @@ +# 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...) + isdir(dir) || mkpath(dir) + 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 + for line in eachline(io) + isempty(strip(line)) && continue + push!(objs, JSON3.read(line)) + 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/src/BooleanInference.jl b/src/BooleanInference.jl index 0c43bae..c8e1133 100644 --- a/src/BooleanInference.jl +++ b/src/BooleanInference.jl @@ -1,59 +1,81 @@ module BooleanInference -using JuMP -using HiGHS using TropicalNumbers using SparseArrays using OptimalBranchingCore -using OptimalBranchingCore: AbstractProblem,select_variables,apply_branch,reduce_problem,_vec2int,optimal_branching_rule,candidate_clauses +using OptimalBranchingCore: AbstractProblem, select_variables, reduce_problem, _vec2int, candidate_clauses using OptimalBranchingCore.BitBasis using GenericTensorNetworks using GenericTensorNetworks.OMEinsum import ProblemReductions -import ProblemReductions: CircuitSAT,Circuit,Factoring,reduceto -using KaHyPar +import ProblemReductions: CircuitSAT, Circuit, Factoring, reduceto, Satisfiability +using DataStructures +using DataStructures: PriorityQueue # using GenericTensorNetworks: ∧, ∨, ¬ -# status -export BranchingStatus, initialize_branching_status -# stride -export tensor2vec,get_tensor_number,slice_tensor, vec2tensor,vec2lluint,lluint2vec -# types -export BooleanInferenceProblem,BooleanResultBranchCount,NumOfVertices,NumOfClauses,NumOfDegrees +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") -# interface -export cnf2bip,cir2bip,sat2bip,solvebip,solve_factoring,solve_sat,solve_factoring_count +# Export core types +export Variable, EdgeRef, BoolTensor, TNStatic, DomainMask, TNProblem +export DM_BOTH, DM_0, DM_1 +export Region, RegionContraction, RegionCacheEntry, RegionCacheState, RegionCache +export DynamicWorkspace -# reducer -export NoReducer,decide_literal +# Export domain mask functions +export is_fixed, has0, has1, init_doms, get_var_value -# selector -export KNeighborSelector,neighboring,k_neighboring,subhg,Smallest2NeighborSelector,KaHyParSelector +# Export problem setup functions +export setup_problem, setup_from_tensor_network, setup_from_cnf, setup_from_circuit, setup_from_sat -# tablesolver -export TNContractionSolver +# Export problem state functions +export is_solved, cache_branch_solution!, reset_last_branch_problem!, has_last_branch_problem, last_branch_problem -# readcnf -export readcnf,solvecnf +# Export solving functions +export solve, solve_sat_problem, solve_sat_with_assignments, solve_factoring -# tropicalsvd -export tropical_svd +# Export measure types +export NumUnfixedVars -# inspired by boolean matrix factorization -export MEBF -export find_formal_concepts +# Export selector types +export LeastOccurrenceSelector, AbstractSelector + +# 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 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 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 -include("status.jl") -include("stride.jl") -include("types.jl") -include("interface.jl") -include("reducer.jl") -include("selector.jl") -include("tablesolver.jl") -include("branch.jl") -include("readcnf.jl") -include("tropicalsvd.jl") -include("bmf_mebf.jl") -include("bmf_factor.jl") end diff --git a/src/bmf_factor.jl b/src/bmf_factor.jl deleted file mode 100644 index a8f1e8d..0000000 --- a/src/bmf_factor.jl +++ /dev/null @@ -1,75 +0,0 @@ -# This is an exact algorithm for Boolean matrix factorization (BMF) based on the MBF algorithm. -# https://doi.org/10.1016/j.jcss.2009.05.002 - -# compute C↓ (from D to C) -function down_closure(I::BitMatrix, D::Vector{Int}) - C = Int[x for x in 1:size(I, 1) if all(I[x, y] == 1 for y in D)] - return C -end - -# compute D↑ (from C to D) -function up_closure(I::BitMatrix, C::Vector{Int}) - D = Int[y for y in 1:size(I, 2) if all(I[x, y] == 1 for x in C)] - return D -end - -function restore_matrix(F::Set{Tuple{Vector{Int}, Vector{Int}}}, m::Int, n::Int) - # Initialize C and D matrices with Bool type - C = zeros(Bool, m, length(F)) # C matrix: m rows, length of F columns - D = zeros(Bool, length(F), n) # D matrix: length of F rows, n columns - - # Iterate over the set F - for (idx, (C_cols, D_rows)) in enumerate(F) - # Set values in C matrix - for col in 1:length(C_cols) - for row in C_cols[col] - C[row, idx] = true # Set corresponding position in C to true - end - end - - # Set values in D matrix - for col in 1:length(D_rows) - for row in D_rows[col] - D[idx, row] = true # Set corresponding position in D to true - end - end - end - return Matrix{TropicalAndOr}(C), Matrix{TropicalAndOr}(D) -end - - -function find_formal_concepts(I::Matrix{TropicalAndOr}) - # Convert the input matrix to a BitMatrix - I = getproperty.(I, :n) - U = I .== 1 - F = Set{Tuple{Vector{Int}, Vector{Int}}}() - while sum(U) > 0 - j_list = collect(1: size(I, 2)) - best_V = 0 - best_C, best_D, best_j = Int[], Int[], nothing - - for j in j_list - if j ∉ best_D - C_temp = down_closure(I, vcat(best_D, Int[j])) - D_temp = up_closure(I, C_temp) - V_temp = compute_coverage_matrix((C_temp, D_temp), U) - if V_temp > best_V - best_V, best_C, best_D, best_j = V_temp, C_temp, D_temp, j - end - end - end - push!(F, (best_C, best_D)) - U[best_C, best_D] .= false - end - - return restore_matrix(F, size(I, 1), size(I, 2)) -end - -function compute_coverage_matrix(concepts::Tuple{Vector{Int}, Vector{Int}}, U::BitMatrix) - covered = falses(size(U)) - - covered[concepts[1], concepts[2]] .= true - - coverage_ratio = sum(covered .& U) - return coverage_ratio -end diff --git a/src/bmf_mebf.jl b/src/bmf_mebf.jl deleted file mode 100644 index 812bf7e..0000000 --- a/src/bmf_mebf.jl +++ /dev/null @@ -1,116 +0,0 @@ -# This is an approximate algorithm for Boolean matrix factorization (BMF) based on the MBF algorithm. -function MEBF(MATₜ::Matrix{TropicalAndOr}; DIM::Int=1000, Thres::Float64=0.5, P::Float64=0.05) - # This function is directly migrated from the R language code provided in the paper. https://doi.org/10.1609/aaai.v34i04.6072 - # Source: https://github.com/clwan/MEBF.git - # Dim: maximum number of patterns (k) - # Thres: between 0 and 1. A smaller t could achieve higher coverage with less number of patterns, while a larger t enables a more sparse decomposition of the input matrix with greater number of patterns. - # P: percentage of uncovered 1s - MAT = getproperty.(MATₜ, :n) - if min(size(MAT)...) < DIM - DIM = min(size(MAT)...) - 1 - end - - m1 = copy(MAT) - SUM = sum(MAT) - MAT_B = Vector{Int}[] - MAT_C = Vector{Int}[] - - for _ in 1:DIM - if sum(m1) <= P * SUM - @info "Uncovered percentage is lower than $P. No more growth possible." - break - - end - - C1 = 0 - B1_use = zeros(Int, size(m1, 1)) - B2_use = zeros(Int, size(m1, 2)) - - COL = vec(sum(m1, dims=1)) # column sum - ROW = vec(sum(m1, dims=2)) # row sum - - if median(filter(x -> x > 0, COL)) > 1 - TEMP = findall(==(minimum(filter(x -> x >= median(filter(x -> x > 0, COL)), COL))), COL) - for j in TEMP - B1 = m1[:, j] - B2 = zeros(Int, size(m1, 2)) - B2[findall(x -> x >= min(Thres * sum(B1) + 1, sum(B1)), vec(sum(m1[B1 .== 1, :], dims=1)))] .= 1 - C2 = (sum(B1) - 1) * (sum(B2) - 1) - sum(m1[B1 .== 1, B2 .== 1] .== 0) - if C2 > C1 - C1 = C2 - B1_use = B1 - B2_use = B2 - end - end - end - - if median(filter(x -> x > 0, ROW)) > 1 - TEMP = findall(==(minimum(filter(x -> x >= median(filter(x -> x > 0, ROW)), ROW))), ROW) - for j in TEMP - B2 = m1[j, :] - B1 = zeros(Int, size(m1, 1)) - B1[findall(x -> x >= min(Thres * sum(B2) + 1, sum(B2)), vec(sum(m1[:, B2 .== 1], dims=2)))] .= 1 - C2 = (sum(B1) - 1) * (sum(B2) - 1) - sum(m1[B1 .== 1, B2 .== 1] .== 0) - if C2 > C1 - C1 = C2 - B1_use = B1 - B2_use = B2 - end - end - end - - if C1 == 0 - ROW_ORDER = sortperm(ROW, rev=true) - COL_ORDER = sortperm(COL, rev=true) - - B1_1 = zeros(Int, size(m1, 1)) - if count(x -> x == 2, sum(m1[:, COL_ORDER[1:2]], dims=2)) > 1 - B1_1[findall(x -> x == 2, sum(m1[:, COL_ORDER[1:2]], dims=2))] .= 1 - B1_2 = zeros(Int, size(m1, 2)) - B1_2[findall(x -> x >= min(Thres * sum(m1[B1_1 .== 1, COL_ORDER[1]]) + 1, sum(m1[B1_1 .== 1, COL_ORDER[1]])), vec(sum(m1[B1_1 .== 1, :], dims=1)))] .= 1 - C1 = (sum(B1_1) - 1) * (sum(B1_2) - 1) - sum(m1[B1_1 .== 1, B1_2 .== 1] .== 0) - else - C1 = -Inf - end - - B2_2 = zeros(Int, size(m1, 2)) - if count(x -> x == 2, sum(m1[ROW_ORDER[1:2], :], dims=1)) > 1 - B2_2[findall(x -> x == 2, vec(sum(m1[ROW_ORDER[1:2], :], dims=1)))] .= 1 - B2_1 = zeros(Int, size(m1, 1)) - B2_1[findall(x -> x >= min(Thres * sum(m1[ROW_ORDER[1], B2_2 .== 1]) + 1, sum(m1[ROW_ORDER[1], B2_2 .== 1])), sum(m1[:, B2_2 .== 1], dims=2))] .= 1 - C2 = (sum(B2_1) - 1) * (sum(B2_2) - 1) - sum(m1[B2_1 .== 1, B2_2 .== 1] .== 0) - else - C2 = -Inf - end - - if C1 == -Inf && C2 == -Inf - break - elseif C1 > C2 - B1_use = B1_1 - B2_use = B1_2 - m1[B1_1 .== 1, B1_2 .== 1] .= 0 - else - B1_use = B2_1 - B2_use = B2_2 - m1[B2_1 .== 1, B2_2 .== 1] .= 0 - end - end - MAT_B = push!(MAT_B, B1_use) - MAT_C = push!(MAT_C, B2_use) - m1[B1_use .== 1, B2_use .== 1] .= 0 - end - @info "MEBF finished. Number of patterns: $(length(MAT_B))" - return (Matrix{TropicalAndOr}(BitMatrix(hcat(MAT_B...))), Matrix{TropicalAndOr}(BitMatrix(hcat(MAT_C...)'))) -end - -function median(arr) - sorted_arr = sort(arr) - n = length(sorted_arr) - if n == 0 - error("Cannot compute median of an empty array") - elseif isodd(n) - return sorted_arr[(n ÷ 2) + 1] - else - return (sorted_arr[n ÷ 2] + sorted_arr[(n ÷ 2) + 1]) / 2 - end -end diff --git a/src/branch.jl b/src/branch.jl index 649cd86..f9246f6 100644 --- a/src/branch.jl +++ b/src/branch.jl @@ -1,110 +1,169 @@ -function OptimalBranchingCore.apply_branch(p::BooleanInferenceProblem, bs::AbstractBranchingStatus,clause::Clause{INT}, vertices::Vector{T}) where {INT<:Integer, T<:Integer} - bs_new, aedges = decide_literal(bs,p, vertices, clause) - return deduction_reduce(p, bs_new, aedges) +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.branch_and_reduce(problem::BooleanInferenceProblem, bs::AbstractBranchingStatus,config::BranchingStrategy, reducer::AbstractReducer) - stopped, res,count_num = check_stopped(bs) - stopped && return res, bs,count_num +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 - # branch the problem - subbip = select_variables(problem, bs, config.measure, config.selector) # select a subset of variables - tbl = branching_table(problem,bs, config.table_solver, subbip) # compute the BranchingTable - iszero(tbl.bit_length) && return false,bs,1 - result = optimal_branching_rule(tbl, subbip.vs, bs,problem, config.measure, config.set_cover_solver) # compute the optimal branching rule - for branch in OptimalBranchingCore.get_clauses(result) - res, bs_new ,count_num1= branch_and_reduce(problem, apply_branch(problem,bs, branch, subbip.vs), config, reducer) - count_num += count_num1 - if res - return res, bs_new,count_num +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 - return false, bs,count_num -end - -function check_stopped(bs::AbstractBranchingStatus) - # global BRANCHNUMBER - # if there is no clause, then the problem is solved. - if all(bs.undecided_literals .== -1) - # BRANCHNUMBER += 1 - return true,true,1 + + @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 - if any(bs.undecided_literals .== 0) - # BRANCHNUMBER += 1 - return true,false,1 + + # 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 - return false,false,0 -end -function OptimalBranchingCore.optimal_branching_rule(table::BranchingTable, variables::Vector, bs::AbstractBranchingStatus,p::BooleanInferenceProblem, m::AbstractMeasure, solver::AbstractSetCoverSolver) - candidates = collect(candidate_clauses(table)) - size_reductions = [measure(bs, m) - measure((apply_branch(p,bs, candidate, variables)), m) for candidate in candidates] - return minimize_γ(table, candidates, size_reductions, solver; γ0=2.0) + # 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 mybranch_and_reduce(problem::BooleanInferenceProblem, bs::AbstractBranchingStatus,config::BranchingStrategy, reducer::AbstractReducer) - stopped, res = check_stopped(bs) - stopped && return res, bs +function OptimalBranchingCore.reduce_problem(::Type{T}, problem::TNProblem, ::OptimalBranchingCore.NoReducer) where T + return (problem, one(T)) +end - v = findfirst(x -> readbit(bs.decided_mask,x) == 0, 1:problem.literal_num) +# function reduce_problem(::Type{T}, problem::TNProblem, reducer::UnitPropagationReducer) where T +# propagated = propagate(problem.static, problem.doms) - res, bs_new = branch_and_reduce(problem, apply_branch(problem,bs, Clause(0b1, 0b1), [v]), config, reducer) - if res - return res, bs_new - end +# doms = problem.doms +# changed = false +# n_unfixed::Int = 0 - res, bs_new = branch_and_reduce(problem, apply_branch(problem,bs, Clause(0b1, 0b0), [v]), config, reducer) - if res - return res, bs_new - end +# @inbounds for i in eachindex(doms) +# dm_new = propagated[i] +# bits = dm_new.bits - return false, bs -end +# if bits == 0x00 +# clear_region_cache!(problem) +# return (TNProblem(problem.static, propagated, 0, problem.ws), zero(T)) +# end -function OptimalBranchingCore.optimal_branching_rule(table::BranchingTable, variables::Vector, bs::AbstractBranchingStatus,p::BooleanInferenceProblem, m::AbstractMeasure, solver::AbstractSetCoverSolver) - candidates = OptimalBranchingCore.bit_clauses(table) - return OptimalBranchingCore.greedymerge(candidates, p,bs, variables, m) -end +# changed |= bits != doms[i].bits +# n_unfixed += is_fixed(dm_new) ? 0 : 1 +# end -# TODO: NOT DRY! -function OptimalBranchingCore.greedymerge(cls::Vector{Vector{Clause{INT}}}, problem::AbstractProblem,bs::AbstractBranchingStatus, variables::Vector, m::AbstractMeasure) where {INT} - active_cls = collect(1:length(cls)) - cls = copy(cls) - merging_pairs = [(i, j) for i in active_cls, j in active_cls if i < j] - n = length(variables) - size_reductions = [OptimalBranchingCore.size_reduction(problem, m,bs, candidate[1], variables) for candidate in cls] - γ = OptimalBranchingCore.complexity_bv(size_reductions) - while !isempty(merging_pairs) - i, j = popfirst!(merging_pairs) - if i in active_cls && j in active_cls - for ii in 1:length(cls[i]), jj in 1:length(cls[j]) - if OptimalBranchingCore.bdistance(cls[i][ii], cls[j][jj]) == 1 - cl12 = OptimalBranchingCore.gather2(n, cls[i][ii], cls[j][jj]) - if cl12.mask == 0 - continue - end - l12 = OptimalBranchingCore.size_reduction(problem, m,bs, cl12, variables) - if γ^(-size_reductions[i]) + γ^(-size_reductions[j]) >= γ^(-l12) + 1e-8 - push!(cls, [cl12]) - k = length(cls) - deleteat!(active_cls, findfirst(==(i), active_cls)) - deleteat!(active_cls, findfirst(==(j), active_cls)) - for ii in active_cls - push!(merging_pairs, (ii, k)) - end - push!(active_cls, k) - push!(size_reductions, l12) - γ = OptimalBranchingCore.complexity_bv(size_reductions[active_cls]) - break - end - end - end - end - end - return [cl[1] for cl in cls[active_cls]] -end +# if !changed +# return (problem, one(T)) +# end -function OptimalBranchingCore.size_reduction(p::AbstractProblem, m::AbstractMeasure,bs::AbstractBranchingStatus, cl::Clause{INT}, variables::Vector) where {INT} - return measure(bs, m) - measure(apply_branch(p,bs, cl, variables), m) -end +# clear_region_cache!(problem) +# return (TNProblem(problem.static, propagated, n_unfixed, problem.ws), one(T)) +# end diff --git a/src/branchtable.jl b/src/branchtable.jl new file mode 100644 index 0000000..2f8fa0a --- /dev/null +++ b/src/branchtable.jl @@ -0,0 +1,223 @@ +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/contraction.jl b/src/contraction.jl new file mode 100644 index 0000000..8535772 --- /dev/null +++ b/src/contraction.jl @@ -0,0 +1,116 @@ +function contract_region(tn::TNStatic, region::Region, doms::Vector{DomainMask}) + # @show region.id + # println(union(region.boundary_vars, region.inner_vars)) + n_tensors = length(region.tensors) + + sliced_tensors = Vector{Vector{Tropical{Float64}}}(undef, n_tensors) + tensor_indices = Vector{Vector{Int}}(undef, n_tensors) + + @inbounds for (i, tensor_id) in enumerate(region.tensors) + # Get original tensor and its variable axes + original_tensor = tn.tensors[tensor_id] + var_axes = original_tensor.var_axes + tensor_data = original_tensor.tensor + + # Slice tensor according to fixed variables + sliced = slicing(tensor_data, doms, var_axes) + sliced_tensors[i] = sliced + + # Determine remaining (unfixed) variable indices + # After slicing, only unfixed variables remain as indices + remaining_vars = Int[] + for var_id in var_axes + if !is_fixed(doms[var_id]) + push!(remaining_vars, var_id) + end + end + tensor_indices[i] = remaining_vars + end + + # Output: only unfixed variables in the region + # (fixed variables have already been sliced out) + output_vars = Int[] + + # Collect unfixed variables + for var_id in region.boundary_vars + if !is_fixed(doms[var_id]) + push!(output_vars, var_id) + end + end + for var_id in region.inner_vars + if !is_fixed(doms[var_id]) + push!(output_vars, var_id) + end + end + + if isempty(output_vars) + # if all variables are fixed, output is scalar + contracted = contract_tensors(sliced_tensors, tensor_indices, Int[]) + @assert length(contracted) == 1 + return contracted, Int[] + else + contracted = contract_tensors(sliced_tensors, tensor_indices, output_vars) + return contracted, output_vars + end +end + +function contract_tensors(tensors::Vector{Vector{T}}, ixs::Vector{Vector{Int}}, iy::Vector{Int}) where T + eincode = EinCode(ixs, iy) + optcode = optimize_code(eincode, uniformsize(eincode, 2), GreedyMethod()) + unwrapped_tensors = [tensor_unwrapping(t) for t in tensors] + result = optcode(unwrapped_tensors...) + return result +end + + +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 + 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) + 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) + + # Iterate over all assignments to free variables + @inbounds for free_idx in 0:(out_len-1) + # Build full index + full_idx = 0 + + # 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 + + +function tensor_unwrapping(vec::Vector{T}) where T + len = length(vec) + @assert len > 0 + k = trailing_zeros(len) + @assert (1 << k) == len "vector length is not power-of-two" + dims = ntuple(_->2, k) + return reshape(vec, dims) +end \ No newline at end of file diff --git a/src/interface.jl b/src/interface.jl index 9228a5e..423c6de 100644 --- a/src/interface.jl +++ b/src/interface.jl @@ -1,45 +1,74 @@ -function cnf2bip(cnf::CNF) - return sat2bip(Satisfiability(cnf)) +function setup_from_cnf(cnf::CNF) + return setup_from_sat(Satisfiability(cnf; use_constraints=true)) end -function sat2bip(sat::ConstraintSatisfactionProblem) - problem = GenericTensorNetwork(sat) - he2v = getixsv(problem.code) - tensors = GenericTensorNetworks.generate_tensors(Tropical(1.0), problem) - vec_tensors = [vec(t) for t in tensors] - new_tensors = [replace(t,Tropical(1.0) => zero(Tropical{Float64})) for t in vec_tensors] - return BooleanInferenceProblem(new_tensors, he2v, length(problem.problem.symbols)), problem.problem.symbols +function setup_from_circuit(cir::Circuit) + return setup_from_sat(CircuitSAT(cir; use_constraints=true)) end -function cir2bip(cir::Circuit) - return sat2bip(CircuitSAT(cir)) +function setup_from_sat(sat::ConstraintSatisfactionProblem) + tn = GenericTensorNetwork(sat) + static = setup_from_tensor_network(tn) + TNProblem(static) end -function solvesat(sat::ConstraintSatisfactionProblem; bsconfig::BranchingStrategy = BranchingStrategy(table_solver = TNContractionSolver(), selector = KNeighborSelector(1,1), measure=NumOfVertices()), reducer=NoReducer()) - p,syms = sat2bip(sat) - return solvebip(p; bsconfig, reducer) +function solve(problem::TNProblem, bsconfig::BranchingStrategy, reducer::AbstractReducer) + reset_branching_stats!(problem) # Reset stats before solving + depth = branch_and_reduce(problem, bsconfig, reducer, Tropical{Float64}; show_progress=false) + res = last_branch_problem(problem) + stats = get_branching_stats(problem) + return (res, depth, stats) end -function solve_factoring(n::Int, m::Int, N::Int; bsconfig::BranchingStrategy = BranchingStrategy(table_solver = TNContractionSolver(), selector =KNeighborSelector(1,1), measure=NumOfVertices()), reducer=NoReducer()) - # global BRANCHNUMBER = 0 - fproblem = Factoring(m, n, N) - res = reduceto(CircuitSAT,fproblem) - ans,vals = solvesat(res.circuit; bsconfig, reducer) - a, b = ProblemReductions.read_solution(fproblem, [vals[res.p]...,vals[res.q]...]) - @show ans - # @show BRANCHNUMBER - return a,b +function solve_sat_problem( + sat::ConstraintSatisfactionProblem; + bsconfig::BranchingStrategy=BranchingStrategy( + table_solver=TNContractionSolver(), + selector=LeastOccurrenceSelector(1, 2), + measure=NumUnfixedVars() + ), + reducer::AbstractReducer=NoReducer() +) + tn_problem = setup_from_sat(sat) + result, depth, stats = solve(tn_problem, bsconfig, reducer) + satisfiable = !isnothing(result) + return satisfiable, result, depth, stats end -function solve_sat(sat::ConstraintSatisfactionProblem) - res,vals = solvesat(sat) - return res, Dict(zip(sat.symbols,vals)) +function solve_sat_with_assignments( + sat::ConstraintSatisfactionProblem; + bsconfig::BranchingStrategy=BranchingStrategy( + table_solver=TNContractionSolver(), + selector=LeastOccurrenceSelector(1, 2), + measure=NumUnfixedVars() + ), + reducer::AbstractReducer=NoReducer() +) + satisfiable, result, depth, stats = solve_sat_problem(sat; bsconfig, reducer) + if satisfiable && !isnothing(result) + # Convert TNProblem result to variable assignments + assignments = Dict{Symbol, Int}() + for (i, symbol) in enumerate(sat.symbols) + assignments[symbol] = get_var_value(result, i) + end + return satisfiable, assignments, depth, stats + else + return false, Dict{Symbol, Int}(), depth, stats + end end -function solvebip(bip::BooleanInferenceProblem; bsconfig::BranchingStrategy = BranchingStrategy(table_solver = TNContractionSolver(), selector = KNeighborSelector(1,1), measure=NumOfVertices()), reducer=NoReducer()) - bs = initialize_branching_status(bip) - bs = deduction_reduce(bip,bs,collect(1:length(bip.he2v))) - ns,res,count_num = branch_and_reduce(bip,bs, bsconfig, reducer) - @show count_num - return ns,get_answer(res,bip.literal_num) -end \ No newline at end of file +function solve_factoring( + n::Int, m::Int, N::Int; + bsconfig::BranchingStrategy=BranchingStrategy(table_solver=TNContractionSolver(), selector=LeastOccurrenceSelector(1, 2), measure=NumUnfixedVars()), + reducer::AbstractReducer=NoReducer() +) + 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) + 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 diff --git a/src/knn.jl b/src/knn.jl new file mode 100644 index 0000000..5c9f5c2 --- /dev/null +++ b/src/knn.jl @@ -0,0 +1,124 @@ +# Local workspace for k-neighboring algorithm +mutable struct KNNWorkspace + epoch::Int + visited_vars::Vector{Int} + visited_tensors::Vector{Int} + frontier::Vector{Int} + next_frontier::Vector{Int} + collected_vars::Vector{Int} + collected_tensors::Vector{Int} +end + +function KNNWorkspace(var_num::Int, tensor_num::Int) + return KNNWorkspace( + 1, + zeros(Int, var_num), + zeros(Int, tensor_num), + Int[], + Int[], + Int[], + Int[] + ) +end + +@inline function _mark_var!(ws::KNNWorkspace, var::Int) + ws.visited_vars[var] = ws.epoch +end +@inline function _mark_tensor!(ws::KNNWorkspace, tensor::Int) + ws.visited_tensors[tensor] = ws.epoch +end +@inline function _check_seen_var(ws::KNNWorkspace, var::Int) + return ws.visited_vars[var] == ws.epoch +end +@inline function _check_seen_tensor(ws::KNNWorkspace, tensor::Int) + return ws.visited_tensors[tensor] == ws.epoch +end + +@inline function _bump_epoch!(ws::KNNWorkspace) + ws.epoch += 1 +end + +function expand_one_var!( + tn::TNStatic, + ws::KNNWorkspace, + doms::Vector{DomainMask}, + max_tensors::Int +)::Bool # returns `true` if STOPPED EARLY due to limits; `false` otherwise + empty!(ws.next_frontier) + @inbounds for var in ws.frontier + for tensor_id in tn.v2t[var] + if !_check_seen_tensor(ws, tensor_id) + _mark_tensor!(ws, tensor_id) + push!(ws.collected_tensors, tensor_id) + if length(ws.collected_tensors) >= max_tensors + return true + end + end + 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]) + _mark_var!(ws, var_id) + push!(ws.collected_vars, var_id) + push!(ws.next_frontier, var_id) + end + end + end + end + return false +end + +function classify_inner_boundary!(tn::TNStatic, ws::KNNWorkspace, vars::Vector{Int}) + inner = Int[] + boundary = Int[] + @inbounds for vid in vars + deg_total = tn.vars[vid].deg + deg_in = 0 + for tensor_id in tn.v2t[vid] + _check_seen_tensor(ws, tensor_id) && (deg_in += 1) + end + if deg_in == deg_total + push!(inner, vid) + else + push!(boundary, vid) + end + end + return inner, boundary +end + + +function k_neighboring( + tn::TNStatic, + doms::Vector{DomainMask}, + focus_var::Int; + max_tensors::Int, + k::Int = 2, +) + @debug "k_neighboring: focus_var = $focus_var" + @assert k ≥ 0 + @assert !is_fixed(doms[focus_var]) "Focus variable must be unfixed" + + # Create local workspace for this call + ws = KNNWorkspace(length(tn.vars), length(tn.tensors)) + + # Initialize with focus variable + _mark_var!(ws, focus_var) + push!(ws.frontier, focus_var) + push!(ws.collected_vars, focus_var) + + # k var-hops (only expanding to unfixed variables) + stopped = false + @inbounds for _ in 1:k + stopped = expand_one_var!(tn, ws, doms, max_tensors) + ws.frontier, ws.next_frontier = ws.next_frontier, ws.frontier + empty!(ws.next_frontier) + (stopped || isempty(ws.frontier)) && break + end + + # classify inner/boundary using current epoch marks for tensor-membership + inner, boundary = classify_inner_boundary!(tn, ws, ws.collected_vars) + + # deterministic ordering + sort!(inner); sort!(boundary); sort!(ws.collected_tensors) + return Region(focus_var, ws.collected_tensors, inner, boundary) +end \ No newline at end of file diff --git a/src/measure.jl b/src/measure.jl new file mode 100644 index 0000000..fb4cd52 --- /dev/null +++ b/src/measure.jl @@ -0,0 +1,16 @@ +struct NumUnfixedVars <: AbstractMeasure end +function OptimalBranchingCore.measure(problem::TNProblem, ::NumUnfixedVars) + # problem.n_unfixed == 0 && println("n_unfixed is 0") + return Int(problem.n_unfixed) +end + +# struct NumUnfixedTensors <: AbstractMeasure end +# function OptimalBranchingCore.measure(problem::TNProblem, ::NumUnfixedTensors) +# problem.doms +# end + +# struct NumUnfixedTensors <: AbstractMeasure end +# function OptimalBranchingCore.measure(problem::TNProblem, ::NumUnfixedTensors) +# length(problem.ws.active_tensors) +# end + diff --git a/src/problems.jl b/src/problems.jl new file mode 100644 index 0000000..8bf45e7 --- /dev/null +++ b/src/problems.jl @@ -0,0 +1,198 @@ +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 new file mode 100644 index 0000000..c0d27fc --- /dev/null +++ b/src/propagate.jl @@ -0,0 +1,330 @@ +""" + 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/readcnf.jl b/src/readcnf.jl deleted file mode 100644 index 8d1d946..0000000 --- a/src/readcnf.jl +++ /dev/null @@ -1,31 +0,0 @@ -function readcnf(filename::String) - f = open(filename) - vals = split(readline(f)) - tensors = Vector{Vector{Tropical{Float64}}}() - he2v = Vector{Vector{Int}}() - for line in eachline(f) - cl = filter(!iszero, parse.(Int, split(line))) - tensor, literals = clause2tensors(cl) - push!(tensors, tensor) - push!(he2v, literals) - end - @assert length(he2v) == parse(Int, vals[4]) - return BooleanInferenceProblem(tensors, he2v, parse(Int, vals[3])) -end - -function clause2tensors(clause::Vector{Int}) - tensor = ones(Tropical{Float64}, 2^length(clause)) - literals = Int[] - sum = 0 - for j in 1:length(clause) - sum += (1 << (j - 1)) * (sign(clause[j]) == 1 ? 0 : 1) - push!(literals, abs(clause[j])) - end - tensor[sum+1] = Tropical(-Inf) - return tensor,literals -end - -function solvecnf(filename::String; bsconfig::BranchingStrategy = BranchingStrategy(table_solver = TNContractionSolver(), selector = KNeighborSelector(1,1), measure=NumOfVertices()), reducer=NoReducer()) - p = readcnf(filename) - return solvebip(p; bsconfig, reducer) -end \ No newline at end of file diff --git a/src/reducer.jl b/src/reducer.jl deleted file mode 100644 index 6df53ed..0000000 --- a/src/reducer.jl +++ /dev/null @@ -1,105 +0,0 @@ -struct NoReducer <: AbstractReducer end -struct RuleReducer <: AbstractReducer # Rules from paper "Fast Exact Algorithms for the SAT Problem with Bounded Occurrences of Variables" - rules::Vector{Int} -end -function OptimalBranchingCore.reduce_problem(p::BooleanInferenceProblem, bs::AbstractBranchingStatus,reducer::NoReducer) - return bs -end -function OptimalBranchingCore.reduce_problem(p::BooleanInferenceProblem, bs::AbstractBranchingStatus,reducer::RuleReducer) - if 4 ∈ reducer.rules - bs = reduction_4(p, bs) - end - return bs -end -function reduction_4(p::BooleanInferenceProblem, bs::AbstractBranchingStatus) - # Rule 4: If x is true in all clauses, set x to true. - for (i, undecided_literal_num) in enumerate(bs.undecided_literals) - if undecided_literal_num == 1 - edge = p.he2v[i] - undecided_literal = findfirst(x -> readbit(bs.decided_mask, x) == 0, edge) - bs.decided_mask = bs.decided_mask | LongLongUInt(1) << (undecided_literal - 1) - bs.config = bs.config & (~(LongLongUInt(1) << (undecided_literal - 1))) - bs.undecided_literals[i] = -1 - end - end - return bs -end - -function deduction_reduce(p::BooleanInferenceProblem, bs::AbstractBranchingStatus, reducing_queue::Vector{Int}) - while !isempty(reducing_queue) - isempty(reducing_queue) && break - edge_num = popfirst!(reducing_queue) - (bs.undecided_literals[edge_num] <= 0) && continue - - zerocount, sumpos = check_reduce(p.he2v[edge_num], bs.decided_mask,bs.config, p.tensors[edge_num]) - (zerocount == 1) || continue - - bs, aedges = decide_literal(bs, p, p.he2v[edge_num], Clause(2^length(p.he2v[edge_num]) - 1, sumpos - 1)) - - reducing_queue = reducing_queue ∪ aedges - end - return bs -end - -function check_reduce(he2vi, mask, config, tensor) - count = 0 - sum = 0 - decided_literal_num = 0 - for j in 1:length(he2vi) - if readbit(mask, he2vi[j]) == 1 - sum += Int(readbit(config, he2vi[j])) * (1 << (j - 1)) - decided_literal_num += 1 - end - end - sumpos = 0 - for i in 0:2^(length(he2vi)-decided_literal_num)-1 - sum1 = sum - counti = 0 - for j in 1:length(he2vi) - if !(readbit(mask, he2vi[j]) == 1) - counti += 1 - sum1 += (1 << (j - 1)) * readbit(i, counti) - end - end - if tensor[sum1+1] == Tropical(0.0) - count += 1 - sumpos = sum1 + 1 - end - end - return count, sumpos -end - -function decide_literal(bs::AbstractBranchingStatus{C}, p::BooleanInferenceProblem, vertices::Vector{Int}, clause::Clause{N}) where {N,C} - config = copy(bs.config) - # mask = copy(bs.decided_mask) - dls = Int[] - for (k, v) in enumerate(vertices) - if readbit(clause.mask, k) == 1 && (readbit(bs.decided_mask, v) == 0) - push!(dls, v) - # mask = mask | LongLongUInt{C}(1) << (v - 1) - if readbit(clause.val, k) == 1 - config = config | LongLongUInt{C}(1) << (v - 1) - end - end - end - mask = bs.decided_mask | vec2lluint(dls, typeof(config)) - undecided_literals = copy(bs.undecided_literals) - aedges = Int[] # Edges that have been changed - for edge_num in mapreduce(v -> p.v2he[v], ∪, dls) - if bs.undecided_literals[edge_num] > 0 - zerocount, _ = check_reduce(p.he2v[edge_num], mask, config, p.tensors[edge_num]) - - decided_num = count(x -> x in dls, p.he2v[edge_num]) - # edge_mask = vec2lluint(p.he2v[edge_num], typeof(config)) - # decided_num = count_ones(edge_mask & mask) - count_ones(edge_mask & bs.decided_mask) - - if (zerocount == 2^(undecided_literals[edge_num] - decided_num)) - undecided_literals[edge_num] = -1 - else - undecided_literals[edge_num] -= decided_num - push!(aedges, edge_num) - end - end - end - return BranchingStatus(config, mask, undecided_literals), aedges -end \ No newline at end of file diff --git a/src/region.jl b/src/region.jl new file mode 100644 index 0000000..acdcbbe --- /dev/null +++ b/src/region.jl @@ -0,0 +1,80 @@ +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 index ed71ad4..ada49c2 100644 --- a/src/selector.jl +++ b/src/selector.jl @@ -1,112 +1,32 @@ -struct SubBIP{N} - vs::Vector{Int} - edges::Vector{Int} - outside_vs_ind::Vector{Int} - sub_tensors::Array{Tropical{Float64}, N} +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 -function SubBIP(p::BooleanInferenceProblem, bs::AbstractBranchingStatus, vs::Vector{Int}) - he2v, edge_list, decided_v = subhg(p, bs) - edges = sort([i for i in 1:length(he2v) if (he2v[i] ⊆ vs)]) - outside_vs_ind = [v for v in vs if any(i -> !(i in edges) && (v in he2v[i]), 1:length(he2v))] - return SubBIP{length(vs)}(vs, edges, outside_vs_ind, gen_sub_tensor(p, bs, vs, edges, he2v, edge_list)) -end - -struct KNeighborSelector <: AbstractSelector - k::Int - initial_vertex_strategy::Int # 1: maximum, 2: minimum,3: minimum weight -end - -struct Smallest2NeighborSelector <: AbstractSelector end -struct KaHyParSelector <: AbstractSelector - app_domain_size::Int -end - -function OptimalBranchingCore.select_variables(p::BooleanInferenceProblem, bs::AbstractBranchingStatus, m::M, selector::KaHyParSelector) where {M <: AbstractMeasure} - he2v, edge_list, decided_v = subhg(p, bs) - h = KaHyPar.HyperGraph(he2v2sparse(he2v)) - imbalance = 1-2*selector.app_domain_size/p.literal_num - @show imbalance - parts = KaHyPar.partition(h, 2; configuration=:edge_cut, imbalance) - - zero_num = count(x-> x ≈ 0,parts) - one_num = length(parts)-zero_num - - part0 = abs(zero_num-selector.app_domain_size) < abs(one_num-selector.app_domain_size) ? findall(iszero,parts) : findall(!iszero,parts) - return SubBIP(p,bs,part0) -end - -function OptimalBranchingCore.select_variables(p::BooleanInferenceProblem, bs::AbstractBranchingStatus, m::M, selector::KNeighborSelector) where {M <: AbstractMeasure} - he2v, edge_list, decided_v = subhg(p, bs) - undecided_literals = setdiff(1:p.literal_num, decided_v) - - v2he = [count(x -> i ∈ x, he2v) for i in undecided_literals] - index = argmin(x -> iszero(v2he[x]) ? Inf : v2he[x], 1:length(v2he)) - initial_v = undecided_literals[index] - # initial_v = selector.initial_vertex_strategy == 1 ? maximum(undecided_literals) : minimum(undecided_literals) - - vs, edges, outside_vs_ind = k_neighboring(he2v, initial_v, selector.k) - - return SubBIP{length(vs)}(vs, edge_list[edges], outside_vs_ind, gen_sub_tensor(p, bs, vs, edges, he2v, edge_list)) -end - -function OptimalBranchingCore.select_variables(p::BooleanInferenceProblem, bs::AbstractBranchingStatus, m::M, selector::Smallest2NeighborSelector) where {M <: AbstractMeasure} - he2v, edge_list, decided_v = subhg(p, bs) - undecided_literals = setdiff(1:p.literal_num, decided_v) - - minval = -Inf - local min_vs, min_edges, min_outside_vs_ind - for v in undecided_literals - vs, edges, outside_vs_ind = k_neighboring(he2v, v, 2) - if length(outside_vs_ind) / length(vs) > minval - minval = length(outside_vs_ind) / length(vs) - min_vs = vs - min_edges = edges - min_outside_vs_ind = outside_vs_ind - end - end - - return SubBIP{length(min_vs)}(min_vs, edge_list[min_edges], min_outside_vs_ind, gen_sub_tensor(p, bs, min_vs, min_edges, he2v, edge_list)) -end - - -function k_neighboring(he2v::Vector{Vector{Int}}, vs, k::Int) - for _ in 1:k-1 - vs = first(_neighboring(he2v, vs)) - end - vs, edges = _neighboring(he2v, vs) - - outside_vs_ind = [ind for ind in 1:length(vs) if any([vs[ind] ∈ v for v in he2v[setdiff(1:length(he2v), edges)]])] - return vs, edges, outside_vs_ind -end -_neighboring(he2v::Vector{Vector{Int}}, vs::Int) = _neighboring(he2v, [vs]) -function _neighboring(he2v::Vector{Vector{Int}}, vs::Vector{Int}) - edges = sort([i for i in 1:length(he2v) if !isempty(he2v[i] ∩ vs)]) - vs = sort(reduce(∪, he2v[edges])) - return vs, edges -end - -function subhg(bip::BooleanInferenceProblem, bs::AbstractBranchingStatus) - decided_v = [i for i in 1:bip.literal_num if readbit(bs.decided_mask, i) == 1] - return [setdiff(bip.he2v[e], decided_v) for e in 1:length(bip.he2v) if bs.undecided_literals[e] > 0], [e for e in 1:length(bip.he2v) if bs.undecided_literals[e] > 0], decided_v -end - -function gen_sub_tensor(p::BooleanInferenceProblem, bs::AbstractBranchingStatus, vs::Vector{Int}, edges::Vector{Int}, he2v::Vector{Vector{Int}}, edge_list::Vector{Int}) - eincode = EinCode(he2v[edges], vs) - optcode = optimize_code(eincode, uniformsize(eincode, 2), GreedyMethod()) - - sub_tensors = optcode([vec2tensor(slice_tensor(p.tensors[e], bs.decided_mask, bs.config, p.he2v[e])) for e in edge_list[edges]]...) - return sub_tensors -end - -function he2v2sparse(he2v::Vector{Vector{Int}}) - I = Int[] - J = Int[] - for i in 1:length(he2v) - for e in 1:length(he2v[i]) - push!(I, he2v[i][e]) - push!(J, i) - end - end - return sparse(I, J, ones(length(I))) -end diff --git a/src/status.jl b/src/status.jl deleted file mode 100644 index ac29029..0000000 --- a/src/status.jl +++ /dev/null @@ -1,20 +0,0 @@ -abstract type AbstractBranchingStatus{C} end -struct BranchingStatus{C} <: AbstractBranchingStatus{C} - config::LongLongUInt{C} - decided_mask::LongLongUInt{C} - undecided_literals::Vector{Int} # undecided literals in each clause, 0 means unsatisfiable, -1 means already satisfied -end - -Base.:+(bs1::BranchingStatus, bs2::BranchingStatus) = BranchingStatus(bs1.config, bs1.decided_mask, bs1.undecided_literals) - -struct BranchingStatusBranchCount{C} <: AbstractBranchingStatus{C} - config::LongLongUInt{C} - decided_mask::LongLongUInt{C} - undecided_literals::Vector{Int} # undecided literals in each clause, 0 means unsatisfiable, -1 means already satisfied - count::Int -end - -Base.:+(bs1::BranchingStatusBranchCount, bs2::BranchingStatusBranchCount) = BranchingStatusBranchCount(bs1.config, bs1.decided_mask, bs1.undecided_literals, bs1.count + bs2.count) -function get_answer(bs::BranchingStatus,n::Int) - return [Int(readbit(bs.config,i)) for i in 1:n] -end \ No newline at end of file diff --git a/src/stride.jl b/src/stride.jl deleted file mode 100644 index b64a023..0000000 --- a/src/stride.jl +++ /dev/null @@ -1,51 +0,0 @@ -function get_tensor_number(tensor::AbstractVector, vals::Vector{Int}) - return tensor[sum([(1<<(i-1))*(vals[i]-1) for i in 1:length(vals)])+1] -end - -bin_to_vec(x,n) = [ (x >> i) & 1 for i in 0:(n-1)] -function get_vals!(vals::Vector{Int},float_pos::Vector{Int},i::Int) - vals[float_pos] = bin_to_vec(i,length(float_pos)).+ 1 - return vals -end -function slice_tensor(tensor::AbstractVector, pos::Vector{Int}, pos_vals::Vector{Int},n::Int) - vals = fill(0,n) - vals[pos] = pos_vals - float_pos = setdiff(1:n,pos) - return [get_tensor_number(tensor,get_vals!(vals,float_pos,i)) for i in 0:2^(n-length(pos))-1] -end - -function slice_tensor(tensor::AbstractVector, mask::LongLongUInt, config::LongLongUInt, he2vi::Vector{Int}) - decided_v, decided_vals = lluint2vec(mask, config, he2vi) - return slice_tensor(tensor, _vertex_in_edge(he2vi, decided_v), decided_vals .+ 1, length(he2vi)) -end - - -function vec2tensor(vec::AbstractVector) - return reshape(vec, fill(2, Int(log2(length(vec))))...) -end - -function vec2lluint(vec::Vector{Int}, ::Type{T}) where T - return sum(i->T(1)<<(i-1), vec) -end -function lluint2vec(x::LongLongUInt,vals::LongLongUInt, vec::Vector{Int}) - pos = Int[] - pos_vals = Int[] - for i in vec - if readbit(x,i) == 1 - push!(pos, i) - push!(pos_vals, readbit(vals,i)) - end - end - return pos,pos_vals -end - -function _vertex_in_edge(he2vi, dls::Vector{Int}) - pos = Int[] - for i in 1:length(dls) - v = dls[i] - if v in he2vi - push!(pos, findfirst(==(v), he2vi)) - end - end - return pos -end \ No newline at end of file diff --git a/src/tablesolver.jl b/src/tablesolver.jl deleted file mode 100644 index 4970c19..0000000 --- a/src/tablesolver.jl +++ /dev/null @@ -1,22 +0,0 @@ -struct TNContractionSolver <: AbstractTableSolver end - -function OptimalBranchingCore.branching_table(bip::BooleanInferenceProblem, bs::AbstractBranchingStatus,solver::TNContractionSolver, subbip::SubBIP) - - out_vs_num = length(subbip.outside_vs_ind) - vs_num = length(subbip.vs) - ind_pos = [i ∈ subbip.outside_vs_ind ? findfirst(==(i), subbip.outside_vs_ind) : findfirst(==(i), setdiff(1:vs_num, subbip.outside_vs_ind)) for i in 1:vs_num] - possible_configurations = Vector{Vector{Bool}}[] - for i in 0:(2^out_vs_num-1) - answer = [i & (1 << j) != 0 for j in 0:out_vs_num-1] - out_index = [i ? 2 : 1 for i in answer] - vec = [i ∈ subbip.outside_vs_ind ? out_index[ind_pos[i]] : (:) for i in 1:vs_num] - in_indies = findall(==(Tropical(0.0)), subbip.sub_tensors[vec...]) - if length(in_indies) == 0 - continue - end - pcs = [[i ∈ subbip.outside_vs_ind ? out_index[ind_pos[i]] : in_index[ind_pos[i]] for i in 1:vs_num] .== fill(2, vs_num) for in_index in in_indies] - push!(possible_configurations, pcs) - end - # @show possible_configurations - isempty(possible_configurations) ? (return BranchingTable(0,[Int[]])) : (return BranchingTable(vs_num, possible_configurations)) -end \ No newline at end of file diff --git a/src/tropicalsvd.jl b/src/tropicalsvd.jl deleted file mode 100644 index bcded10..0000000 --- a/src/tropicalsvd.jl +++ /dev/null @@ -1,34 +0,0 @@ -function tropical_svd(C::Matrix{TropicalAndOr},k::Int;solver = IPSolver()) - # C = A * B - # A: m x k, B: k x n - # IP by JuMP - - m,n = size(C) - model = Model(solver.optimizer) - !solver.verbose && set_silent(model) - - @variable(model, 0 <= a[i = 1:m*k] <= 1, Int) # a[i,l] = a[(l-1)*m+i] - @variable(model, 0 <= b[i = 1:k*n] <= 1, Int) # b[l,j] = b[(j-1)*k+l] - @variable(model, 0 <= d[i = 1:m*k*n] <= 1, Int) # d[i,l,j] = d[(j-1)*k*m+(l-1)*m+i] - - for i in 1:m - for j in 1:n - if C[i,j].n - for l in 1:k - @constraint(model, d[(j-1)*k*m+(l-1)*m+i] <= a[(l-1)*m+i]) - @constraint(model, d[(j-1)*k*m+(l-1)*m+i] <= b[(j-1)*k+l]) - end - @constraint(model, sum(d[(j-1)*k*m+(l-1)*m+i] for l in 1:k) >= 1) - else - for l in 1:k - @constraint(model, d[(j-1)*k*m+(l-1)*m+i] + 1 >= a[(l-1)*m+i] + b[(j-1)*k+l]) - @constraint(model, d[(j-1)*k*m+(l-1)*m+i] <= 0) - end - end - end - end - @objective(model, Max, 1) - - optimize!(model) - return is_solved_and_feasible(model),reshape([TropicalAndOr(v ≈ 1.0) for v in value.(a)],m,k), reshape([TropicalAndOr(v ≈ 1.0) for v in value.(b)],k,n) -end \ No newline at end of file diff --git a/src/types.jl b/src/types.jl deleted file mode 100644 index c219f95..0000000 --- a/src/types.jl +++ /dev/null @@ -1,24 +0,0 @@ -struct BooleanInferenceProblem <: AbstractProblem - tensors::Vector{Vector{Tropical{Float64}}} - he2v::Vector{Vector{Int}} - v2he::Vector{Vector{Int}} - literal_num::Int -end - -BooleanInferenceProblem(tensor::Vector{Vector{Tropical{Float64}}}, he2v::Vector{Vector{Int}},literal_num::Int)= BooleanInferenceProblem(tensor,he2v,[findall(x->i in x, he2v) for i in 1:literal_num], literal_num) - -Base.copy(p::BooleanInferenceProblem) = BooleanInferenceProblem(copy(p.tensors), copy(p.he2v),copy(p.v2he), p.literal_num) - -struct NumOfVertices <: AbstractMeasure end -OptimalBranchingCore.measure(bs::AbstractBranchingStatus, ::NumOfVertices) = - count_ones(bs.decided_mask) - -struct NumOfClauses <: AbstractMeasure end -OptimalBranchingCore.measure(bs::AbstractBranchingStatus, ::NumOfClauses) = count( >=(0) ,bs.undecided_literals) - -function initialize_branching_status(p::BooleanInferenceProblem) - t = (p.literal_num-1) ÷ 64 + 1 - return BranchingStatus{t}(LongLongUInt{t}(0), LongLongUInt{t}(0), length.(p.he2v)) -end - -struct NumOfDegrees <: AbstractMeasure end -OptimalBranchingCore.measure(bs::AbstractBranchingStatus, ::NumOfDegrees) = sum(x -> x > 0 ? x : 0 ,bs.undecided_literals) \ No newline at end of file diff --git a/src/utils.jl b/src/utils.jl new file mode 100644 index 0000000..d63bf86 --- /dev/null +++ b/src/utils.jl @@ -0,0 +1,22 @@ +@inline function get_unfixed_vars(doms::Vector{DomainMask})::Vector{Int} + unfixed_vars = Int[] + @inbounds for (i, dm) in enumerate(doms) + if !is_fixed(dm) + push!(unfixed_vars, i) + end + end + return unfixed_vars +end + +# Convenience: compute number of unfixed variables quickly. +@inline function count_unfixed(doms::Vector{DomainMask})::Int + c::Int = 0 + @inbounds for dm in doms + if !is_fixed(dm) + c += 1 + end + end + return c +end + +bits_to_int(v::Vector{Bool}) = sum(b << (i - 1) for (i, b) in enumerate(v)) \ No newline at end of file diff --git a/test/bmf.jl b/test/bmf.jl deleted file mode 100644 index 9fcf879..0000000 --- a/test/bmf.jl +++ /dev/null @@ -1,19 +0,0 @@ -using TropicalNumbers -using Test -using BooleanInference - -@testset "MEBF" begin - # This is an approximate algorithm for Boolean matrix factorization (BMF) based on the MBF algorithm. - A = Matrix{TropicalAndOr}(BitMatrix(rand(0:1, 100, 100))) - # MBF function contains 3 keyword arguments: - # Dim: maximum number of patterns (k) - # Thres: between 0 and 1. A smaller t could achieve higher coverage with less number of patterns, while a larger t enables a more sparse decomposition of the input matrix with greater number of patterns. - # P: percentage of uncovered 1s - B, C = MEBF(A) -end - -@testset "factor analysis" begin - # This is an exact algorithm for Boolean matrix factorization (BMF) based on the MBF algorithm. - I = Matrix{TropicalAndOr}(BitMatrix(rand(0:1, 50, 50))) - C, D = find_formal_concepts(I) -end diff --git a/test/branch.jl b/test/branch.jl index e7bd7b7..a8dc47f 100644 --- a/test/branch.jl +++ b/test/branch.jl @@ -1,17 +1,45 @@ -using BooleanInference -using BooleanInference.GenericTensorNetworks -using BooleanInference.GenericTensorNetworks: ∧, ∨, ¬ using Test -using BooleanInference.GenericTensorNetworks.ProblemReductions -using BooleanInference.OptimalBranchingCore: apply_branch, branch_and_reduce, Clause, BranchingStrategy +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: BranchingStrategy, NoReducer +using OptimalBranchingCore: Clause +using OptimalBranchingCore: branching_table, branch_and_reduce, apply_branch +using ProblemReductions: Factoring, reduceto, CircuitSAT, read_solution +using GenericTensorNetworks +using TropicalNumbers: Tropical +# using Logging + +# ENV["JULIA_DEBUG"] = "BooleanInference" + +@testset "branch" begin + fproblem = Factoring(12, 12, 10395529) + # fproblem = Factoring(2, 2, 4) + circuit_sat = reduceto(CircuitSAT, fproblem) + problem = CircuitSAT(circuit_sat.circuit.circuit; use_constraints=true) + + tn = GenericTensorNetwork(problem) + 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) + if res != Tropical(-Inf) + @test has_last_branch_problem(tn_problem) + res = last_branch_problem(tn_problem) + @test res.n_unfixed == 0 + # @show circuit_sat.q, circuit_sat.p + + a = get_var_value(res, circuit_sat.q) + b = get_var_value(res, circuit_sat.p) + # @show a, b + @show bits_to_int(a), bits_to_int(b) + end +end @testset "apply_branch" begin @bools a b c d e f g cnf = ∧(∨(a, b, ¬d, ¬e), ∨(¬a, d, e, ¬f), ∨(f, g), ∨(¬b, c)) - bip, syms = cnf2bip(cnf) - bs = initialize_branching_status(bip) - bs = apply_branch(bip, bs, Clause(0b110, 0b100), [1, 2, 3]) - @test bs.undecided_literals == [2, -1, 2, -1] - @test bs.config == 4 - @test bs.decided_mask == 6 -end \ No newline at end of file + tn_problem = BooleanInference.setup_from_cnf(cnf) + new_problem, _ = apply_branch(tn_problem, Clause(0b110, 0b100), [1, 2, 3]) + @test new_problem.n_unfixed == 5 +end diff --git a/test/contraction.jl b/test/contraction.jl new file mode 100644 index 0000000..211c091 --- /dev/null +++ b/test/contraction.jl @@ -0,0 +1,414 @@ +using Test +using BooleanInference +using BooleanInference: setup_from_tensor_network, TNProblem, DynamicWorkspace, setup_problem, select_variables, get_cached_region, LeastOccurrenceSelector, 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 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, + [1, 2], + [3, 4], + [1, 2]) + + @test region.id == 1 + @test length(region.tensors) == 2 + @test length(region.inner_vars) == 2 + @test length(region.boundary_vars) == 2 +end + +@testset "tensor_unwrapping" begin + # Test 2x2 tensor (2^1 = 2) + @test_throws AssertionError tensor_unwrapping([1.0, 2.0, 3.0]) # not power of 2 + + # Test 2-element vector (2^1) + vec2 = [1.0, 2.0] + t2 = tensor_unwrapping(vec2) + @test size(t2) == (2,) + @test t2[1] == 1.0 + @test t2[2] == 2.0 + + # Test 4-element vector (2^2) + vec4 = [1.0, 2.0, 3.0, 4.0] + t4 = tensor_unwrapping(vec4) + @test size(t4) == (2, 2) + @test t4[1,1] == 1.0 + @test t4[2,1] == 2.0 + @test t4[1,2] == 3.0 + @test t4[2,2] == 4.0 + + # Test 8-element vector (2^3) + vec8 = collect(1.0:8.0) + t8 = tensor_unwrapping(vec8) + @test size(t8) == (2, 2, 2) + @test length(t8) == 8 +end + +@testset "slice_tensor basic" begin + # Test 1D tensor (single boolean variable) + # one(Tropical{Float64}) = 0.0ₜ = satisfied + # zero(Tropical{Float64}) = -Infₜ = unsatisfied + T1 = one(Tropical{Float64}) + T0 = zero(Tropical{Float64}) + tensor1d = [T1, T0] # Allows x=0, forbids x=1 + axis_vars = [1] + + # Allow both values - shape unchanged + doms_both = [DM_BOTH] + result = slicing(tensor1d, doms_both, axis_vars) + @test length(result) == 2 + @test result == [T1, T0] + + # Allow only 0 - slice to single element + doms_0 = [DM_0] + result = slicing(tensor1d, doms_0, axis_vars) + @test length(result) == 1 # Only 1 free variable value + @test result[1] == T1 # x=0, keeps original value + + # Allow only 1 - slice to single element + doms_1 = [DM_1] + result = slicing(tensor1d, doms_1, axis_vars) + @test length(result) == 1 # Only 1 free variable value + @test result[1] == T0 # x=1, original was unsatisfied +end + +@testset "slice_tensor 2D" begin + # Test 2D tensor (two boolean variables) + # tensor[i,j] represents variable assignment (i-1, j-1) + # Using Tropical: one=satisfied, zero=unsatisfied + T1 = one(Tropical{Float64}) + T0 = zero(Tropical{Float64}) + tensor2d = [T1, T1, T1, T0] # (0,0), (1,0), (0,1), (1,1) + axis_vars = [1, 2] + + # Allow both variables to take both values - shape unchanged + doms = [DM_BOTH, DM_BOTH] + result = slicing(tensor2d, doms, axis_vars) + @test length(result) == 4 + @test result == [T1, T1, T1, T0] + + # Fix first variable to 0, allow second to vary - becomes 1D + doms = [DM_0, DM_BOTH] + result = slicing(tensor2d, doms, axis_vars) + @test length(result) == 2 # Only x2 free: x2=0, x2=1 + @test result[1] == T1 # x1=0, x2=0 + @test result[2] == T1 # x1=0, x2=1 + + # Fix second variable to 1, allow first to vary - becomes 1D + doms = [DM_BOTH, DM_1] + result = slicing(tensor2d, doms, axis_vars) + @test length(result) == 2 # Only x1 free: x1=0, x1=1 + @test result[1] == T1 # x1=0, x2=1 + @test result[2] == T0 # x1=1, x2=1 + + # Fix both variables to (0,1) - becomes scalar + doms = [DM_0, DM_1] + result = slicing(tensor2d, doms, axis_vars) + @test length(result) == 1 # No free variables + @test result[1] == T1 # x1=0, x2=1 +end + +@testset "slice_tensor 3D" begin + # Test 3D tensor (three boolean variables) + # Only use one and zero of Tropical + T1 = one(Tropical{Float64}) + T0 = zero(Tropical{Float64}) + # Pattern: allow most, forbid some specific combinations + # Index order: (x1,x2,x3) = (0,0,0), (1,0,0), (0,1,0), (1,1,0), (0,0,1), (1,0,1), (0,1,1), (1,1,1) + tensor3d = [T1, T1, T0, T1, T1, T0, T1, T0] + axis_vars = [1, 2, 3] + + # Allow all - shape unchanged + doms = [DM_BOTH, DM_BOTH, DM_BOTH] + result = slicing(tensor3d, doms, axis_vars) + @test length(result) == 8 + @test result == tensor3d + + # Fix x1=1 and x3=0, allow x2 to vary - becomes 1D + doms = [DM_1, DM_BOTH, DM_0] + result = slicing(tensor3d, doms, axis_vars) + @test length(result) == 2 # Only x2 free: x2=0, x2=1 + # x1=1, x2=0, x3=0 -> index 010 binary (bit pattern) = index 2 (0-indexed) = tensor3d[2] = T1 + # x1=1, x2=1, x3=0 -> index 110 binary = index 6 (0-indexed) = tensor3d[7] = T1 + @test result[1] == T1 # x1=1, x2=0, x3=0 + @test result[2] == T1 # x1=1, x2=1, x3=0 + + # Fix all variables - becomes scalar + doms = [DM_0, DM_1, DM_1] + result = slicing(tensor3d, doms, axis_vars) + @test length(result) == 1 # No free variables + # x1=0, x2=1, x3=1 -> index 011 reversed = 110 binary = index 6 (0-indexed) = tensor3d[7] = T1 + @test result[1] == T1 +end + +@testset "DomainMask helpers" begin + @test has0(DM_BOTH) == true + @test has1(DM_BOTH) == true + @test has0(DM_0) == true + @test has1(DM_0) == false + @test has0(DM_1) == false + @test has1(DM_1) == true +end + +function AND_test() + T1 = one(Tropical{Float64}) + T0 = zero(Tropical{Float64}) + + # T[x1, x2, y] + T_and = Array{Tropical{Float64}}(undef, 2, 2, 2) + + for x1 in 0:1, x2 in 0:1, y in 0:1 + if y == (x1 & x2) + T_and[x1+1, x2+1, y+1] = T1 + else + T_and[x1+1, x2+1, y+1] = T0 + end + end + return T_and +end + +function NOT_test() + T0 = zero(Tropical{Float64}) + T1 = one(Tropical{Float64}) + T_not = Array{Tropical{Float64}}(undef, 2, 2) + for x in 0:1, y in 0:1 + if y != x + T_not[x+1, y+1] = T1 + else + T_not[x+1, y+1] = T0 + end + end + return T_not +end + +@testset "contract_tensors" begin + tensor1 = AND_test() + vector1 = vec(tensor1) + DOMs = DomainMask[DM_BOTH, DM_BOTH, DM_0] + sliced_tensor1 = slicing(vector1, DOMs, [1, 2, 3]) + reshaped_tensor1 = tensor_unwrapping(sliced_tensor1) + + tensor2 = NOT_test() + vector2 = vec(tensor2) + DOMs = DomainMask[DM_BOTH, DM_BOTH] + sliced_tensor2 = slicing(vector2, DOMs, [1, 2]) + @test size(sliced_tensor2) == size(vector2) + + result = contract_tensors([sliced_tensor1, sliced_tensor2], Vector{Int}[Int[1,2], Int[4,2]], Int[1,2,4]) + @test result[2,2,1] == zero(Tropical{Float64}) +end + +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 +end + +@testset "contract_region" begin + problem = generate_example_problem() + 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)) + + contracted, _ = contract_region(tn_static, region, tn_problem.doms) + @test contracted != nothing + @test length(size(contracted)) == length(region.boundary_vars) + length(region.inner_vars) + + table = branching_table(tn_problem, TNContractionSolver(), vcat([v for v in region.boundary_vars], [v for v in region.inner_vars])) + @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 + +@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}) + + contracted_1d = reshape([T1, T0], 2) + configs = extract_inner_configs(contracted_1d, 1) + @test length(configs) == 1 + @test configs[1] == [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) + + 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 + +@testset "combine_configs" begin + boundary_config = [true, false] + inner_configs = [[false, true], [true, false]] + + full_configs = combine_configs(boundary_config, inner_configs) + @test length(full_configs) == 2 + @test full_configs[1] == [true, false, false, true] + @test full_configs[2] == [true, false, true, false] + + boundary_config = [true] + inner_configs = [[false]] + + full_configs = combine_configs(boundary_config, inner_configs) + @test length(full_configs) == 1 + @test full_configs[1] == [true, false] + + boundary_config = Bool[] + inner_configs = [[true, true], [false, false]] + + full_configs = combine_configs(boundary_config, inner_configs) + @test length(full_configs) == 2 + @test full_configs[1] == [true, true] + @test full_configs[2] == [false, false] +end diff --git a/test/interface.jl b/test/interface.jl index 45bc6a8..1ba4f79 100644 --- a/test/interface.jl +++ b/test/interface.jl @@ -5,62 +5,127 @@ using Test using BooleanInference.GenericTensorNetworks.ProblemReductions using BooleanInference.OptimalBranchingCore: BranchingStrategy -@testset "cnf2bip" begin +@testset "setup_from_cnf" begin @bools a b c d e f g cnf = ∧(∨(a, b, ¬d, ¬e), ∨(¬a, d, e, ¬f), ∨(f, g), ∨(¬b, c), ∨(¬a)) - bip, syms = cnf2bip(cnf) - @test bip.he2v == [[1, 2, 3, 4], [1, 3, 4, 5], [5, 6], [2, 7], [1]] - @test bip.tensors[3][1] == zero(Tropical{Float64}) - @test bip.literal_num == 7 + + he2v = [] + tnproblem = setup_from_cnf(cnf) + for tensors in tnproblem.static.t2v + list = [] + for item in tensors + push!(list, item.var) + end + push!(he2v, list) + end + @test he2v == [[1, 2, 3, 4], [1, 3, 4, 5], [5, 6], [2, 7], [1]] + @show tnproblem.static.tensors[3].tensor[1] == zero(Tropical{Float64}) + @test tnproblem.n_unfixed == 6 end -@testset "cir2bip" begin +@testset "convert_circuit_to_bip" begin circuit = @circuit begin c = x ∧ y end push!(circuit.exprs, Assignment([:c],BooleanExpr(true))) - bip, syms = cir2bip(circuit) - @test bip.he2v == [[1, 2, 3],[1]] - @test bip.tensors == [vec(Tropical.([0.0 0.0; -Inf -Inf;;; 0.0 -Inf; -Inf 0.0])), vec(Tropical.([-Inf, 0.0]))] - @test bip.literal_num == 3 + tnproblem = setup_from_circuit(circuit) + he2v = [] + for tensors in tnproblem.static.t2v + list = [] + for item in tensors + push!(list, item.var) + end + push!(he2v, list) + end + @test he2v == [[1, 2, 3],[1]] + @test tnproblem.static.tensors[1].tensor == vec(Tropical.([0.0 0.0; -Inf -Inf;;; 0.0 -Inf; -Inf 0.0])) + @test tnproblem.static.tensors[2].tensor == [Tropical(-Inf), Tropical(0.0)] + # After initial propagation, all variables are fixed (problem is solved) + @test tnproblem.n_unfixed == 0 end -@testset "solve_sat" begin +@testset "solve_sat_with_assignments" begin @bools a b c d e f g cnf = ∧(∨(a, b, ¬d, ¬e), ∨(¬a, d, e, ¬f), ∨(f, g), ∨(¬b, c), ∨(¬a)) - sat = Satisfiability(cnf) - res,dict = solve_sat(sat) + sat = Satisfiability(cnf; use_constraints=true) + res, dict, _, stats = solve_sat_with_assignments(sat) @test res == true @test satisfiable(cnf, dict) == true + # Test that stats are recorded + @test stats.total_branches >= 0 + @test stats.max_depth >= 0 cnf = ∧(∨(a), ∨(a,¬c), ∨(d,¬b), ∨(¬c,¬d), ∨(a,e), ∨(a,e,¬c), ∨(¬a)) - sat = Satisfiability(cnf) - res,dict = solve_sat(sat) + sat = Satisfiability(cnf; use_constraints=true) + res, dict, _, stats = solve_sat_with_assignments(sat) @test res == false - @test satisfiable(cnf, dict) == false + # @test satisfiable(cnf, dict) == false + @test isempty(dict) + @test stats.total_branches >= 0 end @testset "solve_factoring" begin - a,b = solve_factoring(5,5,31*29) + a, b, stats = solve_factoring(5, 5, 31*29) @test a*b == 31*29 + @test stats.total_branches >= 0 + @test stats.total_subproblems >= 0 + println("Factoring stats: branches=$(stats.total_branches), subproblems=$(stats.total_subproblems), max_depth=$(stats.max_depth)") end -@testset "benchmark" begin - table_solver = TNContractionSolver() - reducer = NoReducer() - for selector in [KNeighborSelector(1, 1), KNeighborSelector(2, 1), KNeighborSelector(1, 2), KNeighborSelector(2, 2)] - for measure in [NumOfVertices(), NumOfClauses(), NumOfDegrees()] - println("$measure,$selector") - solve_factoring(8, 8, 1019 * 1021; bsconfig = BranchingStrategy(; table_solver, selector, measure), reducer) - end - end +@testset "branching_statistics" begin + # Test with a simple SAT problem + @bools a b c d + cnf = ∧(∨(a, b), ∨(¬a, c), ∨(¬b, d)) + sat = Satisfiability(cnf; use_constraints=true) + tn_problem = setup_from_sat(sat) + + # Test initial stats are zero + initial_stats = get_branching_stats(tn_problem) + @test initial_stats.total_branches == 0 + @test initial_stats.total_subproblems == 0 + @test initial_stats.max_depth == 0 + + # Solve and check stats are recorded + result, depth, stats = BooleanInference.solve(tn_problem, + BranchingStrategy(table_solver=TNContractionSolver(), + selector=LeastOccurrenceSelector(1, 2), + measure=NumUnfixedVars()), + NoReducer()) + + # Stats should have been recorded + @test stats.total_branches > 0 || result !== nothing # Either branched or solved immediately + @test stats.total_subproblems >= 0 # Each branch creates at least 2 subproblems + @test stats.max_depth >= 0 + @test stats.avg_branching_factor >= 0.0 + + # Print stats for debugging + println("\nBranching Statistics:") + print_branching_stats(tn_problem) + + # Test reset functionality + reset_branching_stats!(tn_problem) + reset_stats = get_branching_stats(tn_problem) + @test reset_stats.total_branches == 0 + @test reset_stats.total_subproblems == 0 + @test reset_stats.max_depth == 0 end -@testset "interface" begin - solve_factoring(8, 8, 1019 * 1021; bsconfig = BranchingStrategy(; table_solver= TNContractionSolver(), selector=KNeighborSelector(1, 1), measure=NumOfDegrees()), reducer= NoReducer()) - solve_factoring(5, 5, 899; bsconfig = BranchingStrategy(; table_solver= TNContractionSolver(), selector=KNeighborSelector(1, 1), measure=NumOfDegrees()), reducer= NoReducer()) +# @testset "benchmark" begin +# table_solver = TNContractionSolver() +# reducer = NoReducer() +# for selector in [] +# for measure in [NumOfVertices(), NumOfClauses(), NumOfDegrees()] +# println("$measure,$selector") +# solve_factoring(8, 8, 1019 * 1021; bsconfig = BranchingStrategy(; table_solver, selector, measure), reducer) +# end +# end +# end + +# @testset "interface" begin +# solve_factoring(8, 8, 1019 * 1021; bsconfig = BranchingStrategy(; table_solver= TNContractionSolver(), selector=KNeighborSelector(1, 1), measure=NumOfDegrees()), reducer= NoReducer()) +# solve_factoring(5, 5, 899; bsconfig = BranchingStrategy(; table_solver= TNContractionSolver(), selector=KNeighborSelector(1, 1), measure=NumOfDegrees()), reducer= NoReducer()) - bs = BranchingStrategy(table_solver = TNContractionSolver(), selector = KNeighborSelector(1, 1), set_cover_solver = BooleanInference.OptimalBranchingCore.GreedyMerge(),measure=NumOfDegrees()) +# bs = BranchingStrategy(table_solver = TNContractionSolver(), selector = KNeighborSelector(1, 1), set_cover_solver = BooleanInference.OptimalBranchingCore.GreedyMerge(),measure=NumOfDegrees()) - solve_factoring(8, 8, 1019 * 1021; bsconfig = bs, reducer= NoReducer()) -end \ No newline at end of file +# solve_factoring(8, 8, 1019 * 1021; bsconfig = bs, reducer= NoReducer()) +# end \ No newline at end of file diff --git a/test/knn.jl b/test/knn.jl new file mode 100644 index 0000000..c2b30fd --- /dev/null +++ b/test/knn.jl @@ -0,0 +1,55 @@ +using Test +using BooleanInference +using ProblemReductions: Factoring, reduceto, CircuitSAT +using GenericTensorNetworks +using BooleanInference: setup_from_cnf, k_neighboring +using BooleanInference.GenericTensorNetworks: ∧, ∨, ¬ + +function generate_example_problem() + fproblem = Factoring(12, 12, 10395529) + res = reduceto(CircuitSAT, fproblem) + problem = CircuitSAT(res.circuit.circuit; use_constraints=true) + return problem +end + +@testset "knn" begin + problem = generate_example_problem() + tn = GenericTensorNetwork(problem) + tn = BooleanInference.setup_from_tensor_network(tn) + doms = BooleanInference.init_doms(tn) + region = BooleanInference.k_neighboring(tn, doms, 1; k=2, max_tensors=10) + @show region +end + +@testset "original_test_for_neighboring" begin + @bools a b c d e f g + cnf = ∧(∨(b), ∨(a,¬c), ∨(d,¬b), ∨(¬c,¬d), ∨(a,e), ∨(a,e,¬c)) + problem = setup_from_cnf(cnf) + # Use unpropagated doms for testing k_neighboring + doms = BooleanInference.init_doms(problem.static) + region = k_neighboring(problem.static, doms, 1; max_tensors=100, k=1) + @test region.tensors == [1, 3] + @test Set(union(region.inner_vars, region.boundary_vars)) == Set([1, 4]) + + region = k_neighboring(problem.static, doms, 2; max_tensors=100, k=1) + @test region.tensors == [2, 5, 6] + @test Set(union(region.inner_vars, region.boundary_vars)) == Set([2, 3, 5]) +end + + +@testset "original_test_for_k_neighboring" begin + @bools a b c d e + cnf = ∧(∨(b), ∨(a,¬c), ∨(d,¬b), ∨(¬c,¬d), ∨(a,e), ∨(a,e,¬c)) + problem = setup_from_cnf(cnf) + # Use unpropagated doms for testing k_neighboring + doms = BooleanInference.init_doms(problem.static) + region = k_neighboring(problem.static, doms, 1; max_tensors=100, k=2) + @test region.tensors == [1, 3, 4] + @test Set(union(region.inner_vars, region.boundary_vars)) == Set([1, 3, 4]) + @test region.boundary_vars == [3] + + region = k_neighboring(problem.static, doms, 2; max_tensors=100, k=2) + @test region.tensors == [2, 4, 5, 6] + @test Set(union(region.inner_vars, region.boundary_vars)) == Set([2, 3, 4, 5]) + @test region.boundary_vars == [4] +end \ No newline at end of file diff --git a/test/problems.jl b/test/problems.jl new file mode 100644 index 0000000..3b7cdc1 --- /dev/null +++ b/test/problems.jl @@ -0,0 +1,67 @@ +using Test +using BooleanInference +using ProblemReductions: Factoring, reduceto, CircuitSAT +using GenericTensorNetworks +using BooleanInference: setup_from_tensor_network, TNProblem, DynamicWorkspace, setup_problem +using TropicalNumbers: Tropical + +function generate_example_problem() + fproblem = Factoring(2, 2, 6) + res = reduceto(CircuitSAT, fproblem) + problem = CircuitSAT(res.circuit.circuit; use_constraints=true) + return problem +end + +@testset "generate_example_problem" begin + problem = generate_example_problem() + tn = GenericTensorNetwork(problem) + problem = setup_from_tensor_network(tn) + @test length(problem.vars) == length(tn.problem.symbols) + @test length(problem.tensors) > 0 + @test length(problem.v2t) == length(problem.vars) + @test length(problem.t2v) == length(problem.tensors) + @show problem +end + +@testset "ids are just Int" begin + var_id = 1 + @test var_id == 1 + tensor_id = 1 + @test tensor_id == 1 +end + +@testset "setup_problem" begin + var_num = 2 + tensors_to_vars = [[1, 2], [2]] + tensor_data = [ + [Tropical(0.0), Tropical(0.0), Tropical(0.0), Tropical(1.0)], # AND: [0,0,0,1] + [Tropical(1.0), Tropical(0.0)] # NOT: [1,0] + ] + + tn = setup_problem(var_num, tensors_to_vars, tensor_data) + + @test length(tn.vars) == 2 + @test all(v.dom_size == 2 for v in tn.vars) + + @test length(tn.tensors) == 2 + @test length(tn.tensors[1].var_axes) == 2 + @test length(tn.tensors[2].var_axes) == 1 + + @test length(tn.v2t) == 2 + @test length(tn.v2t[1]) == 1 + @test length(tn.v2t[2]) == 2 + + @test length(tn.t2v) == 2 + @test length(tn.t2v[1]) == 2 + @test length(tn.t2v[2]) == 1 + + @test haskey(tn.axis_of_t, (1, 1)) + @test haskey(tn.axis_of_t, (1, 2)) + @test haskey(tn.axis_of_t, (2, 2)) +end + +@testset "setup_from_tensor_network" begin + tn = GenericTensorNetwork(generate_example_problem()) + tn_static = setup_from_tensor_network(tn) + tn_problem = TNProblem(tn_static) +end diff --git a/test/propagate.jl b/test/propagate.jl new file mode 100644 index 0000000..08603d1 --- /dev/null +++ b/test/propagate.jl @@ -0,0 +1,174 @@ +using Test +using BooleanInference +using TropicalNumbers +using BooleanInference: setup_from_cnf, propagate, has1, is_fixed, has0,setup_from_circuit +using ProblemReductions: @circuit, Assignment, BooleanExpr, Factoring, reduceto +using GenericTensorNetworks + +@testset "propagate" begin + # Test 1: Simple AND gate (x1 ∧ x2 = 1) + # Only one feasible config: (1, 1) + T1 = one(Tropical{Float64}) + T0 = zero(Tropical{Float64}) + @testset "Simple unit propagation - AND gate" begin + # Create a 2-variable AND constraint + # Tensor encodes: only (1,1) is feasible (Tropical(0.0)) + tensor_data = [ + T0, # (0,0) - infeasible + T0, # (1,0) - infeasible + T0, # (0,1) - infeasible + T1 # (1,1) - feasible + ] + + static = BooleanInference.setup_problem( + 2, + [[1, 2]], + [tensor_data] + ) + + # Initially both variables are unfixed + doms = BooleanInference.init_doms(static) + @test doms[1] == BooleanInference.DM_BOTH + @test doms[2] == BooleanInference.DM_BOTH + + # After propagation, both should be fixed to 1 + propagated = BooleanInference.propagate(static, doms) + @test propagated[1] == BooleanInference.DM_1 + @test propagated[2] == BooleanInference.DM_1 + end + + # Test 2: No unit propagation possible + @testset "No propagation - multiple solutions" begin + # OR gate: (0,0) is infeasible, others are feasible + tensor_data = [ + T0, # (0,0) - infeasible + T1, # (1,0) - feasible + T1, # (0,1) - feasible + T1 # (1,1) - feasible + ] + + static = BooleanInference.setup_problem( + 2, + [[1, 2]], + [tensor_data] + ) + + doms = BooleanInference.init_doms(static) + + # No unit propagation should occur + propagated = BooleanInference.propagate(static, doms) + @test propagated[1] == BooleanInference.DM_BOTH + @test propagated[2] == BooleanInference.DM_BOTH + end + + # Test 3: Partial assignment leading to unit propagation + @testset "Propagation after partial assignment" begin + # AND gate again + tensor_data = [ + T0, # (0,0) - infeasible + T0, # (1,0) - infeasible + T0, # (0,1) - infeasible + T1 # (1,1) - feasible + ] + + static = BooleanInference.setup_problem( + 2, + [[1, 2]], + [tensor_data] + ) + + # Fix x1 = 1 + doms = BooleanInference.init_doms(static) + doms[1] = BooleanInference.DM_1 + + # Propagation should fix x2 = 1 + propagated = BooleanInference.propagate(static, doms) + @test propagated[1] == BooleanInference.DM_1 + @test propagated[2] == BooleanInference.DM_1 + end + + # Test 4: Contradiction detection + @testset "Contradiction detection" begin + # AND gate with x1 = 0 should lead to contradiction + tensor_data = [ + T0, # (0,0) - infeasible + T0, # (1,0) - infeasible + T0, # (0,1) - infeasible + T1 # (1,1) - feasible + ] + + static = BooleanInference.setup_problem( + 2, + [[1, 2]], + [tensor_data] + ) + + # Fix x1 = 0 (contradicts the AND constraint) + doms = BooleanInference.init_doms(static) + doms[1] = BooleanInference.DM_0 + + # Propagation should detect contradiction + propagated = BooleanInference.propagate(static, doms) + @test all(d -> d.bits == 0x00, propagated) + end + + # Test 5: Chain propagation + @testset "Chain propagation" begin + # Two constraints: x1 = x2, x2 = x3 + # x1 = x2: (0,0) and (1,1) are feasible + tensor1 = [ + T1, # (0,0) - feasible + T0, # (1,0) - infeasible + T0, # (0,1) - infeasible + T1 # (1,1) - feasible + ] + + # x2 = x3: (0,0) and (1,1) are feasible + tensor2 = [ + T1, # (0,0) - feasible + T0, # (1,0) - infeasible + T0, # (0,1) - infeasible + T1 # (1,1) - feasible + ] + + static = BooleanInference.setup_problem( + 3, + [[1, 2], [2, 3]], + [tensor1, tensor2] + ) + + # Fix x1 = 1 + doms = BooleanInference.init_doms(static) + doms[1] = BooleanInference.DM_1 + + # Propagation should fix x2 = 1 and x3 = 1 + propagated = BooleanInference.propagate(static, doms) + @test propagated[1] == BooleanInference.DM_1 + @test propagated[2] == BooleanInference.DM_1 + @test propagated[3] == BooleanInference.DM_1 + end +end + +@testset "reduce_problem" begin + @bools a b c d e f g + cnf = ∧(∨(a, b, ¬d, ¬e), ∨(¬a, d, e, ¬f), ∨(f, g), ∨(¬b, c), ∨(¬a)) + problem = setup_from_cnf(cnf) + new_doms = propagate(problem.static, problem.doms) + @show new_doms[1] + @test has0(new_doms[1]) && is_fixed(new_doms[1]) == true + # TODO: undecided_literal has not been refactored yet + + @bools x1 x2 x3 x4 x5 + cnf = ∧(∨(x1), ∨(x2, ¬x3), ∨(x4, ¬x1), ∨(¬x3, ¬x4), ∨(x2, x5), ∨(x2, x5, ¬x3)) + problem = setup_from_cnf(cnf) + new_doms = propagate(problem.static, problem.doms) + # TODO: undecided_literal has not been refactored yet + + circuit = @circuit begin + c = x ∧ y + end + push!(circuit.exprs, Assignment([:c], BooleanExpr(true))) + problem = setup_from_circuit(circuit) + new_doms = propagate(problem.static, problem.doms) + @show new_doms +end \ No newline at end of file diff --git a/test/readcnf.jl b/test/readcnf.jl deleted file mode 100644 index 13b1d52..0000000 --- a/test/readcnf.jl +++ /dev/null @@ -1,36 +0,0 @@ -using BooleanInference -using BooleanInference.GenericTensorNetworks -using BooleanInference.GenericTensorNetworks: ∧, ∨, ¬ -using Test -using BooleanInference.GenericTensorNetworks.ProblemReductions - -@testset "clause2tensors" begin - @bools a b c d - cnf = ∧(∨(a, b, ¬c)) - bip, syms = cnf2bip(cnf) - tensor,lits = BooleanInference.clause2tensors([1,2,-3]) - @test bip.tensors[1]== tensor - - cnf = ∧(∨(¬a, b, ¬c, d)) - bip, syms = cnf2bip(cnf) - tensor,lits = BooleanInference.clause2tensors([-1,2,-3,4]) - @test bip.tensors[1]== tensor -end - -@testset "readcnf" begin - @bools a b c d e f - cnf = ∧(∨(a, b, c), ∨(d, e), ∨(a,c,¬e), ∨(¬b, c, e,f)) - bip, syms = cnf2bip(cnf) - bip2 = readcnf("datas/test.cnf") - @test bip.he2v == bip2.he2v - @test bip.tensors == bip2.tensors - @test bip.literal_num == bip2.literal_num - @test bip.v2he == bip2.v2he -end - -@testset "solvecnf" begin - solvecnf("datas/test.cnf") - - bs = BranchingStrategy(table_solver = TNContractionSolver(), selector = KNeighborSelector(1, 1), set_cover_solver = BooleanInference.OptimalBranchingCore.GreedyMerge(),measure=NumOfDegrees()) - solvecnf("datas/1.cnf"; bsconfig = bs) -end \ No newline at end of file diff --git a/test/reducer.jl b/test/reducer.jl deleted file mode 100644 index 2367922..0000000 --- a/test/reducer.jl +++ /dev/null @@ -1,60 +0,0 @@ -using BooleanInference -using BooleanInference.GenericTensorNetworks -using Test -using BooleanInference.GenericTensorNetworks: ∧, ∨, ¬ -using BooleanInference.OptimalBranchingCore: reduce_problem,apply_branch,Clause -import BooleanInference.ProblemReductions: @circuit, Assignment,BooleanExpr,Factoring,reduceto,CircuitSAT -using BooleanInference.OptimalBranchingCore.BitBasis - - -@testset "reduce_problem" begin - @bools a b c d e f g - cnf = ∧(∨(a, b, ¬d, ¬e), ∨(¬a, d, e, ¬f), ∨(f, g), ∨(¬b, c), ∨(¬a)) - bip, syms = cnf2bip(cnf) - bs = initialize_branching_status(bip) - bsnew = BooleanInference.deduction_reduce(bip, bs, collect(1:5)) - @test bsnew.decided_mask == 1 - @test bsnew.config == 0 - @test bsnew.undecided_literals == [3,-1,2,2,-1] - - @bools x1 x2 x3 x4 x5 - cnf = ∧(∨(x1), ∨(x2, ¬x3), ∨(x4, ¬x1), ∨(¬x3, ¬x4), ∨(x2, x5), ∨(x2, x5, ¬x3)) - bip, syms = cnf2bip(cnf) - bs = initialize_branching_status(bip) - bsnew = BooleanInference.deduction_reduce(bip, bs, collect(1:6)) - @test bsnew.decided_mask == 13 - @test bsnew.config == 9 - @test bsnew.undecided_literals == [-1,-1,-1,-1,2,-1] - - circuit = @circuit begin - c = x ∧ y - end - push!(circuit.exprs, Assignment([:c], BooleanExpr(true))) - bip, syms = cir2bip(circuit) - bs = initialize_branching_status(bip) - bsnew = BooleanInference.deduction_reduce(bip, bs, collect(1:2)) - @test bsnew.decided_mask == 7 - @test bsnew.config == 7 - @test bsnew.undecided_literals == [-1,-1] -end - -@testset "decide_literal" begin - @bools a b c d e f g - cnf = ∧(∨(a, b, ¬d, ¬e), ∨(¬a, d, e, ¬f), ∨(f, g), ∨(¬b, c)) - bip, syms = cnf2bip(cnf) - bs = initialize_branching_status(bip) - bs_new ,aedges = decide_literal(bs,bip,[1,2],Clause(0b11, 0b10)) - @test bs_new.undecided_literals == [-1,-1,2,1] - @test bs_new.config == 2 - @test bs_new.decided_mask == 3 - @test aedges == [4] -end - -@testset "benchmark" begin - fproblem = Factoring(3, 5,7* 17) - sat = reduceto(CircuitSAT,fproblem) - p,syms = BooleanInference.sat2bip(sat.circuit) - bs = initialize_branching_status(p) - - bs = apply_branch(p, bs, Clause(0b001, 0b000),[1,2,3]) -end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index df4015a..ca1ea47 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,31 +1,34 @@ using BooleanInference using Test +@testset "branch.jl" begin + include("branch.jl") +end + +@testset "contraction.jl" begin + include("contraction.jl") +end + @testset "interface.jl" begin include("interface.jl") end -@testset "reducer.jl" begin - include("reducer.jl") +@testset "knn.jl" begin + include("knn.jl") end -@testset "selector.jl" begin - include("selector.jl") +@testset "problems.jl" begin + include("problems.jl") end - -@testset "branch.jl" begin - include("branch.jl") +@testset "propagate.jl" begin + include("propagate.jl") end -@testset "status.jl" begin - include("status.jl") +@testset "selector.jl" begin + include("selector.jl") end @testset "tablesolver.jl" begin include("tablesolver.jl") -end - -@testset "stride.jl" begin - include("stride.jl") end \ No newline at end of file diff --git a/test/selector.jl b/test/selector.jl index caac323..8acd6bf 100644 --- a/test/selector.jl +++ b/test/selector.jl @@ -1,70 +1,92 @@ -using BooleanInference using Test -using BooleanInference.ProblemReductions -using BooleanInference.GenericTensorNetworks +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.OptimalBranchingCore: select_variables, apply_branch, Clause -@testset "subhg" begin - @bools a b c d e f g - cnf = ∧(∨(a, b, ¬d, ¬e), ∨(¬a, d, e, ¬f), ∨(f, g), ∨(¬b, c), ∨(¬a)) - bip, syms = cnf2bip(cnf) - bs = initialize_branching_status(bip) - bsnew = BooleanInference.deduction_reduce(bip, bs, collect(1:5)) - subhe2v,elist = subhg(bip,bsnew) - - @test subhe2v == [[2, 3, 4],[5, 6],[2, 7]] - @test elist == [1,3,4] +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 end -@testset "_neighboring" begin - @bools a b c d e - cnf = ∧(∨(b), ∨(a,¬c), ∨(d,¬b), ∨(¬c,¬d), ∨(a,e), ∨(a,e,¬c)) - bip,syms = cnf2bip(cnf) +@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 +end - subbip,edges = BooleanInference._neighboring(bip.he2v,1) - @test subbip == [1, 4] - @test edges == [1, 3] +@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) - subbip, edges = BooleanInference._neighboring(bip.he2v,2) - @test subbip == [2, 3, 5] - @test edges == [2, 5, 6] + 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]) + end + # Check that focus variable is properly in inner vars + @test first_var in region.inner_vars end -@testset "k_neighboring" begin - @bools a b c d e - cnf = ∧(∨(b), ∨(a,¬c), ∨(d,¬b), ∨(¬c,¬d), ∨(a,e), ∨(a,e,¬c)) - bip,syms = cnf2bip(cnf) - - subbip, edges, outside_vs_ind = k_neighboring(bip.he2v,1,2) - @test subbip == [1, 3, 4] - @test edges == [1, 3, 4] - @test outside_vs_ind == [2] +@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 - subbip, edges, outside_vs_ind = k_neighboring(bip.he2v,2,2) - @test subbip == [2, 3, 4, 5] - @test edges == [2, 4, 5, 6] - @test outside_vs_ind == [3] + 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 end -@testset "KNeighborSelector,neighbor_subbip" begin +@testset "KNeighborSelector" begin @bools a b c d e f g cnf = ∧(∨(a, b, ¬d, ¬e), ∨(¬a, d, e, ¬f), ∨(f, g), ∨(¬b, c)) - sat = Satisfiability(cnf) - bip,syms = sat2bip(sat) - bs = initialize_branching_status(bip) + 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 - subbip = select_variables(bip,bs, BooleanInference.NumOfVertices(),KNeighborSelector(4,3)) - @test subbip.vs == collect(1:7) - @test subbip.edges == collect(1:4) - @test subbip.outside_vs_ind == [] -end - -@testset "KaHyParSelector" begin - fproblem = Factoring(8, 8, 251*249) - res = reduceto(CircuitSAT,fproblem) - - bip,syms = sat2bip(res.circuit) - bs = initialize_branching_status(bip) - subbip = select_variables(bip,bs, BooleanInference.NumOfVertices(),KaHyParSelector(10)) + 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 diff --git a/test/status.jl b/test/status.jl deleted file mode 100644 index 145e42b..0000000 --- a/test/status.jl +++ /dev/null @@ -1,15 +0,0 @@ -using BooleanInference -using BooleanInference.GenericTensorNetworks -using BooleanInference.GenericTensorNetworks: ∧, ∨, ¬ -using Test -using BooleanInference.GenericTensorNetworks.ProblemReductions - -@testset "initialize_branching_status" begin - @bools a b c d e f g - cnf = ∧(∨(a, b, ¬d, ¬e), ∨(¬a, d, e, ¬f), ∨(f, g), ∨(¬b, c), ∨(¬a)) - bip, syms = cnf2bip(cnf) - branching_status = initialize_branching_status(bip) - @test branching_status.undecided_literals == [4,4,2,2,1] - @test branching_status.config == 0 - @test branching_status.decided_mask == 0 -end \ No newline at end of file diff --git a/test/stride.jl b/test/stride.jl deleted file mode 100644 index f65e65b..0000000 --- a/test/stride.jl +++ /dev/null @@ -1,39 +0,0 @@ -using BooleanInference -using Test - -@testset "vec2tensor" begin - a = rand(2,2,2,2) - b = vec(a) - c = vec2tensor(b) - @test a == c -end - -@testset "get_tensor_number" begin - a = rand(2,2,2,2) - b = vec(a) - test_tag = true - for i in CartesianIndices(a) - test_tag = (get_tensor_number(b, collect(i.I)) == a[i]) && test_tag - end - @test test_tag -end - -@testset "slice_tensor" begin - a = reshape(1:32, 2,2,2,2,2) - b = vec(a) - subb = slice_tensor(b, [2,4,5], [2,1,1], 5) - suba = a[:,2,:,1,1] - @test subb == vec(suba) -end - -@testset "vec2lluint" begin - a = LongLongUInt{1}(11) - avals = LongLongUInt{1}(10) - b,bvals = lluint2vec(a,avals, [2,3,4]) - @test b == [2,4] - @test bvals == [1,1] - - c,cvals = lluint2vec(a,avals, [1,2,3,4,5]) - @test a == vec2lluint(c, typeof(a)) - @test cvals == [0,1,1] -end \ No newline at end of file diff --git a/test/tablesolver.jl b/test/tablesolver.jl index eaff32a..fb01e12 100644 --- a/test/tablesolver.jl +++ b/test/tablesolver.jl @@ -2,26 +2,26 @@ using BooleanInference using Test using BooleanInference.GenericTensorNetworks using BooleanInference.GenericTensorNetworks: ∧, ∨, ¬ -using BooleanInference.OptimalBranchingCore: branching_table,select_variables +using BooleanInference.OptimalBranchingCore: select_variables using BooleanInference.OptimalBranchingCore.BitBasis -using BooleanInference:SubBIP +# 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 = cnf2bip(cnf) - bs = initialize_branching_status(bip) +# @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 = branching_table(bip, bs,BooleanInference.TNContractionSolver(), subbip) +# 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 +# @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 diff --git a/test/tropicalsvd.jl b/test/tropicalsvd.jl deleted file mode 100644 index 298e351..0000000 --- a/test/tropicalsvd.jl +++ /dev/null @@ -1,28 +0,0 @@ -using Test -using TropicalNumbers -using Random -using BooleanInference -using BooleanInference.OptimalBranchingCore - -@testset "tropical_svd" begin - m = 5 - n = 5 - Random.seed!(1234) - c = [TropicalAndOr(rand()>0.5) for i in 1:m, j in 1:n] - ans,a,b = tropical_svd(c,4) - @test ans == true - @test c == a*b -end - -@testset "tropical_svd" begin - m = 50 - k= 6 - n = 50 - Random.seed!(1234) - a = [TropicalAndOr(rand()>0.5) for i in 1:m, j in 1:k] - b = [TropicalAndOr(rand()>0.5) for i in 1:k, j in 1:n] - c = a*b - ans,a,b = tropical_svd(c,k) - @test ans == true - @test c == a*b -end \ No newline at end of file