Colgen algorithm style, params & tests

src/algorithms/colgen.jl

@@ -1,33 +1,61 @@
 Base.@kwdef struct ColumnGeneration <: AbstractAlgorithm
-    option::Bool = false
+    max_nb_iterations::Int = 1000
+    optimality_tol::Float64 = 1e-5
 end
 
-mutable struct ColumnGenTmpRecord
+# Data stored while algorithm is running
+mutable struct ColumnGenerationAuxRec
     incumbents::Incumbents
     has_converged::Bool
     is_feasible::Bool
+    params::ColumnGeneration
 end
 
-function ColumnGenTmpRecord(S::Type{<:AbstractObjSense}, node_inc::Incumbents)
-    i = Incumbents(S)
-    update_ip_primal_sol!(i, get_ip_primal_sol(node_inc))
-    return ColumnGenTmpRecord(i, false, true)
+function ColumnGenerationAuxRec(
+    algparams::ColumnGeneration, form::Reformulation, node::Node
+)
+    inc = Incumbents(form.master.obj_sense)
+    update_ip_primal_sol!(inc, get_ip_primal_sol(node.incumbents))
+    return ColumnGenerationAuxRec(inc, false, true, algparams)
 end
 
 # Data needed for another round of column generation
-mutable struct ColumnGenerationRecord <: AbstractAlgorithmResult
+mutable struct ColumnGenerationResult <: AbstractAlgorithmResult
     incumbents::Incumbents
     proven_infeasible::Bool
 end
 
 # Overload of the algorithm's prepare function
-function prepare!(algo::ColumnGeneration, form, node)
+function prepare!(alg::ColumnGeneration, form, node)
     @logmsg LogLevel(-1) "Prepare ColumnGeneration."
     return
 end
 
-function should_do_ph_1(cg_rec::ColumnGenerationRecord)
-    primal_lp_sol = getsol(get_lp_primal_sol(cg_rec.incumbents))
+# Overload of the algorithm's run function
+function run!(alg::ColumnGeneration, form::Reformulation, node::Node)
+    @logmsg LogLevel(-1) "Run ColumnGeneration."
+    algdata = ColumnGenerationAuxRec(alg, form, node)
+    result = cg_main_loop(algdata, form, 2)
+    if should_do_ph_1(result)
+        record!(form, node)
+        set_ph_one(form.master)
+        result = cg_main_loop(algdata, form, 1)
+    end
+    if result.proven_infeasible
+        result.incumbents = Incumbents(getsense(result.incumbents))
+    end
+    if result.proven_infeasible
+        @logmsg LogLevel(-1) "ColumnGeneration terminated with status INFEASIBLE."
+    else
+        @logmsg LogLevel(-1) "ColumnGeneration terminated with status FEASIBLE."
+    end
+    update!(node.incumbents, result.incumbents)
+    return result
+end
+
+# Internal methods to the column generation
+function should_do_ph_1(result::ColumnGenerationResult)
+    primal_lp_sol = getsol(get_lp_primal_sol(result.incumbents))
     art_vars = filter(x->(getduty(x) isa ArtificialDuty), primal_lp_sol)
     if !isempty(art_vars)
         @logmsg LogLevel(-2) "Artificial variables in lp solution, need to do phase one"
@@ -45,63 +73,37 @@ function set_ph_one(master::Formulation)
     return
 end
 
-function run!(algo::ColumnGeneration, form, node)
-    @logmsg LogLevel(-1) "Run ColumnGeneration."
-    algdata = ColumnGenTmpRecord(form.master.obj_sense, node.incumbents)
-    cg_rec = cg_main_loop(algdata, form, 2)
-    if should_do_ph_1(cg_rec)
-        record!(form, node)
-        set_ph_one(form.master)
-        cg_rec = cg_main_loop(algdata, form, 1)
-    end
-    if cg_rec.proven_infeasible
-        cg_rec.incumbents = Incumbents(getsense(cg_rec.incumbents))
-    end
-    if cg_rec.proven_infeasible
-        @logmsg LogLevel(-1) "ColumnGeneration terminated with status INFEASIBLE."
-    else
-        @logmsg LogLevel(-1) "ColumnGeneration terminated with status FEASIBLE."
-    end
-    update!(node.incumbents, cg_rec.incumbents)
-    return cg_rec
-end
-
-# Internal methods to the column generation
 function update_pricing_problem!(spform::Formulation, dual_sol::DualSolution)
-
-    masterform = spform.parent_formulation
-
+    masterform = getmaster(spform)
     for (var_id, var) in Iterators.filter(_active_pricing_sp_var_ , getvars(spform))
         setcurcost!(spform, var, computereducedcost(masterform, var_id, dual_sol))
     end
