From 857970acd4ba84d769d291781b3b657fcf993fbf Mon Sep 17 00:00:00 2001 From: GiggleLiu Date: Mon, 14 Mar 2022 16:20:26 -0400 Subject: [PATCH 1/4] save --- examples/IndependentSet.jl | 2 +- examples/spectrumcomputing.jl | 47 ++++++++++++++++++ src/arithematics.jl | 90 +++++++++++++++++++++++++++++++++++ src/bounding.jl | 2 +- test/arithematics.jl | 11 +++++ 5 files changed, 150 insertions(+), 2 deletions(-) create mode 100644 examples/spectrumcomputing.jl diff --git a/examples/IndependentSet.jl b/examples/IndependentSet.jl index 65a5dcf9..3d488553 100644 --- a/examples/IndependentSet.jl +++ b/examples/IndependentSet.jl @@ -32,7 +32,7 @@ show_graph(graph; locs=locations) # x_i^{w_i} # \end{matrix}\right), # ``` -# where ``W(x_i)_0=1`` is the first element associated with ``s_i=0`` and ``W(x_i)_1=x_i^{w_i}`` is the second element associated with ``s_i=1``, and `w_i` is the weight of vertex ``i``. +# where ``W(x_i)_0=1`` is the first element associated with ``s_i=0`` and ``W(x_i)_1=x_i^{w_i}`` is the second element associated with ``s_i=1``, and ``w_i`` is the weight of vertex ``i``. # Similarly, on each edge ``(u, v)``, we define a matrix ``B`` indexed by ``s_u`` and ``s_v`` as # ```math # B = \left(\begin{matrix} diff --git a/examples/spectrumcomputing.jl b/examples/spectrumcomputing.jl new file mode 100644 index 00000000..87ed64d5 --- /dev/null +++ b/examples/spectrumcomputing.jl @@ -0,0 +1,47 @@ +using GraphTensorNetworks, Graphs + +function k23() + g = SimpleGraph(5) + for (i, j) in [(1,5), (1,4), (1,2), (3,2), (3,4), (3, 5)] + add_edge!(g, i, j) + end + return g +end + +function mapped_k23() + g = SimpleGraph(9) + for (i, j) in [(1,5), (1,4), (1,7), (2,8), (2,9), + (3, 5), (3, 4), (3, 6), (4,8), (4, 9), (6, 8), + (7,9), + ] + add_edge!(g, i, j) + end + return g +end + +function compute_spectrums(problem, nlevel, weights_list) + res = zeros(nlevel, length(weights_list)) + for i=1:length(weights_list) + weights = copy(problem.weights) + weights[1:length(weights_list[i])] .+= weights_list[i] + prob = IndependentSet(problem.code, length(weights), weights) + res[:,i] .= solve(prob, SizeMax(nlevel))[].orders + end + return res +end + +source = Independence(k23(), weights=zeros(5)) +mapped_weight0 = [1.0, 0, 1, 0, 0, 2, 2, 1, 1] +mapped = Independence(mapped_k23(), weights=mapped_weight0) +nlevel = Int(solve(source, CountingAll())[]) + +weights_list = [[0.01+0.01*sin(2π*t/(1<<(k-1))) for k=1:5] for t=0.0:0.1:16.0] +spectrum1 = compute_spectrums(source, nlevel, weights_list) +spectrum2 = compute_spectrums(mapped, nlevel, weights_list) + +using PyPlot +PyPlot.subplot(211) +PyPlot.plot(spectrum1') +PyPlot.subplot(212) +PyPlot.plot(spectrum2') +PyPlot.xlabel("weights") \ No newline at end of file diff --git a/src/arithematics.jl b/src/arithematics.jl index 35b7d25f..4ae8851d 100644 --- a/src/arithematics.jl +++ b/src/arithematics.jl @@ -222,6 +222,96 @@ function Base.:*(a::ExtendedTropical{K,TO}, b::ExtendedTropical{K,TO}) where {K, end function sorted_sum_combination!(res::AbstractVector{TO}, A::AbstractVector{TO}, B::AbstractVector{TO}) where TO + K = length(res) + @assert length(B) == length(A) == K + @inbounds high = A[K] * B[K] + + mA = findfirst(!iszero, A) + mB = findfirst(!iszero, B) + if mA === nothing || mB === nothing + res .= Ref(zero(TO)) + return res + end + @inbounds low = A[mA] * B[mB] + # count number bigger than x + c = count_geq(A, B, mB, low) + @inbounds if c <= K # return + res[K-c+1:K] .= sort!(collect_geq!(view(res,1:c), A, B, mB, low)) + if c < K + res[1:K-c] .= zero(TO) + end + return res + end + # calculate by bisection for at most 30 times. + for _ = 1:30 + mid = mid_point(high, low) + c = count_geq(A, B, mB, mid) + if c > K + low = mid + elseif c == K # return + return sort!(collect_geq!(res, A, B, mB, mid)) + else + high = mid + end + end + clow = count_geq(A, B, mB, low) + res .= sort!(collect_geq!(similar(res, clow), A, B, mB, low))[end-K+1:end] + return res +end + +function count_geq(A, B, mB, low) + K = length(A) + k = 1 # TODO: we should tighten mA, mB later! + @inbounds Ak = A[K-k+1] + @inbounds Bq = B[K-mB+1] + c = 0 + @inbounds for q = K-mB+1:-1:1 + Bq = B[K-q+1] + while k < K && Ak * Bq >= low + k += 1 + Ak = A[K-k+1] + end + if Ak * Bq >= low + c += k + else + c += (k-1) + end + #if c > K + #return c + #end + end + return c +end + +function collect_geq!(res, A, B, mB, low) + K = length(A) + k = 1 # TODO: we should tighten mA, mB later! + @inbounds Ak = A[K-k+1] + @inbounds Bq = B[K-mB+1] + l = 0 + @inbounds for q = K-mB+1:-1:1 + Bq = B[K-q+1] + while k < K && Ak * Bq >= low + k += 1 + Ak = A[K-k+1] + end + # push data + ck = Ak * Bq >= low ? k : k-1 + for j=1:ck + l += 1 + res[l] = Bq * A[end-j+1] + end + end + return res +end + +# for bisection +mid_point(a::Tropical{T}, b::Tropical{T}) where T = Tropical{T}((a.n + b.n) / 2) +mid_point(a::CountingTropical{T,CT}, b::CountingTropical{T,CT}) where {T,CT} = CountingTropical{T,CT}((a.n + b.n) / 2, a.c) +mid_point(a::Tropical{T}, b::Tropical{T}) where T<:Integer = Tropical{T}((a.n + b.n) ÷ 2) +mid_point(a::CountingTropical{T,CT}, b::CountingTropical{T,CT}) where {T<:Integer,CT} = CountingTropical{T,CT}((a.n + b.n) ÷ 2, a.c) + +function sorted_sum_combination1!(res::AbstractVector{TO}, A::AbstractVector{TO}, B::AbstractVector{TO}) where TO K = length(res) @assert length(B) == length(A) == K @inbounds maxval = A[K] * B[K] diff --git a/src/bounding.jl b/src/bounding.jl index 99426924..be66b553 100644 --- a/src/bounding.jl +++ b/src/bounding.jl @@ -60,7 +60,7 @@ struct CacheTree{T} end function cached_einsum(se::SlicedEinsum, @nospecialize(xs), size_dict) if length(se.slicing) != 0 - @warn "Slicing is not supported for caching! Fallback to `NestedEinsum`." + @warn "Slicing is not supported for caching, got nslices = $(length(se.slicing))! Fallback to `NestedEinsum`." end return cached_einsum(se.eins, xs, size_dict) end diff --git a/test/arithematics.jl b/test/arithematics.jl index b7059fee..5b978c43 100644 --- a/test/arithematics.jl +++ b/test/arithematics.jl @@ -153,6 +153,17 @@ end end end +@testset "count geq" begin + A = collect(1:10) + B = collect(2:2:20) + low = 20 + c = GraphTensorNetworks.count_geq(A, B, 1, low) + @test c == count(x->x>=low, vec([a*b for a in A, b in B])) + res = similar(A, c) + @test sort!(GraphTensorNetworks.collect_geq!(res, A, B, 1, low)) == sort!(filter(x->x>=low, vec([a*b for a in A, b in B]))) +end + + # check the correctness of sampling @testset "generate samples" begin Random.seed!(2) From fbd3781a5ad5d21c42d46d0e5a2e637750fc2462 Mon Sep 17 00:00:00 2001 From: GiggleLiu Date: Mon, 14 Mar 2022 18:29:44 -0400 Subject: [PATCH 2/4] update --- src/arithematics.jl | 68 +++++++++++++------------------------------- test/arithematics.jl | 2 +- 2 files changed, 20 insertions(+), 50 deletions(-) diff --git a/src/arithematics.jl b/src/arithematics.jl index 4ae8851d..6438e20c 100644 --- a/src/arithematics.jl +++ b/src/arithematics.jl @@ -234,7 +234,7 @@ function sorted_sum_combination!(res::AbstractVector{TO}, A::AbstractVector{TO}, end @inbounds low = A[mA] * B[mB] # count number bigger than x - c = count_geq(A, B, mB, low) + c, _ = count_geq(A, B, mB, low, true) @inbounds if c <= K # return res[K-c+1:K] .= sort!(collect_geq!(view(res,1:c), A, B, mB, low)) if c < K @@ -243,28 +243,31 @@ function sorted_sum_combination!(res::AbstractVector{TO}, A::AbstractVector{TO}, return res end # calculate by bisection for at most 30 times. - for _ = 1:30 + @inbounds for _ = 1:30 mid = mid_point(high, low) - c = count_geq(A, B, mB, mid) + c, nB = count_geq(A, B, mB, mid, true) if c > K low = mid + mB = nB elseif c == K # return + # NOTE: this is the bottleneck return sort!(collect_geq!(res, A, B, mB, mid)) else high = mid end end - clow = count_geq(A, B, mB, low) - res .= sort!(collect_geq!(similar(res, clow), A, B, mB, low))[end-K+1:end] + clow, _ = count_geq(A, B, mB, low, false) + @inbounds res .= sort!(collect_geq!(similar(res, clow), A, B, mB, low))[end-K+1:end] return res end -function count_geq(A, B, mB, low) +function count_geq(A, B, mB, low, earlybreak) K = length(A) k = 1 # TODO: we should tighten mA, mB later! @inbounds Ak = A[K-k+1] @inbounds Bq = B[K-mB+1] c = 0 + nB = mB @inbounds for q = K-mB+1:-1:1 Bq = B[K-q+1] while k < K && Ak * Bq >= low @@ -275,21 +278,24 @@ function count_geq(A, B, mB, low) c += k else c += (k-1) + if k==1 + nB += 1 + end + end + if earlybreak && c > K + return c, nB end - #if c > K - #return c - #end end - return c + return c, nB end function collect_geq!(res, A, B, mB, low) K = length(A) k = 1 # TODO: we should tighten mA, mB later! - @inbounds Ak = A[K-k+1] - @inbounds Bq = B[K-mB+1] + Ak = A[K-k+1] + Bq = B[K-mB+1] l = 0 - @inbounds for q = K-mB+1:-1:1 + for q = K-mB+1:-1:1 Bq = B[K-q+1] while k < K && Ak * Bq >= low k += 1 @@ -311,42 +317,6 @@ mid_point(a::CountingTropical{T,CT}, b::CountingTropical{T,CT}) where {T,CT} = C mid_point(a::Tropical{T}, b::Tropical{T}) where T<:Integer = Tropical{T}((a.n + b.n) ÷ 2) mid_point(a::CountingTropical{T,CT}, b::CountingTropical{T,CT}) where {T<:Integer,CT} = CountingTropical{T,CT}((a.n + b.n) ÷ 2, a.c) -function sorted_sum_combination1!(res::AbstractVector{TO}, A::AbstractVector{TO}, B::AbstractVector{TO}) where TO - K = length(res) - @assert length(B) == length(A) == K - @inbounds maxval = A[K] * B[K] - ptr = K - @inbounds res[ptr] = maxval - @inbounds queue = [(K,K-1,A[K]*B[K-1]), (K-1,K,A[K-1]*B[K])] - for k = 1:K-1 - @inbounds (i, j, res[K-k]) = _pop_max_sum!(queue) # TODO: do not enumerate, use better data structures - _push_if_not_exists!(queue, i, j-1, A, B) - _push_if_not_exists!(queue, i-1, j, A, B) - end - return res -end - -function _push_if_not_exists!(queue, i, j, A, B) - @inbounds if j>=1 && i>=1 && !any(x->x[1] >= i && x[2] >= j, queue) - push!(queue, (i, j, A[i]*B[j])) - end -end - -function _pop_max_sum!(queue) - maxsum = first(queue)[3] - maxloc = 1 - @inbounds for i=2:length(queue) - m = queue[i][3] - if m > maxsum - maxsum = m - maxloc = i - end - end - @inbounds data = queue[maxloc] - deleteat!(queue, maxloc) - return data -end - function Base.:+(a::ExtendedTropical{K,TO}, b::ExtendedTropical{K,TO}) where {K,TO} res = Vector{TO}(undef, K) ptr1, ptr2 = K, K diff --git a/test/arithematics.jl b/test/arithematics.jl index 5b978c43..036f413f 100644 --- a/test/arithematics.jl +++ b/test/arithematics.jl @@ -157,7 +157,7 @@ end A = collect(1:10) B = collect(2:2:20) low = 20 - c = GraphTensorNetworks.count_geq(A, B, 1, low) + c, _ = GraphTensorNetworks.count_geq(A, B, 1, low, false) @test c == count(x->x>=low, vec([a*b for a in A, b in B])) res = similar(A, c) @test sort!(GraphTensorNetworks.collect_geq!(res, A, B, 1, low)) == sort!(filter(x->x>=low, vec([a*b for a in A, b in B]))) From fc072253264e3c0dce96709613e5de9de4abcf03 Mon Sep 17 00:00:00 2001 From: GiggleLiu Date: Mon, 14 Mar 2022 18:30:09 -0400 Subject: [PATCH 3/4] clean up --- examples/spectrumcomputing.jl | 47 ----------------------------------- 1 file changed, 47 deletions(-) delete mode 100644 examples/spectrumcomputing.jl diff --git a/examples/spectrumcomputing.jl b/examples/spectrumcomputing.jl deleted file mode 100644 index 87ed64d5..00000000 --- a/examples/spectrumcomputing.jl +++ /dev/null @@ -1,47 +0,0 @@ -using GraphTensorNetworks, Graphs - -function k23() - g = SimpleGraph(5) - for (i, j) in [(1,5), (1,4), (1,2), (3,2), (3,4), (3, 5)] - add_edge!(g, i, j) - end - return g -end - -function mapped_k23() - g = SimpleGraph(9) - for (i, j) in [(1,5), (1,4), (1,7), (2,8), (2,9), - (3, 5), (3, 4), (3, 6), (4,8), (4, 9), (6, 8), - (7,9), - ] - add_edge!(g, i, j) - end - return g -end - -function compute_spectrums(problem, nlevel, weights_list) - res = zeros(nlevel, length(weights_list)) - for i=1:length(weights_list) - weights = copy(problem.weights) - weights[1:length(weights_list[i])] .+= weights_list[i] - prob = IndependentSet(problem.code, length(weights), weights) - res[:,i] .= solve(prob, SizeMax(nlevel))[].orders - end - return res -end - -source = Independence(k23(), weights=zeros(5)) -mapped_weight0 = [1.0, 0, 1, 0, 0, 2, 2, 1, 1] -mapped = Independence(mapped_k23(), weights=mapped_weight0) -nlevel = Int(solve(source, CountingAll())[]) - -weights_list = [[0.01+0.01*sin(2π*t/(1<<(k-1))) for k=1:5] for t=0.0:0.1:16.0] -spectrum1 = compute_spectrums(source, nlevel, weights_list) -spectrum2 = compute_spectrums(mapped, nlevel, weights_list) - -using PyPlot -PyPlot.subplot(211) -PyPlot.plot(spectrum1') -PyPlot.subplot(212) -PyPlot.plot(spectrum2') -PyPlot.xlabel("weights") \ No newline at end of file From 6fb0c35ff38d3ce304730a7ecc0c8e5e86f01e93 Mon Sep 17 00:00:00 2001 From: GiggleLiu Date: Mon, 14 Mar 2022 18:51:16 -0400 Subject: [PATCH 4/4] add one more comment --- src/arithematics.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/arithematics.jl b/src/arithematics.jl index 6438e20c..15924622 100644 --- a/src/arithematics.jl +++ b/src/arithematics.jl @@ -221,6 +221,9 @@ function Base.:*(a::ExtendedTropical{K,TO}, b::ExtendedTropical{K,TO}) where {K, return ExtendedTropical{K,TO}(sorted_sum_combination!(res, a.orders, b.orders)) end +# 1. bisect over summed value and find the critical value `c`, +# 2. collect the values with sum combination `≥ c`, +# 3. sort the collected values function sorted_sum_combination!(res::AbstractVector{TO}, A::AbstractVector{TO}, B::AbstractVector{TO}) where TO K = length(res) @assert length(B) == length(A) == K @@ -261,6 +264,7 @@ function sorted_sum_combination!(res::AbstractVector{TO}, A::AbstractVector{TO}, return res end +# count the number of sum-combinations with the sum >= low function count_geq(A, B, mB, low, earlybreak) K = length(A) k = 1 # TODO: we should tighten mA, mB later!