sample tree config (#25)

* samplg tree config

* hamming dist test

* fix doctest

Author: GiggleLiu

Project.toml

@@ -17,6 +17,7 @@ OMEinsum = "ebe7aa44-baf0-506c-a96f-8464559b3922"
 OMEinsumContractionOrders = "6f22d1fd-8eed-4bb7-9776-e7d684900715"
 Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 Primes = "27ebfcd6-29c5-5fa9-bf4b-fb8fc14df3ae"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 SIMDTypes = "94e857df-77ce-4151-89e5-788b33177be4"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
@@ -44,8 +45,7 @@ julia = "1"
 
 [extras]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "Documenter", "Random"]
+test = ["Test", "Documenter"]

docs/src/performancetips.md

@@ -93,6 +93,48 @@ Benchmark shows the performance of `TropicalGEMM` is close to the theoretical op
 It is a sum-product expression tree to store [`ConfigEnumerator`](@ref) in a lazy style, configurations can be extracted by depth first searching the tree with the `Base.collect` method. Although it is space efficient, it is in general not easy to extract information from it.
 This tree structure supports directed sampling so that one can get some statistic properties from it with an intermediate effort.
 
+For example, if we want to check some property of an intermediate scale graph, one can type
+```julia
+julia> graph = random_regular_graph(70, 3)
+
+julia> problem = IndependentSet(graph; optimizer=TreeSA());
+
+julia> tree = solve(problem, ConfigsAll(; tree_storage=true))[];
+16633909006371
+```
+If one wants to store these configurations, he will need a hard disk of size 256 TB!
+However, this sum-product binary tree structure supports efficient and unbiased direct sampling.
+
+```julia
+samples = generate_samples(tree, 1000);
+```
+
+With these samples, one can already compute useful properties like distribution of hamming distance (see [`hamming_distribution`](@ref)).
+
+```julia
+julia> using UnicodePlots
+
+julia> lineplot(hamming_distribution(samples, samples))
+          ┌────────────────────────────────────────┐ 
+   100000 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢰⠹⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
+          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡎⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
+          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⡇⠀⢣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
+          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
+          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸⠀⠀⠸⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
+          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
+          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⣇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
+          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢰⠃⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
+          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸⠀⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
+          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡞⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
+          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⣇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
+          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⢰⠁⠀⠀⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
+          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⡼⠀⠀⠀⠀⠀⠀⠈⡆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
+          │⠀⠀⠀⠀⠀⠀⠀⠀⢠⠇⠀⠀⠀⠀⠀⠀⠀⢳⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
+        0 │⢀⣀⣀⣀⣀⣀⣀⣀⠎⠀⠀⠀⠀⠀⠀⠀⠀⠀⠓⢄⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⠀⠀⠀⠀│ 
+          └────────────────────────────────────────┘ 
+          ⠀0⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀80⠀ 
+```
+
 (To be written.)
 
 ## Make use of GPUs

docs/src/ref.md

@@ -97,6 +97,9 @@ save_configs
 load_configs
 @bv_str
 onehotv
