diff --git a/Project.toml b/Project.toml index 113c91d..4d536f5 100644 --- a/Project.toml +++ b/Project.toml @@ -8,21 +8,31 @@ BitBasis = "50ba71b6-fa0f-514d-ae9a-0916efc90dcf" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MLStyle = "d8e11817-5142-5d16-987a-aa16d5891078" PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d" +[weakdeps] +JuMP = "4076af6c-e467-56ae-b986-b466b2749572" + +[extensions] +IPSolverExt = "JuMP" + [compat] BitBasis = "0.9" DocStringExtensions = "0.9" Graphs = "1" InteractiveUtils = "1" +JuMP = "1" +LinearAlgebra = "1" MLStyle = "0.4" PrettyTables = "2" julia = "1.10" [extras] Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +SCIP = "82193955-e24f-5292-bf16-6f2c5261a85f" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Test", "Documenter"] +test = ["Test", "Documenter", "JuMP", "SCIP"] diff --git a/ext/IPSolverExt.jl b/ext/IPSolverExt.jl new file mode 100644 index 0000000..6fdd868 --- /dev/null +++ b/ext/IPSolverExt.jl @@ -0,0 +1,106 @@ +module IPSolverExt + +import JuMP +using ProblemReductions +using LinearAlgebra + +function Base.findmin(problem::AbstractProblem, solver::IPSolver) + return _find(problem, solver,true) +end +function Base.findmax(problem::AbstractProblem, solver::IPSolver) + return _find(problem, solver,false) +end + +# tag: true for min, false for max +function _find(problem::AbstractProblem, solver::IPSolver,tag::Bool) + @assert num_flavors(problem) == 2 "IPSolver only supports boolean variables" + 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) + @assert all(length(obj.variables) <= 1 for obj in objs) "IPSolver only supports objectives with at most 1 variables" + + # IP by JuMP + model = JuMP.Model(solver.optimizer) + !solver.verbose && JuMP.set_silent(model) + + JuMP.@variable(model, 0 <= x[i = 1:nsc] <= 1, Int) + + for con in cons + f_vec = findall(!,con.specification) + num_vars = length(con.variables) + for f in f_vec + JuMP.@constraint(model, sum(j-> iszero(combs[num_vars][f][j]) ? (1 - x[con.variables[j]]) : x[con.variables[j]], 1:num_vars) <= num_vars -1) + end + end + if isempty(objs) + JuMP.@objective(model, Min, 0) + else + obj_sum = sum(objs) do obj + (1-x[obj.variables[1]])*obj.specification[1] + x[obj.variables[1]]*obj.specification[2] + end + tag ? JuMP.@objective(model, Min, obj_sum) : JuMP.@objective(model, Max, obj_sum) + end + + JuMP.optimize!(model) + @assert JuMP.is_solved_and_feasible(model) "The problem is infeasible" + return round.(Int, JuMP.value.(x)) +end + +""" + minimal_set_cover(coverset::Vector{Int}, subsets::Vector{Vector{Int}}, optimizer) + +Solve the set cover problem: all elements in the coverset must be covered by the subsets. +The objective is to minimize the number of subsets used. + +# Arguments +- `coverset::Vector{Int}`: The set of all elements to be covered. +- `subsets::Vector{Vector{Int}}`: The set of subsets to choose from. +- `optimizer`: The optimizer to use, e.g. SCIP.Optimizer or HiGHS.Optimizer + +# Returns +- `Vector{Int}`: The indices of the subsets to choose. +""" +function minimal_set_cover(coverset::Vector{Int}, subsets::Vector{Vector{Int}}, weights::AbstractVector, optimizer, verbose::Bool=false) + # Remove subsets that cover not existing elements in the coverset + @assert all(set -> issubset(set, coverset), subsets) "subsets ($subsets) must not cover any elements absent in the coverset ($coverset)" + + # Create a JuMP model for exact set cover + model = JuMP.Model(optimizer) + !verbose && JuMP.set_silent(model) + + # Define binary variables for each subset (1 if selected, 0 if not) + n = length(subsets) + JuMP.@variable(model, x[1:n], Bin) + + # Each element in the coverset must be covered exactly once + for element in coverset + # Find all subsets containing this element + covering_subsets = [i for i in 1:n if element in subsets[i]] + + # Each element must be covered exactly once + JuMP.@constraint(model, sum(x[i] for i in covering_subsets) >= 1) + end + + # Minimize the number of subsets used (optional objective) + JuMP.@objective(model, Min, sum(weights[i]*x[i] for i in 1:n)) + # Solve the model + JuMP.optimize!(model) + + # Return the solution if feasible + @assert JuMP.is_solved_and_feasible(model) "The problem is infeasible" + return [i for i in 1:n if JuMP.value(x[i]) > 0.5] +end +minimal_set_cover(coverset::Vector{Int}, subsets::Vector{Vector{Int}}, optimizer, verbose::Bool=false) = minimal_set_cover(coverset, subsets, UnitWeight(length(subsets)), optimizer, verbose) + +function Base.findmin(problem::SetCovering, solver::IPSolver) + return minimal_set_cover(problem.elements, problem.sets, problem.weights, solver.optimizer, solver.verbose) +end + +function Base.findmax(problem::SetCovering, solver::IPSolver) + return minimal_set_cover(problem.elements, problem.sets, - problem.weights, solver.optimizer, solver.verbose) +end + +end \ No newline at end of file diff --git a/src/ProblemReductions.jl b/src/ProblemReductions.jl index 0ee2489..8876ccb 100644 --- a/src/ProblemReductions.jl +++ b/src/ProblemReductions.jl @@ -41,7 +41,7 @@ export LogicGadget, truth_table export ReductionSATTo3SAT export ReductionCircuitToSpinGlass, ReductionMaxCutToSpinGlass, ReductionSpinGlassToMaxCut, ReductionVertexCoveringToSetCovering, ReductionSatToColoring, ReductionSpinGlassToQUBO, ReductionQUBOToSpinGlass -export findbest, BruteForce +export findbest, BruteForce, AbstractSolver, IPSolver export CNF export ReductionSATToIndependentSet, ReductionSATToDominatingSet export ReductionIndependentSetToSetPacking, ReductionSetPackingToIndependentSet @@ -57,7 +57,7 @@ include("topology.jl") include("bitvector.jl") include("models/models.jl") include("rules/rules.jl") -include("bruteforce.jl") +include("solvers.jl") include("reduction_path.jl") include("deprecated.jl") diff --git a/src/bruteforce.jl b/src/solvers.jl similarity index 86% rename from src/bruteforce.jl rename to src/solvers.jl index b602b02..7954a2f 100644 --- a/src/bruteforce.jl +++ b/src/solvers.jl @@ -1,9 +1,11 @@ +abstract type AbstractSolver end + """ $TYPEDEF A brute force method to find the best configuration of a problem. """ -Base.@kwdef struct BruteForce{T} +Base.@kwdef struct BruteForce{T} <: AbstractSolver atol::T=eps(Float64) rtol::T=eps(Float64) end @@ -38,3 +40,10 @@ function _find!(compare, best_configs, configs, sizes, initial, atol, rtol) end end end + +# Interface for IPSolver +Base.@kwdef struct IPSolver <: AbstractSolver + optimizer # e.g. HiGHS.Optimizer + max_itr::Int = 20 + verbose::Bool = false +end \ No newline at end of file diff --git a/test/IPSolverExt.jl b/test/IPSolverExt.jl new file mode 100644 index 0000000..217458b --- /dev/null +++ b/test/IPSolverExt.jl @@ -0,0 +1,56 @@ +using Test +using JuMP +using ProblemReductions +using SCIP +using Graphs + +@testset "IPSolverExt" begin + # Test exact_set_cover with HiGHS optimizer + optimizer = SCIP.Optimizer + nflavor = 5 + subsets = [[1, 2], [2, 3], [3, 4], [4, 5]] + coverset = [1, 2, 3, 4, 5] + + # Test exact_set_cover with HiGHS optimizer + Ext = Base.get_extension(ProblemReductions, :IPSolverExt) + result = Ext.minimal_set_cover(coverset, subsets, optimizer) + @test result == [1, 2, 4] || result == [1, 3, 4] +end + +@testset "SetCovering" begin + problem = SetCovering([[1, 2], [2, 3], [3, 4], [4, 5]], [1, 2, 3, 4]) + @test findmin(problem, IPSolver(SCIP.Optimizer,20,false)) == [1, 2, 4] + @test findmax(problem, IPSolver(SCIP.Optimizer,20,false)) == [1, 2, 3, 4] +end + +@testset "IPSolver" begin + graph = smallgraph(:petersen) + problem = MaximalIS(graph) + @test findmin(problem, IPSolver(SCIP.Optimizer,20,false)) ∈ findmin(problem, BruteForce()) + + problem = IndependentSet(graph) + @test findmax(problem, IPSolver(SCIP.Optimizer,20,false)) ∈ findmax(problem, BruteForce()) + + fact3 = Factoring(2, 1, 3) + res3 = reduceto(CircuitSAT, fact3) + problem = CircuitSAT(res3.circuit.circuit; use_constraints=true) + @test findmin(problem, IPSolver(SCIP.Optimizer,20,false)) ∈ findmin(problem, BruteForce()) + best_config3 = findmin(problem, IPSolver(SCIP.Optimizer,20,false)) + assignment3 = Dict(zip(res3.circuit.symbols, best_config3)) + @test (2* assignment3[:p2]+ assignment3[:p1]) * assignment3[:q1] == 3 + + m1 = Matching(graph) + @test findmax(m1, IPSolver(SCIP.Optimizer,20,false)) ∈ findbest(m1, BruteForce()) +end + +@testset "Factoring" begin + function factoring(m,n,N,solver) + fact3 = Factoring(m, n, N) + res3 = reduceto(CircuitSAT, fact3) + problem = CircuitSAT(res3.circuit.circuit; use_constraints=true) + vals = findmin(problem, IPSolver(solver,20,true)) + return ProblemReductions.read_solution(fact3, [vals[res3.p]...,vals[res3.q]...]) + end + a,b = factoring(5,5,899,SCIP.Optimizer) + @test a*b == 899 +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index fbb41fb..7a01bfd 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,6 +6,10 @@ using Documenter include("bitvector.jl") end +@testset "solvers" begin + include("solvers.jl") +end + @testset "models" begin include("models/models.jl") end @@ -30,5 +34,8 @@ end include("deprecated.jl") end +@testset "IPSolverExt" begin + include("IPSolverExt.jl") +end DocMeta.setdocmeta!(ProblemReductions, :DocTestSetup, :(using ProblemReductions); recursive=true) Documenter.doctest(ProblemReductions; manual=false, fix=false) \ No newline at end of file diff --git a/test/solvers.jl b/test/solvers.jl new file mode 100644 index 0000000..67aa086 --- /dev/null +++ b/test/solvers.jl @@ -0,0 +1,12 @@ +using Test, ProblemReductions, Graphs + +@testset "BruteForce" begin + graph = smallgraph(:petersen) + problem = IndependentSet(graph) + solver = BruteForce() + res = findbest(problem, solver) + @test res == [[0, 0, 1, 0, 1, 1, 1, 0, 0, 0], [1, 0, 0, 1, 0, 0, 1, 1, 0, 0], [0, 1, 0, 0, 1, 0, 0, 1, 1, 0], [0, 1, 0, 1, 0, 1, 0, 0, 0, 1], [1, 0, 1, 0, 0, 0, 0, 0, 1, 1]] + solver = BruteForce() + res = findbest(problem, solver) + @test res == [[0, 0, 1, 0, 1, 1, 1, 0, 0, 0], [1, 0, 0, 1, 0, 0, 1, 1, 0, 0], [0, 1, 0, 0, 1, 0, 0, 1, 1, 0], [0, 1, 0, 1, 0, 1, 0, 0, 0, 1], [1, 0, 1, 0, 0, 0, 0, 0, 1, 1]] +end \ No newline at end of file