Fix ColGen iteration test for gap with pure master vars

src/Algorithm/colgen/default.jl

@@ -176,8 +176,8 @@ function ColGen.check_primal_ip_feasibility(master_lp_primal_sol, phase, reform)
 end
 
 # Reduced costs calculation
-ColGen.get_orig_costs(ctx::ColGenContext) = ctx.reduced_cost_helper.c
-ColGen.get_coef_matrix(ctx::ColGenContext) = ctx.reduced_cost_helper.A
+ColGen.get_subprob_var_orig_costs(ctx::ColGenContext) = ctx.reduced_cost_helper.dw_subprob_c
+ColGen.get_subprob_var_coef_matrix(ctx::ColGenContext) = ctx.reduced_cost_helper.dw_subprob_A
 
 function ColGen.update_sp_vars_red_costs!(ctx::ColGenContext, sp::Formulation{DwSp}, red_costs)
     for (var_id, _) in getvars(sp)
@@ -234,6 +234,7 @@ end
 struct ColGenPricingResult{F,S}
     result::OptimizationState{F,S}
     columns::Vector{GeneratedColumn}
+    best_red_cost::Float64
 end
 
 function ColGen.is_infeasible(pricing_res::ColGenPricingResult)
@@ -247,9 +248,10 @@ function ColGen.is_unbounded(pricing_res::ColGenPricingResult)
 end
 
 ColGen.get_primal_sols(pricing_res) = pricing_res.columns
-ColGen.get_dual_bound(pricing_res) = get_ip_dual_bound(pricing_res.result)
+ColGen.get_dual_bound(pricing_res) = pricing_res.best_red_cost
 
-has_improving_red_cost(column::GeneratedColumn) = column.red_cost < 0
+is_improving_red_cost(red_cost) = red_cost < 0
+has_improving_red_cost(column::GeneratedColumn) = is_improving_red_cost(column.red_cost)
 # In our implementation of `push_in_set!`, we keep only columns that have improving reduced 
 # cost.
 function ColGen.push_in_set!(pool, column)