-
     return false
 end
 
 function update_pricing_target!(spform::Formulation)
     # println("pricing target will only be needed after automating convexity constraints")
 end
 
-function insert_cols_in_master!(masterform::Formulation,
-                               spform::Formulation,
-                               sp_sols::Vector{PrimalSolution{S}}) where {S}
-
+function insert_cols_in_master!(
+    masterform::Formulation, spform::Formulation, 
+    spsols::Vector{PrimalSolution{S}}
+) where {S}
     sp_uid = getuid(spform)
     nb_of_gen_col = 0
-
-    for sp_sol in sp_sols
-        if contrib_improves_mlp(getbound(sp_sol))
+    for spsol in spsols
+        if contrib_improves_mlp(getbound(spsol))
             nb_of_gen_col += 1
             ref = getvarcounter(masterform) + 1
             name = string("MC", sp_uid, "_", ref)
-            resetsolvalue!(masterform, sp_sol)
+            resetsolvalue!(masterform, spsol)
             lb = 0.0
             ub = Inf
             kind = Continuous
             duty = MasterCol
             sense = Positive
             mc = setprimaldwspsol!(
-                masterform, name, sp_sol, duty; lb = lb, ub = ub,
+                masterform, name, spsol, duty; lb = lb, ub = ub,
                 kind = kind, sense = sense
             )
             @logmsg LogLevel(-2) string("Generated column : ", name)
