Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Better support of single variable constraints #590

Merged
merged 22 commits into from
Oct 13, 2021
Merged
26 changes: 12 additions & 14 deletions src/Algorithm/benders.jl
Original file line number Diff line number Diff line change
@@ -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
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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
Expand Down
27 changes: 17 additions & 10 deletions src/Algorithm/colgen.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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)
Expand All @@ -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]
Expand All @@ -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)
Expand Down Expand Up @@ -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
Expand Down
21 changes: 15 additions & 6 deletions src/Algorithm/colgenstabilization.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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!(
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down
2 changes: 1 addition & 1 deletion src/Algorithm/node.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This can be deleted

branchdescription::String
recordids::RecordsVector
conquerwasrun::Bool
Expand Down
7 changes: 2 additions & 5 deletions src/Algorithm/utilities/optimizationstate.jl
Original file line number Diff line number Diff line change
@@ -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}
Expand All @@ -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
Expand All @@ -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
Expand Down
4 changes: 2 additions & 2 deletions src/ColunaBase/ColunaBase.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
43 changes: 18 additions & 25 deletions src/ColunaBase/solsandbounds.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand All @@ -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}
Expand Down
5 changes: 4 additions & 1 deletion src/MOIcallbacks.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Loading