Merge branch 'master' into compathelper/new_version/2022-04-01-00-53-…

…18-789-4156876630

Project.toml

@@ -1,7 +1,7 @@
 name = "Coluna"
 uuid = "88b4ec78-b192-11e8-04aa-4d367dd96a64"
 authors = ["François Vanderbeck", "Guillaume Marques", "Vitor Nesello", "Ruslan Sadykov"]
-version = "0.4.0"
+version = "0.4.1"
 
 [deps]
 BlockDecomposition = "6cde8614-403a-11e9-12f1-c10d0f0caca0"
@@ -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/cuts.jl

@@ -5,12 +5,12 @@
 
 # Consider the following instance:
 
-J = 1:15
-M = 1:6
-c = [10.13 15.6 15.54 13.41 17.08 17.78 16.36; 19.58 16.83 10.75 15.8 14.89 11.91 15.01; 14.23 17.36 16.05 14.49 18.96 13.58 12.52; 16.47 16.38 18.14 15.46 11.64 10.56 12.65; 17.87 18.25 13.12 19.16 16.33 11.28 17.71; 11.09 16.76 15.5 12.08 13.06 19.03 14.89; 15.19 13.86 16.08 19.47 15.79 13.05 11.06; 10.79 18.96 16.11 19.78 15.57 12.68 14.76; 12.03 19.03 16.01 14.46 12.77 19.57 19.83; 14.48 11.75 16.97 19.95 18.32 13.43 12.7; 10.01 10.8 16.22 16.69 17.26 11.73 13.55; 17.84 11.91 10.47 10.91 17.87 14.63 11.52; 17.93 11.24 12.54 10.48 18.69 19.58 17.67; 16.05 15.83 10.31 19.81 18.72 10.2 14.95; 17.1 15.53 12.34 18.34 13.64 12.3 12.79];
-w = [5, 4, 5, 6, 8, 9, 5, 8, 10, 7, 7, 7, 7, 7, 8];
-Q = [ 25,  24,  31,  28,  24,  30];
-f = [105, 103, 109, 112, 100, 104];
+J = 1:10
+M = 1:5
+c = [10.13 15.6 15.54 13.41 17.08;19.58 16.83 10.75 15.8 14.89;14.23 17.36 16.05 14.49 18.96;16.47 16.38 18.14 15.46 11.64;17.87 18.25 13.12 19.16 16.33;11.09 16.76 15.5 12.08 13.06;15.19 13.86 16.08 19.47 15.79;10.79 18.96 16.11 19.78 15.55;12.03 19.03 16.01 14.46 12.77;14.48 11.75 16.97 19.95 18.32];
+w = [5, 4, 5, 6, 8, 9, 5, 8, 10, 7];
+Q = [25,  24,  31,  28,  24];
+f = [105, 103, 109, 112, 100];
 
 # We define the dependencies:
 

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)

docs/src/start/pricing.jl

@@ -13,11 +13,11 @@ const MOI = MathOptInterface;
 
 # Let us see an example with the following generalized assignment problem :
 
