test angle; fix subgradient calculation (#935)

src/Algorithm/colgen/stabilization.jl

@@ -141,33 +141,43 @@ function _primal_solution(master::Formulation, generated_columns, is_minimizatio
     return sparsevec(var_ids, var_vals)
 end
 
-function _dynamic_alpha_schedule(
-    stab::ColGenStab, smooth_dual_sol, h, primal_solution, is_minimization
-)
+function _increase(smooth_dual_sol, cur_stab_center, h, primal_solution, is_minimization)
     # Calculate the in-sep direction.
-    in_sep_direction = smooth_dual_sol - stab.cur_stab_center
+    in_sep_direction = smooth_dual_sol - cur_stab_center
     in_sep_dir_norm = norm(in_sep_direction)
 
+    # if in & sep are the same point, we need to decrease α becase it is the weight of the
+    # stability center (in) in the formula to compute the sep point.
+    if iszero(in_sep_dir_norm)
+        return false
+    end
+
     # Calculate the subgradient
     subgradient = h.a - h.A * primal_solution
     subgradient_norm = norm(subgradient)
 
     # we now calculate the angle between the in-sep direction and the subgradient 
-    angle = (transpose(in_sep_direction) * subgradient) / (in_sep_dir_norm * subgradient_norm)
+    cos_angle = (transpose(in_sep_direction) * subgradient) / (in_sep_dir_norm * subgradient_norm)
+
     if !is_minimization
-        angle *= -1
+        cos_angle *= -1
     end
+    return cos_angle < 1e-12
+end
 
+function _dynamic_alpha_schedule(
+    α, smooth_dual_sol, cur_stab_center, h, primal_solution, is_minimization
+)
+    increase = _increase(smooth_dual_sol, cur_stab_center, h, primal_solution, is_minimization)
     # we modify the alpha parameter based on the calculated angle
-    α = angle > 1e-12 ? f_decr(stab.base_α) : f_incr(stab.base_α)
-    return α
+    return increase ? f_incr(α) : f_decr(α)
 end
 
 function ColGen.update_stabilization_after_iter!(stab::ColGenStab, ctx, master, generated_columns, mast_dual_sol)
     if stab.smooth_factor == 1
         is_min = ColGen.is_minimization(ctx)
         primal_sol = _primal_solution(master, generated_columns, is_min)
-        α = _dynamic_alpha_schedule(stab, mast_dual_sol, subgradient_helper(ctx), primal_sol, is_min)
+        α = _dynamic_alpha_schedule(stab.base_α, mast_dual_sol, stab.cur_stab_center, subgradient_helper(ctx), primal_sol, is_min)
         stab.base_α = α
     end
 

src/Algorithm/colgen/utils.jl

@@ -47,17 +47,20 @@ function _submatrix(
     form::Formulation, 
     keep_constr::Function, 
     keep_var::Function,
+    m::Function = (form, is_min, constr_id, var_id) -> 1.0
 )
+    is_min = getobjsense(form) == MinSense
     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)
+                c = m(form, is_min, constr_id, var_id)
                 push!(constr_ids, constr_id)
                 push!(var_ids, var_id)
-                push!(nz, coeff)
+                push!(nz, c * coeff)
             end
         end
     end
@@ -154,18 +157,34 @@ end
 function SubgradientCalculationHelper(master)
     constr_ids = ConstrId[]
     constr_rhs = Float64[]
+
+    m_rhs = (master, is_min, constr_id) -> begin
+        constr_sense = getcursense(master, constr_id)
+        if is_min
+            return constr_sense == Less ? -1.0 : 1.0
+        else
+            return constr_sense == Greater ? -1.0 : 1.0
+        end
+    end
+    m_submatrix = (master, is_min, constr_id, var_id) -> begin
+        m_rhs(master, is_min, constr_id)
+    end
+
+    is_min = getobjsense(master) == MinSense
     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))
+            push!(constr_rhs, m_rhs(master, is_min, constr_id) * 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
+        var_id -> getduty(var_id) <= MasterPureVar || getduty(var_id) <= MasterRepPricingVar,
+        m_submatrix
     )
     return SubgradientCalculationHelper(a, A)
 end

test/unit/ColGen/colgen_default.jl

@@ -155,9 +155,9 @@ function test_subgradient_calculation_helper()
 
     helper = ClA.SubgradientCalculationHelper(master)
     @test helper.a[cids["c1"]] == 10
-    @test helper.a[cids["c2"]] == 100
+    @test helper.a[cids["c2"]] == -100
     @test helper.a[cids["c3"]] == 100
