diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 53ccc5119..ac21fd3a2 100644
--- a/.github/workflows/ci.yml
+++ b/.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
diff --git a/src/Algorithm/Algorithm.jl b/src/Algorithm/Algorithm.jl
index a147b9c98..502454334 100644
--- a/src/Algorithm/Algorithm.jl
+++ b/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
 
diff --git a/src/Algorithm/benders.jl b/src/Algorithm/benders.jl
index 46aad6178..70e992cf7 100644
--- a/src/Algorithm/benders.jl
+++ b/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
diff --git a/src/Algorithm/branching/branchingalgo.jl b/src/Algorithm/branching/branchingalgo.jl
index 243526fa0..a55b58d01 100644
--- a/src/Algorithm/branching/branchingalgo.jl
+++ b/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)
diff --git a/src/Algorithm/branching/branchinggroup.jl b/src/Algorithm/branching/branchinggroup.jl
index ec9e480ce..cd7c9ad34 100644
--- a/src/Algorithm/branching/branchinggroup.jl
+++ b/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
diff --git a/src/Algorithm/colgen.jl b/src/Algorithm/colgen.jl
index 126854d90..f01d73f2a 100644
--- a/src/Algorithm/colgen.jl
+++ b/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)
diff --git a/src/Algorithm/conquer.jl b/src/Algorithm/conquer.jl
index 37a66ea29..e11f0e52c 100644
--- a/src/Algorithm/conquer.jl
+++ b/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 ?
diff --git a/src/Algorithm/treesearch.jl b/src/Algorithm/treesearch.jl
index 4c0e6b1ab..228958226 100644
--- a/src/Algorithm/treesearch.jl
+++ b/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
diff --git a/src/Algorithm/utilities/optimizationstate.jl b/src/Algorithm/utilities/optimizationstate.jl
index 5c46310de..f6e20c58b 100644
--- a/src/Algorithm/utilities/optimizationstate.jl
+++ b/src/Algorithm/utilities/optimizationstate.jl
@@ -10,9 +10,9 @@ mutable struct OptimizationState{F<:AbstractFormulation,S<:Coluna.AbstractSense}
     insert_function_ip_primal_sols::Function
     insert_function_lp_primal_sols::Function
     insert_function_lp_dual_sols::Function
-    ip_primal_sols::Union{Nothing, Vector{PrimalSolution{F}}}
-    lp_primal_sols::Union{Nothing, Vector{PrimalSolution{F}}}
-    lp_dual_sols::Union{Nothing, Vector{DualSolution{F}}}
+    ip_primal_sols::Vector{PrimalSolution{F}}
+    lp_primal_sols::Vector{PrimalSolution{F}}
+    lp_dual_sols::Vector{DualSolution{F}}
 end
 
 _sort!(sols::Vector{PrimalSolution{F}}, f::Function) where {F} = sort!(sols, by = x -> f(PrimalBound(x.model, x.bound)))
@@ -27,7 +27,7 @@ function bestbound!(solutions::Vector{Sol}, max_len::Int, new_sol::Sol) where {S
     return
 end
 
-function set!(solutions::Vector{Sol}, max_len::Int, new_sol::Sol) where {Sol<:Solution}
+function set!(solutions::Vector{Sol}, ::Int, new_sol::Sol) where {Sol<:Solution}
     empty!(solutions)
     push!(solutions, new_sol)
     return
@@ -47,10 +47,10 @@ end
         insert_function_ip_primal_sols = bestbound!, 
         insert_function_lp_primal_sols = bestbound!, 
         insert_function_lp_dual_sols = bestbound!
-        )
+    )
 
 A convenient structure to maintain and return solutions and bounds of a formulation `form` during an
-optimization process. The termination statuses is considered as
+optimization process. The termination status is considered as
 unknown by default. You can define the initial incumbent bounds using `ip_primal_bound`,
 `ip_dual_bound`, `lp_primal_bound`, and `lp_primal_bound` keyword arguments. Incumbent
 bounds are set to infinite (according to formulation objective sense) by default.
