update docs

Author: GiggleLiu

Project.toml

@@ -18,11 +18,13 @@ Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 Primes = "27ebfcd6-29c5-5fa9-bf4b-fb8fc14df3ae"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 SIMDTypes = "94e857df-77ce-4151-89e5-788b33177be4"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 TropicalNumbers = "b3a74e9c-7526-4576-a4eb-79c0d4c32334"
 Viznet = "52a3aca4-6234-47fd-b74a-806bdf78ede9"
 
 [compat]
 CUDA = "3.5"
+Cairo = "1.0"
 Compose = "0.9"
 FFTW = "1.4"
 Graphs = "1.4"
@@ -35,7 +37,7 @@ Requires = "1"
 SIMDTypes = "0.1"
 TropicalNumbers = "0.4, 0.5"
 Viznet = "0.3"
-Cairo = "1.0"
+StatsBase = "0.33"
 julia = "1"
 
 [extras]

README.md

@@ -58,7 +58,7 @@ time/space complexity is (7.965784284662086, 4.0)
 
 Here, the `problem` is a `Independence` instance, it contains the tensor network contraction tree for the target graph.
 Here, we choose the `TreeSA` optimizer to optimize the tensor network contraciton tree, it is a local search based algorithm, check [arXiv: 2108.05665](https://arxiv.org/abs/2108.05665). You will see some warnings, do not panic, this is because we set `sc_target` (target space complex) to 1 for agressive optimization of space complexity. Type `?TreeSA` in a Julia REPL for more information about the key word arguments.
-Similarly, one can select tensor network structures for solving other problems like `MaximalIndependence`, `MaxCut`, `Matching`, `Coloring{K}` and `set_packing`.
+Similarly, one can select tensor network structures for solving other problems like `MaximalIndependence`, `MaxCut`, `Matching`, `Coloring{K}`, `PaintShop` and `set_packing`.
 
 #### 1. find MIS size, count MISs and count ISs
 ```julia

docs/make.jl

@@ -1,6 +1,6 @@
 using Pkg
 using GraphTensorNetworks
-using GraphTensorNetworks: TropicalNumbers, Polynomials, Mods
+using GraphTensorNetworks: TropicalNumbers, Polynomials, Mods, OMEinsum, OMEinsumContractionOrders
 using Documenter
 using DocThemeIndigo
 using Literate
@@ -19,7 +19,7 @@ indigo = DocThemeIndigo.install(GraphTensorNetworks)
 DocMeta.setdocmeta!(GraphTensorNetworks, :DocTestSetup, :(using GraphTensorNetworks); recursive=true)
 
 makedocs(;
-    modules=[GraphTensorNetworks, TropicalNumbers, Polynomials, Mods],
+    modules=[GraphTensorNetworks, TropicalNumbers, Polynomials, Mods, OMEinsum, OMEinsumContractionOrders],
     authors="Jinguo Liu",
     repo="https://github.com/Happy-Diode/GraphTensorNetworks.jl/blob/{commit}{path}#{line}",
     sitename="GraphTensorNetworks.jl",

docs/src/ref.md

@@ -1,26 +1,86 @@
-```@meta
-DocTestSetup = quote
-    import TropicalNumbers, Mods, Polynomials
-end
-```
 # References
+## Graph problems
+```@docs
+solve
+Independence
+MaximalIndependence
+Matching
+Coloring
+MaxCut
+PaintShop
+```
 
-```@autodocs
-Modules = [GraphTensorNetworks]
+```@docs
+set_packing
 ```
 
+## Properties
 ```@docs
-TropicalNumbers.Tropical
+SizeMax
+CountingAll
+CountingMax
+GraphPolynomial
+SingleConfigMax
+ConfigsAll
+ConfigsMax
 ```
 
+## Element Algebras
 ```@docs
+TropicalNumbers.Tropical
 TropicalNumbers.CountingTropical
+Mods.Mod
+Polynomials.Polynomial
+TruncatedPoly
+Max2Poly
+ConfigEnumerator
+ConfigSampler
 ```
 
 ```@docs
-Mods.Mod
+StaticBitVector
+StaticElementVector
+save_configs
+load_configs
+@bv_str
+onehotv
+is_commutative_semiring
 ```
 
+## Tensor Network
 ```@docs
-Polynomials.Polynomial
+optimize_code
+timespace_complexity
+@ein_str
+GreedyMethod
+TreeSA
+SABipartite
+KaHyParBipartite
+MergeVectors
+MergeGreedy
+```
+
+## Others
+#### Graph
+```@docs
+is_independent_set
+mis_compactify!
+show_graph
+
+diagonal_coupled_graph
+square_lattice_graph
+unitdisk_graph
+```
+
+One can also use `random_regular_graph` and `smallgraph` in [Graphs](https://github.com/JuliaGraphs/Graphs.jl) to build special graphs.
+
+#### Lower level APIs
+```@docs
+best_solutions
+best2_solutions
+solutions
+all_solutions
+graph_polynomial
+max_size
+max_size_count
 ```

examples/IndependentSet/main.jl

@@ -1,11 +1,12 @@
 # # Independent set problem
 
-# ## Problem definition
+# ## Introduction
+using GraphTensorNetworks, Graphs
 
 # Please check the docstring of [`Independence`](@ref) for the definition of independence problem.
-# In the following, we are going to defined an independent set problem for the Petersen graph.
+@doc Independence
 
-using GraphTensorNetworks, Graphs
+# In the following, we are going to defined an independent set problem for the Petersen graph.
 
 graph = Graphs.smallgraph(:petersen)
 
@@ -18,13 +19,13 @@ show_graph(graph; locs=locations)
 
 # Let us contruct the problem instance with optimized tensor network contraction order as bellow.
 problem = Independence(graph; optimizer=TreeSA(sc_weight=1.0, ntrials=10,
-                         βs=0.01:0.1:15.0, niters=20, rw_weight=0.2));
+                         βs=0.01:0.1:15.0, niters=20, rw_weight=0.2),
+                         simplifier=MergeGreedy());
 
 
 # ## Solving properties
 
