Fix ColGen algorithm

src/Algorithm/colgen.jl

@@ -9,14 +9,15 @@ end
 # Data stored while algorithm is running
 mutable struct ColGenRuntimeData
     optstate::OptimizationState
+    nb_iterations::Int64
     phase::Int64
 end
 
 function ColGenRuntimeData(
     algo::ColumnGeneration, reform::Reformulation, nodestate::OptimizationState
 )
     optstate = CopyBoundsAndStatusesFromOptState(getmaster(reform), nodestate, false)
-    return ColGenRuntimeData(optstate, 2)
+    return ColGenRuntimeData(optstate, 0, 3)
 end
 
 struct ReducedCostsVector
@@ -33,7 +34,7 @@ function ReducedCostsVector(varids::Vector{VarId}, form::Vector{Formulation})
     permute!(varids, p)
     permute!(form, p)
     for i in 1:len
-        perenecosts[i] = getperenecost(form[i], varids[i])
+        perenecosts[i] = getcurcost(getmaster(form[i]), varids[i])
     end
     return ReducedCostsVector(len, varids, perenecosts, form)
 end
@@ -43,24 +44,26 @@ getoptstate(data::ColGenRuntimeData) = data.optstate
 function run!(algo::ColumnGeneration, reform::Reformulation, input::OptimizationInput)::OptimizationOutput    
     data = ColGenRuntimeData(algo, reform, getoptstate(input))
     optstate = getoptstate(data)
-
-    cg_main_loop!(algo, data, reform)
     masterform = getmaster(reform)
-    if should_do_ph_1(masterform, data)
-        set_ph_one(masterform, data)        
-        cg_main_loop!(algo, data, reform)
-        # TO DO : to implement unsetting phase one !!
-        # TO DO : to implement repeating of phase two in the case phase 1 succeded
+
+    set_ph3!(masterform, data) # mixed ph1 & ph2
+    stop = cg_main_loop!(algo, data, reform)
+
+    if !stop && should_do_ph_1(masterform, data) 
+        set_ph1!(masterform, data) # pure ph1 (mixed phase infeasible or artifical variables not enough expensive)
+        stop = cg_main_loop!(algo, data, reform)
+        if !stop
+            set_ph2!(masterform, data) # pure ph2
+            cg_main_loop!(algo, data, reform)
+        end
     end
 
     @logmsg LogLevel(-1) string("ColumnGeneration terminated with status ", getfeasibilitystatus(optstate))
 
     return OptimizationOutput(optstate)
 end
 
-# Internal methods to the column generation
 function should_do_ph_1(master::Formulation, data::ColGenRuntimeData)
-    ip_gap(getoptstate(data)) <= 0.00001 && return false
     primal_lp_sol = get_lp_primal_sols(getoptstate(data))[1]
     if contains(primal_lp_sol, vid -> isanArtificialDuty(getduty(vid)))
         @logmsg LogLevel(-2) "Artificial variables in lp solution, need to do phase one"
@@ -71,12 +74,40 @@ function should_do_ph_1(master::Formulation, data::ColGenRuntimeData)
     end
 end
 
-function set_ph_one(master::Formulation, data::ColGenRuntimeData)
+function set_ph1!(master::Formulation, data::ColGenRuntimeData)
     for (varid, var) in getvars(master)
-        isanArtificialDuty(getduty(varid)) && continue
-        setcurcost!(master, varid, 0.0)
+        if !isanArtificialDuty(getduty(varid))
+            setcurcost!(master, varid, 0.0)
+        end
     end
     data.phase = 1
+    set_lp_dual_bound!(data.optstate, DualBound(master))
+    set_ip_dual_bound!(data.optstate, DualBound(master))
+    return
+end
+
+function set_ph2!(master::Formulation, data::ColGenRuntimeData)
+    for (varid, var) in getvars(master)
+        if isanArtificialDuty(getduty(varid))
+            deactivate!(master, varid)
+        else
+            setcurcost!(master, varid, getperenecost(master, var))
+        end
+    end
+    set_lp_dual_bound!(data.optstate, DualBound(master))
+    set_ip_dual_bound!(data.optstate, DualBound(master))
+    data.phase = 2
+end
+
+function set_ph3!(master::Formulation, data::ColGenRuntimeData)
+    for (varid, var) in getvars(master)
+        if isanArtificialDuty(getduty(varid))
+            activate!(master, varid)
+        else
+            setcurcost!(master, varid, getperenecost(master, var))
+        end
+    end
+    data.phase = 3
     return
 end
 
