Skip to content

Commit

Permalink
Refactor Interfaces (#81)
Browse files Browse the repository at this point in the history
* update

* save

* NoWeight->UnitWeight

* update

* update

* update

* update

* update

* update docs

* update

* fix CUDA test

* fix document

* fix document

* fix doc build
  • Loading branch information
GiggleLiu authored Apr 7, 2024
1 parent 43e64c2 commit 0d23e37
Show file tree
Hide file tree
Showing 63 changed files with 783 additions and 1,092 deletions.
28 changes: 28 additions & 0 deletions Makefile
Original file line number Diff line number Diff line change
@@ -0,0 +1,28 @@
JL = julia --project

default: init test

init:
$(JL) -e 'using Pkg; Pkg.precompile()'
init-docs:
$(JL) -e 'using Pkg; Pkg.activate("docs"); Pkg.develop(path="."), Pkg.precompile()'

update:
$(JL) -e 'using Pkg; Pkg.update(); Pkg.precompile()'
update-docs:
$(JL) -e 'using Pkg; Pkg.activate("docs"); Pkg.update(); Pkg.precompile()'

test:
$(JL) -e 'using Pkg; Pkg.test("GenericTensorNetworks")'

coverage:
$(JL) -e 'using Pkg; Pkg.test("GenericTensorNetworks"; coverage=true)'

serve:
$(JL) -e 'using Pkg; Pkg.activate("docs"); using LiveServer; servedocs(;skip_dirs=["docs/src/assets", "docs/src/generated"], literate_dir="examples")'

clean:
rm -rf docs/build
find . -name "*.cov" -type f -print0 | xargs -0 /bin/rm -f

.PHONY: init test coverage serve clean init-docs update update-docs
1 change: 0 additions & 1 deletion docs/make.jl
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,6 @@ makedocs(;
"Satisfiability problem" => "generated/Satisfiability.md",
"Set covering problem" => "generated/SetCovering.md",
"Set packing problem" => "generated/SetPacking.md",
"Hyper Spin glass problem" => "generated/HyperSpinGlass.md",
#"Other problems" => "generated/Others.md",
],
"Topics" => [
Expand Down
25 changes: 11 additions & 14 deletions docs/src/index.md
Original file line number Diff line number Diff line change
Expand Up @@ -23,32 +23,29 @@ If you find our paper or software useful in your work, we would be grateful if y
You can find a set up guide in our [README](https://github.com/QuEraComputing/GenericTensorNetworks.jl).
To get started, open a Julia REPL and type the following code.

```julia
julia> using GenericTensorNetworks, Graphs

julia> # using CUDA

julia> solve(
IndependentSet(
Graphs.random_regular_graph(20, 3);
```@repl
using GenericTensorNetworks, Graphs#, CUDA
solve(
GenericTensorNetwork(IndependentSet(
Graphs.random_regular_graph(20, 3),
UnitWeight()); # default: uniform weight 1
optimizer = TreeSA(),
weights = NoWeight(), # default: uniform weight 1
openvertices = (), # default: no open vertices
fixedvertices = Dict() # default: no fixed vertices
),
GraphPolynomial();
usecuda=false # the default value
)
0-dimensional Array{Polynomial{BigInt, :x}, 0}:
Polynomial(1 + 20*x + 160*x^2 + 659*x^3 + 1500*x^4 + 1883*x^5 + 1223*x^6 + 347*x^7 + 25*x^8)
```

Here the main function [`solve`](@ref) takes three input arguments, the problem instance of type [`IndependentSet`](@ref), the property instance of type [`GraphPolynomial`](@ref) and an optional key word argument `usecuda` to decide use GPU or not.
If one wants to use GPU to accelerate the computation, the `using CUDA` statement must uncommented.
If one wants to use GPU to accelerate the computation, the `, CUDA` should be uncommented.

An [`IndependentSet`](@ref) instance takes two positional arguments to initialize, the graph instance that one wants to solve and the weights for each vertex. Here, we use a random regular graph with 20 vertices and degree 3, and the default uniform weight 1.

The problem instance takes four arguments to initialize, the only positional argument is the graph instance that one wants to solve, the key word argument `optimizer` is for specifying the tensor network optimization algorithm, the key word argument `weights` is for specifying the weights of vertices as either a vector or `NoWeight()`.
The [`GenericTensorNetwork`](@ref) function is a constructor for the problem instance, which takes the problem instance as the first argument and optional key word arguments. The key word argument `optimizer` is for specifying the tensor network optimization algorithm.
The keyword argument `openvertices` is a tuple of labels for specifying the degrees of freedom not summed over, and `fixedvertices` is a label-value dictionary for specifying the fixed values of the degree of freedoms.
Here, we use [`TreeSA`](@ref) method as the tensor network optimizer, and leave `weights` and `openvertices` the default values.
Here, we use [`TreeSA`](@ref) method as the tensor network optimizer, and leave `openvertices` the default values.
The [`TreeSA`](@ref) method finds the best contraction order in most of our applications, while the default [`GreedyMethod`](@ref) runs the fastest.

The first execution of this function will be a bit slow due to Julia's just in time compiling.
Expand Down
80 changes: 25 additions & 55 deletions docs/src/performancetips.md
Original file line number Diff line number Diff line change
Expand Up @@ -2,17 +2,15 @@
## Optimize contraction orders

Let us use the independent set problem on 3-regular graphs as an example.
```julia
julia> using GenericTensorNetworks, Graphs, Random

julia> graph = random_regular_graph(120, 3)
{120, 180} undirected simple Int64 graph

julia> problem = IndependentSet(graph; optimizer=TreeSA(
sc_target=20, sc_weight=1.0, rw_weight=3.0, ntrials=10, βs=0.01:0.1:15.0, niters=20), simplifier=MergeGreedy());
```@repl performancetips
using GenericTensorNetworks, Graphs, Random
graph = random_regular_graph(120, 3)
iset = IndependentSet(graph)
problem = GenericTensorNetwork(iset; optimizer=TreeSA(
sc_target=20, sc_weight=1.0, rw_weight=3.0, ntrials=10, βs=0.01:0.1:15.0, niters=20));
```

The [`IndependentSet`](@ref) constructor maps an independent set problem to a tensor network with optimized contraction order.
The [`GenericTensorNetwork`](@ref) constructor maps an independent set problem to a tensor network with optimized contraction order.
The key word argument `optimizer` specifies the contraction order optimizer of the tensor network.
Here, we choose the local search based [`TreeSA`](@ref) algorithm, which often finds the smallest time/space complexity and supports slicing.
One can type `?TreeSA` in a Julia REPL for more information about how to configure the hyper-parameters of the [`TreeSA`](@ref) method,
Expand All @@ -22,46 +20,32 @@ Alternative tensor network contraction order optimizers include
* [`KaHyParBipartite`](@ref)
* [`SABipartite`](@ref)

The keyword argument `simplifier` specifies the preprocessor to improve the searching speed of the contraction order finding.
For example, the `MergeGreedy()` here "contracts" tensors greedily whenever the contraction result has a smaller space complexity.
It can remove all vertex tensors (vectors) before entering the contraction order optimization algorithm.

The returned object `problem` contains a field `code` that specifies the tensor network with optimized contraction order.
For an independent set problem, the optimal contraction time/space complexity is ``\sim 2^{{\rm tw}(G)}``, where ``{\rm tw(G)}`` is the [tree-width](https://en.wikipedia.org/wiki/Treewidth) of ``G``.
One can check the time, space and read-write complexity with the [`contraction_complexity`](@ref) function.

```julia
julia> contraction_complexity(problem)
Time complexity (number of element-wise multiplications) = 2^20.568850503058382
Space complexity (number of elements in the largest intermediate tensor) = 2^16.0
Read-write complexity (number of element-wise read and write) = 2^18.70474460304404
```@repl performancetips
contraction_complexity(problem)
```

The return values are `log2` values of the number of multiplications, the number elements in the largest tensor during contraction and the number of read-write operations to tensor elements.
In this example, the number `*` operations is ``\sim 2^{21.9}``, the number of read-write operations are ``\sim 2^{20}``, and the largest tensor size is ``2^{17}``.
One can check the element size by typing
```julia
julia> sizeof(TropicalF64)
8

julia> sizeof(TropicalF32)
4

julia> sizeof(StaticBitVector{200,4})
32

julia> sizeof(TruncatedPoly{5,Float64,Float64})
48
```@repl performancetips
sizeof(TropicalF64)
sizeof(TropicalF32)
sizeof(StaticBitVector{200,4})
sizeof(TruncatedPoly{5,Float64,Float64})
```

One can use [`estimate_memory`](@ref) to get a good estimation of peak memory in bytes.
For example, to compute the graph polynomial, the peak memory can be estimated as follows.
```julia
julia> estimate_memory(problem, GraphPolynomial(; method=:finitefield))
297616

julia> estimate_memory(problem, GraphPolynomial(; method=:polynomial))
71427840
```@repl performancetips
estimate_memory(problem, GraphPolynomial(; method=:finitefield))
estimate_memory(problem, GraphPolynomial(; method=:polynomial))
```
The finite field approach only requires 298 KB memory, while using the [`Polynomial`](https://juliamath.github.io/Polynomials.jl/stable/polynomials/polynomial/#Polynomial-2) number type requires 71 MB memory.

Expand All @@ -75,25 +59,11 @@ For large scale applications, it is also possible to slice over certain degrees
loop and accumulate over certain degrees of freedom so that one can have a smaller tensor network inside the loop due to the removal of these degrees of freedom.
In the [`TreeSA`](@ref) optimizer, one can set `nslices` to a value larger than zero to turn on this feature.

```julia
julia> using GenericTensorNetworks, Graphs, Random

julia> graph = random_regular_graph(120, 3)
{120, 180} undirected simple Int64 graph

julia> problem = IndependentSet(graph; optimizer=TreeSA(βs=0.01:0.1:25.0, ntrials=10, niters=10));

julia> contraction_complexity(problem)
Time complexity (number of element-wise multiplications) = 2^20.53277253647484
Space complexity (number of elements in the largest intermediate tensor) = 2^16.0
Read-write complexity (number of element-wise read and write) = 2^19.34699193791874

julia> problem = IndependentSet(graph; optimizer=TreeSA(βs=0.01:0.1:25.0, ntrials=10, niters=10, nslices=5));

julia> contraction_complexity(problem)
Time complexity (number of element-wise multiplications) = 2^21.117277836449578
Space complexity (number of elements in the largest intermediate tensor) = 2^11.0
Read-write complexity (number of element-wise read and write) = 2^19.854965754099602
```@repl performancetips
problem = GenericTensorNetwork(iset; optimizer=TreeSA(βs=0.01:0.1:25.0, ntrials=10, niters=10));
contraction_complexity(problem)
problem = GenericTensorNetwork(iset; optimizer=TreeSA(βs=0.01:0.1:25.0, ntrials=10, niters=10, nslices=5));
contraction_complexity(problem)
```

In the second `IndependentSet` constructor, we slice over 5 degrees of freedom, which can reduce the space complexity by at most 5.
Expand All @@ -103,7 +73,7 @@ i.e. the peak memory usage is reduced by a factor ``32``, while the (theoretical
## GEMM for Tropical numbers
One can speed up the Tropical number matrix multiplication when computing the solution space property [`SizeMax`](@ref)`()` by using the Tropical GEMM routines implemented in package [`TropicalGEMM`](https://github.com/TensorBFS/TropicalGEMM.jl/).

```julia
```julia-repl
julia> using BenchmarkTools
julia> @btime solve(problem, SizeMax())
Expand Down Expand Up @@ -139,7 +109,7 @@ results = multiprocess_run(collect(1:10)) do seed
n = 10
@info "Graph size $n x $n, seed= $seed"
g = random_diagonal_coupled_graph(n, n, 0.8)
gp = Independence(g; optimizer=TreeSA(), simplifier=MergeGreedy())
gp = GenericTensorNetwork(IndependentSet(g); optimizer=TreeSA())
res = solve(gp, GraphPolynomial())[]
return res
end
Expand All @@ -165,7 +135,7 @@ You will see a vector of polynomials printed out.
## Make use of GPUs
To upload the computation to GPU, you just add `using CUDA` before calling the `solve` function, and set the keyword argument `usecuda` to `true`.
```julia
```julia-repl
julia> using CUDA
[ Info: OMEinsum loaded the CUDA module successfully

Expand Down
7 changes: 2 additions & 5 deletions docs/src/ref.md
Original file line number Diff line number Diff line change
Expand Up @@ -2,14 +2,14 @@
## Graph problems
```@docs
solve
GenericTensorNetwork
GraphProblem
IndependentSet
MaximalIS
Matching
Coloring
DominatingSet
SpinGlass
HyperSpinGlass
MaxCut
PaintShop
Satisfiability
Expand All @@ -27,14 +27,12 @@ Interfaces [`GenericTensorNetworks.generate_tensors`](@ref), [`labels`](@ref), [
```@docs
GenericTensorNetworks.generate_tensors
labels
terms
energy_terms
flavors
get_weights
chweights
nflavor
fixedvertices
extract_result
```

#### Graph Problem Utilities
Expand All @@ -49,7 +47,6 @@ is_set_packing
cut_size
spinglass_energy
hyperspinglass_energy
num_paint_shop_color_switch
paint_shop_coloring_from_config
mis_compactify!
Expand Down
40 changes: 9 additions & 31 deletions docs/src/sumproduct.md
Original file line number Diff line number Diff line change
Expand Up @@ -3,45 +3,23 @@
It is a sum-product expression tree to store [`ConfigEnumerator`](@ref) in a lazy style, where configurations can be extracted by depth first searching the tree with the `Base.collect` method.
Although it is space efficient, it is in general not easy to extract information from it due to the exponential large configuration space.
Directed sampling is one of its most important operations, with which one can get some statistic properties from it with an intermediate effort. For example, if we want to check some property of an intermediate scale graph, one can type
```julia
julia> graph = random_regular_graph(70, 3)

julia> problem = IndependentSet(graph; optimizer=TreeSA());

julia> tree = solve(problem, ConfigsAll(; tree_storage=true))[];
16633909006371
```@repl sumproduct
graph = random_regular_graph(70, 3)
problem = GenericTensorNetwork(IndependentSet(graph); optimizer=TreeSA());
tree = solve(problem, ConfigsAll(; tree_storage=true))[]
```
If one wants to store these configurations, he will need a hard disk of size 256 TB!
However, this sum-product binary tree structure supports efficient and unbiased direct sampling.

```julia
samples = generate_samples(tree, 1000);
```@repl sumproduct
samples = generate_samples(tree, 1000)
```

With these samples, one can already compute useful properties like Hamming distance (see [`hamming_distribution`](@ref)) distribution.

```julia
julia> using UnicodePlots

julia> lineplot(hamming_distribution(samples, samples))
┌────────────────────────────────────────┐
100000 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢰⠹⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
│⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡎⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
│⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⡇⠀⢣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
│⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
│⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸⠀⠀⠸⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
│⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
│⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⣇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
│⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢰⠃⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
│⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸⠀⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
│⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡞⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
│⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⣇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
│⠀⠀⠀⠀⠀⠀⠀⠀⠀⢰⠁⠀⠀⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
│⠀⠀⠀⠀⠀⠀⠀⠀⠀⡼⠀⠀⠀⠀⠀⠀⠈⡆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
│⠀⠀⠀⠀⠀⠀⠀⠀⢠⠇⠀⠀⠀⠀⠀⠀⠀⢳⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
0 │⢀⣀⣀⣀⣀⣀⣀⣀⠎⠀⠀⠀⠀⠀⠀⠀⠀⠀⠓⢄⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⠀⠀⠀⠀│
└────────────────────────────────────────┘
⠀0⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀80⠀
```@repl sumproduct
using UnicodePlots
lineplot(hamming_distribution(samples, samples))
```

Here, the ``x``-axis is the Hamming distance and the ``y``-axis is the counting of pair-wise Hamming distances.
10 changes: 6 additions & 4 deletions examples/Coloring.jl
Original file line number Diff line number Diff line change
Expand Up @@ -19,9 +19,11 @@ locations = [[rot15(0.0, 2.0, i) for i=0:4]..., [rot15(0.0, 1.0, i) for i=0:4]..
show_graph(graph; locs=locations, format=:svg)

# ## Generic tensor network representation
#
# We construct the tensor network for the 3-coloring problem as
problem = Coloring{3}(graph);
# We can define the 3-coloring problem with the [`Coloring`](@ref) type as
coloring = Coloring{3}(graph)

# The tensor network representation of the 3-coloring problem can be obtained by
problem = GenericTensorNetwork(coloring)

# ### Theory (can skip)
# Type [`Coloring`](@ref) can be used for constructing the tensor network with optimized contraction order for a coloring problem.
Expand Down Expand Up @@ -67,7 +69,7 @@ show_graph(linegraph; locs=[0.5 .* (locations[e.src] .+ locations[e.dst])
# Let us construct the tensor network and see if there are solutions.
lineproblem = Coloring{3}(linegraph);

num_of_coloring = solve(lineproblem, CountingMax())[]
num_of_coloring = solve(GenericTensorNetwork(lineproblem), CountingMax())[]

# You will see the maximum size 28 is smaller than the number of edges in the `linegraph`,
# meaning no solution for the 3-coloring on edges of a Petersen graph.
6 changes: 5 additions & 1 deletion examples/DominatingSet.jl
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,11 @@ show_graph(graph; locs=locations, format=:svg)

# ## Generic tensor network representation
# We can use [`DominatingSet`](@ref) to construct the tensor network for solving the dominating set problem as
problem = DominatingSet(graph; optimizer=TreeSA());
dom = DominatingSet(graph)

# The tensor network representation of the dominating set problem can be obtained by
problem = GenericTensorNetwork(dom; optimizer=TreeSA())
# where the key word argument `optimizer` specifies the tensor network contraction order optimizer as a local search based optimizer [`TreeSA`](@ref).

# ### Theory (can skip)
# Let ``G=(V,E)`` be the target graph that we want to solve.
Expand Down
Loading

0 comments on commit 0d23e37

Please sign in to comment.