-# ### The decision problem
-# ##### maximum independent set size ``\alpha(G)``
+# ### Maximum independent set size ``\alpha(G)``
 maximum_independent_set_size = solve(problem, SizeMax())
 
 # ### Counting properties
@@ -34,12 +35,18 @@ count_all_independent_sets = solve(problem, CountingAll())
 # ##### counting independent sets with sizes ``\alpha(G)`` and ``\alpha(G)-1``
 count_max2_independent_sets = solve(problem, CountingMax(2))
 
-# ##### computing the independence polynomial
+# ##### independence polynomial
 # For the definition of independence polynomial, please check the docstring of [`Independence`](@ref) or this [wiki page](https://mathworld.wolfram.com/IndependencePolynomial.html).
+# There are 3 methods to compute a graph polynomial, `:finitefield`, `:fft` and `:polynomial`.
+# These methods are introduced in the docstring of [`GraphPolynomial`](@ref).
 independence_polynomial = solve(problem, GraphPolynomial(; method=:finitefield))
 
 # ### Configuration properties
 # ##### finding one maximum independent set (MIS)
+# There are two approaches to find one of the best solution.
+# The unbounded (default) version uses [`ConfigSampler`](@ref) to sample one of the best solutions directly.
+# The bounded version uses the binary gradient back-propagation (see our paper) to compute the gradients.
+# It requires caching intermediate states, but is often faster on CPU because it can use [`TropicalGEMM`](https://github.com/TensorBFS/TropicalGEMM.jl).
 max_config = solve(problem, SingleConfigMax(; bounded=false))[]
 
 # The return value contains a bit string, and one should read this bit string from left to right.
@@ -49,7 +56,9 @@ show_graph(graph; locs=locations, colors=[iszero(max_config.c.data[i]) ? "white"
                                  for i=1:nv(graph)])
 
 # ##### enumeration of all MISs
-all_max_configs = solve(problem, ConfigsMax(; bounded=true))[]
+# There are two approaches to enumerate all best-K solutions.
+# The bounded (default) version is always prefered because it can significantly use the memory usage.
+all_max_configs = solve(problem, ConfigsMax(1; bounded=true))[]
 
 using Compose
 

src/GraphTensorNetworks.jl

@@ -23,7 +23,7 @@ export bestk_solutions
 export contractx, contractf, graph_polynomial, max_size, max_size_count
 
 # Graphs