@@ -122,17 +124,16 @@ function insert_cols_in_master!(masterform::Formulation,
             ==#
         end
     end
-
     return nb_of_gen_col
 end
 
 contrib_improves_mlp(sp_primal_bound::PrimalBound{MinSense}) = (sp_primal_bound < 0.0 - 1e-8)
 contrib_improves_mlp(sp_primal_bound::PrimalBound{MaxSense}) = (sp_primal_bound > 0.0 + 1e-8)
 
-function compute_pricing_db_contrib(spform::Formulation,
-                                    sp_sol_primal_bound::PrimalBound{S},
-                                    sp_lb::Float64,
-                                    sp_ub::Float64) where {S}
+function compute_pricing_db_contrib(
+    spform::Formulation, sp_sol_primal_bound::PrimalBound{S}, sp_lb::Float64,
+    sp_ub::Float64
+) where {S}
     # Since convexity constraints are not automated and there is no stab
     # the pricing_dual_bound_contrib is just the reduced cost * multiplicty
     if contrib_improves_mlp(sp_sol_primal_bound)
@@ -143,12 +144,10 @@ function compute_pricing_db_contrib(spform::Formulation,
     return contrib
 end
 
-function solve_sp_to_gencol!(masterform::Formulation,
-                     spform::Formulation,
-                     dual_sol::DualSolution,
-                     sp_lb::Float64,
-                     sp_ub::Float64)
-
+function solve_sp_to_gencol!(
+    masterform::Formulation, spform::Formulation, dual_sol::DualSolution,
+    sp_lb::Float64, sp_ub::Float64
+)
     #flag_need_not_generate_more_col = 0 # Not used
     flag_is_sp_infeasible = -1
     #flag_cannot_generate_more_col = -2 # Not used
@@ -198,19 +197,19 @@ function solve_sp_to_gencol!(masterform::Formulation,
     return insertion_status, pricing_db_contrib
 end
 
-function solve_sps_to_gencols!(reformulation::Reformulation,
-                  dual_sol::DualSolution{S},
-                  sp_lbs::Dict{FormId, Float64},
-                  sp_ubs::Dict{FormId, Float64}) where {S}
-
+function solve_sps_to_gencols!(
+    reform::Reformulation, dual_sol::DualSolution{S},
+    sp_lbs::Dict{FormId, Float64}, sp_ubs::Dict{FormId, Float64}
+) where {S}
     nb_new_cols = 0
     dual_bound_contrib = DualBound{S}(0.0)
-    masterform = getmaster(reformulation)
-    sps = get_dw_pricing_sp(reformulation)
+    masterform = getmaster(reform)
+    sps = get_dw_pricing_sp(reform)
     for spform in sps
         sp_uid = getuid(spform)
-        gen_status, contrib = solve_sp_to_gencol!(masterform, spform, dual_sol, sp_lbs[sp_uid], sp_ubs[sp_uid])
-
+        gen_status, contrib = solve_sp_to_gencol!(
+            masterform, spform, dual_sol, sp_lbs[sp_uid], sp_ubs[sp_uid]
+        )
         if gen_status > 0
             nb_new_cols += gen_status
             dual_bound_contrib += float(contrib)
@@ -221,15 +220,17 @@ function solve_sps_to_gencols!(reformulation::Reformulation,
     return (nb_new_cols, dual_bound_contrib)
 end
 
-function compute_master_db_contrib(algdata::ColumnGenTmpRecord,
-                                   restricted_master_sol_value::PrimalBound{S}) where {S}
+function compute_master_db_contrib(
+    algdata::ColumnGenerationAuxRec, restricted_master_sol_value::PrimalBound{S}
+) where {S}
     # TODO: will change with stabilization
     return DualBound{S}(restricted_master_sol_value)
 end
 
-function update_lagrangian_db!(algdata::ColumnGenTmpRecord,
-                               restricted_master_sol_value::PrimalBound{S},
-                               pricing_sp_dual_bound_contrib::DualBound{S}) where {S}
+function update_lagrangian_db!(
+    algdata::ColumnGenerationAuxRec, restricted_master_sol_value::PrimalBound{S},
+    pricing_sp_dual_bound_contrib::DualBound{S}
+) where {S}
     lagran_bnd = DualBound{S}(0.0)
     lagran_bnd += compute_master_db_contrib(algdata, restricted_master_sol_value)
     lagran_bnd += pricing_sp_dual_bound_contrib
@@ -245,8 +246,10 @@ function solve_restricted_master!(master::Formulation)
     getprimalsols(opt_result), getdualsols(opt_result), elapsed_time)
 end
 
-function generatecolumns!(algdata::ColumnGenTmpRecord, reform::Reformulation,
-                          master_val, dual_sol, sp_lbs, sp_ubs)
+function generatecolumns!(
+    algdata::ColumnGenerationAuxRec, reform::Reformulation, master_val, 
+    dual_sol, sp_lbs, sp_ubs
+)
     nb_new_columns = 0
     while true # TODO Replace this condition when starting implement stabilization
         nb_new_col, sp_db_contrib =  solve_sps_to_gencols!(reform, dual_sol, sp_lbs, sp_ubs)
@@ -264,21 +267,21 @@ end
 ph_one_infeasible_db(db::DualBound{MinSense}) = getvalue(db) > (0.0 + 1e-5)
 ph_one_infeasible_db(db::DualBound{MaxSense}) = getvalue(db) < (0.0 - 1e-5)
 
-function cg_main_loop(algdata::ColumnGenTmpRecord,
-                      reformulation::Reformulation, 
-                      phase::Int)::ColumnGenerationRecord
+function cg_main_loop(
+    algdata::ColumnGenerationAuxRec, reform::Reformulation, phase::Int
+)::ColumnGenerationResult
     nb_cg_iterations = 0
     # Phase II loop: Iterate while can generate new columns and
     # termination by bound does not apply
-    masterform = reformulation.master
+    masterform = reform.master
     sp_lbs = Dict{FormId, Float64}()
     sp_ubs = Dict{FormId, Float64}()
 
     # collect multiplicity current bounds for each sp
-    for spform in reformulation.dw_pricing_subprs
+    for spform in reform.dw_pricing_subprs
         sp_uid = getuid(spform)
-        lb_convexity_constr_id = reformulation.dw_pricing_sp_lb[sp_uid]
-        ub_convexity_constr_id = reformulation.dw_pricing_sp_ub[sp_uid]
+        lb_convexity_constr_id = reform.dw_pricing_sp_lb[sp_uid]
+        ub_convexity_constr_id = reform.dw_pricing_sp_ub[sp_uid]
         sp_lbs[sp_uid] = getcurrhs(getconstr(masterform, lb_convexity_constr_id))
         sp_ubs[sp_uid] = getcurrhs(getconstr(masterform, ub_convexity_constr_id))
     end
@@ -290,7 +293,7 @@ function cg_main_loop(algdata::ColumnGenTmpRecord,
         if (phase != 1 && (master_status == MOI.INFEASIBLE
             || master_status == MOI.INFEASIBLE_OR_UNBOUNDED))
             @error "Solver returned that restricted master LP is infeasible or unbounded (status = $master_status) during phase != 1."
-            return ColumnGenerationRecord(algdata.incumbents, true)
+            return ColumnGenerationResult(algdata.incumbents, true)
         end
 
         update_lp_primal_sol!(algdata.incumbents, primal_sols[1])
@@ -306,13 +309,13 @@ function cg_main_loop(algdata::ColumnGenTmpRecord,
         # generate new columns by solving the subproblems
         sp_time = @elapsed begin
             nb_new_col = generatecolumns!(
-                algdata, reformulation, master_val, dual_sols[1], sp_lbs, sp_ubs
+                algdata, reform, master_val, dual_sols[1], sp_lbs, sp_ubs
             )
         end
 
         if nb_new_col < 0
             @error "Infeasible subproblem."
-            return ColumnGenerationRecord(algdata.incumbents, true)
+            return ColumnGenerationResult(algdata.incumbents, true)
         end
 
         print_intermediate_statistics(
@@ -323,42 +326,41 @@ function cg_main_loop(algdata::ColumnGenTmpRecord,
 
         dual_bound = get_ip_dual_bound(algdata.incumbents)
         primal_bound = get_lp_primal_bound(algdata.incumbents)     
-        cur_gap = gap(primal_bound, dual_bound)
-
         ip_primal_bound = get_ip_primal_bound(algdata.incumbents)
 
-        if diff(dual_bound, ip_primal_bound) < -1e-6
+        if diff(dual_bound, ip_primal_bound) < algdata.params.optimality_tol
             algdata.has_converged = false # ???
             @logmsg LogLevel(1) "Dual bound reached primal bound."
-            return ColumnGenerationRecord(algdata.incumbents, false)
+            return ColumnGenerationResult(algdata.incumbents, false)
         end
         if phase == 1 && ph_one_infeasible_db(dual_bound)
             algdata.is_feasible = false
             @logmsg LogLevel(1) "Phase one determines infeasibility."
-            return ColumnGenerationRecord(algdata.incumbents, true)
+            return ColumnGenerationResult(algdata.incumbents, true)
         end
-        if nb_new_col == 0 || cur_gap < 0.00001 #_params_.relative_optimality_tolerance
-            @logmsg LogLevel(1) "Column Generation Algorithm has converged." #nb_new_col cur_gap
+        if nb_new_col == 0 || gap(primal_bound, dual_bound) < algdata.params.optimality_tol
+            @logmsg LogLevel(1) "Column Generation Algorithm has converged."
             algdata.has_converged = true
-            return ColumnGenerationRecord(algdata.incumbents, false)
+            return ColumnGenerationResult(algdata.incumbents, false)
         end
-        if nb_cg_iterations > 1000 ##TDalgdata.max_nb_cg_iterations
+        if nb_cg_iterations > algdata.params.max_nb_iterations
             @warn "Maximum number of column generation iteration is reached."
-            return ColumnGenerationRecord(algdata.incumbents, false)
+            return ColumnGenerationResult(algdata.incumbents, false)
         end
     end
-    return ColumnGenerationRecord(algdata.incumbents, false)
+    return ColumnGenerationResult(algdata.incumbents, false)
 end
 
-function print_intermediate_statistics(algdata::ColumnGenTmpRecord,
-                                       nb_new_col::Int,
-                                       nb_cg_iterations::Int,
-                                       mst_time::Float64, sp_time::Float64)
+function print_intermediate_statistics(
+    algdata::ColumnGenerationAuxRec, nb_new_col::Int, nb_cg_iterations::Int,
+    mst_time::Float64, sp_time::Float64
+)
     mlp = getvalue(get_lp_primal_bound(algdata.incumbents))
     db = getvalue(get_ip_dual_bound(algdata.incumbents))
     pb = getvalue(get_ip_primal_bound(algdata.incumbents))
     @printf(
-            "<it=%i> <et=%i> <mst=%.3f> <sp=%.3f> <cols=%i> <mlp=%.4f> <DB=%.4f> <PB=%.4f>\n",
-            nb_cg_iterations, _elapsed_solve_time(), mst_time, sp_time, nb_new_col, mlp, db, pb
+        "<it=%i> <et=%i> <mst=%.3f> <sp=%.3f> <cols=%i> <mlp=%.4f> <DB=%.4f> <PB=%.4f>\n",
+        nb_cg_iterations, _elapsed_solve_time(), mst_time, sp_time, nb_new_col, mlp, db, pb
     )
+    return
 end

test/full_instances_tests.jl

@@ -43,6 +43,32 @@ function generalized_assignment_tests()
         @test CLD.GeneralizedAssignment.print_and_check_sol(data, problem, x)
     end
 
+    @testset "gap - ColGen max nb iterations" begin
+        data = CLD.GeneralizedAssignment.data("smallgap3.txt")
+
+        coluna = JuMP.with_optimizer(
+            Coluna.Optimizer, params = CL.Params(
+                global_strategy = CL.GlobalStrategy(
+                    CL.SimpleBnP(
+                        colgen = CL.ColumnGeneration(
+                            max_nb_iterations = 8
+                        )
+                    ),
+                    CL.SimpleBranching(), 
+                    CL.DepthFirst()
+                )
+            ),
+            default_optimizer = with_optimizer(GLPK.Optimizer)
+        )
+
+        problem, x, dec = CLD.GeneralizedAssignment.model(data, coluna)
+
+        JuMP.optimize!(problem)
+        @test abs(JuMP.objective_value(problem) - 438.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
+
     @testset "gap with penalties - pure master variables" begin
         data = CLD.GeneralizedAssignment.data("smallgap3.txt")