@@ -270,16 +272,20 @@ function ColGen.optimize_pricing_problem!(ctx::ColGenContext, sp::Formulation{Dw
     # (C) the contribution of the pure master variables.
 
     # Master convexity constraints contribution.
+    # TODO: talk with fv & Ruslan because this way to take into account convexity constraints has
+    # drawbacks (numerical stability).
     lb_dual = master_dual_sol[sp.duty_data.lower_multiplicity_constr_id]
     ub_dual = master_dual_sol[sp.duty_data.upper_multiplicity_constr_id]
 
     # Pure master variables contribution.
-    # TODO
+    # TODO (only when stabilization is used otherwise already taken into account by master obj val)
 
     generated_columns = GeneratedColumn[]
     for col in get_ip_primal_sols(opt_state)
         red_cost = getvalue(col) - lb_dual - ub_dual
         push!(generated_columns, GeneratedColumn(col, red_cost))
     end
-    return ColGenPricingResult(opt_state, generated_columns)
+
+    best_red_cost = getvalue(get_ip_dual_bound(opt_state)) - lb_dual - ub_dual
+    return ColGenPricingResult(opt_state, generated_columns, best_red_cost)
 end

src/Algorithm/colgen/utils.jl

@@ -67,43 +67,67 @@ function _submatrix(
 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 
+Extracted information to speed-up calculation of reduced costs of subproblem representatives
+and pure master variables.
+We extract from the master 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
+- `dw_subprob_c` contains the perenial cost of DW subproblem representative variables
+- `dw_subprob_A` is a submatrix of the master coefficient matrix that involves only DW subproblem
   representative variables.
+We also extract from the master the information we need to compute the reduced cost of pure
+master variables:
+- `pure_master_c` contains the perenial cost of pure master variables
+- `pure_master_A` is a submatrix of the master coefficient matrix that involves only pure master
+  variables.
 
 Calculation is `c - transpose(A) * master_lp_dual_solution`.
 
 This information is given to the generic implementation of the column generation algorithm
 through methods:
-- ColGen.get_orig_costs 
+- ColGen.get_subprob_var_orig_costs 
 - ColGen.get_orig_coefmatrix
 """
 struct ReducedCostsCalculationHelper
-    c::SparseVector{Float64,VarId}
-    A::DynamicSparseMatrix{ConstrId,VarId,Float64}
+    dw_subprob_c::SparseVector{Float64,VarId}
+    dw_subprob_A::DynamicSparseMatrix{ConstrId,VarId,Float64}
+    master_c::SparseVector{Float64,VarId}
+    master_A::DynamicSparseMatrix{ConstrId,VarId,Float64}
 end
 
-function ReducedCostsCalculationHelper(master)
-    dwspvar_ids = VarId[]
+"""
+Function `var_duty_func(var_id)` returns `true` if we want to keep the variable `var_id`; `false` otherwise.
+Same for `constr_duty_func(constr_id)`.
+"""
+function _get_costs_and_coeffs(master, var_duty_func, constr_duty_func)
+    var_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)
+        if iscuractive(master, var_id) && var_duty_func(var_id)
+            push!(var_ids, var_id)
             push!(peren_costs, getperencost(master, var_id))
         end
     end
 
-    dwsprep_costs = sparsevec(dwspvar_ids, peren_costs, Coluna.MAX_NB_ELEMS)
-    dwsprep_coefmatrix = _submatrix(
+    costs = sparsevec(var_ids, peren_costs, Coluna.MAX_NB_ELEMS)
+    coef_matrix = _submatrix(master, constr_duty_func, var_duty_func)
+    return costs, coef_matrix 
+end
+
+function ReducedCostsCalculationHelper(master)
+    dw_subprob_c, dw_subprob_A = _get_costs_and_coeffs(
         master, 
-        constr_id -> !(getduty(constr_id) <= MasterConvexityConstr),
-        var_id -> getduty(var_id) <= AbstractMasterRepDwSpVar
+        var_id -> getduty(var_id) <= AbstractMasterRepDwSpVar,
+        constr_id -> !(getduty(constr_id) <= MasterConvexityConstr)
+    )
+
+    master_c, master_A = _get_costs_and_coeffs(
+        master, 
+        var_id -> getduty(var_id) <= AbstractOriginMasterVar,
+        constr_id -> !(getduty(constr_id) <= MasterConvexityConstr)
     )
-    return ReducedCostsCalculationHelper(dwsprep_costs, dwsprep_coefmatrix)
+
+    return ReducedCostsCalculationHelper(dw_subprob_c, dw_subprob_A, master_c, master_A)
 end
 
 """

src/ColGen/interface.jl

@@ -149,13 +149,13 @@ LP solution is integer feasible; `nothing` otherwise.
 Returns the original cost `c` of subproblems variables.
 to compute reduced cost `̄c = c - transpose(A) * π`.
 """
-@mustimplement "ColGenReducedCosts " get_orig_costs(ctx::AbstractColGenContext)
+@mustimplement "ColGenReducedCosts " get_subprob_var_orig_costs(ctx::AbstractColGenContext)
 
 """
 Returns the coefficient matrix `A` of subproblem variables in the master
 to compute reduced cost `̄c = c - transpose(A) * π`.
 """
-@mustimplement "ColGenReducedCosts" get_coef_matrix(ctx::AbstractColGenContext)
+@mustimplement "ColGenReducedCosts" get_subprob_var_coef_matrix(ctx::AbstractColGenContext)
 
 "Updates reduced costs of variables of a given subproblem."
 @mustimplement "ColGenReducedCosts" update_sp_vars_red_costs!(ctx::AbstractColGenContext, sp, red_costs)
@@ -184,6 +184,8 @@ end
 
 function compute_dual_bound(ctx, phase, master_lp_obj_val, master_dbs)
     # TODO pure master variables are missing.
+    @show master_lp_obj_val
+    @show master_dbs
     return master_lp_obj_val + mapreduce(((id, val),) -> val, +, master_dbs)
 end
 
@@ -242,15 +244,14 @@ function run_colgen_iteration!(context, phase, env)
     # stabcenter is master_dual_sol
     # return alpha * stab_center + (1 - alpha) * lp_dual_sol
 
-
     # With stabilization, you solve several times the suproblem because you can have misprice
     # loop:
     #   - solve all subproblems 
     #   - check if misprice 
 
     # Compute reduced cost (generic operation) by you must support math operations.
-    c = get_orig_costs(context)
-    A = get_coef_matrix(context)
+    c = get_subprob_var_orig_costs(context)
+    A = get_subprob_var_coef_matrix(context)
     red_costs = c - transpose(A) * mast_dual_sol
 
     # Updates subproblems reduced costs.

test/unit/ColGen/colgen_default.jl

@@ -2,11 +2,12 @@
 form1() = """
 master
     min
-    3x1 + 2x2 + 5x3 + 4y1 + 3y2 + 5y3
+    3x1 + 2x2 + 5x3 + 4y1 + 3y2 + 5y3 + z
     s.t.
-    x1 + x2 + x3 + y1 + y2 + y3  >= 10
-    x1 + 2x2     + y1 + 2y2      <= 100
-    x1 +     3x3 + y1 +    + 3y3 == 100
+    x1 + x2 + x3 + y1 + y2 + y3 + 2z >= 10
+    x1 + 2x2     + y1 + 2y2     + z <= 100
+    x1 +     3x3 + y1 +    + 3y3    == 100
+                                  z <= 5   
 
 dw_sp
     min
@@ -17,6 +18,9 @@ dw_sp
     integer
         representatives
             x1, x2, x3, y1, y2, y3
+
+        pure
+            z
 
     bounds
         x1 >= 0
@@ -25,6 +29,7 @@ dw_sp
         y1 >= 0
         y2 >= 0
         y3 >= 0
+        z >= 0
 """
 
 function get_name_to_varids(form)
@@ -46,33 +51,52 @@ end
 # Simple case with only subproblem representatives variables.
 function test_reduced_costs_calculation_helper()
     _, master, _, _, _ = reformfromstring(form1())
+    @show master
     vids = get_name_to_varids(master)
     cids = get_name_to_constrids(master)
 
     helper = ClA.ReducedCostsCalculationHelper(master)
-    @test helper.c[vids["x1"]] == 3
-    @test helper.c[vids["x2"]] == 2
-    @test helper.c[vids["x3"]] == 5
-    @test helper.c[vids["y1"]] == 4
-    @test helper.c[vids["y2"]] == 3
-    @test helper.c[vids["y3"]] == 5
-
-    @test helper.A[cids["c1"], vids["x1"]] == 1
-    @test helper.A[cids["c1"], vids["x2"]] == 1
-    @test helper.A[cids["c1"], vids["x3"]] == 1
-    @test helper.A[cids["c1"], vids["y1"]] == 1
-    @test helper.A[cids["c1"], vids["y2"]] == 1
-    @test helper.A[cids["c1"], vids["y3"]] == 1
-
-    @test helper.A[cids["c2"], vids["x1"]] == 1
-    @test helper.A[cids["c2"], vids["x2"]] == 2
-    @test helper.A[cids["c2"], vids["y1"]] == 1
-    @test helper.A[cids["c2"], vids["y2"]] == 2
-
-    @test helper.A[cids["c3"], vids["x1"]] == 1
-    @test helper.A[cids["c3"], vids["x3"]] == 3
-    @test helper.A[cids["c3"], vids["y1"]] == 1
-    @test helper.A[cids["c3"], vids["y3"]] == 3
+    @test helper.dw_subprob_c[vids["x1"]] == 3
+    @test helper.dw_subprob_c[vids["x2"]] == 2
+    @test helper.dw_subprob_c[vids["x3"]] == 5
+    @test helper.dw_subprob_c[vids["y1"]] == 4
+    @test helper.dw_subprob_c[vids["y2"]] == 3
+    @test helper.dw_subprob_c[vids["y3"]] == 5
+    @test helper.dw_subprob_c[vids["z"]] == 0
+
+    @test helper.dw_subprob_A[cids["c1"], vids["x1"]] == 1
+    @test helper.dw_subprob_A[cids["c1"], vids["x2"]] == 1
+    @test helper.dw_subprob_A[cids["c1"], vids["x3"]] == 1
+    @test helper.dw_subprob_A[cids["c1"], vids["y1"]] == 1
+    @test helper.dw_subprob_A[cids["c1"], vids["y2"]] == 1
+    @test helper.dw_subprob_A[cids["c1"], vids["y3"]] == 1
+    @test helper.dw_subprob_A[cids["c1"], vids["z"]] == 0  # z is not in the subproblem.
+
+    @test helper.dw_subprob_A[cids["c2"], vids["x1"]] == 1
+    @test helper.dw_subprob_A[cids["c2"], vids["x2"]] == 2
+    @test helper.dw_subprob_A[cids["c2"], vids["y1"]] == 1
+    @test helper.dw_subprob_A[cids["c2"], vids["y2"]] == 2
+    @test helper.dw_subprob_A[cids["c2"], vids["z"]] == 0 # z is not in the subproblem.
+
+    @test helper.dw_subprob_A[cids["c3"], vids["x1"]] == 1
+    @test helper.dw_subprob_A[cids["c3"], vids["x3"]] == 3
+    @test helper.dw_subprob_A[cids["c3"], vids["y1"]] == 1
+    @test helper.dw_subprob_A[cids["c3"], vids["y3"]] == 3
+    @test helper.dw_subprob_A[cids["c3"], vids["z"]] == 0 # z is not in the subproblem.
+
+    @test helper.master_c[vids["x1"]] == 0 # x1 is not in the master.
+    @test helper.master_c[vids["x2"]] == 0 # x2 is not in the master.
+    @test helper.master_c[vids["x3"]] == 0 # x3 is not in the master.
+    @test helper.master_c[vids["y1"]] == 0 # y1 is not in the master.
+    @test helper.master_c[vids["y2"]] == 0 # y2 is not in the master.
+    @test helper.master_c[vids["y3"]] == 0 # y3 is not in the master.
+    @test helper.master_c[vids["z"]] == 1
+
+    @test helper.master_A[cids["c1"], vids["x1"]] == 0 # x1 is not in the master.
+    @test helper.master_A[cids["c1"], vids["z"]] == 2
+    @test helper.master_A[cids["c2"], vids["z"]] == 1
+    @test helper.master_A[cids["c3"], vids["z"]] == 0
+    @test helper.master_A[cids["c4"], vids["z"]] == 1
 end
 register!(unit_tests, "colgen_default", test_reduced_costs_calculation_helper)
 
@@ -486,7 +510,6 @@ function ColGen.optimize_master_lp_problem!(master, ctx::TestColGenIterationCont
     end
 
     dual_sol = ColGen.get_dual_sol(output)
-    @show dual_sol
     for (constr_id, constr) in ClMP.getconstrs(master)
         name = ClMP.getname(master, constr)
         if !haskey(ctx.master_lp_dual_sol, name)
@@ -501,17 +524,13 @@ end
 ColGen.is_unbounded(ctx::TestColGenIterationContext) = ColGen.is_unbounded(ctx.context)
 ColGen.is_infeasible(ctx::TestColGenIterationContext) = ColGen.is_infeasible(ctx.context)
 ColGen.update_master_constrs_dual_vals!(ctx::TestColGenIterationContext, phase, reform, master_lp_dual_sol) = ColGen.update_master_constrs_dual_vals!(ctx.context, phase, reform, master_lp_dual_sol)
-ColGen.get_orig_costs(ctx::TestColGenIterationContext) = ColGen.get_orig_costs(ctx.context)
-ColGen.get_coef_matrix(ctx::TestColGenIterationContext) = ColGen.get_coef_matrix(ctx.context)
+ColGen.get_subprob_var_orig_costs(ctx::TestColGenIterationContext) = ColGen.get_subprob_var_orig_costs(ctx.context)
+ColGen.get_subprob_var_coef_matrix(ctx::TestColGenIterationContext) = ColGen.get_subprob_var_coef_matrix(ctx.context)
 
 function ColGen.update_sp_vars_red_costs!(ctx::TestColGenIterationContext, sp::Formulation{DwSp}, red_costs)
-    for i in 1:5
-        println("\e[34m ***************** \e[00m")
-    end
     ColGen.update_sp_vars_red_costs!(ctx.context, sp, red_costs)
     for (_, var) in ClMP.getvars(sp)
         name = ClMP.getname(sp, var)
-        println(" ---- name = $(name) ---- expected : $(ctx.pricing_var_reduced_costs[name]) ---- actual : $(ClMP.getcurcost(sp, var)) --- cur_cost = $(ClMP.getcurcost(sp, var)))")
         @test ctx.pricing_var_reduced_costs[name] ≈ ClMP.getcurcost(sp, var)
     end
     return
@@ -597,7 +616,7 @@ function test_colgen_iteration_min_gap()
     @test output.infeasible_subproblem == false
     @test output.unbounded_subproblem == false
 end
-#register!(unit_tests, "colgen_default", test_colgen_iteration_min_gap)
+register!(unit_tests, "colgen_default", test_colgen_iteration_min_gap)
 
 function test_colgen_iteration_max_gap()
     env, master, sps, reform = max_toy_gap()
@@ -659,10 +678,8 @@ function test_colgen_iteration_max_gap()
     @test output.unbounded_master == false
     @test output.infeasible_subproblem == false
     @test output.unbounded_subproblem == false
-
-
 end
-#register!(unit_tests, "colgen_default", test_colgen_iteration_max_gap)
+register!(unit_tests, "colgen_default", test_colgen_iteration_max_gap)
 
 function test_colgen_iteration_pure_master_vars()
     env, master, sps, reform = toy_gap_with_penalties()
@@ -681,13 +698,13 @@ function test_colgen_iteration_pure_master_vars()
         "c6" => 15.41666667,
         "c7" => 19.26666667,  # fixed
         "c8" => -10.86666667,
-        "c10" =>  -22.55,
+        "c10" => -22.55,
         "c12" => -30.83333334
     )
     master_obj_val = 52.95 
 
     pricing_var_reduced_costs = Dict(
-        "x_11" => 0.26666666999999933,
+        "x_11" => - 0.26666666999999933,
         "x_12" => - 12.13333333,
         "x_13" => - 7.56666667,
         "x_14" => 0.0,

test/unit/ColGen/colgen_iteration.jl

@@ -153,8 +153,8 @@ function ColGen.optimize_pricing_problem!(ctx::ColGenIterationTestContext, form,
 end
 
 # Reduced costs
-ColGen.get_orig_costs(::ColGenIterationTestContext) = [7, 2, 1, 5, 3, 6, 8, 4, 2, 9]
-ColGen.get_coef_matrix(::ColGenIterationTestContext) = [
+ColGen.get_subprob_var_orig_costs(::ColGenIterationTestContext) = [7, 2, 1, 5, 3, 6, 8, 4, 2, 9]
+ColGen.get_subprob_var_coef_matrix(::ColGenIterationTestContext) = [
     1 1 1 1 0 0 0 0 0 0; 
     1 0 0 0 1 1 1 0 0 0; 
     0 1 0 0 1 0 0 1 1 0;