-    @test helper.a[cids["c4"]] == 5
+    @test helper.a[cids["c4"]] == -5
 
     @test helper.A[cids["c1"], vids["x1"]] == 1
     @test helper.A[cids["c1"], vids["x2"]] == 1
@@ -167,22 +167,22 @@ function test_subgradient_calculation_helper()
     @test helper.A[cids["c1"], vids["y3"]] == 1
     @test helper.A[cids["c1"], vids["z1"]] == 2
     @test helper.A[cids["c1"], vids["z2"]] == 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["c2"], vids["z1"]] == 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["c2"], vids["z1"]] == -1
     @test helper.A[cids["c2"], vids["z2"]] == 0
     @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.A[cids["c3"], vids["z1"]] == 0
     @test helper.A[cids["c3"], vids["z2"]] == 0
-    @test helper.A[cids["c4"], vids["z1"]] == 1
-    @test helper.A[cids["c4"], vids["z2"]] == 1
+    @test helper.A[cids["c4"], vids["z1"]] == -1
+    @test helper.A[cids["c4"], vids["z2"]] == -1
 end
-register!(unit_tests, "colgen_default", test_subgradient_calculation_helper)
+register!(unit_tests, "colgen_default", test_subgradient_calculation_helper; f = true)
 
 # All the tests are based on the Generalized Assignment problem.
 # x_mj = 1 if job j is assigned to machine m

test/unit/ColGen/colgen_stabilization.jl

@@ -64,7 +64,6 @@ form_primal_solution() = """
             z2 >= 3
 """
 
-
 function test_primal_solution()
     _, master, sps, _, _ = reformfromstring(form_primal_solution())
 
@@ -113,11 +112,180 @@ function test_primal_solution()
 end
 register!(unit_tests, "colgen_stabilization", test_primal_solution; f = true)
 
+form_primal_solution2() = """
+    master
+        min
+        3x_11 + 2x_12 + 2x_21 + x_22 + z1 + z2 + 0s1 + 0s2
+        s.t.
+        x_11  + x_12  +  x_21 +         2z1 + z2  >= 10
+        x_11  + 2x_12 +  x_21 + 2x_22 +  z1       <= 100
+        x_11  +                  x_22 +           == 100
+                                         z1  + z2 >= 5  
+        s1                                        >= 1 {MasterConvexityConstr}
+        s1                                        <= 2 {MasterConvexityConstr}
+        s2                                        >= 0 {MasterConvexityConstr}
+        s2                                        <= 3 {MasterConvexityConstr}
+
+    dw_sp
+        min
+        x_11 + x_12 + 0s1
+        s.t.
+        x_11 + x_12  >= 10
+
+    dw_sp
+        min
+        x_21 + x_22 + 0s2
+        s.t.
+        x_21 + x_22  >= 10
+        
+        integer
+            representatives
+                x_11, x_12, x_21, x_22
+
+            pure
+                z1, z2
+
+            pricing_setup
+                s1, s2
+        
+        bounds
+            x_11 >= 0
+            x_12 >= 0
+            x_21 >= 0
+            x_22 >= 1
+            z1 >= 0
+            z2 >= 3
+"""
+
+# We consider the master with the following coefficient matrix:
+# master_coeff_matrix = [
+#      1 1 1 0 1 1;
+#      -1 -2 -1 -2 -1 0;
+#      1 0 0 1 0 0;      # is it correct to handle an "==" constraint like this in subgradient computation?
+#      0 0 0 0 1 1;
+# ]
+# the following rhs: rhs = [10, -100, 100, 5]
+
+# We consider the primal solution: primal = [1, 2, 12, 12, 0, 0]
+# The subgradient is therefore: rhs - master_coeff_matrix * primal = [-5, -59, 87, 5]
+# We use the following stability center:  stab = [1, 2, 0, 1]
+
+function _test_angle_primal_sol(master, sp1, sp2)
+    vids = get_name_to_varids(master)
+    pool = Coluna.Algorithm.ColumnsSet()
+
+    sol1 = Coluna.MathProg.PrimalSolution(
+        sp1,
+        [vids["x_11"], vids["x_12"]],
+        [1.0, 2.0],
+        11.0,
+        Coluna.MathProg.FEASIBLE_SOL
+    )
+    Coluna.Algorithm.add_primal_sol!(pool.subprob_primal_sols, sol1, false) # improving = true
+
+    sol2 = Coluna.MathProg.PrimalSolution(
+        sp2,
+        [vids["x_21"], vids["x_22"]],
+        [4.0, 4.0],
+        13.0,
+        Coluna.MathProg.FEASIBLE_SOL
+    )
+    Coluna.Algorithm.add_primal_sol!(pool.subprob_primal_sols, sol2, true)
+    is_minimization = true
+    primal_sol = Coluna.Algorithm._primal_solution(master, pool, is_minimization)
+    return primal_sol
+end
+
+function _test_angle_stab_center(master)
+    cids = get_name_to_constrids(master)
+    return Coluna.MathProg.DualSolution(
+        master,
+        [cids["c1"], cids["c2"], cids["c4"]],
+        [1.0, 2.0, 1.0],
+        Coluna.MathProg.VarId[], Float64[], Coluna.MathProg.ActiveBound[],
+        0.0,
+        Coluna.MathProg.FEASIBLE_SOL
+    )
+end
+
 # Make sure the angle is well computed.
