Merge branch 'master' into doc

atoptima · Aug 10, 2021 · b6a9ecb · b6a9ecb
2 parents a7956e0 + b4a66a6
commit b6a9ecb
Show file tree

Hide file tree

Showing 13 changed files with 398 additions and 104 deletions.
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
@@ -29,6 +29,7 @@ jobs:
       - run: julia --project=@. -e 'using Pkg; Pkg.add(PackageSpec(url="https://github.com/atoptima/ColunaDemos.jl.git"));'
       - run: julia --project=@. -e 'using Pkg; Pkg.add(PackageSpec(url="https://github.com/rafaelmartinelli/KnapsackLib.jl.git"));'
       - run: julia --project=@. -e 'using Pkg; Pkg.add(PackageSpec(url="https://github.com/atoptima/BlockDecomposition.jl.git"));'
+      - run: julia --project=@. -e 'using Pkg; Pkg.add(PackageSpec(url="https://github.com/atoptima/DynamicSparseArrays.jl.git"));'
       - uses: actions/cache@v1
         env:
           cache-name: cache-artifacts

diff --git a/src/Algorithm/benders.jl b/src/Algorithm/benders.jl
@@ -24,7 +24,6 @@ function get_units_usage(algo::BendersCutGeneration, reform::Reformulation)
     end
     return units_usage
 end
-
 mutable struct BendersCutGenRuntimeData
     optstate::OptimizationState
     spform_phase::Dict{FormId, FormulationPhase}
@@ -33,13 +32,6 @@ mutable struct BendersCutGenRuntimeData
     #slack_cost_increase_applied::Bool
 end
 
-function all_sp_in_phase2(algdata::BendersCutGenRuntimeData)
-    for (key, phase) in algdata.spform_phase
-        phase != PurePhase2 && return false
-    end
-    return true
-end
-
 function BendersCutGenRuntimeData(form::Reformulation, init_optstate::OptimizationState)
     optstate = OptimizationState(getmaster(form))
     best_ip_primal_sol = get_best_ip_primal_sol(init_optstate)
