Replace dynamic sparse vector by a static one in Solution. (#710)

* wip

* ok

* better comment

* Update Project.toml

* missing SparseArrays dep

Project.toml

@@ -16,6 +16,7 @@ Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RandomNumbers = "e6cf234a-135c-5ec9-84dd-332b85af5143"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
@@ -24,7 +25,7 @@ TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 BlockDecomposition = "1.10"
 Crayons = "4.1"
 DataStructures = "0.17, 0.18"
-DynamicSparseArrays = "0.5.3"
+DynamicSparseArrays = "0.6"
 MathOptInterface = "0.10, 1"
 Parameters = "0.12"
 RandomNumbers = "1.5"

src/Algorithm/benders.jl

@@ -105,7 +105,7 @@ function update_benders_sp_problem!(
     for (varid, var) in getvars(spform)
         iscuractive(spform, varid) || continue
         getduty(varid) <= BendSpSlackFirstStageVar || continue
-        haskey(master_primal_sol, varid) || continue
+        !iszero(master_primal_sol[varid]) || continue
         setcurub!(spform, var, getperenub(spform, var) - master_primal_sol[varid])
     end
 
@@ -428,12 +428,14 @@ function generatecuts!(
 )::Tuple{Int, Bool, PrimalBound}
     masterform = getmaster(reform)
     S = getobjsense(masterform)
-    filtered_dual_sol = filter(elem -> getduty(elem[1]) == MasterPureConstr, master_dual_sol)
+
+    # following variable is not used:
+    #filtered_dual_sol = filter(elem -> getduty(elem[1]) == MasterPureConstr, master_dual_sol)
 
     ## TODO stabilization : move the following code inside a loop
     nb_new_cuts, spsols_relaxed, pb_correction, sp_pb_contrib =
         solve_sps_to_gencuts!(
-            algo, env, algdata, reform, master_primal_sol, filtered_dual_sol, phase
+            algo, env, algdata, reform, master_primal_sol, master_dual_sol, phase
         )
     update_lagrangian_pb!(algdata, reform, master_dual_sol, sp_pb_contrib)
     if nb_new_cuts < 0

src/ColunaBase/ColunaBase.jl

@@ -2,7 +2,7 @@ module ColunaBase
 
 using ..Coluna
 
-using DynamicSparseArrays, MathOptInterface, TimerOutputs, RandomNumbers, Random
+using DynamicSparseArrays, MathOptInterface, TimerOutputs, RandomNumbers, Random, SparseArrays
 
 const MOI = MathOptInterface
 const TO = TimerOutputs

src/ColunaBase/hashtable.jl

@@ -5,35 +5,57 @@ const MT_MASK = 0x0ffff # hash keys from 1 to 65536
 This datastructure allows us to quickly find solution that shares the same members:
 variables for primal solutions and constraints for dual solutions.
 """
-struct HashTable{VarConstrId}
+struct HashTable{MemberId,SolId}
     rng::MersenneTwisters.MT19937
-    memberid_to_hash::Dict{VarConstrId, UInt32} # members of the primal/dual solution -> hash
-    hash_to_solids::Vector{Vector{VarConstrId}} # hash of the primal/dual solution -> solution id
+    memberid_to_hash::Dict{MemberId, UInt32} # members of the primal/dual solution -> hash
+    hash_to_solids::Vector{Vector{SolId}} # hash of the primal/dual solution -> solution id
 
-    HashTable{VarConstrId}() where {VarConstrId} = new(
+    HashTable{MemberId,SolId}() where {MemberId,SolId} = new(
         MersenneTwisters.MT19937(MT_SEED),
-        Dict{VarConstrId, UInt32}(),
-        [VarConstrId[] for _ in 0:MT_MASK]
+        Dict{MemberId, UInt32}(),
+        [SolId[] for _ in 0:MT_MASK]
     )
 end
 
+function _gethash!(
+    hashtable::HashTable{MemberId,SolId}, id::MemberId, bad_hash = Int(MT_MASK) + 2
+) where {MemberId,SolId}
+    hash = UInt32(get(hashtable.memberid_to_hash, id, bad_hash) - 1)
+    if hash > MT_MASK
+        hash = MersenneTwisters.mt_get(hashtable.rng) & MT_MASK
+        hashtable.memberid_to_hash[id] = Int(hash) + 1
+    end
+    return hash
+end
+
+_gethash!(hashtable, entry::Tuple, bad_hash = Int(MT_MASK) + 2) = 
+    _gethash!(hashtable, first(entry), bad_hash)
+
+# By default, we consider that the iterator of the `sol` argument returns a tuple that 
+# contains the id as first element.
 function gethash(hashtable::HashTable, sol)
-    bad_hash = Int(MT_MASK) + 2
     acum_hash = UInt32(0)
-    for (varconstrid, _) in sol
-        hash = UInt32(get(hashtable.memberid_to_hash, varconstrid, bad_hash) - 1)
-        if hash > MT_MASK
-            hash = MersenneTwisters.mt_get(hashtable.rng) & MT_MASK
-            hashtable.memberid_to_hash[varconstrid] = Int(hash) + 1
-        end
-        acum_hash ⊻= hash
+    for entry in sol
+        acum_hash ⊻= _gethash!(hashtable, entry)
+    end
+    return Int(acum_hash) + 1
+end
+
+# If the solution is in a sparse vector, we just want to check indices associated to non-zero
+# values.
+function gethash(hashtable::HashTable, sol::SparseVector)
+    acum_hash = UInt32(0)
+    for nzid in SparseArrays.nonzeroinds(sol)
+        acum_hash ⊻= _gethash!(hashtable, nzid)
     end
     return Int(acum_hash) + 1
 end
 
-savesolid!(hashtable::HashTable, solid, sol) = push!(getsolids(hashtable, sol), solid)
+savesolid!(hashtable::HashTable, solid, sol) =
+    push!(getsolids(hashtable, sol), solid)
 
-getsolids(hashtable::HashTable, sol) = hashtable.hash_to_solids[gethash(hashtable, sol)]
+getsolids(hashtable::HashTable, sol) =
+    hashtable.hash_to_solids[gethash(hashtable, sol)]
 
 function Base.show(io::IO, ht::HashTable)
     println(io, typeof(ht), ":")

src/ColunaBase/solsandbounds.jl

@@ -257,11 +257,11 @@ function convert_status(coluna_status::SolutionStatus)
 end
 
 # Basic structure of a solution
-struct Solution{Model<:AbstractModel,Decision,Value} <: AbstractDict{Decision,Value}
+struct Solution{Model<:AbstractModel,Decision<:Integer,Value} <: AbstractSparseVector{Decision,Value}
     model::Model
     bound::Float64
     status::SolutionStatus
-    sol::DynamicSparseArrays.PackedMemoryArray{Decision,Value}
+    sol::SparseVector{Value,Decision}
 end
 
 """
@@ -283,10 +283,9 @@ Create a solution to the `model`. Other arguments are:
 - `status` is the solution status.
 """
 function Solution{Mo,De,Va}(
-    model::Mo, decisions::Vector{De}, values::Vector{Va}, solution_value::Float64, 
-    status::SolutionStatus
+    model::Mo, decisions::Vector{De}, values::Vector{Va}, solution_value::Float64, status::SolutionStatus
 ) where {Mo<:AbstractModel,De,Va}
-    sol = DynamicSparseArrays.dynamicsparsevec(decisions, values)
+    sol = sparsevec(decisions, values, typemax(De))
     return Solution(model, solution_value, status, sol)
 end
 
@@ -302,22 +301,56 @@ getvalue(s::Solution) = float(s.bound)
 "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)
+# implementing indexing interface
+Base.getindex(s::Solution, i::Integer) = getindex(s.sol, i)
+Base.setindex!(s::Solution, v, i::Integer) = setindex!(s.sol, v, i)
+Base.firstindex(s::Solution) = firstindex(s.sol)
+Base.lastindex(s::Solution) = lastindex(s.sol)
+
+# implementing abstract array interface
+Base.size(s::Solution) = size(s.sol)
 Base.length(s::Solution) = length(s.sol)
-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
+Base.IndexStyle(::Type{<:Solution{Mo,De,Va}}) where {Mo,De,Va} =
+    IndexStyle(SparseVector{Va,De})
+SparseArrays.nnz(s::Solution) = nnz(s.sol)
+
+# It iterates only on non-zero values because:
+# - we use indices (`Id`) that behaves like an Int with additional information and given a 
+#   indice, we cannot deduce the additional information for the next one (i.e. impossible to
+#   create an Id for next integer);
+# - we don't know the length of the vector (it depends on the number of variables & 
+#   constraints that varies over time).
+function Base.iterate(s::Solution)
+    iterator = Iterators.zip(findnz(s.sol)...)
+    next = iterate(iterator)
+    isnothing(next) && return nothing
+    (item, zip_state) = next
+    return (item, (zip_state, iterator))
+end
+
+function Base.iterate(::Solution, state)
+    (zip_state, iterator) = state
+    next = iterate(iterator, zip_state)
+    isnothing(next) && return nothing
+    (next_item, next_zip_state) = next
+    return (next_item, (next_zip_state, iterator))
+end
+
+# # implementing sparse array interface
+# SparseArrays.nnz(s::Solution) = nnz(s.sol)
+# SparseArrays.nonzeroinds(s::Solution) = SparseArrays.nonzeroinds(s.sol)
+# SparseArrays.nonzeros(s::Solution) = nonzeros(s.sol)
+
+function _eq_sparse_vec(a::SparseVector, b::SparseVector)
+    a_ids, a_vals = findnz(a)
+    b_ids, b_vals = findnz(b)
+    return a_ids == b_ids && a_vals == b_vals
+end
 
 Base.:(==)(::Solution, ::Solution) = false
 function Base.:(==)(a::S, b::S) where {S<:Solution}
     return a.model == b.model && a.bound == b.bound && a.status == b.status &&
-            a.sol == b.sol
-end
-
-# 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))
+        _eq_sparse_vec(a.sol, b.sol)
 end
 
 function Base.in(p::Tuple{De,Va}, a::Solution{Mo,De,Va}, valcmp=(==)) where {Mo,De,Va}
@@ -328,12 +361,23 @@ function Base.in(p::Tuple{De,Va}, a::Solution{Mo,De,Va}, valcmp=(==)) where {Mo,
     return false
 end
 
-function Base.show(io::IO, solution::Solution{Mo,De,Va}) where {Mo,De,Va}
+function Base.show(io::IOContext, solution::Solution{Mo,De,Va}) where {Mo,De,Va}
     println(io, "Solution")
     for (decision, value) in solution
         println(io, "| ", decision, " = ", value)
     end
     Printf.@printf(io, "└ value = %.2f \n", getvalue(solution))
 end
+
 # Todo : revise method
 Base.copy(s::S) where {S<:Solution} = S(s.bound, copy(s.sol))
+
+# Implementing comparison between solution & dynamic matrix col view for solution comparison
+function Base.:(==)(v1::DynamicMatrixColView, v2::Solution)
+    for ((i1,j1), (i2,j2)) in Iterators.zip(v1,v2)
+        if !(i1 == i2 && j1 == j2)
+            return false
+        end
+    end
+    return true
+end

src/MathProg/MathProg.jl

@@ -9,7 +9,7 @@ using ..ColunaBase
 
 import Base: haskey, length, iterate, diff, delete!, contains, setindex!, getindex, view
 
-using DynamicSparseArrays, Logging, Printf
+using DynamicSparseArrays, SparseArrays, Logging, Printf
 
 const BD = BlockDecomposition
 const ClB = ColunaBase

src/MathProg/duties.jl

@@ -28,7 +28,7 @@ mutable struct DwSp <: AbstractSpDuty
     ## [colid, varid] = value
     primalsols_pool::VarVarMatrix
     # Hash table to quickly find identical solutions
-    hashtable_primalsols_pool::HashTable{VarId}
+    hashtable_primalsols_pool::HashTable{VarId,VarId}
     ## Perennial cost of solutions
     costs_primalsols_pool::Dict{VarId, Float64}
     ## Custom representation of solutions
@@ -40,7 +40,7 @@ function DwSp(setup_var, lower_multiplicity, upper_multiplicity, column_var_kind
     return DwSp(
         setup_var, lower_multiplicity, upper_multiplicity, column_var_kind,
         dynamicsparse(VarId, VarId, Float64; fill_mode = false),
-        HashTable{VarId}(),
+        HashTable{VarId, VarId}(),
         Dict{VarId, Float64}(),
         Dict{VarId, BD.AbstractCustomData}()
     )

src/MathProg/formulation.jl

@@ -255,8 +255,7 @@ get_primal_sol_pool_hash_table(form::Formulation{DwSp}) = form.duty_data.hashtab
 function _get_same_sol_in_pool(pool_sols, pool_hashtable, sol)
     sols_with_same_members = getsolids(pool_hashtable, sol)
     for existing_sol_id in sols_with_same_members
-        # TODO: implement comparison between view & dynamicsparsevec
-        existing_sol = pool_sols[existing_sol_id,:]
+        existing_sol = @view pool_sols[existing_sol_id,:]
         if existing_sol == sol
             return existing_sol_id
         end
@@ -316,7 +315,7 @@ function get_column_from_pool(primal_sol::PrimalSolution{Formulation{DwSp}})
     spform = primal_sol.solution.model
     pool = get_primal_sol_pool(spform)
     pool_hashtable = get_primal_sol_pool_hash_table(spform)
-    return _get_same_sol_in_pool(pool, pool_hashtable, primal_sol.solution.sol)
+    return _get_same_sol_in_pool(pool, pool_hashtable, primal_sol)
 end
 
 """

src/MathProg/manager.jl

@@ -1,13 +1,11 @@
-const DynSparseVector{I} = DynamicSparseArrays.PackedMemoryArray{I, Float64}
-
 const VarMembership = Dict{VarId, Float64}
 const ConstrMembership = Dict{ConstrId, Float64}
 const ConstrConstrMatrix = DynamicSparseArrays.DynamicSparseMatrix{ConstrId,ConstrId,Float64}
 const VarConstrDualSolMatrix = DynamicSparseArrays.DynamicSparseMatrix{VarId,ConstrId,Tuple{Float64,ActiveBound}}
 const VarVarMatrix = DynamicSparseArrays.DynamicSparseMatrix{VarId,VarId,Float64}
 
 # Define the semaphore of the dynamic sparse matrix using MathProg.Id as index
-DynamicSparseArrays.semaphore_key(::Type{I}) where {I <: Id} = zero(I)
+DynamicSparseArrays.semaphore_key(I::Type{Id{VC}}) where VC = I(Duty{VC}(0), -1, -1, -1, -1)
 
 # We wrap the coefficient matrix because we need to buffer the changes.
 struct CoefficientMatrix{C,V,T}
@@ -51,7 +49,7 @@ mutable struct FormulationManager
     coefficients::ConstrVarMatrix # rows = constraints, cols = variables
     dual_sols::ConstrConstrMatrix # cols = dual solutions with constrid, rows = constrs
     dual_sols_varbounds::VarConstrDualSolMatrix # cols = dual solutions with constrid, rows = variables
-    dual_sol_rhss::DynSparseVector{ConstrId} # dual solutions with constrid map to their rhs
+    dual_sol_rhss::DynamicSparseVector{ConstrId} # dual solutions with constrid map to their rhs
     robust_constr_generators::Vector{RobustConstraintsGenerator}
     custom_families_id::Dict{DataType,Int}
 end

src/MathProg/solutions.jl

@@ -102,10 +102,11 @@ ColunaBase.getmodel(s::AbstractSolution) = getmodel(s.solution)
 ColunaBase.getvalue(s::AbstractSolution) = getvalue(s.solution)
 ColunaBase.getbound(s::AbstractSolution) = getbound(s.solution)
 ColunaBase.getstatus(s::AbstractSolution) = getstatus(s.solution)
+
 Base.length(s::AbstractSolution) = length(s.solution)
-Base.get(s::AbstractSolution, id, default) = Base.get(s.solution, id, default)
-Base.getindex(s::AbstractSolution, id) = Base.getindex(s.solution, id)
-Base.setindex!(s::AbstractSolution, val, id) = Base.setindex!(s.solution, val, id)
+Base.get(s::AbstractSolution, id, default) = get(s.solution, id, default)
+Base.getindex(s::AbstractSolution, id) = getindex(s.solution, id)
+Base.setindex!(s::AbstractSolution, val, id) = setindex!(s.solution, val, id)
 
 # Iterating over a PrimalSolution or a DualSolution is similar to iterating over
 # ColunaBase.Solution
@@ -191,13 +192,11 @@ function Base.show(io::IO, solution::PrimalSolution{M}) where {M}
     Printf.@printf(io, "└ value = %.2f \n", getvalue(solution))
 end
 
-# Following methods are needed by Benders
-# TODO : check if we can remove them during refactoring of Benders
-# not performant
-Base.haskey(s::AbstractSolution, key) = haskey(s.solution, key)
-# we can't filter the constraints, the variables, and the custom data.
-function Base.filter(f::Function, s::DualSolution)
-    return DualSolution(
-        filter(f, s.solution), s.var_redcosts, s.custom_data
-    )
-end
+
+# To check if a solution is part of solutions from the pool.
+Base.:(==)(v1::DynamicMatrixColView, v2::AbstractSolution) = v1 == v2.solution
+
+# To allocate an array with size equals to the number of non-zero elements when using
+# "generation" syntax.
+Base.length(gen::Base.Generator{<:AbstractSolution}) = nnz(gen.iter.solution)
+

src/MathProg/vcids.jl

@@ -49,12 +49,25 @@ function Id{VC}(
 end
 
 Base.hash(a::Id, h::UInt) = hash(a.uid, h)
-Base.zero(I::Type{Id{VC}}) where {VC} = I(Duty{VC}(0), -1, -1, -1, -1) # semaphore in the PMA (DynamicSparseArrays).
-Base.isequal(a::Id{VC}, b::Id{VC}) where {VC} = Base.isequal(a.uid, b.uid)
+Base.zero(I::Type{Id{VC}}) where {VC} = I(Duty{VC}(0), zero(Int32), -1, -1, -1)
+Base.zero(::Id{VC}) where {VC} = Id{VC}(Duty{VC}(0), zero(Int32), -1, -1, -1)
+Base.one(I::Type{Id{VC}}) where {VC} = I(Duty{VC}(0), one(Int32), -1, -1, -1)
+Base.typemax(I::Type{Id{VC}}) where {VC} = I(Duty{VC}(0), typemax(Int32), -1, -1, -1)
+Base.isequal(a::Id{VC}, b::Id{VC}) where {VC} = isequal(a.uid, b.uid)
 
 Base.promote_rule(::Type{T}, ::Type{<:Id}) where {T<:Integer} = T
 Base.promote_rule(::Type{<:Id}, ::Type{T}) where {T<:Integer} = T
-Base.promote_rule(::Type{I}, ::Type{I}) where {I<:Id} = Int32
+Base.promote_rule(::Type{<:Id}, ::Type{<:Id}) = Int32
+
+# Promotion mechanism will never call the following rule:
+#   Base.promote_rule(::Type{I}, ::Type{I}) where {I<:Id} = Int32
+#
+# The problem is that an Id is an integer with additional information and we
+# cannot generate additional information of a new id from the operation of two
+# existing ids.
+# As we want that all operations on ids results on operations on the uid,
+# we redefine the promotion mechanism for Ids so that operations on Ids return integer:
+Base.promote_type(::Type{I}, ::Type{I}) where {I<:Id} = Int32
 
 Base.convert(::Type{Int}, id::I) where {I<:Id} = Int(id.uid)
 Base.convert(::Type{Int32}, id::I) where {I<:Id} = id.uid

test/ColunaTests.jl

@@ -1,6 +1,6 @@
 module ColunaTests
     using Base.CoreLogging: error
-    using DynamicSparseArrays, Coluna, TOML
+    using DynamicSparseArrays, SparseArrays, Coluna, TOML
 
     using ReTest, GLPK, ColunaDemos, JuMP, BlockDecomposition, Random, MathOptInterface, MathOptInterface.Utilities, Base.CoreLogging, Logging
     global_logger(ConsoleLogger(stderr, LogLevel(0)))