Projection and mapping for identical dw subproblem

src/Algorithm/basic/cutcallback.jl

@@ -40,7 +40,7 @@ function run!(algo::CutCallbacks, env::Env, form::Formulation, input::CutCallbac
     if length(robust_generators) > 0 && (algo.call_robust_facultative || algo.call_robust_essential)
         !projection_is_possible(form) && error("Cannot do projection on original variables. Open an issue.")
 
-        projsol1 = proj_cols_on_rep(input.primalsol, form)
+        projsol1 = proj_cols_on_rep(input.primalsol)
         projsol2 = Dict{VarId, Float64}(varid => val for (varid, val) in projsol1)
         viol_vals = Float64[]
 

src/Algorithm/branching/interface.jl

@@ -16,7 +16,7 @@ function get_original_sol(::Branching.AbstractDivideContext, reform, opt_state)
     original_sol = nothing
     if !isnothing(extended_sol)
         original_sol = if projection_is_possible(master)
-            proj_cols_on_rep(extended_sol, master)
+            proj_cols_on_rep(extended_sol)
         else
             get_best_lp_primal_sol(opt_state) # it means original_sol equals extended_sol(requires discussion)
         end

src/Algorithm/colgen.jl

@@ -558,7 +558,7 @@ end
 
 function _is_feasible_and_integer(rm_lp_primal_sol)
     return !contains(rm_lp_primal_sol, varid -> isanArtificialDuty(getduty(varid))) &&
-        isinteger(proj_cols_on_rep(rm_lp_primal_sol, getmodel(rm_lp_primal_sol)))
+        isinteger(proj_cols_on_rep(rm_lp_primal_sol))
 end
 
 function _assert_has_lp_dual_sol(rm_optstate, master, rm_optimizer_id)

src/Algorithm/colgen/default.jl

@@ -17,6 +17,8 @@ mutable struct ColGenContext <: ColGen.AbstractColGenContext
     opt_rtol::Float64
     opt_atol::Float64
 
+    incumbent_primal_solution::PrimalSolution
+
     # # Information to solve the master
     # master_solve_alg
     # master_optimizer_id
@@ -216,8 +218,17 @@ function ColGen.update_master_constrs_dual_vals!(ctx::ColGenContext, phase, refo
 
 end
 
-function ColGen.check_primal_ip_feasibility(master_lp_primal_sol, phase, reform)
+function ColGen.check_primal_ip_feasibility(master_lp_primal_sol, ::ColGenContext, phase, reform)
+    primal_sol_is_integer = MathProg.proj_cols_is_integer(master_lp_primal_sol)
+    if !primal_sol_is_integer
+        return nothing
+    end
+    return MathProg.proj_cols_on_rep(master_lp_primal_sol)
+end
 
+function ColGen.update_inc_primal_sol!(ctx::ColGenContext, ip_primal_sol)
+    ctx.incumbent_primal_solution = ip_primal_sol
+    return
 end
 
 # Reduced costs calculation

src/Algorithm/colgen/printer.jl

@@ -36,6 +36,9 @@ function ColGen.update_master_constrs_dual_vals!(ctx::ColGenPrinterContext, phas
     return ColGen.update_master_constrs_dual_vals!(ctx.inner, phase, reform, master_lp_dual_sol)
 end
 
+ColGen.check_primal_ip_feasibility(mast_primal_sol, ctx::ColGenPrinterContext, phase, reform) = ColGen.check_primal_ip_feasibility(mast_primal_sol, ctx.inner, phase, reform)
+ColGen.update_inc_primal_sol!(ctx::ColGenPrinterContext, ip_primal_sol) = ColGen.update_inc_primal_sol!(ctx.inner, ip_primal_sol)
+
 ColGen.get_subprob_var_orig_costs(ctx::ColGenPrinterContext) = ColGen.get_subprob_var_orig_costs(ctx.inner)
 ColGen.get_subprob_var_coef_matrix(ctx::ColGenPrinterContext) = ColGen.get_subprob_var_coef_matrix(ctx.inner)
 

src/ColGen/interface.jl

@@ -169,7 +169,9 @@ information at two different places.
 Returns a primal solution expressed in the original problem variables if the current master
 LP solution is integer feasible; `nothing` otherwise.
 """
-@mustimplement "ColGenMaster" check_primal_ip_feasibility(mast_lp_primal_sol, phase, reform) = nothing
+@mustimplement "ColGenMaster" check_primal_ip_feasibility(mast_lp_primal_sol, ::AbstractColGenContext, phase, reform) = nothing
+
+@mustimplement "ColGen" update_inc_primal_sol!(ctx::AbstractColGenContext, ip_primal_sol) = nothing
 
 ############################################################################################
 # Reduced costs calculation.
@@ -262,7 +264,7 @@ function run_colgen_iteration!(context, phase, env)
     if !isnothing(mast_primal_sol)
         # If the master LP problem has a primal solution, we can try to find a integer feasible
         # solution.
-        ip_primal_sol = check_primal_ip_feasibility(phase, mast_primal_sol, get_reform(context))
+        ip_primal_sol = check_primal_ip_feasibility(mast_primal_sol, context, phase, get_reform(context))
         if !isnothing(ip_primal_sol)
             update_inc_primal_sol!(context, ip_primal_sol)
         end

src/MathProg/projection.jl

@@ -1,19 +1,124 @@
+"Returns `true` if we can project a solution of `form` to the original formulation."
+projection_is_possible(form) = false
+
+############################################################################################
+# Projection of Dantzig-Wolfe master on original formulation.
+############################################################################################
+
 projection_is_possible(master::Formulation{DwMaster}) = true
-function proj_cols_on_rep(sol::PrimalSolution, master::Formulation{DwMaster})
-    projected_sol_vars = Vector{VarId}()
-    projected_sol_vals = Vector{Float64}()
+
+Base.isless(A::DynamicMatrixColView{VarId, VarId, Float64}, B::DynamicMatrixColView{VarId, VarId, Float64}) = cmp(A, B) < 0
+
+function Base.cmp(A::DynamicMatrixColView{VarId, VarId, Float64}, B::DynamicMatrixColView{VarId, VarId, Float64})
+    for (a, b) in zip(A, B)
+        if !isequal(a, b)
+            return isless(a, b) ? -1 : 1
+        end
+    end
+    return 0 # no length for dynamic sparse vectors
+end
+
+function _assign_width!(cur_roll, col::Vector, width_to_assign)
+    for i in 1:length(col)
+        cur_roll[i] += col[i] * width_to_assign
+    end
+    return
+end
+
+function _assign_width!(cur_roll::Dict, col::DynamicMatrixColView, width_to_assign)
+    for (id, val) in col
+        if !haskey(cur_roll, id)
+            cur_roll[id] = 0.0
+        end
+        cur_roll[id] += val * width_to_assign
+    end
+    return
+end
+
+_new_set_of_rolls(::Type{Vector{E}}) where {E} = Vector{Float64}[]
+_new_roll(::Type{Vector{E}}, col_len) where {E} = zeros(Float64, col_len)
+_roll_is_integer(roll::Vector{Float64}) = all(isinteger.(roll))
+
+_new_set_of_rolls(::Type{DynamicMatrixColView{VarId, VarId, Float64}}) = Dict{VarId, Float64}[]
+_new_roll(::Type{DynamicMatrixColView{VarId, VarId, Float64}}, _) = Dict{VarId, Float64}()
+_roll_is_integer(roll::Dict{VarId, Float64}) = all(isinteger.(values(roll)))
+
+function _mapping(columns::Vector{A}, values::Vector{B}; col_len::Int = 10) where {A,B}
+    p = sortperm(columns, rev=true)
+    columns = columns[p]
+    values = values[p]
+
+    rolls = _new_set_of_rolls(eltype(columns))
+    total_width_assigned = 0 
+    nb_roll_opened = 1 # roll is width 1
+    cur_roll = _new_roll(eltype(columns), col_len)
+
+    for (val, col) in zip(values, columns)
+        cur_unassigned_width = val
+        while cur_unassigned_width > 0
+            width_to_assign = min(cur_unassigned_width, nb_roll_opened - total_width_assigned)
+            _assign_width!(cur_roll, col, width_to_assign)
+            cur_unassigned_width -= width_to_assign
+            total_width_assigned += width_to_assign
+            if total_width_assigned == nb_roll_opened
+                push!(rolls, cur_roll)
+                cur_roll = _new_roll(eltype(columns), col_len)
+                nb_roll_opened += 1
+            end
+        end
+    end
+    return rolls
+end
+
+function _mapping_by_subproblem(columns::Dict{Int, Vector{A}}, values::Dict{Int, Vector{B}}) where {A,B}
+    return Dict(
+        uid =>  _mapping(cols, values[uid]) for (uid, cols) in columns
+    )
+end
+
+_rolls_are_integer(rolls) = all(_roll_is_integer.(rolls))
+_subproblem_rolls_are_integer(rolls_by_sp::Dict) = all(_rolls_are_integer.(values(rolls_by_sp)))
+
+function _extract_data_for_mapping(sol::PrimalSolution{Formulation{DwMaster}})
+    columns = Dict{Int, Vector{DynamicMatrixColView{VarId, VarId, Float64}}}()
+    values = Dict{Int, Vector{Float64}}()
+    master = getmodel(sol)
+
+    for (varid, val) in sol
+        duty = getduty(varid)
+        if duty <= MasterCol
+            origin_form_uid = getoriginformuid(varid)
+            spform = get_dw_pricing_sps(master.parent_formulation)[origin_form_uid]
+            column = @view get_primal_sol_pool(spform)[varid,:]
+            if !haskey(columns, origin_form_uid)
+                columns[origin_form_uid] = DynamicMatrixColView{VarId, VarId, Float64}[]
+                values[origin_form_uid] = Float64[]
+            end
+            push!(columns[origin_form_uid], column)
+            push!(values[origin_form_uid], val)
+        end
+    end
+    return columns, values
+end
+
+function _proj_cols_on_rep(sol::PrimalSolution{Formulation{DwMaster}}, extracted_cols, extracted_vals)
+    projected_sol_vars = VarId[]
+    projected_sol_vals = Float64[]
 
     for (varid, val) in sol
         duty = getduty(varid)
         if duty <= MasterPureVar
             push!(projected_sol_vars, varid)
             push!(projected_sol_vals, val)
-        elseif duty <= MasterCol
-            origin_form_uid = getoriginformuid(varid)
-            spform = get_dw_pricing_sps(master.parent_formulation)[origin_form_uid]
+        end
+    end
 
-            for (repid, repval) in @view get_primal_sol_pool(spform)[varid,:]
-                if getduty(repid) <= DwSpPricingVar || getduty(repid) <= DwSpSetupVar
+    master = getmodel(sol)
+    for spid in keys(extracted_cols)
+        for (column, val) in Iterators.zip(extracted_cols[spid], extracted_vals[spid])
+            for (repid, repval) in column
+                if getduty(repid) <= DwSpPricingVar || getduty(repid) <= DwSpSetupVar || 
+                    getduty(repid) <= MasterRepPricingVar || getduty(repid) <= MasterRepPricingSetupVar
                     mastrepid = getid(getvar(master, repid))
                     push!(projected_sol_vars, mastrepid)
                     push!(projected_sol_vals, repval * val)
@@ -24,6 +129,24 @@ function proj_cols_on_rep(sol::PrimalSolution, master::Formulation{DwMaster})
     return PrimalSolution(master, projected_sol_vars, projected_sol_vals, getvalue(sol), FEASIBLE_SOL)
 end
 
+function proj_cols_on_rep(sol::PrimalSolution{Formulation{DwMaster}})
+    columns, values = _extract_data_for_mapping(sol)
+    projected_sol = _proj_cols_on_rep(sol, columns, values)
+    return projected_sol
+end
+
+function proj_cols_is_integer(sol::PrimalSolution{Formulation{DwMaster}})
+    columns, values = _extract_data_for_mapping(sol)
+    projected_sol = _proj_cols_on_rep(sol, columns, values)
+    rolls = _mapping_by_subproblem(columns, values)
+    integer_rolls = _subproblem_rolls_are_integer(rolls)
+    return isinteger(projected_sol) && integer_rolls
+end
+
+############################################################################################
+# Porjection of Benders master on original formulation.
+############################################################################################
+
 projection_is_possible(master::Formulation{BendersMaster}) = false
 
 function proj_cols_on_rep(sol::PrimalSolution, master::Formulation{BendersMaster})

src/optimize.jl

@@ -102,14 +102,14 @@ function optimize!(
     ip_primal_sols = get_ip_primal_sols(algstate)
     if !isnothing(ip_primal_sols)
         for sol in ip_primal_sols
-            add_ip_primal_sol!(outstate, proj_cols_on_rep(sol, master))
+            add_ip_primal_sol!(outstate, proj_cols_on_rep(sol))
         end
     end
 
     # lp_primal_sol may also be of interest, for example when solving the relaxation
     lp_primal_sol = get_best_lp_primal_sol(algstate)
     if !isnothing(lp_primal_sol)
-        add_lp_primal_sol!(outstate, proj_cols_on_rep(lp_primal_sol, master))
+        add_lp_primal_sol!(outstate, proj_cols_on_rep(lp_primal_sol))
     end
 
     # lp_dual_sol to retrieve, for instance, the dual value of generated cuts

test/unit/ColGen/colgen_default.jl

@@ -674,6 +674,7 @@ function ColGen.optimize_master_lp_problem!(master, ctx::TestColGenIterationCont
     return output
 end
 
+ColGen.check_primal_ip_feasibility(master_lp_primal_sol, ::TestColGenIterationContext, phase, reform) = nothing
 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)

test/unit/ColGen/colgen_iteration.jl

@@ -174,7 +174,7 @@ function ColGen.update_sp_vars_red_costs!(::ColGenIterationTestContext, subprob,
     return
 end
 
-ColGen.check_primal_ip_feasibility(::ColGenIterationTestPhase, sol, reform) = nothing
+ColGen.check_primal_ip_feasibility(sol, ::ColGenIterationTestContext, ::ColGenIterationTestPhase, reform) = nothing
 ColGen.update_master_constrs_dual_vals!(::ColGenIterationTestContext, ::ColGenIterationTestPhase, reform, dual_mast_sol) = nothing
 
 function ColGen.insert_columns!(reform, ::ColGenIterationTestContext, phase, generated_columns)