Add JuMP interface (#1)

* Add JuMP interface

Project.toml

@@ -7,6 +7,7 @@ Cbc = "9961bab8-2fa3-5c5a-9d89-47fab24efd76"
 Clp = "e2554f3b-3117-50c0-817c-e040a3ddf72d"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"
+JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Scratch = "6c6a2e73-6563-6170-7368-637461726353"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
@@ -17,6 +18,7 @@ Clp = "1"
 GoGameGraphs = "1"
 Graphs = "1.4"
 HiGHS = "1.10"
+JuMP = "1.23.4"
 Scratch = "1.2"
 julia = "1.6"
 

src/FeedbackArcSets.jl

@@ -27,6 +27,7 @@ include("edge_subgraph.jl")
 include("optimization.jl")
 include("cbc.jl")
 include("highs.jl")
+include("jump.jl")
 include("dfs_fas.jl")
 include("greedy_fas.jl")
 include("pagerank_fas.jl")
@@ -140,6 +141,9 @@ function find_feedback_arc_set_ip(graph::EdgeSubGraph;
                                   max_iterations = typemax(Int),
                                   time_limit = typemax(Int),
                                   solver = "cbc",
+                                  solver_options = Dict(
+                                    "allowableGap" => 0,
+                                    "relativeGap" => 0),
                                   solver_time_limit = 10,
                                   log_level = 1,
                                   iteration_callback = print_iteration_data)
@@ -169,9 +173,11 @@ function find_feedback_arc_set_ip(graph::EdgeSubGraph;
         end
 
         solution = solve_IP(O; solver,
-                            seconds = solver_time,
-                            allowableGap = 0,
-                            logLevel = max(0, log_level - 1))
+                            solver_options = merge(Dict(solver_options),
+                                Dict(
+                                    "seconds" => solver_time,
+                                    "logLevel" => max(0, log_level - 1)),
+                                ))
 
         objbound = solution.attrs[:objbound]
 

src/highs.jl

@@ -5,15 +5,16 @@ using HiGHS: Highs_create, Highs_setStringOptionValue, Highs_setSolution,
              Highs_getModelStatus, HighsInt
 
 function solve_IP(::Val{:highs}, O::OptProblem, initial_solution = Int[],
-                  use_warmstart = true; options...)
+                  use_warmstart = true; solver_options, options...)
     model = Highs_create()
-    Highs_setStringOptionValue(model, "mip_rel_gap", "0")
-    for (name, value) in options
-        name = get(Dict(:seconds => "time_limit",
-                        :allowableGap => "mip_abs_gap"),
+    for (name, value) in solver_options
+        name = get(Dict("seconds" => "time_limit",
+                        "allowableGap" => "mip_abs_gap",
+                        "relativeGap"  => "mip_rel_gap",
+                        ),
                    name, name)
-        if name == :logLevel
-            name = :log_to_console
+        if name == "logLevel"
+            name = "log_to_console"
             value = value > 0
         end
         Highs_setStringOptionValue(model, string(name), string(value))

src/jump.jl

@@ -0,0 +1,80 @@
+import JuMP
+using JuMP: MOI, @variable, @constraint, @objective
+
+function solve_IP(::Val{:jump}, O::OptProblem, initial_solution = Int[],
+                  use_warmstart = true; solver_options, options...)
+
+    optimizer = pop!(solver_options, "optimizer")
+    model = JuMP.direct_model(MOI.instantiate(optimizer));
+    JuMP.set_string_names_on_creation(model, false)
+
+    # Pass through solver options:
+    set_JuMP_options!(model, Dict(solver_options))
+
+    A, lb, ub = add_cycle_constraints_to_formulation(O)
+    JuMP_loadproblem!(model, A, O.l, O.u, O.c, lb, ub)
+
+    if !isempty(initial_solution) && use_warmstart
+        JuMP.set_start_value.(model[:x], initial_solution)
+    end
+    JuMP.optimize!(model)
+
+    if !JuMP.is_solved_and_feasible(model)
+        error("Problem not solved or feasible.")
+    end
+
+    attrs = Dict()
+    attrs[:objbound] = JuMP.objective_bound(model)
+    attrs[:solver] = :ip
+    attrs[:model] = model
+    status = JuMP.termination_status(model) 
+    objval = JuMP.objective_value(model)
+    sol = JuMP.value.(model[:x])
+    solution = Solution(status, objval, sol, attrs)
+    return solution
+end
+
+function set_JuMP_options!(model::JuMP.Model, options_dict::Dict)
+
+    if haskey(options_dict, "seconds")
+        JuMP.set_time_limit_sec(model, Float64(pop!(options_dict, "seconds")))
+    end
+
+    if haskey(options_dict, "allowableGap")
+        JuMP.set_attribute(model, MOI.AbsoluteGapTolerance(), Float64(pop!(options_dict, "allowableGap")))
+    end
+
+    if haskey(options_dict, "relativeGap")
+        JuMP.set_attribute(model, MOI.RelativeGapTolerance(), Float64(pop!(options_dict, "relativeGap")))
+    end
+
+    if haskey(options_dict, "logLevel")
+        log_value = pop!(options_dict, "logLevel")
+        if isinteger(log_value) && (log_value < 0 || iszero(log_value))
+            @info "Solver logging set to silent."
+            JuMP.set_silent(model)
+        end
+    end
+
+    for (name, value) in options_dict
+        JuMP.set_attribute(model, string(name), value)
+    end
+
+    # for (name, value) in options
+    #     MOI.set(model, MOI.RawOptimizerAttribute(string(name)), value)
+    # end
+
+end
+
+function JuMP_loadproblem!(model, A, l, u, c, lb, ub)
+    m, n = size(A)  # row/constraints, columns/variables
+    @variable(model, x[i=1:n], Bin, lower_bound = l[i], upper_bound = u[i])
+    @objective(model, Min, c' * x)
+    @constraint(model, A*x .>= lb)
+    @constraint(model, A*x .<= ub)
+
+    # # If the wrapper supports Interval constraints, we can do:
+    # @constraint(model, lb .<= A*x .<= ub)
+
+    return nothing
+end

src/optimization.jl

@@ -34,10 +34,10 @@ mutable struct Solution
     attrs
 end
 
-function solve_IP(O::OptProblem; solver, initial_solution = Int[],
+function solve_IP(O::OptProblem; solver, solver_options, initial_solution = Int[],
                   use_warmstart = true, options...)
     solve_IP(Val(Symbol(solver)), O, initial_solution, use_warmstart;
-             options...)
+             solver_options, options...)
 end
 
 function solve_IP(solver::Val{T}, args...; kwargs...) where T