From e7a5ac4c20cc286c53ba51d2cde425ae1e85e6f6 Mon Sep 17 00:00:00 2001 From: Guillaume Marques Date: Fri, 7 Apr 2023 15:53:58 +0200 Subject: [PATCH 1/5] implementing projection and mapping --- src/MathProg/projection.jl | 132 +++++++++++++++++++++--- test/unit/MathProg/projection.jl | 166 +++++++++++++++++++++++++++++++ test/unit/MathProg/solutions.jl | 2 +- 3 files changed, 287 insertions(+), 13 deletions(-) create mode 100644 test/unit/MathProg/projection.jl diff --git a/src/MathProg/projection.jl b/src/MathProg/projection.jl index 528cdc53d..faa2c13e9 100644 --- a/src/MathProg/projection.jl +++ b/src/MathProg/projection.jl @@ -1,29 +1,137 @@ +"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) + +_new_set_of_rolls(::Type{DynamicMatrixColView{VarId, VarId, Float64}}) = Dict{VarId, Float64}[] +_new_roll(::Type{DynamicMatrixColView{VarId, VarId, Float64}}, col_len) = Dict{VarId, Float64}() + +function _mapping(columns::Vector{A}, values::Vector{B}, col_len::Int) 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 _extract_data_for_mapping(sol::PrimalSolution{Formulation{DwMaster}}) + columns = DynamicMatrixColView{VarId, VarId, Float64}[] + values = 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,:] + push!(columns, column) + push!(values, val) + end + end + + col_len = 0 + for (var_id, _) in getvars(master) + if getduty(var_id) <= DwSpPricingVar || getduty(var_id) <= DwSpSetupVar + col_len += 1 + end + end + return columns, values, col_len +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 - mastrepid = getid(getvar(master, repid)) - push!(projected_sol_vars, mastrepid) - push!(projected_sol_vals, repval * val) - end + master = getmodel(sol) + for (column, val) in Iterators.zip(extracted_cols, extracted_vals) + for (repid, repval) in column + if getduty(repid) <= DwSpPricingVar || getduty(repid) <= DwSpSetupVar + mastrepid = getid(getvar(master, repid)) + push!(projected_sol_vars, mastrepid) + push!(projected_sol_vals, repval * val) end end end 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, col_len = _extract_data_for_mapping(sol) + projected_sol = _proj_cols_on_rep(sol, columns, values) + println("\e[34m ~~~~~world starts here ~~~~~~~ \e[00m]") + rolls = _mapping(columns, values, col_len) + @show rolls + return projected_sol +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}) diff --git a/test/unit/MathProg/projection.jl b/test/unit/MathProg/projection.jl new file mode 100644 index 000000000..5df8420ca --- /dev/null +++ b/test/unit/MathProg/projection.jl @@ -0,0 +1,166 @@ +function test_mapping_operator_1() + G = Vector{Float64}[ + [0, 0, 1, 1, 0, 0, 1], + [1, 0, 0, 1, 1, 0, 1], + [1, 1, 0, 0, 1, 1, 0], + [0, 1, 1, 0, 0, 1, 1] + ] + + v = Float64[3, 2, 0, 0] + + result = Coluna.MathProg._mapping(G, v, 7)#, 6) + @test result == [ + [1, 0, 0, 1, 1, 0, 1], + [1, 0, 0, 1, 1, 0, 1], + [0, 0, 1, 1, 0, 0, 1], + [0, 0, 1, 1, 0, 0, 1], + [0, 0, 1, 1, 0, 0, 1], + ] +end +register!(unit_tests, "projection", test_mapping_operator_1; f= true) + +function test_mapping_operator_2() + # Example from the paper: + # Branching in Branch-and-Price: a Generic Scheme + # François Vanderbeck + + # A = [ + # 1 0 1 0; + # 0 1 0 1; + # 1 1 0 2; + # ] + # a = [5 5 10] + + G = Vector{Float64}[ + [1, 1, 1, 0], + [1, 1, 0, 1], + [1, 1, 0, 0], + [1, 0, 1, 1], + [1, 0, 1, 0], + [1, 0, 0, 1], + [1, 0, 0, 0], + [0, 1, 1, 1], + [0, 1, 1, 0], + [0, 1, 0, 1], + [0, 1, 0, 0], + [0, 0, 1, 1], + [0, 0, 1, 0], + [0, 0, 0, 1], + ] + + v = Float64[0, 1/2, 1, 1/2, 0, 0, 1, 1, 0, 0, 1/2, 0, 1/2, 0, 0] + + result = Coluna.MathProg._mapping(G, v, 4)#, 5) + @test result == [ + [1, 1, 0, 1/2], + [1, 1/2, 1/2, 1/2], + [1, 0, 0, 0], + [0, 1, 1, 1], + [0, 1/2, 1/2, 0] + ] +end +register!(unit_tests, "projection", test_mapping_operator_2; f=true) + +function identical_subproblems_vrp() + # We consider a vrp problem (with fake subproblem) where routes are: + # - MC1 : 1 -> 2 -> 3 + # - MC2 : 2 -> 3 -> 4 + # - MC4 : 3 -> 4 -> 1 + # - MC3 : 4 -> 1 -> 3 + # At three vehicles are available to visit all customers. + form = """ + master + min + x_12 + x_13 + x_14 + x_23 + x_24 + x_34 + MC1 + MC2 + MC3 + MC4 + 0.0 PricingSetupVar_sp_5 + s.t. + x_12 + x_13 + x_14 + MC1 + MC3 + MC4 >= 1.0 + x_12 + x_23 + x_24 + MC1 + MC2 + MC4 >= 1.0 + x_13 + x_23 + x_34 + MC1 + MC2 + MC3 >= 1.0 + x_14 + x_24 + x_34 + MC2 + MC3 + MC4 >= 1.0 + PricingSetupVar_sp_5 >= 0.0 {MasterConvexityConstr} + PricingSetupVar_sp_5 <= 3.0 {MasterConvexityConstr} + + dw_sp + min + x_12 + x_13 + x_14 + x_23 + x_24 + x_34 + 0.0 PricingSetupVar_sp_5 + s.t. + x_12 + x_13 + x_14 + x_23 + x_24 + x_34 >= 0 + + continuous + columns + MC1, MC2, MC3, MC4 + + integer + pricing_setup + PricingSetupVar_sp_5 + + binary + representatives + x_12, x_13, x_14, x_23, x_24, x_34 + + bounds + 0.0 <= x_12 <= 1.0 + 0.0 <= x_13 <= 1.0 + 0.0 <= x_14 <= 1.0 + 0.0 <= x_23 <= 1.0 + 0.0 <= x_24 <= 1.0 + 0.0 <= x_34 <= 1.0 + MC1 >= 0 + MC2 >= 0 + MC3 >= 0 + MC4 >= 0 + 1.0 <= PricingSetupVar_sp_5 <= 1.0 + """ + env, master, sps, _, reform = reformfromstring(form) + return env, master, sps, reform +end + +function projection_from_dw_reform_to_master_1() + env, master, sps, reform = identical_subproblems_vrp() + mastervarids = Dict(CL.getname(master, var) => varid for (varid, var) in CL.getvars(master)) + + # Register column in the pool + spform = first(sps) + pool = ClMP.get_primal_sol_pool(spform) + pool_hashtable = ClMP._get_primal_sol_pool_hash_table(spform) + costs_pool = spform.duty_data.costs_primalsols_pool + custom_pool = spform.duty_data.custom_primalsols_pool + + var_ids = map(n -> mastervarids[n], ["x_12", "x_13", "x_14", "x_23", "x_24", "x_34"]) + for (name, vals) in Iterators.zip( + ["MC1", "MC2", "MC3", "MC4"], + [ + [1.0, 0.0, 0.0, 1.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 1.0, 0.0, 1.0], + [0.0, 0.0, 1.0, 0.0, 0.0, 1.0], + [1.0, 0.0, 1.0, 0.0, 0.0, 0.0] + ] + ) + col_id = ClMP.VarId(mastervarids[name]; duty = DwSpPrimalSol) + addrow!(pool, col_id, var_ids, vals) + costs_pool[col_id] = 1.0 + ClMP.savesolid!(pool_hashtable, col_id, ClMP.PrimalSolution(spform, var_ids, vals, 1.0, ClMP.FEASIBLE_SOL)) + end + + # Create primal solution where each route is used 1/2 time. + # This solution is integer feasible. + solution = Coluna.MathProg.PrimalSolution( + master, + map(n -> ClMP.VarId(mastervarids[n]; origin_form_uid = 2), ["MC1", "MC2", "MC3", "MC4"]), + [1/2, 1/2, 1/2, 1/2], + 2.0, + ClB.FEASIBLE_SOL + ) + + @show mastervarids["MC1"].uid + @show mastervarids["MC1"].origin_form_uid + @show mastervarids["MC1"].assigned_form_uid + + proj_cols_on_rep(solution) +end +register!(unit_tests, "projection", projection_from_dw_reform_to_master_1; f = true) + +function projection_from_dw_reform_to_master_2() + +end +register!(unit_tests, "projection", projection_from_dw_reform_to_master_2) \ No newline at end of file diff --git a/test/unit/MathProg/solutions.jl b/test/unit/MathProg/solutions.jl index 5fc2a9663..107ec814b 100644 --- a/test/unit/MathProg/solutions.jl +++ b/test/unit/MathProg/solutions.jl @@ -169,4 +169,4 @@ function lin_alg_3_spMv() g = transpose(dyn_mat) * vec_len5 @test e == f == g end -register!(unit_tests, "solutions", lin_alg_3_spMv) \ No newline at end of file +register!(unit_tests, "solutions", lin_alg_3_spMv) From 72847e48c92e3d640c21d44870607cbc735e078d Mon Sep 17 00:00:00 2001 From: Guillaume Marques Date: Fri, 7 Apr 2023 21:39:08 +0200 Subject: [PATCH 2/5] wip --- src/Algorithm/colgen/default.jl | 2 +- src/MathProg/projection.jl | 60 ++++++++++++++++------------ test/unit/MathProg/projection.jl | 68 ++++++++++++++++++++++---------- 3 files changed, 83 insertions(+), 47 deletions(-) diff --git a/src/Algorithm/colgen/default.jl b/src/Algorithm/colgen/default.jl index db2f9a862..17ab49768 100644 --- a/src/Algorithm/colgen/default.jl +++ b/src/Algorithm/colgen/default.jl @@ -213,7 +213,7 @@ function ColGen.update_master_constrs_dual_vals!(ctx::ColGenContext, phase, refo end function ColGen.check_primal_ip_feasibility(master_lp_primal_sol, phase, reform) - + end # Reduced costs calculation diff --git a/src/MathProg/projection.jl b/src/MathProg/projection.jl index faa2c13e9..c34934ce4 100644 --- a/src/MathProg/projection.jl +++ b/src/MathProg/projection.jl @@ -39,9 +39,9 @@ _new_set_of_rolls(::Type{Vector{E}}) where {E} = Vector{Float64}[] _new_roll(::Type{Vector{E}}, col_len) where {E} = zeros(Float64, col_len) _new_set_of_rolls(::Type{DynamicMatrixColView{VarId, VarId, Float64}}) = Dict{VarId, Float64}[] -_new_roll(::Type{DynamicMatrixColView{VarId, VarId, Float64}}, col_len) = Dict{VarId, Float64}() +_new_roll(::Type{DynamicMatrixColView{VarId, VarId, Float64}}, _) = Dict{VarId, Float64}() -function _mapping(columns::Vector{A}, values::Vector{B}, col_len::Int) where {A,B} +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] @@ -68,30 +68,32 @@ function _mapping(columns::Vector{A}, values::Vector{B}, col_len::Int) where {A, 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 function _extract_data_for_mapping(sol::PrimalSolution{Formulation{DwMaster}}) - columns = DynamicMatrixColView{VarId, VarId, Float64}[] - values = Float64[] + 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,:] - push!(columns, column) - push!(values, val) - end - end - - col_len = 0 - for (var_id, _) in getvars(master) - if getduty(var_id) <= DwSpPricingVar || getduty(var_id) <= DwSpSetupVar - col_len += 1 + 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, col_len + return columns, values end function _proj_cols_on_rep(sol::PrimalSolution{Formulation{DwMaster}}, extracted_cols, extracted_vals) @@ -107,12 +109,16 @@ function _proj_cols_on_rep(sol::PrimalSolution{Formulation{DwMaster}}, extracted end master = getmodel(sol) - for (column, val) in Iterators.zip(extracted_cols, extracted_vals) - for (repid, repval) in column - if getduty(repid) <= DwSpPricingVar || getduty(repid) <= DwSpSetupVar - mastrepid = getid(getvar(master, repid)) - push!(projected_sol_vars, mastrepid) - push!(projected_sol_vals, repval * val) + for spid in keys(extracted_cols) + for (column, val) in Iterators.zip(extracted_cols[spid], extracted_vals[spid]) + @show column, val + for (repid, repval) in column + @show getduty(repid) <= DwSpPricingVar + if getduty(repid) <= DwSpPricingVar || getduty(repid) <= DwSpSetupVar + mastrepid = getid(getvar(master, repid)) + push!(projected_sol_vars, mastrepid) + push!(projected_sol_vals, repval * val) + end end end end @@ -120,14 +126,18 @@ function _proj_cols_on_rep(sol::PrimalSolution{Formulation{DwMaster}}, extracted end function proj_cols_on_rep(sol::PrimalSolution{Formulation{DwMaster}}) - columns, values, col_len = _extract_data_for_mapping(sol) + columns, values = _extract_data_for_mapping(sol) projected_sol = _proj_cols_on_rep(sol, columns, values) - println("\e[34m ~~~~~world starts here ~~~~~~~ \e[00m]") - rolls = _mapping(columns, values, col_len) - @show rolls 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) + return isinteger(projected_sol) && all(isinteger, rolls) +end + ############################################################################################ # Porjection of Benders master on original formulation. ############################################################################################ diff --git a/test/unit/MathProg/projection.jl b/test/unit/MathProg/projection.jl index 5df8420ca..ab65c0566 100644 --- a/test/unit/MathProg/projection.jl +++ b/test/unit/MathProg/projection.jl @@ -8,7 +8,7 @@ function test_mapping_operator_1() v = Float64[3, 2, 0, 0] - result = Coluna.MathProg._mapping(G, v, 7)#, 6) + result = Coluna.MathProg._mapping(G, v; col_len = 7)#, 6) @test result == [ [1, 0, 0, 1, 1, 0, 1], [1, 0, 0, 1, 1, 0, 1], @@ -50,24 +50,38 @@ function test_mapping_operator_2() v = Float64[0, 1/2, 1, 1/2, 0, 0, 1, 1, 0, 0, 1/2, 0, 1/2, 0, 0] - result = Coluna.MathProg._mapping(G, v, 4)#, 5) - @test result == [ - [1, 1, 0, 1/2], - [1, 1/2, 1/2, 1/2], - [1, 0, 0, 0], - [0, 1, 1, 1], - [0, 1/2, 1/2, 0] - ] + result = Coluna.MathProg._mapping(G, v; col_len = 4) end register!(unit_tests, "projection", test_mapping_operator_2; f=true) +# function test_mapping_operator_3() +# G = Vector{Float64}[ +# #x_12, x_13, x_14, x_15, x_23, x_24, x_25, x_34, x_35, x_45 +# [1, 0, 1, 0, 0, 1, 0, 0, 0, 0], +# [1, 0, 0, 1, 1, 0, 0, 0, 1, 0], +# [0, 1, 1, 0, 0, 0, 0, 1, 0, 0], +# [0, 0, 0, 2, 0, 0, 0, 0, 0, 0], +# [1, 0, 1, 0, 1, 0, 0, 1, 0, 0], +# [0, 1, 0, 1, 0, 0, 0, 0, 1, 0] +# ] + +# v = Float64[2/3, 1/3, 1/3, 2/3, 1/3, 1/3] + +# result = Coluna.MathProg._mapping(G, v, 10) +# @show result + +# end +# register!(unit_tests, "projection", test_mapping_operator_3; f= true) + function identical_subproblems_vrp() # We consider a vrp problem (with fake subproblem) where routes are: # - MC1 : 1 -> 2 -> 3 # - MC2 : 2 -> 3 -> 4 # - MC4 : 3 -> 4 -> 1 - # - MC3 : 4 -> 1 -> 3 - # At three vehicles are available to visit all customers. + # - MC3 : 4 -> 1 -> 2 + # At most, three vehicles are available to visit all customers. + # We can visit a customer multiple times. + # Fractional solution is 1/2 for all columns form = """ master min @@ -152,15 +166,27 @@ function projection_from_dw_reform_to_master_1() ClB.FEASIBLE_SOL ) - @show mastervarids["MC1"].uid - @show mastervarids["MC1"].origin_form_uid - @show mastervarids["MC1"].assigned_form_uid - - proj_cols_on_rep(solution) -end -register!(unit_tests, "projection", projection_from_dw_reform_to_master_1; f = true) - -function projection_from_dw_reform_to_master_2() + # Test integration + columns, values = Coluna.MathProg._extract_data_for_mapping(solution) + rolls = Coluna.MathProg._mapping_by_subproblem(columns, values) + + # Expected: + # | 1/2 of [1.0, 0.0, 1.0, 0.0, 0.0, 0.0] + # | 1/2 of [1.0, 0.0, 0.0, 1.0, 0.0, 0.0] + # -----> [1.0, 0.0, 0.5, 0.5, 0.0, 0.0] + # | 1/2 of [0.0, 0.0, 1.0, 0.0, 0.0, 1.0] + # | 1/2 of [0.0, 0.0, 0.0, 1.0, 0.0, 1.0] + # -----> [0.0, 0.0, 0.5, 0.5, 0.0, 1.0] + @test rolls == Dict(2 => [ + Dict(mastervarids["x_14"] => 0.5, mastervarids["x_23"] => 0.5, mastervarids["x_34"] => 1.0) + Dict(mastervarids["x_12"] => 1.0, mastervarids["x_14"] => 0.5, mastervarids["x_23"] => 0.5) + ]) + + proj = Coluna.MathProg.proj_cols_on_rep(solution) + @test proj[mastervarids["x_12"]] == 1.0 + @test proj[mastervarids["x_14"]] == 1.0 + @test proj[mastervarids["x_23"]] == 1.0 + @test proj[mastervarids["x_34"]] == 1.0 end -register!(unit_tests, "projection", projection_from_dw_reform_to_master_2) \ No newline at end of file +register!(unit_tests, "projection", projection_from_dw_reform_to_master_1; f = true) \ No newline at end of file From 612143244e5b95ff3e70d7dd9e83c9caabbfb5fb Mon Sep 17 00:00:00 2001 From: Guillaume Marques Date: Tue, 11 Apr 2023 10:24:56 +0200 Subject: [PATCH 3/5] wip --- src/MathProg/projection.jl | 7 +++---- test/unit/MathProg/projection.jl | 11 ++++++++++- 2 files changed, 13 insertions(+), 5 deletions(-) diff --git a/src/MathProg/projection.jl b/src/MathProg/projection.jl index c34934ce4..c546cbf35 100644 --- a/src/MathProg/projection.jl +++ b/src/MathProg/projection.jl @@ -111,10 +111,9 @@ function _proj_cols_on_rep(sol::PrimalSolution{Formulation{DwMaster}}, extracted master = getmodel(sol) for spid in keys(extracted_cols) for (column, val) in Iterators.zip(extracted_cols[spid], extracted_vals[spid]) - @show column, val for (repid, repval) in column - @show getduty(repid) <= DwSpPricingVar - if getduty(repid) <= DwSpPricingVar || getduty(repid) <= DwSpSetupVar + 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) @@ -135,7 +134,7 @@ 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) - return isinteger(projected_sol) && all(isinteger, rolls) + return isinteger(projected_sol) && all(isinteger, map(r -> values(r), rolls)) end ############################################################################################ diff --git a/test/unit/MathProg/projection.jl b/test/unit/MathProg/projection.jl index ab65c0566..cb22fbeed 100644 --- a/test/unit/MathProg/projection.jl +++ b/test/unit/MathProg/projection.jl @@ -51,6 +51,13 @@ function test_mapping_operator_2() v = Float64[0, 1/2, 1, 1/2, 0, 0, 1, 1, 0, 0, 1/2, 0, 1/2, 0, 0] result = Coluna.MathProg._mapping(G, v; col_len = 4) + @test result == [ + [1.0, 1.0, 0.0, 0.5], + [1.0, 0.5, 0.5, 0.5], + [1.0, 0.0, 0.0, 0.0], + [0.0, 1.0, 1.0, 1.0], + [0.0, 0.5, 0.5, 0.0] + ] end register!(unit_tests, "projection", test_mapping_operator_2; f=true) @@ -141,6 +148,7 @@ function projection_from_dw_reform_to_master_1() custom_pool = spform.duty_data.custom_primalsols_pool var_ids = map(n -> mastervarids[n], ["x_12", "x_13", "x_14", "x_23", "x_24", "x_34"]) + # VarId[Variableu2, Variableu1, Variableu3, Variableu4, Variableu5, Variableu6] for (name, vals) in Iterators.zip( ["MC1", "MC2", "MC3", "MC4"], [ @@ -188,5 +196,6 @@ function projection_from_dw_reform_to_master_1() @test proj[mastervarids["x_23"]] == 1.0 @test proj[mastervarids["x_34"]] == 1.0 + @test Coluna.MathProg.proj_cols_is_integer(solution) == false end -register!(unit_tests, "projection", projection_from_dw_reform_to_master_1; f = true) \ No newline at end of file +register!(unit_tests, "projection", projection_from_dw_reform_to_master_1; f = true) From ea783564e4fede9f0e798c2224aa348a24fd2ac6 Mon Sep 17 00:00:00 2001 From: Guillaume Marques Date: Tue, 11 Apr 2023 16:28:10 +0200 Subject: [PATCH 4/5] ok --- src/Algorithm/colgen/default.jl | 15 +++++++++++++-- src/Algorithm/colgen/printer.jl | 2 ++ src/ColGen/interface.jl | 6 ++++-- src/MathProg/projection.jl | 8 +++++++- test/unit/ColGen/colgen_default.jl | 1 + test/unit/ColGen/colgen_iteration.jl | 2 +- test/unit/MathProg/projection.jl | 10 +++++++--- 7 files changed, 35 insertions(+), 9 deletions(-) diff --git a/src/Algorithm/colgen/default.jl b/src/Algorithm/colgen/default.jl index fdca13f8c..26380e84e 100644 --- a/src/Algorithm/colgen/default.jl +++ b/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 diff --git a/src/Algorithm/colgen/printer.jl b/src/Algorithm/colgen/printer.jl index be2ff99c8..2bcbe6184 100644 --- a/src/Algorithm/colgen/printer.jl +++ b/src/Algorithm/colgen/printer.jl @@ -36,6 +36,8 @@ 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.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) diff --git a/src/ColGen/interface.jl b/src/ColGen/interface.jl index eeb920ccb..fd1161828 100644 --- a/src/ColGen/interface.jl +++ b/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 diff --git a/src/MathProg/projection.jl b/src/MathProg/projection.jl index c546cbf35..69e54acf9 100644 --- a/src/MathProg/projection.jl +++ b/src/MathProg/projection.jl @@ -37,9 +37,11 @@ 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) @@ -74,6 +76,9 @@ function _mapping_by_subproblem(columns::Dict{Int, Vector{A}}, values::Dict{Int, ) 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}}() @@ -134,7 +139,8 @@ 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) - return isinteger(projected_sol) && all(isinteger, map(r -> values(r), rolls)) + integer_rolls = _subproblem_rolls_are_integer(rolls) + return isinteger(projected_sol) && integer_rolls end ############################################################################################ diff --git a/test/unit/ColGen/colgen_default.jl b/test/unit/ColGen/colgen_default.jl index f7b52e89a..e8c41da09 100644 --- a/test/unit/ColGen/colgen_default.jl +++ b/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) diff --git a/test/unit/ColGen/colgen_iteration.jl b/test/unit/ColGen/colgen_iteration.jl index 16582943f..dea22de8f 100644 --- a/test/unit/ColGen/colgen_iteration.jl +++ b/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) diff --git a/test/unit/MathProg/projection.jl b/test/unit/MathProg/projection.jl index cb22fbeed..209f446b6 100644 --- a/test/unit/MathProg/projection.jl +++ b/test/unit/MathProg/projection.jl @@ -16,8 +16,10 @@ function test_mapping_operator_1() [0, 0, 1, 1, 0, 0, 1], [0, 0, 1, 1, 0, 0, 1], ] + @test Coluna.MathProg._rolls_are_integer(result) + return end -register!(unit_tests, "projection", test_mapping_operator_1; f= true) +register!(unit_tests, "projection", test_mapping_operator_1) function test_mapping_operator_2() # Example from the paper: @@ -58,8 +60,10 @@ function test_mapping_operator_2() [0.0, 1.0, 1.0, 1.0], [0.0, 0.5, 0.5, 0.0] ] + @test Coluna.MathProg._rolls_are_integer(result) == false + return end -register!(unit_tests, "projection", test_mapping_operator_2; f=true) +register!(unit_tests, "projection", test_mapping_operator_2) # function test_mapping_operator_3() # G = Vector{Float64}[ @@ -198,4 +202,4 @@ function projection_from_dw_reform_to_master_1() @test Coluna.MathProg.proj_cols_is_integer(solution) == false end -register!(unit_tests, "projection", projection_from_dw_reform_to_master_1; f = true) +register!(unit_tests, "projection", projection_from_dw_reform_to_master_1) From 69dc653d4a32c1c92acdae7e6fcd6fc83daf9a83 Mon Sep 17 00:00:00 2001 From: Guillaume Marques Date: Tue, 11 Apr 2023 16:46:18 +0200 Subject: [PATCH 5/5] ok --- src/Algorithm/basic/cutcallback.jl | 2 +- src/Algorithm/branching/interface.jl | 2 +- src/Algorithm/colgen.jl | 2 +- src/Algorithm/colgen/printer.jl | 1 + src/optimize.jl | 4 ++-- 5 files changed, 6 insertions(+), 5 deletions(-) diff --git a/src/Algorithm/basic/cutcallback.jl b/src/Algorithm/basic/cutcallback.jl index dc41d896b..89c25f819 100644 --- a/src/Algorithm/basic/cutcallback.jl +++ b/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[] diff --git a/src/Algorithm/branching/interface.jl b/src/Algorithm/branching/interface.jl index 004895f84..2c292d7b8 100644 --- a/src/Algorithm/branching/interface.jl +++ b/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 diff --git a/src/Algorithm/colgen.jl b/src/Algorithm/colgen.jl index 0b437c606..c9a6b5de8 100644 --- a/src/Algorithm/colgen.jl +++ b/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) diff --git a/src/Algorithm/colgen/printer.jl b/src/Algorithm/colgen/printer.jl index 2bcbe6184..ee45f9971 100644 --- a/src/Algorithm/colgen/printer.jl +++ b/src/Algorithm/colgen/printer.jl @@ -36,6 +36,7 @@ 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) diff --git a/src/optimize.jl b/src/optimize.jl index ff5bde5df..fce18fb52 100644 --- a/src/optimize.jl +++ b/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