Integer programming solver

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"]

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

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

src/bruteforce.jl → 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

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

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)

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