ColGen Iteration & Phases API

src/Algorithm/Algorithm.jl

@@ -4,7 +4,7 @@ using DataStructures
 import MathOptInterface
 import TimerOutputs
 
-using ..Coluna, ..ColunaBase, ..MathProg, ..MustImplement
+using ..Coluna, ..ColunaBase, ..MathProg, ..MustImplement, ..ColGen
 
 using Crayons, DynamicSparseArrays, Logging, Parameters, Printf, Random, Statistics, SparseArrays, LinearAlgebra
 
@@ -41,6 +41,12 @@ include("basic/cutcallback.jl")
 
 # Child algorithms used by conquer algorithms
 include("colgenstabilization.jl")
+
+# Column generation
+include("colgen/utils.jl")
+include("colgen/default.jl")
+
+
 include("colgen.jl")
 include("benders.jl")
 

src/Algorithm/colgen.jl

@@ -112,148 +112,8 @@ function get_units_usage(algo::ColumnGeneration, reform::Reformulation)
     return units_usage
 end
 
-############################################################################################
-# Errors and warnings
-############################################################################################
-
-"""
-Error thrown when when a subproblem generates a column with negative (resp. positive) 
-reduced cost in min (resp. max) problem that already exists in the master 
-and that is already active. 
-An active master column cannot have a negative reduced cost.
-"""
-struct ColumnAlreadyInsertedColGenWarning
-    column_in_master::Bool
-    column_is_active::Bool
-    column_reduced_cost::Float64
-    column_id::VarId
-    master::Formulation{DwMaster}
-    subproblem::Formulation{DwSp}
-end
-
-function Base.show(io::IO, err::ColumnAlreadyInsertedColGenWarning)
-    msg = """
-    Unexpected variable state during column insertion.
-    ======
-    Column id: $(err.column_id).
-    Reduced cost of the column: $(err.column_reduced_cost).
-    The column is in the master ? $(err.column_in_master).
-    The column is active ? $(err.column_is_active).
-    ======
-    If the column is in the master and active, it means a subproblem found a solution with
-    negative (minimization) / positive (maximization) reduced cost that is already active in
-    the master. This should not happen.
-    ======
-    If you are using a pricing callback, make sure there is no bug in your code.
-    If you are using a solver (e.g. GLPK, Gurobi...), check the reduced cost tolerance 
-    `redcost_tol` parameter of `ColumnGeneration`.
-    If you find a bug in Coluna, please open an issue at https://github.com/atoptima/Coluna.jl/issues with an example
-    that reproduces the bug.
-    ======
-    """
-    println(io, msg)
-end
-
-############################################################################################
-# Information extracted to speed-up some computations.
-############################################################################################
-function _submatrix(
-    form::Formulation, 
-    keep_constr::Function, 
-    keep_var::Function,
-)
-    matrix = getcoefmatrix(form)
-    constr_ids = ConstrId[]
-    var_ids = VarId[]
-    nz = Float64[]
-    for constr_id in Iterators.filter(keep_constr, Iterators.keys(getconstrs(form)))
-        for (var_id, coeff) in @view matrix[constr_id, :]
-            if keep_var(var_id)
-                push!(constr_ids, constr_id)
-                push!(var_ids, var_id)
-                push!(nz, coeff)
-            end
-        end
-    end
-    return dynamicsparse(
-        constr_ids, var_ids, nz, ConstrId(Coluna.MAX_NB_ELEMS), VarId(Coluna.MAX_NB_ELEMS)
-    )
-end
-
-"""
-Extracted information to speed-up calculation of reduced costs of master variables.
-We extract from the master only the information we need to compute the reduced cost of DW 
-subproblem variables:
-- `c` contains the perenial cost of DW subproblem representative variables
-- `A` is a submatrix of the master coefficient matrix that involves only DW subproblem
-  representative variables.
-
-Calculation is `c - transpose(A) * master_lp_dual_solution`.
-"""
-struct ReducedCostsCalculationHelper
-    c::SparseVector{Float64,VarId}
-    A::DynamicSparseMatrix{ConstrId,VarId,Float64}
-end
-
-function ReducedCostsCalculationHelper(master)
-    dwspvar_ids = VarId[]
-    peren_costs = Float64[]
-
-    for var_id in Iterators.keys(getvars(master))
-        if iscuractive(master, var_id) && getduty(var_id) <= AbstractMasterRepDwSpVar
-            push!(dwspvar_ids, var_id)
-            push!(peren_costs, getcurcost(master, var_id))
-        end
-    end
-
-    dwsprep_costs = sparsevec(dwspvar_ids, peren_costs, Coluna.MAX_NB_ELEMS)
-    dwsprep_coefmatrix = _submatrix(
-        master, 
-        constr_id -> !(getduty(constr_id) <= MasterConvexityConstr),
-        var_id -> getduty(var_id) <= AbstractMasterRepDwSpVar
-    )
-    return ReducedCostsCalculationHelper(dwsprep_costs, dwsprep_coefmatrix)
-end
-
-"""
-Precompute information to speed-up calculation of subgradient of master variables.
-We extract from the master follwowing information:
-- `a` contains the perenial rhs of all master constraints except convexity constraints;
-- `A` is a submatrix of the master coefficient matrix that involves only representative of
-  original variables (pure master vars + DW subproblem represtative vars) 
-
-Calculation is `a - A * (m .* z)`
-where :
- - `m` contains a multiplicity factor for each variable involved in the calculation
-       (lower or upper sp multiplicity depending on variable reduced cost);
- - `z` is the concatenation of the solution to the master (for pure master vars) and pricing
-       subproblems (for DW subproblem represtative vars).
-
-Operation `m .* z` "mimics" a solution in the original space.
-"""
-struct SubgradientCalculationHelper
-    a::SparseVector{Float64,ConstrId}
-    A::DynamicSparseMatrix{ConstrId,VarId,Float64}
-end
-
-function SubgradientCalculationHelper(master)
-    constr_ids = ConstrId[]
-    constr_rhs = Float64[]
-    for (constr_id, constr) in getconstrs(master)
-        if !(getduty(constr_id) <= MasterConvexityConstr) && 
-           iscuractive(master, constr) && isexplicit(master, constr)
-            push!(constr_ids, constr_id)
-            push!(constr_rhs, getcurrhs(master, constr_id))
-        end 
-    end
-    a = sparsevec(constr_ids, constr_rhs, Coluna.MAX_NB_ELEMS)
-    A = _submatrix(
-        master, 
-        constr_id -> !(getduty(constr_id) <= MasterConvexityConstr),
-        var_id -> getduty(var_id) <= MasterPureVar || getduty(var_id) <= MasterRepPricingVar
-    )
-    return SubgradientCalculationHelper(a, A)
-end
+#####
+include("colgen/utils.jl")
 
 ############################################################################################
 # Column generation algorithm.

