Refactoring ObjValues & OptimizationState

.github/workflows/ci.yml

@@ -27,6 +27,8 @@ jobs:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
       - 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"));'
       - uses: actions/cache@v1
         env:
           cache-name: cache-artifacts

src/Algorithm/Algorithm.jl

@@ -14,13 +14,6 @@ const MOI = MathOptInterface
 
 import Base: push!
 
-# Import to extend methods to OptimizationState
-import ..MathProg:
-    get_ip_primal_bound, get_ip_dual_bound,
-    get_lp_primal_bound, get_lp_dual_bound, update_ip_primal_bound!, update_ip_dual_bound!,
-    update_lp_primal_bound!, update_lp_dual_bound!, set_ip_primal_bound!,
-    set_ip_dual_bound!, set_lp_primal_bound!, set_lp_dual_bound!, ip_gap, lp_gap, ip_gap_closed, lp_gap_closed
-
 # Utilities to build algorithms
 include("utilities/optimizationstate.jl")
 
@@ -60,17 +53,19 @@ include("treesearch.jl")
 
 # Utilities
 export getterminationstatus, setterminationstatus!,
-    nb_ip_primal_sols, nb_lp_primal_sols, nb_lp_dual_sols,
     get_ip_primal_sols, get_lp_primal_sols, get_lp_dual_sols, get_best_ip_primal_sol,
     get_best_lp_primal_sol, get_best_lp_dual_sol, update_ip_primal_sol!,
-    update_lp_primal_sol!, update_lp_dual_sol!, add_ip_primal_sol!, add_lp_primal_sol!,
-    add_lp_dual_sol!, set_ip_primal_sol!, set_lp_primal_sol!, set_lp_dual_sol!,
-    get_ip_dual_bound, set_ip_dual_bound!, update_all_ip_primal_solutions!, getreform,
+    update_lp_primal_sol!, update_lp_dual_sol!, add_ip_primal_sol!, add_ip_primal_sols!,
+    add_lp_primal_sol!, add_lp_dual_sol!, set_ip_primal_sol!, set_lp_primal_sol!, set_lp_dual_sol!,
+    empty_ip_primal_sols!, empty_lp_primal_sols!, empty_lp_dual_sols!,
+    get_ip_primal_bound, get_lp_primal_bound, get_lp_dual_bound, get_ip_dual_bound, 
+    set_ip_primal_bound!, set_lp_primal_bound!, set_lp_dual_bound!, set_ip_dual_bound!,
+    update_ip_primal_bound!, update_lp_primal_bound!, update_lp_dual_bound!, update_ip_dual_bound!,
     getoptstate, run!, isinfeasible
 
 # Algorithm's types
 export AbstractOptimizationAlgorithm, TreeSearchAlgorithm, ColCutGenConquer, ColumnGeneration,
-       BendersConquer, BendersCutGeneration, SolveIpForm, RestrictedMasterHeuristic,
+       BendersConquer, BendersCutGeneration, SolveIpForm, RestrictedMasterIPHeuristic,
        SolveLpForm, ExactBranchingPhase, OnlyRestrictedMasterBranchingPhase, PreprocessAlgorithm, 
        OptimizationInput, OptimizationOutput, OptimizationState, EmptyInput
 

src/Algorithm/benders.jl

@@ -42,8 +42,9 @@ end
 
 function BendersCutGenRuntimeData(form::Reformulation, init_optstate::OptimizationState)
     optstate = OptimizationState(getmaster(form))
-    if nb_ip_primal_sols(init_optstate) > 0
-        add_ip_primal_sol!(optstate, get_best_ip_primal_sol(init_optstate))
+    best_ip_primal_sol = get_best_ip_primal_sol(init_optstate)
+    if best_ip_primal_sol !== nothing
+        add_ip_primal_sol!(optstate, best_ip_primal_sol)
     end
     return BendersCutGenRuntimeData(optstate, Dict{FormId, FormulationPhase}(), Dict{FormId, Bool}())#0.0, true)
 end

src/Algorithm/branching/branchingalgo.jl

