Add implementation of textbook branch-and-bound method

Manifest.toml

@@ -126,9 +126,9 @@ version = "1.1.1"
 
 [[DataStructures]]
 deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
-git-tree-sha1 = "4437b64df1e0adccc3e5d1adbc3ac741095e4677"
+git-tree-sha1 = "7d9d316f04214f7efdbb6398d545446e246eff02"
 uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
-version = "0.18.9"
+version = "0.18.10"
 
 [[DataValueInterfaces]]
 git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"

Project.toml

@@ -9,6 +9,7 @@ Cbc = "9961bab8-2fa3-5c5a-9d89-47fab24efd76"
 Clp = "e2554f3b-3117-50c0-817c-e040a3ddf72d"
 Conda = "8f4d0f93-b110-5947-807f-2305c1781a2d"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
@@ -18,7 +19,9 @@ MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
 PackageCompiler = "9b87118b-4619-50d2-8e1e-99f35a4d4d9d"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 
 [compat]

src/MIPLearn.jl

@@ -33,6 +33,8 @@ include("solvers/macros.jl")
 include("utils/benchmark.jl")
 include("utils/parse.jl")
 
+include("bb/BB.jl")
+
 function __init__()
     copy!(miplearn, pyimport("miplearn"))
     copy!(traceback, pyimport("traceback"))

src/bb/BB.jl

@@ -0,0 +1,22 @@
+#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
+#  Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
+#  Released under the modified BSD license. See COPYING.md for more details.
+
+module BB
+
+frac(x) = x - floor(x)
+
+include("structs.jl")
+
+include("nodepool.jl")
+include("optimize.jl")
+include("log.jl")
+include("lp.jl")
+include("varbranch/hybrid.jl")
+include("varbranch/infeasibility.jl")
+include("varbranch/pseudocost.jl")
+include("varbranch/random.jl")
+include("varbranch/reliability.jl")
+include("varbranch/strong.jl")
+
+end # module

src/bb/log.jl

@@ -0,0 +1,79 @@
+#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
+#  Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
+#  Released under the modified BSD license. See COPYING.md for more details.
+
+using Printf
+
+function print_progress_header(; detailed_output::Bool)
+    @printf(
+        "%8s  %9s %9s %13s %13s %13s %9s %8s",
+        "time",
+        "visited",
+        "pending",
+        "obj",
+        "primal-bound",
+        "dual-bound",
+        "gap",
+        "lp-iter"
+    )
+    if detailed_output
+        @printf(
+            " %6s %6s %-24s %6s %6s %6s",
+            "node",
+            "parent",
+            "branch-var",
+            "b-val",
+            "depth",
+            "iinfes"
+        )
+    end
+    println()
+    flush(stdout)
+end
+
+function print_progress(
+    pool::NodePool,
+    node::Node;
+    time_elapsed::Float64,
+    print_interval::Int,
+    detailed_output::Bool,
+    primal_update::Bool,
+)::Nothing
+    prefix = primal_update ? "*" : " "
+    if (pool.processed % print_interval == 0) || isempty(pool.pending) || primal_update
+        @printf(
+            "%8.2f %1s%9d %9d %13.6e %13.6e %13.6e %9.2e %8d",
+            time_elapsed,
+            prefix,
+            pool.processed,
+            length(pool.processing) + length(pool.pending),
+            node.obj * node.mip.sense,
+            pool.primal_bound * node.mip.sense,
+            pool.dual_bound * node.mip.sense,
+            pool.gap,
+            pool.mip.lp_iterations,
+        )
+        if detailed_output
+            if isempty(node.branch_variables)
+                branch_var_name = "---"
+                branch_value = "---"
+            else
+                branch_var_name = name(node.mip, last(node.branch_variables))
+                L = min(24, length(branch_var_name))
+                branch_var_name = branch_var_name[1:L]
+                branch_value = @sprintf("%.2f", last(node.branch_values))
+            end
+            @printf(
+                " %6d %6s %-24s %6s %6d %6d",
+                node.index,
+                node.parent === nothing ? "---" : @sprintf("%d", node.parent.index),
+                branch_var_name,
+                branch_value,
+                length(node.branch_variables),
+                length(node.fractional_variables)
+            )
+        end
+        println()
+        flush(stdout)
+    end
+end

src/bb/lp.jl