-export random_regular_graph, diagonal_coupled_graph, isindependentset
+export random_regular_graph, diagonal_coupled_graph, is_independent_set
 export square_lattice_graph, unitdisk_graph
 
 # Tensor Networks (Graph problems)
@@ -33,7 +33,6 @@ export Independence, MaximalIndependence, Matching, Coloring, optimize_code, set
 export solve, SizeMax, CountingAll, CountingMax, GraphPolynomial, SingleConfigMax, ConfigsAll, ConfigsMax
 
 # Utilities
-export is_independent_set
 export mis_compactify!
 export save_configs, load_configs
 

src/arithematics.jl

@@ -8,7 +8,27 @@ Base.abs(x::Mod) = x
 Base.isless(x::Mod{N}, y::Mod{N}) where N = mod(x.val, N) < mod(y.val, N)
 
 
-# this function is used for testing
+"""
+    is_commutative_semiring(a::T, b::T, c::T) where T
+
+Check if elements `a`, `b` and `c` satisfied the commutative semiring requirements.
+```math
+\\begin{align*}
+(a \\oplus b) \\oplus c = a \\oplus (b \\oplus c) & \\hspace{5em}\\triangleright\\text{commutative monoid \$\\oplus\$ with identity \$\\mathbb{0}\$}\\\\
+a \\oplus \\mathbb{0} = \\mathbb{0} \\oplus a = a &\\\\
+a \\oplus b = b \\oplus a &\\\\
+&\\\\
+(a \\odot b) \\odot c = a \\odot (b \\odot c)  &   \\hspace{5em}\\triangleright \\text{commutative monoid \$\\odot\$ with identity \$\\mathbb{1}\$}\\\\
+a \\odot  \\mathbb{1} =  \\mathbb{1} \\odot a = a &\\\\
+a \\odot b = b \\odot a &\\\\
+&\\\\
+a \\odot (b\\oplus c) = a\\odot b \\oplus a\\odot c  &  \\hspace{5em}\\triangleright \\text{left and right distributive}\\\\
+(a\\oplus b) \\odot c = a\\odot c \\oplus b\\odot c &\\\\
+&\\\\
+a \\odot \\mathbb{0} = \\mathbb{0} \\odot a = \\mathbb{0}
+\\end{align*}
+```
+"""
 function is_commutative_semiring(a::T, b::T, c::T) where T
     # +
     if (a + b) + c != a + (b + c)
@@ -70,6 +90,10 @@ struct TruncatedPoly{K,T,TO} <: Number
     coeffs::NTuple{K,T}
     maxorder::TO
 end
+
+"""
+    Max2Poly{T,TO} = TruncatedPoly{2,T,TO}
+"""
 const Max2Poly{T,TO} = TruncatedPoly{2,T,TO}
 Max2Poly(a, b, maxorder) = TruncatedPoly((a, b), maxorder)
 Max2Poly{T,TO}(a, b, maxorder) where {T,TO} = TruncatedPoly{2,T,TO}((a, b), maxorder)

src/bitvector.jl