+
+generate_samples
+hamming_distribution
 ```
 
 ## Tensor Network

src/GraphTensorNetworks.jl

@@ -5,7 +5,7 @@ using Core: Argument
 using TropicalNumbers
 using OMEinsum
 using OMEinsum: timespace_complexity, getixsv
-using Graphs
+using Graphs, Random
 
 # OMEinsum
 export timespace_complexity, timespacereadwrite_complexity, @ein_str, getixsv, getiyv
@@ -19,6 +19,7 @@ export StaticBitVector, StaticElementVector, @bv_str
 export is_commutative_semiring
 export Max2Poly, TruncatedPoly, Polynomial, Tropical, CountingTropical, StaticElementVector, Mod, ConfigEnumerator, onehotv, ConfigSampler, TreeConfigEnumerator
 export CountingTropicalF64, CountingTropicalF32, TropicalF64, TropicalF32, ExtendedTropical
+export generate_samples
 
 # Lower level APIs
 export AllConfigs, SingleConfig
@@ -45,7 +46,7 @@ export is_dominating_set
 export solve, SizeMax, SizeMin, CountingAll, CountingMax, CountingMin, GraphPolynomial, SingleConfigMax, SingleConfigMin, ConfigsAll, ConfigsMax, ConfigsMin
 
 # Utilities
-export save_configs, load_configs
+export save_configs, load_configs, hamming_distribution
 
 # Visualization
 export show_graph, spring_layout

src/arithematics.jl

@@ -4,7 +4,7 @@ using Mods, Primes
 using Base.Cartesian
 import AbstractTrees: children, printnode, print_tree
 
-@enum TreeTag LEAF SUM PROD ZERO
+@enum TreeTag LEAF SUM PROD ZERO ONE
 
 # pirate
 Base.abs(x::Mod) = x
@@ -415,6 +415,7 @@ Base.one(::ConfigSampler{N,S,C}) where {N,S,C} = one(ConfigSampler{N,S,C})
 
 Configuration enumerator encoded in a tree, it is the most natural representation given by a sum-product network
 and is often more memory efficient than putting the configurations in a vector.
+One can use [`generate_samples`](@ref) to sample configurations from this tree structure efficiently.
 `N`, `S` and `C` are type parameters from the [`StaticElementVector`](@ref){N,S,C}.
 
 Fields
@@ -478,41 +479,43 @@ mutable struct TreeConfigEnumerator{N,S,C} <: AbstractSetNumber
         new{N,S,C}(LEAF, data)
     end
     function TreeConfigEnumerator{N,S,C}(tag::TreeTag) where {N,S,C}
-        @assert  tag === ZERO
+        @assert  tag === ZERO || tag === ONE
         return new{N,S,C}(tag)
     end
 end
 
 # AbstractTree APIs
 function children(t::TreeConfigEnumerator)
-    if t.tag == ZERO || t.tag == LEAF
+    if t.tag == ZERO || t.tag == LEAF || t.tag == ONE
         return typeof(t)[]
     else
         return [t.left, t.right]
     end
 end
-function printnode(io::IO, t::TreeConfigEnumerator)
+function printnode(io::IO, t::TreeConfigEnumerator{N,S,C}) where {N,S,C}
     if t.tag === LEAF
         print(io, t.data)
     elseif t.tag === ZERO
         print(io, "∅")
+    elseif t.tag === ONE
+        print(io, zero(StaticElementVector{N,S,C}))
     elseif t.tag === SUM
         print(io, "+")
     else  # PROD
         print(io, "*")
     end
 end
 
-Base.length(x::TreeConfigEnumerator) = _length(x, IdDict{typeof(x), Int}())
+Base.length(x::TreeConfigEnumerator) = _length!(x, IdDict{typeof(x), Int}())
 
-function _length(x, d)
+function _length!(x, d)
     haskey(d, x) && return d[x]
     if x.tag === SUM
-        l = _length(x.left, d) + _length(x.right, d)
+        l = _length!(x.left, d) + _length!(x.right, d)
         d[x] = l
         return l
     elseif x.tag === PROD
-        l = _length(x.left, d) * _length(x.right, d)
+        l = _length!(x.left, d) * _length!(x.right, d)
         d[x] = l
         return l
     elseif x.tag === ZERO
@@ -525,7 +528,7 @@ end
 num_nodes(x::TreeConfigEnumerator) = _num_nodes(x, IdDict{typeof(x), Int}())
 function _num_nodes(x, d)
     haskey(d, x) && return 0
-    if x.tag == ZERO
+    if x.tag == ZERO || x.tag == ONE
         res = 1
     elseif x.tag == LEAF
         res = 1
@@ -545,6 +548,8 @@ Base.show(io::IO, t::TreeConfigEnumerator) = print_tree(io, t)
 function Base.collect(x::TreeConfigEnumerator{N,S,C}) where {N,S,C}
     if x.tag == ZERO
         return StaticElementVector{N,S,C}[]
+    elseif x.tag == ONE
+        return StaticElementVector{N,S,C}[zero(StaticElementVector{N,S,C})]
     elseif x.tag == LEAF
         return StaticElementVector{N,S,C}[x.data]
     elseif x.tag == SUM
@@ -555,15 +560,33 @@ function Base.collect(x::TreeConfigEnumerator{N,S,C}) where {N,S,C}
 end
 
 function Base.:+(x::TreeConfigEnumerator{N,S,C}, y::TreeConfigEnumerator{N,S,C}) where {N,S,C}
-    TreeConfigEnumerator(SUM, x, y)
+    if x.tag == ZERO
+        return y
+    elseif y.tag == ZERO
+        return x
+    else
+        return TreeConfigEnumerator(SUM, x, y)
+    end
 end
 
 function Base.:*(x::TreeConfigEnumerator{L,S,C}, y::TreeConfigEnumerator{L,S,C}) where {L,S,C}
-    TreeConfigEnumerator(PROD, x, y)
+    if x.tag == ONE
+        return y
+    elseif y.tag == ONE
+        return x
+    elseif x.tag == ZERO
+        return x
+    elseif y.tag == ZERO
+        return y
+    elseif x.tag == LEAF && y.tag == LEAF
+        return TreeConfigEnumerator(x.data | y.data)
+    else
+        return TreeConfigEnumerator(PROD, x, y)
+    end
 end
 
 Base.zero(::Type{TreeConfigEnumerator{N,S,C}}) where {N,S,C} = TreeConfigEnumerator{N,S,C}(ZERO)
-Base.one(::Type{TreeConfigEnumerator{N,S,C}}) where {N,S,C} = TreeConfigEnumerator(zero(StaticElementVector{N,S,C}))
+Base.one(::Type{TreeConfigEnumerator{N,S,C}}) where {N,S,C} = TreeConfigEnumerator{N,S,C}(ONE)
 Base.zero(::TreeConfigEnumerator{N,S,C}) where {N,S,C} = zero(TreeConfigEnumerator{N,S,C})
 Base.one(::TreeConfigEnumerator{N,S,C}) where {N,S,C} = one(TreeConfigEnumerator{N,S,C})
 # todo, check siblings too?
@@ -572,13 +595,68 @@ function Base.iszero(t::TreeConfigEnumerator)
         iszero(t.left) && iszero(t.right)
     elseif t.tag == ZERO
         true
-    elseif t.tag == LEAF
+    elseif t.tag == LEAF || t.tag == ONE
         false
     else
         iszero(t.left) || iszero(t.right)
     end
 end
 
+"""
+    generate_samples(t::TreeConfigEnumerator, nsamples::Int)
+
+Direct sampling configurations from a [`TreeConfigEnumerator`](@ref) instance.
+
+Example
+-----------------------------
+```jldoctest; setup=:(using GraphTensorNetworks)
+julia> using Graphs
+
+julia> g= smallgraph(:petersen)
+{10, 15} undirected simple Int64 graph
+
+julia> t = solve(IndependentSet(g), ConfigsAll(; tree_storage=true))[];
+
+julia> samples = generate_samples(t, 1000);
+
+julia> all(s->is_independent_set(g, s), samples)
+true
+```
+"""
+function generate_samples(t::TreeConfigEnumerator{N,S,C}, nsamples::Int) where {N,S,C}
+    # get length dict
+    res = fill(zero(StaticElementVector{N,S,C}), nsamples)
+    d = IdDict{typeof(t), Int}()
+    sample_descend!(res, t, d)
+    return res
+end
+
+function sample_descend!(res::AbstractVector, t::TreeConfigEnumerator, d::IdDict)
+    length(res) == 0 && return res
+    if t.tag == LEAF
+        res .|= Ref(t.data)
+    elseif t.tag == SUM
+        ratio = _length!(t.left, d)/_length!(t, d)
+        nleft = 0
+        for _ = 1:length(res)
+            if rand() < ratio
+                nleft += 1
+            end
+        end
+        shuffle!(res)  # shuffle the `res` to avoid biased sampling, very important.
+        sample_descend!(view(res,1:nleft), t.left, d)
+        sample_descend!(view(res,nleft+1:length(res)), t.right, d)
+    elseif t.tag == PROD
+        sample_descend!(res, t.right, d)
+        sample_descend!(res, t.left, d)
+    elseif t.tag == ZERO
+        error("Meet zero when descending.")
+    else
+        # pass for 1
+    end
+    return res
+end
+
 # A patch to make `Polynomial{ConfigEnumerator}` work
 function Base.:*(a::Int, y::AbstractSetNumber)
     a == 0 && return zero(y)
@@ -611,7 +689,8 @@ end
 
 # utilities for creating onehot vectors
 onehotv(::Type{ConfigEnumerator{N,S,C}}, i::Integer, v) where {N,S,C} = ConfigEnumerator([onehotv(StaticElementVector{N,S,C}, i, v)])
-onehotv(::Type{TreeConfigEnumerator{N,S,C}}, i::Integer, v) where {N,S,C} = TreeConfigEnumerator(onehotv(StaticElementVector{N,S,C}, i, v))
+# we treat `v == 0` specially because we want the final result not containing one leaves.
+onehotv(::Type{TreeConfigEnumerator{N,S,C}}, i::Integer, v) where {N,S,C} = v == 0 ? one(TreeConfigEnumerator{N,S,C}) : TreeConfigEnumerator(onehotv(StaticElementVector{N,S,C}, i, v))
 onehotv(::Type{ConfigSampler{N,S,C}}, i::Integer, v) where {N,S,C} = ConfigSampler(onehotv(StaticElementVector{N,S,C}, i, v))
 # just to make matrix transpose work
 Base.transpose(c::ConfigEnumerator) = c
@@ -620,7 +699,7 @@ Base.transpose(c::TreeConfigEnumerator) = c
 function Base.copy(c::TreeConfigEnumerator{N,S,C}) where {N,S,C}
     if c.tag == LEAF
         TreeConfigEnumerator(c.data)
-    elseif c.tag == ZERO
+    elseif c.tag == ZERO || c.tag == ONE
         TreeConfigEnumerator{N,S,C}(c.tag)
     else
         TreeConfigEnumerator(c.tag, c.left, c.right)
@@ -635,9 +714,9 @@ end
 
 # to handle power of polynomials
 function Base.:^(x::TreeConfigEnumerator, y::Real)
-    if y <= 0
+    if y == 0
         return one(x)
-    elseif x.tag == LEAF
+    elseif x.tag == LEAF || x.tag == ONE || x.tag == ZERO
         return x
     else
         error("pow of non-leaf nodes is forbidden!")
@@ -688,6 +767,19 @@ function _x(::Type{T}; invert) where {T<:AbstractSetNumber}
     invert ? pre_invert_exponent(ret) : ret
 end
 
+function _onehotv(::Type{Polynomial{BS,X}}, x, v) where {BS,X}
+    Polynomial{BS,X}([onehotv(BS, x, v)])
+end
+function _onehotv(::Type{TruncatedPoly{K,BS,OS}}, x, v) where {K,BS,OS}
+    TruncatedPoly{K,BS,OS}(ntuple(i->i != K ? zero(BS) : onehotv(BS, x, v), K),zero(OS))
+end
+function _onehotv(::Type{CountingTropical{TV,BS}}, x, v) where {TV,BS}
+    CountingTropical{TV,BS}(zero(TV), onehotv(BS, x, v))
+end
+function _onehotv(::Type{BS}, x, v) where {BS<:AbstractSetNumber}
+    onehotv(BS, x, v)
+end
+
 # negate the exponents before entering the solver
 pre_invert_exponent(t::TruncatedPoly{K}) where K = TruncatedPoly(t.coeffs, -t.maxorder)
 pre_invert_exponent(t::TropicalNumbers.TropicalTypes) = inv(t)

src/configurations.jl

@@ -90,15 +90,30 @@ e.g. when the problem is [`MaximalIS`](@ref), it computes all maximal independen
 """
 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)])
-end
-function _onehotv(::Type{TruncatedPoly{K,BS,OS}}, x, v) where {K,BS,OS}
-    TruncatedPoly{K,BS,OS}(ntuple(i->i != K ? zero(BS) : onehotv(BS, x, v), K),zero(OS))
+# NOTE: do we have more efficient way to compute it?
+# NOTE: doing pair-wise Hamming distance might be biased?
+"""
+    hamming_distribution(S, T)
+
+Compute the distribution of pair-wise Hamming distances, which is defined as:
+```math
+c(k) := \\sum_{\\sigma\\in S, \\tau\\in T} \\delta({\\rm dist}(\\sigma, \\tau), k)
+```
+where ``\\delta`` is a function that returns 1 if two arguments are equivalent, 0 otherwise,
+``{\\rm dist}`` is the Hamming distance function.
+
+Returns the counting as a vector.
+"""
+function hamming_distribution(t1::ConfigEnumerator, t2::ConfigEnumerator)
+    return hamming_distribution(t1.data, t2.data)
 end
-function _onehotv(::Type{CountingTropical{TV,BS}}, x, v) where {TV,BS}
-    CountingTropical{TV,BS}(zero(TV), onehotv(BS, x, v))
+function hamming_distribution(s1::AbstractVector{StaticElementVector{N,S,C}}, s2::AbstractVector{StaticElementVector{N,S,C}}) where {N,S,C}
+    return hamming_distribution!(zeros(Int, N+1), s1, s2)
 end
-function _onehotv(::Type{BS}, x, v) where {BS<:AbstractSetNumber}
-    onehotv(BS, x, v)
+function hamming_distribution!(out::AbstractVector, s1::AbstractVector{StaticElementVector{N,S,C}}, s2::AbstractVector{StaticElementVector{N,S,C}}) where {N,S,C}
+    @assert length(out) == N+1
+    @inbounds for a in s1, b in s2
+        out[count_ones(a ⊻ b)+1] += 1
+    end
+    return out
 end