Feature ad-hoc pricing solver

Manifest.toml

@@ -23,7 +23,7 @@ version = "0.5.4"
 
 [[BlockDecomposition]]
 deps = ["DataStructures", "JuMP", "MathOptInterface", "Test"]
-git-tree-sha1 = "c1b57e5ed4b891b2fe3326dba7f558b2e4893654"
+git-tree-sha1 = "caed35511b5d1dd20fbb53a7bac5bc8f63530973"
 repo-rev = "master"
 repo-url = "https://github.com/atoptima/BlockDecomposition.jl.git"
 uuid = "6cde8614-403a-11e9-12f1-c10d0f0caca0"

src/Coluna.jl

@@ -3,7 +3,6 @@ module Coluna
 import MathOptInterface
 import MathOptInterface.Utilities
 import DataStructures
-import GLPK
 import JuMP
 import BlockDecomposition
 import Distributed

src/MOIwrapper.jl

@@ -23,11 +23,13 @@ end
 
 setinnerprob!(o::Optimizer, prob::Problem) = o.inner = prob
 
-function Optimizer(;master_factory =
-        JuMP.with_optimizer(GLPK.Optimizer), pricing_factory =
-        JuMP.with_optimizer(GLPK.Optimizer), benders_sep_factory =
-        JuMP.with_optimizer(GLPK.Optimizer), params = Params())
-    prob = Problem(master_factory, pricing_factory, benders_sep_factory)
+function Optimizer(;default_optimizer = nothing,
+                   params = Params())
+    b = no_optimizer_builder
+    if default_optimizer != nothing
+        b = ()->MoiOptimizer(default_optimizer())
+    end
+    prob = Problem(b)
     return Optimizer(
         prob, MOIU.IndexMap(), params, Annotations(),
         Dict{MOI.VariableIndex,Id{Variable}}(), OptimizationResult{MinSense}()
@@ -257,4 +259,4 @@ function MOI.get(optimizer::Optimizer, object::MOI.TerminationStatus)
     result = optimizer.result
     isfeasible(result) && return MOI.OPTIMAL
     !isfeasible(result) && return MOI.INFEASIBLE
-end
+end

src/decomposition.jl

@@ -63,7 +63,8 @@ function build_dw_master!(prob::Problem,
                           reformulation::Reformulation,
                           master_form::Formulation,
                           vars_in_form::VarDict,
-                          constrs_in_form::ConstrDict)
+                          constrs_in_form::ConstrDict,
+                          opt_builder::Function)
 
     orig_form = get_original_formulation(prob)
     reformulation.dw_pricing_sp_lb = Dict{FormId, Id}()
@@ -136,21 +137,24 @@ function build_dw_master!(prob::Problem,
     # add artificial var 
     initialize_artificial_variables(master_form, constrs_in_form)
     initialize_local_art_vars(master_form, convexity_constrs)
+    initialize_optimizer!(master_form, opt_builder)
     return
 end
 
 function build_dw_pricing_sp!(prob::Problem,
                               annotation_id::Int,
                               sp_form::Formulation,
                               vars_in_form::VarDict,
-                              constrs_in_form::ConstrDict)
+                              constrs_in_form::ConstrDict,
+                              opt_builder::Function)
 
     orig_form = get_original_formulation(prob)
     master_form = sp_form.parent_formulation
     reformulation = master_form.parent_formulation
     ## Create Pure Pricing Sp Var & constr
     clone_in_formulation!(sp_form, orig_form, vars_in_form, PricingSpVar)
     clone_in_formulation!(sp_form, orig_form, constrs_in_form, PricingSpPureConstr)
+    initialize_optimizer!(sp_form, opt_builder)
     return
 end
 
@@ -188,9 +192,13 @@ function reformulate!(prob::Problem, annotations::Annotations,
     ann_sorted_by_uid = sort(collect(annotation_set), by = ann -> ann.unique_id)
 
     formulations = Dict{Int, Formulation}()
+    optimizer_builders = Dict{Int, Function}()
     master_unique_id = -1
 
     for annotation in ann_sorted_by_uid
+        if BD.getoptimizerbuilder(annotation) != nothing
+            optimizer_builders[BD.getid(annotation)] = BD.getoptimizerbuilder(annotation)
+        end
         if BD.getformulation(annotation) == BD.Master
             master_unique_id = BD.getid(annotation)
             formulations[BD.getid(annotation)] = master_form
@@ -213,18 +221,20 @@ function reformulate!(prob::Problem, annotations::Annotations,
             vars, constrs = find_vcs_in_block(
                 BD.getid(annotation), vars_per_block, constrs_per_block
             )
+            opt_builder = get(optimizer_builders, BD.getid(annotation), prob.default_optimizer_builder)
             build_dw_pricing_sp!(prob, BD.getid(annotation),
                                  formulations[BD.getid(annotation)],
-                                 vars, constrs)
+                                 vars, constrs, opt_builder)
         end
     end
 
     # Build Master
     vars, constrs = find_vcs_in_block(
         master_unique_id, vars_per_block, constrs_per_block
     )
+    opt_builder = get(optimizer_builders, master_unique_id, prob.default_optimizer_builder)
     build_dw_master!(prob, master_unique_id, reformulation,
-                     master_form, vars, constrs)
+                     master_form, vars, constrs, opt_builder)
 
     @debug "\e[1;34m Master formulation \e[00m" master_form
     for sp_form in reformulation.dw_pricing_subprs

