Subproblem local bound propagation

src/Algorithm/presolve/helpers.jl

@@ -3,7 +3,7 @@ const PRECISION_DIGITS = 6 # floating point numbers have between 6 and 9 signifi
 """
 Temporary data structure where we store a representation of the formulation that we presolve.
 """
-struct PresolveFormRepr
+mutable struct PresolveFormRepr
     nb_vars::Int
     nb_constrs::Int
     col_major_coef_matrix::SparseMatrixCSC{Float64,Int64} # col major

src/Algorithm/presolve/interface.jl

@@ -275,9 +275,8 @@ struct PresolveAlgorithm <: AlgoAPI.AbstractAlgorithm
 end
 
 # PresolveAlgorithm does not have child algorithms, therefore get_child_algorithms() is not defined
-
 function get_units_usage(algo::PresolveAlgorithm, reform::Reformulation) 
-    units_usage = Tuple{AbstractModel, UnitType, UnitPermission}[] 
+    units_usage = Tuple{AbstractModel, UnitType, UnitPermission}[]
     master = getmaster(reform)
     push!(units_usage, (master, StaticVarConstrUnit, READ_AND_WRITE))
     push!(units_usage, (master, MasterBranchConstrsUnit, READ_AND_WRITE))

src/Algorithm/presolve/propagation.jl

