new opt method auto

Author: GiggleLiu

Project.toml

@@ -24,18 +24,18 @@ Viznet = "52a3aca4-6234-47fd-b74a-806bdf78ede9"
 [compat]
 CUDA = "3.3"
 Compose = "0.9"
+FFTW = "1.4"
 LightGraphs = "1.3"
 Mods = "1.3"
 OMEinsum = "0.4"
-OMEinsumContractionOrders = "0.1"
+OMEinsumContractionOrders = "0.2"
 Polynomials = "2.0"
 Primes = "0.5"
 Requires = "1"
 TropicalGEMM = "0.1"
 TropicalNumbers = "0.4"
 TupleTools = "1.2"
 Viznet = "0.3"
-FFTW = "1.4"
 julia = "1"
 
 [extras]

src/arithematics.jl

@@ -1,5 +1,6 @@
 export is_commutative_semiring
 export Max2Poly, Polynomial, Tropical, CountingTropical, StaticBitVector, Mod, ConfigEnumerator, onehotv
+export bitstringset_type, bitstringsampler_type
 
 using Polynomials: Polynomial
 using TropicalNumbers: Tropical, CountingTropical, StaticBitVector

src/bounding.jl

@@ -1,7 +1,5 @@
 using TupleTools
 
-export bounding_contract
-
 """
     backward_tropical(mode, ixs, xs, iy, y, ymask, size_dict)
 
@@ -85,6 +83,15 @@ function masked_einsum(code::NestedEinsum, @nospecialize(xs), masks, size_dict)
     y[OMEinsum.asarray(.!masks.content)] .= Ref(zero(eltype(y))); y
 end
 
+"""
+    bounding_contract(code, xsa, ymask, xsb; size_info=nothing)
+
+Contraction method with bounding.
+
+    * `xsa` are input tensors for bounding, e.g. tropical tensors,
+    * `xsb` are input tensors for computing, e.g. tensors elements are counting tropical with set algebra,
+    * `ymask` is the initial gradient mask for the output tensor.
+"""
 function bounding_contract(@nospecialize(code::EinCode), @nospecialize(xsa), ymask, @nospecialize(xsb); size_info=nothing)
     bounding_contract(NestedEinsum((1:length(xsa)), code), xsa, ymask, xsb; size_info=size_info)
 end
@@ -100,11 +107,12 @@ function bounding_contract(code::NestedEinsum, @nospecialize(xsa), ymask, @nospe
     masked_einsum(code, xsb, mt, size_dict)
 end
 
-function mis_config_ad(@nospecialize(code::EinCode), @nospecialize(xsa), ymask; size_info=nothing)
-    mis_config_ad(NestedEinsum((1:length(xsa)), code), xsa, ymask; size_info=size_info)
+# get the optimal solution with automatic differentiation.
+function solution_ad(@nospecialize(code::EinCode), @nospecialize(xsa), ymask; size_info=nothing)
+    solution_ad(NestedEinsum((1:length(xsa)), code), xsa, ymask; size_info=size_info)
 end
 
-function mis_config_ad(code::NestedEinsum, @nospecialize(xsa), ymask; size_info=nothing)
+function solution_ad(code::NestedEinsum, @nospecialize(xsa), ymask; size_info=nothing)
     size_dict = OMEinsum.get_size_dict(getixs(flatten(code)), xsa, size_info)
     # compute intermediate tensors
     @debug "caching einsum..."

src/configurations.jl

@@ -1,6 +1,14 @@
-export getconfigs_bounded, getconfigs_direct
+export optimalsolutions, solutions
 
-function getconfigs_bounded(gp::GraphProblem; all=false, usecuda=false)
+"""
+    optimalsolutions(problem; all=false, usecuda=false)
+    
+Find optimal solutions with bounding.
+
+* When `all` is true, the program will use set for enumerate all possible solutions, otherwise, it will return one solution for each size.
+* `usecuda` can not be true if you want to use set to enumerate all possible solutions.
+"""
+function optimalsolutions(gp::GraphProblem; all=false, usecuda=false)
     if all && usecuda
         throw(ArgumentError("ConfigEnumerator can not be computed on GPU!"))
     end
@@ -18,22 +26,34 @@ function getconfigs_bounded(gp::GraphProblem; all=false, usecuda=false)
         return bounding_contract(gp.code, xst, ymask, xs)
     else
         @assert ndims(ymask) == 0
-        t, res = mis_config_ad(gp.code, xst, ymask)
+        t, res = solution_ad(gp.code, xst, ymask)
         N = length(vertex_index)
         return fill(CountingTropical(asscalar(t).n, ConfigSampler(StaticBitVector(map(l->res[l], 1:N)))))
     end
 end
 
-function getconfigs_direct(gp::GraphProblem; all=false, usecuda=false)
+"""
+    solutions(problem, basetype; all=false, usecuda=false)
+    
+General routine to find solutions without bounding,
+
+* `basetype` can be a type with counting field,
+    * `CountingTropical{Float64,Float64}` for finding optimal solutions,
+    * `Polynomial{Float64, 1.0}` for enumerating all solutions,
+    * `Max2Poly{Float64, 1.0}` for optimal and suboptimal solutions.
+* When `all` is true, the program will use set for enumerate all possible solutions, otherwise, it will return one solution for each size.
+* `usecuda` can not be true if you want to use set to enumerate all possible solutions.
+"""
+function solutions(gp::GraphProblem, ::Type{BT}; all=false, usecuda=false) where BT
     if all && usecuda
         throw(ArgumentError("ConfigEnumerator can not be computed on GPU!"))
     end
-    T = (all ? bitstringset_type : bitstringsampler_type)(CountingTropical{Int64}, length(labels(gp.code)))
-    syms = labels(gp.code)
-    vertex_index = Dict([s=>i for (i, s) in enumerate(syms)])
-    xs = generate_tensors(l->onehotv(T, vertex_index[l]), gp)
-    if usecuda
-        xs = CuArray.(xs)
-    end
-    dynamic_einsum(gp.code, xs)
+    return contractf(fx_solutions(gp, BT, all), gp; usecuda=usecuda)
+end
+
+# return a mapping from label to variable `x`
+function fx_solutions(gp::GraphProblem, ::Type{BT}, all::Bool) where BT
+    T = (all ? bitstringset_type : bitstringsampler_type)(BT, length(labels(gp.code)))
+    vertex_index = Dict([s=>i for (i, s) in enumerate(labels(gp.code))])
+    return l->onehotv(T, vertex_index[l])
 end

src/graph_polynomials.jl

@@ -3,20 +3,47 @@ using OMEinsum: NestedEinsum, getixs, getiy
 using FFTW
 using LightGraphs
 
-export contractx, graph_polynomial, optimize_code
+export contractx, contractf, graph_polynomial, optimize_code
 export Independence, MaximalIndependence, Matching, Coloring
 const EinTypes = Union{EinCode,NestedEinsum}
 
 abstract type GraphProblem end
+"""
+    Independence{CT<:EinTypes} <: GraphProblem
+    Independence(graph; kwargs...)
+
+Independent set problem. For `kwargs`, check `optimize_code` API.
+"""
 struct Independence{CT<:EinTypes} <: GraphProblem
     code::CT
 end
+
+"""
+    Independence{CT<:EinTypes} <: GraphProblem
+    Independence(graph; kwargs...)
+
+Independent set problem. For `kwargs`, check `optimize_code` API.
+"""
 struct MaximalIndependence{CT<:EinTypes} <: GraphProblem
     code::CT
 end
+
+"""
+    Matching{CT<:EinTypes} <: GraphProblem
+    Matching(graph; kwargs...)
+
+Vertex matching problem. For `kwargs`, check `optimize_code` API.
+"""
 struct Matching{CT<:EinTypes} <: GraphProblem
     code::CT
 end
+
+"""
+    Coloring{K,CT<:EinTypes} <: GraphProblem
+    Coloring{K}(graph; kwargs...)
+
+K-Coloring problem. For `kwargs`, check `optimize_code` API.
+"""
 struct Coloring{K,CT<:EinTypes} <: GraphProblem
     code::CT
 end
@@ -39,7 +66,31 @@ function labels(code::EinTypes)
     return res
 end
 
-function graph_polynomial(gp::GraphProblem, ::Val{:fft}; usecuda=false, maxorder=graph_polynomial_maxorder(gp; usecuda=usecuda), r=1.0)
+"""
+    graph_polynomial(problem, method; usecuda=false, kwargs...)
+
+Computing the graph polynomial for specific problem.
+
+* `problem` can be one of the following instances,
+    * `Independence` for the independence polynomial,
+    * `MaximalIndependence` for the maximal independence polynomial,
+    * `Matching` for the matching polynomial,
+
+* `method` can be one of the following inputs,
+    * `Val(:finitefield)`, compute exactly with the finite field method.
+        It consumes additional kwargs [`max_iter`, `maxorder`], where `maxorder` is maximum order of polynomial
+        and `max_iter` is the maximum number of primes numbers to use in the finitefield algebra.
+        `max_iter` can be determined automatically in most cases.
+    * `Val(:polynomial)`, compute directly with `Polynomial` number type,
+    * `Val(:fft)`, compute with the fast fourier transformation approach, fast but needs to tune the hyperparameter `r`.
+        It Consumes additional kwargs [`maxorder`, `r`]. The larger `r` is,
+        the more accurate the factors of high order terms, and the less accurate the factors of low order terms.
+    * `Val(:fitting)`, compute with the polynomial fitting approach, fast but inaccurate for large graphs.
+"""
+function graph_polynomial end
+
+function graph_polynomial(gp::GraphProblem, ::Val{:fft}; usecuda=false, 
+        maxorder=graph_polynomial_maxorder(gp; usecuda=usecuda), r=1.0)
 	ω = exp(-2im*π/(maxorder+1))
 	xs = r .* collect(ω .^ (0:maxorder))
 	ys = [asscalar(contractx(gp, x; usecuda=usecuda)) for x in xs]
@@ -68,7 +119,8 @@ function _polynomial_single(gp::GraphProblem, ::Type{T}; usecuda, maxorder) wher
 	A \ T.(ys)
 end
 
-function graph_polynomial(gp::GraphProblem, ::Val{:finitefield}; usecuda=false, maxorder=graph_polynomial_maxorder(gp; usecuda=usecuda), max_iter=100)
+function graph_polynomial(gp::GraphProblem, ::Val{:finitefield}; usecuda=false,
+        maxorder=graph_polynomial_maxorder(gp; usecuda=usecuda), max_iter=100)
     TI = Int32  # Int 32 is faster
     N = typemax(TI)
     YS = []
@@ -93,21 +145,38 @@ function improved_counting(sequences)
     map(yi->Mods.CRT(yi...), zip(sequences...))
 end
 
-function optimize_code(code; method=:kahypar, sc_target=17, max_group_size=40, nrepeat=10, imbalances=0.0:0.001:0.8)
+"""
+    optimize_code(code; optmethod=:kahypar, sc_target=17, max_group_size=40, nrepeat=10, imbalances=0.0:0.001:0.8)
+
+Optimize the contraction order.
+
+* `optmethod` can be one of
+    * `:kahypar`, the kahypar + greedy approach, takes kwargs [`sc_target`, `max_group_size`, `imbalances`].
+    Check `optimize_kahypar` method in package `OMEinsumContractionOrders`.
+    * `:auto`, also the kahypar + greedy approach, but determines `sc_target` automatically. It is slower!
+    * `:greedy`, the greedy approach. Check in `optimize_greedy` in package `OMEinsum`.
+    * `:raw`, do nothing and return the raw EinCode.
+"""
+function optimize_code(code::EinTypes; optmethod=:kahypar, sc_target=17, max_group_size=40, nrepeat=10, imbalances=0.0:0.001:0.8)
     size_dict = Dict([s=>2 for s in labels(code)])
-    optcode = if method == :kahypar
+    optcode = if optmethod == :kahypar
         optimize_kahypar(code, size_dict; sc_target=sc_target, max_group_size=max_group_size, imbalances=imbalances)
-    elseif method == :greedy
+    elseif optmethod == :greedy
         optimize_greedy(code, size_dict; nrepeat=nrepeat)
+    elseif optmethod == :auto
+        optimize_kahypar_auto(code, size_dict; max_group_size=max_group_size, effort=500)
+    elseif optmethod == :raw
+        code
     else
-        ArgumentError("optimizer `$method` not defined.")
+        ArgumentError("optimizer `$optmethod` not defined.")
     end
     println("time/space complexity is $(OMEinsum.timespace_complexity(optcode, size_dict))")
     return optcode
 end
 
-function contractx(gp::GraphProblem, x::T; usecuda=false) where T
-    xs = generate_tensors(_->x, gp)
+contractx(gp::GraphProblem, x; usecuda=false) = contractf(_->x, gp; usecuda=usecuda)
+function contractf(f, gp::GraphProblem; usecuda=false)
+    xs = generate_tensors(f, gp)
     if usecuda
         xs = CuArray.(xs)
     end
@@ -116,9 +185,10 @@ end
 
 ############### Problem specific implementations ################
 ### independent set ###
-function Independence(g::SimpleGraph)
-    Independence(EinCode(([(i,) for i in LightGraphs.vertices(g)]..., # labels for vertex tensors
-                    [minmax(e.src,e.dst) for e in LightGraphs.edges(g)]...), ()))  # labels for edge tensors
+function Independence(g::SimpleGraph; kwargs...)
+    rawcode = EinCode(([(i,) for i in LightGraphs.vertices(g)]..., # labels for vertex tensors
+                    [minmax(e.src,e.dst) for e in LightGraphs.edges(g)]...), ())  # labels for edge tensors
+    Independence(optimize_code(rawcode; kwargs...))
 end
 
 function generate_tensors(fx, gp::Independence)
@@ -136,7 +206,7 @@ misv(val::T) where T = [one(T), val]
 graph_polynomial_maxorder(gp::Independence; usecuda) = Int(sum(contractx(gp, TropicalF64(1.0); usecuda=usecuda)).n)
 
 ### coloring ###
-Coloring(g::SimpleGraph) = Coloring(Independence(g).code)
+Coloring(g::SimpleGraph; kwargs...) = Coloring(Independence(g; kwargs...).code)
 function generate_tensors(fx, c::Coloring{K}) where K
     ixs = getixs(flatten(code))
     T = eltype(fx(ixs[1][1]))
@@ -159,9 +229,10 @@ end
 coloringv(vals::Vector{T}) where T = vals
 
 ### matching ###
-function Matching(g::SimpleGraph)
-    Matching(EinCode(([(minmax(e.src,e.dst),) for e in LightGraphs.edges(g)]..., # labels for edge tensors
-                    [([minmax(i,j) for j in neighbors(g, i)]...,) for i in LightGraphs.vertices(g)]...,), ()))        # labels for vertex tensors
+function Matching(g::SimpleGraph; kwargs...)
+    rawcode = EinCode(([(minmax(e.src,e.dst),) for e in LightGraphs.edges(g)]..., # labels for edge tensors
+                    [([minmax(i,j) for j in neighbors(g, i)]...,) for i in LightGraphs.vertices(g)]...,), ())       # labels for vertex tensors
+    Matching(optimize_code(rawcode; kwargs...))
 end
 
 function generate_tensors(fx, m::Matching)
@@ -193,8 +264,9 @@ end
 graph_polynomial_maxorder(m::Matching; usecuda) = Int(sum(contractx(m, TropicalF64(1.0); usecuda=usecuda)).n)
 
 ### maximal independent set ###
-function MaximalIndependence(g::SimpleGraph)
-    MaximalIndependence(EinCode(([(LightGraphs.neighbors(g, v)..., v) for v in LightGraphs.vertices(g)]...,), ()))
+function MaximalIndependence(g::SimpleGraph; kwargs...)
+    rawcode = EinCode(([(LightGraphs.neighbors(g, v)..., v) for v in LightGraphs.vertices(g)]...,), ())
+    MaximalIndependence(optimize_code(rawcode; kwargs...))
 end
 
 function generate_tensors(fx, mi::MaximalIndependence)
@@ -213,8 +285,4 @@ function neighbortensor(x::T, d::Int) where T
     return t
 end
 
-graph_polynomial_maxorder(mi::MaximalIndependence; usecuda) = Int(sum(contractx(mi, TropicalF64(1.0); usecuda=usecuda)).n)
-
-for GP in [:Independence, :MaximalIndependence, :Matching, :Coloring]
-    @eval optimize_code(gp::$GP; kwargs...) = $GP(optimize_code(gp.code; kwargs...))
-end
+graph_polynomial_maxorder(mi::MaximalIndependence; usecuda) = Int(sum(contractx(mi, TropicalF64(1.0); usecuda=usecuda)).n)