@@ -61,10 +53,12 @@ function run!(
 end
 
 function update_benders_sp_slackvar_cost_for_ph1!(spform::Formulation)
+    slack_cost = getobjsense(spform) === MinSense ? 1.0 : -1.0
     for (varid, var) in getvars(spform)
         iscuractive(spform, varid) || continue
-        if getduty(varid) == BendSpSlackFirstStageVar
-            setcurcost!(spform, var, 1.0)
+        if getduty(varid) <= BendSpSlackFirstStageVar
+            setcurcost!(spform, var, slack_cost)
+            setcurub!(spform, var, getperenub(spform, var))
         else
             setcurcost!(spform, var, 0.0)
         end
@@ -76,7 +70,7 @@ end
 function update_benders_sp_slackvar_cost_for_ph2!(spform::Formulation) 
     for (varid, var) in getvars(spform)
         iscuractive(spform, varid) || continue
-        if getduty(varid) == BendSpSlackFirstStageVar
+        if getduty(varid) <= BendSpSlackFirstStageVar
             setcurcost!(spform, var, 0.0)
             setcurub!(spform, var, 0.0)
         else
@@ -169,9 +163,10 @@ function record_solutions!(
 )::Vector{ConstrId}
 
     recorded_dual_solution_ids = Vector{ConstrId}()
+    sense = getobjsense(spform) === MinSense ? 1.0 : -1.0
 
     for dual_sol in get_lp_dual_sols(spresult)
-        if getvalue(dual_sol) > algo.feasibility_tol 
+        if sense * getvalue(dual_sol) > algo.feasibility_tol 
             (insertion_status, dual_sol_id) = setdualsol!(spform, dual_sol)
             if insertion_status
                 push!(recorded_dual_solution_ids, dual_sol_id)
@@ -187,9 +182,8 @@ end
 function insert_cuts_in_master!(
     masterform::Formulation, spform::Formulation, sp_dualsol_ids::Vector{ConstrId},
 )
-    sp_uid = getuid(spform)
     nb_of_gen_cuts = 0
-    sense = (getobjsense(masterform) == MinSense ? Greater : Less)
+    sense = getobjsense(masterform) == MinSense ? Greater : Less
 
     for dual_sol_id in sp_dualsol_ids
         nb_of_gen_cuts += 1
@@ -265,10 +259,14 @@ function solve_sp_to_gencut!(
         # Solve sub-problem and insert generated cuts in master
         # @logmsg LogLevel(-3) "optimizing benders_sp prob"
         TO.@timeit Coluna._to "Bender Sep SubProblem" begin
-            optstate = run!(SolveLpForm(get_dual_solution = true), env, spform, OptimizationInput(OptimizationState(spform)))
+            optstate = run!(
+                SolveLpForm(get_dual_solution = true, relax_integrality = true), 
+                env, spform, OptimizationInput(OptimizationState(spform))
+            )
         end
 
         optresult = getoptstate(optstate)
+
         if getterminationstatus(optresult) == INFEASIBLE # if status != MOI.OPTIMAL
             sp_is_feasible = false 
             # @logmsg LogLevel(-3) "benders_sp prob is infeasible"
@@ -279,7 +277,7 @@ function solve_sp_to_gencut!(
         benders_sp_lagrangian_bound_contrib = compute_benders_sp_lagrangian_bound_contrib(algdata, spform, optresult)
 
         primalsol = get_best_lp_primal_sol(optresult)
-        spsol_relaxed = contains(primalsol, varid -> getduty(varid) == BendSpSlackFirstStageVar)
+        spsol_relaxed = contains(primalsol, varid -> getduty(varid) <= BendSpSlackFirstStageVar)
 
         benders_sp_primal_bound_contrib = 0.0
         # compute benders_sp_primal_bound_contrib which stands for the sum of nu var,
@@ -294,12 +292,12 @@ function solve_sp_to_gencut!(
                 end
             end
         end
-        
+
         if - algo.feasibility_tol <= get_lp_primal_bound(optresult) <= algo.feasibility_tol
         # no cuts are generated since there is no violation 
             if spsol_relaxed
                 if algdata.spform_phase[spform_uid] == PurePhase2
-                    error("In PurePhase2, art var were not supposed to be in sp forlumation ")
+                    error("In PurePhase2, art var were not supposed to be in sp formulation ")
                 end
                 if algdata.spform_phase[spform_uid] == PurePhase1
                     error("In PurePhase1, if art var were in sol, the objective should be strictly positive.")
@@ -326,7 +324,6 @@ function solve_sp_to_gencut!(
                     continue
                 end             
             end
-
         else # a cut can be generated since there is a violation
             recorded_dual_solution_ids = record_solutions!(algo, algdata, spform, optresult)
             if spsol_relaxed && algo.option_increase_cost_in_hybrid_phase
@@ -418,9 +415,14 @@ end
 
 function solve_relaxed_master!(master::Formulation, env::Env)
     elapsed_time = @elapsed begin
-        optresult = TO.@timeit Coluna._to "relaxed master" run!(SolveLpForm(get_dual_solution = true), env, master, OptimizationInput(OptimizationState(master)))
+        optstate = TO.@timeit Coluna._to "relaxed master" begin
+            run!(
+                SolveLpForm(get_dual_solution = true, relax_integrality = true),
+                env, master, OptimizationInput(OptimizationState(master))
+            )
+        end
     end
-    return optresult, elapsed_time
+    return optstate, elapsed_time
 end
 
 function generatecuts!(
@@ -470,8 +472,8 @@ function bend_cutting_plane_main_loop!(
         cur_gap = 0.0
 
         optoutput, master_time = solve_relaxed_master!(masterform, env)
-        optresult = getoptstate(optoutput)
 
+        optresult = getoptstate(optoutput)
         if getterminationstatus(optresult) == INFEASIBLE
             db = - getvalue(DualBound(masterform))
             pb = - getvalue(PrimalBound(masterform))

diff --git a/src/MOIwrapper.jl b/src/MOIwrapper.jl
@@ -216,10 +216,10 @@ end
 
 function MOI.get(model::Coluna.Optimizer, ::MOI.ListOfVariableIndices)
     indices = Vector{MathOptInterface.VariableIndex}()
-    for (key,value) in model.moi_varids
+    for (_, value) in model.moi_varids
         push!(indices, value)
     end
-    return sort!(indices)
+    return sort!(indices, by = x -> x.value)
 end
 
 ############################################################################################
@@ -286,11 +286,10 @@ end
 function MOI.get(
     model::Coluna.Optimizer, ::MOI.ListOfConstraintIndices{F, S}
 ) where {F<:MOI.SingleVariable, S}
-    orig_form = get_original_formulation(model.inner)
     indices = MOI.ConstraintIndex{F,S}[]
-    for (id, var) in model.vars
-        if S == MathProg.convert_coluna_kind_to_moi(getperenkind(orig_form, var))
-            push!(indices, MOI.ConstraintIndex{F,S}(id.value))
+    for (id, _) in model.constrs_on_single_var_to_vars
+        if S == typeof(MOI.get(model, MOI.ConstraintSet(), id))
+            push!(indices, id)
         end
     end
     return sort!(indices, by = x -> x.value)
@@ -302,8 +301,8 @@ function MOI.get(
     orig_form = get_original_formulation(model.inner)
     constrid = getid(model.constrs[index])
     terms = MOI.ScalarAffineTerm{Float64}[]
-    for (varid, coeff) in @view getcoefmatrix(orig_form)[constrid, :]
-        push!(terms, MOI.ScalarAffineTerm(coeff, model.moi_varids[varid]))
+    for (varid, coef) in @view getcoefmatrix(orig_form)[constrid, :]
+        push!(terms, MOI.ScalarAffineTerm(coef, model.moi_varids[varid]))
     end
     return MOI.ScalarAffineFunction(terms, 0.0)
 end
@@ -425,6 +424,15 @@ end
 ############################################################################################
 # Attributes of constraints
 ############################################################################################
+function MOI.set(
+    model::MOI.Bridges.LazyBridgeOptimizer{Coluna.Optimizer}, attr::MOI.AbstractConstraintAttribute,
+    bridge::MOI.Bridges.Constraint.SplitIntervalBridge, value
+)
+    MOI.set(model.model, attr, bridge.lower, value)
+    MOI.set(model.model, attr, bridge.upper, value)
+    return
+end
+
 function MOI.set(
     model::Coluna.Optimizer, ::BD.ConstraintDecomposition, constrid::MOI.ConstraintIndex,
     annotation::BD.Annotation
@@ -610,6 +618,8 @@ function MOI.empty!(model::Coluna.Optimizer)
     model.env.varids = CleverDicts.CleverDict{MOI.VariableIndex, VarId}()
     model.moi_varids = Dict{VarId, MOI.VariableIndex}()
     model.constrs = Dict{MOI.ConstraintIndex, Constraint}()
+    model.constrs_on_single_var_to_vars = Dict{MOI.ConstraintIndex, VarId}()
+    model.constrs_on_single_var_to_names = Dict{MOI.ConstraintIndex, String}()
     if model.default_optimizer_builder !== nothing
         set_default_optimizer_builder!(model.inner, model.default_optimizer_builder)
     end
@@ -625,11 +635,13 @@ mutable struct ColumnInfo <: BD.AbstractColumnInfo
 end
 
 function BD.getsolutions(model::Coluna.Optimizer, k)
-    ip_primal_sol = model.disagg_result.ip_primal_sols[k]
+    ip_primal_sol = get_best_ip_primal_sol(model.disagg_result)
     sp_columns_info = Vector{ColumnInfo}()
     for (varid, val) in ip_primal_sol
         if getduty(varid) <= MasterCol
-            push!(sp_columns_info, ColumnInfo(model, varid, val))
+            if  model.annotations.ann_per_form[getoriginformuid(varid)].axis_index_value == k
+                push!(sp_columns_info, ColumnInfo(model, varid, val))
+            end
         end
     end
     return sp_columns_info
@@ -694,18 +706,21 @@ function MOI.get(optimizer::Optimizer, ::MOI.ObjectiveValue)
     return getvalue(get_ip_primal_bound(optimizer.result))
 end
 
+function MOI.get(optimizer::Optimizer, ::MOI.DualObjectiveValue)
+    return getvalue(get_lp_dual_bound(optimizer.result))
+end
+
 function MOI.get(optimizer::Optimizer, ::MOI.RelativeGap)
     return ip_gap(optimizer.result)
 end
 
-function MOI.get(optimizer::Optimizer, ::MOI.VariablePrimal, ref::MOI.VariableIndex)
+function MOI.get(optimizer::Optimizer, attr::MOI.VariablePrimal, ref::MOI.VariableIndex)
     id = getid(optimizer.vars[ref]) # This gets a coluna Id{Variable}
-    best_primal_sol = get_best_ip_primal_sol(optimizer.result)
-    if best_primal_sol === nothing
-        @warn "Coluna did not find a primal feasible solution."
-        return NaN
+    primalsols = get_ip_primal_sols(optimizer.result)
+    if 1 <= attr.N <= length(primalsols)
+        return get(primalsols[attr.N], id, 0.0)
     end
-    return get(best_primal_sol, id, 0.0)
+    return error("Invalid result index.")
 end
 
 function MOI.get(optimizer::Optimizer, ::MOI.VariablePrimal, refs::Vector{MOI.VariableIndex})
@@ -766,7 +781,20 @@ end
 MOI.get(optimizer::Optimizer, ::MOI.NodeCount) = optimizer.env.kpis.node_count
 MOI.get(optimizer::Optimizer, ::MOI.SolveTime) = optimizer.env.kpis.elapsed_optimization_time
 
-# function MOI.get(optimizer::Optimizer, ::MOI.ConstraintDual, index::MOI.ConstraintIndex)
+function MOI.get(
+    optimizer::Optimizer, attr::MOI.ConstraintDual, 
+    index::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}}
+)
+    dualsols = get_lp_dual_sols(optimizer.result)
+    if 1 <= attr.N <= length(dualsols)
+        return get(dualsols[attr.N], getid(optimizer.constrs[index]), 0.0)
+    end
+    return error("Invalid result index.")
+end
+
+# function MOI.get(
+#     optimizer::Optimizer, attr::MOI.ConstraintDual, index::MOI.ConstraintIndex{F,S}
+# ) where {F<:MOI.SingleVariable,S}
 #     return 0.0
 # end
 

diff --git a/src/MathProg/MOIinterface.jl b/src/MathProg/MOIinterface.jl
@@ -275,7 +275,6 @@ function get_dual_solutions(form::F, optimizer::MoiOptimizer) where {F <: Formul
             solcost += val * getcurrhs(form, id)
             val = round(val, digits = Coluna.TOL_DIGITS)
             if abs(val) > Coluna.TOL
-                @logmsg LogLevel(-4) string("Constr ", constr.name, " = ", val)
                 push!(solconstrs, id)
                 push!(solvals, val)      
             end

diff --git a/src/MathProg/clone.jl b/src/MathProg/clone.jl
@@ -1,3 +1,4 @@
+# TODO : these methods should not be part of MathProg.
 function clonevar!(
     originform::Formulation,
     destform::Formulation,

diff --git a/src/MathProg/duties.jl b/src/MathProg/duties.jl
@@ -23,7 +23,11 @@ mutable struct DwSp <: AbstractSpDuty
 end
 
 "A Benders subproblem of a formulation decomposed using Benders."
-struct BendersSp <: AbstractSpDuty end
+struct BendersSp <: AbstractSpDuty 
+    slack_to_first_stage::Dict{VarId, VarId}
+end
+
+BendersSp() = BendersSp(Dict{VarId, VarId}())
 
 #
 # Duties tree for a Variable
@@ -53,6 +57,8 @@ struct BendersSp <: AbstractSpDuty end
         AbstractBendSpVar <= Duty{Variable}
             AbstractBendSpSlackMastVar <= AbstractBendSpVar
                 BendSpSlackFirstStageVar <= AbstractBendSpSlackMastVar
+                    BendSpPosSlackFirstStageVar <= BendSpSlackFirstStageVar
+                    BendSpNegSlackFirstStageVar <= BendSpSlackFirstStageVar
                 BendSpSlackSecondStageCostVar <= AbstractBendSpSlackMastVar
             BendSpSepVar <= AbstractBendSpVar
             #BendSpPureVar <= AbstractBendSpVar