src/formulation.jl

@@ -351,8 +351,14 @@ function optimize!(form::Formulation)
     return res
 end
 
-function initialize_optimizer(form::Formulation, factory::JuMP.OptimizerFactory)
-    form.optimizer = create_optimizer(factory, form.obj_sense)
+function initialize_optimizer!(form::Formulation, builder::Union{Function})
+    form.optimizer = builder()
+    if form.optimizer isa MoiOptimizer
+        f = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{Float64}[], 0.0)
+        MOI.set(form.optimizer.inner, MoiObjective(), f)
+        set_obj_sense!(form.optimizer, getobjsense(form))
+    end
+    return
 end
 
 function _show_obj_fun(io::IO, f::Formulation)

src/optimizerwrappers.jl

@@ -5,6 +5,8 @@ Wrapper to indicate that no optimizer is assigned to a `Formulation`
 """
 struct NoOptimizer <: AbstractOptimizer end
 
+no_optimizer_builder(args...) = NoOptimizer()
+
 """
     UserOptimizer <: AbstractOptimizer
 
@@ -14,13 +16,9 @@ mutable struct UserOptimizer <: AbstractOptimizer
     optimize_function::Function
 end
 
-function create_optimizer(optimize_function::Function)
-    return UserOptimizer(optimize_function)
-end
-
 function optimize!(form::Formulation, optimizer::UserOptimizer)
-    println("Calling user-defined optimization function.")
-    return OptimizationResult{getobjsense(form)}()
+    @logmsg LogLevel(-2) "Calling user-defined optimization function."
+    return optimizer.optimize_function(form)
 end
 
 """
@@ -35,16 +33,6 @@ end
 
 getinner(optimizer::MoiOptimizer) = optimizer.inner
 
-function create_optimizer(factory::JuMP.OptimizerFactory,
-                          sense::Type{<:AbstractObjSense})
-    moi_optimizer = factory()
-    f = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{Float64}[], 0.0)
-    MOI.set(moi_optimizer, MoiObjective(),f)
-    optimizer = MoiOptimizer(moi_optimizer)
-    set_obj_sense!(optimizer, sense)
-    return optimizer
-end
-
 function retrieve_result(form::Formulation, optimizer::MoiOptimizer)
     result = OptimizationResult{getobjsense(form)}()
     if MOI.get(form.optimizer.inner, MOI.ResultCount()) >= 1
@@ -147,7 +135,3 @@ end
 optimize!(::S) where {S<:AbstractOptimizer} = error(
     string("Function `optimize!` is not defined for object of type ", S)
 )
-
-create_optimizer(::S) where {S<:AbstractOptimizer} = error(
-    string("Function `create_optimizer` is not defined for object of type ", S)
-)

src/problem.jl

@@ -20,35 +20,22 @@ mutable struct Problem <: AbstractProblem
     original_formulation::Union{Nothing, Formulation}
     re_formulation::Union{Nothing, Reformulation}
     form_counter::Counter # 0 is for original form
-    master_factory::Union{Nothing, JuMP.OptimizerFactory}
-    pricing_factory::Union{Nothing, JuMP.OptimizerFactory}
-    benders_sep_factory::Union{Nothing, JuMP.OptimizerFactory}
+    default_optimizer_builder::Function
 end
 
 """
-    Problem(params::Params, master_factory, pricing_factory, benders_sep_factory)
+    Problem(b::Function)
 
 Constructs an empty `Problem`.
 """
-function Problem(master_factory, pricing_factory, benders_sep_factory)
-    return Problem(
-        nothing, nothing, Counter(-1),
-        master_factory, pricing_factory, benders_sep_factory
-    )
-end
+Problem(b::Function) = Problem(nothing, nothing, Counter(-1), b)
 
 set_original_formulation!(m::Problem, of::Formulation) = m.original_formulation = of
 set_re_formulation!(m::Problem, r::Reformulation) = m.re_formulation = r
 
 get_original_formulation(m::Problem) = m.original_formulation
 get_re_formulation(m::Problem) = m.re_formulation
 
-function initialize_optimizer(prob::Problem)
-    initialize_optimizer(
-        prob.re_formulation, prob.master_factory, prob.pricing_factory
-    )
-end
-
 function _welcome_message()
     welcome = """
     Coluna