@@ -85,7 +116,7 @@ function update_pricing_target!(spform::Formulation)
 end
 
 function insert_cols_in_master!(
-    masterform::Formulation, spform::Formulation, sp_solution_ids::Vector{VarId}
+    data::ColGenRuntimeData, masterform::Formulation, spform::Formulation, sp_solution_ids::Vector{VarId}
 ) 
     sp_uid = getuid(spform)
     nb_of_gen_col = 0
@@ -102,6 +133,9 @@ function insert_cols_in_master!(
             masterform, spform, sol_id, name, duty; lb = lb, ub = ub, 
             kind = kind, sense = sense
         )
+        if data.phase == 1
+            setcurcost!(masterform, mc, 0.0)
+        end
         @logmsg LogLevel(-2) string("Generated column : ", name)
     end
 
@@ -265,8 +299,8 @@ function updatereducedcosts!(reform::Reformulation, redcostsvec::ReducedCostsVec
 end
 
 function solve_sps_to_gencols!(
-    algo::ColumnGeneration, reform::Reformulation, redcostsvec::ReducedCostsVector, dual_sol::DualSolution, 
-    sp_lbs::Dict{FormId, Float64}, sp_ubs::Dict{FormId, Float64}
+    algo::ColumnGeneration, data::ColGenRuntimeData, reform::Reformulation, redcostsvec::ReducedCostsVector, 
+    dual_sol::DualSolution, sp_lbs::Dict{FormId, Float64}, sp_ubs::Dict{FormId, Float64}
 )
     masterform = getmaster(reform)
     nb_new_cols = 0
@@ -296,7 +330,7 @@ function solve_sps_to_gencols!(
     nb_new_cols = 0
     for (spuid, spform) in sps
         dual_bound_contrib += sp_dual_bound_contribs[spuid]
-        nb_new_cols += insert_cols_in_master!(masterform, spform, recorded_sp_solution_ids[spuid])
+        nb_new_cols += insert_cols_in_master!(data, masterform, spform, recorded_sp_solution_ids[spuid])
         for colid in sp_solution_to_activate[spuid]
             activate!(masterform, colid)
             nb_new_cols += 1
@@ -330,7 +364,7 @@ function generatecolumns!(
     cg_optstate = getoptstate(data)
     nb_new_columns = 0
     while true # TODO Replace this condition when starting implement stabilization
-        nb_new_col, sp_db_contrib =  solve_sps_to_gencols!(algo, reform, redcostsvec, dual_sol, sp_lbs, sp_ubs)
+        nb_new_col, sp_db_contrib =  solve_sps_to_gencols!(algo, data, reform, redcostsvec, dual_sol, sp_lbs, sp_ubs)
         nb_new_columns += nb_new_col
         lagran_bnd = calculate_lagrangian_db(data, master_val, sp_db_contrib)
         update_ip_dual_bound!(cg_optstate, lagran_bnd)
@@ -344,13 +378,13 @@ function generatecolumns!(
     return nb_new_columns
 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)
+ph_one_infeasible_db(algo, db::DualBound{MinSense}) = getvalue(db) > algo.optimality_tol
+ph_one_infeasible_db(algo, db::DualBound{MaxSense}) = getvalue(db) < - algo.optimality_tol
 
 function cg_main_loop!(algo::ColumnGeneration, data::ColGenRuntimeData, reform::Reformulation)
-    nb_cg_iterations = 0
     # Phase II loop: Iterate while can generate new columns and
     # termination by bound does not apply
+
     masterform = getmaster(reform)
     sp_lbs = Dict{FormId, Float64}()
     sp_ubs = Dict{FormId, Float64}()
@@ -391,7 +425,7 @@ function cg_main_loop!(algo::ColumnGeneration, data::ColGenRuntimeData, reform::
             @warn string("Solver returned that LP restricted master is infeasible or unbounded ",
             "(feasibility status = " , status, ") during phase != 1.")
             setfeasibilitystatus!(cg_optstate, status) 
-            return
+            return true
         end
 
         lp_dual_sol = DualSolution(masterform)
@@ -416,7 +450,7 @@ function cg_main_loop!(algo::ColumnGeneration, data::ColGenRuntimeData, reform::
 
         # TODO: cleanup restricted master columns        
 
-        nb_cg_iterations += 1
+        data.nb_iterations += 1
 
         # generate new columns by solving the subproblems
         sp_time = @elapsed begin
@@ -428,10 +462,10 @@ function cg_main_loop!(algo::ColumnGeneration, data::ColGenRuntimeData, reform::
         if nb_new_col < 0
             @error "Infeasible subproblem."
             setfeasibilitystatus!(cg_optstate, INFEASIBLE) 
-            return 
+            return true
         end
 
-        print_colgen_statistics(cg_optstate, nb_new_col, nb_cg_iterations, rm_time, sp_time)
+        print_colgen_statistics(data, cg_optstate, nb_new_col, rm_time, sp_time)
 
         # TODO: update colgen stabilization
 
@@ -442,41 +476,47 @@ function cg_main_loop!(algo::ColumnGeneration, data::ColGenRuntimeData, reform::
         if ip_gap(cg_optstate) < algo.optimality_tol
             setterminationstatus!(cg_optstate, OPTIMAL) 
             @logmsg LogLevel(0) "Dual bound reached primal bound."
-            return 
+            return true
         end
-        if data.phase == 1 && ph_one_infeasible_db(dual_bound)
+        if data.phase == 1 && ph_one_infeasible_db(algo, dual_bound)
             db = - getvalue(DualBound(reform))
             pb = - getvalue(PrimalBound(reform))
             set_lp_dual_bound!(cg_optstate, DualBound(reform, db))
             set_lp_primal_bound!(cg_optstate, PrimalBound(reform, pb))
             setfeasibilitystatus!(cg_optstate, INFEASIBLE) 
             @logmsg LogLevel(0) "Phase one determines infeasibility."
-            return 
+            return true
         end
         if nb_new_col == 0 || lp_gap(cg_optstate) < algo.optimality_tol
             @logmsg LogLevel(0) "Column Generation Algorithm has converged."
             setterminationstatus!(cg_optstate, OPTIMAL) 
-            return 
+            return false
         end
-        if nb_cg_iterations > algo.max_nb_iterations
+        if data.nb_iterations > algo.max_nb_iterations
             setterminationstatus!(cg_optstate, OTHER_LIMIT)
             @warn "Maximum number of column generation iteration is reached."
-            return 
+            return true
         end
     end
-    return 
+    return false
 end
 
 function print_colgen_statistics(
-    optstate::OptimizationState, nb_new_col::Int, nb_cg_iterations::Int,
-    mst_time::Float64, sp_time::Float64
+    data::ColGenRuntimeData, optstate::OptimizationState, nb_new_col::Int, mst_time::Float64, sp_time::Float64
 )
     mlp = getvalue(get_lp_primal_bound(optstate))
     db = getvalue(get_lp_dual_bound(optstate))
     pb = getvalue(get_ip_primal_bound(optstate))
+    phase = "  "
+    if data.phase == 1
+        phase = "# "
+    elseif data.phase == 2
+        phase = "##"
+    end
+
     @printf(
-        "<it=%3i> <et=%5.2f> <mst=%5.2f> <sp=%5.2f> <cols=%2i> <DB=%10.4f> <mlp=%10.4f> <PB=%.4f>\n",
-        nb_cg_iterations, Coluna._elapsed_solve_time(), mst_time, sp_time, nb_new_col, db, mlp, pb
+        "%s<it=%3i> <et=%5.2f> <mst=%5.2f> <sp=%5.2f> <cols=%2i> <DB=%10.4f> <mlp=%10.4f> <PB=%.4f>\n",
+        phase, data.nb_iterations, Coluna._elapsed_solve_time(), mst_time, sp_time, nb_new_col, db, mlp, pb
     )
     return
 end

src/MathProg/MOIinterface.jl

@@ -165,9 +165,8 @@ function call_moi_optimize_with_silence(optimizer::MoiOptimizer)
     return
 end
 
-function remove_from_optimizer!(optimizer::MoiOptimizer,
-                                var::Variable)
-    inner = getinner(optimizer)
+function remove_from_optimizer!(form::Formulation, var::Variable)                       
+    inner = getinner(form.optimizer)
     moirecord = getmoirecord(var)
     @assert getindex(moirecord).value != -1
     MOI.delete(inner, getbounds(moirecord))
@@ -179,11 +178,10 @@ function remove_from_optimizer!(optimizer::MoiOptimizer,
     return
 end
 
-function remove_from_optimizer!(optimizer::MoiOptimizer,
-                                constr::Constraint)
+function remove_from_optimizer!(form::Formulation, constr::Constraint)
     moirecord = getmoirecord(constr)
     @assert getindex(moirecord).value != -1
-    MOI.delete(getinner(optimizer), getindex(moirecord))
+    MOI.delete(getinner(form.optimizer), getindex(moirecord))
     setindex!(moirecord, MoiConstrIndex())
     return
 end

src/MathProg/decomposition.jl

@@ -4,7 +4,7 @@ function set_glob_art_var(form::Formulation, is_pos::Bool)
     cost *= getobjsense(form) == MinSense ? 1.0 : -1.0
     return setvar!(
         form, name, MasterArtVar; 
-        cost = cost, lb = 0.0, ub = Inf, kind = Continuous, sense = Positive
+        cost = cost, lb = 0.0, ub = Inf, kind = Continuous, sense = Positive 
     )
 end
 

src/MathProg/formulation.jl

@@ -409,7 +409,7 @@ function remove_from_optimizer!(ids::Set{Id{T}}, form::Formulation) where {
     for id in ids
         vc = getelem(form, id)
         @logmsg LogLevel(-3) string("Removing varconstr of name ", getname(form, vc))
-        remove_from_optimizer!(form.optimizer, vc)
+        remove_from_optimizer!(form, vc)
     end
     return
 end
@@ -448,7 +448,7 @@ function optimize!(form::Formulation)
     @logmsg LogLevel(-1) string("Optimizing formulation ", getuid(form))
     @logmsg LogLevel(-3) form
     res = optimize!(form, getoptimizer(form))
-    @logmsg LogLevel(-2) "Optimization finished with result:" print(form, res)
+    @logmsg LogLevel(-3) "Optimization finished with result:" print(form, res)
     return res
 end
 

src/MathProg/new_varconstr.jl

@@ -21,7 +21,7 @@ Set the current cost of variable `var` with id `id` to `cost` in formulation
 """
 function setcurcost!(form::Formulation, varid::VarId, cost::Float64)
     form.manager.var_datas[varid].cost = cost
-    if iscurexplicit(form, varid) 
+    if iscurexplicit(form, varid) && iscuractive(form, varid)
         change_cost!(form.buffer, varid)
     end
     return
@@ -54,7 +54,7 @@ according to `new_lb`. Change on `f.optimizer` will be buffered.
 """
 function setcurlb!(form::Formulation, varid::VarId, lb::Float64)
     form.manager.var_datas[varid].lb = lb
-    if iscurexplicit(form, varid) 
+    if iscurexplicit(form, varid) && iscuractive(form, varid)
         change_bound!(form.buffer, varid)
     end
     return
@@ -85,7 +85,7 @@ according to `new_ub`. Change on `f.optimizer` will be buffered.
 """
 function setcurub!(form::Formulation, varid::VarId, ub::Float64)
     form.manager.var_datas[varid].ub = ub
-    if iscurexplicit(form, varid) 
+    if iscurexplicit(form, varid) && iscuractive(form, varid)
         change_bound!(form.buffer, varid)
     end
     return
@@ -105,7 +105,7 @@ getcurrhs(form::Formulation, constrid::ConstrId) = form.manager.constr_datas[con
 getcurrhs(form::Formulation, constr::Constraint) = getcurrhs(form, getid(constr))
 function setcurrhs!(form::Formulation, constrid::ConstrId, rhs::Float64) 
     form.manager.constr_datas[constrid].rhs = rhs
-    if iscurexplicit(form, constrid) 
+    if iscurexplicit(form, constrid) && iscuractive(form, constrid)
         change_rhs!(form.buffer, constrid)
     end
     return
@@ -140,15 +140,15 @@ according to `new_kind`. Change on `f.optimizer` will be buffered.
 """
 function setcurkind!(form::Formulation, varid::VarId, kind::VarKind)
     form.manager.var_datas[varid].kind = kind
-    if iscurexplicit(form, varid) 
+    if iscurexplicit(form, varid) && iscuractive(form, varid)
         change_kind!(form.buffer, varid)
     end
     return
 end
 setcurkind!(form::Formulation, var::Variable, kind::VarKind) = setcurkind!(form, getid(var), kind)
 function setcurkind!(form::Formulation, constrid::ConstrId, kind::ConstrKind)
     form.manager.constr_datas[constrid].kind = kind
-    if iscurexplicit(form, constrid) 
+    if iscurexplicit(form, constrid) && iscuractive(form, constrid)
         change_kind!(form.buffer, constrid)
     end
     return
@@ -285,7 +285,7 @@ end
 Deactivate a variable or a constraint in the formulation
 """
 function deactivate!(form::Formulation, varconstrid::Id{VC}) where {VC<:AbstractVarConstr}
-    if iscurexplicit(form, varconstrid)
+    if iscurexplicit(form, varconstrid) && iscuractive(form, varconstrid)
         remove!(form.buffer, varconstrid)
     end
     _setiscuractive!(form, varconstrid, false)