@@ -54,13 +54,22 @@ Base.:(&)(x::StaticElementVector{N,S,C}, y::StaticElementVector{N,S,C}) where {N
 # difference of two element sets
 Base.:(⊻)(x::StaticElementVector{N,S,C}, y::StaticElementVector{N,S,C}) where {N,S,C} = StaticElementVector{N,S,C}(x.data .⊻ y.data)
 
+"""
+    onehotv(::Type{<:StaticElementVector}, i, v)
+    onehotv(::Type{<:StaticBitVector, i)
+
+Returns a static element vector, with the value at location `i` being `v` (1 if not specified).
+"""
 function onehotv(::Type{StaticElementVector{N,S,C}}, i, v) where {N,S,C}
     x = zeros(Int,N)
     x[i] = v
     return convert(StaticElementVector{N,S,C}, x)
 end
 
 ##### BitVectors
+"""
+    StaticBitVector{N,C} = StaticElementVector{N,1,C}
+"""
 const StaticBitVector{N,C} = StaticElementVector{N,1,C}
 @inline function Base.getindex(x::StaticBitVector{N,C}, i::Integer) where {N,C}
     @boundscheck (i <= N || throw(BoundsError(x, i)))  # NOTE: still checks bounds in global scope, why?
@@ -86,6 +95,7 @@ staticfalses(::Type{StaticBitVector{N,C}}) where {N,C} = zero(StaticBitVector{N,
 @generated function statictrues(::Type{StaticBitVector{N,C}}) where {N,C}
     Expr(:call, :(StaticBitVector{$N,$C}), Expr(:tuple, fill(typemax(UInt64), C)...))
 end
+
 onehotv(::Type{StaticBitVector{N,C}}, i, v) where {N,C} = v > 0 ? onehotv(StaticBitVector{N,C}, i) : zero(StaticBitVector{N,C})
 function onehotv(::Type{StaticBitVector{N,C}}, i) where {N,C}
     x = falses(N)

src/graphs.jl

@@ -1,21 +1,52 @@
-using Graphs, OMEinsumContractionOrders
+using Graphs, OMEinsumContractionOrders, StatsBase
 
-isindependentset(g, v) = !any(e->v[e.src] == 1 && v[e.dst] == 1, edges(g))
+"""
+    is_independent_set(g::SimpleGraph, vertices)
 
+Return true if `vertices` is an independent set of graph `g`.
+"""
+is_independent_set(g::SimpleGraph, v) = !any(e->v[e.src] == 1 && v[e.dst] == 1, edges(g))
+
+"""
+    square_lattice_graph(mask::AbstractMatrix{Bool})
+
+Create a masked square lattice graph.
+"""
 function square_lattice_graph(mask::AbstractMatrix{Bool})
     locs = [(i, j) for i=1:size(mask, 1), j=1:size(mask, 2) if mask[i,j]]
     unitdisk_graph(locs, 1.1)
 end
 
+"""
+    random_diagonal_coupled_graph(m::Int, n::Int, ρ::Real)
+
+Create a random masked diagonal coupled square lattice graph, with number of vertices fixed to ``\\lfloor mn\\rho \\rceil``.
+"""
 function random_diagonal_coupled_graph(m::Int, n::Int, ρ::Real)
-    diagonal_coupled_graph(rand(m, n) .< ρ)
+    @assert ρ >=0 && ρ <= 1
+    diagonal_coupled_graph(generate_mask(m, n, round(Int,m*n*ρ)))
 end
+function generate_mask(Nx::Int, Ny::Int, natoms::Int)
+    mask = zeros(Bool, Nx, Ny)
+    mask[sample(1:Nx*Ny, natoms; replace=false)] .= true
+    return mask
+end
+
+"""
+    diagonal_coupled_graph(mask::AbstractMatrix{Bool})
 
+Create a masked diagonal coupled square lattice graph from a specified `mask`.
+"""
 function diagonal_coupled_graph(mask::AbstractMatrix{Bool})
     locs = [(i, j) for i=1:size(mask, 1), j=1:size(mask, 2) if mask[i,j]]
     unitdisk_graph(locs, 1.5)
 end
 
+"""
+    unitdisk_graph(locs::AbstractVector, unit::Real)
+
+Create a unit disk graph with locations specified by `locs` and unit distance `unit`.
+"""
 function unitdisk_graph(locs::AbstractVector, unit::Real)
     n = length(locs)
     g = SimpleGraph(n)

src/interfaces.jl

@@ -180,6 +180,20 @@ function solve(gp::GraphProblem, property::AbstractProperty; T=Float64, usecuda=
 end
 _has_weight(gp::GraphProblem) = hasfield(typeof(gp), :weights) && gp.weights isa Vector # ugly but makes life easier
 
+"""
+    max_size(problem; usecuda=false)
+
+Returns the maximum size of the graph problem. 
+A shorthand of `solve(problem, SizeMax(); usecuda=false)`.
+"""
+function max_size end
+"""
+    max_size_count(problem; usecuda=false)
+
+Returns the maximum size and the counting of the graph problem.
+It is a shorthand of `solve(problem, CountingMax(); usecuda=false)`.
+"""
+function max_size_count end
 for TP in [:MaximalIndependence, :Independence, :Matching, :MaxCut, :PaintShop]
     @eval max_size(m::$TP; usecuda=false) = Int(sum(solve(m, SizeMax(); usecuda=usecuda)).n)  # floating point number is faster (BLAS)
     @eval max_size_count(m::$TP; usecuda=false) = (r = sum(solve(m, CountingMax(); usecuda=usecuda)); (Int(r.n), Int(r.c)))