@@ -0,0 +1,212 @@
+#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
+#  Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
+#  Released under the modified BSD license. See COPYING.md for more details.
+
+import Base: values, convert
+using Base.Threads
+import Base.Threads: @threads, nthreads, threadid
+using JuMP
+using MathOptInterface
+const MOI = MathOptInterface
+
+function init(constructor)::MIP
+    return MIP(constructor, Any[nothing for t = 1:nthreads()], Variable[], 1.0, 0)
+end
+
+function read!(mip::MIP, filename::AbstractString)::Nothing
+    @threads for t = 1:nthreads()
+        model = read_from_file(filename)
+        mip.optimizers[t] = backend(model)
+        _replace_zero_one!(mip.optimizers[t])
+        if t == 1
+            _assert_supported(mip.optimizers[t])
+            mip.binary_variables = _get_binary_variables(mip.optimizers[t])
+            mip.sense = _get_objective_sense(mip.optimizers[t])
+        end
+        _relax_integrality!(mip.optimizers[t])
+        set_optimizer(model, mip.constructor)
+        set_silent(model)
+    end
+    return
+end
+
+function _assert_supported(optimizer::MOI.AbstractOptimizer)::Nothing
+    types = MOI.get(optimizer, MOI.ListOfConstraints())
+    for (F, S) in types
+        _assert_supported(F, S)
+    end
+end
+
+function _assert_supported(F::DataType, S::DataType)::Nothing
+    if F in [MOI.ScalarAffineFunction{Float64}, MOI.SingleVariable] && S in [
+        MOI.LessThan{Float64},
+        MOI.GreaterThan{Float64},
+        MOI.EqualTo{Float64},
+        MOI.Interval{Float64},
+    ]
+        return
+    end
+    if F in [MOI.SingleVariable] && S in [MOI.Integer, MOI.ZeroOne]
+        return
+    end
+    error("MOI constraint not supported: $F in $S")
+end
+
+function _get_objective_sense(optimizer::MOI.AbstractOptimizer)::Float64
+    sense = MOI.get(optimizer, MOI.ObjectiveSense())
+    if sense == MOI.MIN_SENSE
+        return 1.0
+    elseif sense == MOI.MAX_SENSE
+        return -1.0
+    else
+        error("objective sense not supported: $sense")
+    end
+end
+
+_bounds_constraint(v::Variable) =
+    MOI.ConstraintIndex{MOI.SingleVariable,MOI.Interval{Float64}}(v.index)
+
+function _replace_zero_one!(optimizer::MOI.AbstractOptimizer)::Nothing
+    constrs_to_delete = MOI.ConstraintIndex[]
+    funcs = MOI.SingleVariable[]
+    sets = Union{MOI.Interval,MOI.Integer}[]
+    for ci in
+        MOI.get(optimizer, MOI.ListOfConstraintIndices{MOI.SingleVariable,MOI.ZeroOne}())
+        func = MOI.get(optimizer, MOI.ConstraintFunction(), ci)
+        var = func.variable
+        push!(constrs_to_delete, ci)
+        push!(funcs, MOI.SingleVariable(var))
+        push!(funcs, MOI.SingleVariable(var))
+        push!(sets, MOI.Interval{Float64}(0.0, 1.0))
+        push!(sets, MOI.Integer())
+    end
+    MOI.delete(optimizer, constrs_to_delete)
+    MOI.add_constraints(optimizer, funcs, sets)
+    return
+end
+
+function _get_binary_variables(optimizer::MOI.AbstractOptimizer)::Vector{Variable}
+    vars = Variable[]
+    for ci in
+        MOI.get(optimizer, MOI.ListOfConstraintIndices{MOI.SingleVariable,MOI.Integer}())
+        func = MOI.get(optimizer, MOI.ConstraintFunction(), ci)
+        var = Variable(func.variable.value)
+
+        MOI.is_valid(optimizer, _bounds_constraint(var)) ||
+            error("$var is not interval-constrained")
+        interval = MOI.get(optimizer, MOI.ConstraintSet(), _bounds_constraint(var))
+        interval.lower == 0.0 || error("$var has lb != 0")
+        interval.upper == 1.0 || error("$var has ub != 1")
+
+        push!(vars, var)
+    end
+    return vars
+end
+
+function _relax_integrality!(optimizer::MOI.AbstractOptimizer)::Nothing
+    indices =
+        MOI.get(optimizer, MOI.ListOfConstraintIndices{MOI.SingleVariable,MOI.Integer}())
+    MOI.delete(optimizer, indices)
+end
+
+"""
+    solve_relaxation(mip::MIP)::Tuple{Symbol, Float64}
+
+Solve the linear relaxation of `mip` and returns a tuple containing the
+solution status (either :Optimal or :Infeasible) and the optimal objective
+value. If the problem is infeasible, the optimal value is Inf for minimization
+problems and -Inf for maximization problems..
+"""
+function solve_relaxation!(mip::MIP)::Tuple{Symbol,Float64}
+    t = threadid()
+    MOI.optimize!(mip.optimizers[t])
+    try
+        mip.lp_iterations += MOI.get(mip.optimizers[t], MOI.SimplexIterations())
+    catch
+        # ignore
+    end
+    status = MOI.get(mip.optimizers[t], MOI.TerminationStatus())
+    if status == MOI.OPTIMAL
+        obj = MOI.get(mip.optimizers[t], MOI.ObjectiveValue())
+        return :Optimal, obj * mip.sense
+    elseif status in [MOI.INFEASIBLE, MOI.INFEASIBLE_OR_UNBOUNDED]
+        return :Infeasible, Inf * mip.sense
+    end
+    error("unknown status: $status")
+end
+
+"""
+    values(mip::MIP, vars::Vector{Variable})::Array{Float64}
+
+Returns a vector `vals` which describes the current primal values for the
+decision variables `vars`. More specifically, `vals[j]` is the current
+primal value of `vars[j]`.
+"""
+function values(mip::MIP, vars::Vector{Variable})::Array{Float64}
+    return MOI.get(
+        mip.optimizers[threadid()],
+        MOI.VariablePrimal(),
+        convert.(MOI.VariableIndex, vars),
+    )
+end
+
+values(mip::MIP) =
+    values(mip, MOI.get(mip.optimizers[threadid()], MOI.ListOfVariableIndices()))
+
+"""
+    set_bounds!(mip::MIP,
+                vars::Vector{Variable},
+                lb::Array{Float64},
+                ub::Array{Float64})
+
+Modify the bounds of the given variables. More specifically, sets
+upper and lower bounds of `vars[j]` to `ub[j]` and `lb[j]`, respectively.
+"""
+function set_bounds!(
+    mip::MIP,
+    vars::Vector{Variable},
+    lb::Array{Float64},
+    ub::Array{Float64},
+)::Nothing
+    t = threadid()
+    MOI.delete(mip.optimizers[t], _bounds_constraint.(vars))
+    funcs = MOI.SingleVariable[]
+    sets = MOI.Interval[]
+    for j = 1:length(vars)
+        push!(funcs, MOI.SingleVariable(vars[j]))
+        push!(sets, MOI.Interval(lb[j], ub[j]))
+    end
+    MOI.add_constraints(mip.optimizers[t], funcs, sets)
+    return
+end
+
+"""
+    name(mip::MIP, var::Variable)::String
+
+Return the name of the decision variable `var`.
+"""
+function name(mip::MIP, var::Variable)::String
+    t = threadid()
+    return MOI.get(mip.optimizers[t], MOI.VariableName(), convert(MOI.VariableIndex, var))
+end
+
+convert(::Type{MOI.VariableIndex}, v::Variable) = MOI.VariableIndex(v.index)
+
+"""
+    probe(mip::MIP, var)::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,
+respectively, to the the optimal values of the LP relaxations having the
+constraints var=1 and var=0 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)::Tuple{Float64,Float64}
+    set_bounds!(mip, [var], [1.0], [1.0])
+    status_up, obj_up = solve_relaxation!(mip)
+    set_bounds!(mip, [var], [0.0], [0.0])
+    status_down, obj_down = solve_relaxation!(mip)
+    set_bounds!(mip, [var], [0.0], [1.0])
+    return obj_up * mip.sense, obj_down * mip.sense
+end

src/bb/nodepool.jl

@@ -0,0 +1,188 @@
+#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
+#  Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
+#  Released under the modified BSD license. See COPYING.md for more details.
+
+using Statistics
+using DataStructures
+import Base.Threads: threadid
+
+function take(
+    pool::NodePool;
+    suggestions::Array{Node} = [],
+    time_remaining::Float64,
+    gap_limit::Float64,
+    node_limit::Int,
+)::Union{Symbol,Node}
+    t = threadid()
+    lock(pool.lock) do
+        n_processing = length(pool.processing)
+        if (
+            (pool.gap < gap_limit) ||
+            (n_processing + pool.processed >= node_limit) ||
+            (time_remaining < 0)
+        )
+            return :END
+        end
+        if isempty(pool.pending)
+            if isempty(pool.processing)
+                return :END
+            else
+                return :WAIT
+            end
+        else
+            # If one of the suggested nodes is still pending, return it.
+            # This is known in the literature as plunging.
+            for s in suggestions
+                if s in keys(pool.pending)
+                    delete!(pool.pending, s)
+                    pool.processing[s] = s.obj
+                    return s
+                end
+            end
+
+            # If all suggestions have already been processed
+            # or pruned, find another node based on best bound.
+            node = dequeue!(pool.pending)
+            pool.processing[node] = node.obj
+            return node
+        end
+    end
+end
+
+function offer(
+    pool::NodePool;
+    parent_node::Union{Nothing,Node},
+    child_nodes::Vector{Node},
+    time_elapsed::Float64 = 0.0,
+    print_interval::Int = 100,
+    detailed_output::Bool = false,
+)::Nothing
+    lock(pool.lock) do
+        primal_update = false
+
+        # Update node.processing and node.processed
+        pool.processed += 1
+        if parent_node !== nothing
+            delete!(pool.processing, parent_node)
+        end
+
+        # Queue child nodes
+        for node in child_nodes
+            if node.status == :Infeasible
+                continue
+            end
+            if node.obj >= pool.primal_bound
+                continue
+            end
+            if isempty(node.fractional_variables)
+                pool.primal_bound = min(pool.primal_bound, node.obj)
+                primal_update = true
+                continue
+            end
+            pool.pending[node] = node.obj
+        end
+
+        # Update dual bound
+        pool.dual_bound = pool.primal_bound
+        if !isempty(pool.pending)
+            pool.dual_bound = min(pool.dual_bound, peek(pool.pending)[2])
+        end
+        if !isempty(pool.processing)
+            pool.dual_bound = min(pool.dual_bound, peek(pool.processing)[2])
+        end
+
+        # Update gap
+        if pool.primal_bound == pool.dual_bound
+            pool.gap = 0
+        else
+            pool.gap = abs((pool.primal_bound - pool.dual_bound) / pool.primal_bound)
+        end
+
+        if parent_node !== nothing
+            # Update branching variable history
+            branch_var = child_nodes[1].branch_variables[end]
+            offset = findfirst(isequal(branch_var), parent_node.fractional_variables)
+            x = parent_node.fractional_values[offset]
+            obj_change_up = child_nodes[1].obj - parent_node.obj
+            obj_change_down = child_nodes[2].obj - parent_node.obj
+            _update_var_history(
+                pool = pool,
+                var = branch_var,
+                x = x,
+                obj_change_down = obj_change_down,
+                obj_change_up = obj_change_up,
+            )
+
+            # Update global history
+            pool.history.avg_pseudocost_up =
+                mean(vh.pseudocost_up for vh in values(pool.var_history))
+            pool.history.avg_pseudocost_down =
+                mean(vh.pseudocost_down for vh in values(pool.var_history))
+
+            # Print progress
+            print_progress(
+                pool,
+                parent_node,
+                time_elapsed = time_elapsed,
+                print_interval = print_interval,
+                detailed_output = detailed_output,
+                primal_update = primal_update,
+            )
+        end
+    end
+    return
+end
+
+function _update_var_history(;
+    pool::NodePool,
+    var::Variable,
+    x::Float64,
+    obj_change_down::Float64,
+    obj_change_up::Float64,
+)::Nothing
+    # Create new history entry
+    if var ∉ keys(pool.var_history)
+        pool.var_history[var] = VariableHistory()
+    end
+    varhist = pool.var_history[var]
+
+    # Push fractional value
+    push!(varhist.fractional_values, x)
+
+    # Push objective value changes
+    push!(varhist.obj_change_up, obj_change_up)
+    push!(varhist.obj_change_down, obj_change_down)
+
+    # Push objective change ratios
+    f_up = x - floor(x)
+    f_down = ceil(x) - x
+    if isfinite(obj_change_up)
+        push!(varhist.obj_ratio_up, obj_change_up / f_up)
+    end
+    if isfinite(obj_change_down)
+        push!(varhist.obj_ratio_down, obj_change_down / f_down)
+    end
+
+    # Update variable pseudocosts
+    varhist.pseudocost_up = 0.0
+    varhist.pseudocost_down = 0.0
+    if !isempty(varhist.obj_ratio_up)
+        varhist.pseudocost_up = sum(varhist.obj_ratio_up) / length(varhist.obj_ratio_up)
+    end
+    if !isempty(varhist.obj_ratio_down)
+        varhist.pseudocost_down =
+            sum(varhist.obj_ratio_down) / length(varhist.obj_ratio_down)
+    end
+    return
+end
+
+function generate_indices(pool::NodePool, n::Int)::Vector{Int}
+    lock(pool.lock) do
+        result = Int[]
+        for i = 1:n
+            push!(result, pool.next_index)
+            pool.next_index += 1
+        end
+        return result
+    end
+end

src/bb/optimize.jl

@@ -0,0 +1,125 @@
+#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
+#  Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
+#  Released under the modified BSD license. See COPYING.md for more details.
+
+using Printf
+using Base.Threads
+import Base.Threads: @threads, nthreads, threadid
+
+function solve!(
+    mip::MIP;
+    time_limit::Float64 = Inf,
+    node_limit::Int = typemax(Int),
+    gap_limit::Float64 = 1e-4,
+    print_interval::Int = 5,
+    initial_primal_bound::Float64 = Inf,
+    detailed_output::Bool = false,
+    branch_rule::VariableBranchingRule = HybridBranching(),
+)::NodePool
+    time_initial = time()
+    pool = NodePool(mip = mip)
+    pool.primal_bound = initial_primal_bound
+    print_progress_header(detailed_output = detailed_output)
+
+    root_node = _create_node(mip)
+    if isempty(root_node.fractional_variables)
+        println("root relaxation is integer feasible")
+        pool.dual_bound = root_node.obj
+        pool.primal_bound = root_node.obj
+        return pool
+    end
+
+    offer(pool, parent_node = nothing, child_nodes = [root_node])
+    @threads for t = 1:nthreads()
+        child_one, child_zero, suggestions = nothing, nothing, Node[]
+        while true
+            time_elapsed = time() - time_initial
+            if child_one !== nothing
+                suggestions = Node[child_one, child_zero]
+            end
+            node = take(
+                pool,
+                suggestions = suggestions,
+                time_remaining = time_limit - time_elapsed,
+                node_limit = node_limit,
+                gap_limit = gap_limit,
+            )
+            if node == :END
+                break
+            elseif node == :WAIT
+                sleep(0.1)
+                continue
+            else
+                ids = generate_indices(pool, 2)
+                branch_variable = find_branching_var(branch_rule, node, pool)
+                child_zero = _create_node(
+                    mip,
+                    index = ids[1],
+                    parent = node,
+                    branch_variable = branch_variable,
+                    branch_value = 0.0,
+                )
+                child_one = _create_node(
+                    mip,
+                    index = ids[2],
+                    parent = node,
+                    branch_variable = branch_variable,
+                    branch_value = 1.0,
+                )
+                offer(
+                    pool,
+                    parent_node = node,
+                    child_nodes = [child_one, child_zero],
+                    time_elapsed = time_elapsed,
+                    print_interval = print_interval,
+                    detailed_output = detailed_output,
+                )
+            end
+        end
+    end
+    return pool
+end
+
+function _create_node(
+    mip;
+    index::Int = 0,
+    parent::Union{Nothing,Node} = nothing,
+    branch_variable::Union{Nothing,Variable} = nothing,
+    branch_value::Union{Nothing,Float64} = nothing,
+)::Node
+    if parent === nothing
+        branch_variables = Variable[]
+        branch_values = Float64[]
+        depth = 1
+    else
+        branch_variables = [parent.branch_variables; branch_variable]
+        branch_values = [parent.branch_values; branch_value]
+        depth = parent.depth + 1
+    end
+    set_bounds!(mip, branch_variables, branch_values, branch_values)
+    status, obj = solve_relaxation!(mip)
+    if status == :Optimal
+        vals = values(mip, mip.binary_variables)
+        fractional_indices =
+            [j for j in 1:length(mip.binary_variables) if 1e-6 < vals[j] < 1 - 1e-6]
+        fractional_values = vals[fractional_indices]
+        fractional_variables = mip.binary_variables[fractional_indices]
+    else
+        fractional_variables = Variable[]
+        fractional_values = Float64[]
+    end
+    n_branch = length(branch_variables)
+    set_bounds!(mip, branch_variables, zeros(n_branch), ones(n_branch))
+    return Node(
+        mip,
+        index,
+        depth,
+        obj,
+        status,
+        branch_variables,
+        branch_values,
+        fractional_variables,
+        fractional_values,
+        parent,
+    )
+end

src/bb/structs.jl

@@ -0,0 +1,70 @@
+#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
+#  Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
+#  Released under the modified BSD license. See COPYING.md for more details.
+
+using DataStructures
+abstract type VariableBranchingRule end
+
+struct Variable
+    index::Any
+end
+
+mutable struct MIP
+    constructor
+    optimizers::Vector
+    binary_variables::Vector{Variable}
+    sense::Float64
+    lp_iterations::Int64
+end
+
+struct Node
+    mip::MIP
+    index::Int
+    depth::Int
+    obj::Float64
+    status::Symbol
+    branch_variables::Array{Variable}
+    branch_values::Array{Float64}
+    fractional_variables::Array{Variable}
+    fractional_values::Array{Float64}
+    parent::Union{Nothing,Node}
+end
+
+Base.@kwdef mutable struct History
+    avg_pseudocost_up::Float64 = 1.0
+    avg_pseudocost_down::Float64 = 1.0
+end
+
+mutable struct VariableHistory
+    fractional_values::Array{Float64}
+    obj_change_up::Array{Float64}
+    obj_change_down::Array{Float64}
+    obj_ratio_up::Array{Float64}
+    obj_ratio_down::Array{Float64}
+    pseudocost_up::Float64
+    pseudocost_down::Float64
+
+    VariableHistory() = new(
+        Float64[], # fractional_values
+        Float64[], # obj_change_up
+        Float64[], # obj_change_down
+        Float64[], # obj_ratio_up
+        Float64[], # obj_ratio_up
+        0.0, # pseudocost_up
+        0.0, # pseudocost_down
+    )
+end
+
+Base.@kwdef mutable struct NodePool
+    mip::MIP
+    pending::PriorityQueue{Node,Float64} = PriorityQueue{Node,Float64}()
+    processing::PriorityQueue{Node,Float64} = PriorityQueue{Node,Float64}()
+    processed::Int = 0
+    next_index::Int = 1
+    lock::ReentrantLock = ReentrantLock()
+    primal_bound::Float64 = Inf
+    dual_bound::Float64 = Inf
+    gap::Float64 = Inf
+    history::History = History()
+    var_history::Dict{Variable,VariableHistory} = Dict()
+end

src/bb/varbranch/hybrid.jl

@@ -0,0 +1,31 @@
+#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
+#  Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
+#  Released under the modified BSD license. See COPYING.md for more details.
+
+"""
+    HybridBranching(depth_cutoff::Int,
+                    shallow_rule::VariableBranchingRule,
+                    deep_rule::::VariableBranchingRule)
+
+Branching strategy that switches between two variable branching strategies,
+according to the depth of the node.
+
+More specifically, if `node.depth <= depth_cutoff`, then `shallow_rule` is
+applied. Otherwise, `deep_rule` is applied.
+"""
+mutable struct HybridBranching <: VariableBranchingRule
+    depth_cutoff::Int
+    shallow_rule::VariableBranchingRule
+    deep_rule::VariableBranchingRule
+end
+
+HybridBranching(depth_cutoff::Int = 10) =
+    HybridBranching(depth_cutoff, StrongBranching(), PseudocostBranching())
+
+function find_branching_var(rule::HybridBranching, node::Node, pool::NodePool)::Variable
+    if node.depth <= rule.depth_cutoff
+        return find_branching_var(rule.shallow_rule, node, pool)
+    else
+        return find_branching_var(rule.deep_rule, node, pool)
+    end
+end

src/bb/varbranch/infeasibility.jl

@@ -0,0 +1,54 @@
+#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
+#  Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
+#  Released under the modified BSD license. See COPYING.md for more details.
+
+"""
+    FirstInfeasibleBranching()
+
+Branching rule that always selects the first fractional variable.
+"""
+struct FirstInfeasibleBranching <: VariableBranchingRule end
+
+function find_branching_var(
+    rule::FirstInfeasibleBranching,
+    node::Node,
+    pool::NodePool,
+)::Variable
+    return node.fractional_variables[1]
+end
+
+"""
+    LeastInfeasibleBranching()
+
+Branching strategy that select the fractional variable whose value is the closest
+to an integral value.
+"""
+struct LeastInfeasibleBranching <: VariableBranchingRule end
+
+function find_branching_var(
+    rule::LeastInfeasibleBranching,
+    node::Node,
+    pool::NodePool,
+)::Variable
+    scores = [max(v - floor(v), ceil(v) - v) for v in node.fractional_values]
+    _, max_offset = findmax(scores)
+    return node.fractional_variables[max_offset]
+end
+
+"""
+    MostInfeasibleBranching()
+
+Branching strategy that selects the fractional variable whose value is closest
+to 1/2.
+"""
+struct MostInfeasibleBranching <: VariableBranchingRule end
+
+function find_branching_var(
+    rule::MostInfeasibleBranching,
+    node::Node,
+    pool::NodePool,
+)::Variable
+    scores = [min(v - floor(v), ceil(v) - v) for v in node.fractional_values]
+    _, max_offset = findmax(scores)
+    return node.fractional_variables[max_offset]
+end

src/bb/varbranch/pseudocost.jl

@@ -0,0 +1,45 @@
+#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
+#  Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
+#  Released under the modified BSD license. See COPYING.md for more details.
+
+"""
+    PseudocostBranching()
+
+Branching strategy that uses historical changes in objective value to estimate
+strong branching scores at lower computational cost.
+"""
+struct PseudocostBranching <: VariableBranchingRule end
+
+function find_branching_var(rule::PseudocostBranching, node::Node, pool::NodePool)::Variable
+    scores = [
+        _pseudocost_score(
+            node,
+            pool,
+            node.fractional_variables[j],
+            node.fractional_values[j],
+        ) for j = 1:length(node.fractional_variables)
+    ]
+    _, max_offset = findmax(scores)
+    return node.fractional_variables[max_offset]
+end
+
+function _pseudocost_score(
+    node::Node,
+    pool::NodePool,
+    var::Variable,
+    x::Float64,
+)::Tuple{Float64,Int}
+    f_up = x - floor(x)
+    f_down = ceil(x) - x
+    pseudo_up = pool.history.avg_pseudocost_up * f_up
+    pseudo_down = pool.history.avg_pseudocost_down * f_down
+    if var in keys(pool.var_history)
+        if isfinite(pool.var_history[var].pseudocost_up)
+            pseudo_up = pool.var_history[var].pseudocost_up * f_up
+        end
+        if isfinite(pool.var_history[var].pseudocost_down)
+            pseudo_down = pool.var_history[var].pseudocost_down * f_down
+        end
+    end
+    return (pseudo_up * f_up * pseudo_down * f_down, var.index)
+end

src/bb/varbranch/random.jl

@@ -0,0 +1,17 @@
+#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
+#  Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
+#  Released under the modified BSD license. See COPYING.md for more details.
+
+using Random
+
+
+"""
+    RandomBranching()
+
+Branching strategy that picks a fractional variable randomly.
+"""
+struct RandomBranching <: VariableBranchingRule end
+
+function find_branching_var(rule::RandomBranching, node::Node, pool::NodePool)::Variable
+    return shuffle(node.fractional_variables)[1]
+end

src/bb/varbranch/reliability.jl

@@ -0,0 +1,78 @@
+#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
+#  Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
+#  Released under the modified BSD license. See COPYING.md for more details.
+
+"""
+    ReliabilityBranching
+
+Branching strategy that uses pseudocosts if there are enough observations
+to make an accurate prediction of strong branching scores, or runs the
+actual strong branching routine if not enough observations are available.
+"""
+Base.@kwdef mutable struct ReliabilityBranching <: VariableBranchingRule
+    min_samples::Int = 8
+    max_sb_calls::Int = 100
+    look_ahead::Int = 10
+    n_sb_calls::Int = 0
+    side_effect::Bool = true
+end
+
+function find_branching_var(
+    rule::ReliabilityBranching,
+    node::Node,
+    pool::NodePool,
+)::Variable
+    nfrac = length(node.fractional_variables)
+    pseudocost_scores = [
+        _pseudocost_score(
+            node,
+            pool,
+            node.fractional_variables[j],
+            node.fractional_values[j],
+        ) for j = 1:nfrac
+    ]
+    σ = sortperm(pseudocost_scores, rev = true)
+    sorted_vars = node.fractional_variables[σ]
+    _strong_branch_start(node)
+    no_improv_count, n_sb_calls = 0, 0
+    max_score, max_var = pseudocost_scores[σ[1]], sorted_vars[1]
+    for (i, var) in enumerate(sorted_vars)
+        use_strong_branch = true
+        if n_sb_calls >= rule.max_sb_calls
+            use_strong_branch = false
+        else
+            if var in keys(pool.var_history)
+                varhist = pool.var_history[var]
+                hlength = min(
+                    length(varhist.obj_ratio_up),
+                    length(varhist.obj_ratio_down),
+                )
+                if hlength >= rule.min_samples
+                    use_strong_branch = false
+                end
+            end
+        end
+        if use_strong_branch
+            n_sb_calls += 1
+            rule.n_sb_calls += 1
+            score = _strong_branch_score(
+                node = node,
+                pool = pool,
+                var = var,
+                x = node.fractional_values[σ[i]],
+                side_effect = rule.side_effect,
+            )
+        else
+            score = pseudocost_scores[σ[i]]
+        end
+        if score > max_score
+            max_score, max_var = score, var
+            no_improv_count = 0
+        else
+            no_improv_count += 1
+        end
+        no_improv_count <= rule.look_ahead || break
+    end
+    _strong_branch_end(node)
+    return max_var
+end

src/bb/varbranch/strong.jl

@@ -0,0 +1,93 @@
+#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
+#  Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
+#  Released under the modified BSD license. See COPYING.md for more details.
+
+using Random
+
+"""
+    StrongBranching(look_ahead::Int, max_calls::Int)
+
+Branching strategy that selects a subset of fractional variables
+as candidates (according to pseudocosts) the solves two linear
+programming problems for each candidate.
+
+
+"""
+Base.@kwdef struct StrongBranching <: VariableBranchingRule
+    look_ahead::Int = 10
+    max_calls::Int = 100
+    side_effect::Bool = true
+end
+
+function find_branching_var(rule::StrongBranching, node::Node, pool::NodePool)::Variable
+    nfrac = length(node.fractional_variables)
+    pseudocost_scores = [
+        _pseudocost_score(
+            node,
+            pool,
+            node.fractional_variables[j],
+            node.fractional_values[j],
+        ) for j = 1:nfrac
+    ]
+    σ = sortperm(pseudocost_scores, rev = true)
+    sorted_vars = node.fractional_variables[σ]
+    _strong_branch_start(node)
+    no_improv_count, call_count = 0, 0
+    max_score, max_var = pseudocost_scores[σ[1]], sorted_vars[1]
+    for (i, var) in enumerate(sorted_vars)
+        call_count += 1
+        score = _strong_branch_score(
+            node = node,
+            pool = pool,
+            var = var,
+            x = node.fractional_values[σ[i]],
+            side_effect = rule.side_effect,
+        )
+        if score > max_score
+            max_score, max_var = score, var
+            no_improv_count = 0
+        else
+            no_improv_count += 1
+        end
+        no_improv_count <= rule.look_ahead || break
+        call_count <= rule.max_calls || break
+    end
+    _strong_branch_end(node)
+    return max_var
+end
+
+function _strong_branch_score(;
+    node::Node,
+    pool::NodePool,
+    var::Variable,
+    x::Float64,
+    side_effect::Bool,
+)::Tuple{Float64,Int}
+    obj_up, obj_down = 0, 0
+    try
+        obj_up, obj_down = probe(node.mip, var)
+    catch
+        @warn "strong branch error" var=var
+    end
+    obj_change_up = obj_up - node.obj
+    obj_change_down = obj_down - node.obj
+    if side_effect
+        _update_var_history(
+            pool = pool,
+            var = var,
+            x = x,
+            obj_change_down = obj_change_down,
+            obj_change_up = obj_change_up,
+        )
+    end
+    return (obj_change_up * obj_change_down, var.index)
+end
+
+function _strong_branch_start(node::Node)
+    set_bounds!(node.mip, node.branch_variables, node.branch_values, node.branch_values)
+end
+
+function _strong_branch_end(node::Node)
+    nbranch = length(node.branch_variables)
+    set_bounds!(node.mip, node.branch_variables, zeros(nbranch), ones(nbranch))
+end

test/bb/lp_test.jl

@@ -0,0 +1,115 @@
+#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
+#  Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
+#  Released under the modified BSD license. See COPYING.md for more details.
+
+using Clp
+using JuMP
+using Test
+using MIPLearn.BB
+
+basepath = @__DIR__
+
+function runtests(optimizer_name, optimizer; large = true)
+    @testset "Solve ($optimizer_name)" begin
+        @testset "interface" begin
+            filename = "$basepath/../fixtures/danoint.mps.gz"
+
+            mip = BB.init(optimizer)
+            BB.read!(mip, filename)
+
+            @test mip.sense == 1.0
+            @test length(mip.binary_variables) == 56
+
+            status, obj = BB.solve_relaxation!(mip)
+            @test status == :Optimal
+            @test round(obj, digits = 6) == 62.637280
+
+            @test BB.name(mip, mip.binary_variables[1]) == "xab"
+            @test BB.name(mip, mip.binary_variables[2]) == "xac"
+            @test BB.name(mip, mip.binary_variables[3]) == "xad"
+
+            vals = BB.values(mip, mip.binary_variables)
+            @test round(vals[1], digits = 6) == 0.046933
+            @test round(vals[2], digits = 6) == 0.000841
+            @test round(vals[3], digits = 6) == 0.248696
+
+            # Probe (up and down are feasible)
+            probe_up, probe_down = BB.probe(mip, mip.binary_variables[1])
+            @test round(probe_down, digits = 6) == 62.690000
+            @test round(probe_up, digits = 6) == 62.714100
+
+            # Fix one variable to zero
+            BB.set_bounds!(mip, mip.binary_variables[1:1], [0.0], [0.0])
+            status, obj = BB.solve_relaxation!(mip)
+            @test status == :Optimal
+            @test round(obj, digits = 6) == 62.690000
+
+            # Fix one variable to one and another variable variable to zero
+            BB.set_bounds!(mip, mip.binary_variables[1:2], [1.0, 0.0], [1.0, 0.0])
+            status, obj = BB.solve_relaxation!(mip)
+            @test status == :Optimal
+            @test round(obj, digits = 6) == 62.714777
+
+            # Probe (up is infeasible, down is feasible)
+            BB.set_bounds!(
+                mip,
+                mip.binary_variables[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.binary_variables[3])
+            @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.binary_variables)
+            BB.set_bounds!(mip, mip.binary_variables, ones(N), ones(N))
+            status, obj = BB.solve_relaxation!(mip)
+            @test status == :Infeasible
+            @test obj == Inf
+
+            # Restore original problem
+            N = length(mip.binary_variables)
+            BB.set_bounds!(mip, mip.binary_variables, zeros(N), ones(N))
+            status, obj = BB.solve_relaxation!(mip)
+            @test status == :Optimal
+            @test round(obj, digits = 6) == 62.637280
+        end
+
+        @testset "varbranch" begin
+            branch_rules = [
+                BB.RandomBranching(),
+                BB.FirstInfeasibleBranching(),
+                BB.LeastInfeasibleBranching(),
+                BB.MostInfeasibleBranching(),
+                BB.PseudocostBranching(),
+                BB.StrongBranching(),
+                BB.ReliabilityBranching(),
+                BB.HybridBranching(),
+            ]
+            for branch_rule in branch_rules
+                filename = "$basepath/../fixtures/vpm2.mps.gz"
+                mip = BB.init(optimizer)
+                BB.read!(mip, filename)
+                @info optimizer_name, branch_rule
+                @time BB.solve!(
+                    mip,
+                    initial_primal_bound = 13.75,
+                    print_interval = 10,
+                    node_limit = 100,
+                    branch_rule = branch_rule,
+                )
+            end
+        end
+    end
+end
+
+@testset "BB" begin
+    @time runtests("Clp", Clp.Optimizer)
+    if is_gurobi_available
+        using Gurobi
+        @time runtests("Gurobi", Gurobi.Optimizer)
+    end
+end

test/runtests.jl

@@ -6,6 +6,7 @@ using Test
 using MIPLearn
 
 MIPLearn.setup_logger()
+const is_gurobi_available = ("GUROBI_HOME" in keys(ENV))
 
 @testset "MIPLearn" begin
     include("fixtures/knapsack.jl")
@@ -15,4 +16,5 @@ MIPLearn.setup_logger()
     include("solvers/learning_solver_test.jl")
     # include("utils/benchmark_test.jl")
     include("utils/parse_test.jl")
+    include("bb/lp_test.jl")
 end

test/solvers/jump_solver_test.jl

@@ -8,7 +8,6 @@ using MIPLearn
 using PyCall
 using Test
 
-const is_gurobi_available = ("GUROBI_HOME" in keys(ENV))
 if is_gurobi_available
     using Gurobi
 end