Dynamic sparse vector for solutions (#247)

* draft

* tol

* fix tolerance in fill_primal_sol

* rename tol params

src/Algorithm/benders.jl

@@ -611,3 +611,4 @@ function print_intermediate_statistics(algdata::BendersCutGenData,
             nb_bc_iterations, Coluna._elapsed_solve_time(), mst_time, sp_time, nb_new_cut, mlp, db, pb
     )
 end
+

src/Containers/Containers.jl

@@ -4,7 +4,7 @@ using DynamicSparseArrays
 
 import ..Coluna
 
-import DataStructures
+import DynamicSparseArrays
 import Primes
 import Printf
 import Base: <=, setindex!, get, getindex, haskey, keys, values, iterate, 

src/Containers/solsandbounds.jl

@@ -123,7 +123,7 @@ Base.:>(b1::B, b2::B) where {B<:Bound} = b1.value > b2.value
 # Solution
 struct Solution{Space<:Coluna.AbstractSpace,Sense<:Coluna.AbstractSense,Decision,Value} <: AbstractDict{Decision,Value}
     bound::Bound{Space,Sense}
-    sol::DataStructures.SortedDict{Decision,Value}
+    sol::DynamicSparseArrays.PackedMemoryArray{Decision,Value}
 end
 
 """
@@ -134,18 +134,18 @@ doc todo
 """
 function Solution{Sp,Se,De,Va}() where {Sp<:Coluna.AbstractSpace,Se<:Coluna.AbstractSense,De,Va}
     bound = Bound{Sp,Se}()
-    sol = DataStructures.SortedDict{De,Va}()
+    sol = DynamicSparseArrays.dynamicsparsevec(De[], Va[])
     return Solution(bound, sol)
 end
 
-function Solution{Sp,Se,De,Va}(solution::Dict{De,Va}, value::Float64) where {Sp<:Coluna.AbstractSpace,Se<:Coluna.AbstractSense,De,Va}
+function Solution{Sp,Se,De,Va}(decisions::Vector{De}, vals::Vector{Va}, value::Float64) where {Sp<:Coluna.AbstractSpace,Se<:Coluna.AbstractSense,De,Va}
     bound = Bound{Sp,Se}(value)
-    sol = DataStructures.SortedDict{De,Va}(solution)
+    sol = DynamicSparseArrays.dynamicsparsevec(decisions, vals)
     return Solution(bound, sol)
 end
 
-function Solution{Sp,Se,De,Va}(solution::Dict{De,Va}, bound::Bound{Sp,Se}) where {Sp<:Coluna.AbstractSpace,Se<:Coluna.AbstractSense,De,Va}
-    sol = DataStructures.SortedDict{De,Va}(solution)
+function Solution{Sp,Se,De,Va}(decisions::Vector{De}, vals::Vector{Va}, bound::Bound{Sp,Se}) where {Sp<:Coluna.AbstractSpace,Se<:Coluna.AbstractSense,De,Va}
+    sol = DynamicSparseArrays.dynamicsparsevec(decisions, vals)
     return Solution(bound, sol)
 end
 
@@ -156,8 +156,21 @@ setvalue!(s::Solution, v::Float64) = s.bound.value = v
 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{Sp,Se,De,Va}, id::De, default) where {Sp,Se,De,Va} = Base.get(s.sol, id, default)
-Base.setindex!(s::Solution{Sp,Se,De,Va}, val::Va, id::De) where {Sp,Se,De,Va} = Base.setindex!(s.sol, val, id)
+Base.get(s::Solution{Sp,Se,De,Va}, id::De, default) where {Sp,Se,De,Va} = s.sol[id]
+Base.setindex!(s::Solution{Sp,Se,De,Va}, val::Va, id::De) where {Sp,Se,De,Va} = s.sol[id] = val
+
+# todo: move in DynamicSparseArrays or avoid using filter ?
+function Base.filter(f::Function, pma::DynamicSparseArrays.PackedMemoryArray{K,T,P}) where {K,T,P}
+    ids = Vector{K}()
+    vals = Vector{T}()
+    for e in pma
+        if f(e)
+            push!(ids, e[1])
+            push!(vals, e[2])
+        end
+    end
+    return DynamicSparseArrays.dynamicsparsevec(ids, vals)
+end
 
 function Base.filter(f::Function, s::S) where {S <: Solution}
     return S(s.bound, filter(f, s.sol))

