Finding minimum sets (#15)

* finding minimum sets

* fix tests

Author: GiggleLiu

docs/src/ref.md

@@ -34,6 +34,8 @@ nflavor
 is_independent_set
 mis_compactify!
 
+is_maximal_independent_set
+
 cut_size
 cut_assign
 
@@ -48,10 +50,13 @@ is_good_vertex_coloring
 SizeMax
 CountingAll
 CountingMax
+CountingMin
 GraphPolynomial
 SingleConfigMax
+SingleConfigMin
 ConfigsAll
 ConfigsMax
+ConfigsMin
 ```
 
 ## Element Algebras

examples/IndependentSet.jl

@@ -62,15 +62,21 @@ timespacereadwrite_complexity(problem)
 # ### Maximum independent set size ``\alpha(G)``
 maximum_independent_set_size = solve(problem, SizeMax())[]
 
+# Function [`solve`](@ref) takes two positional arguments, the problem instance and the wanted property.
+# Here [`SizeMax`](@ref) means finding the solution with maximum set size.
+
 # ### Counting properties
 # ##### counting all independent sets
 count_all_independent_sets = solve(problem, CountingAll())[]
 
-# Function [`solve`](@ref) takes two positional arguments, the problem instance and the wanted property.
+# Here [`CountingAll`](@ref) means finding all possible solutions.
 
 # ##### counting independent sets with sizes ``\alpha(G)`` and ``\alpha(G)-1``
 count_max2_independent_sets = solve(problem, CountingMax(2))[]
 
+# Here [`CountingMax`](@ref) means finding solutions with largest-2 sizes,
+# the positional argument has default value `1`.
+
 # ##### independence polynomial
 # The graph polynomial for the independent set problem is known as the independence polynomial.
 # ```math
@@ -85,7 +91,7 @@ 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.
+# One can use [`SingleConfigMax`](@ref) to find one of the solution with largest set size, and it has two implementations.
 # 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) (see [Performance Tips](@ref)).
@@ -98,7 +104,8 @@ show_graph(graph; locs=locations, vertex_colors=
     [iszero(max_config.c.data[i]) ? "white" : "red" for i=1:nv(graph)])
 
 # ##### enumeration of all MISs
-# There are two approaches to enumerate all best-K solutions.
+# One can use [`ConfigsMax`](@ref) to find all solutions with largest-K set sizes.
+# It also has two implementations.
 # The bounded (default) version is always prefered because it can significantly use the memory usage.
 all_max_configs = solve(problem, ConfigsMax(1; bounded=true))[]
 
@@ -117,6 +124,7 @@ Compose.compose(context(),
      ntuple(k->(context((k-1)/m, 0.0, 1.2/m, 1.0), imgs[k]), m)...)
 
 # ##### enumeration of all IS configurations
+# One can use [`ConfigsAll`](@ref) to enumerate all sets satisfying the problem constraint.
 all_independent_sets = solve(problem, ConfigsAll())[]
 
 # To save/read a set of configuration to disk, one can type the following
@@ -130,6 +138,14 @@ loaded_sets = load_configs(filename; format=:binary, bitlength=10)
 #     When loading data, one needs to provide the `bitlength` if the data is saved in binary format.
 #     Because the bitstring length is not stored.
 
+# !!! note
+#     Check section [Maximal independent set problem](@ref) for examples of finding graph properties related to minimum sizes:
+#     * [`SizeMin`](@ref) for enumerating minimum set size,
+#     * [`CountingMin`](@ref) for enumerating minimum set size,
+#     * [`SingleConfigMin`](@ref) for enumerating minimum set size,
+#     * [`ConfigsMin`](@ref) for enumerating minimum set size,
+
+
 # ## Weights and open vertices
 # [`IndependentSet`](@ref) accepts weights as a key word argument.
 # The following code computes the weighted MIS problem.

examples/MaximalIS.jl

@@ -57,14 +57,24 @@ timespacereadwrite_complexity(problem)
 # ```
 # where ``b_k`` is the number of maximal independent sets of size ``k`` in graph ``G=(V, E)``.
 
-max_config = solve(problem, GraphPolynomial())[]
+maximal_indenpendence_polynomial = solve(problem, GraphPolynomial())[]
 
 # One can see the first several coefficients are 0, because it only counts the maximal independent sets, 
+# The minimum maximal independent set size is also known as the independent domination number.
+# It can be computed with the [`SizeMin`](@ref) property:
+independent_domination_number = solve(problem, SizeMin())[]
+
+# Similarly, we have its counting [`CountingMin`](@ref):
+counting_min_maximal_independent_set = solve(problem, CountingMin())[]
 
 # ### Configuration properties
 # ##### finding all maximal independent set
 maximal_configs = solve(problem, ConfigsAll())[]
 
+all(c->is_maximal_independent_set(g, i), maximal_configs)
+
+#
+
 imgs = ntuple(k->show_graph(graph;
                 locs=locations, scale=0.25,
                 vertex_colors=[iszero(maximal_configs[k][i]) ? "white" : "red"
@@ -76,4 +86,19 @@ Compose.set_default_graphic_size(18cm, 12cm); Compose.compose(context(),
 # This result should be consistent with that given by the [Bron Kerbosch algorithm](https://en.wikipedia.org/wiki/Bron%E2%80%93Kerbosch_algorithm) on the complement of Petersen graph.
 cliques = maximal_cliques(complement(graph))
 
-# For sparse graphs, the generic tensor network approach is usually much faster and memory efficient than the Bron Kerbosch algorithm.
+# For sparse graphs, the generic tensor network approach is usually much faster and memory efficient than the Bron Kerbosch algorithm.
+
+# ##### finding minimum maximal independent set
+# It is the [`ConfigsMin`](@ref) property in the program.
+minimum_maximal_configs = solve(problem, ConfigsMin())[]
+
+imgs2 = ntuple(k->show_graph(graph;
+                locs=locations, scale=0.25,
+                vertex_colors=[iszero(minimum_maximal_configs[k][i]) ? "white" : "red"
+                for i=1:nv(graph)]), length(minimum_maximal_configs));
+
+Compose.set_default_graphic_size(15cm, 12cm); Compose.compose(context(),
+     ntuple(k->(context((mod1(k,4)-1)/4, ((k-1)÷5)/3, 1.2/4, 1.0/3), imgs2[k]), 10)...)
+
+# Similarly, if one is only interested in computing one of the minimum sets,
+# one can use the graph property [`SingleConfigMin`](@ref).

src/GraphTensorNetworks.jl

@@ -24,7 +24,7 @@ export bestk_solutions
 export contractx, contractf, graph_polynomial, max_size, max_size_count
 
 # Graphs
-export random_regular_graph, diagonal_coupled_graph, is_independent_set
+export random_regular_graph, diagonal_coupled_graph, is_independent_set, is_maximal_independent_set
 export square_lattice_graph, unit_disk_graph, random_diagonal_coupled_graph, random_square_lattice_graph
 export line_graph
 
@@ -35,7 +35,7 @@ export mis_compactify!, cut_assign, cut_size, num_paint_shop_color_switch, paint
 export is_good_vertex_coloring
 
 # Interfaces
-export solve, SizeMax, CountingAll, CountingMax, GraphPolynomial, SingleConfigMax, ConfigsAll, ConfigsMax
+export solve, SizeMax, SizeMin, CountingAll, CountingMax, CountingMin, GraphPolynomial, SingleConfigMax, SingleConfigMin, ConfigsAll, ConfigsMax, ConfigsMin
 
 # Utilities
 export save_configs, load_configs

src/arithematics.jl

@@ -193,6 +193,8 @@ struct ConfigEnumerator{N,S,C}
 end
 
 Base.length(x::ConfigEnumerator{N}) where N = length(x.data)
+Base.iterate(x::ConfigEnumerator{N}) where N = iterate(x.data)
+Base.iterate(x::ConfigEnumerator{N}, state) where N = iterate(x.data, state)
 Base.getindex(x::ConfigEnumerator, i) = x.data[i]
 Base.:(==)(x::ConfigEnumerator{N,S,C}, y::ConfigEnumerator{N,S,C}) where {N,S,C} = Set(x.data) == Set(y.data)
 

src/bounding.jl

@@ -155,29 +155,36 @@ function solution_ad(code::Union{NestedEinsum,SlicedEinsum}, @nospecialize(xsa),
     # compute masks from cached tensors
     @debug "generating masked tree..."
     mt = generate_masktree(SingleConfig(), code, c, ymask, size_dict)
-    n, read_config!(code, mt, Dict())
+    config = read_config!(code, mt, Dict())
+    if length(config) !== length(labels(code))  # equal to the # of degree of freedoms
+        error("configuration `$(config)` is not fully determined!")
+    end
+    n, config
 end
 
 # get the solution configuration from gradients.
 function read_config!(code::SlicedEinsum, mt, out)
     read_config!(code.eins, mt, out)
 end
+
 function read_config!(code::NestedEinsum, mt, out)
     for (arg, ix, sibling) in zip(code.args, getixs(code.eins), mt.siblings)
         if OMEinsum.isleaf(arg)
-            assign = convert(Array, sibling.content)  # note: the content can be CuArray
-            if length(ix) == 1
-                if !assign[1] && assign[2]
-                    out[ix[1]] = 1
-                elseif !assign[2] && assign[1]
-                    out[ix[1]] = 0
-                else
-                    error("invalid assign $(assign)")
+            mask = convert(Array, sibling.content)  # note: the content can be CuArray
+            for ci in CartesianIndices(mask)
+                if mask[ci]
+                    for k in 1:ndims(mask)
+                        if haskey(out, ix[k])
+                            @assert out[ix[k]] == ci.I[k] - 1
+                        else
+                            out[ix[k]] = ci.I[k] - 1
+                        end
+                    end
                 end
             end
         else  # nested
             read_config!(arg, sibling, out)
         end
     end
     return out
-end
+end

src/configurations.jl

@@ -1,34 +1,37 @@
 """
-    best_solutions(problem; all=false, usecuda=false)
+    best_solutions(problem; all=false, usecuda=false, invert=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.
+* If `invert` is true, find the minimum.
 """
-function best_solutions(gp::GraphProblem; all=false, usecuda=false)
+function best_solutions(gp::GraphProblem; all=false, usecuda=false, invert=false)
     if all && usecuda
         throw(ArgumentError("ConfigEnumerator can not be computed on GPU!"))
     end
-    xst = generate_tensors(l->TropicalF64.(get_weights(gp, l)), gp)
+    xst = generate_tensors(l->(vs = TropicalF64.(get_weights(gp, l)); invert ? pre_invert_exponent.(vs) : vs), gp)
     ymask = trues(fill(2, length(getiyv(gp.code)))...)
     if usecuda
         xst = CuArray.(xst)
         ymask = CuArray(ymask)
     end
     if all
         # we use `Float64` types because we want to support weighted graphs.
-        xs = generate_tensors(fx_solutions(gp, CountingTropical{Float64,Float64}, all), gp)
-        return bounding_contract(AllConfigs{1}(), gp.code, xst, ymask, xs)
+        xs = generate_tensors(fx_solutions(gp, CountingTropical{Float64,Float64}, all, invert), gp)
+        ret = bounding_contract(AllConfigs{1}(), gp.code, xst, ymask, xs)
+        return invert ? post_invert_exponent.(ret) : ret
     else
         @assert ndims(ymask) == 0
         t, res = solution_ad(gp.code, xst, ymask)
-        return fill(CountingTropical(asscalar(t).n, ConfigSampler(StaticBitVector(map(l->res[l], 1:length(res))))))
+        ret = fill(CountingTropical(asscalar(t).n, ConfigSampler(StaticBitVector(map(l->res[l], 1:length(res))))))
+        return invert ? post_invert_exponent.(ret) : ret
     end
 end
 
 """
-    solutions(problem, basetype; all, usecuda=false)
+    solutions(problem, basetype; all, usecuda=false, invert=false)
     
 General routine to find solutions without bounding,
 
@@ -39,25 +42,27 @@ General routine to find solutions without 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 solutions(gp::GraphProblem, ::Type{BT}; all, usecuda=false) where BT
+function solutions(gp::GraphProblem, ::Type{BT}; all::Bool, usecuda::Bool=false, invert::Bool=false) where BT
     if all && usecuda
         throw(ArgumentError("ConfigEnumerator can not be computed on GPU!"))
     end
-    return contractf(fx_solutions(gp, BT, all), gp; usecuda=usecuda)
+    ret = contractf(fx_solutions(gp, BT, all, invert), gp; usecuda=usecuda)
+    return invert ? post_invert_exponent.(ret) : ret
 end
 
 """
-    best2_solutions(problem; all=true, usecuda=false)
+    best2_solutions(problem; all=true, usecuda=false, invert=false)
 
 Finding optimal and suboptimal solutions.
 """
-best2_solutions(gp::GraphProblem; all=true, usecuda=false) = solutions(gp, Max2Poly{Float64,Float64}; all=all, usecuda=usecuda)
+best2_solutions(gp::GraphProblem; all=true, usecuda=false, invert::Bool=false) = solutions(gp, Max2Poly{Float64,Float64}; all, usecuda, invert)
 
-function bestk_solutions(gp::GraphProblem, k::Int)
+function bestk_solutions(gp::GraphProblem, k::Int; invert::Bool=false)
     xst = generate_tensors(l->TropicalF64.(get_weights(gp, l)), gp)
     ymask = trues(fill(2, length(getiyv(gp.code)))...)
-    xs = generate_tensors(fx_solutions(gp, TruncatedPoly{k,Float64,Float64}, true), gp)
-    return bounding_contract(AllConfigs{k}(), gp.code, xst, ymask, xs)
+    xs = generate_tensors(fx_solutions(gp, TruncatedPoly{k,Float64,Float64}, true, invert), gp)
+    ret = bounding_contract(AllConfigs{k}(), gp.code, xst, ymask, xs)
+    return invert ? post_invert_exponent.(ret) : ret
 end
 
 """
@@ -69,12 +74,13 @@ e.g. when the problem is [`MaximalIS`](@ref), it computes all maximal independen
 all_solutions(gp::GraphProblem) = solutions(gp, Polynomial{Float64,:x}, all=true, usecuda=false)
 
 # return a mapping from label to onehot bitstrings (degree of freedoms).
-function fx_solutions(gp::GraphProblem, ::Type{BT}, all::Bool) where {BT}
+function fx_solutions(gp::GraphProblem, ::Type{BT}, all::Bool, invert::Bool) where {BT}
     syms = symbols(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)
-        _onehotv.(Ref(T), vertex_index[l], flavors(gp), get_weights(gp, l))
+        ret = _onehotv.(Ref(T), vertex_index[l], flavors(gp), get_weights(gp, l))
+        invert ? pre_invert_exponent.(ret) : ret
     end
 end
 function _onehotv(::Type{Polynomial{BS,X}}, x, v, w) where {BS,X}