-M = 1:4;
-J = 1:30;
-c = [12.7 22.5 8.9 20.8 13.6 12.4 24.8 19.1 11.5 17.4 24.7 6.8 21.7 14.3 10.5 15.2 14.3 12.6 9.2 20.8 11.7 17.3 9.2 20.3 11.4 6.2 13.8 10.0 20.9 20.6;  19.1 24.8 24.4 23.6 16.1 20.6 15.0 9.5 7.9 11.3 22.6 8.0 21.5 14.7 23.2 19.7 19.5 7.2 6.4 23.2 8.1 13.6 24.6 15.6 22.3 8.8 19.1 18.4 22.9 8.0;  18.6 14.1 22.7 9.9 24.2 24.5 20.8 12.9 17.7 11.9 18.7 10.1 9.1 8.9 7.7 16.6 8.3 15.9 24.3 18.6 21.1 7.5 16.8 20.9 8.9 15.2 15.7 12.7 20.8 10.4;  13.1 16.2 16.8 16.7 9.0 16.9 17.9 12.1 17.5 22.0 19.9 14.6 18.2 19.6 24.2 12.9 11.3 7.5 6.5 11.3 7.8 13.8 20.7 16.8 23.6 19.1 16.8 19.3 12.5 11.0];
-w = [61 70 57 82 51 74 98 64 86 80 69 79 60 76 78 71 50 99 92 83 53 91 68 61 63 97 91 77 68 80; 50 57 61 83 81 79 63 99 82 59 83 91 59 99 91 75 66 100 69 60 87 98 78 62 90 89 67 87 65 100; 91 81 66 63 59 81 87 90 65 55 57 68 92 91 86 74 80 89 95 57 55 96 77 60 55 57 56 67 81 52;  62 79 73 60 75 66 68 99 69 60 56 100 67 68 54 66 50 56 70 56 72 62 85 70 100 57 96 69 65 50];
-Q = [1020 1460 1530 1190];
+M = 1:3;
+J = 1:15;
+c = [12.7 22.5 8.9 20.8 13.6 12.4 24.8 19.1 11.5 17.4 24.7 6.8 21.7 14.3 10.5; 19.1 24.8 24.4 23.6 16.1 20.6 15.0 9.5 7.9 11.3 22.6 8.0 21.5 14.7 23.2;  18.6 14.1 22.7 9.9 24.2 24.5 20.8 12.9 17.7 11.9 18.7 10.1 9.1 8.9 7.7; 13.1 16.2 16.8 16.7 9.0 16.9 17.9 12.1 17.5 22.0 19.9 14.6 18.2 19.6 24.2];
+w = [61 70 57 82 51 74 98 64 86 80 69 79 60 76 78; 50 57 61 83 81 79 63 99 82 59 83 91 59 99 91;91 81 66 63 59 81 87 90 65 55 57 68 92 91 86; 62 79 73 60 75 66 68 99 69 60 56 100 67 68 54];
+Q = [1020 1460 1530];
 
 # with the following Coluna configuration
 
@@ -68,7 +68,9 @@ end
 function my_pricing_callback(cbdata)
     ## Retrieve the index of the subproblem (it will be one of the values in M_axis)
     cur_machine = BD.callback_spid(cbdata, model)
-    println("Pricing callback for machine $(cur_machine).")
+
+    ## Uncomment to see that the pricing callback is called.
+    ## println("Pricing callback for machine $(cur_machine).")
 
     ## Retrieve reduced costs of subproblem variables
     red_costs = [BD.callback_reduced_cost(cbdata, x[cur_machine, j]) for j in J]

docs/src/start/start.jl

@@ -26,11 +26,11 @@
 
 # Let us consider the following instance. 
 
-M = 1:4;
-J = 1:30;
-c = [12.7 22.5 8.9 20.8 13.6 12.4 24.8 19.1 11.5 17.4 24.7 6.8 21.7 14.3 10.5 15.2 14.3 12.6 9.2 20.8 11.7 17.3 9.2 20.3 11.4 6.2 13.8 10.0 20.9 20.6;  19.1 24.8 24.4 23.6 16.1 20.6 15.0 9.5 7.9 11.3 22.6 8.0 21.5 14.7 23.2 19.7 19.5 7.2 6.4 23.2 8.1 13.6 24.6 15.6 22.3 8.8 19.1 18.4 22.9 8.0;  18.6 14.1 22.7 9.9 24.2 24.5 20.8 12.9 17.7 11.9 18.7 10.1 9.1 8.9 7.7 16.6 8.3 15.9 24.3 18.6 21.1 7.5 16.8 20.9 8.9 15.2 15.7 12.7 20.8 10.4;  13.1 16.2 16.8 16.7 9.0 16.9 17.9 12.1 17.5 22.0 19.9 14.6 18.2 19.6 24.2 12.9 11.3 7.5 6.5 11.3 7.8 13.8 20.7 16.8 23.6 19.1 16.8 19.3 12.5 11.0];
-w = [61 70 57 82 51 74 98 64 86 80 69 79 60 76 78 71 50 99 92 83 53 91 68 61 63 97 91 77 68 80; 50 57 61 83 81 79 63 99 82 59 83 91 59 99 91 75 66 100 69 60 87 98 78 62 90 89 67 87 65 100; 91 81 66 63 59 81 87 90 65 55 57 68 92 91 86 74 80 89 95 57 55 96 77 60 55 57 56 67 81 52;  62 79 73 60 75 66 68 99 69 60 56 100 67 68 54 66 50 56 70 56 72 62 85 70 100 57 96 69 65 50];
-Q = [1020 1460 1530 1190];
+M = 1:3;
+J = 1:15;
+c = [12.7 22.5 8.9 20.8 13.6 12.4 24.8 19.1 11.5 17.4 24.7 6.8 21.7 14.3 10.5; 19.1 24.8 24.4 23.6 16.1 20.6 15.0 9.5 7.9 11.3 22.6 8.0 21.5 14.7 23.2;  18.6 14.1 22.7 9.9 24.2 24.5 20.8 12.9 17.7 11.9 18.7 10.1 9.1 8.9 7.7; 13.1 16.2 16.8 16.7 9.0 16.9 17.9 12.1 17.5 22.0 19.9 14.6 18.2 19.6 24.2];
+w = [61 70 57 82 51 74 98 64 86 80 69 79 60 76 78; 50 57 61 83 81 79 63 99 82 59 83 91 59 99 91;91 81 66 63 59 81 87 90 65 55 57 68 92 91 86; 62 79 73 60 75 66 68 99 69 60 56 100 67 68 54];
+Q = [1020 1460 1530];
 
 # We write the model with [JuMP](https://github.com/jump-dev/JuMP.jl), a domain-specific modeling
 # language for mathematical optimization embedded in Julia. We optimize with GLPK.

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/MOIinterface.jl

