Refactor the graph problems

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,

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

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")

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->i<K ? zero(BS) : onehotv(BS, x, v), K),one(OS))
+    TruncatedPoly{K,BS,OS}(ntuple(i->i<K ? zero(BS) : onehotv(BS, x, v), K),one(OS))^w
 end
 function _onehotv(::Type{CountingTropical{TV,BS}}, x, v, w) where {TV,BS}
     CountingTropical{TV,BS}(TV(w), onehotv(BS, x, v))
 end
 function _onehotv(::Type{BS}, x, v, w) where {BS<:ConfigEnumerator}
     onehotv(BS, x, v)
-end
-
-for GP in [:Independence, :Matching, :MaximalIndependence, :Coloring]
-    @eval symbols(gp::$GP) = labels(gp.code)
-end
-symbols(gp::MaxCut) = getixsv(gp.code)
-symbols(gp::PaintShop) = getixsv(gp.code)
+end

src/graph_polynomials.jl

@@ -96,130 +96,6 @@ function improved_counting(ys::AbstractArray...)
 end
 improved_counting(ys::Mod...) = Mods.CRT(ys...)
 
-contractx(gp::GraphProblem, x; usecuda=false) = contractf(_->x, 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

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)