fix isinteger + continuous columns if sp has only continuous vars (#…

…551)

* fix #550

* correct

src/Algorithm/colgen.jl

@@ -235,19 +235,15 @@ set_bestcol_id!(spinfo::SubprobInfo, varid::VarId) = spinfo.bestcol_id = varid
 function insert_cols_in_master!(
     masterform::Formulation, spinfo::SubprobInfo, phase::Int64, spform::Formulation,
 )
-    sp_uid = getuid(spform)
     nb_of_gen_col = 0
 
     for sol_id in spinfo.recorded_sol_ids
         nb_of_gen_col += 1
         name = string("MC_", getsortuid(sol_id))
         lb = 0.0
         ub = Inf
-        kind = Continuous
         duty = MasterCol
-        mc = setcol_from_sp_primalsol!(
-            masterform, spform, sol_id, name, duty; lb=lb, ub=ub, kind=kind
-        )
+        mc = setcol_from_sp_primalsol!(masterform, spform, sol_id, name, duty; lb=lb, ub=ub)
         if phase == 1
             setcurcost!(masterform, mc, 0.0)
         end

src/MathProg/duties.jl

@@ -19,6 +19,7 @@ mutable struct DwSp <: AbstractSpDuty
     setup_var::Union{VarId, Nothing}
     lower_multiplicity::Int
     upper_multiplicity::Int
+    column_var_kind::VarKind
 end
 
 "A Benders subproblem of a formulation decomposed using Benders."
@@ -47,14 +48,14 @@ struct BendersSp <: AbstractSpDuty end
         AbstractDwSpVar <= Duty{Variable}
             DwSpPricingVar <= AbstractDwSpVar
             DwSpSetupVar <= AbstractDwSpVar
-            DwSpPureVar <= AbstractDwSpVar
+            #DwSpPureVar <= AbstractDwSpVar
             DwSpPrimalSol <= AbstractDwSpVar
         AbstractBendSpVar <= Duty{Variable}
             AbstractBendSpSlackMastVar <= AbstractBendSpVar
                 BendSpSlackFirstStageVar <= AbstractBendSpSlackMastVar
                 BendSpSlackSecondStageCostVar <= AbstractBendSpSlackMastVar
             BendSpSepVar <= AbstractBendSpVar
-            BendSpPureVar <= AbstractBendSpVar
+            #BendSpPureVar <= AbstractBendSpVar
             BendSpPrimalSol <= AbstractBendSpVar
 end
 
@@ -108,11 +109,11 @@ function isaStaticDuty(duty::NestedEnum)
     duty <= MasterRepPricingSetupVar ||
     duty <= DwSpPricingVar ||
     duty <= DwSpSetupVar ||
-    duty <= DwSpPureVar ||
+    #duty <= DwSpPureVar ||
     duty <= DwSpPrimalSol ||
     duty <= DwSpDualSol ||
     duty <= BendSpSepVar ||
-    duty <= BendSpPureVar ||
+    #duty <= BendSpPureVar ||
     duty <= BendSpSlackFirstStageVar  ||
     duty <= BendSpSlackSecondStageCostVar ||
     duty <= OriginalConstr ||

src/MathProg/formulation.jl

@@ -301,7 +301,7 @@ end
 
 function setcol_from_sp_primalsol!(
     masterform::Formulation, spform::Formulation, sol_id::VarId, name::String,
-    duty::Duty{Variable}; lb::Float64 = 0.0, ub::Float64 = Inf, kind::VarKind = Continuous,
+    duty::Duty{Variable}; lb::Float64 = 0.0, ub::Float64 = Inf,
     inc_val::Float64 = 0.0, is_active::Bool = true, is_explicit::Bool = true,
     moi_index::MoiVarIndex = MoiVarIndex(), custom_data::Union{Nothing, AbstractCustomData} = nothing
 )
