Memory estimation (#22)

* update

* estimate memory

Author: GiggleLiu

README.md

@@ -52,7 +52,7 @@ end
 
 ## Examples
 
-You can find many examples in the documentation, a good one to start with is [solving the independent set problem](https://psychic-meme-f4d866f8.pages.github.io/dev/tutorials/IndependentSet/)
+You can find many examples in the documentation, a good one to start with is [solving the independent set problem](https://psychic-meme-f4d866f8.pages.github.io/dev/tutorials/IndependentSet/).
 
 ## Supporting and Citing
 

docs/src/performancetips.md

@@ -49,6 +49,17 @@ julia> sizeof(TruncatedPoly{5,Float64,Float64})
 48
 ```
 
+One can use [`estimate_memory`](@ref) to get a good estimation of peak memory in bytes.
+```julia
+julia> estimate_memory(problem, GraphPolynomial(; method=:finitefield))
+297616
+
+julia> estimate_memory(problem, GraphPolynomial(; method=:polynomial))
+71427840
+```
+It means one only needs 298 KB memory to find the graph polynomial with the finite field approach,
+but needs 71 MB memory to find the graph polynomial using the [`Polynomial`](@ref) type.
+
 !!! note
     * The actual run time memory can be several times larger than the size of the maximum tensor.
     There is no constant bound for the factor, an empirical value for it is 3x.

docs/src/ref.md

@@ -105,6 +105,7 @@ getixsv
 getiyv
 timespace_complexity
 timespacereadwrite_complexity
+estimate_memory
 @ein_str
 GreedyMethod
 TreeSA

src/GraphTensorNetworks.jl

@@ -11,6 +11,9 @@ using Graphs
 export timespace_complexity, timespacereadwrite_complexity, @ein_str, getixsv, getiyv
 export GreedyMethod, TreeSA, SABipartite, KaHyParBipartite, MergeVectors, MergeGreedy
 
+# estimate memory
+export estimate_memory
+
 # Algebras
 export StaticBitVector, StaticElementVector, @bv_str
 export is_commutative_semiring

src/configurations.jl

@@ -20,19 +20,19 @@ Find optimal solutions with bounding.
 * If `invert` is true, find the minimum.
 * If `tree_storage` is true, use [`TreeConfigEnumerator`](@ref) as the storage of solutions.
 """
-function best_solutions(gp::GraphProblem; all=false, usecuda=false, invert=false, tree_storage::Bool=false)
+function best_solutions(gp::GraphProblem; all=false, usecuda=false, invert=false, tree_storage::Bool=false, T=Float64)
     if all && usecuda
         throw(ArgumentError("ConfigEnumerator can not be computed on GPU!"))
     end
-    xst = generate_tensors(_x(TropicalF64; invert), gp)
+    xst = generate_tensors(_x(Tropical{T}; invert), gp)
     ymask = trues(fill(2, length(getiyv(gp.code)))...)
     if usecuda
         xst = CuArray.(xst)
         ymask = CuArray(ymask)
     end
     if all
-        # we use `Float64` types because we want to support weighted graphs.
-        T = config_type(CountingTropical{Float64,Float64}, length(labels(gp)), nflavor(gp); all, tree_storage)
+        # we use `Float64` as default because we want to support weighted graphs.
+        T = config_type(CountingTropical{T,T}, length(labels(gp)), nflavor(gp); all, tree_storage)
         xs = generate_tensors(_x(T; invert), gp)
         ret = bounding_contract(AllConfigs{1}(), gp.code, xst, ymask, xs)
         return invert ? post_invert_exponent.(ret) : ret
@@ -71,12 +71,12 @@ end
 
 Finding optimal and suboptimal solutions.
 """
-best2_solutions(gp::GraphProblem; all=true, usecuda=false, invert::Bool=false) = solutions(gp, Max2Poly{Float64,Float64}; all, usecuda, invert)
+best2_solutions(gp::GraphProblem; all=true, usecuda=false, invert::Bool=false, T=Float64) = solutions(gp, Max2Poly{T,T}; all, usecuda, invert)
 
-function bestk_solutions(gp::GraphProblem, k::Int; invert::Bool=false, tree_storage::Bool=false)
-    xst = generate_tensors(_x(TropicalF64; invert), gp)
+function bestk_solutions(gp::GraphProblem, k::Int; invert::Bool=false, tree_storage::Bool=false, T=Float64)
+    xst = generate_tensors(_x(Tropical{T}; invert), gp)
     ymask = trues(fill(2, length(getiyv(gp.code)))...)
-    T = config_type(TruncatedPoly{k,Float64,Float64}, length(labels(gp)), nflavor(gp); all=true, tree_storage)
+    T = config_type(TruncatedPoly{k,T,T}, length(labels(gp)), nflavor(gp); all=true, tree_storage)
     xs = generate_tensors(_x(T; invert), gp)
     ret = bounding_contract(AllConfigs{k}(), gp.code, xst, ymask, xs)
     return invert ? post_invert_exponent.(ret) : ret
@@ -88,7 +88,7 @@ end
 Finding all solutions grouped by size.
 e.g. when the problem is [`MaximalIS`](@ref), it computes all maximal independent sets, or the maximal cliques of it complement.
 """
-all_solutions(gp::GraphProblem) = solutions(gp, Polynomial{Float64,:x}, all=true, usecuda=false, tree_storage=false)
+all_solutions(gp::GraphProblem; T=Float64) = solutions(gp, Polynomial{T,:x}, all=true, usecuda=false, tree_storage=false)
 
 function _onehotv(::Type{Polynomial{BS,X}}, x, v) where {BS,X}
     Polynomial{BS,X}([onehotv(BS, x, v)])

src/graph_polynomials.jl

@@ -4,10 +4,12 @@ using FFTW
 using Graphs
 
 """
-    graph_polynomial(problem, method; usecuda=false, kwargs...)
+    graph_polynomial(problem, method; usecuda=false, T=Float64, kwargs...)
 
 Computing the graph polynomial for specific problem.
 
+Positional Arguments
+----------------------------------
 * `problem` can be one of the following instances,
     * `IndependentSet` for the independence polynomial,
     * `MaximalIS` for the maximal independence polynomial,
@@ -23,27 +25,32 @@ Computing the graph polynomial for specific problem.
         It Consumes additional kwargs [`maxorder`, `r`]. The larger `r` is,
         the more accurate the factors of high order terms, and the less accurate the factors of low order terms.
     * `Val(:fitting)`, compute with the polynomial fitting approach, fast but inaccurate for large graphs.
+
+Keyword Arguments
+----------------------------------
+* `usecuda` is true if one wants to compute with CUDA arrays,
+* `T` is the base type
 """
 function graph_polynomial end
 
-function graph_polynomial(gp::GraphProblem, ::Val{:fft}; usecuda=false, 
+function graph_polynomial(gp::GraphProblem, ::Val{:fft}; usecuda=false, T=Float64,
         maxorder=max_size(gp; usecuda=usecuda), r=1.0)
     ω = exp(-2im*π/(maxorder+1))
     xs = r .* collect(ω .^ (0:maxorder))
-    ys = [Array(contractx(gp, x; usecuda=usecuda)) for x in xs]
+    ys = [Array(contractx(gp, Complex{T}(x); usecuda=usecuda)) for x in xs]
     map(ci->Polynomial(ifft(getindex.(ys, Ref(ci))) ./ (r .^ (0:maxorder))), CartesianIndices(ys[1]))
 end
 
-function graph_polynomial(gp::GraphProblem, ::Val{:fitting}; usecuda=false,
+function graph_polynomial(gp::GraphProblem, ::Val{:fitting}; usecuda=false, T=Float64,
         maxorder = max_size(gp; usecuda=usecuda))
     xs = (0:maxorder)
-    ys = [Array(contractx(gp, x; usecuda=usecuda)) for x in xs]
+    ys = [Array(contractx(gp, T(x); usecuda=usecuda)) for x in xs]
     map(ci->fit(xs, getindex.(ys, Ref(ci))), CartesianIndices(ys[1]))
 end
 
-function graph_polynomial(gp::GraphProblem, ::Val{:polynomial}; usecuda=false)
+function graph_polynomial(gp::GraphProblem, ::Val{:polynomial}; usecuda=false, T=Float64)
     @assert !usecuda "Polynomial type can not be computed on GPU!"
-    contractx(gp::GraphProblem, Polynomial([0, 1.0]))
+    contractx(gp::GraphProblem, Polynomial(T[0, 1]))
 end
 
 function _polynomial_single(gp::GraphProblem, ::Type{T}; usecuda, maxorder) where T
@@ -60,7 +67,8 @@ function _polynomial_single(gp::GraphProblem, ::Type{T}; usecuda, maxorder) wher
     return res
 end
 
-function graph_polynomial(gp::GraphProblem, ::Val{:finitefield}; usecuda=false,
+# T is not used in finitefield approach
+function graph_polynomial(gp::GraphProblem, ::Val{:finitefield}; usecuda=false, T=Float64,
         maxorder=max_size(gp; usecuda=usecuda), max_iter=100)
     return map(Polynomial, big_integer_solve(T->_polynomial_single(gp, T; usecuda=usecuda, maxorder=maxorder), Int32, max_iter))
 end

src/interfaces.jl

@@ -90,6 +90,10 @@ Method Argument
     * It accepts keyword arguments `maxorder` (optional) and `r`,
         if `r > 1`, one has better precision for coefficients of large order, if `r < 1`,
         one has better precision for coefficients of small order.
+* `:fitting`, fit the polynomial directly.
+    * The corresponding tensor element type is floating point numbers like `Base.Float64`.
+    * It has round-off error.
+    * BLAS and GPU are supported, it is the fastest among all methods.
 
 Graph polynomials are not defined for weighted graph problems.
 """
@@ -237,7 +241,7 @@ function solve(gp::GraphProblem, property::AbstractProperty; T=Float64, usecuda=
     elseif property isa CountingMin
         return post_invert_exponent.(contractx(gp, pre_invert_exponent(TruncatedPoly(ntuple(i->i == min_k(property) ? one(T) : zero(T), min_k(property)), one(T))); usecuda=usecuda))
     elseif property isa GraphPolynomial
-        return graph_polynomial(gp, Val(graph_polynomial_method(property)); usecuda=usecuda, property.kwargs...)
+        return graph_polynomial(gp, Val(graph_polynomial_method(property)); usecuda=usecuda, T=T, property.kwargs...)
     elseif property isa SingleConfigMax{false}
         return solutions(gp, CountingTropical{T,T}; all=false, usecuda=usecuda, )
     elseif property isa SingleConfigMin{false}
@@ -253,19 +257,19 @@ function solve(gp::GraphProblem, property::AbstractProperty; T=Float64, usecuda=
     elseif property isa ConfigsAll
         return solutions(gp, Real; all=true, usecuda=usecuda, tree_storage=tree_storage(property))
     elseif property isa SingleConfigMax{true}
-        return best_solutions(gp; all=false, usecuda=usecuda)
+        return best_solutions(gp; all=false, usecuda=usecuda, T=T)
     elseif property isa SingleConfigMin{true}
-        return best_solutions(gp; all=false, usecuda=usecuda, invert=true)
+        return best_solutions(gp; all=false, usecuda=usecuda, invert=true, T=T)
     elseif property isa ConfigsMax{1,true}
-        return best_solutions(gp; all=true, usecuda=usecuda, tree_storage=tree_storage(property))
+        return best_solutions(gp; all=true, usecuda=usecuda, tree_storage=tree_storage(property), T=T)
     elseif property isa ConfigsMin{1,true}
-        return best_solutions(gp; all=true, usecuda=usecuda, invert=true, tree_storage=tree_storage(property))
+        return best_solutions(gp; all=true, usecuda=usecuda, invert=true, tree_storage=tree_storage(property), T=T)
     elseif property isa (ConfigsMax{K,true} where K)
-        return bestk_solutions(gp, max_k(property), tree_storage=tree_storage(property))
+        return bestk_solutions(gp, max_k(property), tree_storage=tree_storage(property), T=T)
     elseif property isa (ConfigsMin{K,true} where K)
-        return bestk_solutions(gp, min_k(property), invert=true, tree_storage=tree_storage(property))
+        return bestk_solutions(gp, min_k(property), invert=true, tree_storage=tree_storage(property), T=T)
     else
-        error("unknown property $property.")
+        error("unknown property: `$property`.")
     end
 end
 
@@ -379,3 +383,61 @@ end
 # convert to Matrix
 Base.Matrix(ce::ConfigEnumerator) = plain_matrix(ce)
 Base.Vector(ce::StaticElementVector) = collect(ce)
+
+########## memory estimation ###############
+"""
+    estimate_memory(problem, property; T=Float64)
+
+Memory estimation in number of bytes to compute certain `property` of a `problem`.
+`T` is the base type.
+"""
+function estimate_memory(problem::GraphProblem, property::AbstractProperty; T=Float64)
+    _estimate_memory(tensor_element_type(T, length(labels(problem)), nflavor(problem), property), problem)
+end
+function estimate_memory(problem::GraphProblem, ::Union{SingleConfigMax{true},SingleConfigMin{true}}; T=Float64)
+    tc, sc, rw = timespacereadwrite_complexity(problem.code, _size_dict(problem))
+    # caching all tensors is equivalent to counting the total number of writes
+    return ceil(Int, exp2(rw - 1)) * sizeof(Tropical{T})
+end
+function estimate_memory(problem::GraphProblem, ::GraphPolynomial{:polynomial}; T=Float64)
+    # this is the upper bound
+    return peak_memory(problem.code, _size_dict(problem)) * (sizeof(T) * length(labels(problem)))
+end
+
+function _size_dict(problem)
+    lbs = labels(problem)
+    nf = nflavor(problem)
+    return Dict([lb=>nf for lb in lbs])
+end
+
+function _estimate_memory(::Type{ET}, problem::GraphProblem) where ET
+    if !isbitstype(ET) && !(ET <: Mod)
+        @warn "Target tensor element type `$ET` is not a bits type, the estimation of memory might be unreliable."
+    end
+    return peak_memory(problem.code, _size_dict(problem)) * sizeof(ET)
+end
+
+for (PROP, ET) in [(:SizeMax, :(Tropical{T})), (:SizeMin, :(Tropical{T})),
+        (:(SingleConfigMax{true}), :(Tropical{T})), (:(SingleConfigMin{true}), :(Tropical{T})),
+        (:(CountingAll), :T), (:(CountingMax{1}), :(CountingTropical{T,T})), (:(CountingMin{1}), :(CountingTropical{T,T})),
+        (:(CountingMax{K}), :(TruncatedPoly{K,T,T})), (:(CountingMin{K}), :(TruncatedPoly{K,T,T})),
+        (:(GraphPolynomial{:finitefield}), :(Mod{N,Int32} where N)), (:(GraphPolynomial{:fft}), :(Complex{T})), 
+        (:(GraphPolynomial{:polynomial}), :(Polynomial{T, :x})), (:(GraphPolynomial{:fitting}), :T),
+    ]
+    @eval tensor_element_type(::Type{T}, n::Int, nflavor::Int, ::$PROP) where {T,K} = $ET
+end
+
+for (PROP, ET) in [(:(SingleConfigMax{false}), :(CountingTropical{T,T})), (:(SingleConfigMin{false}), :(CountingTropical{T,T}))]
+    @eval function tensor_element_type(::Type{T}, n::Int, nflavor::Int, ::$PROP) where {T}
+        sampler_type($ET, n, nflavor)
+    end
+end
+for (PROP, ET) in [
+        (:(ConfigsMax{1}), :(CountingTropical{T,T})), (:(ConfigsMin{1}), :(CountingTropical{T,T})),
+        (:(ConfigsMax{K}), :(TruncatedPoly{K,T,T})), (:(ConfigsMin{K}), :(TruncatedPoly{K,T,T})),
+        (:(ConfigsAll), :(Real))
+    ]
+    @eval function tensor_element_type(::Type{T}, n::Int, nflavor::Int, ::$PROP) where {T,K}
+        set_type($ET, n, nflavor)
+    end
+end

test/interfaces.jl

@@ -179,4 +179,31 @@ end
         @test length(res1) == length(res2)
         @test Set(res2 |> collect) == Set(res1 |> collect)
     end
+end
+
+@testset "memory estimation" begin
+    gp = IndependentSet(smallgraph(:petersen))
+    for property in [
+            SizeMax(), SizeMin(), CountingMax(), CountingMin(), CountingMax(2), CountingMin(2),
+            ConfigsMax(;bounded=true), ConfigsMin(;bounded=true), ConfigsMax(2;bounded=true), ConfigsMin(2;bounded=true), 
+            ConfigsMax(;bounded=false), ConfigsMin(;bounded=false), ConfigsMax(2;bounded=false), ConfigsMin(2;bounded=false), SingleConfigMax(;bounded=false), SingleConfigMin(;bounded=false),
+            CountingAll(), ConfigsAll(),
+        ]
+        @show property
+        ET = GraphTensorNetworks.tensor_element_type(Float32, 10, 2, property)
+        @test eltype(solve(gp, property, T=Float32)) <: ET
+        @test estimate_memory(gp, property) isa Integer
+    end
+    @test GraphTensorNetworks.tensor_element_type(Float32, 10, 2, GraphPolynomial(method=:polynomial)) == Polynomial{Float32, :x}
+    @test sizeof(GraphTensorNetworks.tensor_element_type(Float32, 10, 2, GraphPolynomial(method=:fitting))) == 4
+    @test sizeof(GraphTensorNetworks.tensor_element_type(Float32, 10, 2, GraphPolynomial(method=:fft))) == 8
+    @test sizeof(GraphTensorNetworks.tensor_element_type(Float64, 10, 2, GraphPolynomial(method=:finitefield))) == 4
+    @test GraphTensorNetworks.tensor_element_type(Float32, 10, 2, SingleConfigMax(;bounded=true)) == Tropical{Float32}
+    @test GraphTensorNetworks.tensor_element_type(Float32, 10, 2, SingleConfigMin(;bounded=true)) == Tropical{Float32}
+
+    @test estimate_memory(gp, SizeMax()) * 2 == estimate_memory(gp, CountingMax())
+    @test estimate_memory(gp, SingleConfigMax(bounded=true)) > estimate_memory(gp, SingleConfigMax(bounded=false))
+    @test estimate_memory(gp, ConfigsMax(bounded=true)) == estimate_memory(gp, SingleConfigMax(bounded=false))
+    @test estimate_memory(gp, GraphPolynomial(method=:fitting); T=Float32) * 4 == estimate_memory(gp, GraphPolynomial(method=:fft))
+    @test estimate_memory(gp, GraphPolynomial(method=:finitefield)) * 10 == estimate_memory(gp, GraphPolynomial(method=:polynomial); T=Float32)
 end