Initial columns callback (#656)

* wip

* initial columns callback

* update BlockDecomposition dep

Project.toml

@@ -18,7 +18,7 @@ TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 
 [compat]
-BlockDecomposition = "1.7"
+BlockDecomposition = "1.8"
 DataStructures = "0.17, 0.18"
 DynamicSparseArrays = "0.5.3"
 MathOptInterface = "0.10, 1"

docs/make.jl

@@ -3,14 +3,17 @@ using Documenter, Coluna, Literate, BlockDecomposition
 TUTORIAL_GAP = joinpath(@__DIR__, "src", "start", "start.jl")
 TUTORIAL_CUTS = joinpath(@__DIR__, "src", "start", "cuts.jl")
 TUTORIAL_PRICING = joinpath(@__DIR__, "src", "start", "pricing.jl")
+TUTORIAL_INITCOLS = joinpath(@__DIR__, "src", "start", "initial_columns.jl")
 
 OUTPUT_GAP = joinpath(@__DIR__, "src", "start")
 OUTPUT_CUTS = joinpath(@__DIR__, "src", "start")
 OUTPUT_PRICING = joinpath(@__DIR__, "src", "start")
+OUTPUT_INITCOLS = joinpath(@__DIR__, "src", "start")
 
 Literate.markdown(TUTORIAL_GAP, OUTPUT_GAP, documenter=true)
 Literate.markdown(TUTORIAL_CUTS, OUTPUT_CUTS, documenter=true)
 Literate.markdown(TUTORIAL_PRICING, OUTPUT_PRICING, documenter=true)
+Literate.markdown(TUTORIAL_INITCOLS, OUTPUT_INITCOLS, documenter=true)
 
 makedocs(
     modules = [Coluna, BlockDecomposition],
@@ -28,6 +31,7 @@ makedocs(
             "Column generation" => joinpath("start", "start.md"),
             "Valid inequalities" => joinpath("start", "cuts.md"),
             "Pricing callback" => joinpath("start", "pricing.md"),
+            "Initial columns callback" => joinpath("start", "initial_columns.md")
         ],
         "Manual" => Any[
             "Decomposition" => joinpath("man", "decomposition.md"),

docs/src/start/initial_columns.jl

@@ -0,0 +1,83 @@
+# # Initial columns 
+
+# The initial columns callback let you provide initial columns associated to each problem
+# ahead the optimization.
+# This callback is useful when you have an efficient heuristic that finds feasible solutions
+# to the problem. You can then extract columns from the solutions and give them to Coluna
+# through the callback.
+# You have to make sure the columns you provide are feasible because Coluna won't check their 
+# feasibility.
+# The cost of the columns will be computed using the perennial cost of subproblem variables.
+
+# Let us see an example with the following generalized assignment problem :
+
+M = 1:3;
+J = 1:5;
+c = [1 1 1 1 1; 1.2 1.2 1.1 1.1 1; 1.3 1.3 1.1 1.2 1.4];
+Q = [3, 2, 3];
+
+# with the following Coluna configuration
+
+using JuMP, GLPK, BlockDecomposition, Coluna;
+
+coluna = optimizer_with_attributes(
+    Coluna.Optimizer,
+    "params" => Coluna.Params(
+        solver = Coluna.Algorithm.TreeSearchAlgorithm() # default branch-cut-and-price
+    ),
+    "default_optimizer" => GLPK.Optimizer # GLPK for the master & the subproblems
+);
+
+# for which the JuMP model takes the form:
+
+@axis(M_axis, M);
+model = BlockModel(coluna); 
+
+@variable(model, x[m in M_axis, j in J], Bin);
+@constraint(model, cov[j in J], sum(x[m, j] for m in M_axis) >= 1);
+@constraint(model, knp[m in M_axis], sum(x[m, j] for j in J) <= Q[m]);
+@objective(model, Min, sum(c[m, j] * x[m, j] for m in M_axis, j in J));
+
+@dantzig_wolfe_decomposition(model, decomposition, M_axis)
+
+subproblems = getsubproblems(decomposition)
+specify!.(subproblems, lower_multiplicity = 0, upper_multiplicity = 1)
+
+
+# Let's consider that the following assignement patterns are good candidates:
+
+machine1 = [[1,2,4], [1,3,4], [2,3,4], [2,3,5]];
+machine2 = [[1,2], [1,5], [2,5], [3,4]];
+machine3 = [[1,2,3], [1,3,4], [1,3,5], [2,3,4]];
+
+initial_columns = [machine1, machine2, machine3];
+
+# We can write the initial columns callback:
+
+function initial_columns_callback(cbdata)
+    ## Retrieve the index of the subproblem (it will be one of the values in M_axis)
+    spid = BlockDecomposition.callback_spid(cbdata, model)
+    println("initial columns callback $spid")
+
+    ## Retrieve assignment patterns of a given machine
+    for col in initial_columns[spid]
+        ## Create the column in the good representation
+        vars = [x[spid, j] for j in col]
+        vals = [1.0 for _ in col]
+
+        ## Submit the column
+        MOI.submit(model, BlockDecomposition.InitialColumn(cbdata), vars, vals)
+    end
+end
+
+# The initial columns callback is a function.
+# It takes as argument `cbdata` which is a data structure
+# that allows the user to interact with Coluna within the callback.
+
+# We provide the initial columns callback to Coluna through the following method:
+
+MOI.set(model, BlockDecomposition.InitialColumnsCallback(), initial_columns_callback)
+
+# You can then optimize:
+
+optimize!(model)

src/MOIcallbacks.jl

@@ -21,6 +21,13 @@ function MOI.set(model::Coluna.Optimizer, ::BD.PricingCallback, ::Nothing)
     return
 end
 
+function MOI.set(model::Coluna.Optimizer, ::BD.InitialColumnsCallback, callback_function::Function)
+    model.has_initialcol_cb = true
+    problem = model.inner
+    MathProg._register_initcols_callback!(problem, callback_function)
+    return
+end
+
 ############################################################################################
 #  Pricing Callback                                                                        #
 ############################################################################################
@@ -174,3 +181,31 @@ end
 
 MOI.supports(::Optimizer, ::MOI.UserCutCallback) = true
 MOI.supports(::Optimizer, ::MOI.LazyConstraintCallback) = true
+
+############################################################################################
+#  Initial columns Callback                                                                        #
+############################################################################################
+function _submit_initial_solution(env, cbdata, variables, values, custom_data)
+    @assert length(variables) == length(values)
+    form = cbdata.form
+    colunavarids = [_get_varid_of_origvar_in_form(env, form, v) for v in variables]
+    cost = sum(value * getperencost(form, varid) for (varid, value) in Iterators.zip(colunavarids, values))
+    return _submit_pricing_solution(env, cbdata, cost, variables, values, custom_data)
+end
+
+function MOI.submit(
+    model::Optimizer,
+    cb::BD.InitialColumn{MathProg.InitialColumnsCallbackData},
+    variables::Vector{MOI.VariableIndex},
+    values::Vector{Float64},
+    custom_data::Union{Nothing, BD.AbstractCustomData} = nothing
+)
+    return _submit_initial_solution(model.env, cb.callback_data, variables, values, custom_data)
+end
+
+function MOI.get(model::Optimizer, spid::BD.PricingSubproblemId{MathProg.InitialColumnsCallbackData})
+    callback_data = spid.callback_data
+    uid = getuid(callback_data.form)
+    axis_index_value = model.annotations.ann_per_form[uid].axis_index_value
+    return axis_index_value
+end

src/MOIwrapper.jl

@@ -48,6 +48,7 @@ mutable struct Optimizer <: MOI.AbstractOptimizer
     has_pricing_cb::Bool
     has_usercut_cb::Bool
     has_lazyconstraint_cb::Bool
+    has_initialcol_cb::Bool
 
     function Optimizer()
         model = new()
@@ -70,6 +71,7 @@ mutable struct Optimizer <: MOI.AbstractOptimizer
         model.has_pricing_cb = false
         model.has_usercut_cb = false
         model.has_lazyconstraint_cb = false
+        model.has_initialcol_cb = false
         return model
     end
 end
@@ -101,6 +103,7 @@ function MOI.empty!(model::Optimizer)
     model.has_pricing_cb = false
     model.has_usercut_cb = false
     model.has_lazyconstraint_cb = false
+    model.has_initialcol_cb = false
     return
 end
 
@@ -1076,6 +1079,9 @@ function MOI.get(model::Optimizer, ::MOI.ListOfModelAttributesSet)
     if model.has_pricing_cb
         push!(attributes, BD.PricingCallback())
     end
+    if model.has_initialcol_cb
+        push!(attributes, BD.InitialColumnsCallback())
+    end
     return attributes
 end
 

src/MathProg/formulation.jl

@@ -278,6 +278,16 @@ function _sol_repr_for_pool(primal_sol::PrimalSolution, ::DwSp)
     return var_ids, vals
 end
 
+function initialize_solution_pool!(form::Formulation{DwSp}, initial_columns_callback::Function)
+    master = getmaster(form)
+    cbdata = InitialColumnsCallbackData(form, PrimalSolution[])
+    initial_columns_callback(cbdata)
+    for sol in cbdata.primal_solutions
+        insert_column!(master, sol, "iMC")
+    end
+    return
+end
+
 ############################################################################################
 # Insertion of a column in the master
 ############################################################################################

src/MathProg/problem.jl

@@ -4,6 +4,7 @@ mutable struct Problem <: AbstractProblem
     original_formulation::Formulation
     re_formulation::Union{Nothing, Reformulation}
     default_optimizer_builder::Function
+    initial_columns_callback::Union{Nothing, Function}
 end
 
 """
@@ -15,7 +16,7 @@ function Problem(env)
     original_formulation = create_formulation!(env, Original())
     return Problem(
         nothing, nothing, original_formulation, nothing,
-        no_optimizer_builder
+        no_optimizer_builder, nothing
     )
 end
 
@@ -63,3 +64,8 @@ function get_optimization_target(p::Problem)
     end
     return p.re_formulation
 end
+
+function _register_initcols_callback!(problem::Problem, callback_function::Function)
+    problem.initial_columns_callback = callback_function
+    return
+end

src/MathProg/reformulation.jl

@@ -100,6 +100,23 @@ subproblem with id `spid`.
 """
 get_dw_pricing_sp_lb_constrid(r::Reformulation, spid::FormId) = r.dw_pricing_sp_lb[spid]
 
+############################################################################################
+# Initial columns callback
+############################################################################################
+struct InitialColumnsCallbackData
+    form::Formulation
+    primal_solutions::Vector{PrimalSolution}
+end
+
+# Method to initial the solution pools of the subproblems
+function initialize_solution_pools!(reform::Reformulation, initial_columns_callback::Function)
+    for (_, sp) in get_dw_pricing_sps(reform)
+        initialize_solution_pool!(sp, initial_columns_callback)
+    end
+    return
+end
+
+initialize_solution_pools!(::Reformulation, ::Nothing) = nothing # fallback
 
 # Following two functions are temporary, we must store a pointer to the vc
 # being represented by a representative vc

src/optimize.jl

@@ -39,6 +39,7 @@ function optimize!(env::Env, prob::MathProg.Problem, annotations::Annotations)
     ## Retrieve initial bounds on the objective given by the user
     init_pb = get_initial_primal_bound(prob)
     init_db = get_initial_dual_bound(prob)
+    init_cols = prob.initial_columns_callback
     _adjust_params(env.params, init_pb)
 
     # Apply decomposition
@@ -50,7 +51,7 @@ function optimize!(env::Env, prob::MathProg.Problem, annotations::Annotations)
     @logmsg LogLevel(-1) env.params
 
     TO.@timeit _to "Coluna" begin
-        outstate, algstate = optimize!(get_optimization_target(prob), env, init_pb, init_db)
+        outstate, algstate = optimize!(get_optimization_target(prob), env, init_pb, init_db, init_cols)
     end
 
     env.kpis.elapsed_optimization_time = elapsed_optim_time(env)
@@ -67,7 +68,8 @@ function optimize!(env::Env, prob::MathProg.Problem, annotations::Annotations)
 end
 
 function optimize!(
-    reform::MathProg.Reformulation, env::Env, initial_primal_bound, initial_dual_bound
+    reform::MathProg.Reformulation, env::Env, initial_primal_bound, initial_dual_bound,
+    initial_columns
 )
     master = getmaster(reform)
     initstate = OptimizationState(
@@ -79,6 +81,9 @@ function optimize!(
 
     algorithm = env.params.solver
 
+    # retrieve initial columns
+    MathProg.initialize_solution_pools!(reform, initial_columns)
+
     # initialize all the units used by the algorithm and its child algorithms
     Algorithm.initialize_storage_units!(reform, algorithm)
 
@@ -124,7 +129,7 @@ function optimize!(
 end
 
 function optimize!(
-    form::MathProg.Formulation, env::Env, initial_primal_bound, initial_dual_bound
+    form::MathProg.Formulation, env::Env, initial_primal_bound, initial_dual_bound, _
 )
     initstate = OptimizationState(
         form,
@@ -140,14 +145,14 @@ end
 """
 Fallback if no solver provided by the user.
 """
-function optimize!(::MathProg.Reformulation, ::Nothing, ::Real, ::Real)
+function optimize!(::MathProg.Reformulation, ::Nothing, ::Real, ::Real, _)
     error("""
         No solver to optimize the reformulation. You should provide a solver through Coluna parameters. 
         Please, check the starting guide of Coluna.
     """)
 end
 
-function optimize!(::MathProg.Formulation, ::Nothing, ::Real, ::Real)
+function optimize!(::MathProg.Formulation, ::Nothing, ::Real, ::Real, _)
     error("""
         No solver to optimize the formulation. You should provide a solver through Coluna parameters. 
         Please, check the starting guide of Coluna.