@@ -322,7 +322,7 @@ function setcol_from_sp_primalsol!(
         cost = cost,
         lb = lb,
         ub = ub,
-        kind = kind,
+        kind = spform.duty_data.column_var_kind,
         inc_val = inc_val,
         is_active = is_active,
         is_explicit = is_explicit,
@@ -331,6 +331,9 @@ function setcol_from_sp_primalsol!(
         id = sol_id,
         custom_data = get(spform.manager.primal_sols_custom_data, sol_id, nothing)
     )
+
+    setcurkind!(masterform, mast_col, Continuous)
+
     return mast_col
 end
 
@@ -481,29 +484,36 @@ function _addlocalartvar!(form::Formulation, constr::Constraint, abs_cost::Float
     return
 end
 
+"""
+    enforce_integrality!(formulation)
+
+Set the current kind of each active & explicit variable of the formulation to its perennial kind.
+"""
 function enforce_integrality!(form::Formulation)
-    @logmsg LogLevel(-1) string("Enforcing integrality of formulation ", getuid(form))
-    for (varid, var) in getvars(form)
-        !iscuractive(form, varid) && continue
-        !isexplicit(form, varid) && continue
-        getcurkind(form, varid) == Integ && continue
-        getcurkind(form, varid) == Binary && continue
-        if getduty(varid) <= MasterCol || getperenkind(form, varid) != Continuous
-            @logmsg LogLevel(-3) string("Setting kind of var ", getname(form, var), " to Integer")
-            setcurkind!(form, varid, Integ)
+    for (_, var) in getvars(form)
+        enforce = iscuractive(form, var) && isexplicit(form, var)
+        enforce &= getcurkind(form, var) === Continuous
+        enforce &= getperenkind(form, var) !== Continuous
+        if enforce
+            setcurkind!(form, var, getperenkind(form, var))
         end
     end
     return
 end
 
-function relax_integrality!(form::Formulation) # TODO remove : should be in Algorithm
-    @logmsg LogLevel(-1) string("Relaxing integrality of formulation ", getuid(form))
-    for (varid, var) in getvars(form)
-        !iscuractive(form, varid) && continue
-        !isexplicit(form, varid) && continue
-        getcurkind(form, var) == Continuous && continue
-        @logmsg LogLevel(-3) string("Setting kind of var ", getname(form, var), " to continuous")
-        setcurkind!(form, varid, Continuous)
+"""
+    relax_integrality!(formulation)
+
+Set the current kind of each active & explicit integer or binary variable of the formulation
+to continuous.
+"""
+function relax_integrality!(form::Formulation)
+    for (_, var) in getvars(form)
+        relax = iscuractive(form, var) && isexplicit(form, var)
+        relax &= getcurkind(form, var) !== Continuous
+        if relax
+            setcurkind!(form, var, Continuous)
+        end
     end
     return
 end

src/MathProg/solutions.jl

@@ -20,15 +20,13 @@ end
 
 function Base.isinteger(sol::Solution)
     for (vc_id, val) in sol
-        #if getperenkind(sol.model, vc_id) != Continuous
-            abs(round(val) - val) <= 1e-5 || return false
-        #end
+        if getperenkind(sol.model, vc_id) !== Continuous && abs(round(val) - val) > 1e-5
+            return false
+        end
     end
     return true
 end
 
-isfractional(sol::Solution) = !Base.isinteger(sol)
-
 function contains(sol::PrimalSolution, f::Function)
     for (varid, val) in sol
         f(varid) && return true

src/decomposition.jl

@@ -60,7 +60,7 @@ function instantiatesp!(
 )
     spform = create_formulation!(
         env,
-        MathProg.DwSp(nothing, 1, 1);
+        MathProg.DwSp(nothing, 1, 1, Integ);
         parent_formulation = masterform,
         obj_sense = getobjsense(masterform)
     )
@@ -138,7 +138,9 @@ function create_side_vars_constrs!(
         @assert length(setupvars) == 1
         setupvar = collect(values(setupvars))[1]
         setuprepvar = clonevar!(origform, masterform, spform, setupvar, MasterRepPricingSetupVar, is_explicit = false)
-        # create convexity constraint
+
+        # create convexity constraint & storing information about the convexity constraint
+        # in the duty data of the formulation
         lb_mult = Float64(BD.getlowermultiplicity(ann))
         name = string("sp_lb_", spuid)
         lb_conv_constr = setconstr!(
@@ -156,11 +158,24 @@ function create_side_vars_constrs!(
             kind = Essential, sense = Less, inc_val = 100.0, 
             loc_art_var_abs_cost = env.params.local_art_var_cost
         )
+        masterform.parent_formulation.dw_pricing_sp_ub[spuid] = getid(ub_conv_constr)
+        coefmatrix[getid(ub_conv_constr), getid(setuprepvar)] = 1.0
+
         spform.duty_data.lower_multiplicity = lb_mult
         spform.duty_data.upper_multiplicity = ub_mult
         spform.duty_data.setup_var = getid(setupvar)
-        masterform.parent_formulation.dw_pricing_sp_ub[spuid] = getid(ub_conv_constr)
-        coefmatrix[getid(ub_conv_constr), getid(setuprepvar)] = 1.0
+
+        # If pricing subproblem variables are continuous, the master columns generated by
+        # the subproblem must have a continuous perenkind.
+        # This piece of information is stored in the duty data of the formulation.
+        continuous_columns = true
+        for (varid, var) in getvars(spform)
+            if getduty(varid) <= DwSpPricingVar && getperenkind(spform, var) !== Continuous
+                continuous_columns = false
+                break
+            end
+        end
+        spform.duty_data.column_var_kind = continuous_columns ? Continuous : Integ
     end
     return
 end

test/issues_tests.jl

@@ -234,6 +234,78 @@ function branching_file_completion()
     @test JuMP.termination_status(model) == MOI.OPTIMAL
 end
 
+# Issue https://github.com/atoptima/Coluna.jl/issues/550
+function continuous_vars_in_sp()
+    # Simple min cost flow problem
+    #   
+    # n1 ------ f[1] -------> n3
+    #  \                     ^
+    #   \                   /
+    #    -f[2]-> n2 --f[3]--
+    # 
+    #  n1: demand = -10.1
+    #  n2: demand = 0
+    #  n3: demand = 10.1 
+    #  f[1]: cost = 0, capacity = 8.5, in mip model only integer flow allowed
+    #  f[2]: cost = 50, (capacity = 5 can be activated by removing comment at constraint, line 93)
+    #  f[3]: cost = 50
+    #
+    #  Correct solution for non-integer f[1]
+    #           f[1] = 8.5, f[2] = f[3] = 1.6, cost = 8.5*0 + 1.6*2*50 = 160
+    #  Correct solution for integer f[1]
+    #           f[1] = 8, f[2] = f[3] = 2.1, cost = 8.5*0 + 2.1*2*50 = 210
+    #
+    function solve_flow_model(f1_integer, coluna)
+        @axis(M, 1:1)
+        model = BlockDecomposition.BlockModel(coluna, direct_model=true)
+        @variable(model, f[1:3, m in M])
+        if f1_integer
+            JuMP.set_integer(f[1, 1])
+        end
+        @constraint(model, n1[m in M], f[1,m] + f[2,m] == 10.1)
+        @constraint(model, n2[m in M], f[2,m] == f[3,m])
+        @constraint(model, n3[m in M], f[1,m] + f[3,m] == 10.1)
+        @constraint(model, cap1, sum(f[1,m] for m in M) <= 8.5)
+        #@JuMP.constraint(model, cap2, sum(f[2,m] for m in M) <= 5)
+        @objective(model, Min, 50 * f[2,1] + 50 * f[3,1])
+
+        @dantzig_wolfe_decomposition(model, decomposition, M)
+
+        subproblems = BlockDecomposition.getsubproblems(decomposition)
+        BlockDecomposition.specify!.(subproblems, lower_multiplicity=1, upper_multiplicity=1)
+
+        optimize!(model)
+
+        if f1_integer
+            @test termination_status(model) == MOI.OPTIMAL
+            @test primal_status(model) == MOI.FEASIBLE_POINT
+            @test objective_value(model) ≈ 210
+            @test value(f[1,1]) ≈ 8
+            @test value(f[2,1]) ≈ 2.1
+            @test value(f[3,1]) ≈ 2.1
+        else
+            @test termination_status(model) == MOI.OPTIMAL
+            @test primal_status(model) == MOI.FEASIBLE_POINT
+            @test objective_value(model) ≈ 160
+            @test value(f[1,1]) ≈ 8.5
+            @test value(f[2,1]) ≈ 1.6
+            @test value(f[3,1]) ≈ 1.6
+        end
+    end
+
+    coluna = JuMP.optimizer_with_attributes(
+        Coluna.Optimizer,
+        "params" => Coluna.Params(
+            solver=Coluna.Algorithm.TreeSearchAlgorithm(),
+        ),
+        "default_optimizer" => GLPK.Optimizer 
+    );
+
+    solve_flow_model(false, coluna)
+    solve_flow_model(true, coluna)
+    return
+end
+
 function test_issues_fixed()
     @testset "no_decomposition" begin
         solve_with_no_decomposition()
@@ -262,6 +334,10 @@ function test_issues_fixed()
     @testset "branching_file_completion" begin
         branching_file_completion()
     end
+
+    @testset "continuous_vars_in_sp" begin
+        continuous_vars_in_sp()
+    end
 end
 
 test_issues_fixed()