@@ -148,13 +148,19 @@ function perform_strong_branching_with_phases!(
                     @printf "), value = %6.2f\n" getvalue(get_lp_primal_bound(getoptstate(node)))
                 end
 
-                update_ip_primal!(getoptstate(node), sbstate, exploitsprimalsolutions)
+                nodestate = getoptstate(node)
+
+                update_ip_primal_bound!(nodestate, get_ip_primal_bound(sbstate))
+                best_ip_primal_sol = get_best_ip_primal_sol(sbstate)
+                if exploitsprimalsolutions && best_ip_primal_sol !== nothing
+                    set_ip_primal_sol!(nodestate, best_ip_primal_sol)
+                end                
 
                 apply_conquer_alg_to_node!(
                     node, current_phase.conquer_algo, env, reform, conquer_units_to_restore
                 )        
 
-                update_all_ip_primal_solutions!(sbstate, getoptstate(node))
+                add_ip_primal_sols!(sbstate, get_ip_primal_sols(nodestate)...)
 
                 if to_be_pruned(node) 
                     if isverbose(current_phase.conquer_algo)
@@ -209,10 +215,9 @@ function run!(algo::StrongBranching, env::Env, reform::Reformulation, input::Div
     # we obtain the original and extended solutions
     master = getmaster(reform)
     original_solution = nothing
-    extended_solution = nothing
-    if nb_lp_primal_sols(optstate) > 0
+    extended_solution = get_best_lp_primal_sol(optstate)
+    if extended_solution !== nothing
         if projection_is_possible(master)
-            extended_solution = get_best_lp_primal_sol(optstate)
             original_solution = proj_cols_on_rep(extended_solution, master)
         else
             original_solution = get_best_lp_primal_sol(optstate)

src/Algorithm/branching/branchinggroup.jl

@@ -59,10 +59,9 @@ function regenerate_children!(group::BranchingGroup, parent::Node)
 end
 
 function product_score(group::BranchingGroup, parent_optstate::OptimizationState)
-    parent_inc = getincumbents(parent_optstate)
     # TO DO : we need to mesure the gap to the cut-off value
-    parent_lp_dual_bound = get_lp_dual_bound(parent_inc)
-    parent_delta = diff(get_ip_primal_bound(parent_inc), parent_lp_dual_bound)
+    parent_lp_dual_bound = get_lp_dual_bound(parent_optstate)
+    parent_delta = diff(get_ip_primal_bound(parent_optstate), parent_lp_dual_bound)
 
     all_branches_above_delta = true
     deltas = zeros(Float64, length(group.children))
@@ -126,11 +125,9 @@ function tree_depth_score(group::BranchingGroup, parent_optstate::OptimizationSt
         return 0.0
     end
 
-    parent_inc = getincumbents(parent_optstate)
-
     # TO DO : we need to mesure the gap to the cut-off value
-    parent_lp_dual_bound = get_lp_dual_bound(parent_inc)
-    parent_delta = diff(get_ip_primal_bound(parent_inc), parent_lp_dual_bound)
+    parent_lp_dual_bound = get_lp_dual_bound(parent_optstate)
+    parent_delta = diff(get_ip_primal_bound(parent_optstate), parent_lp_dual_bound)
 
     deltas = zeros(Float64, length(group.children))
     nb_zero_deltas = 0

src/Algorithm/colgen.jl

@@ -319,27 +319,24 @@ function solve_sp_to_gencol!(
     compute_db_contributions!(spinfo, get_ip_dual_bound(sp_optstate), sp_sol_value)
 
     sense = getobjsense(masterform)
-    if nb_ip_primal_sols(sp_optstate) > 0
-        bestsol = get_best_ip_primal_sol(sp_optstate)
-
-        if bestsol.status == FEASIBLE_SOL
-            spinfo.bestsol = bestsol
-            spinfo.isfeasible = true
-            for sol in get_ip_primal_sols(sp_optstate)
-                if improving_red_cost(compute_red_cost(algo, masterform, spinfo, sol, dualsol), algo, sense)
-                    insertion_status, col_id = setprimalsol!(spform, sol)
-                    if insertion_status
-                        push!(spinfo.recorded_sol_ids, col_id)
-                    elseif !insertion_status && haskey(masterform, col_id) && !iscuractive(masterform, col_id)
-                        push!(spinfo.sol_ids_to_activate, col_id)
-                    elseif !insertion_status && !haskey(masterform, col_id)
-                        @warn "The pricing problem generated the column with id $(col_id) twice."
-                    else
-                        msg = """
-                        Column already exists as $(getname(masterform, col_id)) and is already active.
-                        """
-                        @warn string(msg)
-                    end
+    bestsol = get_best_ip_primal_sol(sp_optstate)
+    if bestsol !== nothing && bestsol.status == FEASIBLE_SOL
+        spinfo.bestsol = bestsol
+        spinfo.isfeasible = true
+        for sol in get_ip_primal_sols(sp_optstate)
+            if improving_red_cost(compute_red_cost(algo, masterform, spinfo, sol, dualsol), algo, sense)
+                insertion_status, col_id = setprimalsol!(spform, sol)
+                if insertion_status
+                    push!(spinfo.recorded_sol_ids, col_id)
+                elseif !insertion_status && haskey(masterform, col_id) && !iscuractive(masterform, col_id)
+                    push!(spinfo.sol_ids_to_activate, col_id)
+                elseif !insertion_status && !haskey(masterform, col_id)
+                    @warn "The pricing problem generated the column with id $(col_id) twice."
+                else
+                    msg = """
+                    Column already exists as $(getname(masterform, col_id)) and is already active.
+                    """
+                    @warn string(msg)
                 end
             end
         end
@@ -351,7 +348,6 @@ function updatereducedcosts!(
     reform::Reformulation, redcostshelper::ReducedCostsCalculationHelper, dualsol::DualSolution
 )
     redcosts = deepcopy(redcostshelper.perencosts)
-    master = getmaster(reform)
 
     result = transpose(redcostshelper.dwsprep_coefmatrix) * getsol(dualsol)
 
@@ -642,10 +638,8 @@ function cg_main_loop!(
             return true, false
         end
 
-        lp_dual_sol = nothing
-        if nb_lp_dual_sols(rm_optstate) > 0
-            lp_dual_sol = get_best_lp_dual_sol(rm_optstate)
-        else
+        lp_dual_sol = get_best_lp_dual_sol(rm_optstate)
+        if lp_dual_sol === nothing
             error("Solver returned that the LP restricted master is feasible but ",
             "did not return a dual solution. ",
             "Please open an issue (https://github.com/atoptima/Coluna.jl/issues).")
@@ -657,9 +651,8 @@ function cg_main_loop!(
         lp_dual_sol = move_convexity_constrs_dual_values!(spinfos, lp_dual_sol)
 
         TO.@timeit Coluna._to "Getting primal solution" begin
-        if nb_lp_primal_sols(rm_optstate) > 0
-            rm_sol = get_best_lp_primal_sol(rm_optstate)
-
+        rm_sol = get_best_lp_primal_sol(rm_optstate)
+        if rm_sol !== nothing
             set_lp_primal_sol!(cg_optstate, rm_sol)
             lp_bound = get_lp_primal_bound(rm_optstate) + getvalue(partial_solution)
             set_lp_primal_bound!(cg_optstate, lp_bound)

src/Algorithm/conquer.jl

@@ -54,7 +54,7 @@ function apply_conquer_alg_to_node!(
         @logmsg LogLevel(-1) string("Tree IP PB: ", get_ip_primal_bound(nodestate))
     end
     if ip_gap_closed(nodestate, rtol = opt_rtol, atol = opt_atol)
-        @info "IP Gap is closed: $(ip_gap(getincumbents(nodestate))). Abort treatment."
+        @info "IP Gap is closed: $(ip_gap(nodestate)). Abort treatment."
     else
         isverbose(algo) && @logmsg LogLevel(-1) string("IP Gap is positive. Need to treat node.")
 
@@ -69,6 +69,13 @@ end
 #                      ParameterisedHeuristic
 ####################################################################
 
+"""
+    Coluna.Algorithm.RestrictedMasterHeuristic()
+
+This algorithm enforces integrality of column variables in the master formulation and then solves the master formulation with its optimizer.
+"""
+RestrictedMasterIPHeuristic() = SolveIpForm(moi_params = MoiOptimize(get_dual_bound = false))
+
 struct ParameterisedHeuristic
     algorithm::AbstractOptimizationAlgorithm
     root_priority::Float64
@@ -78,14 +85,9 @@ struct ParameterisedHeuristic
     name::String
 end
 
-"""
-    Coluna.Algorithm.RestrictedMasterHeuristic()
-
-This algorithm enforces integrality of column variables in the master formulation and then solves the master formulation with its optimizer.
-"""
-RestrictedMasterHeuristic() = 
+ParamRestrictedMasterHeuristic() = 
     ParameterisedHeuristic(
-        SolveIpForm(moi_params = MoiOptimize(get_dual_bound = false)), 
+        RestrictedMasterIPHeuristic(), 
         1.0, 1.0, 1, 1000, "Restricted Master IP"
     )
 
@@ -124,7 +126,7 @@ end
 """
     Coluna.Algorithm.ColCutGenConquer(
         stages::Vector{ColumnGeneration} = [ColumnGeneration()]
-        primal_heuristics::Vector{ParameterisedHeuristic} = [RestrictedMasterHeuristic()]
+        primal_heuristics::Vector{ParameterisedHeuristic} = [ParamRestrictedMasterHeuristic()]
         preprocess = PreprocessAlgorithm()
         cutgen = CutCallbacks()
         run_preprocessing::Bool = false
@@ -139,7 +141,7 @@ several primal heuristics to more efficiently find feasible solutions.
 """
 @with_kw struct ColCutGenConquer <: AbstractConquerAlgorithm 
     stages::Vector{ColumnGeneration} = [ColumnGeneration()]
-    primal_heuristics::Vector{ParameterisedHeuristic} = [RestrictedMasterHeuristic()]
+    primal_heuristics::Vector{ParameterisedHeuristic} = [ParamRestrictedMasterHeuristic()]
     preprocess = PreprocessAlgorithm()
     cutgen = CutCallbacks()
     max_nb_cut_rounds::Int = 3 # TODO : tailing-off ?

src/Algorithm/treesearch.jl

@@ -244,7 +244,12 @@ function run_conquer_algorithm!(
 
     treestate = getoptstate(tsdata)
     nodestate = getoptstate(node)
-    update_ip_primal!(nodestate, treestate, tsdata.exploitsprimalsolutions)
+
+    update_ip_primal_bound!(nodestate, get_ip_primal_bound(treestate))
+    best_ip_primal_sol = get_best_ip_primal_sol(nodestate)
+    if tsdata.exploitsprimalsolutions && best_ip_primal_sol !== nothing
+        set_ip_primal_sol!(treestate, best_ip_primal_sol)
+    end
 
     # in the case the conquer was already run (in strong branching),
     # we still need to update the node IP primal bound before exiting 
@@ -256,10 +261,11 @@ function run_conquer_algorithm!(
         algo.opt_rtol, algo.opt_atol
     )        
 
-    update_all_ip_primal_solutions!(treestate, nodestate)
+    add_ip_primal_sols!(treestate, get_ip_primal_sols(nodestate)...)
 
-    if algo.storelpsolution && isrootnode(node) && nb_lp_primal_sols(nodestate) > 0
-        set_lp_primal_sol!(treestate, get_best_lp_primal_sol(nodestate)) 
+    best_lp_primal_sol = get_best_lp_primal_sol(nodestate)
+    if algo.storelpsolution && isrootnode(node) && best_lp_primal_sol !== nothing
+        set_lp_primal_sol!(treestate, best_lp_primal_sol) 
     end 
     return
 end
@@ -271,7 +277,7 @@ function run_divide_algorithm!(
     treestate = getoptstate(tsdata)
     output = run!(algo.dividealg, env, reform, DivideInput(node, treestate))
 
-    update_all_ip_primal_solutions!(treestate, getoptstate(output))
+    add_ip_primal_sols!(treestate, get_ip_primal_sols(getoptstate(output))...)
 
     @logmsg LogLevel(-1) string("Updating tree.")
 
@@ -348,7 +354,10 @@ function run!(
 
         remove_records!(node.recordids)
         # we delete solutions from the node optimization state, as they are not needed anymore
-        clear_solutions!(getoptstate(node))
+        nodestate = getoptstate(node)
+        empty_ip_primal_sols!(nodestate)
+        empty_lp_primal_sols!(nodestate)
+        empty_lp_dual_sols!(nodestate)
 
         if nodestatus == TIME_LIMIT
             println("Time limit is reached. Tree search is interrupted")
@@ -358,7 +367,7 @@ function run!(
     finish_branching_tree_file(algo)
 
     if treeisempty(tsdata) # it means that the BB tree has been fully explored
-        if nb_ip_primal_sols(tsdata.optstate) >= 1
+        if length(get_ip_primal_sols(tsdata.optstate)) >= 1
             if ip_gap_closed(tsdata.optstate, rtol = algo.opt_rtol, atol = algo.opt_atol)
                 setterminationstatus!(tsdata.optstate, OPTIMAL)
             else
@@ -375,7 +384,6 @@ function run!(
     while !treeisempty(tsdata)
         node = popnode!(tsdata)
         remove_records!(node.recordids)
-        clear_solutions!(node.optstate)
     end
 
     env.kpis.node_count = get_tree_order(tsdata) - 1 # TODO : check why we need to remove 1