@@ -63,7 +50,6 @@ function coluna_initialization(prob::Problem, annotations::Annotations,
     _set_global_params(params)
     reformulate!(prob, annotations, params.global_strategy)
     relax_integrality!(prob.re_formulation.master)
-    initialize_optimizer(prob)
     @info "Coluna initialized."
 end
 

src/reformulation.jl

@@ -30,15 +30,6 @@ getmaster(r::Reformulation) = r.master
 setmaster!(r::Reformulation, f) = r.master = f
 add_dw_pricing_sp!(r::Reformulation, f) = push!(r.dw_pricing_subprs, f)
 
-function initialize_optimizer(reformulation::Reformulation,
-                                  master_factory::JuMP.OptimizerFactory,
-                                  pricing_factory::JuMP.OptimizerFactory)
-    initialize_optimizer(reformulation.master, master_factory)
-    for problem in reformulation.dw_pricing_subprs
-        initialize_optimizer(problem, pricing_factory)
-    end
-end
-
 function optimize!(reformulation::Reformulation)
     opt_result = apply!(GlobalStrategy, reformulation)
     opt_result.primal_sols = [proj_cols_on_rep(

test/full_instances_tests.jl

@@ -2,9 +2,9 @@ function full_instances_tests()
     @testset "play gap" begin
         data = CLD.GeneralizedAssignment.data("play2.txt")
 
-        coluna = JuMP.with_optimizer(Coluna.Optimizer,# #params = params,
-            master_factory = with_optimizer(GLPK.Optimizer),
-            pricing_factory = with_optimizer(GLPK.Optimizer))
+        coluna = JuMP.with_optimizer(Coluna.Optimizer,
+            default_optimizer = with_optimizer(GLPK.Optimizer)
+        )
 
         problem, x, dec = CLD.GeneralizedAssignment.model(data, coluna)
 
@@ -17,9 +17,9 @@ function full_instances_tests()
     @testset "gap - JuMP/MOI modeling" begin
         data = CLD.GeneralizedAssignment.data("smallgap3.txt")
 
-        coluna = JuMP.with_optimizer(Coluna.Optimizer,# #params = params,
-            master_factory = with_optimizer(GLPK.Optimizer),
-            pricing_factory = with_optimizer(GLPK.Optimizer))
+        coluna = JuMP.with_optimizer(Coluna.Optimizer,
+            default_optimizer = with_optimizer(GLPK.Optimizer)
+        )
 
         problem, x, dec = CLD.GeneralizedAssignment.model(data, coluna)
 
@@ -32,9 +32,8 @@ function full_instances_tests()
     @testset "gap with penalties - pure master variables" begin
         data = CLD.GeneralizedAssignment.data("smallgap3.txt")
 
-        coluna = JuMP.with_optimizer(Coluna.Optimizer,# #params = params,
-            master_factory = with_optimizer(GLPK.Optimizer),
-            pricing_factory = with_optimizer(GLPK.Optimizer)
+        coluna = JuMP.with_optimizer(Coluna.Optimizer,
+            default_optimizer = with_optimizer(GLPK.Optimizer)
         )
 
         problem, x, y, dec = CLD.GeneralizedAssignment.model_with_penalties(data, coluna)
@@ -46,9 +45,8 @@ function full_instances_tests()
     @testset "gap with maximisation objective function" begin
         data = CLD.GeneralizedAssignment.data("smallgap3.txt")
 
-        coluna = JuMP.with_optimizer(Coluna.Optimizer,# #params = params,
-            master_factory = with_optimizer(GLPK.Optimizer),
-            pricing_factory = with_optimizer(GLPK.Optimizer)
+        coluna = JuMP.with_optimizer(Coluna.Optimizer,
+            default_optimizer = with_optimizer(GLPK.Optimizer)
         )
 
         problem, x, dec = CLD.GeneralizedAssignment.model_max(data, coluna)
@@ -60,9 +58,9 @@ function full_instances_tests()
     @testset "gap with infeasible subproblem" begin
     data = CLD.GeneralizedAssignment.data("root_infeas.txt")
 
-    coluna = JuMP.with_optimizer(Coluna.Optimizer,# #params = params,
-        master_factory = with_optimizer(GLPK.Optimizer),
-        pricing_factory = with_optimizer(GLPK.Optimizer))
+    coluna = JuMP.with_optimizer(Coluna.Optimizer,
+        default_optimizer = with_optimizer(GLPK.Optimizer)
+    )
 
     problem, x, dec = CLD.GeneralizedAssignment.model(data, coluna)
 
@@ -73,10 +71,9 @@ function full_instances_tests()
     # @testset "gap BIG instance" begin
     #     data = CLD.GeneralizedAssignment.data("gapC-5-100.txt")
 
-    #     coluna = JuMP.with_optimizer(Coluna.Optimizer,# #params = params,
-    #         master_factory = with_optimizer(GLPK.Optimizer),
-    #         pricing_factory = with_optimizer(GLPK.Optimizer)
-    #     )
+    # coluna = JuMP.with_optimizer(Coluna.Optimizer,
+    #     default_optimizer = with_optimizer(GLPK.Optimizer)
+    # )
 
     #     problem, x, dec = CLD.GeneralizedAssignment.model(data, coluna)
     #     JuMP.optimize!(problem)
@@ -91,9 +88,8 @@ function full_instances_tests()
     @testset "play gap" begin
         data = CLD.GeneralizedAssignment.data("play2.txt")
 
-        coluna = JuMP.with_optimizer(Coluna.Optimizer,# #params = params,
-            master_factory = with_optimizer(GLPK.Optimizer),
-            pricing_factory = with_optimizer(GLPK.Optimizer)
+        coluna = JuMP.with_optimizer(Coluna.Optimizer,
+            default_optimizer = with_optimizer(GLPK.Optimizer)
         )
 
         problem, x, dec = CLD.GeneralizedAssignment.model(data, coluna)

test/pricing_callback_tests.jl

@@ -0,0 +1,50 @@
+function mycallback(form::CL.Formulation)
+    vars = [v for (id,v) in filter(CL._active_explicit_, CL.getvars(form))]
+    constr = [c for (id,c) in filter(CL._active_explicit_, CL.getconstrs(form))][1]
+    matrix = CL.getcoefmatrix(form)
+    m = JuMP.Model(with_optimizer(GLPK.Optimizer))
+    @variable(m, CL.getcurlb(vars[i]) <= x[i=1:length(vars)] <= CL.getcurub(vars[i]), Int)
+    @objective(m, Min, sum(CL.getcurcost(vars[j]) * x[j] for j in 1:length(vars)))
+    @constraint(m, knp, 
+        sum(matrix[CL.getid(constr),CL.getid(vars[j])] * x[j]
+        for j in 1:length(vars)) <= CL.getcurrhs(constr)
+    )
+    optimize!(m)
+    result = CL.OptimizationResult{CL.MinSense}()
+    result.primal_bound = CL.PrimalBound{CL.MinSense}(JuMP.objective_value(m))
+    sol = CL.MembersVector{Float64}(CL.getvars(form))
+    for i in 1:length(x)
+        val = JuMP.value(x[i])
+        if val > 0.000001  || val < - 0.000001 # todo use a tolerance
+            sol[CL.getid(vars[i])] = val
+        end
+    end
+    push!(result.primal_sols, CL.PrimalSolution{CL.MinSense}(result.primal_bound, sol))
+    return result
+end
+
+build_sp_moi_optimizer() = CL.UserOptimizer(mycallback)
+build_master_moi_optimizer() = CL.MoiOptimizer(with_optimizer(GLPK.Optimizer)())
+
+function pricing_callback_tests()
+
+    @testset "GAP with ad-hoc pricing callback " begin
+        data = CLD.GeneralizedAssignment.data("play2.txt")
+
+        coluna = JuMP.with_optimizer(Coluna.Optimizer)
+
+        problem, x, dec = CLD.GeneralizedAssignment.model(data, coluna)
+
+        BD.assignsolver(dec, build_master_moi_optimizer)
+        BD.assignsolver(dec[1:2], build_sp_moi_optimizer)
+        @test BD.getoptimizerbuilder(dec) == build_master_moi_optimizer
+        @test BD.getoptimizerbuilder(dec[1]) == build_sp_moi_optimizer
+        @test BD.getoptimizerbuilder(dec[2]) == build_sp_moi_optimizer
+
+        JuMP.optimize!(problem)
+        @test abs(JuMP.objective_value(problem) - 75.0) <= 0.00001
+        @test MOI.get(problem.moi_backend.optimizer, MOI.TerminationStatus()) == MOI.OPTIMAL
+        @test CLD.GeneralizedAssignment.print_and_check_sol(data, problem, x)
+    end
+
+end