BB: Use CPXstrongbranch if optimizer is CPLEX

Project.toml

@@ -22,6 +22,7 @@ Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 ProgressBars = "49802e3a-d2f1-5c88-81d8-b72133a6f568"
 PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"

src/MIPLearn.jl

@@ -5,6 +5,7 @@
 module MIPLearn
 
 using PyCall
+using Requires
 
 global DynamicLazyConstraintsComponent = PyNULL()
 global JuMPSolver = PyNULL()

src/bb/BB.jl

@@ -4,6 +4,8 @@
 
 module BB
 
+using Requires
+
 frac(x) = x - floor(x)
 
 include("structs.jl")
@@ -19,4 +21,8 @@ include("varbranch/random.jl")
 include("varbranch/reliability.jl")
 include("varbranch/strong.jl")
 
+function __init__()
+    @require CPLEX = "a076750e-1247-5638-91d2-ce28b192dca0" include("cplex.jl")
+end
+
 end # module

src/bb/cplex.jl

@@ -0,0 +1,41 @@
+#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
+#  Copyright (C) 2020-2022, UChicago Argonne, LLC. All rights reserved.
+#  Released under the modified BSD license. See COPYING.md for more details.
+
+using CPLEX
+
+function _probe(
+    mip::MIP,
+    cpx::CPLEX.Optimizer,
+    var::Variable,
+    ::Float64,
+    ::Float64,
+    ::Float64,
+    itlim::Int,
+)::Tuple{Float64,Float64}
+    indices = [var.index - Cint(1)]
+    downobj, upobj, cnt = [0.0], [0.0], 1
+
+    status = CPXlpopt(cpx.env, cpx.lp)
+    status == 0 || error("CPXlpopt failed ($status)")
+
+    status = CPXstrongbranch(
+        cpx.env,
+        cpx.lp,
+        indices,
+        cnt,
+        downobj,
+        upobj,
+        itlim,
+    )
+    status == 0 || error("CPXstrongbranch failed ($status)")
+
+    return upobj[1] * mip.sense, downobj[1] * mip.sense
+end
+
+
+function _relax_integrality!(cpx::CPLEX.Optimizer)::Nothing
+    status = CPXchgprobtype(cpx.env, cpx.lp, CPLEX.CPXPROB_LP)
+    status == 0 || error("CPXchgprobtype failed ($status)")
+    return
+end

src/bb/lp.jl

@@ -28,7 +28,7 @@ end
 
 function load!(mip::MIP, prototype::JuMP.Model)
     @threads for t = 1:nthreads()
-        model = Model()
+        model = direct_model(mip.constructor)
         MOI.copy_to(model, backend(prototype))
         _replace_zero_one!(backend(model))
         if t == 1
@@ -38,7 +38,6 @@ function load!(mip::MIP, prototype::JuMP.Model)
             mip.sense = _get_objective_sense(backend(model))
         end
         _relax_integrality!(backend(model))
-        set_optimizer(model, mip.constructor)
         mip.optimizers[t] = backend(model)
         set_silent(model)
     end
@@ -236,23 +235,40 @@ function name(mip::MIP, var::Variable)::String
     return MOI.get(mip.optimizers[t], MOI.VariableName(), MOI.VariableIndex(var.index))
 end
 
-# convert(::Type{MOI.VariableIndex}, v::Variable) = MOI.VariableIndex(v.index)
-
 """
-    probe(mip::MIP, var, frac, lb, ub)::Tuple{Float64, Float64}
+    probe(mip::MIP, var, x, lb, ub, max_iterations)::Tuple{Float64, Float64}
 
 Suppose that the LP relaxation of `mip` has been solved and that `var` holds
-a fractional value `f`. This function returns two numbers corresponding,
+a fractional value `x`. This function returns two numbers corresponding,
 respectively, to the the optimal values of the LP relaxations having the
-constraints `ceil(frac) <= var <= ub` and `lb <= var <= floor(frac)` enforced.
+constraints `ceil(x) <= var <= ub` and `lb <= var <= floor(x)` enforced.
 If any branch is infeasible, the optimal value for that branch is Inf for
 minimization problems and -Inf for maximization problems.
 """
-function probe(mip::MIP, var::Variable, frac::Float64, lb::Float64, ub::Float64)::Tuple{Float64,Float64}
-    set_bounds!(mip, [var], [ceil(frac)], [ub])
-    status_up, obj_up = solve_relaxation!(mip)
-    set_bounds!(mip, [var], [lb], [floor(frac)])
-    status_down, obj_down = solve_relaxation!(mip)
+function probe(
+    mip::MIP,
+    var::Variable,
+    x::Float64,
+    lb::Float64,
+    ub::Float64,
+    max_iterations::Int,
+)::Tuple{Float64,Float64}
+    return _probe(mip, mip.optimizers[threadid()], var, x, lb, ub, max_iterations)
+end
+
+function _probe(
+    mip::MIP,
+    _,
+    var::Variable,
+    x::Float64,
+    lb::Float64,
+    ub::Float64,
+    ::Int,
+)::Tuple{Float64,Float64}
+    set_bounds!(mip, [var], [ceil(x)], [ceil(x)])
+    _, obj_up = solve_relaxation!(mip)
+    set_bounds!(mip, [var], [floor(x)], [floor(x)])
+    _, obj_down = solve_relaxation!(mip)
     set_bounds!(mip, [var], [lb], [ub])
     return obj_up * mip.sense, obj_down * mip.sense
 end