src/Algorithm/colgen/default.jl

@@ -0,0 +1,164 @@
+struct ColGenContext <: ColGen.AbstractColGenContext
+    # Information to solve the master
+    master_solve_alg
+    master_optimizer_id
+
+    # Memoization to compute reduced costs (this is a precompute)
+    redcost_mem
+end
+
+###############################################################################
+# Sequence of phases
+###############################################################################
+struct ColunaColGenPhaseIterator <: ColGen.AbstractColGenPhaseIterator end
+
+"""
+Phase 1 sets the cost of variables to 0 except for artifical variables.
+The goal is to find a solution to the master LP problem that has no artificial variables.
+"""
+struct ColGenPhase1 <: ColGen.AbstractColGenPhase end
+
+"""
+Phase 2 solves the master LP without artificial variables.
+To starts, it requires a set of columns that forms a feasible solution to the LP master.
+This set is found with phase 1.
+"""
+struct ColGenPhase2 <: ColGen.AbstractColGenPhase end
+
+"""
+Phase 3 is a mix of phase 1 and phase 2.
+It sets a very large cost to artifical variable to force them to be removed from the master 
+LP solution.
+If the final master LP solution contains artifical variables either the master is infeasible
+or the cost of artificial variables is not large enough. Phase 1 must be run.
+"""
+struct ColGenPhase3 <: ColGen.AbstractColGenPhase end
+
+# Implementation of ColGenPhase interface
+## Implementation of `initial_phase`.
+ColGen.initial_phase(::ColunaColGenPhaseIterator) = ColGenPhase3()
+
+## Implementation of `next_phase`.
+function ColGen.next_phase(::ColunaColGenPhaseIterator, ::ColGenPhase1, ctx)
+    # If master LP solution has no artificial vars, it means that the phase 1 has succeeded.
+    # We have a set of columns that forms a feasible solution to the LP master and we can 
+    # thus start phase 2.
+    if !colgen_mast_lp_sol_has_art_vars(ctx)
+        return ColGenPhase2()
+    end
+    return nothing
+end
+
+function ColGen.next_phase(::ColunaColGenPhaseIterator, ::ColGenPhase2, ctx)
+    # The phase 2 is always the last phase of the column generation algorithm.
+    # It means the algorithm converged or hit a user limit.
+    return nothing
+end
+
+function ColGen.next_phase(::ColunaColGenPhaseIterator, ::ColGenPhase3, ctx)
+    # Master LP solution has artificial vars.
+    if colgen_mast_lp_sol_has_art_vars(ctx)
+        return ColGenPhase1()
+    end
+    return nothing
+end
+
+## Methods used in the implementation and that we should mock in tests.
+function colgen_mast_lp_sol_has_art_vars(ctx::ColGenContext)
+
+end
+
+## Implementatation of `setup_reformulation!`
+## Phase 1 => non-artifical variables have cost equal to 0
+function ColGen.setup_reformulation!(reform, ::ColGenPhase1)
+    master = getmaster(reform)
+    for (varid, _) in getvars(master)
+        if !isanArtificialDuty(getduty(varid))
+            setcurcost!(master, varid, 0.0)
+        end
+    end
+    return
+end
+
+## Phase 2 => deactivate artifical variables and make sure that the cost of non-artifical
+## variables is correct.
+function ColGen.setup_reformulation!(reform, ::ColGenPhase2)
+    master = getmaster(reform)
+    for (varid, var) in getvars(master)
+        if isanArtificialDuty(getduty(varid))
+            deactivate!(master, varid)
+        else
+            setcurcost!(master, varid, getperencost(master, var))
+        end
+    end
+    return
+end
+
+## Phase 3 => make sure artifical variables are active and cost is correct.
+function ColGen.setup_reformulation!(reform, ::ColGenPhase3)
+    master = getmaster(reform)
+    for (varid, var) in getvars(master)
+        if isanArtificialDuty(getduty(varid))
+            activate!(master, varid)
+        end
+        setcurcost!(master, varid, getperencost(master, var))
+    end
+    return
+end
+
+######### Pricing strategy
+struct ClassicPricingStrategy <: ColGen.AbstractPricingStrategy 
+    subprobs::Dict{FormId, Formulation{DwSp}}
+end
+
+function ColGen.get_pricing_strategy(ctx::ColGen.AbstractColGenContext, _)
+    ClassicPricingStrategy(Dict(i => sp for (i, sp) in ColGen.get_pricing_subprobs(ctx)))
+end
+
+ColGen.pricing_strategy_iterate(ps::ClassicPricingStrategy) = iterate(ps.subprobs)
+
+
+######### Column generation
+
+# Placeholder methods:  
+ColGen.before_colgen_iteration(::ColGenContext, _, _) = nothing
+ColGen.after_colgen_iteration(::ColGenContext, _, _, _) = nothing
+
+######### Column generation iteration
+function ColGen.optimize_master_lp_problem!(master, context, env)
+    println("\e[31m optimize master lp problem \e[00m")
+    input = OptimizationState(master, ip_primal_bound=0.0) # TODO : ip_primal_bound=get_ip_primal_bound(cg_optstate)
+    return run!(context.master_solve_alg, env, master, input, context.master_optimizer_id)
+end
+
+#get_primal_sol(mast_result)
+
+function ColGen.check_primal_ip_feasibility(ctx, mast_lp_primal_sol)
+    println("\e[31m check primal ip feasibility \e[00m")
+    return !contains(mast_lp_primal_sol, varid -> isanArtificialDuty(getduty(varid))) &&
+        isinteger(proj_cols_on_rep(mast_lp_primal_sol, getmodel(mast_lp_primal_sol)))
+end
+
+#update_inc_primal_sol!
+
+#get_dual_sol(mast_result)
+
+function ColGen.update_master_constrs_dual_vals!(ctx, master, smooth_dual_sol)
+    println("\e[32m update_master_constrs_dual_vals \e[00m")
+    # Set all dual value of all constraints to 0.
+    for constr in Iterators.values(getconstrs(master))
+        setcurincval!(master, constr, 0.0)
+    end
+    # Update constraints that have non-zero dual values.
+    for (constr_id, val) in smooth_dual_sol
+        setcurincval!(master, constr_id, val)
+    end
+end
+
+function ColGen.update_sp_vars_red_costs!(ctx, sp, red_costs)
+    println("\e[34m update_sp_vars_red_costs \e[00m")
+    for (var_id, _) in getvars(sp)
+        setcurcost!(sp, var_id, red_costs[var_id])
+    end
+    return
+end