src/MathProg/MOIinterface.jl

@@ -183,25 +183,48 @@ function remove_from_optimizer!(optimizer::MoiOptimizer,
     return
 end
 
-function fill_primal_result!(optimizer::MoiOptimizer, 
+function _getcolunakind(record::MoiVarRecord)
+    record.kind.value == -1 && return Continuous
+    record.kind isa MoiBinary && return Binary
+    return Integ
+end
+
+function fill_primal_result!(form::Formulation, optimizer::MoiOptimizer, 
                              result::OptimizationResult{S},
-                             vars::VarDict) where {S<:Coluna.AbstractSense}
+                             ) where {S<:Coluna.AbstractSense}
     inner = getinner(optimizer)
     for res_idx in 1:MOI.get(inner, MOI.ResultCount())
         if MOI.get(inner, MOI.PrimalStatus(res_idx)) != MOI.FEASIBLE_POINT
             continue
         end
         pb = PrimalBound{S}(MOI.get(inner, MOI.ObjectiveValue()))
-        sol = Dict{Id{Variable}, Float64}()
-        for (id, var) in vars
-            moi_index = getindex(getmoirecord(var))
+        solvars = Vector{VarId}()
+        solvals = Vector{Float64}()
+        for (id, var) in getvars(form)
+            getcurisactive(form, id) && getcurisexplicit(form, id) || continue
+            moirec = getmoirecord(var)
+            moi_index = getindex(moirec)
+            kind = _getcolunakind(moirec)
             val = MOI.get(inner, MOI.VariablePrimal(res_idx), moi_index)
-            if val > 0.000001  || val < - 0.000001 # todo use a tolerance
+            if kind == Integ || kind == Binary
+                val = round(val, digits = Coluna._params_.tol_digits)
+            else # Continuous
+                moi_bounds = getbounds(moirec)
+                moiset = MOI.get(inner, MOI.ConstraintSet(), getbounds(moirec))
+                if val < moiset.lower
+                    val = moiset.lower
+                end
+                if val > moiset.upper
+                    val = moiset.upper
+                end
+            end
+            if abs(val) > Coluna._params_.tol
                 @logmsg LogLevel(-4) string("Var ", var.name , " = ", val)
-                sol[id] = val
+                push!(solvars, id)
+                push!(solvals, val)
             end
         end
-        push!(result.primal_sols, PrimalSolution{S}(sol, pb))
+        push!(result.primal_sols, PrimalSolution{S}(solvars, solvals, pb))
     end
     result.primal_bound = PrimalBound{S}()
     if nbprimalsols(result) > 0
@@ -220,18 +243,20 @@ function fill_dual_result!(optimizer::MoiOptimizer,
             continue
         end
         db = DualBound{S}(MOI.get(inner, MOI.ObjectiveValue()))
-        sol = Dict{Id{Constraint}, Float64}()
+        solconstrs = Vector{ConstrId}()
+        solvals = Vector{Float64}()
         # Getting dual bound is not stable in some solvers. 
         # Getting primal bound instead, which will work for lps
         for (id, constr) in constrs
             moi_index = getindex(getmoirecord(constr))
             val = MOI.get(inner, MOI.ConstraintDual(res_idx), moi_index)
             if val > 0.000001 || val < - 0.000001 # todo use a tolerance
                 @logmsg LogLevel(-4) string("Constr ", constr.name, " = ", val)
-                sol[id] = val           
+                push!(solconstrs, id)
+                push!(solvals, val)      
             end
         end
-        push!(result.dual_sols, DualSolution{S}(sol, db))
+        push!(result.dual_sols, DualSolution{S}(solconstrs, solvals, db))
     end
     result.dual_bound = DualBound{S}()
     if nbdualsols(result) > 0

src/MathProg/formulation.jl

@@ -193,7 +193,7 @@ function setprimalsol!(
 
         is_identical = true
         for (var_id, var_val) in getrecords(sol)
-            if !haskey(newprimalsol.sol, var_id)
+            if newprimalsol[var_id] == 0
                 is_identical = false
                 break
             end
@@ -203,7 +203,7 @@ function setprimalsol!(
         end
 
         for (var_id, var_val) in getrecords(newprimalsol.sol)
-            if !haskey(sol, var_id)
+            if sol[var_id] == 0
                 is_identical = false
                 break
             end
@@ -243,7 +243,6 @@ function adddualsol!(
     return dualsol_id
 end
 
-
 function setdualsol!(
     form::Formulation,
     new_dual_sol::DualSolution{S}
@@ -266,7 +265,7 @@ function setdualsol!(
 
         is_identical = true
         for (constr_id, constr_val) in getrecords(prev_dual_sol)
-            if !haskey(new_dual_sol.sol, constr_id)
+            if new_dual_sol[constr_id] == 0
                 is_identical = false
                 break
             end
@@ -276,7 +275,7 @@ function setdualsol!(
         end
 
         for (constr_id, constr_val) in new_dual_sol
-            if !haskey(prev_dual_sol, constr_id)
+            if prev_dual_sol[constr_id] == 0
                 is_identical = false
                 break
             end
@@ -291,14 +290,12 @@ function setdualsol!(
         end    
     end
 
-
     ### else not identical to any existing dual sol
     new_dual_sol_id = Id{Constraint}(generateconstrid(form), getuid(form))
     adddualsol!(form, new_dual_sol, new_dual_sol_id)
     return (true, new_dual_sol_id)
 end
 
-
 function setcol_from_sp_primalsol!(
     masterform::Formulation, spform::Formulation, sol_id::VarId,
     name::String, duty::Duty{Variable}; lb::Float64 = 0.0,

src/MathProg/incumbents.jl

@@ -10,17 +10,17 @@ function PrimalBound(form::AbstractFormulation, val::Float64)
 end
 
 function PrimalSolution(
-    form::AbstractFormulation, sol::Dict{De,Va}, val::Float64
+    form::AbstractFormulation, decisions::Vector{De}, vals::Vector{Va}, val::Float64
 ) where {De,Va}
     Se = getobjsense(form)
-    return Coluna.Containers.Solution{Primal,Se,De,Va}(sol, val)
+    return Coluna.Containers.Solution{Primal,Se,De,Va}(decisions, vals, val)
 end
 
 function PrimalSolution(
-    form::AbstractFormulation, sol::Dict{De,Va}, bound::Coluna.Containers.Bound{Primal,Se}
+    form::AbstractFormulation, decisions::Vector{De}, vals::Vector{Va}, bound::Coluna.Containers.Bound{Primal,Se}
 ) where {Se,De,Va}
     @assert Se == getobjsense(form)
-    return Coluna.Containers.Solution{Primal,Se,De,Va}(sol, bound)
+    return Coluna.Containers.Solution{Primal,Se,De,Va}(decisions, vals, bound)
 end
 
 function DualBound(form::AbstractFormulation)
@@ -34,17 +34,17 @@ function DualBound(form::AbstractFormulation, val::Float64)
 end
 
 function DualSolution(
-    form::AbstractFormulation, sol::Dict{De,Va}, val::Float64
+    form::AbstractFormulation, decisions::Vector{De}, vals::Vector{Va}, val::Float64
 ) where {De,Va}
     Se = getobjsense(form)
-    return Coluna.Containers.Solution{Dual,Se,De,Va}(sol, val)
+    return Coluna.Containers.Solution{Dual,Se,De,Va}(decisions, vals, val)
 end
 
 function DualSolution(
-    form::AbstractFormulation, sol::Dict{De,Va}, bound::Coluna.Containers.Bound{Dual,Se}
+    form::AbstractFormulation, decisions::Vector{De}, vals::Vector{Va}, bound::Coluna.Containers.Bound{Dual,Se}
 ) where {Se,De,Va}
     @assert Se == getobjsense(form)
-    return Coluna.Containers.Solution{Dual,Se,De,Va}(sol, bound)
+    return Coluna.Containers.Solution{Dual,Se,De,Va}(decisions, vals, bound)
 end
 
 valueinminsense(b::PrimalBound{MinSense}) = b.value

src/MathProg/optimizerwrappers.jl

@@ -40,12 +40,7 @@ function retrieve_result(form::Formulation, optimizer::MoiOptimizer)
             terminationstatus != MOI.DUAL_INFEASIBLE &&
             terminationstatus != MOI.INFEASIBLE_OR_UNBOUNDED &&
             terminationstatus != MOI.OPTIMIZE_NOT_CALLED
-        fill_primal_result!(
-            optimizer, result, filter(
-                vc -> getcurisactive(form,vc) && getcurisexplicit(form,vc), 
-                getvars(form)
-            )
-        )
+        fill_primal_result!(form, optimizer, result)
         fill_dual_result!(
             optimizer, result, filter(
                 vc -> getcurisactive(form,vc) && getcurisexplicit(form,vc), 

src/MathProg/projection.jl

@@ -1,7 +1,8 @@
 projection_is_possible(master::Formulation{DwMaster}) = true
 
 function proj_cols_on_rep(sol::PrimalSolution{Sense}, master::Formulation{DwMaster}) where {Sense}
-    projected_sol = Dict{VarId, Float64}()
+    projected_sol_vars = Vector{VarId}()
+    projected_sol_vals = Vector{Float64}()
     for (mc_id, mc_val) in sol
         origin_form_uid = getoriginformuid(mc_id)
         # TODO : enhance following
@@ -14,10 +15,11 @@ function proj_cols_on_rep(sol::PrimalSolution{Sense}, master::Formulation{DwMast
         for (rep_id, rep_val) in Iterators.filter(
             v -> getduty(v) <= DwSpPricingVar || getduty(v) <= DwSpSetupVar,
             col)
-            projected_sol[rep_id] = (get!(projected_sol, rep_id, 0.0)) + rep_val * mc_val
+            push!(projected_sol_vars, rep_id)
+            push!(projected_sol_vals, rep_val * mc_val)
         end
     end
-    return PrimalSolution(master, projected_sol, float(getbound(sol)))
+    return PrimalSolution(master, projected_sol_vars, projected_sol_vals, float(getbound(sol)))
 end
 
 projection_is_possible(master::Formulation{BendersMaster}) = false

src/MathProg/vcids.jl

@@ -50,10 +50,10 @@ function Id{VC}(id::Id{VC}) where {VC}
 end
 
 Base.hash(a::Id, h::UInt) = hash(a._hash, h)
-Base.isequal(a::Id, b::Id) = Base.isequal(a._hash, b._hash)
+Base.isequal(a::Id{VC}, b::Id{VC}) where {VC} = Base.isequal(a._hash, b._hash)
 Base.isequal(a::Int, b::Id) = Base.isequal(a, b._hash)
 Base.isequal(a::Id, b::Int) = Base.isequal(a._hash, b)
-Base.isless(a::Id, b::Id) = Base.isless(a.uid, b.uid)
+Base.isless(a::Id{VC}, b::Id{VC}) where {VC} = Base.isless(a._hash, b._hash)
 Base.zero(I::Type{<:Id}) = I(-1, -1, -1, -1, -1)
 
 Base.:(<)(a::Id{VC}, b::Id{VC}) where {VC} = a._hash < b._hash

src/parameters.jl

@@ -1,9 +1,11 @@
 Base.@kwdef mutable struct Params
     max_num_nodes::Int = 10000
     open_nodes_limit::Int = 100000
-    integrality_tolerance::Float64 = 1e-5
-    absolute_optimality_tolerance::Float64 = 1e-5
-    relative_optimality_tolerance::Float64 = 1e-5
+    #integrality_tolerance::Float64 = 1e-5
+    #absolute_optimality_tolerance::Float64 = 1e-5
+    #relative_optimality_tolerance::Float64 = 1e-5
+    tol::Float64 = 1e-6 # if - ϵ_tol < val < ϵ_tol, we consider val = 0
+    tol_digits::Int = 6 # because round(val, digits = 6)
     cut_up::Float64 = Inf
     cut_lo::Float64 = -Inf
     force_copy_names::Bool = false

test/full_instances_tests.jl

@@ -73,7 +73,6 @@ function generalized_assignment_tests()
         @test CLD.GeneralizedAssignment.print_and_check_sol(data, problem, x)
     end
 
-
     @testset "gap - ColGen max nb iterations" begin
         data = CLD.GeneralizedAssignment.data("smallgap3.txt")
 

test/pricing_callback_tests.jl

@@ -17,14 +17,16 @@ function mycallback(form::CL.Formulation)
     optimize!(m)
     result = CL.OptimizationResult{CL.MinSense}()
     result.primal_bound = CL.PrimalBound(form, JuMP.objective_value(m))
-    sol = Dict{CL.VarId, Float64}()
+    solvarids = Vector{CL.VarId}()
+    solvarvals = Vector{CL.Float64}()
     for i in 1:length(x)
         val = JuMP.value(x[i])
         if val > 0.000001  || val < - 0.000001 # todo use a tolerance
-            sol[CL.getid(vars[i])] = val
+            push!(solvarids, CL.getid(vars[i]))
+            push!(solvarvals, val)
         end
     end
-    push!(result.primal_sols, CL.PrimalSolution(form, sol, result.primal_bound))
+    push!(result.primal_sols, CL.PrimalSolution(form, solvarids, solvarvals, result.primal_bound))
     CL.setfeasibilitystatus!(result, CL.FEASIBLE)
     CL.setterminationstatus!(result, CL.OPTIMAL)
     return result

test/runtests.jl

@@ -31,6 +31,7 @@ unit_tests()
 @testset "Full instances " begin
     full_instances_tests()
 end
+
 @testset "Preprocessing " begin
     preprocessing_tests()
 end

test/unit/containers/solsandbounds.jl

@@ -179,11 +179,16 @@ function fake_solution_factory(nbdecisions)
             i += 1
         end
     end
-    solution = Dict{Int, Float64}()
+    dict = Dict{Int, Float64}()
+    soldecisions = Vector{Int}()
+    solvals = Vector{Float64}()
     for d in decisions
-        solution[d] = rand(rng, 0:0.0001:1000)
+        val = rand(rng, 0:0.0001:1000)
+        dict[d] = val
+        push!(soldecisions, d)
+        push!(solvals, val)
     end
-    return solution
+    return dict, soldecisions, solvals
 end
 
 function test_solution_iterations(solution::Coluna.Containers.Solution, dict::Dict)
@@ -208,16 +213,16 @@ function solution_unit()
     PrimalSolution{S} = Coluna.Containers.Solution{Primal,S,Int,Float64}
     DualSolution{S} = Coluna.Containers.Solution{Dual,S,Int,Float64}
 
-    dict_sol = fake_solution_factory(100)
-    primal_sol = PrimalSolution{MinSense}(dict_sol, 12.3)
+    dict_sol, soldecs, solvals = fake_solution_factory(100)
+    primal_sol = PrimalSolution{MinSense}(soldecs, solvals, 12.3)
     test_solution_iterations(primal_sol, dict_sol)
     @test Coluna.Containers.getvalue(primal_sol) == 12.3
     Coluna.Containers.setvalue!(primal_sol, 123.4)
     @test Coluna.Containers.getvalue(primal_sol) == 123.4
     @test typeof(Coluna.Containers.getbound(primal_sol)) == Coluna.Containers.Bound{Primal,MinSense}
 
-    dict_sol = fake_solution_factory(100)
-    dual_sol = DualSolution{MaxSense}(dict_sol, 32.1)
+    dict_sol, soldecs, solvals = fake_solution_factory(100)
+    dual_sol = DualSolution{MaxSense}(soldecs, solvals, 32.1)
     test_solution_iterations(dual_sol, dict_sol)
     @test Coluna.Containers.getvalue(dual_sol) == 32.1
     Coluna.Containers.setvalue!(dual_sol, 432.1)