Hyper Spin-glass model and Partition function (#55)

* hyper spin glass

* new hyper spin glass model doc

* fix doc make

Author: GiggleLiu

Project.toml

@@ -8,6 +8,7 @@ AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
+DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -25,8 +26,9 @@ StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 TropicalNumbers = "b3a74e9c-7526-4576-a4eb-79c0d4c32334"
 
 [compat]
-AbstractTrees = "0.4"
+AbstractTrees = "0.3, 0.4"
 CUDA = "3"
+DocStringExtensions = "0.8, 0.9"
 FFTW = "1.4"
 Graphs = "1.7"
 LuxorGraphPlot = "0.1"

docs/make.jl

@@ -57,6 +57,7 @@ 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" => [

docs/serve.jl

@@ -9,10 +9,10 @@ function serve(;host::String="0.0.0.0", port::Int=8000)
     LiveServer.servedocs(;
         doc_env=true,
         skip_dirs=[
-            joinpath("docs", "src", "assets"),
-            joinpath("docs", "src", "generated"),
-            joinpath("docs", "src", "notebooks"),
+            #joinpath("docs", "src", "generated"),
+            joinpath("docs", "src"),
             joinpath("docs", "build"),
+            joinpath("docs", "Manifest.toml"),
         ],
         literate="examples",
         host=\"$host\",

docs/src/ref.md

@@ -9,6 +9,7 @@ Matching
 Coloring
 DominatingSet
 SpinGlass
+HyperSpinGlass
 MaxCut
 PaintShop
 Satisfiability
@@ -46,6 +47,7 @@ is_set_packing
 
 cut_size
 spinglass_energy
+hyperspinglass_energy
 num_paint_shop_color_switch
 paint_shop_coloring_from_config
 mis_compactify!
@@ -65,6 +67,7 @@ print_mining
 
 ## Properties
 ```@docs
+PartitionFunction
 SizeMax
 SizeMin
 CountingAll
@@ -117,6 +120,7 @@ hamming_distribution
 optimize_code
 getixsv
 getiyv
+contraction_complexity
 timespace_complexity
 timespacereadwrite_complexity
 estimate_memory

examples/HyperSpinGlass.jl

@@ -0,0 +1,95 @@
+# # Hyper-Spin-glass problem
+
+# !!! note
+#     It is highly recommended to read the [Independent set problem](@ref) chapter before reading this one.
+
+# ## Problem definition
+# The hyper-spin-glass problem is a natural generalization of the spin-glass problem from simple graph to hypergraph.
+# A hyper-spin-glass problem of hypergraph ``H = (V, E)`` can be characterized by the following energy function
+# ```math
+# E = \sum_{c \in E} w_{c} \prod_{v\in c} S_v
+# ```
+# where ``S_v \in \{-1, 1\}``, ``w_c`` is coupling strength associated with hyperedge ``c``.
+# In the program, we use boolean variable ``s_v = \frac{1-S_v}{2}`` to represent a spin configuration.
+
+using GenericTensorNetworks
+
+# In the following, we are going to defined an spin glass problem for the following hypergraph.
+num_vertices = 15
+
+hyperedges = [[1,3,4,6,7], [4,7,8,12], [2,5,9,11,13],
+    [1,2,14,15], [3,6,10,12,14], [8,14,15], 
+    [1,2,6,11], [1,2,4,6,8,12]]
+
+weights = [-1, 1, -1, 1, -1, 1, -1, 1];
+
+# The energy function is
+# ```math
+# \begin{align*}
+# E = &-s_1s_3s_4s_6s_7 + s_4s_7s_8s_{12} - s_2s_5s_9s_{11}s_{13} +\\
+#    &s_1s_2s_{14}s_{15} - s_3s_6s_{10}s_{12}s_{14} + s_8s_{14}s_{15} +\\
+#    &s_1s_2s_6s_{11} - s_1s_2s_4s_6s_8s_{12}
+# \end{align*}
+# ```
+
+# ## Generic tensor network representation
+# We define an anti-ferromagnetic spin glass problem as
+problem = HyperSpinGlass(num_vertices, hyperedges; weights);
+
+# ### Theory (can skip)
+# Let ``H = (V, E)`` be a hypergraph. The tensor network for the partition function of the energy model for ``H`` can be defined as the following tiple of (alphabet of labels, input tensors, output labels).
+# ```math
+# \begin{cases}
+# \Lambda &= \{s_v \mid v \in V\}\\
+# \mathcal{T} &= \{B^{(c)}_{s_{N(c, 1),N(c, 2),\ldots,N(c, d(c))}} \mid c \in E\} \cup \{W^{(v)}_{s_v} \mid v \in V\}\\
+# \sigma_o &= \varepsilon
+# \end{cases}
+# ```
+# where ``s_v \in \{0, 1\}`` is the boolean degreen associated to vertex ``v``,
+# ``N(c, k)`` is the ``k``th vertex of hyperedge ``c``, and ``d(c)`` is the degree of ``c``.
+# The edge tensor ``B^{(c)}`` is defined as
+# ```math
+# B^{(c)} = \begin{cases}
+# x^{w_c} &  (\sum_{v\in c} s_v) \;is\; even, \\
+# x^{-w_c}  &  otherwise.
+# \end{cases}
+# ```
+# and the vertex tensor ``W^{(v)}`` (used to carry labels) is defined as
+# ```math
+# W^{(v)} = \left(\begin{matrix}1_v\\ 1_v\end{matrix}\right)
+# ```
+
+# ## Solving properties
+# ### Minimum and maximum energies
+# Its ground state energy is -8.
+Emin = solve(problem, SizeMin())[]
+# While the state correspond to the highest energy has the ferromagnetic order.
+Emax = solve(problem, SizeMax())[]
+
+# In this example, the spin configurations can be choosen to make all hyperedges having even or odd spin parity.
+
+# ### Counting properties
+# ##### partition function and graph polynomial
+# The graph polynomial defined for the hyper-spin-glass problem is a Laurent polynomial
+# ```math
+# Z(G, w, x) = \sum_{k=E_{\rm min}}^{E_{\rm max}} c_k x^k,
+# ```
+# where ``E_{\rm min}`` and ``E_{\rm max}`` are minimum and maximum energies,
+# ``c_k`` is the number of spin configurations with energy ``k``.
+# Let the inverse temperature ``\beta = 2``, the partition function is
+β = 2.0
+Z = solve(problem, PartitionFunction(β))[]
+
+# The infinite temperature partition function is the counting of all feasible configurations
+solve(problem, PartitionFunction(0.0))[]
+
+# Let ``x = e^\beta``, it corresponds to the partition function of a spin glass at temperature ``\beta^{-1}``.
+poly = solve(problem, GraphPolynomial())[]
+
+# ### Configuration properties
+# ##### finding a ground state
+ground_state = solve(problem, SingleConfigMin())[].c.data
+
+Emin_verify = hyperspinglass_energy(hyperedges, ground_state; weights)
+
+# You should see a consistent result as above `Emin`.

examples/IndependentSet.jl

@@ -71,7 +71,8 @@ maximum_independent_set_size.n
 # We can count all independent sets with the [`CountingAll`](@ref) property.
 count_all_independent_sets = solve(problem, CountingAll())[]
 
-# The return value has type `Float64`.
+# The return value has type `Float64`. The counting of all independent sets is equivalent to the infinite temperature partition function
+solve(problem, PartitionFunction(0.0))[]
 
 # We can count the maximum independent sets with [`CountingMax`](@ref).
 count_maximum_independent_sets = solve(problem, CountingMax())[]

src/GenericTensorNetworks.jl

@@ -11,6 +11,7 @@ import Polynomials
 using Polynomials: Polynomial, LaurentPolynomial, printpoly, fit
 using FFTW
 using Mods, Primes
+using DocStringExtensions
 using Base.Cartesian
 import AbstractTrees: children, printnode, print_tree
 import StatsBase
@@ -42,6 +43,7 @@ export IndependentSet, mis_compactify!, is_independent_set
 export MaximalIS, is_maximal_independent_set
 export cut_size, MaxCut
 export spinglass_energy, SpinGlass
+export hyperspinglass_energy, HyperSpinGlass
 export PaintShop, paintshop_from_pairs, num_paint_shop_color_switch, paint_shop_coloring_from_config
 export Coloring, is_vertex_coloring
 export Satisfiability, CNF, CNFClause, BoolVar, satisfiable, @bools, ∨, ¬, ∧
@@ -52,7 +54,7 @@ export SetCovering, is_set_covering
 export OpenPitMining, is_valid_mining, print_mining
 
 # Interfaces
-export solve, SizeMax, SizeMin, CountingAll, CountingMax, CountingMin, GraphPolynomial, SingleConfigMax, SingleConfigMin, ConfigsAll, ConfigsMax, ConfigsMin, Single
+export solve, SizeMax, SizeMin, PartitionFunction, CountingAll, CountingMax, CountingMin, GraphPolynomial, SingleConfigMax, SingleConfigMin, ConfigsAll, ConfigsMax, ConfigsMin, Single
 
 # Utilities
 export save_configs, load_configs, hamming_distribution, save_sumproduct, load_sumproduct

src/arithematics.jl

@@ -198,7 +198,7 @@ This algebra maps
 
 Fields
 ------------------------
-* `orders` is a vector of [`Tropical`](@ref) ([`CoutingTropical`](@ref)) numbers as the largest-K solution sizes (solutions).
+* `orders` is a vector of [`Tropical`](@ref) ([`CountingTropical`](@ref)) numbers as the largest-K solution sizes (solutions).
 
 Examples
 ------------------------------

src/interfaces.jl

@@ -44,6 +44,18 @@ Counting the total number of sets. e.g. for the [`IndependentSet`](@ref) problem
 """
 struct CountingAll <: AbstractProperty end
 
+"""
+$TYPEDEF
+$FIELDS
+
+Compute the partition function for the target problem.
+
+* The corresponding tensor element type is `T`.
+"""
+struct PartitionFunction{T} <: AbstractProperty
+    beta::T
+end
+
 """
     CountingMax{K} <: AbstractProperty
     CountingMax(K=Single)
@@ -213,12 +225,15 @@ Positional Arguments
 * `problem` is the graph problem with tensor network information,
 * `property` is string specifying the task. Using the maximum independent set problem as an example, it can be one of
 
+    * [`PartitionFunction`](@ref)`()` for computing the partition function,
+
     * [`SizeMax`](@ref)`(k=Single)` for finding maximum-``k`` set sizes,
     * [`SizeMin`](@ref)`(k=Single)` for finding minimum-``k`` set sizes,
 
     * [`CountingMax`](@ref)`(k=Single)` for counting configurations with maximum-``k`` sizes,
     * [`CountingMin`](@ref)`(k=Single)` for counting configurations with minimum-``k`` sizes,
     * [`CountingAll`](@ref)`()` for counting all configurations,
+    * [`PartitionFunction`](@ref)`()` for counting all configurations,
     * [`GraphPolynomial`](@ref)`(; method=:finitefield, kwargs...)` for evaluating the graph polynomial,
 
     * [`SingleConfigMax`](@ref)`(k=Single; bounded=false)` for finding one maximum-``k`` configuration for each size,
@@ -247,6 +262,8 @@ function solve(gp::GraphProblem, property::AbstractProperty; T=Float64, usecuda=
         return asarray(post_invert_exponent.(res), res)
     elseif property isa CountingAll
         return contractx(gp, one(T); usecuda=usecuda)
+    elseif property isa PartitionFunction
+        return contractx(gp, exp(property.beta); usecuda=usecuda)
     elseif property isa CountingMax{Single}
         return contractx(gp, _x(CountingTropical{T,T}; invert=false); usecuda=usecuda)
     elseif property isa CountingMin{Single}
@@ -353,29 +370,29 @@ function has_noninteger_weights(problem::GraphProblem)
 end
 
 """
-    max_size(problem; usecuda=false)
+$TYPEDSIGNATURES
 
 Returns the maximum size of the graph problem. 
 A shorthand of `solve(problem, SizeMax(); usecuda=false)`.
 """
-max_size(m::GraphProblem; usecuda=false) = Int(sum(solve(m, SizeMax(); usecuda=usecuda)).n)  # floating point number is faster (BLAS)
+max_size(m::GraphProblem; usecuda=false)::Int = Int(sum(solve(m, SizeMax(); usecuda=usecuda)).n)  # floating point number is faster (BLAS)
 
 """
-    max_size_count(problem; usecuda=false)
+$TYPEDSIGNATURES
 
 Returns the maximum size and the counting of the graph problem.
 It is a shorthand of `solve(problem, CountingMax(); usecuda=false)`.
 """
-max_size_count(m::GraphProblem; usecuda=false) = (r = sum(solve(m, CountingMax(); usecuda=usecuda)); (Int(r.n), Int(r.c)))
+max_size_count(m::GraphProblem; usecuda=false)::Tuple{Int,Int} = (r = sum(solve(m, CountingMax(); usecuda=usecuda)); (Int(r.n), Int(r.c)))
 
 ########## memory estimation ###############
 """
-    estimate_memory(problem, property; T=Float64)
+$TYPEDSIGNATURES
 
 Memory estimation in number of bytes to compute certain `property` of a `problem`.
 `T` is the base type.
 """
-function estimate_memory(problem::GraphProblem, property::AbstractProperty; T=Float64)
+function estimate_memory(problem::GraphProblem, property::AbstractProperty; T=Float64)::Real
     _estimate_memory(tensor_element_type(T, length(labels(problem)), nflavor(problem), property), problem)
 end
 function estimate_memory(problem::GraphProblem, property::Union{SingleConfigMax{K,BOUNDED},SingleConfigMin{K,BOUNDED}}; T=Float64) where {K, BOUNDED}
@@ -419,6 +436,7 @@ function _estimate_memory(::Type{ET}, problem::GraphProblem) where ET
 end
 
 for (PROP, ET) in [
+        (:(PartitionFunction{T}), :(T)),
         (:(SizeMax{Single}), :(Tropical{T})), (:(SizeMin{Single}), :(Tropical{T})),
         (:(SizeMax{K}), :(ExtendedTropical{K,Tropical{T}})), (:(SizeMin{K}), :(ExtendedTropical{K,Tropical{T}})),
         (:(CountingAll), :T), (:(CountingMax{Single}), :(CountingTropical{T,T})), (:(CountingMin{Single}), :(CountingTropical{T,T})),

src/networks/Coloring.jl

@@ -5,7 +5,7 @@
             fixedvertices=Dict()
         )
 
-The [Vertex Coloring](https://psychic-meme-f4d866f8.pages.github.io/dev/generated/Coloring.html) problem.
+The [Vertex Coloring](https://queracomputing.github.io/GenericTensorNetworks.jl/dev/generated/Coloring/) problem.
 
 Positional arguments
 -------------------------------