@@ -47,3 +47,68 @@ function propagate_var_bounds_from!(dest::PresolveFormulation, src::PresolveForm
     end
     return
 end
+
+function propagate_local_and_global_bounds!(reform::Reformulation, presolve_form_repr::DwPresolveReform)
+    master = getmaster(reform)
+    presolve_master = presolve_form_repr.restricted_master
+    for (spid, spform) in get_dw_pricing_sps(reform)
+        presolve_sp = presolve_form_repr.dw_sps[spid]
+        propagate_local_bounds!(presolve_sp, presolve_master, spform, master)
+    end
+    return
+end
+
+function propagate_global_bounds!(presolve_repr_master::PresolveFormulation, presolve_restricted_master::PresolveFormulation, master::Formulation)
+
+end
+
+function propagate_local_bounds!(presolve_sp::PresolveFormulation, presolve_master::PresolveFormulation, spform::Formulation, master::Formulation)
+    partial_solution = zeros(Float64, presolve_sp.form.nb_vars)
+
+    # Get new columns in partial solution.
+    pool = get_primal_sol_pool(spform)
+
+    # Get new columns in partial solution.
+    for (col, partial_sol_value) in enumerate(presolve_master.form.partial_solution)
+        if abs(partial_sol_value) > Coluna.TOL
+            var = presolve_master.col_to_var[col]
+            varid = getid(var)
+            column = @view pool.solutions[varid,:]
+            for (varid, val) in column
+                sp_var_col = presolve_sp.var_to_col[varid]
+                partial_solution[sp_var_col] += val * partial_sol_value
+            end
+        end
+    end
+
+    new_lbs = zeros(Float64, presolve_sp.form.nb_vars)
+    new_ubs = zeros(Float64, presolve_sp.form.nb_vars)
+
+    # when there is a partial solution, we update the bound so that the lower bound
+    # on the absolute value of the variable is improved.
+    # Examples:
+    # -3 <= x <= 3  &  x = 2    ->    2 <= x <= 3
+    # -3 <= x <= 3  &  x = -1   ->   -3 <= x <= -1
+    #  2 <= x <= 4  &  x = 3    ->    3 <= x <= 4
+    # -5 <= x <= 0  &  x = -2   ->   -5 <= x <= -2
+    for (col, (val, lb, ub)) in enumerate(Iterators.zip(partial_solution, presolve_sp.form.lbs, presolve_sp.form.ubs))
+        if val < 0
+            new_lbs[col] = lb - val
+            new_ubs[col] = min(0, ub - val)
+        elseif val > 0
+            new_lbs[col] = max(0, lb - val)
+            new_ubs[col] = ub - val
+        else
+            new_lbs[col] = lb
+            new_ubs[col] = ub
+        end
+
+        if lb > ub
+            error("Infeasible.")
+        end
+    end
+
+    presolve_sp.form.lbs = new_lbs
+    presolve_sp.form.ubs = new_ubs
+    return
+end

test/unit/Presolve/propagation.jl

@@ -336,6 +336,273 @@ function test_var_bound_propagation_within_restricted_master()
 end
 register!(unit_tests, "presolve_propagation", test_var_bound_propagation_within_restricted_master)
 
+function test_col_bounds_propagation_from_restricted_master()
+    # Original Master
+    # min x1 + x2 
+    # s.t. x1 + x2  >= 10
+    #      1 <= x1 <= 2 (repr)
+    #      1 <= x2 <= 3 (repr)
+
+    # Restricted master
+    # min 2MC1 + 3MC2 + 3MC3 + 1000a1
+    # s.t.  2MC1 + 3MC2 + 3MC3 + a1 >= 10
+    #       MC1 + MC2 + MC3 >= 0 (convexity)
+    #       MC1 + MC2 + MC3 <= 2 (convexity)
+    #       MC1, MC2, MC3 >= 0
+
+    # subproblem variables:
+    #   1 <= x1 <= 2, 
+    #   1 <= x2 <= 3
+
+    env = Coluna.Env{Coluna.MathProg.VarId}(Coluna.Params())
+
+    master_form, master_name_to_var, master_name_to_constr = _mathprog_formulation!(
+        env,
+        Coluna.MathProg.DwMaster(),
+        [
+            # name, duty, cost, lb, ub, id
+            ("x1", Coluna.MathProg.MasterRepPricingVar, 1.0, 1.0, 2.0, nothing),
+            ("x2", Coluna.MathProg.MasterRepPricingVar, 1.0, 1.0, 3.0, nothing),
+            ("MC1", Coluna.MathProg.MasterCol, 2.0, 0.0, Inf, nothing),
+            ("MC2", Coluna.MathProg.MasterCol, 3.0, 0.0, Inf, nothing),
+            ("MC3", Coluna.MathProg.MasterCol, 3.0, 0.0, Inf, nothing),
+            ("a", Coluna.MathProg.MasterArtVar, 1000.0, 0.0, Inf, nothing),
+            ("pricing_setup", Coluna.MathProg.MasterRepPricingSetupVar, 0.0, 0.0, 1.0, nothing)
+        ],
+        [
+            # name, duty, rhs, sense, id
+            ("c1", Coluna.MathProg.MasterMixedConstr, 10.0, ClMP.Greater, nothing),
+            ("c2", Coluna.MathProg.MasterConvexityConstr, 0.0, ClMP.Greater, nothing),
+            ("c3", Coluna.MathProg.MasterConvexityConstr, 2.0, ClMP.Less, nothing)
+        ]
+    )
+
+    spform, sp_name_to_var, sp_name_to_constr = _mathprog_formulation!(
+        env,
+        Coluna.MathProg.DwSp(
+            Coluna.MathProg.getid(master_name_to_var["pricing_setup"]),
+            Coluna.MathProg.getid(master_name_to_constr["c2"]),
+            Coluna.MathProg.getid(master_name_to_constr["c3"]),
+            Coluna.MathProg.Integ
+        ),
+        [
+            ("x1", Coluna.MathProg.DwSpPricingVar, 1.0, 1.0, 2.0, Coluna.Algorithm.getid(master_name_to_var["x1"])),
+            ("x2", Coluna.MathProg.DwSpPricingVar, 1.0, 1.0, 3.0, Coluna.Algorithm.getid(master_name_to_var["x2"]))
+        ],
+        []
+    )
+
+    var_ids = [Coluna.MathProg.getid(sp_name_to_var["x1"]), Coluna.MathProg.getid(sp_name_to_var["x2"])]
+    pool = Coluna.MathProg.get_primal_sol_pool(spform)
+
+    for (name, vals) in Iterators.zip(
+        ["MC1", "MC2", "MC3"],
+        [
+            #x1, x2,
+            Float64[1.0, 1.0],
+            Float64[1.0, 2.0],
+            Float64[2.0, 2.0]
+        ]
+    )
+        col_id = Coluna.MathProg.VarId(Coluna.MathProg.getid(master_name_to_var[name]); duty = Coluna.MathProg.DwSpPrimalSol)
+        Coluna.MathProg.push_in_pool!(
+            pool,
+            Coluna.MathProg.PrimalSolution(spform, var_ids, vals, 1.0, Coluna.MathProg.FEASIBLE_SOL),
+            col_id,
+            1.0
+        )
+    end
+
+    ## MC3 >= 1
+    Coluna.MathProg.setcurlb!(master_form, master_name_to_var["MC3"], 1.0)
+
+    ## We create the presolve formulations
+    master_presolve_form = _presolve_formulation(
+        ["x1", "x2"], ["c1"], [1 1], master_form, master_name_to_var, master_name_to_constr
+    )
+
+    restricted_presolve_form = _presolve_formulation(
+        ["MC1", "MC2", "MC3", "a", "pricing_setup"], ["c1", "c2", "c3"],
+        [2 3 3 1 0; 1 1 1 0 0; 1 1 1 0 1], master_form, master_name_to_var, master_name_to_constr
+    )
+
+    sp_presolve_form = _presolve_formulation(
+        ["x1", "x2"], [], [1 1], spform, sp_name_to_var, sp_name_to_constr
+    )
+
+    ## We run the presolve on the restricted master
+    tightened_bounds = Coluna.Algorithm.bounds_tightening(restricted_presolve_form.form)
+    new_restricted_master = Coluna.Algorithm.propagate_in_presolve_form(
+        restricted_presolve_form, Int[], tightened_bounds
+    )
+
+    Coluna.Algorithm.update_form_from_presolve!(master_form, new_restricted_master)
+    Coluna.Algorithm.propagate_local_bounds!(sp_presolve_form, new_restricted_master, spform, master_form)
+    Coluna.Algorithm.propagate_global_bounds!(master_presolve_form, new_restricted_master, master_form)
+    Coluna.Algorithm.update_form_from_presolve!(spform, sp_presolve_form)
+
+    # Fix MC3 [x1 = 2, x2 = 2]
+    # we should have 
+    #   0 <= x1 <= 0,
+    #   0 <= x2 <= 1
+
+    # Partial solution: MC3 = 1
+    @test Coluna.MathProg.get_value_in_partial_sol(master_form, master_name_to_var["MC1"]) ≈ 0
+    @test Coluna.MathProg.get_value_in_partial_sol(master_form, master_name_to_var["MC2"]) ≈ 0
+    @test Coluna.MathProg.get_value_in_partial_sol(master_form, master_name_to_var["MC3"]) ≈ 1
+
+    @test Coluna.MathProg.getcurlb(master_form, master_name_to_var["MC1"]) ≈ 0
+    @test Coluna.MathProg.getcurub(master_form, master_name_to_var["MC1"]) ≈ Inf
+    @test Coluna.MathProg.getcurlb(master_form, master_name_to_var["MC2"]) ≈ 0
+    @test Coluna.MathProg.getcurub(master_form, master_name_to_var["MC2"]) ≈ Inf
+    @test Coluna.MathProg.getcurlb(master_form, master_name_to_var["MC3"]) ≈ 0
+    @test Coluna.MathProg.getcurub(master_form, master_name_to_var["MC3"]) ≈ Inf
+
+    # Partial solution: x1 = 2, x2 = 2 |||| x1 = 0, 0 <= x2 <= 1
+    #@test Coluna.MathProg.get_value_in_partial_sol(spform, sp_name_to_var["x1"]) ≈ 2
+    #@test Coluna.MathProg.get_value_in_partial_sol(spform, sp_name_to_var["x2"]) ≈ 2
+
+    @test Coluna.MathProg.getcurlb(spform, sp_name_to_var["x1"]) ≈ 0
+    @test Coluna.MathProg.getcurub(spform, sp_name_to_var["x1"]) ≈ 0
+    @test Coluna.MathProg.getcurlb(spform, sp_name_to_var["x2"]) ≈ 0
+    @test Coluna.MathProg.getcurub(spform, sp_name_to_var["x2"]) ≈ 1
+end
+register!(unit_tests, "presolve_propagation", test_col_bounds_propagation_from_restricted_master)
+
+function test_col_bounds_propagation_from_restricted_master2()
+    # Original Master
+    # min x1 + x2 
+    # s.t. x1 + x2  >= 10
+    #      1 <= x1 <= 2 (repr)
+    #      1 <= x2 <= 3 (repr)
+
+    # Restricted master
+    # min 2MC1 + 3MC2 + 3MC3 + 1000a1
+    # s.t.  2MC1 + 3MC2 + 3MC3 + a1 >= 10
+    #       MC1 + MC2 + MC3 >= 0 (convexity)
+    #       MC1 + MC2 + MC3 <= 2 (convexity)
+    #       MC1, MC2, MC3 >= 0
+
+    # subproblem variables:
+    #   -3 <= x1 <=  3, 
+    #  -10 <= x2 <= -1
+
+    env = Coluna.Env{Coluna.MathProg.VarId}(Coluna.Params())
+
+    master_form, master_name_to_var, master_name_to_constr = _mathprog_formulation!(
+        env,
+        Coluna.MathProg.DwMaster(),
+        [
+            # name, duty, cost, lb, ub, id
+            ("x1", Coluna.MathProg.MasterRepPricingVar, 1.0, -3.0, 3.0, nothing),
+            ("x2", Coluna.MathProg.MasterRepPricingVar, 1.0, -10.0, -1.0, nothing),
+            ("MC1", Coluna.MathProg.MasterCol, 2.0, 0.0, Inf, nothing),
+            ("MC2", Coluna.MathProg.MasterCol, 3.0, 0.0, Inf, nothing),
+            ("MC3", Coluna.MathProg.MasterCol, 3.0, 0.0, Inf, nothing),
+            ("a", Coluna.MathProg.MasterArtVar, 1000.0, 0.0, Inf, nothing),
+            ("pricing_setup", Coluna.MathProg.MasterRepPricingSetupVar, 0.0, 0.0, 1.0, nothing)
+        ],
+        [
+            # name, duty, rhs, sense, id
+            ("c1", Coluna.MathProg.MasterMixedConstr, 10.0, ClMP.Greater, nothing),
+            ("c2", Coluna.MathProg.MasterConvexityConstr, 0.0, ClMP.Greater, nothing),
+            ("c3", Coluna.MathProg.MasterConvexityConstr, 2.0, ClMP.Less, nothing)
+        ]
+    )
+
+    spform, sp_name_to_var, sp_name_to_constr = _mathprog_formulation!(
+        env,
+        Coluna.MathProg.DwSp(
+            Coluna.MathProg.getid(master_name_to_var["pricing_setup"]),
+            Coluna.MathProg.getid(master_name_to_constr["c2"]),
+            Coluna.MathProg.getid(master_name_to_constr["c3"]),
+            Coluna.MathProg.Integ
+        ),
+        [
+            ("x1", Coluna.MathProg.DwSpPricingVar, 1.0, -3.0, 3.0, Coluna.Algorithm.getid(master_name_to_var["x1"])),
+            ("x2", Coluna.MathProg.DwSpPricingVar, 1.0, -10.0, -1.0, Coluna.Algorithm.getid(master_name_to_var["x2"]))
+        ],
+        []
+    )
+
+    var_ids = [Coluna.MathProg.getid(sp_name_to_var["x1"]), Coluna.MathProg.getid(sp_name_to_var["x2"])]
+    pool = Coluna.MathProg.get_primal_sol_pool(spform)
+
+    for (name, vals) in Iterators.zip(
+        ["MC1", "MC2", "MC3"],
+        [
+            #x1, x2,
+            Float64[-1.0, -1.0],
+            Float64[-1.0, -2.0],
+            Float64[-2.0, -2.0]
+        ]
+    )
+        col_id = Coluna.MathProg.VarId(Coluna.MathProg.getid(master_name_to_var[name]); duty = Coluna.MathProg.DwSpPrimalSol)
+        Coluna.MathProg.push_in_pool!(
+            pool,
+            Coluna.MathProg.PrimalSolution(spform, var_ids, vals, 1.0, Coluna.MathProg.FEASIBLE_SOL),
+            col_id,
+            1.0
+        )
+    end
+
+    ## MC3 >= 1
+    Coluna.MathProg.setcurlb!(master_form, master_name_to_var["MC2"], 1.0)
+
+    ## We create the presolve formulations
+    master_presolve_form = _presolve_formulation(
+        ["x1", "x2"], ["c1"], [1 1], master_form, master_name_to_var, master_name_to_constr
+    )
+
+    restricted_presolve_form = _presolve_formulation(
+        ["MC1", "MC2", "MC3", "a", "pricing_setup"], ["c1", "c2", "c3"],
+        [2 3 3 1 0; 1 1 1 0 0; 1 1 1 0 1], master_form, master_name_to_var, master_name_to_constr
+    )
+
+    sp_presolve_form = _presolve_formulation(
+        ["x1", "x2"], [], [1 1], spform, sp_name_to_var, sp_name_to_constr
+    )
+
+    ## We run the presolve on the restricted master
+    tightened_bounds = Coluna.Algorithm.bounds_tightening(restricted_presolve_form.form)
+    new_restricted_master = Coluna.Algorithm.propagate_in_presolve_form(
+        restricted_presolve_form, Int[], tightened_bounds
+    )
+
+    Coluna.Algorithm.update_form_from_presolve!(master_form, new_restricted_master)
+    Coluna.Algorithm.propagate_local_bounds!(sp_presolve_form, new_restricted_master, spform, master_form)
+    Coluna.Algorithm.propagate_global_bounds!(master_presolve_form, new_restricted_master, master_form)
+    Coluna.Algorithm.update_form_from_presolve!(spform, sp_presolve_form)
+
+    # Fix MC2 [x1 = -1, x2 = -2]
+    # we should have 
+    #   -2 <= x1 <= 0,
+    #   -8 <= x2 <= 0
+
+
+    # Partial solution: MC3 = 1
+    @test Coluna.MathProg.get_value_in_partial_sol(master_form, master_name_to_var["MC1"]) ≈ 0
+    @test Coluna.MathProg.get_value_in_partial_sol(master_form, master_name_to_var["MC2"]) ≈ 1
+    @test Coluna.MathProg.get_value_in_partial_sol(master_form, master_name_to_var["MC3"]) ≈ 0
+
+    @test Coluna.MathProg.getcurlb(master_form, master_name_to_var["MC1"]) ≈ 0
+    @test Coluna.MathProg.getcurub(master_form, master_name_to_var["MC1"]) ≈ Inf
+    @test Coluna.MathProg.getcurlb(master_form, master_name_to_var["MC2"]) ≈ 0
+    @test Coluna.MathProg.getcurub(master_form, master_name_to_var["MC2"]) ≈ Inf
+    @test Coluna.MathProg.getcurlb(master_form, master_name_to_var["MC3"]) ≈ 0
+    @test Coluna.MathProg.getcurub(master_form, master_name_to_var["MC3"]) ≈ Inf
+
+    # Partial solution: x1 = 2, x2 = 2 |||| x1 = 0, 0 <= x2 <= 1
+    # @test Coluna.MathProg.get_value_in_partial_sol(spform, sp_name_to_var["x1"]) ≈ 2
+    # @test Coluna.MathProg.get_value_in_partial_sol(spform, sp_name_to_var["x2"]) ≈ 2
+
+    @test Coluna.MathProg.getcurlb(spform, sp_name_to_var["x1"]) ≈ -2
+    @test Coluna.MathProg.getcurub(spform, sp_name_to_var["x1"]) ≈ 0
+    @test Coluna.MathProg.getcurlb(spform, sp_name_to_var["x2"]) ≈ -8
+    @test Coluna.MathProg.getcurub(spform, sp_name_to_var["x2"]) ≈ 0
+end
+register!(unit_tests, "presolve_propagation", test_col_bounds_propagation_from_restricted_master2)
+
 ## OriginalVar -> DwSpPricingVar (mapping exists)
 ## otherwise no propagation
 function test_var_bound_propagation_from_original_to_subproblem()