src/bb/varbranch/reliability.jl

@@ -15,6 +15,7 @@ Base.@kwdef mutable struct ReliabilityBranching <: VariableBranchingRule
     look_ahead::Int = 10
     n_sb_calls::Int = 0
     side_effect::Bool = true
+    max_iterations::Int = 1_000_000
 end
 
 function find_branching_var(
@@ -58,6 +59,7 @@ function find_branching_var(
                 var = var,
                 x = node.fractional_values[σ[i]],
                 side_effect = rule.side_effect,
+                max_iterations = rule.max_iterations
             )
         else
             score = pseudocost_scores[σ[i]]

src/bb/varbranch/strong.jl

@@ -17,6 +17,7 @@ Base.@kwdef struct StrongBranching <: VariableBranchingRule
     look_ahead::Int = 10
     max_calls::Int = 100
     side_effect::Bool = true
+    max_iterations::Int = 1_000_000
 end
 
 function find_branching_var(rule::StrongBranching, node::Node, pool::NodePool)::Variable
@@ -42,6 +43,7 @@ function find_branching_var(rule::StrongBranching, node::Node, pool::NodePool)::
             var = var,
             x = node.fractional_values[σ[i]],
             side_effect = rule.side_effect,
+            max_iterations = rule.max_iterations,
         )
         # @show name(node.mip, var), round(score[1], digits=2)
         if score > max_score
@@ -63,6 +65,7 @@ function _strong_branch_score(;
     var::Variable,
     x::Float64,
     side_effect::Bool,
+    max_iterations::Int,
 )::Tuple{Float64,Int}
 
     # Find current variable lower and upper bounds
@@ -77,11 +80,7 @@ function _strong_branch_score(;
     end
 
     obj_up, obj_down = 0, 0
-    try
-        obj_up, obj_down = probe(node.mip, var, x, var_lb, var_ub)
-    catch
-        @warn "strong branch error" var = var
-    end
+    obj_up, obj_down = probe(node.mip, var, x, var_lb, var_ub, max_iterations)
     obj_change_up = obj_up - node.obj
     obj_change_down = obj_down - node.obj
     if side_effect

test/Project.toml

@@ -2,6 +2,7 @@
 Cbc = "9961bab8-2fa3-5c5a-9d89-47fab24efd76"
 Clp = "e2554f3b-3117-50c0-817c-e040a3ddf72d"
 Conda = "8f4d0f93-b110-5947-807f-2305c1781a2d"
+CPLEX = "a076750e-1247-5638-91d2-ce28b192dca0"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Gurobi = "2e9cd046-0924-5485-92f1-d5272153d98b"

test/bb/lp_test.jl

@@ -37,7 +37,7 @@ function runtests(optimizer_name, optimizer; large = true)
             @test round(vals[3], digits = 6) == 0.248696
 
             # Probe (up and down are feasible)
-            probe_up, probe_down = BB.probe(mip, mip.int_vars[1], 0.5, 0.0, 1.0)
+            probe_up, probe_down = BB.probe(mip, mip.int_vars[1], 0.5, 0.0, 1.0, 1_000_000)
             @test round(probe_down, digits = 6) == 62.690000
             @test round(probe_up, digits = 6) == 62.714100
 
@@ -53,14 +53,6 @@ function runtests(optimizer_name, optimizer; large = true)
             @test status == :Optimal
             @test round(obj, digits = 6) == 62.714777
 
-            # Probe (up is infeasible, down is feasible)
-            BB.set_bounds!(mip, mip.int_vars[1:3], [1.0, 1.0, 0.0], [1.0, 1.0, 1.0])
-            status, obj = BB.solve_relaxation!(mip)
-            @test status == :Optimal
-            probe_up, probe_down = BB.probe(mip, mip.int_vars[3], 0.5, 0.0, 1.0)
-            @test round(probe_up, digits = 6) == Inf
-            @test round(probe_down, digits = 6) == 63.073992
-
             # Fix all binary variables to one, making problem infeasible
             N = length(mip.int_vars)
             BB.set_bounds!(mip, mip.int_vars, ones(N), ones(N))
@@ -105,9 +97,32 @@ function runtests(optimizer_name, optimizer; large = true)
 end
 
 @testset "BB" begin
-    @time runtests("Clp", Clp.Optimizer)
+    # @time runtests(
+    #     "Clp",
+    #     optimizer_with_attributes(
+    #         Clp.Optimizer,
+    #     ),
+    # )
+
     if is_gurobi_available
         using Gurobi
-        @time runtests("Gurobi", Gurobi.Optimizer)
+        @time runtests(
+            "Gurobi",
+            optimizer_with_attributes(
+                Gurobi.Optimizer,
+                "Threads" => 1,
+            )
+        )
+    end
+
+    if is_cplex_available
+        using CPLEX
+        @time runtests(
+            "CPLEX",
+            optimizer_with_attributes(
+                CPLEX.Optimizer,
+                "CPXPARAM_Threads" => 1,
+            ),
+        )
     end
 end

test/runtests.jl

@@ -7,6 +7,7 @@ using MIPLearn
 
 MIPLearn.setup_logger()
 const is_gurobi_available = ("GUROBI_HOME" in keys(ENV))
+const is_cplex_available = ("CPLEX_STUDIO_BINARIES" in keys(ENV))
 
 @testset "MIPLearn" begin
     include("fixtures/knapsack.jl")