diff --git a/docs/make.jl b/docs/make.jl index 3e21363a..b399192a 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -35,7 +35,6 @@ makedocs(; "Max-Cut Problem" => "tutorials/MaxCut.md", "Other Problems" => "tutorials/Others.md", ], - "Method Selection Guide" => "methodselection.md", "References" => "ref.md", ], doctest=false, diff --git a/docs/src/methodselection.md b/docs/src/methodselection.md deleted file mode 100644 index a04b48dd..00000000 --- a/docs/src/methodselection.md +++ /dev/null @@ -1,11 +0,0 @@ -```@meta -CurrentModule = GraphTensorNetworks -``` - -# Method Selection Guide - -## Tensor network contraction tree optimizer - -## Independence polynomial - -## Configuration enumeration diff --git a/docs/src/ref.md b/docs/src/ref.md index a0df597b..1814b618 100644 --- a/docs/src/ref.md +++ b/docs/src/ref.md @@ -2,6 +2,7 @@ ## Graph problems ```@docs solve +GraphProblem Independence MaximalIndependence Matching @@ -14,6 +15,20 @@ PaintShop set_packing ``` +#### Graph Problem Interfaces +```@docs +generate_tensors +symbols +flavors +get_weights +nflavor +``` + +To subtype [`GraphProblem`](@ref), the new type must contain a `code` field to represent the (optimized) tensor network. +Interfaces [`generate_tensors`](@ref), [`symbols`](@ref), [`flavors`](@ref) and [`get_weights`] are required. +[`nflavor`] is optimal. + + ## Properties ```@docs SizeMax @@ -50,7 +65,10 @@ is_commutative_semiring ## Tensor Network ```@docs optimize_code +getixsv +getiyv timespace_complexity +timespacereadwrite_complexity @ein_str GreedyMethod TreeSA diff --git a/src/GraphTensorNetworks.jl b/src/GraphTensorNetworks.jl index 94bc6c5d..64bdfdbe 100644 --- a/src/GraphTensorNetworks.jl +++ b/src/GraphTensorNetworks.jl @@ -8,7 +8,7 @@ using OMEinsum: timespace_complexity, getixsv using Graphs # OMEinsum -export timespace_complexity, @ein_str +export timespace_complexity, timespacereadwrite_complexity, @ein_str, getixsv, getiyv export GreedyMethod, TreeSA, SABipartite, KaHyParBipartite, MergeVectors, MergeGreedy # Algebras @@ -27,7 +27,8 @@ export random_regular_graph, diagonal_coupled_graph, is_independent_set export square_lattice_graph, unitdisk_graph # Tensor Networks (Graph problems) -export Independence, MaximalIndependence, Matching, Coloring, optimize_code, set_packing, MaxCut, PaintShop, paintshop_from_pairs, UnWeighted +export GraphProblem, Independence, MaximalIndependence, Matching, Coloring, optimize_code, set_packing, MaxCut, PaintShop, paintshop_from_pairs, UnWeighted +export flavors, symbols, nflavor, get_weights # Interfaces export solve, SizeMax, CountingAll, CountingMax, GraphPolynomial, SingleConfigMax, ConfigsAll, ConfigsMax @@ -44,7 +45,7 @@ project_relative_path(xs...) = normpath(joinpath(dirname(dirname(pathof(@__MODUL include("utils.jl") include("bitvector.jl") include("arithematics.jl") -include("networks.jl") +include("networks/networks.jl") include("graph_polynomials.jl") include("configurations.jl") include("graphs.jl") diff --git a/src/configurations.jl b/src/configurations.jl index d0b9a3fb..886bb459 100644 --- a/src/configurations.jl +++ b/src/configurations.jl @@ -11,16 +11,16 @@ function best_solutions(gp::GraphProblem; all=false, usecuda=false) throw(ArgumentError("ConfigEnumerator can not be computed on GPU!")) end syms = symbols(gp) - T = (all ? set_type : sampler_type)(CountingTropical{Int64}, length(syms), bondsize(gp)) + T = (all ? set_type : sampler_type)(CountingTropical{Int64}, length(syms), nflavor(gp)) vertex_index = Dict([s=>i for (i, s) in enumerate(syms)]) - xst = generate_tensors(l->TropicalF64(get_weight(gp, vertex_index[l])), gp) + xst = generate_tensors(l->TropicalF64.(get_weights(gp, l)), gp) ymask = trues(fill(2, length(getiyv(gp.code)))...) if usecuda xst = CuArray.(xst) ymask = CuArray(ymask) end if all - xs = generate_tensors(l->_onehotv(T, vertex_index[l], 1, get_weight(gp, vertex_index[l])), gp) + xs = generate_tensors(l->_onehotv.(Ref(T), vertex_index[l], flavors(gp), get_weights(gp, l)), gp) return bounding_contract(AllConfigs{1}(), gp.code, xst, ymask, xs) else @assert ndims(ymask) == 0 @@ -59,57 +59,39 @@ best2_solutions(gp::GraphProblem; all=true, usecuda=false) = solutions(gp, Max2P function bestk_solutions(gp::GraphProblem, k::Int) syms = symbols(gp) vertex_index = Dict([s=>i for (i, s) in enumerate(syms)]) - xst = generate_tensors(l->TropicalF64(get_weight(gp, vertex_index[l])), gp) + xst = generate_tensors(l->TropicalF64.(get_weights(gp, l)), gp) ymask = trues(fill(2, length(getiyv(gp.code)))...) - T = set_type(TruncatedPoly{k,Float64,Float64}, length(syms), bondsize(gp)) - xs = generate_tensors(l->_onehotv(T, vertex_index[l], 1, get_weight(gp, vertex_index[l])), gp) + T = set_type(TruncatedPoly{k,Float64,Float64}, length(syms), nflavor(gp)) + xs = generate_tensors(l->_onehotv.(Ref(T), vertex_index[l], flavors(gp), get_weights(gp, l)), gp) return bounding_contract(AllConfigs{k}(), gp.code, xst, ymask, xs) end """ all_solutions(problem) -Finding all solutions. +Finding all solutions grouped by size. e.g. when the problem is `MaximalIndependence`, 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) -# return a mapping from label to variable `x` -for GP in [:Independence, :Matching, :MaximalIndependence, :MaxCut, :PaintShop] - @eval function fx_solutions(gp::$GP, ::Type{BT}, all::Bool) where BT - syms = symbols(gp) - T = (all ? set_type : sampler_type)(BT, length(syms), bondsize(gp)) - vertex_index = Dict([s=>i for (i, s) in enumerate(syms)]) - return l->_onehotv(T, vertex_index[l], 1, get_weight(gp, vertex_index[l])) - end -end -function fx_solutions(gp::Coloring{K}, ::Type{BT}, all::Bool) where {K,BT} +# return a mapping from label to onehot bitstrings (degree of freedoms). +function fx_solutions(gp::GraphProblem, ::Type{BT}, all::Bool) where {K,BT} syms = symbols(gp) - T = (all ? set_type : sampler_type)(BT, length(syms), bondsize(gp)) + T = (all ? set_type : sampler_type)(BT, length(syms), nflavor(gp)) vertex_index = Dict([s=>i for (i, s) in enumerate(syms)]) return function (l) - map(1:K) do k - _onehotv(T, vertex_index[l], k, get_weight(gp, vertex_index[l])) - end + _onehotv.(Ref(T), vertex_index[l], flavors(gp), get_weights(gp, l)) end end function _onehotv(::Type{Polynomial{BS,X}}, x, v, w) where {BS,X} - @assert isone(w) - Polynomial{BS,X}([zero(BS), onehotv(BS, x, v)]) + Polynomial{BS,X}([zero(BS), onehotv(BS, x, v)])^w end function _onehotv(::Type{TruncatedPoly{K,BS,OS}}, x, v, w) where {K,BS,OS} - @assert isone(w) - TruncatedPoly{K,BS,OS}(ntuple(i->iix, gp; usecuda=usecuda) -function contractf(f, gp::GraphProblem; usecuda=false) - @debug "generating tensors ..." - xs = generate_tensors(f, gp) - @debug "contracting tensors ..." - if usecuda - gp.code([CuArray(x) for x in xs]...) - else - gp.code(xs...) - end -end - -############### Problem specific implementations ################ -### independent set ### -function generate_tensors(fx, gp::Independence) - ixs = getixsv(gp.code) - n = length(unique!(vcat(ixs...))) - T = typeof(fx(ixs[1][1])) - return map(enumerate(ixs)) do (i, ix) - if i <= n - misv(fx(ix[1])) - else - misb(T, length(ix)) # if n!=2, it corresponds to set packing problem. - end - end -end - -function misb(::Type{T}, n::Integer=2) where T - res = zeros(T, fill(2, n)...) - res[1] = one(T) - for i=1:n - res[1+1<<(i-1)] = one(T) - end - return res -end -misv(val::T) where T = [one(T), val] - -### coloring ### -function generate_tensors(fx, c::Coloring{K}) where K - ixs = getixsv(c.code) - T = eltype(fx(ixs[1][1])) - return map(ixs) do ix - # if the tensor rank is 1, create a vertex tensor. - # otherwise the tensor rank must be 2, create a bond tensor. - length(ix)==1 ? coloringv(fx(ix[1])) : coloringb(T, K) - end -end - -# coloring bond tensor -function coloringb(::Type{T}, k::Int) where T - x = ones(T, k, k) - for i=1:k - x[i,i] = zero(T) - end - return x -end -# coloring vertex tensor -coloringv(vals::Vector{T}) where T = vals - -### matching ### -function generate_tensors(fx, m::Matching) - ixs = getixsv(m.code) - T = typeof(fx(ixs[1][1])) - n = length(unique!(vcat(ixs...))) # number of vertices - tensors = [] - for i=1:length(ixs) - if i<=n - @assert length(ixs[i]) == 1 - t = T[one(T), fx(ixs[i][1])] # fx is defined on edges. - else - t = match_tensor(T, length(ixs[i])) - end - push!(tensors, t) - end - return tensors -end -function match_tensor(::Type{T}, n::Int) where T - t = zeros(T, fill(2, n)...) - for ci in CartesianIndices(t) - if sum(ci.I .- 1) <= 1 - t[ci] = one(T) - end - end - return t -end - -### maximal independent set ### -function generate_tensors(fx, mi::MaximalIndependence) - ixs = getixsv(mi.code) - T = eltype(fx(ixs[1][end])) - return map(ixs) do ix - neighbortensor(fx(ix[end]), length(ix)) - end -end -function neighbortensor(x::T, d::Int) where T - t = zeros(T, fill(2, d)...) - for i = 2:1<<(d-1) - t[i] = one(T) - end - t[1<<(d-1)+1] = x - return t -end - -### max cut/spin glass problem ### -function generate_tensors(fx, gp::MaxCut) - ixs = getixsv(gp.code) - return map(enumerate(ixs)) do (i, ix) - maxcutb(fx(ix)) - end -end -function maxcutb(expJ::T) where T - return T[one(T) expJ; expJ one(T)] -end - -### paint shop ### -function generate_tensors(fx, c::PaintShop) - ixs = getixsv(c.code) - T = eltype(fx(ixs[1])) - [paintshop_bond_tensor(fx(ixs[i]), c.isfirst[i], c.isfirst[i+1]) for i=1:length(ixs)] -end - -function paintshop_bond_tensor(x::T, if1::Bool, if2::Bool) where T - m = T[x one(T); one(T) x] - m = if1 ? m : m[[2,1],:] - m = if2 ? m : m[:,[2,1]] - return m +for GP in [:Independence, :Matching, :MaxCut, :MaximalIndependence, :PaintShop] + @eval contractx(gp::$GP, x::T; usecuda=false) where T = contractf(l->Ref(x) .^ get_weights(gp, l), gp; usecuda=usecuda) end \ No newline at end of file diff --git a/src/interfaces.jl b/src/interfaces.jl index 95db966a..c527e891 100644 --- a/src/interfaces.jl +++ b/src/interfaces.jl @@ -1,5 +1,4 @@ abstract type AbstractProperty end -_support_weight(::AbstractProperty) = false """ SizeMax <: AbstractProperty @@ -12,7 +11,6 @@ The maximum independent set size. * BLAS (on CPU) and GPU are supported, """ struct SizeMax <: AbstractProperty end -_support_weight(::SizeMax) = true """ CountingAll <: AbstractProperty @@ -25,7 +23,6 @@ Counting the total number of sets. e.g. for [`Independence`](@ref) problem, it c * BLAS (GPU and CPU) and GPU are supported, """ struct CountingAll <: AbstractProperty end -_support_weight(::CountingAll) = true """ CountingMax{K} <: AbstractProperty @@ -41,7 +38,6 @@ it counts independent sets of size ``\\alpha(G), \\alpha(G)-1, \\ldots, \\alpha( struct CountingMax{K} <: AbstractProperty end CountingMax(K::Int=1) = CountingMax{K}() max_k(::CountingMax{K}) where K = K -_support_weight(::CountingMax{1}) = true """ GraphPolynomial{METHOD} <: AbstractProperty @@ -87,7 +83,6 @@ Finding single best solution, e.g. for [`Independence`](@ref) problem, it is one """ struct SingleConfigMax{BOUNDED} <:AbstractProperty end SingleConfigMax(; bounded::Bool=false) = SingleConfigMax{bounded}() -_support_weight(::SingleConfigMax) = true """ ConfigsAll <:AbstractProperty @@ -99,7 +94,6 @@ Find all valid configurations, e.g. for [`Independence`](@ref) problem, it is fi * Weights do not take effect. """ struct ConfigsAll <:AbstractProperty end -_support_weight(::ConfigsAll) = true """ ConfigsMax{K, BOUNDED} <:AbstractProperty @@ -114,7 +108,6 @@ it is finding all independent sets of sizes ``\\alpha(G), \\alpha(G)-1, \\ldots, struct ConfigsMax{K, BOUNDED} <:AbstractProperty end ConfigsMax(K::Int=1; bounded::Bool=true) = ConfigsMax{K,bounded}() max_k(::ConfigsMax{K}) where K = K -_support_weight(::ConfigsMax{1}) = true """ solve(problem, property; usecuda=false, T=Float64) @@ -143,19 +136,16 @@ Keyword arguments * `T` is the "base" element type, sometimes can be used to reduce the memory cost. """ function solve(gp::GraphProblem, property::AbstractProperty; T=Float64, usecuda=false) - if !_support_weight(property) && _has_weight(gp) - throw(ArgumentError("weighted instance of type $(typeof(gp)) is not supported in computing $(property)")) - end if property isa SizeMax syms = symbols(gp) vertex_index = Dict([s=>i for (i, s) in enumerate(syms)]) - return contractf(x->Tropical(T(get_weight(gp, vertex_index[x]))), gp; usecuda=usecuda) + return contractf(x->Tropical{T}.(get_weights(gp, vertex_index[x])), gp; usecuda=usecuda) elseif property isa CountingAll return contractx(gp, one(T); usecuda=usecuda) elseif property isa CountingMax{1} syms = symbols(gp) vertex_index = Dict([s=>i for (i, s) in enumerate(syms)]) - return contractf(x->CountingTropical(T(get_weight(gp, vertex_index[x])), one(T)), gp; usecuda=usecuda) + return contractf(x->CountingTropical{T,T}.(get_weights(gp, vertex_index[x])), gp; usecuda=usecuda) elseif property isa CountingMax return contractx(gp, TruncatedPoly(ntuple(i->i == max_k(property) ? one(T) : zero(T), max_k(property)), one(T)); usecuda=usecuda) elseif property isa GraphPolynomial @@ -178,7 +168,6 @@ function solve(gp::GraphProblem, property::AbstractProperty; T=Float64, usecuda= error("unknown property $property.") end end -_has_weight(gp::GraphProblem) = hasfield(typeof(gp), :weights) && gp.weights isa Vector # ugly but makes life easier """ max_size(problem; usecuda=false) diff --git a/src/networks.jl b/src/networks.jl deleted file mode 100644 index 06677eb7..00000000 --- a/src/networks.jl +++ /dev/null @@ -1,374 +0,0 @@ -const EinTypes = Union{EinCode,NestedEinsum,SlicedEinsum} - -abstract type GraphProblem end - -struct UnWeighted end - -""" - Independence{CT<:EinTypes,WT<:Union{UnWeighted, Vector}} <: GraphProblem - Independence(graph; weights=UnWeighted(), openvertices=(), - optimizer=GreedyMethod(), simplifier=nothing) - -The [Independent set problem](https://en.wikipedia.org/wiki/Independent_set_(graph_theory)). -In the constructor, `weights` are the weights of vertices. -`openvertices` specifies labels for the output tensor. -`optimizer` and `simplifier` are for tensor network optimization, check `optimize_code` for details. - -Problem definition ---------------------------- -An independent set is defined in the [monadic second order logic](https://digitalcommons.oberlin.edu/cgi/viewcontent.cgi?article=1678&context=honors) as -```math -\\exists x_i,\\ldots,x_M\\left[\\bigwedge_{i\\neq j} (x_i\\neq x_j \\wedge \\neg \\textbf{adj}(x_i, x_j))\\right] -``` - -Graph polynomial ---------------------------- -The graph polynomial defined for the independence problem is known as the independence polynomial. -```math -I(G, x) = \\sum_{k=0}^{\\alpha(G)} a_k x^k, -``` -where ``\\alpha(G)`` is the maximum independent set size, -``a_k`` is the number of independent sets of size ``k`` in graph ``G=(V,E)``. -The total number of independent sets is thus equal to ``I(G, 1)``. - -Tensor network ---------------------------- -In tensor network representation of the independent set problem, -we map a vertex ``i\\in V`` to a label ``s_i \\in \\{0, 1\\}`` of dimension 2, -where we use 0 (1) to denote a vertex is absent (present) in the set. -For each label ``s_i``, we defined a parametrized rank-one vertex tensor ``W(x_i)`` as -```math -W(x_i)_{s_i} = \\left(\\begin{matrix} - 1 \\\\ - x_i -\\end{matrix}\\right)_{s_i} -``` -We use subscripts to index tensor elements, e.g.``W(x_i)_0=1`` is the first element associated -with ``s_i=0`` and ``W(x_i)_1=x_i`` is the second element associated with ``s_i=1``. -Similarly, on each edge ``(u, v)``, we define a matrix ``B`` indexed by ``s_u`` and ``s_v`` as -```math -B_{s_i s_j} = \\left(\\begin{matrix} - 1 & 1\\\\ - 1 & 0 -\\end{matrix}\\right)_{s_is_j} -``` - -Its contraction time space complexity is ``2^{{\\rm tw}(G)}``, where ``{\\rm tw}`` is the [tree-width](https://en.wikipedia.org/wiki/Treewidth) of ``G``. -""" -struct Independence{CT<:EinTypes,WT<:Union{UnWeighted, Vector}} <: GraphProblem - code::CT - weights::WT -end - -function Independence(g::SimpleGraph; weights=UnWeighted(), openvertices=(), optimizer=GreedyMethod(), simplifier=nothing) - @assert weights isa UnWeighted || length(weights) == nv(g) - rawcode = EinCode(([[i] for i in Graphs.vertices(g)]..., # labels for vertex tensors - [[minmax(e.src,e.dst)...] for e in Graphs.edges(g)]...), collect(Int, openvertices)) # labels for edge tensors - code = _optimize_code(rawcode, uniformsize(rawcode, 2), optimizer, simplifier) - Independence(code, weights) -end - -""" - MaxCut{CT<:EinTypes,WT<:Union{UnWeighted, Vector}} <: GraphProblem - MaxCut(graph; weights=UnWeighted(), openvertices=(), - optimizer=GreedyMethod(), simplifier=nothing) - -[Cut](https://en.wikipedia.org/wiki/Maximum_cut) problem (or spin glass problem). -In the constructor, `weights` are the weights of edges. -`optimizer` and `simplifier` are for tensor network optimization, check `optimize_code` for details. - -Problem definition ---------------------------- -In graph theory, a cut is a partition of the vertices of a graph into two disjoint subsets. -A maximum cut is a cut whose size is at least the size of any other cut, -where the size of a cut is the number of edges (or the sum of weights on edges) crossing the cut. - -Graph polynomial ---------------------------- -The graph polynomial defined for the cut problem is -```math -C(G, x) = \\sum_{k=0}^{\\gamma(G)} c_k x^k, -``` -where ``\\alpha(G)`` is the maximum independent set size, -``c_k/2`` is the number of cuts of size ``k`` in graph ``G=(V,E)``. - -Tensor network ---------------------------- -For a vertex ``v\\in V``, we define a boolean degree of freedom ``s_v\\in\\{0, 1\\}``. -Then the maximum cut problem can be encoded to tensor networks by mapping an edge ``(i,j)\\in E`` to an edge matrix labelled by ``s_is_j`` -```math -B(x_{\\langle i, j\\rangle}) = \\left(\\begin{matrix} - 1 & x_{\\langle i, j\\rangle}\\\\ - x_{\\langle i, j\\rangle} & 1 -\\end{matrix}\\right), -``` -where variable ``x_{\\langle i, j\\rangle}`` represents a cut on edge ``(i, j)`` or a domain wall of an Ising spin glass. -Similar to other problems, we can define a polynomial about edges variables by setting ``x_{\\langle i, j\\rangle} = x``, -where its k-th coefficient is two times the number of configurations of cut size k. - -Its contraction time space complexity is ``2^{{\\rm tw}(G)}``, where ``{\\rm tw}`` is the [tree-width](https://en.wikipedia.org/wiki/Treewidth) of ``G``. -""" -struct MaxCut{CT<:EinTypes,WT<:Union{UnWeighted, Vector}} <: GraphProblem - code::CT - weights::WT -end -function MaxCut(g::SimpleGraph; weights=UnWeighted(), openvertices=(), optimizer=GreedyMethod(), simplifier=nothing) - @assert weights isa UnWeighted || length(weights) == ne(g) - rawcode = EinCode([[minmax(e.src,e.dst)...] for e in Graphs.edges(g)], collect(Int, openvertices)) # labels for edge tensors - MaxCut(_optimize_code(rawcode, uniformsize(rawcode, 2), optimizer, simplifier), weights) -end - -""" - MaximalIndependence{CT<:EinTypes,WT<:Union{UnWeighted, Vector}} <: GraphProblem - MaximalIndependence(graph; weights=UnWeighted(), openvertices=(), - optimizer=GreedyMethod(), simplifier=nothing) - -[Maximal independent set](https://en.wikipedia.org/wiki/Maximal_independent_set) problem. In the constructor, `weights` are the weights of vertices. -`optimizer` and `simplifier` are for tensor network optimization, check `optimize_code` for details. - -Problem definition ---------------------------- -In graph theory, a maximal independent set is an independent set that is not a subset of any other independent set. -It is different from maximum independent set because it does not require the set to have the max size. - -Graph polynomial ---------------------------- -The graph polynomial defined for the maximal independent set problem is -```math -I_{\\rm max}(G, x) = \\sum_{k=0}^{\\alpha(G)} b_k x^k, -``` -where ``b_k`` is the number of maximal independent sets of size ``k`` in graph ``G=(V, E)``. - -Tensor network ---------------------------- -For a vertex ``v\\in V``, we define a boolean degree of freedom ``s_v\\in\\{0, 1\\}``. -We defined the restriction on its neighbourhood ``N[v]``: -```math -T(x_v)_{s_1,s_2,\\ldots,s_{|N(v)|},s_v} = \\begin{cases} - s_vx_v & s_1=s_2=\\ldots=s_{|N(v)|}=0,\\\\ - 1-s_v& \\text{otherwise}.\\ -\\end{cases} -``` -Intuitively, it means if all the neighbourhood vertices are not in ``I_{m}``, i.e., ``s_1=s_2=\\ldots=s_{|N(v)|}=0``, then ``v`` should be in ``I_{m}`` and contribute a factor ``x_{v}``, -otherwise, if any of the neighbourhood vertices is in ``I_{m}``, then ``v`` cannot be in ``I_{m}``. - -Its contraction time space complexity is no longer determined by the tree-width of the original graph ``G``. -It is often harder to contract this tensor network than to contract the one for regular independent set problem. -""" -struct MaximalIndependence{CT<:EinTypes,WT<:Union{UnWeighted, Vector}} <: GraphProblem - code::CT - weights::WT -end - -function MaximalIndependence(g::SimpleGraph; weights=UnWeighted(), openvertices=(), optimizer=GreedyMethod(), simplifier=nothing) - @assert weights isa UnWeighted || length(weights) == nv(g) - rawcode = EinCode(([[Graphs.neighbors(g, v)..., v] for v in Graphs.vertices(g)]...,), collect(Int, openvertices)) - MaximalIndependence(_optimize_code(rawcode, uniformsize(rawcode, 2), optimizer, simplifier), weights) -end - -""" - Matching{CT<:EinTypes} <: GraphProblem - Matching(graph; openvertices=(), optimizer=GreedyMethod(), simplifier=nothing) - -[Vertex matching](https://mathworld.wolfram.com/Matching.html) problem. -`optimizer` and `simplifier` are for tensor network optimization, check `optimize_code` for details. - -Problem definition ---------------------------- -A ``k``-matching in a graph ``G`` is a set of k edges, no two of which have a vertex in common. - -Graph polynomial ---------------------------- -The matching polynomial adopts the first definition in [wiki page](https://en.wikipedia.org/wiki/Matching_polynomial) -```math -M(G, x) = \\sum\\limits_{k=1}^{|V|/2} c_k x^k, -``` -where ``k`` is the number of matches, and coefficients ``c_k`` are the corresponding counting. - -Tensor network ---------------------------- -We map an edge ``(u, v) \\in E`` to a label ``\\langle u, v\\rangle \\in \\{0, 1\\}`` in a tensor network, -where 1 means two vertices of an edge are matched, 0 means otherwise. -Then we define a tensor of rank ``d(v) = |N(v)|`` on vertex ``v`` such that, -```math -W_{\\langle v, n_1\\rangle, \\langle v, n_2 \\rangle, \\ldots, \\langle v, n_{d(v)}\\rangle} = \\begin{cases} - 1, & \\sum_{i=1}^{d(v)} \\langle v, n_i \\rangle \\leq 1,\\\\ - 0, & \\text{otherwise}, -\\end{cases} -``` -and a tensor of rank 1 on the bond -```math -B_{\\langle v, w\\rangle} = \\begin{cases} -1, & \\langle v, w \\rangle = 0 \\\\ -x, & \\langle v, w \\rangle = 1, -\\end{cases} -``` -where label ``\\langle v, w \\rangle`` is equivalent to ``\\langle w,v\\rangle``. -""" -struct Matching{CT<:EinTypes} <: GraphProblem - code::CT -end - -function Matching(g::SimpleGraph; openvertices=(), optimizer=GreedyMethod(), simplifier=nothing) - rawcode = EinCode(vcat([[minmax(e.src,e.dst)] for e in Graphs.edges(g)], # labels for edge tensors - [[minmax(i,j) for j in neighbors(g, i)] for i in Graphs.vertices(g)]), - collect(Tuple{Int,Int}, openvertices)) - Matching(_optimize_code(rawcode, uniformsize(rawcode, 2), optimizer, simplifier)) -end - -""" - Coloring{K,CT<:EinTypes} <: GraphProblem - Coloring{K}(graph; openvertices=(), optimizer=GreedyMethod(), simplifier=nothing) - -[Vertex Coloring](https://en.wikipedia.org/wiki/Graph_coloring) problem. -`optimizer` and `simplifier` are for tensor network optimization, check `optimize_code` for details. - -Problem definition ---------------------------- -A vertex coloring is an assignment of labels or colors to each vertex of a graph such that no edge connects two identically colored vertices. - -Tensor network ---------------------------- -Let us use 3-colouring problem defined on vertices as an example. -For a vertex ``v``, we define the degree of freedoms ``c_v\\in\\{1,2,3\\}`` and a vertex tensor labelled by it as -```math -W(v) = \\left(\\begin{matrix} - r_v\\\\ - g_v\\\\ - b_v -\\end{matrix}\\right). -``` -For an edge ``(u, v)``, we define an edge tensor as a matrix labelled by ``(c_u, c_v)`` to specify the constraint -```math -B = \\left(\\begin{matrix} - 0 & 1 & 1\\\\ - 1 & 0 & 1\\\\ - 1 & 1 & 0 -\\end{matrix}\\right). -``` -The number of possible colouring can be obtained by contracting this tensor network by setting vertex tensor elements ``r_v, g_v`` and ``b_v`` to 1. -""" -struct Coloring{K,CT<:EinTypes} <: GraphProblem - code::CT -end -Coloring{K}(code::ET) where {K,ET<:EinTypes} = Coloring{K,ET}(code) -# same network layout as independent set. -Coloring{K}(g::SimpleGraph; openvertices=(), optimizer=GreedyMethod(), simplifier=nothing) where K = Coloring{K}(Independence(g; openvertices=openvertices, optimizer=optimizer, simplifier=simplifier).code) - -""" - PaintShop{CT<:EinTypes} <: GraphProblem - PaintShop(labels::AbstractVector; openvertices=(), - optimizer=GreedyMethod(), simplifier=nothing) - -The [binary paint shop problem](http://m-hikari.com/ams/ams-2012/ams-93-96-2012/popovAMS93-96-2012-2.pdf). - -Example ------------------------------------------ -One can encode the paint shop problem `abaccb` as the following - -```jldoctest; setup=:(using GraphTensorNetworks) -julia> symbols = collect("abaccb"); - -julia> pb = PaintShop(symbols); - -julia> solve(pb, SizeMax())[] -3.0ₜ - -julia> solve(pb, ConfigsMax())[].c.data -2-element Vector{StaticBitVector{5, 1}}: - 01101 - 01101 -``` -In our definition, we find the maximum number of unchanged color in this sequence, i.e. (n-1) - (minimum number of color changes) -In the output of maximum configurations, the two configurations are defined on 5 bonds i.e. pairs of (i, i+1), `0` means color changed, while `1` means color not changed. -If we denote two "colors" as `r` and `b`, then the optimal painting is `rbbbrr` or `brrrbb`, both change the colors twice. -""" -struct PaintShop{CT<:EinTypes,LT} <: GraphProblem - code::CT - labels::Vector{LT} - isfirst::Vector{Bool} -end - -function paintshop_from_pairs(pairs::AbstractVector{Tuple{Int,Int}}; openvertices=(), optimizer=GreedyMethod(), simplifier=nothing) - n = length(pairs) - @assert sort!(vcat(collect.(pairs)...)) == collect(1:2n) - labels = zeros(Int, 2*n) - @inbounds for i=1:n - labels[pairs[i]] .= i - end - return PaintShop(pairs; openvertices, optimizer, simplifier) -end -function PaintShop(labels::AbstractVector{T}; openvertices=(), optimizer=GreedyMethod(), simplifier=nothing) where T - @assert all(l->count(==(l), labels)==2, labels) - n = length(labels) - isfirst = [findfirst(==(labels[i]), labels) == i for i=1:n] - rawcode = EinCode(vcat( - [[labels[i], labels[i+1]] for i=1:n-1], # labels for edge tensors - ), - collect(T, openvertices)) - PaintShop(_optimize_code(rawcode, uniformsize(rawcode, 2), optimizer, simplifier), labels, isfirst) -end - -""" - labels(code) - -Return a vector of unique labels in an Einsum token. -""" -function labels(code::EinTypes) - res = [] - for ix in getixsv(code) - for l in ix - if l ∉ res - push!(res, l) - end - end - end - return res -end - -OMEinsum.timespace_complexity(gp::GraphProblem) = timespace_complexity(gp.code, uniformsize(gp.code, bondsize(gp))) - -for T in [:Independence, :Matching, :MaximalIndependence, :MaxCut, :PaintShop] - @eval bondsize(gp::$T) = 2 -end -bondsize(::Coloring{K}) where K = K - -get_weight(gp::GraphProblem, x::Int) = 1 -for T in [:Independence, :MaximalIndependence, :MaxCut] - @eval get_weight(gp::$T, x::Int) = gp.weights isa UnWeighted ? 1 : gp.weights[x] -end - -""" -set_packing(sets; openvertices=(), optimizer=GreedyMethod(), simplifier=nothing) - -Set packing is a generalization of independent set problem to hypergraphs. -Calling this function will return you an `Independence` instance. -`sets` are a vector of vectors, each element being a vertex in the independent set problem. -`optimizer` and `simplifier` are for tensor network optimization, check `optimize_code` for details. - -Example ------------------------------------ -```julia -julia> sets = [[1, 2, 5], [1, 3], [2, 4], [3, 6], [2, 3, 6]]; # each set is a vertex - -julia> gp = set_packing(sets); - -julia> res = best_solutions(gp; all=true)[] -(2, {10010, 00110, 01100})ₜ -``` -""" -function set_packing(sets; weights=UnWeighted(), openvertices=(), optimizer=GreedyMethod(), simplifier=nothing) - n = length(sets) - code = EinCode(vcat([[i] for i=1:n], [[i,j] for i=1:n,j=1:n if j>i && !isempty(sets[i] ∩ sets[j])]), collect(Int,openvertices)) - Independence(_optimize_code(code, uniformsize(code, 2), optimizer, simplifier), weights) -end - -_optimize_code(code, size_dict, optimizer::Nothing, simplifier) = code -_optimize_code(code, size_dict, optimizer, simplifier) = optimize_code(code, size_dict, optimizer, simplifier) - -# TODOs: -# 1. Dominating set -# \exists x_i,\ldots,x_K \forall y\left[\bigwedge_{i=1}^{K}(y=x_i\wedge \textbf{adj}(y, x_i))\right] -# 2. Polish reading data -# * consistent configuration assign of max-cut -# 3. Support transverse field in max-cut \ No newline at end of file diff --git a/src/networks/Coloring.jl b/src/networks/Coloring.jl new file mode 100644 index 00000000..87d1ae8f --- /dev/null +++ b/src/networks/Coloring.jl @@ -0,0 +1,64 @@ +""" + Coloring{K,CT<:AbstractEinsum} <: GraphProblem + Coloring{K}(graph; openvertices=(), optimizer=GreedyMethod(), simplifier=nothing) + +[Vertex Coloring](https://en.wikipedia.org/wiki/Graph_coloring) problem. +`optimizer` and `simplifier` are for tensor network optimization, check `optimize_code` for details. + +Problem definition +--------------------------- +A vertex coloring is an assignment of labels or colors to each vertex of a graph such that no edge connects two identically colored vertices. + +Tensor network +--------------------------- +Let us use 3-colouring problem defined on vertices as an example. +For a vertex ``v``, we define the degree of freedoms ``c_v\\in\\{1,2,3\\}`` and a vertex tensor labelled by it as +```math +W(v) = \\left(\\begin{matrix} + r_v\\\\ + g_v\\\\ + b_v +\\end{matrix}\\right). +``` +For an edge ``(u, v)``, we define an edge tensor as a matrix labelled by ``(c_u, c_v)`` to specify the constraint +```math +B = \\left(\\begin{matrix} + 0 & 1 & 1\\\\ + 1 & 0 & 1\\\\ + 1 & 1 & 0 +\\end{matrix}\\right). +``` +The number of possible colouring can be obtained by contracting this tensor network by setting vertex tensor elements ``r_v, g_v`` and ``b_v`` to 1. +""" +struct Coloring{K,CT<:AbstractEinsum} <: GraphProblem + code::CT + nv::Int +end +Coloring{K}(code::ET, nv::Int) where {K,ET<:AbstractEinsum} = Coloring{K,ET}(code, nv) +# same network layout as independent set. +Coloring{K}(g::SimpleGraph; openvertices=(), optimizer=GreedyMethod(), simplifier=nothing) where K = Coloring{K}(Independence(g; openvertices=openvertices, optimizer=optimizer, simplifier=simplifier).code, nv(g)) + +flavors(::Type{<:Coloring{K}}) where K = collect(0:K-1) +symbols(c::Coloring{K}) where K = [i for i=1:c.nv] +get_weights(::Coloring{K}, sym) where K = ones(Int, K) + +# `fx` is a function defined on symbols, it returns a vector of elements, the size of vector is same as the number of flavors (or the bond dimension). +function generate_tensors(fx, c::Coloring{K}) where K + ixs = getixsv(c.code) + T = eltype(fx(ixs[1][1])) + return map(1:length(ixs)) do i + i <= c.nv ? coloringv(fx(ixs[i][1])) : coloringb(T, K) + end +end + +# coloring bond tensor +function coloringb(::Type{T}, k::Int) where T + x = ones(T, k, k) + for i=1:k + x[i,i] = zero(T) + end + return x +end +# coloring vertex tensor +coloringv(vals::Vector{T}) where T = vals + diff --git a/src/networks/Independence.jl b/src/networks/Independence.jl new file mode 100644 index 00000000..0bae3554 --- /dev/null +++ b/src/networks/Independence.jl @@ -0,0 +1,120 @@ +""" + Independence{CT<:AbstractEinsum,WT<:Union{UnWeighted, Vector}} <: GraphProblem + Independence(graph; weights=UnWeighted(), openvertices=(), + optimizer=GreedyMethod(), simplifier=nothing) + +The [Independent set problem](https://en.wikipedia.org/wiki/Independent_set_(graph_theory)). +In the constructor, `weights` are the weights of vertices. +`openvertices` specifies labels for the output tensor. +`optimizer` and `simplifier` are for tensor network optimization, check `optimize_code` for details. + +Problem definition +--------------------------- +An independent set is defined in the [monadic second order logic](https://digitalcommons.oberlin.edu/cgi/viewcontent.cgi?article=1678&context=honors) as +```math +\\exists x_i,\\ldots,x_M\\left[\\bigwedge_{i\\neq j} (x_i\\neq x_j \\wedge \\neg \\textbf{adj}(x_i, x_j))\\right] +``` + +Graph polynomial +--------------------------- +The graph polynomial defined for the independence problem is known as the independence polynomial. +```math +I(G, x) = \\sum_{k=0}^{\\alpha(G)} a_k x^k, +``` +where ``\\alpha(G)`` is the maximum independent set size, +``a_k`` is the number of independent sets of size ``k`` in graph ``G=(V,E)``. +The total number of independent sets is thus equal to ``I(G, 1)``. + +Tensor network +--------------------------- +In tensor network representation of the independent set problem, +we map a vertex ``i\\in V`` to a label ``s_i \\in \\{0, 1\\}`` of dimension 2, +where we use 0 (1) to denote a vertex is absent (present) in the set. +For each label ``s_i``, we defined a parametrized rank-one vertex tensor ``W(x_i)`` as +```math +W(x_i)_{s_i} = \\left(\\begin{matrix} + 1 \\\\ + x_i +\\end{matrix}\\right)_{s_i} +``` +We use subscripts to index tensor elements, e.g.``W(x_i)_0=1`` is the first element associated +with ``s_i=0`` and ``W(x_i)_1=x_i`` is the second element associated with ``s_i=1``. +Similarly, on each edge ``(u, v)``, we define a matrix ``B`` indexed by ``s_u`` and ``s_v`` as +```math +B_{s_i s_j} = \\left(\\begin{matrix} + 1 & 1\\\\ + 1 & 0 +\\end{matrix}\\right)_{s_is_j} +``` + +Its contraction time space complexity is ``2^{{\\rm tw}(G)}``, where ``{\\rm tw(G)}`` is the [tree-width](https://en.wikipedia.org/wiki/Treewidth) of ``G``. +""" +struct Independence{CT<:AbstractEinsum,WT<:Union{UnWeighted, Vector}} <: GraphProblem + code::CT + nv::Int + weights::WT +end + +function Independence(g::SimpleGraph; weights=UnWeighted(), openvertices=(), optimizer=GreedyMethod(), simplifier=nothing) + @assert weights isa UnWeighted || length(weights) == nv(g) + rawcode = EinCode(([[i] for i in Graphs.vertices(g)]..., # labels for vertex tensors + [[minmax(e.src,e.dst)...] for e in Graphs.edges(g)]...), collect(Int, openvertices)) # labels for edge tensors + code = _optimize_code(rawcode, uniformsize(rawcode, 2), optimizer, simplifier) + Independence(code, nv(g), weights) +end + +flavors(::Type{<:Independence}) = [0, 1] +symbols(gp::Independence) = [i for i in 1:gp.nv] +get_weights(gp::Independence, label) = [0, gp.weights isa UnWeighted ? 1 : gp.weights[findfirst(==(label), symbols(gp))]] + +# generate tensors +function generate_tensors(fx, gp::Independence) + syms = symbols(gp) + isempty(syms) && return [] + ixs = getixsv(gp.code) + T = eltype(fx(syms[1])) + return map(enumerate(ixs)) do (i, ix) + if i <= length(syms) + misv(fx(ix[1])) + else + misb(T, length(ix)) # if n!=2, it corresponds to set packing problem. + end + end +end + +function misb(::Type{T}, n::Integer=2) where T + res = zeros(T, fill(2, n)...) + res[1] = one(T) + for i=1:n + res[1+1<<(i-1)] = one(T) + end + return res +end +misv(vals) = vals + +############### set packing ##################### +""" +set_packing(sets; openvertices=(), optimizer=GreedyMethod(), simplifier=nothing) + +Set packing is a generalization of independent set problem to hypergraphs. +Calling this function will return you an `Independence` instance. +`sets` are a vector of vectors, each element being a vertex in the independent set problem. +`optimizer` and `simplifier` are for tensor network optimization, check `optimize_code` for details. + +Example +----------------------------------- +```julia +julia> sets = [[1, 2, 5], [1, 3], [2, 4], [3, 6], [2, 3, 6]]; # each set is a vertex + +julia> gp = set_packing(sets); + +julia> res = best_solutions(gp; all=true)[] +(2, {10010, 00110, 01100})ₜ +``` +""" +function set_packing(sets; weights=UnWeighted(), openvertices=(), optimizer=GreedyMethod(), simplifier=nothing) + n = length(sets) + code = EinCode(vcat([[i] for i=1:n], [[i,j] for i=1:n,j=1:n if j>i && !isempty(sets[i] ∩ sets[j])]), collect(Int,openvertices)) + Independence(_optimize_code(code, uniformsize(code, 2), optimizer, simplifier), n, weights) +end + diff --git a/src/networks/Matching.jl b/src/networks/Matching.jl new file mode 100644 index 00000000..a4768367 --- /dev/null +++ b/src/networks/Matching.jl @@ -0,0 +1,84 @@ +""" + Matching{CT<:AbstractEinsum} <: GraphProblem + Matching(graph; openvertices=(), optimizer=GreedyMethod(), simplifier=nothing) + +[Vertex matching](https://mathworld.wolfram.com/Matching.html) problem. +`optimizer` and `simplifier` are for tensor network optimization, check `optimize_code` for details. + +Problem definition +--------------------------- +A ``k``-matching in a graph ``G`` is a set of k edges, no two of which have a vertex in common. + +Graph polynomial +--------------------------- +The matching polynomial adopts the first definition in [wiki page](https://en.wikipedia.org/wiki/Matching_polynomial) +```math +M(G, x) = \\sum\\limits_{k=1}^{|V|/2} c_k x^k, +``` +where ``k`` is the number of matches, and coefficients ``c_k`` are the corresponding counting. + +Tensor network +--------------------------- +We map an edge ``(u, v) \\in E`` to a label ``\\langle u, v\\rangle \\in \\{0, 1\\}`` in a tensor network, +where 1 means two vertices of an edge are matched, 0 means otherwise. +Then we define a tensor of rank ``d(v) = |N(v)|`` on vertex ``v`` such that, +```math +W_{\\langle v, n_1\\rangle, \\langle v, n_2 \\rangle, \\ldots, \\langle v, n_{d(v)}\\rangle} = \\begin{cases} + 1, & \\sum_{i=1}^{d(v)} \\langle v, n_i \\rangle \\leq 1,\\\\ + 0, & \\text{otherwise}, +\\end{cases} +``` +and a tensor of rank 1 on the bond +```math +B_{\\langle v, w\\rangle} = \\begin{cases} +1, & \\langle v, w \\rangle = 0 \\\\ +x, & \\langle v, w \\rangle = 1, +\\end{cases} +``` +where label ``\\langle v, w \\rangle`` is equivalent to ``\\langle w,v\\rangle``. +""" +struct Matching{CT<:AbstractEinsum} <: GraphProblem + symbols::Vector{Tuple{Int,Int}} + code::CT +end + +function Matching(g::SimpleGraph; openvertices=(), optimizer=GreedyMethod(), simplifier=nothing) + symbols = [minmax(e.src,e.dst) for e in Graphs.edges(g)] + rawcode = EinCode(vcat([[s] for s in symbols], # labels for edge tensors + [[minmax(i,j) for j in neighbors(g, i)] for i in Graphs.vertices(g)]), + collect(Tuple{Int,Int}, openvertices)) + Matching(symbols, _optimize_code(rawcode, uniformsize(rawcode, 2), optimizer, simplifier)) +end + +flavors(::Type{<:Matching}) = [0, 1] +get_weights(::Matching, sym) = [0, 1] +symbols(m::Matching) = m.symbols + +function generate_tensors(fx, m::Matching) + syms = symbols(m) + isempty(syms) && return [] + ixs = getixsv(m.code) + T = eltype(fx(syms[1])) + n = length(syms) # number of vertices + tensors = [] + for i=1:length(ixs) + if i<=n + @assert length(ixs[i]) == 1 + t = fx(syms[i]) # fx is defined on edges. + else + t = match_tensor(T, length(ixs[i])) + end + push!(tensors, t) + end + return tensors +end +function match_tensor(::Type{T}, n::Int) where T + t = zeros(T, fill(2, n)...) + for ci in CartesianIndices(t) + if sum(ci.I .- 1) <= 1 + t[ci] = one(T) + end + end + return t +end + diff --git a/src/networks/MaxCut.jl b/src/networks/MaxCut.jl new file mode 100644 index 00000000..26c20f33 --- /dev/null +++ b/src/networks/MaxCut.jl @@ -0,0 +1,64 @@ +""" + MaxCut{CT<:AbstractEinsum,WT<:Union{UnWeighted, Vector}} <: GraphProblem + MaxCut(graph; weights=UnWeighted(), openvertices=(), + optimizer=GreedyMethod(), simplifier=nothing) + +[Cut](https://en.wikipedia.org/wiki/Maximum_cut) problem (or spin glass problem). +In the constructor, `weights` are the weights of edges. +`optimizer` and `simplifier` are for tensor network optimization, check `optimize_code` for details. + +Problem definition +--------------------------- +In graph theory, a cut is a partition of the vertices of a graph into two disjoint subsets. +A maximum cut is a cut whose size is at least the size of any other cut, +where the size of a cut is the number of edges (or the sum of weights on edges) crossing the cut. + +Graph polynomial +--------------------------- +The graph polynomial defined for the cut problem is +```math +C(G, x) = \\sum_{k=0}^{\\gamma(G)} c_k x^k, +``` +where ``\\alpha(G)`` is the maximum independent set size, +``c_k/2`` is the number of cuts of size ``k`` in graph ``G=(V,E)``. + +Tensor network +--------------------------- +For a vertex ``v\\in V``, we define a boolean degree of freedom ``s_v\\in\\{0, 1\\}``. +Then the maximum cut problem can be encoded to tensor networks by mapping an edge ``(i,j)\\in E`` to an edge matrix labelled by ``s_is_j`` +```math +B(x_{\\langle i, j\\rangle}) = \\left(\\begin{matrix} + 1 & x_{\\langle i, j\\rangle}\\\\ + x_{\\langle i, j\\rangle} & 1 +\\end{matrix}\\right), +``` +where variable ``x_{\\langle i, j\\rangle}`` represents a cut on edge ``(i, j)`` or a domain wall of an Ising spin glass. +Similar to other problems, we can define a polynomial about edges variables by setting ``x_{\\langle i, j\\rangle} = x``, +where its k-th coefficient is two times the number of configurations of cut size k. + +Its contraction time space complexity is ``2^{{\\rm tw}(G)}``, where ``{\\rm tw(G)}`` is the [tree-width](https://en.wikipedia.org/wiki/Treewidth) of ``G``. +""" +struct MaxCut{CT<:AbstractEinsum,WT<:Union{UnWeighted, Vector}} <: GraphProblem + code::CT + weights::WT +end +function MaxCut(g::SimpleGraph; weights=UnWeighted(), openvertices=(), optimizer=GreedyMethod(), simplifier=nothing) + @assert weights isa UnWeighted || length(weights) == ne(g) + rawcode = EinCode([[minmax(e.src,e.dst)...] for e in Graphs.edges(g)], collect(Int, openvertices)) # labels for edge tensors + MaxCut(_optimize_code(rawcode, uniformsize(rawcode, 2), optimizer, simplifier), weights) +end + +flavors(::Type{<:MaxCut}) = [0, 1] +symbols(gp::MaxCut) = getixsv(gp.code) +get_weights(gp::MaxCut, label) = [0, gp.weights isa UnWeighted ? 1 : gp.weights[findfirst(==(label), symbols(gp))]] + +function generate_tensors(fx, gp::MaxCut) + ixs = getixsv(gp.code) + return map(enumerate(ixs)) do (i, ix) + maxcutb(fx(ix)...) + end +end +function maxcutb(a, b) + return [a b; b a] +end + diff --git a/src/networks/MaximalIndependence.jl b/src/networks/MaximalIndependence.jl new file mode 100644 index 00000000..aff0a7a2 --- /dev/null +++ b/src/networks/MaximalIndependence.jl @@ -0,0 +1,70 @@ +""" + MaximalIndependence{CT<:AbstractEinsum,WT<:Union{UnWeighted, Vector}} <: GraphProblem + MaximalIndependence(graph; weights=UnWeighted(), openvertices=(), + optimizer=GreedyMethod(), simplifier=nothing) + +[Maximal independent set](https://en.wikipedia.org/wiki/Maximal_independent_set) problem. In the constructor, `weights` are the weights of vertices. +`optimizer` and `simplifier` are for tensor network optimization, check `optimize_code` for details. + +Problem definition +--------------------------- +In graph theory, a maximal independent set is an independent set that is not a subset of any other independent set. +It is different from maximum independent set because it does not require the set to have the max size. + +Graph polynomial +--------------------------- +The graph polynomial defined for the maximal independent set problem is +```math +I_{\\rm max}(G, x) = \\sum_{k=0}^{\\alpha(G)} b_k x^k, +``` +where ``b_k`` is the number of maximal independent sets of size ``k`` in graph ``G=(V, E)``. + +Tensor network +--------------------------- +For a vertex ``v\\in V``, we define a boolean degree of freedom ``s_v\\in\\{0, 1\\}``. +We defined the restriction on its neighbourhood ``N[v]``: +```math +T(x_v)_{s_1,s_2,\\ldots,s_{|N(v)|},s_v} = \\begin{cases} + s_vx_v & s_1=s_2=\\ldots=s_{|N(v)|}=0,\\\\ + 1-s_v& \\text{otherwise}.\\ +\\end{cases} +``` +Intuitively, it means if all the neighbourhood vertices are not in ``I_{m}``, i.e., ``s_1=s_2=\\ldots=s_{|N(v)|}=0``, then ``v`` should be in ``I_{m}`` and contribute a factor ``x_{v}``, +otherwise, if any of the neighbourhood vertices is in ``I_{m}``, then ``v`` cannot be in ``I_{m}``. + +Its contraction time space complexity is no longer determined by the tree-width of the original graph ``G``. +It is often harder to contract this tensor network than to contract the one for regular independent set problem. +""" +struct MaximalIndependence{CT<:AbstractEinsum,WT<:Union{UnWeighted, Vector}} <: GraphProblem + code::CT + weights::WT +end + +function MaximalIndependence(g::SimpleGraph; weights=UnWeighted(), openvertices=(), optimizer=GreedyMethod(), simplifier=nothing) + @assert weights isa UnWeighted || length(weights) == nv(g) + rawcode = EinCode(([[Graphs.neighbors(g, v)..., v] for v in Graphs.vertices(g)]...,), collect(Int, openvertices)) + MaximalIndependence(_optimize_code(rawcode, uniformsize(rawcode, 2), optimizer, simplifier), weights) +end + +flavors(::Type{<:MaximalIndependence}) = [0, 1] +symbols(gp::MaximalIndependence) = [i for i in 1:length(getixsv(gp.code))] +get_weights(gp::MaximalIndependence, label) = [0, gp.weights isa UnWeighted ? 1 : gp.weights[findfirst(==(label), symbols(gp))]] + +function generate_tensors(fx, mi::MaximalIndependence) + ixs = getixsv(mi.code) + isempty(ixs) && return [] + T = eltype(fx(ixs[1][end])) + return map(ixs) do ix + neighbortensor(fx(ix[end])..., length(ix)) + end +end +function neighbortensor(a::T, b::T, d::Int) where T + t = zeros(T, fill(2, d)...) + for i = 2:1<<(d-1) + t[i] = one(T) + end + t[1<<(d-1)+1] = a + t[1<<(d-1)+1] = b + return t +end + diff --git a/src/networks/PaintShop.jl b/src/networks/PaintShop.jl new file mode 100644 index 00000000..1a3ced7a --- /dev/null +++ b/src/networks/PaintShop.jl @@ -0,0 +1,70 @@ +""" + PaintShop{CT<:AbstractEinsum} <: GraphProblem + PaintShop(labels::AbstractVector; openvertices=(), + optimizer=GreedyMethod(), simplifier=nothing) + +The [binary paint shop problem](http://m-hikari.com/ams/ams-2012/ams-93-96-2012/popovAMS93-96-2012-2.pdf). + +Example +----------------------------------------- +One can encode the paint shop problem `abaccb` as the following + +```jldoctest; setup=:(using GraphTensorNetworks) +julia> symbols = collect("abaccb"); + +julia> pb = PaintShop(symbols); + +julia> solve(pb, SizeMax())[] +3.0ₜ + +julia> solve(pb, ConfigsMax())[].c.data +2-element Vector{StaticBitVector{5, 1}}: + 01101 + 01101 +``` +In our definition, we find the maximum number of unchanged color in this sequence, i.e. (n-1) - (minimum number of color changes) +In the output of maximum configurations, the two configurations are defined on 5 bonds i.e. pairs of (i, i+1), `0` means color changed, while `1` means color not changed. +If we denote two "colors" as `r` and `b`, then the optimal painting is `rbbbrr` or `brrrbb`, both change the colors twice. +""" +struct PaintShop{CT<:AbstractEinsum,LT} <: GraphProblem + code::CT + labels::Vector{LT} + isfirst::Vector{Bool} +end + +function paintshop_from_pairs(pairs::AbstractVector{Tuple{Int,Int}}; openvertices=(), optimizer=GreedyMethod(), simplifier=nothing) + n = length(pairs) + @assert sort!(vcat(collect.(pairs)...)) == collect(1:2n) + labels = zeros(Int, 2*n) + @inbounds for i=1:n + labels[pairs[i]] .= i + end + return PaintShop(pairs; openvertices, optimizer, simplifier) +end +function PaintShop(labels::AbstractVector{T}; openvertices=(), optimizer=GreedyMethod(), simplifier=nothing) where T + @assert all(l->count(==(l), labels)==2, labels) + n = length(labels) + isfirst = [findfirst(==(labels[i]), labels) == i for i=1:n] + rawcode = EinCode(vcat( + [[labels[i], labels[i+1]] for i=1:n-1], # labels for edge tensors + ), + collect(T, openvertices)) + PaintShop(_optimize_code(rawcode, uniformsize(rawcode, 2), optimizer, simplifier), labels, isfirst) +end + +flavors(::Type{<:PaintShop}) = [0, 1] +get_weights(::PaintShop, sym) = [0, 1] +symbols(gp::PaintShop) = getixsv(gp.code) + +function generate_tensors(fx, c::PaintShop) + ixs = getixsv(c.code) + [paintshop_bond_tensor(fx(ixs[i])..., c.isfirst[i], c.isfirst[i+1]) for i=1:length(ixs)] +end + +function paintshop_bond_tensor(a::T, b::T, if1::Bool, if2::Bool) where T + m = T[b a; a b] + m = if1 ? m : m[[2,1],:] + m = if2 ? m : m[:,[2,1]] + return m +end + diff --git a/src/networks/networks.jl b/src/networks/networks.jl new file mode 100644 index 00000000..c9413e80 --- /dev/null +++ b/src/networks/networks.jl @@ -0,0 +1,119 @@ +""" + GraphProblem + +The abstract base type of graph problems. +""" +abstract type GraphProblem end + +struct UnWeighted end + +######## Interfaces for graph problems ########## +""" + get_weights(problem::GraphProblem, sym) + +The weights for the degree of freedom specified by `sym` of the graph problem, where `sym` is a symbol. +In graph polynomial, integer weights are the orders of `x`. +""" +function get_weights end + +""" + symbols(problem::GraphProblem) + +The symbols of a graph problem, they are the degrees of freedoms in the graph. +""" +function symbols end + +""" + flavors(::Type{<:GraphProblem}) + +It returns a vector of integers as the flavors of a degree of freedom. +Its size is the same as the degree of freedom on a single vertex/edge. +""" +flavors(::GT) where GT<:GraphProblem = flavors(GT) + +""" + nflavor(::Type{<:GraphProblem}) + +Bond size is equal to the number of flavors. +""" +nflavor(::Type{GT}) where GT = length(flavors(GT)) +nflavor(::GT) where GT<:GraphProblem = nflavor(GT) + +""" + generate_tensors(func, problem::GraphProblem) + +Generate a vector of tensors as the inputs of the tensor network contraction code `problem.code`. +`func` is a function to customize the tensors. +`func(symbol)` returns a vector of elements, the length of which is same as the number of flavors. + +Example +-------------------------- +The following code gives your the maximum independent set size +```jldoctest +julia> using Graphs, GraphTensorNetworks + +julia> gp = Independence(smallgraph(:petersen)); + +julia> f(x) = Tropical.([0, 1.0]) +f (generic function with 1 method) + +julia> getixsv(gp.code) +25-element Vector{Vector{Int64}}: + [1] + [2] + [3] + [4] + [5] + [6] + [7] + [8] + [9] + [10] + ⋮ + [3, 8] + [4, 5] + [4, 9] + [5, 10] + [6, 8] + [6, 9] + [7, 9] + [7, 10] + [8, 10] + +julia> gp.code(GraphTensorNetworks.generate_tensors(f, gp)...) +0-dimensional Array{TropicalNumbers.TropicalF64, 0}: +4.0ₜ +``` +""" +generate_tensors(::Type{GT}) where GT = length(flavors(GT)) + +# requires field `code` + +include("Independence.jl") +include("MaximalIndependence.jl") +include("MaxCut.jl") +include("Matching.jl") +include("Coloring.jl") +include("PaintShop.jl") + +# forward the time, space and readwrite complexity +OMEinsum.timespacereadwrite_complexity(gp::GraphProblem) = timespacereadwrite_complexity(gp.code, uniformsize(gp.code, nflavor(gp))) + +# contract the graph tensor network +function contractf(f, gp::GraphProblem; usecuda=false) + @debug "generating tensors ..." + xs = generate_tensors(f, gp) + @debug "contracting tensors ..." + if usecuda + gp.code([CuArray(x) for x in xs]...) + else + gp.code(xs...) + end +end + +# TODOs: +# 1. Dominating set +# \exists x_i,\ldots,x_K \forall y\left[\bigwedge_{i=1}^{K}(y=x_i\wedge \textbf{adj}(y, x_i))\right] +# 2. Polish reading data +# * consistent configuration assign of max-cut +# 3. Support transverse field in max-cut \ No newline at end of file diff --git a/src/utils.jl b/src/utils.jl index e69de29b..cdd5ede9 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -0,0 +1,16 @@ +# Return a vector of unique labels in an Einsum token. +function labels(code::AbstractEinsum) + res = [] + for ix in getixsv(code) + for l in ix + if l ∉ res + push!(res, l) + end + end + end + return res +end + +# a unified interface to optimize the contraction code +_optimize_code(code, size_dict, optimizer::Nothing, simplifier) = code +_optimize_code(code, size_dict, optimizer, simplifier) = optimize_code(code, size_dict, optimizer, simplifier) diff --git a/test/bounding.jl b/test/bounding.jl index 5a8a978b..afe575a2 100644 --- a/test/bounding.jl +++ b/test/bounding.jl @@ -25,7 +25,7 @@ end rawcode = Independence(random_regular_graph(10, 3); optimizer=nothing).code optcode = Independence(random_regular_graph(10, 3); optimizer=GreedyMethod()).code xs = map(OMEinsum.getixs(rawcode)) do ix - length(ix)==1 ? GraphTensorNetworks.misv(TropicalF64(1.0)) : GraphTensorNetworks.misb(TropicalF64) + length(ix)==1 ? GraphTensorNetworks.misv([one(TropicalF64), TropicalF64(1.0)]) : GraphTensorNetworks.misb(TropicalF64) end y1 = rawcode(xs...) y2 = bounding_contract(AllConfigs{1}(), rawcode, xs, BitArray(fill(true)), xs) diff --git a/test/configurations.jl b/test/configurations.jl index 830253d7..30ea8cf0 100644 --- a/test/configurations.jl +++ b/test/configurations.jl @@ -32,7 +32,7 @@ end @testset "enumerating" begin rawcode = Independence(random_regular_graph(10, 3); optimizer=nothing) - optcode = Independence(optimize_code(rawcode.code, uniformsize(rawcode.code, 2), GreedyMethod()), UnWeighted()) + optcode = Independence(optimize_code(rawcode.code, uniformsize(rawcode.code, 2), GreedyMethod()), 10, UnWeighted()) for code in [rawcode, optcode] res0 = max_size(code) _, res1 = max_size_count(code) diff --git a/test/graph_polynomials.jl b/test/graph_polynomials.jl index 2c398c0d..f6258894 100644 --- a/test/graph_polynomials.jl +++ b/test/graph_polynomials.jl @@ -5,7 +5,7 @@ using GraphTensorNetworks: StaticBitVector @testset "bond and vertex tensor" begin @test GraphTensorNetworks.misb(TropicalF64) == [TropicalF64(0) TropicalF64(0); TropicalF64(0) TropicalF64(-Inf)] - @test GraphTensorNetworks.misv(TropicalF64(2.0)) == [TropicalF64(0), TropicalF64(2.0)] + @test GraphTensorNetworks.misv([one(TropicalF64), TropicalF64(2.0)]) == [TropicalF64(0), TropicalF64(2.0)] end @testset "graph generator" begin diff --git a/test/interfaces.jl b/test/interfaces.jl index 9f4f3909..6d0fe0a6 100644 --- a/test/interfaces.jl +++ b/test/interfaces.jl @@ -131,7 +131,8 @@ end @test res1 == Tropical(24.0) res2 = solve(gp, CountingMax(1))[] @test res2 == CountingTropical(24.0, 1.0) - @test_throws ArgumentError solve(gp, CountingMax(2))[] + res2b = solve(gp, CountingMax(2))[] + @test res2b == Max2Poly(1.0, 1.0, 24.0) res3 = solve(gp, SingleConfigMax(; bounded=false))[] res3b = solve(gp, SingleConfigMax(; bounded=true))[] @test res3 == res3b diff --git a/test/runtests.jl b/test/runtests.jl index a4370589..34bfd784 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -25,9 +25,11 @@ end include("interfaces.jl") end +@testset "visualize" begin + include("visualize.jl") +end + # -------------- # doctests # -------------- -@static if VERSION >= v"1.1.0" - doctest(GraphTensorNetworks) -end +doctest(GraphTensorNetworks) diff --git a/test/visualize.jl b/test/visualize.jl new file mode 100644 index 00000000..caea55ff --- /dev/null +++ b/test/visualize.jl @@ -0,0 +1,7 @@ +using GraphTensorNetworks, Test, Graphs + +@testset "visualize" begin + locations = [(1.0, 2.0), (2.0, 3.0)] + @test show_graph(locations, [(1, 2)]) isa Any + @test show_graph(smallgraph(:petersen); locs=[(randn(), randn()) for i=1:10]) isa Any +end \ No newline at end of file