@@ -52,14 +52,6 @@ function update_cost_in_optimizer!(form::Formulation, optimizer::MoiOptimizer, v
     return
 end
 
-function update_obj_const_in_optimizer!(form::Formulation, optimizer::MoiOptimizer)
-    MOI.modify(
-        getinner(optimizer), MoiObjective(), 
-        MOI.ScalarConstantChange{Float64}(getobjconst(form))
-    )
-    return
-end
-
 function update_constr_member_in_optimizer!(
     optimizer::MoiOptimizer, c::Constraint, v::Variable, coeff::Float64
 )

src/MathProg/buffer.jl

@@ -54,7 +54,6 @@ Changes must be always buffered.
 """
 mutable struct FormulationBuffer{Vi,V,Ci,C}
     changed_obj_sense::Bool # sense of the objective function
-    changed_obj_const::Bool # constant in the objective function
     changed_cost::Set{Vi} # cost of a variable
     changed_bound::Set{Vi} # bound of a variable
     changed_var_kind::Set{Vi} # kind of a variable
@@ -65,14 +64,13 @@ mutable struct FormulationBuffer{Vi,V,Ci,C}
 end
 
 FormulationBuffer{Vi,V,Ci,C}() where {Vi,V,Ci,C} = FormulationBuffer(
-    false, false, Set{Vi}(), Set{Vi}(), Set{Vi}(), Set{Ci}(),
+    false, Set{Vi}(), Set{Vi}(), Set{Vi}(), Set{Ci}(),
     VarConstrBuffer{Vi, V}(), VarConstrBuffer{Ci, C}(), 
     Dict{Pair{Ci,Vi},Float64}()
 )
 
 function empty!(buffer::FormulationBuffer{Vi,V,Ci,C}) where {Vi,V,Ci,C}
     buffer.changed_obj_sense = false
-    buffer.changed_obj_const = false
     buffer.changed_cost = Set{Vi}()
     buffer.changed_bound = Set{Vi}()
     buffer.changed_var_kind = Set{Vi}()

src/MathProg/formulation.jl

@@ -92,7 +92,6 @@ getobjconst(form::Formulation) = form.manager.objective_constant
 "Sets objective constant of the formulation."
 function setobjconst!(form::Formulation, val::Float64)
     form.manager.objective_constant = val
-    form.buffer.changed_obj_const = true
     return
 end
 
@@ -278,6 +277,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/optimizerwrappers.jl

@@ -77,12 +77,6 @@ function sync_solver!(optimizer::MoiOptimizer, f::Formulation)
         update_cost_in_optimizer!(f, optimizer, getvar(f, id))
     end
 
-    # Update objective constant
-    if buffer.changed_obj_const
-        update_obj_const_in_optimizer!(f, optimizer)
-        buffer.changed_obj_const = false
-    end
-
     # Update objective sense
     if buffer.changed_obj_sense
         set_obj_sense!(optimizer, getobjsense(f))

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