Better support of single variable constraints (#590)

* wip, improve support of singlvar constraints

* wip

* wip

* tests ok

* wip

* wip

* wip

* v1

* clean Project

* wip

* v2

* clean unused methods

* add tests

* wip

* solutions.jl ok

* more tests for constraints

* ok?

* rm show form

* ok

src/Algorithm/benders.jl

@@ -1,5 +1,4 @@
 @with_kw struct BendersCutGeneration <: AbstractOptimizationAlgorithm
-    option_use_reduced_cost::Bool = false
     option_increase_cost_in_hybrid_phase::Bool = false
     feasibility_tol::Float64 = 1e-5
     optimality_tol::Float64 = Coluna.DEF_OPTIMALITY_ATOL
@@ -110,15 +109,16 @@ function update_benders_sp_problem!(
         setcurub!(spform, var, getperenub(spform, var) - master_primal_sol[varid])
     end
 
-    if algo.option_use_reduced_cost
-        for (varid, var) in getvars(spform)
-            iscuractive(spform, varid) || continue
-            getduty(varid) <= BendSpSlackFirstStageVar || continue
-            cost = getcurcost(spform, var)
-            rc = computereducedcost(masterform, varid, master_dual_sol)
-            setcurcost!(spform, var, rc)
-        end
-    end
+    # TODO : it was untested
+    # if algo.option_use_reduced_cost
+    #     for (varid, var) in getvars(spform)
+    #         iscuractive(spform, varid) || continue
+    #         getduty(varid) <= BendSpSlackFirstStageVar || continue
+    #         cost = getcurcost(spform, var)
+    #         rc = computereducedcost(masterform, varid, master_dual_sol)
+    #         setcurcost!(spform, var, rc)
+    #     end
+    # end
     return false
 end
 
@@ -190,7 +190,7 @@ function insert_cuts_in_master!(
         name = string("BC_", getsortuid(dual_sol_id))
         kind = Essential
         duty = MasterBendCutConstr
-        bc = setcut_from_sp_dualsol!(
+        setcut_from_sp_dualsol!(
             masterform,
             spform,
             dual_sol_id,
@@ -233,7 +233,6 @@ function solve_sp_to_gencut!(
     benders_sp_primal_bound_contrib = 0.0
     benders_sp_lagrangian_bound_contrib =  0.0
 
-    insertion_status = 0
     spsol_relaxed = false
 
     # Compute target
@@ -363,7 +362,7 @@ function solve_sps_to_gencuts!(
 
     ### BEGIN LOOP TO BE PARALLELIZED
     for (spuid, spform) in sps
-        recorded_sp_dual_solution_ids[spuid] = Vector{ConstrId}()
+        recorded_sp_dual_solution_ids[spuid] = ConstrId[]
         gen_status, spsol_relaxed, recorded_dual_solution_ids, benders_sp_primal_bound_contrib, benders_sp_lagrangian_bound_contrib = solve_sp_to_gencut!(
             algo, env, algdata, masterform, spform,
             master_primalsol, master_dualsol,
@@ -559,7 +558,6 @@ function bend_cutting_plane_main_loop!(
             setterminationstatus!(bnd_optstate, OPTIMAL)
             break # loop on master lp solution          
         end
-
     end  # loop on master lp solution 
 
     if !one_spsol_is_a_relaxed_sol                

src/Algorithm/colgen.jl

@@ -342,7 +342,7 @@ function solve_sp_to_gencol!(
 
     sense = getobjsense(masterform)
     bestsol = get_best_ip_primal_sol(sp_optstate)
-    if bestsol !== nothing && bestsol.status == FEASIBLE_SOL
+    if bestsol !== nothing && getstatus(bestsol) == FEASIBLE_SOL
         spinfo.bestsol = bestsol
         spinfo.isfeasible = true
         for sol in get_ip_primal_sols(sp_optstate)
@@ -380,7 +380,7 @@ function updatereducedcosts!(
     reform::Reformulation, redcostshelper::ReducedCostsCalculationHelper, masterdualsol::DualSolution
 )
     redcosts = Dict{VarId,Float64}()
-    result = redcostshelper.dwsprep_coefmatrix * getsol(masterdualsol)
+    result = redcostshelper.dwsprep_coefmatrix * masterdualsol.solution.sol
     for (i, varid) in enumerate(redcostshelper.dwspvarids)
         redcosts[varid] = redcostshelper.perencosts[i] - get(result, varid, 0.0)
     end
@@ -524,7 +524,7 @@ function update_lagrangian_dual_bound!(
         end
     end
 
-    valid_lagr_bound = DualBound{S}(puremastvars_contrib + dualsol.bound)
+    valid_lagr_bound = DualBound{S}(puremastvars_contrib + getbound(dualsol))
     for (spuid, spinfo) in spinfos
         valid_lagr_bound += spinfo.valid_dual_bound_contrib
     end
@@ -533,7 +533,7 @@ function update_lagrangian_dual_bound!(
     update_lp_dual_bound!(optstate, valid_lagr_bound)
 
     if stabilization_is_used(algo)
-        pseudo_lagr_bound = DualBound{S}(puremastvars_contrib + dualsol.bound)
+        pseudo_lagr_bound = DualBound{S}(puremastvars_contrib + getbound(dualsol))
         for (spuid, spinfo) in spinfos
             pseudo_lagr_bound += spinfo.pseudo_dual_bound_contrib
         end
@@ -560,9 +560,9 @@ function compute_subgradient_contribution(
             end
         end
 
-        for (spuid, spinfo) in spinfos
+        for (_, spinfo) in spinfos
             iszero(spinfo.ub) && continue
-            mult = improving_red_cost(spinfo.bestsol.bound, algo, sense) ? spinfo.ub : spinfo.lb
+            mult = improving_red_cost(getbound(spinfo.bestsol), algo, sense) ? spinfo.ub : spinfo.lb
             for (sp_var_id, sp_var_val) in spinfo.bestsol
                 for (master_constrid, sp_var_coef) in @view master_coef_matrix[:,sp_var_id]
                     if !(getduty(master_constrid) <= MasterConvexityConstr)
@@ -574,13 +574,16 @@ function compute_subgradient_contribution(
         end
     end
 
-    return DualSolution(master, constrids, constrvals, 0.0, UNKNOWN_SOLUTION_STATUS)
+    return DualSolution(
+        master, constrids, constrvals, VarId[], Float64[], ActiveBound[], 0.0, 
+        UNKNOWN_SOLUTION_STATUS
+    )
 end
 
 function move_convexity_constrs_dual_values!(
     spinfos::Dict{FormId,SubprobInfo}, dualsol::DualSolution
 )
-    newbound = dualsol.bound
+    newbound = getbound(dualsol)
     for (spuid, spinfo) in spinfos
         spinfo.lb_dual = dualsol[spinfo.lb_constr_id]
         spinfo.ub_dual = dualsol[spinfo.ub_constr_id]
@@ -594,7 +597,10 @@ function move_convexity_constrs_dual_values!(
             push!(values, value)
         end
     end
-    return DualSolution(dualsol.model, constrids, values, newbound, FEASIBLE_SOL)
+    return DualSolution(
+        getmodel(dualsol), constrids, values, VarId[], Float64[], ActiveBound[], newbound, 
+        FEASIBLE_SOL
+    )
 end
 
 function get_pure_master_vars(master::Formulation)
@@ -803,8 +809,9 @@ function cg_main_loop!(
         if nb_new_columns == 0 && !essential_cuts_separated
             @logmsg LogLevel(0) "No new column generated by the pricing problems."
             setterminationstatus!(cg_optstate, OTHER_LIMIT)
-            # if no columns are generated and lp gap is not closed then this col.gen. stage
+            # If no columns are generated and lp gap is not closed then this col.gen. stage
             # is a heuristic one, so we do not run phase 1 to save time
+            # Comment by @guimarqu : It may also be a bug
             return true, false
         end
         if iteration > algo.max_nb_iterations

src/Algorithm/colgenstabilization.jl

@@ -49,7 +49,7 @@ function init_stab_before_colgen_loop!(unit::ColGenStabilizationUnit)
 end
 
 function componentwisefunction(in_dual_sol::DualSolution, out_dual_sol::DualSolution, f::Function)
-    form = out_dual_sol.model
+    form = getmodel(out_dual_sol)
     constrids = Vector{ConstrId}()
     constrvals = Vector{Float64}()
     value::Float64 = 0.0
@@ -100,12 +100,15 @@ function linear_combination(in_dual_sol::DualSolution, out_dual_sol::DualSolutio
         (x, y) -> coeff * x + (1.0 - coeff) * y
     )
 
-    form = in_dual_sol.model
+    form = getmodel(in_dual_sol)
     bound = 0.0
     for (i, constrid) in enumerate(constrids)
         bound += constrvals[i] * getcurrhs(form, constrid) 
     end
-    return DualSolution(form, constrids, constrvals, bound, UNKNOWN_FEASIBILITY)
+    return DualSolution(
+        form, constrids, constrvals, VarId[], Float64[], ActiveBound[], bound, 
+        UNKNOWN_FEASIBILITY
+    )
 end
 
 function update_stab_after_rm_solve!(
@@ -148,11 +151,14 @@ function update_alpha_automatically!(
     smooth_dual_sol::DualSolution{M}, subgradient_contribution::DualSolution{M}
 ) where {M}    
 
-    master = lp_dual_sol.model 
+    master = getmodel(lp_dual_sol)
 
     # first we calculate the in-sep direction
     constrids, constrvals = componentwisefunction(smooth_dual_sol, unit.stabcenter, -)
-    in_sep_direction = DualSolution(master, constrids, constrvals, 0.0, UNKNOWN_FEASIBILITY)
+    in_sep_direction = DualSolution(
+        master, constrids, constrvals, VarId[], Float64[], ActiveBound[], 0.0, 
+        UNKNOWN_FEASIBILITY
+    )
     in_sep_dir_norm = norm(in_sep_direction)
 
     # we initialize the subgradient with the right-hand-side of all master constraints
@@ -166,7 +172,10 @@ function update_alpha_automatically!(
             push!(constrrhs, getcurrhs(master, constrid))
         end 
     end
-    subgradient = DualSolution(master, constrids, constrrhs, 0.0, UNKNOWN_FEASIBILITY)
+    subgradient = DualSolution(
+        master, constrids, constrrhs, VarId[], Float64[], ActiveBound[], 0.0, 
+        UNKNOWN_FEASIBILITY
+    )
 
     # we calculate the subgradient at the sep point
     for (constrid, value) in subgradient_contribution

src/Algorithm/node.jl

@@ -8,7 +8,7 @@ mutable struct Node
     depth::Int
     parent::Union{Nothing, Node}
     optstate::OptimizationState
-    #branch::Union{Nothing, Branch} # branch::Id{Constraint}
+    #branch::Union{Nothing, Branch} # branch::ConstrId
     branchdescription::String
     recordids::RecordsVector
     conquerwasrun::Bool

src/Algorithm/utilities/optimizationstate.jl

@@ -1,6 +1,3 @@
-getvalue(sol::Solution) = ColunaBase.getvalue(sol)
-getvalue(bnd::Bound) = ColunaBase.getvalue(bnd)
-
 mutable struct OptimizationState{F<:AbstractFormulation,S<:Coluna.AbstractSense}
     termination_status::TerminationStatus
     incumbents::ObjValues{S}
@@ -15,7 +12,7 @@ mutable struct OptimizationState{F<:AbstractFormulation,S<:Coluna.AbstractSense}
     lp_dual_sols::Vector{DualSolution{F}}
 end
 
-function bestbound!(solutions::Vector{Sol}, max_len::Int, new_sol::Sol) where {Sol<:Solution}
+function bestbound!(solutions::Vector{Sol}, max_len::Int, new_sol::Sol) where {Sol<:AbstractSolution}
     push!(solutions, new_sol)
     sort!(solutions, rev = true)
     while length(solutions) > max_len
@@ -24,7 +21,7 @@ function bestbound!(solutions::Vector{Sol}, max_len::Int, new_sol::Sol) where {S
     return
 end
 
-function set!(solutions::Vector{Sol}, ::Int, new_sol::Sol) where {Sol<:Solution}
+function set!(solutions::Vector{Sol}, ::Int, new_sol::Sol) where {Sol<:AbstractSolution}
     empty!(solutions)
     push!(solutions, new_sol)
     return

src/ColunaBase/ColunaBase.jl

@@ -19,8 +19,8 @@ export AbstractModel, AbstractProblem, AbstractSense, AbstractMinSense, Abstract
 export NestedEnum, @nestedenum, @exported_nestedenum
 
 # solsandbounds.jl
-export Bound, Solution, getvalue, isbetter, diff, gap, printbounds, getsol,
-    getstatus, remove_until_last_point
+export Bound, Solution, getvalue, getbound, isbetter, diff, gap, printbounds,
+    getstatus, remove_until_last_point, getmodel
 
 # Statuses
 export TerminationStatus, SolutionStatus, OPTIMIZE_NOT_CALLED, OPTIMAL,

src/ColunaBase/solsandbounds.jl

@@ -214,23 +214,24 @@ function convert_status(coluna_status::SolutionStatus)
     return MOI.OTHER_RESULT_STATUS
 end
 
-# Solution
+# Basic structure of a solution
 struct Solution{Model<:AbstractModel,Decision,Value} <: AbstractDict{Decision,Value}
     model::Model
     bound::Float64
     status::SolutionStatus
     sol::DynamicSparseArrays.PackedMemoryArray{Decision,Value}
-    custom_data::Union{Nothing, BlockDecomposition.AbstractCustomData}
 end
 
 """
+Solution is an internal data structure of Coluna and should not be used in
+algorithms. See `MathProg.PrimalSolution` & `MathProg.DualSolution` instead.
+
     Solution(
         model::AbstractModel,
         decisions::Vector,
         values::Vector,
-        solution_values::Float64,
-        status::SolutionStatus,
-        [custom_data]::Union{Nothing, BlockDecomposition.AbstractCustomData}
+        solution_value::Float64,
+        status::SolutionStatus
     )
 
 Create a solution to the `model`. Other arguments are: 
@@ -241,42 +242,34 @@ Create a solution to the `model`. Other arguments are:
 """
 function Solution{Mo,De,Va}(
     model::Mo, decisions::Vector{De}, values::Vector{Va}, solution_value::Float64, 
-    status::SolutionStatus, custom_data::Union{Nothing, BlockDecomposition.AbstractCustomData} = nothing
-) where {Mo<:AbstractModel,De,Va}
+    status::SolutionStatus
+) where {Mo<:AbstractModel,De,Va,T}
     sol = DynamicSparseArrays.dynamicsparsevec(decisions, values)
-    return Solution(model, solution_value, status, sol, custom_data)
+    return Solution(model, solution_value, status, sol)
 end
 
-"""
-    getsol(solution)
-    
-Return the dynamic sparse vector that describes `solution`.
-"""
-getsol(s::Solution) = s.sol
+"Return the model of a solution."
+getmodel(s::Solution) = s.model
 
-"""
-    getvalue(solution) -> Float64
+"Return the value (as a Bound) of `solution`"
+getbound(s::Solution) = s.bound
 
-Return the value of `solution`.
-"""
+"Return the value of `solution`."
 getvalue(s::Solution) = float(s.bound)
 
-"""
-    getstatus(solution) -> SolutionStatus
-
-Return the solution status of `solution`.
-"""
+"Return the solution status of `solution`."
 getstatus(s::Solution) = s.status
 
 Base.iterate(s::Solution) = iterate(s.sol)
 Base.iterate(s::Solution, state) = iterate(s.sol, state)
 Base.length(s::Solution) = length(s.sol)
-Base.get(s::Solution{Mo,De,Va}, id::De, default) where {Mo,De,Va} = s.sol[id] # TODO : REMOVE
+Base.get(s::Solution{Mo,De,Va}, id::De, default) where {Mo,De,Va} = s.sol[id]
 Base.getindex(s::Solution{Mo,De,Va}, id::De) where {Mo,De,Va} = Base.getindex(s.sol, id)
 Base.setindex!(s::Solution{Mo,De,Va}, val::Va, id::De) where {Mo,De,Va} = s.sol[id] = val
 
+# TODO : remove when refactoring Benders
 function Base.filter(f::Function, s::S) where {S <: Solution}
-    return S(s.model, s.bound, s.status, filter(f, s.sol), s.custom_data)
+    return S(s.model, s.bound, s.status, filter(f, s.sol))
 end
 
 function Base.in(p::Tuple{De,Va}, a::Solution{Mo,De,Va}, valcmp=(==)) where {Mo,De,Va}

src/MOIcallbacks.jl

@@ -34,7 +34,10 @@ function MOI.submit(
     push!(values, 1.0)
     solval += getcurcost(form, setup_var_id)
 
-    sol = PrimalSolution(form, colunavarids, values, solval, FEASIBLE_SOL, custom_data)
+    sol = PrimalSolution(
+        form, colunavarids, values, solval, FEASIBLE_SOL; 
+        custom_data = custom_data
+    )
     push!(cb.callback_data.primal_solutions, sol)
     return
 end