@@ -61,6 +61,16 @@ size of the lists. Keywords `insert_function_ip_primal_sols`, `insert_function_l
 and `insert_function_lp_dual_sols` let you provide a function to define the way
 you want to insert a new solution in each list. By default, lists are sorted by
 best bound.
+
+You can also create an `OptimizationState` from another one : 
+
+    OptimizationState(
+        form, source_state, copy_ip_primal_sol, copy_lp_primal_sol
+    )
+
+It copies the termination status, all the bounds of `source_state`.
+If copies the best IP primal solution when `copy_ip_primal_sol` equals `true` and
+the best LP primal solution when `copy_lp_primal_sol` equals `true`.
 """
 function OptimizationState(
     form::F;
@@ -93,36 +103,11 @@ function OptimizationState(
         insert_function_ip_primal_sols,
         insert_function_lp_primal_sols,
         insert_function_lp_dual_sols,
-        nothing, nothing, nothing
+        PrimalSolution{F}[], PrimalSolution{F}[], DualSolution{F}[]
     )
     return state
 end
 
-# function OptimizationState(
-#     form::AbstractFormulation, or::OptimizationState
-# )
-#     newor = OptimizationState(
-#         form,
-#         termination_status = getterminationstatus(or),
-#         ip_primal_bound = get_ip_primal_bound(or),
-#         ip_dual_bound = get_ip_dual_bound(or),
-#         lp_primal_bound = get_lp_primal_bound(or),
-#         lp_dual_bound = get_lp_dual_bound(or)
-#     )
-#     if or.ip_primal_sols !== nothing
-#         newor.ip_primal_sols = copy(or.ip_primal_sols)
-#     end
-#     if or.lp_primal_sols !== nothing
-#         newor.lp_primal_sols = copy(or.lp_primal_sols)
-#     end
-#     if or.lp_primal_sols !== nothing
-#         newor.lp_dual_sols = copy(or.lp_primal_sols)
-#     end
-#     return newor
-# end
-
-getform(state::OptimizationState{F,S}) where {F, S} = F
-
 function OptimizationState(
     form::AbstractFormulation, source_state::OptimizationState, 
     copy_ip_primal_sol::Bool, copy_lp_primal_sol::Bool
@@ -135,108 +120,135 @@ function OptimizationState(
         ip_dual_bound = get_ip_dual_bound(source_state),
         lp_dual_bound = get_lp_dual_bound(source_state)
     )
-    if copy_ip_primal_sol && nb_ip_primal_sols(source_state) > 0
-        set_ip_primal_sol!(state, copy(get_best_ip_primal_sol(source_state)))
+
+    best_ip_primal_sol = get_best_ip_primal_sol(source_state)
+    if copy_ip_primal_sol && best_ip_primal_sol !== nothing
+        set_ip_primal_sol!(state, copy(best_ip_primal_sol))
     end
-    if copy_lp_primal_sol && nb_lp_primal_sols(source_state) > 0
-        set_lp_primal_sol!(state, copy(get_best_lp_primal_sol(source_state)))
+
+    best_lp_primal_sol = get_best_lp_primal_sol(source_state)
+    if copy_lp_primal_sol && best_lp_primal_sol !== nothing
+        set_lp_primal_sol!(state, copy(best_lp_primal_sol))
     end
+
     return state
 end
 
 getterminationstatus(state::OptimizationState) = state.termination_status
 setterminationstatus!(state::OptimizationState, status::TerminationStatus) = state.termination_status = status
 
-getincumbents(state::OptimizationState) = state.incumbents
 
-get_ip_primal_bound(state::OptimizationState) = get_ip_primal_bound(state.incumbents)
-get_lp_primal_bound(state::OptimizationState) = get_lp_primal_bound(state.incumbents)
-get_ip_dual_bound(state::OptimizationState) = get_ip_dual_bound(state.incumbents)
-get_lp_dual_bound(state::OptimizationState) = get_lp_dual_bound(state.incumbents)
+"Return the best IP primal bound."
+get_ip_primal_bound(state::OptimizationState) = state.incumbents.ip_primal_bound
 
-update_ip_primal_bound!(state::OptimizationState, val) = update_ip_primal_bound!(state.incumbents, val)
-update_ip_dual_bound!(state::OptimizationState, val) = update_ip_dual_bound!(state.incumbents, val)
-update_lp_primal_bound!(state::OptimizationState, val) = update_lp_primal_bound!(state.incumbents, val)
-update_lp_dual_bound!(state::OptimizationState, val) = update_lp_dual_bound!(state.incumbents, val)
+"Return the best LP primal bound."
+get_lp_primal_bound(state::OptimizationState) = state.incumbents.lp_primal_bound
 
-set_ip_primal_bound!(state::OptimizationState, val) = set_ip_primal_bound!(state.incumbents, val)
-set_lp_primal_bound!(state::OptimizationState, val) = set_lp_primal_bound!(state.incumbents, val)
-set_ip_dual_bound!(state::OptimizationState, val) = set_ip_dual_bound!(state.incumbents, val)
-set_lp_dual_bound!(state::OptimizationState, val) = set_lp_dual_bound!(state.incumbents, val)
+"Return the best IP dual bound."
+get_ip_dual_bound(state::OptimizationState) = state.incumbents.ip_dual_bound
 
-ip_gap(state::OptimizationState) = ip_gap(state.incumbents)
-lp_gap(state::OptimizationState) = lp_gap(state.incumbents)
+"Return the best LP dual bound."
+get_lp_dual_bound(state::OptimizationState) = state.incumbents.lp_dual_bound
 
-ip_gap_closed(state::OptimizationState; kw...) = ip_gap_closed(state.incumbents; kw...)
-lp_gap_closed(state::OptimizationState; kw...) = lp_gap_closed(state.incumbents; kw...)
+"""
+Update the primal bound of the mixed-integer program if the new one is better
+than the current one according to the objective sense.
+"""
+update_ip_primal_bound!(state::OptimizationState, bound) = MathProg._update_ip_primal_bound!(state.incumbents, bound)
 
-function nb_ip_primal_sols(state::OptimizationState)
-    return state.ip_primal_sols === nothing ? 0 : length(state.ip_primal_sols)
-end
+"""
+Update the dual bound of the mixed-integer program if the new one is better than
+the current one according to the objective sense.
+"""
+update_ip_dual_bound!(state::OptimizationState, bound) = MathProg._update_ip_dual_bound!(state.incumbents, bound)
 
-function nb_lp_primal_sols(state::OptimizationState)
-    return state.lp_primal_sols === nothing ? 0 : length(state.lp_primal_sols)
-end
+"""
+Update the primal bound of the linear program if the new one is better than the
+current one according to the objective sense.
+"""
+update_lp_primal_bound!(state::OptimizationState, bound) = MathProg._update_lp_primal_bound!(state.incumbents, bound)
 
-function nb_lp_dual_sols(state::OptimizationState)
-    return state.lp_dual_sols === nothing ? 0 : length(state.lp_dual_sols)
-end
+"""
+Update the dual bound of the linear program if the new one is better than the 
+current one according to the objective sense.
+"""
+update_lp_dual_bound!(state::OptimizationState, bound) = MathProg._update_lp_dual_bound!(state.incumbents, bound)
+
+"Set the best IP primal bound."
+set_ip_primal_bound!(state::OptimizationState, bound) = state.incumbents.ip_primal_bound = bound
+
+"Set the best LP primal bound."
+set_lp_primal_bound!(state::OptimizationState, bound) = state.incumbents.lp_primal_bound = bound
+
+"Set the best IP dual bound."
+set_ip_dual_bound!(state::OptimizationState, bound) = state.incumbents.ip_dual_bound = bound
+
+"Set the best LP dual bound."
+set_lp_dual_bound!(state::OptimizationState, bound) = state.incumbents.lp_dual_bound = bound
 
+"""
+Return the gap between the best primal and dual bounds of the integer program.
+Should not be used to check convergence
+"""
+ip_gap(state::OptimizationState) = MathProg._ip_gap(state.incumbents)
+
+"Return the gap between the best primal and dual bounds of the linear program."
+lp_gap(state::OptimizationState) = MathProg._lp_gap(state.incumbents)
+
+"""
+    ip_gap_closed(optstate; atol = Coluna.DEF_OPTIMALITY_ATOL, rtol = Coluna.DEF_OPTIMALITY_RTOL)
+
+Return true if the gap between the best primal and dual bounds of the integer program is closed
+given optimality tolerances.
+"""
+ip_gap_closed(state::OptimizationState; kw...) = MathProg._ip_gap_closed(state.incumbents; kw...)
+
+"""
+    lp_gap_closed(optstate; atol = Coluna.DEF_OPTIMALITY_ATOL, rtol = Coluna.DEF_OPTIMALITY_RTOL)
+
+Return true if the gap between the best primal and dual bounds of the linear program is closed 
+given optimality tolerances.
+"""
+lp_gap_closed(state::OptimizationState; kw...) = MathProg._lp_gap_closed(state.incumbents; kw...)
+
+"Return all IP primal solutions."
 get_ip_primal_sols(state::OptimizationState) = state.ip_primal_sols
 
+
+"Return the best IP primal solution if it exists; `nothing` otherwise."
 function get_best_ip_primal_sol(state::OptimizationState)
-    nb_ip_primal_sols(state) == 0 && return nothing
-    return get_ip_primal_sols(state)[1]
+    length(state.ip_primal_sols) == 0 && return nothing
+    return state.ip_primal_sols[1]
 end
 
+"Return all LP primal solutions."
 get_lp_primal_sols(state::OptimizationState) = state.lp_primal_sols
 
+"Return the best LP primal solution if it exists; `nothing` otherwise."
 function get_best_lp_primal_sol(state::OptimizationState)
-    nb_lp_primal_sols(state) == 0 && return nothing
-    return get_lp_primal_sols(state)[1]
+    length(state.lp_primal_sols) == 0 && return nothing
+    return state.lp_primal_sols[1]
 end
 
+"Return all LP dual solutions."
 get_lp_dual_sols(state::OptimizationState) = state.lp_dual_sols
 
+"Return the best LP dual solution if it exists; `nothing` otherwise."
 function get_best_lp_dual_sol(state::OptimizationState)
-    nb_lp_dual_sols(state) == 0 && return nothing
-    return get_lp_dual_sols(state)[1]
-end
-
-function clear_solutions!(state::OptimizationState)
-    state.lp_primal_sols = nothing
-    state.ip_primal_sols = nothing
-    state.lp_dual_sols = nothing
-end
-
-function update_ip_primal!(
-    dest_state::OptimizationState, orig_state::OptimizationState, set_solution::Bool
-)
-    update_ip_primal_bound!(dest_state, get_ip_primal_bound(orig_state))    
-    if set_solution && nb_ip_primal_sols(orig_state) > 0
-        set_ip_primal_sol!(dest_state, get_best_ip_primal_sol(orig_state))
-    end
-end
-
-function update_all_ip_primal_solutions!(
-    dest_state::OptimizationState, orig_state::OptimizationState
-)
-    # we do it in reverse order in order not to store the same solution several times
-    if nb_ip_primal_sols(orig_state) > 0
-        for ip_primal_sol in reverse(get_ip_primal_sols(orig_state))
-            update_ip_primal_sol!(dest_state, ip_primal_sol)
-        end    
-    end    
+    length(state.lp_dual_sols) == 0 && return nothing
+    return state.lp_dual_sols[1]
 end
 
+# TODO : refactoring ?
 function update!(dest_state::OptimizationState, orig_state::OptimizationState)
     setterminationstatus!(dest_state, getterminationstatus(orig_state))
-    update_all_ip_primal_solutions!(dest_state, orig_state)
+    add_ip_primal_sols!(dest_state, get_ip_primal_sols(orig_state)...)
     update_ip_dual_bound!(dest_state, get_ip_dual_bound(orig_state))
     update_lp_dual_bound!(dest_state, get_lp_dual_bound(orig_state))
     set_lp_primal_bound!(dest_state, get_lp_primal_bound(orig_state))
-    if nb_lp_primal_sols(orig_state) > 0
-        set_lp_primal_sol!(dest_state, get_best_lp_primal_sol(orig_state))
+    best_lp_primal_sol = get_best_lp_primal_sol(orig_state) 
+    if best_lp_primal_sol !== nothing
+        set_lp_primal_sol!(dest_state, best_lp_primal_sol)
     end        
 end
 
@@ -248,15 +260,16 @@ value of the solution is better than the incumbent. The solution is inserted in
 by the method defined in `insert_function_ip_primal_sols` field of `OptimizationState`.
 If the maximum length of the list is reached, the solution located at the end of the list
 is removed.
+
+Similar methods : 
+
+    update_lp_primal_sol!(optstate, sol)
+    update_lp_dual_sol!(optstate, sol)
 """
 function update_ip_primal_sol!(state::OptimizationState{F, S}, sol::PrimalSolution{F}) where {F, S}
     state.max_length_ip_primal_sols == 0 && return
-
-    if state.ip_primal_sols === nothing
-        state.ip_primal_sols = PrimalSolution{F}[]
-    end 
     b = PrimalBound{S}(getvalue(sol))
-    if update_ip_primal_bound!(state.incumbents, b)
+    if update_ip_primal_bound!(state, b)
         state.insert_function_ip_primal_sols(state.ip_primal_sols, state.max_length_ip_primal_sols, sol)
     end
     return
@@ -264,20 +277,29 @@ end
 
 """
     add_ip_primal_sol!(optstate, sol)
+    add_ip_primal_sols!(optstate, sols...)
 
 Add the solution `sol` at the end of the solution list of `opstate`, sort the solution list,
 remove the worst solution if the solution list size is exceded, and update the incumbent bound if 
 the solution is better.
+
+Similar methods :
+
+    add_lp_primal_sol!(optstate, sol)
+    add_lp_dual_sol!(optstate, sol)
 """
 function add_ip_primal_sol!(state::OptimizationState{F, S}, sol::PrimalSolution{F}) where {F, S}
     state.max_length_ip_primal_sols == 0 && return
-    
-    if state.ip_primal_sols === nothing
-        state.ip_primal_sols = PrimalSolution{F}[]
-    end
     state.insert_function_ip_primal_sols(state.ip_primal_sols, state.max_length_ip_primal_sols, sol)
     b = PrimalBound{S}(getvalue(sol))
-    update_ip_primal_bound!(state.incumbents, b)
+    update_ip_primal_bound!(state, b)
+    return
+end
+
+function add_ip_primal_sols!(state::OptimizationState, sols...)
+    for sol in sols
+        add_ip_primal_sol!(state, sol)
+    end
     return
 end
 
@@ -286,108 +308,107 @@ end
 
 Empties the list of solutions and add solution `sol` in the list.
 The incumbent bound is not updated even if the value of the solution is better.
+
+Similar methods :
+
+    set_lp_primal_sol!(optstate, sol)
+    set_lp_dual_sol!(optstate, sol)
 """
 function set_ip_primal_sol!(state::OptimizationState{F, S}, sol::PrimalSolution{F}) where {F, S}
     state.max_length_ip_primal_sols == 0 && return
-    
-    if state.ip_primal_sols === nothing
-        state.ip_primal_sols = PrimalSolution{F}[]
-    end
     set!(state.ip_primal_sols, state.max_length_ip_primal_sols, sol)
     return
 end
 
+"""
+    empty_ip_primal_sols!(optstate)
+
+Remove all IP primal solutions from `optstate`.
+
+Similar methods :
+
+    empty_lp_primal_sols!(optstate)
+    empty_lp_dual_sols!(optstate)
+"""
+empty_ip_primal_sols!(state::OptimizationState) = empty!(state.ip_primal_sols)
+
+"Similar to [`update_ip_primal_sol!`](@ref)."
 function update_lp_primal_sol!(state::OptimizationState{F, S}, sol::PrimalSolution{F}) where {F, S}
     state.max_length_lp_primal_sols == 0 && return
-
-    if state.lp_primal_sols === nothing
-        state.lp_primal_sols = PrimalSolution{F}[]
-    end
     b = PrimalBound{S}(getvalue(sol))
-    if update_lp_primal_bound!(state.incumbents, b)
+    if update_lp_primal_bound!(state, b)
         state.insert_function_lp_primal_sols(state.lp_primal_sols, state.max_length_lp_primal_sols, sol)
     end
     return
 end
 
+"Similar to [`add_ip_primal_sol!`](@ref)."
 function add_lp_primal_sol!(state::OptimizationState{F, S}, sol::PrimalSolution{F}) where {F, S}
     state.max_length_lp_primal_sols == 0 && return
-    
-    if state.lp_primal_sols === nothing
-        state.lp_primal_sols = PrimalSolution{F}[]
-    end
-    
     state.insert_function_lp_primal_sols(state.lp_primal_sols, state.max_length_lp_primal_sols, sol)
     b = PrimalBound{S}(getvalue(sol))
-    update_lp_primal_bound!(state.incumbents, b)
+    update_lp_primal_bound!(state, b)
     return
 end
 
+"Similar to [`set_ip_primal_sol!`](@ref)."
 function set_lp_primal_sol!(state::OptimizationState{F, S}, sol::PrimalSolution{F}) where {F, S}
     state.max_length_lp_primal_sols == 0 && return
-
-    if state.lp_primal_sols === nothing
-        state.lp_primal_sols = PrimalSolution{F}[]
-    end
     set!(state.lp_primal_sols, state.max_length_lp_primal_sols, sol)
     return
 end
 
+"Similar to [`empty_ip_primal_sols!`](@ref)."
+empty_lp_primal_sols!(state::OptimizationState) = empty!(state.lp_primal_sols)
+
+"Similar to [`update_ip_primal_sol!`](@ref)."
 function update_lp_dual_sol!(state::OptimizationState{F, S}, sol::DualSolution{F}) where {F, S}
     state.max_length_lp_dual_sols == 0 && return
-
-    if state.lp_dual_sols === nothing
-        state.lp_dual_sols = DualSolution{F}[]
-    end
     b = DualBound{S}(getvalue(sol))
-    if update_lp_dual_bound!(state.incumbents, b)
+    if update_lp_dual_bound!(state, b)
         state.insert_function_lp_dual_sols(state.lp_dual_sols, state.max_length_lp_dual_sols, sol)
     end
     return
 end
 
+"Similar to [`add_ip_primal_sol!`](@ref)."
 function add_lp_dual_sol!(state::OptimizationState{F, S}, sol::DualSolution{F}) where {F, S}
     state.max_length_lp_dual_sols == 0 && return
-
-    if state.lp_dual_sols === nothing
-        state.lp_dual_sols = DualSolution{F}[]
-    end
     state.insert_function_lp_dual_sols(state.lp_dual_sols, state.max_length_lp_dual_sols, sol)
     b = DualBound{S}(getvalue(sol))
-    update_lp_dual_bound!(state.incumbents, b)
+    update_lp_dual_bound!(state, b)
     return
 end
 
+"Similar to [`set_ip_primal_sol!`](@ref)."
 function set_lp_dual_sol!(state::OptimizationState{F, S}, sol::DualSolution{F}) where {F, S}
     state.max_length_lp_dual_sols == 0 && return
-
-    if state.lp_dual_sols === nothing
-        state.lp_dual_sols = DualSolution{F}[]
-    end
     set!(state.lp_dual_sols, state.max_length_lp_dual_sols, sol)
     return
 end
 
+"Similar to [`empty_ip_primal_sols!`](@ref)."
+empty_lp_dual_sols!(state::OptimizationState) = empty!(state.lp_dual_sols)
 
 function Base.print(io::IO, form::AbstractFormulation, optstate::OptimizationState)
     println(io, "┌ Optimization state ")
     println(io, "│ Termination status: ", optstate.termination_status)
     println(io, "| Incumbents: ", optstate.incumbents)
-    n = nb_ip_primal_sols(optstate)
+    n = length(optstate.ip_primal_sols)
     println(io, "| IP Primal solutions (",n,")")
     if n > 0
         for sol in optstate.ip_primal_sols
             print(io, form, sol)
         end
     end
-    n = nb_lp_primal_sols(optstate)
+    n = length(optstate.lp_primal_sols)
     println(io, "| LP Primal solutions (",n,"):")
     if n > 0
         for sol in optstate.lp_primal_sols
             print(io, form, sol)
         end
     end
-    n = nb_lp_dual_sols(optstate)
+    n = length(optstate.lp_dual_sols)
     println(io, "| LP Dual solutions (",n,"):")
     if n > 0
         for sol in optstate.lp_dual_sols
diff --git a/src/ColunaBase/solsandbounds.jl b/src/ColunaBase/solsandbounds.jl
index c88ed074a..7a18394e9 100644
--- a/src/ColunaBase/solsandbounds.jl
+++ b/src/ColunaBase/solsandbounds.jl
@@ -40,7 +40,8 @@ isbetter(b1::Bound{Sp,Se}, b2::Bound{Sp,Se}) where {Sp<:Dual,Se<:MinSense} = b1.
 isbetter(b1::Bound{Sp,Se}, b2::Bound{Sp,Se}) where {Sp<:Dual,Se<:MaxSense} = b1.value < b2.value
 
 """
-    diff(b1, b2)
+    diff(pb, db)
+    diff(db, pb)
 
 Distance between a primal bound and a dual bound that have the same objective sense.
 Distance is non-positive if dual bound reached primal bound.
@@ -62,9 +63,10 @@ function Base.diff(db::Bound{<:Dual,<:MaxSense}, pb::Bound{<:Primal,<:MaxSense})
 end
 
 """
-    gap
+    gap(pb, db)
+    gap(db, pb)
 
-relative gap. Gap is non-positive if pb reached db
+Return relative gap. Gap is non-positive if pb reached db.
 """
 function gap(pb::Bound{<:Primal,<:MinSense}, db::Bound{<:Dual,<:MinSense})
     return diff(pb, db) / abs(db.value)
diff --git a/src/MOIwrapper.jl b/src/MOIwrapper.jl
index 75de5f631..2b2025526 100644
--- a/src/MOIwrapper.jl
+++ b/src/MOIwrapper.jl
@@ -683,7 +683,7 @@ function MOI.get(optimizer::Optimizer, ::MOI.RawStatusString)
 end
 
 function MOI.get(optimizer::Optimizer, ::MOI.ResultCount)
-    return nb_ip_primal_sols(optimizer.result)
+    return length(get_ip_primal_sols(optimizer.result))
 end
 
 function MOI.get(
diff --git a/src/MathProg/MathProg.jl b/src/MathProg/MathProg.jl
index f53ef77da..5c08bacca 100644
--- a/src/MathProg/MathProg.jl
+++ b/src/MathProg/MathProg.jl
@@ -100,12 +100,7 @@ export Variable, Constraint, VarId, ConstrId, VarMembership, ConstrMembership,
     isexplicit, getname, getbranchingpriority, reset!, getreducedcost
 
 # Types & methods related to solutions & bounds
-export PrimalBound, DualBound, PrimalSolution, DualSolution, ObjValues,
-    get_ip_primal_bound, get_lp_primal_bound,
-    get_ip_dual_bound, get_lp_dual_bound, update_ip_primal_bound!, update_lp_primal_bound!,
-    update_ip_dual_bound!, update_lp_dual_bound!, set_ip_primal_bound!,
-    set_lp_primal_bound!, set_ip_dual_bound!, set_lp_dual_bound!, ip_gap, lp_gap, ip_gap_closed, 
-    lp_gap_closed
+export PrimalBound, DualBound, PrimalSolution, DualSolution, ObjValues
 
 # Methods related to projections
 export projection_is_possible, proj_cols_on_rep
diff --git a/src/MathProg/bounds.jl b/src/MathProg/bounds.jl
index e772e28a8..6c7a75997 100644
--- a/src/MathProg/bounds.jl
+++ b/src/MathProg/bounds.jl
@@ -1,12 +1,20 @@
 const PrimalBound{S} = Bound{Primal, S}
 const DualBound{S} = Bound{Dual, S}
 
+"""
+    PrimalBound(formulation)
+    PrimalBound(formulation, value)
+    PrimalBound(formualtion, pb)
+
+Create a new primal bound for the formulation `formulation`.
+The value of the primal bound is infinity if you do not specify any initial value.
+"""
 function PrimalBound(form::AbstractFormulation)
     Se = getobjsense(form)
     return Bound{Primal,Se}()
 end
 
-function PrimalBound(form::AbstractFormulation, val::Float64)
+function PrimalBound(form::AbstractFormulation, val)
     Se = getobjsense(form)
     return Bound{Primal,Se}(val)
 end
@@ -23,12 +31,20 @@ function PrimalBound(form::AbstractFormulation, pb::PrimalBound{S}) where {S}
     return Bound{Primal,Se}(getvalue(pb))
 end
 
+"""
+    DualBound(formulation)
+    DualBound(formulation, value)
+    DualBound(formulation, db)
+
+Create a new dual bound for the formulation `formulation`.
+The value of the dual bound is infinity if you do not specify any initial value.
+"""
 function DualBound(form::AbstractFormulation)
     Se = getobjsense(form)
     return Bound{Dual,Se}()
 end
 
-function DualBound(form::AbstractFormulation, val::Float64)
+function DualBound(form::AbstractFormulation, val)
     Se = getobjsense(form)
     return Bound{Dual,Se}(val)
 end
@@ -58,11 +74,7 @@ mutable struct ObjValues{S}
     ip_dual_bound::DualBound{S}
 end
 
-"""
-    ObjValues(form)
-
-todo
-"""
+"A convenient structure to maintain and return incumbent bounds."
 function ObjValues(
     form::M; 
     ip_primal_bound = nothing,
@@ -75,58 +87,37 @@ function ObjValues(
         PrimalBound(form), DualBound(form), PrimalBound(form), DualBound(form)
     )
     if ip_primal_bound !== nothing
-        set_ip_primal_bound!(ov, PrimalBound(form, ip_primal_bound))
+        ov.ip_primal_bound = PrimalBound(form, ip_primal_bound)
     end
     if ip_dual_bound !== nothing
-        set_ip_dual_bound!(ov, DualBound(ip_dual_bound))
+        ov.ip_dual_bound = DualBound(form, ip_dual_bound)
     end
     if lp_primal_bound !== nothing
-        set_lp_primal_bound!(ov, PrimalBound(lp_primal_bound))
+        ov.lp_primal_bound = PrimalBound(form, lp_primal_bound)
     end
     if lp_dual_bound !== nothing
-        set_lp_dual_bound!(ov, DualBound(lp_dual_bound))
+        ov.lp_dual_bound = DualBound(form, lp_dual_bound)
     end
     return ov
 end
 
-## Getters bounds
-"Return the best primal bound of the mixed-integer program."
-get_ip_primal_bound(ov::ObjValues) = ov.ip_primal_bound
-
-"Return the best dual bound of the mixed-integer program."
-get_ip_dual_bound(ov::ObjValues) = ov.ip_dual_bound
-
-"Return the best primal bound of the linear program."
-get_lp_primal_bound(ov::ObjValues) = ov.lp_primal_bound
-
-"Return the best dual bound of the linear program."
-get_lp_dual_bound(ov::ObjValues) = ov.lp_dual_bound
-
 ## Gaps
-"""
-Return the gap between the best primal and dual bounds of the integer program.
-Should not be used to check convergence
-"""
-ip_gap(ov::ObjValues) = gap(get_ip_primal_bound(ov), get_ip_dual_bound(ov))
-
-"Return the gap between the best primal and dual bounds of the linear program."
-lp_gap(ov::ObjValues) = gap(get_lp_primal_bound(ov), get_lp_dual_bound(ov))
+_ip_gap(ov::ObjValues) = gap(ov.ip_primal_bound, ov.ip_dual_bound)
+_lp_gap(ov::ObjValues) = gap(ov.lp_primal_bound, ov.lp_dual_bound)
 
-"Return true if the gap between the best primal and dual bounds of the integer program is closed given optimality tolerances."
-function ip_gap_closed(
+function _ip_gap_closed(
     ov::ObjValues; atol = Coluna.DEF_OPTIMALITY_ATOL, rtol = Coluna.DEF_OPTIMALITY_RTOL
 )
-    return (ip_gap(ov) <= 0) || _gap_closed(
-        get_ip_primal_bound(ov).value, get_ip_dual_bound(ov).value, atol = atol, rtol = rtol
+    return (_ip_gap(ov) <= 0) || _gap_closed(
+        ov.ip_primal_bound.value, ov.ip_dual_bound.value, atol = atol, rtol = rtol
     )
 end
 
-"Return true if the gap between the best primal and dual bounds of the linear program is closed given optimality tolerances."
-function lp_gap_closed(
+function _lp_gap_closed(
     ov::ObjValues; atol = Coluna.DEF_OPTIMALITY_ATOL, rtol = Coluna.DEF_OPTIMALITY_RTOL
 )
-    return (lp_gap(ov) <= 0) || _gap_closed(
-        get_lp_primal_bound(ov).value, get_lp_dual_bound(ov).value, atol = atol, rtol = rtol
+    return (_lp_gap(ov) <= 0) || _gap_closed(
+       ov.lp_primal_bound.value, ov.lp_dual_bound.value, atol = atol, rtol = rtol
     )
 end
 
@@ -137,86 +128,35 @@ function _gap_closed(
     return (x == y) || (isfinite(x) && isfinite(y) && norm(x - y) <= max(atol, rtol*min(norm(x), norm(y))))
 end
 
-function set_lp_primal_bound!(ov::ObjValues{S}, b::PrimalBound{S}) where {S}
-    ov.lp_primal_bound = b
-    return
-end
-
-function set_lp_dual_bound!(ov::ObjValues{S}, b::DualBound{S}) where {S}
-    ov.lp_dual_bound = b
-    return
-end
-
-function set_ip_primal_bound!(ov::ObjValues{S}, b::PrimalBound{S}) where {S}
-    ov.ip_primal_bound = b
-    return
-end
-
-function set_ip_dual_bound!(ov::ObjValues{S}, b::DualBound{S}) where {S}
-    ov.ip_dual_bound = b
-    return
-end
-
-"""
-Update the primal bound of the linear program if the new one is better than the
-current one according to the objective sense.
-"""
-function update_lp_primal_bound!(ov::ObjValues{S}, b::PrimalBound{S}) where {S}
-    if isbetter(b, get_lp_primal_bound(ov))
+## Bound updates
+function _update_lp_primal_bound!(ov::ObjValues{S}, b::PrimalBound{S}) where {S}
+    if isbetter(b, ov.lp_primal_bound)
         ov.lp_primal_bound = b
         return true
     end
     return false
 end
 
-"""
-Update the dual bound of the linear program if the new one is better than the 
-current one according to the objective sense.
-"""
-function update_lp_dual_bound!(ov::ObjValues{S}, b::DualBound{S}) where {S}
-    if isbetter(b, get_lp_dual_bound(ov))
+function _update_lp_dual_bound!(ov::ObjValues{S}, b::DualBound{S}) where {S}
+    if isbetter(b, ov.lp_dual_bound)
         ov.lp_dual_bound = b
         return true
     end
     return false
 end
 
-"""
-Update the primal bound of the mixed-integer program if the new one is better
-than the current one according to the objective sense.
-"""
-function update_ip_primal_bound!(ov::ObjValues{S}, b::PrimalBound{S}) where {S}
-    if isbetter(b, get_ip_primal_bound(ov))
+function _update_ip_primal_bound!(ov::ObjValues{S}, b::PrimalBound{S}) where {S}
+    if isbetter(b, ov.ip_primal_bound)
         ov.ip_primal_bound = b
         return true
     end
     return false
 end
 
-"""
-Update the dual bound of the mixed-integer program if the new one is better than
-the current one according to the objective sense.
-"""
-function update_ip_dual_bound!(ov::ObjValues{S}, b::DualBound{S}) where {S}
-    if isbetter(b, get_ip_dual_bound(ov))
+function _update_ip_dual_bound!(ov::ObjValues{S}, b::DualBound{S}) where {S}
+    if isbetter(b, ov.ip_dual_bound)
         ov.ip_dual_bound = b
         return true
     end
     return false
-end
-
-function update!(dest::ObjValues{S}, src::ObjValues{S}) where {S}
-    update_ip_dual_bound!(dest, get_ip_dual_bound(src))
-    update_ip_primal_bound!(dest, get_ip_primal_bound(src))
-    update_lp_dual_bound!(dest, get_lp_dual_bound(src))
-    update_lp_primal_bound!(dest, get_lp_primal_bound(src))
-    return
-end
-
-function Base.show(io::IO, ov::ObjValues{S}) where {S}
-    println(io, "ObjValues{", S, "}:")
-    println(io, "ip_primal_bound : ", ov.ip_primal_bound)
-    println(io, "ip_dual_bound : ", ov.ip_dual_bound)
-    println(io, "lp_primal_bound : ", ov.lp_primal_bound)
-    println(io, "lp_dual_bound : ", ov.lp_dual_bound)
 end
\ No newline at end of file
diff --git a/test/MathOptInterface/MOI_wrapper.jl b/test/MathOptInterface/MOI_wrapper.jl
index 2160f5fce..9ae5fec59 100644
--- a/test/MathOptInterface/MOI_wrapper.jl
+++ b/test/MathOptInterface/MOI_wrapper.jl
@@ -48,16 +48,16 @@ end
     model = BlockModel(coluna, direct_model=true)
     @variable(model, x)
     
-    @test BlockDecomposition.branchingpriority(model, x) == 1
-    BlockDecomposition.branchingpriority!(model, x, 2)
-    @test BlockDecomposition.branchingpriority(model, x) == 2
+    @test BlockDecomposition.branchingpriority(x) == 1
+    BlockDecomposition.branchingpriority!(x, 2)
+    @test BlockDecomposition.branchingpriority(x) == 2
     
     model2 = BlockModel(coluna)
     @variable(model2, x)
     
-    @test BlockDecomposition.branchingpriority(model2, x) == 1
-    BlockDecomposition.branchingpriority!(model2, x, 2)
-    @test BlockDecomposition.branchingpriority(model2, x) == 2
+    @test BlockDecomposition.branchingpriority(x) == 1
+    BlockDecomposition.branchingpriority!(x, 2)
+    @test BlockDecomposition.branchingpriority(x) == 2
 end
 
 @testset "write_to_file" begin
diff --git a/test/full_instances_tests.jl b/test/full_instances_tests.jl
index 6bbf48175..61712c72d 100644
--- a/test/full_instances_tests.jl
+++ b/test/full_instances_tests.jl
@@ -84,7 +84,7 @@ function generalized_assignment_tests()
         # we increase the branching priority of variables which assign jobs to the first two machines
         for machine in 1:2
             for job in data.jobs
-                BD.branchingpriority!(model, x[machine,job], 2)
+                BD.branchingpriority!(x[machine,job], 2)
             end
         end  
 
diff --git a/test/runtests.jl b/test/runtests.jl
index 7ac7091cf..effa97f17 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,3 +1,4 @@
+using Base.CoreLogging: Error
 using Coluna
 
 using Test, GLPK, ColunaDemos, JuMP, BlockDecomposition
diff --git a/test/unit/MathProg/bounds.jl b/test/unit/MathProg/bounds.jl
new file mode 100644
index 000000000..cb469d1ad
--- /dev/null
+++ b/test/unit/MathProg/bounds.jl
@@ -0,0 +1,68 @@
+function mathprog_bounds()
+    env = Coluna.Env(Coluna.Params())
+
+    min_form = Coluna.MathProg.create_formulation!(
+        env, Coluna.MathProg.Original();
+        obj_sense = Coluna.MathProg.MinSense
+    )
+
+    max_form = Coluna.MathProg.create_formulation!(
+        env, Coluna.MathProg.Original();
+        obj_sense = Coluna.MathProg.MaxSense
+    )
+
+    @testset "Primal bound constructor" begin
+        pb1 = Coluna.MathProg.PrimalBound(min_form)
+        @test pb1 == Inf
+        pb2 = Coluna.MathProg.PrimalBound(max_form)
+        @test pb2 == -Inf
+        pb3 = Coluna.MathProg.PrimalBound(min_form, 10)
+        @test pb3 == 10
+        @test_throws ErrorException Coluna.MathProg.PrimalBound(max_form, pb3)
+    end
+
+    @testset "Dual bound constructor" begin
+        db1 = Coluna.MathProg.DualBound(min_form)
+        @test db1 == -Inf
+        db2 = Coluna.MathProg.DualBound(max_form)
+        @test db2 == Inf
+        db3 = Coluna.MathProg.DualBound(min_form, 150)
+        @test db3 == 150
+        @test_throws ErrorException Coluna.MathProg.DualBound(max_form, db3)
+    end
+
+    @testset "ObjValues constructor & private methods" begin
+        obj = ObjValues(
+            min_form;
+            ip_primal_bound = 15.0,
+            ip_dual_bound = 12.0,
+            lp_primal_bound = 66,
+            lp_dual_bound = π
+        )
+
+        @test obj.ip_primal_bound == 15.0
+        @test obj.ip_dual_bound == 12.0
+        @test obj.lp_primal_bound == 66
+        @test obj.lp_dual_bound == float(π) # precision...
+
+        # Gap methods are already tested in containers/solsandbounds.jl
+
+        @test Coluna.MathProg._update_ip_primal_bound!(obj, PrimalBound(min_form, 16.0)) == false
+        @test Coluna.MathProg._update_ip_primal_bound!(obj, PrimalBound(min_form, 14.0)) == true
+
+        @test Coluna.MathProg._update_lp_primal_bound!(obj, PrimalBound(min_form, 67.0)) == false
+        @test Coluna.MathProg._update_lp_primal_bound!(obj, PrimalBound(min_form, 65.0)) == true
+
+        @test Coluna.MathProg._update_ip_dual_bound!(obj, DualBound(min_form, 11.0)) == false
+        @test Coluna.MathProg._update_ip_dual_bound!(obj, DualBound(min_form, 13.0)) == true
+
+        @test Coluna.MathProg._update_lp_dual_bound!(obj, DualBound(min_form, 3.0)) == false
+        @test Coluna.MathProg._update_lp_dual_bound!(obj, DualBound(min_form, 3.2)) == true
+
+        @test obj.ip_primal_bound == 14.0
+        @test obj.ip_dual_bound == 13.0
+        @test obj.lp_primal_bound == 65
+        @test obj.lp_dual_bound == 3.2
+    end
+    return
+end
\ No newline at end of file
diff --git a/test/unit/containers/solsandbounds.jl b/test/unit/containers/solsandbounds.jl
index 8da6e7153..406ec035d 100644
--- a/test/unit/containers/solsandbounds.jl
+++ b/test/unit/containers/solsandbounds.jl
@@ -5,9 +5,12 @@ MaxSense = Coluna.AbstractMaxSense
 
 CB = Coluna.ColunaBase
 
-function bound_unit()
+struct FakeModel <: CB.AbstractModel end
 
+function bound_unit()
     @testset "Bound" begin
+        # Make sure that Coluna initializes bounds to infinity.
+        # Check that initial value of the bound is correct.
         pb = CB.Bound{Primal,MinSense}()
         @test pb == Inf
         @test CB.getvalue(pb) == Inf
@@ -245,8 +248,6 @@ function test_solution_iterations(solution::CB.Solution, dict::Dict)
     return
 end
 
-struct FakeModel <: CB.AbstractModel end
-
 function solution_unit()
     @testset "MOI Termination Status" begin
         @test CB.convert_status(MOI.OPTIMAL) == CB.OPTIMAL
@@ -279,30 +280,32 @@ function solution_unit()
         @test CB.convert_status(CB.UNCOVERED_SOLUTION_STATUS) == MOI.OTHER_RESULT_STATUS
     end
 
-    model = FakeModel()
+    @testset "Solution" begin
+        model = FakeModel()
 
-    Solution = CB.Solution{FakeModel,Int,Float64}
+        Solution = CB.Solution{FakeModel,Int,Float64}
 
-    dict_sol, soldecs, solvals = fake_solution_factory(100)
-    primal_sol = Solution(model, soldecs, solvals, 12.3, CB.FEASIBLE_SOL)
-    test_solution_iterations(primal_sol, dict_sol)
-    @test CB.getvalue(primal_sol) == 12.3
-    @test CB.getstatus(primal_sol) == CB.FEASIBLE_SOL
-    
-    dict_sol = Dict(1 => 2.0, 2 => 3.0, 3 => 4.0)
-    primal_sol = Solution(model, collect(keys(dict_sol)), collect(values(dict_sol)), 0.0, Coluna.ColunaBase.FEASIBLE_SOL)
-    
-    @test iterate(primal_sol) == iterate(primal_sol.sol)
-    _, state = iterate(primal_sol)
-    @test iterate(primal_sol, state) == iterate(primal_sol.sol, state)
-    @test length(primal_sol) == 3
-    @test primal_sol[1] == 2.0
-    primal_sol[1] = 5.0 # change the value
-    @test primal_sol[1] == 5.0
-    
-    io = IOBuffer()
-    show(io, primal_sol)
+        dict_sol, soldecs, solvals = fake_solution_factory(100)
+        primal_sol = Solution(model, soldecs, solvals, 12.3, CB.FEASIBLE_SOL)
+        test_solution_iterations(primal_sol, dict_sol)
+        @test CB.getvalue(primal_sol) == 12.3
+        @test CB.getstatus(primal_sol) == CB.FEASIBLE_SOL
+        
+        dict_sol = Dict(1 => 2.0, 2 => 3.0, 3 => 4.0)
+        primal_sol = Solution(model, collect(keys(dict_sol)), collect(values(dict_sol)), 0.0, Coluna.ColunaBase.FEASIBLE_SOL)
+        
+        @test iterate(primal_sol) == iterate(primal_sol.sol)
+        _, state = iterate(primal_sol)
+        @test iterate(primal_sol, state) == iterate(primal_sol.sol, state)
+        @test length(primal_sol) == 3
+        @test primal_sol[1] == 2.0
+        primal_sol[1] = 5.0 # change the value
+        @test primal_sol[1] == 5.0
+        
+        io = IOBuffer()
+        show(io, primal_sol)
 
-    @test String(take!(io)) == "Solution\n| 1 = 5.0\n| 2 = 3.0\n| 3 = 4.0\n└ value = 0.00 \n"
+        @test String(take!(io)) == "Solution\n| 1 = 5.0\n| 2 = 3.0\n| 3 = 4.0\n└ value = 0.00 \n"
+    end
     return
 end
\ No newline at end of file
diff --git a/test/unit/unit_tests.jl b/test/unit/unit_tests.jl
index a257f7592..bef022aec 100644
--- a/test/unit/unit_tests.jl
+++ b/test/unit/unit_tests.jl
@@ -5,6 +5,7 @@ include("MathProg/buffer.jl")
 include("MathProg/formulations.jl")
 include("MathProg/types.jl")
 include("MathProg/variables.jl")
+include("MathProg/bounds.jl")
 
 include("variable.jl")
 include("constraint.jl")
@@ -21,11 +22,13 @@ function unit_tests()
         max_nb_form_unit()
         types_unit_tests()
         variables_unit_tests()
+        mathprog_bounds()
     end
 
     @testset "variable.jl" begin
         variable_unit_tests()
     end
+    
     @testset "constraint.jl" begin
         constraint_unit_tests()
     end
diff --git a/test/user_algorithms_tests.jl b/test/user_algorithms_tests.jl
index aad4221b7..dfe6e8187 100644
--- a/test/user_algorithms_tests.jl
+++ b/test/user_algorithms_tests.jl
@@ -9,7 +9,7 @@ end
 @with_kw struct ConsecutiveColGen <: AbstractOptimizationAlgorithm
     colgen = ColumnGeneration(smoothing_stabilization = 1.0)
     preprocess = PreprocessAlgorithm(preprocess_subproblems = true, printing = true)
-    rm_heur = RestrictedMasterHeuristic()
+    rm_heur = RestrictedMasterIPHeuristic()
     num_calls_to_col_gen = 3  
 end
 
@@ -39,7 +39,7 @@ function Coluna.Algorithm.run!(
 
         cg_output = run!(algo.colgen, env, reform, OptimizationInput(optstate))
         cg_optstate = getoptstate(cg_output)
-        update_all_ip_primal_solutions!(optstate, cg_optstate)
+        add_ip_primal_sols!(optstate, get_ip_primal_sols(cg_optstate)...)
 
         primal_sol = get_best_lp_primal_sol(cg_optstate)
 
@@ -65,7 +65,7 @@ function Coluna.Algorithm.run!(
     end
 
     heur_output = run!(algo.rm_heur, env, reform, OptimizationInput(optstate))
-    update_all_ip_primal_solutions!(optstate, getoptstate(heur_output))
+    add_ip_primal_sols!(optstate, get_ip_primal_sols(getoptstate(heur_output))...)
 
     return OptimizationOutput(optstate)
 end