Implementation of column generation stages (#525)

* Implementation of column generation stages (for example, heuristic and exact stage)

* Update after conversation with Guillaume + stabilization correction

* Simplification for ColCutGenConquer

* Some more modifs due to Guillaume comments

docs/src/dev/todo.md

@@ -35,7 +35,7 @@ Base.diff
 Coluna.MathProg.update_ip_primal_bound!
 Coluna.Algorithm.OptimizationState
 Coluna.MathProg.update_lp_dual_bound!
-Coluna.MathProg.getoptimizer
+Coluna.MathProg.getmoioptimizer
 Coluna.MathProg.update_ip_dual_bound!
 Coluna.MathProg.ip_gap_closed
 Coluna.Algorithm.add_ip_primal_sol!

src/Algorithm/Algorithm.jl

@@ -34,8 +34,10 @@ include("interface.jl")
 include("basic/solveipform.jl")
 include("basic/solvelpform.jl")
 include("basic/cutcallback.jl")
+include("basic/pricingcallback.jl")
 
 # Child algorithms used by conquer algorithms
+include("pricing.jl")
 include("colgenstabilization.jl")
 include("colgen.jl")
 include("benders.jl")
@@ -70,9 +72,9 @@ export getterminationstatus, setterminationstatus!,
 
 # Algorithm's types
 export AbstractOptimizationAlgorithm, TreeSearchAlgorithm, ColCutGenConquer, ColumnGeneration,
-       BendersConquer, BendersCutGeneration, SolveIpForm, SolveLpForm, ExactBranchingPhase,
-       OnlyRestrictedMasterBranchingPhase, PreprocessAlgorithm, RestrictedMasterIPHeuristic,
-       OptimizationInput, OptimizationOutput, OptimizationState,
+       DefaultPricing, PricingCallback, BendersConquer, BendersCutGeneration, SolveIpForm, 
+       SolveLpForm, ExactBranchingPhase, OnlyRestrictedMasterBranchingPhase, PreprocessAlgorithm, 
+       RestrictedMasterIPHeuristic, OptimizationInput, OptimizationOutput, OptimizationState,
        EmptyInput       
 
 # Units 

src/Algorithm/basic/pricingcallback.jl

@@ -0,0 +1,50 @@
+"""
+todo
+"""
+@with_kw struct PricingCallback <: AbstractOptimizationAlgorithm
+    stage::Int = 1 # stage 1 is the exact stage by convention,
+                   # any other stage is heuristic
+    max_nb_ip_primal_sols = 50
+end
+
+function get_units_usage(algo::PricingCallback, spform::Formulation{DwSp}) 
+    units_usage = Tuple{AbstractModel, UnitType, UnitAccessMode}[] 
+    push!(units_usage, (spform, StaticVarConstrUnit, READ_ONLY))
+    return units_usage
+end
+
+function run!(algo::PricingCallback, env::Env, spform::Formulation{DwSp}, input::OptimizationInput)::OptimizationOutput
+    result = OptimizationState(
+        spform, 
+        ip_primal_bound = get_ip_primal_bound(getoptstate(input)),
+        max_length_ip_primal_sols = algo.max_nb_ip_primal_sols
+    )
+
+    @logmsg LogLevel(-2) "Calling user-defined optimization function."
+
+    optimizer = getuseroptimizer(spform)
+    cbdata = MathProg.PricingCallbackData(spform, algo.stage)
+    optimizer.user_oracle(cbdata)
+
+    if length(cbdata.primal_solutions) > 0
+        for primal_sol in cbdata.primal_solutions
+            add_ip_primal_sol!(result, primal_sol)
+        end
+
+        if algo.stage == 1 # stage 1 is exact by convention
+            dual_bound = getvalue(get_ip_primal_bound(result))
+            set_ip_dual_bound!(result, DualBound(spform, dual_bound))
+            setterminationstatus!(result, OPTIMAL) 
+        else    
+            setterminationstatus!(result, OTHER_LIMIT) 
+        end
+    else
+        if algo.stage == 1    
+            setterminationstatus!(result, INFEASIBLE) 
+        else
+            setterminationstatus!(result, OTHER_LIMIT) 
+        end 
+    end
+
+    return OptimizationOutput(result)
+end

src/Algorithm/basic/solveipform.jl

@@ -46,7 +46,8 @@ get_units_usage(algo::SolveIpForm, reform::Reformulation) =
 function check_if_optimizer_supports_ip(optimizer::MoiOptimizer)
     return MOI.supports_constraint(optimizer.inner, MOI.SingleVariable, MOI.Integer)
 end
-check_if_optimizer_supports_ip(optimizer::UserOptimizer) = true
+check_if_optimizer_supports_ip(optimizer::UserOptimizer) = false
+check_if_optimizer_supports_ip(optimizer::NoOptimizer) = false
 
 function run!(algo::SolveIpForm, env::Env, form::Formulation, input::OptimizationInput)::OptimizationOutput
     result = OptimizationState(
@@ -55,15 +56,15 @@ function run!(algo::SolveIpForm, env::Env, form::Formulation, input::Optimizatio
         max_length_ip_primal_sols = algo.max_nb_ip_primal_sols
     )
 
-    ip_supported = check_if_optimizer_supports_ip(getoptimizer(form))
+    ip_supported = check_if_optimizer_supports_ip(getmoioptimizer(form))
     if !ip_supported
         @warn "Optimizer of formulation with id =", getuid(form),
               " does not support integer variables. Skip SolveIpForm algorithm."
         setterminationstatus!(result, UNKNOWN_TERMINATION_STATUS)
         return OptimizationOutput(result)
     end
 
-    primal_sols = optimize_ip_form!(algo, getoptimizer(form), form, result)
+    primal_sols = optimize_ip_form!(algo, getmoioptimizer(form), form, result)
 
     partial_sol = nothing
     partial_sol_value = 0.0
@@ -133,7 +134,7 @@ end
 
 function optimize_with_moi!(optimizer::MoiOptimizer, form::Formulation, result::OptimizationState)
     sync_solver!(optimizer, form)
-    nbvars = MOI.get(form.optimizer.inner, MOI.NumberOfVariables())
+    nbvars = MOI.get(form.moioptimizer.inner, MOI.NumberOfVariables())
     if nbvars <= 0
         @warn "No variable in the formulation."
     end
@@ -168,19 +169,3 @@ function optimize_ip_form!(
     end
     return primal_sols
 end
-
-function optimize_ip_form!(
-    ::SolveIpForm, optimizer::UserOptimizer, form::Formulation, result::OptimizationState
-)
-    @logmsg LogLevel(-2) "Calling user-defined optimization function."
-
-    cbdata = MathProg.PricingCallbackData(form)
-    optimizer.user_oracle(cbdata)
-
-    if length(cbdata.primal_solutions) > 0
-        setterminationstatus!(result, OPTIMAL)
-    else
-        setterminationstatus!(result, INFEASIBLE) # TODO : what if no solution found ?
-    end
-    return cbdata.primal_solutions
-end

src/Algorithm/basic/solvelpform.jl

@@ -46,7 +46,7 @@ end
 function optimize_lp_form!(
     algo::SolveLpForm, optimizer::MoiOptimizer, form::Formulation, result::OptimizationState
 )
-    MOI.set(form.optimizer.inner, MOI.Silent(), algo.silent)
+    MOI.set(form.moioptimizer.inner, MOI.Silent(), algo.silent)
     optimize_with_moi!(optimizer, form, result)
     return
 end
@@ -68,7 +68,7 @@ function run!(algo::SolveLpForm, env::Env, form::Formulation, input::Optimizatio
         partial_sol_val = getvalue(partial_sol)
     end
 
-    optimizer = getoptimizer(form)
+    optimizer = getmoioptimizer(form)
     optimize_lp_form!(algo, optimizer, form, result)
     primal_sols = get_primal_solutions(form, optimizer)
 

src/Algorithm/colgen.jl

@@ -1,11 +1,7 @@
 """
     Coluna.Algorithm.ColumnGeneration(
         restr_master_solve_alg = SolveLpForm(get_dual_solution = true)
-        pricing_prob_solve_alg = SolveIpForm(
-            deactivate_artificial_vars = false,
-            enforce_integrality = false,
-            log_level = 2
-        ),
+        pricing_prob_solve_alg = DefaultPricing(),
         essential_cut_gen_alg = CutCallbacks(call_robust_facultative = false)
         max_nb_iterations::Int = 1000
         log_print_frequency::Int = 1
@@ -25,11 +21,7 @@ restricted master and `pricing_prob_solve_alg` to solve the subproblems.
     restr_master_solve_alg = SolveLpForm(get_dual_solution=true)
     # TODO : pricing problem solver may be different depending on the
     #       pricing subproblem
-    pricing_prob_solve_alg = SolveIpForm(
-        deactivate_artificial_vars=false,
-        enforce_integrality=false,
-        log_level=2
-    )
+    pricing_prob_solve_alg = DefaultPricing()
     essential_cut_gen_alg = CutCallbacks(call_robust_facultative=false)
     max_nb_iterations::Int64 = 1000
     log_print_frequency::Int64 = 1
@@ -760,7 +752,9 @@ function cg_main_loop!(
         if nb_new_columns == 0 && !essential_cuts_separated
             @logmsg LogLevel(0) "No new column generated by the pricing problems."
             setterminationstatus!(cg_optstate, OTHER_LIMIT)
-            return false, false
+            # if no columns are generated and lp gap is not closed then this col.gen. stage
+            # is a heuristic one, so we do not run phase 1 to save time
+            return true, false
         end
         if iteration > algo.max_nb_iterations
             setterminationstatus!(cg_optstate, OTHER_LIMIT)

src/Algorithm/colgenstabilization.jl

@@ -124,8 +124,16 @@ function update_stab_after_rm_solve!(
 
     unit.curalpha = smoothparam == 1.0 ? unit.basealpha : smoothparam
 
-    return linear_combination(unit.stabcenter, lp_dual_sol, unit.curalpha)
+    if !smoothing_is_active(unit)
+        # in this case Lagrangian bound calculation is simplified in col.gen.
+        # (we use the fact that the contribution of pure master variables 
+        #  is included in the value of the LP dual solution)
+        # thus, LP dual solution should be retured, as linear_combination() 
+        # does not include pure master variables contribution to the bound
+        return lp_dual_sol
     end
+    return linear_combination(unit.stabcenter, lp_dual_sol, unit.curalpha)
+end
 
 function norm(dualsol::DualSolution)
     product_sum = 0.0
@@ -175,7 +183,7 @@ function update_alpha_automatically!(
 
     # we modify the alpha parameter based on the calculated angle
     if nb_new_col == 0 || angle > 1e-12 
-        unit.basealpha -= 0.1
+        unit.basealpha = max(0.0, unit.basealpha - 0.1)
     elseif angle < -1e-12 && unit.basealpha < 0.999
         unit.basealpha += (1.0 - unit.basealpha) * 0.1
     end   
@@ -210,6 +218,7 @@ function update_stab_after_gencols!(
 
     if unit.nb_misprices > 10 || unit.curalpha <= 0.0
         unit.curalpha = 0.0
+        # stabilization is deactivated, thus we need to return the original LP dual solution
         return lp_dual_sol
     end