-function test_angle()
+# Here we test the can where the in and sep points are the same.
+# In that case, we should decrease the value of α.
+function test_angle_1()
+    _, master, sps, _, _ = reformfromstring(form_primal_solution2())
+    sp1, sp2 = sps[2], sps[1]
+    cids = get_name_to_constrids(master)
 
+    smooth_dual_sol = Coluna.MathProg.DualSolution(
+        master,
+        [cids["c1"], cids["c2"], cids["c4"]],
+        [1.0, 2.0, 1.0],
+        Coluna.MathProg.VarId[], Float64[], Coluna.MathProg.ActiveBound[],
+        0.0,
+        Coluna.MathProg.FEASIBLE_SOL
+    )
+    cur_stab_center = _test_angle_stab_center(master)
+
+    h = Coluna.Algorithm.SubgradientCalculationHelper(master)
+    is_minimization = true
+    primal_sol = _test_angle_primal_sol(master, sp1, sp2)
+    increase = Coluna.Algorithm._increase(smooth_dual_sol, cur_stab_center, h, primal_sol, is_minimization)
+    @test increase == false
+end
+register!(unit_tests, "colgen_stabilization", test_angle_1; f = true)
+
+# Let's consider the following sep point: sep = [5, 7, 0, 3]
+# The direction will be [4, 5, 0, 2] and should lead to a negative cosinus for the angle.
+# In that case, we need to increase the value of α.
+function test_angle_2()
+    _, master, sps, _, _ = reformfromstring(form_primal_solution2())
+    sp1, sp2 = sps[2], sps[1]
+    cids = get_name_to_constrids(master)
+
+    smooth_dual_sol = Coluna.MathProg.DualSolution(
+        master,
+        [cids["c1"], cids["c2"], cids["c4"]],
+        [5.0, 7.0, 3.0],
+        Coluna.MathProg.VarId[], Float64[], Coluna.MathProg.ActiveBound[],
+        0.0,
+        Coluna.MathProg.FEASIBLE_SOL
+    )
+    cur_stab_center = _test_angle_stab_center(master)
+
+    h = Coluna.Algorithm.SubgradientCalculationHelper(master)
+    is_minimization = true
+    primal_sol = _test_angle_primal_sol(master, sp1, sp2)
+    increase = Coluna.Algorithm._increase(smooth_dual_sol, cur_stab_center, h, primal_sol, is_minimization)
+    @test increase == true
 end
-register!(unit_tests, "colgen_stabilization", test_angle; f = true)
+register!(unit_tests, "colgen_stabilization", test_angle_2; f = true)
+
+# Let's consider the following sep point: sep = [5, 1, 10, 3]
+# The direction will be [4, 1, 10, 2] and should lead to a positive cosinus for the angle.
+# In that case, we need to decrease the value of α.
+function test_angle_3()
+    _, master, sps, _, _ = reformfromstring(form_primal_solution2())
+    sp1, sp2 = sps[2], sps[1]
+    cids = get_name_to_constrids(master)
+
+    smooth_dual_sol = Coluna.MathProg.DualSolution(
+        master,
+        [cids["c1"], cids["c2"], cids["c3"], cids["c4"]],
+        [5.0, 1.0, 10.0, 3.0],
+        Coluna.MathProg.VarId[], Float64[], Coluna.MathProg.ActiveBound[],
+        0.0,
+        Coluna.MathProg.FEASIBLE_SOL
+    )
+    cur_stab_center = _test_angle_stab_center(master)
+
+    h = Coluna.Algorithm.SubgradientCalculationHelper(master)
+    is_minimization = true
+    primal_sol = _test_angle_primal_sol(master, sp1, sp2)
+    increase = Coluna.Algorithm._increase(smooth_dual_sol, cur_stab_center, h, primal_sol, is_minimization)
+    @test increase == false
+end
+register!(unit_tests, "colgen_stabilization", test_angle_3; f = true)
+
 
 function test_dynamic_alpha_schedule()