diff --git a/docs/make.jl b/docs/make.jl
index 4bda2da6..1ce117fd 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -23,7 +23,7 @@ makedocs(;
     repo="https://github.com/Happy-Diode/GraphTensorNetworks.jl/blob/{commit}{path}#{line}",
     sitename="GraphTensorNetworks.jl",
     format=Documenter.HTML(;
-        prettyurls=false,
+        prettyurls=get(ENV, "CI", "false") == "true",
         canonical="https://Happy-Diode.github.io/GraphTensorNetworks.jl",
         assets=String[indigo],
     ),
diff --git a/docs/src/ref.md b/docs/src/ref.md
index d2eca105..b2bce4ab 100644
--- a/docs/src/ref.md
+++ b/docs/src/ref.md
@@ -35,6 +35,11 @@ mis_compactify!
 
 cut_size
 cut_assign
+
+num_paint_shop_color_switch
+paint_shop_coloring_from_config
+
+is_good_vertex_coloring
 ```
 
 ## Properties
@@ -90,10 +95,12 @@ MergeGreedy
 #### Graph
 ```@docs
 show_graph
+spring_layout
 
 diagonal_coupled_graph
 square_lattice_graph
 unit_disk_graph
+line_graph
 
 random_diagonal_coupled_graph
 random_square_lattice_graph
diff --git a/examples/Coloring.jl b/examples/Coloring.jl
index a9efce3e..87a91729 100644
--- a/examples/Coloring.jl
+++ b/examples/Coloring.jl
@@ -47,4 +47,27 @@ show_graph(graph; locs=locations)
 # We construct the tensor network for the 3-coloring problem as
 problem = Coloring{3}(graph);
 
-# ## Solving properties
\ No newline at end of file
+# ## Solving properties
+# ##### counting all possible coloring
+num_of_coloring = solve(problem, CountingAll())[]
+
+# ##### finding one best coloring
+single_solution = solve(problem, SingleConfigMax())[]
+
+is_good_vertex_coloring(graph, single_solution.c.data)
+
+vertex_color_map = Dict(0=>"red", 1=>"green", 2=>"blue")
+
+show_graph(graph; locs=locations, vertex_colors=[vertex_color_map[Int(c)] for c in single_solution.c.data])
+
+# Let us try to solve the same issue on its line graph, a graph that generated by mapping an edge to a vertex and two edges sharing a common vertex will be connected.
+linegraph = line_graph(graph)
+
+show_graph(linegraph; locs=[0.5 .* (locations[e.src] .+ locations[e.dst]) for e in edges(graph)])
+
+# Let us construct the tensor network and see if there are solutions.
+lineproblem = Coloring{3}(linegraph);
+
+num_of_coloring = solve(lineproblem, CountingAll())[]
+
+# You will see a zero printed, meaning no solution for the 3-coloring on edges of a Petersen graph.
\ No newline at end of file
diff --git a/examples/IndependentSet.jl b/examples/IndependentSet.jl
index 426b034a..97053233 100644
--- a/examples/IndependentSet.jl
+++ b/examples/IndependentSet.jl
@@ -89,7 +89,7 @@ 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.
 # Having value 1 at i-th bit means vertex ``i`` is in the maximum independent set.
 # One can visualize this MIS with the following function.
-show_graph(graph; locs=locations, colors=[iszero(max_config.c.data[i]) ? "white" : "red"
+show_graph(graph; locs=locations, vertex_colors=[iszero(max_config.c.data[i]) ? "white" : "red"
                                  for i=1:nv(graph)])
 
 # ##### enumeration of all MISs
@@ -101,12 +101,12 @@ using Compose
 
 m = length(all_max_configs.c)
 
-imgs = ntuple(k->(context((k-1)/m, 0.0, 1.2/m, 1.0), show_graph(graph;
+imgs = ntuple(k->show_graph(graph;
                             locs=locations, scale=0.25,
-                            colors=[iszero(all_max_configs.c[k][i]) ? "white" : "red"
-                                 for i=1:nv(graph)])), m)
+                            vertex_colors=[iszero(all_max_configs.c[k][i]) ? "white" : "red"
+                            for i=1:nv(graph)]), m);
 
-Compose.set_default_graphic_size(18cm, 4cm); Compose.compose(context(), imgs...)
+Compose.set_default_graphic_size(18cm, 4cm); Compose.compose(context(), ntuple(k->(context((k-1)/m, 0.0, 1.2/m, 1.0), imgs[k]), m)...)
 
 # ##### enumeration of all IS configurations
 all_independent_sets = solve(problem, ConfigsAll())[]
@@ -129,7 +129,7 @@ problem = IndependentSet(graph; weights=collect(1:10))
 
 max_config_weighted = solve(problem, SingleConfigMax())[]
 
-show_graph(graph; locs=locations, colors=
+show_graph(graph; locs=locations, vertex_colors=
           [iszero(max_config_weighted.c.data[i]) ? "white" : "red" for i=1:nv(graph)])
 
 # The following code computes the MIS tropical tensor (reference to be added) with open vertices 1 and 2.
diff --git a/examples/Matching.jl b/examples/Matching.jl
index c765ab32..1443bf2d 100644
--- a/examples/Matching.jl
+++ b/examples/Matching.jl
@@ -48,9 +48,25 @@ show_graph(graph; locs=locations)
 problem = Matching(graph);
 
 # ## Solving properties
+# ### Maximum matching
+# ### Configuration properties
+max_matching = solve(problem, SizeMax())[]
+# The largest number of matching is 5, which means we have a perfect matching (vertices are all paired).
+
+# ##### matching polynomial
 # The graph polynomial defined for the independence problem is known as the matching polynomial.
 # Here, we adopt the first definition in the [wiki page](https://en.wikipedia.org/wiki/Matching_polynomial).
 # ```math
 # M(G, x) = \sum\limits_{k=1}^{|V|/2} c_k x^k,
 # ```
-# where ``k`` is the number of matches, and coefficients ``c_k`` are the corresponding counting.
\ No newline at end of file
+# where ``k`` is the number of matches, and coefficients ``c_k`` are the corresponding counting.
+
+matching_poly = solve(problem, GraphPolynomial())[]
+
+# ## Configuration properties
+
+# ##### one of the perfect matches
+match_config = solve(problem, SingleConfigMax())[]
+
+# Let us show the result by coloring the matched edges to red
+show_graph(graph; locs=locations, edge_colors=[isone(match_config.c.data[i]) ? "red" : "black" for i=1:ne(graph)])
\ No newline at end of file
diff --git a/examples/MaxCut.jl b/examples/MaxCut.jl
index 449354c1..df707bbe 100644
--- a/examples/MaxCut.jl
+++ b/examples/MaxCut.jl
@@ -71,5 +71,5 @@ max_cut_size_verify = cut_size(graph, max_vertex_config)
 
 # You should see a consistent result as above `max_cut_size`.
 
-show_graph(graph; locs=locations, colors=[
+show_graph(graph; locs=locations, vertex_colors=[
         iszero(max_vertex_config[i]) ? "white" : "red" for i=1:nv(graph)])
\ No newline at end of file
diff --git a/examples/MaximalIS.jl b/examples/MaximalIS.jl
index aa7e572f..b68d1c94 100644
--- a/examples/MaximalIS.jl
+++ b/examples/MaximalIS.jl
@@ -9,7 +9,7 @@
 #     * how to compute weighted graphs and handle open vertices.
 
 # ## Problem definition
-using GraphTensorNetworks, Graphs
+using GraphTensorNetworks, Graphs, Compose
 
 # In graph theory, a [maximal independent set](https://en.wikipedia.org/wiki/Maximal_independent_set) is an independent set that is not a subset of any other independent set.
 # It is different from maximum independent set because it does not require the set to have the max size.
@@ -26,6 +26,7 @@ locations = [[rot15(0.0, 1.0, i) for i=0:4]..., [rot15(0.0, 0.6, i) for i=0:4]..
 show_graph(graph; locs=locations)
 
 # ## Tensor network representation
+# Type [`MaximalIS`](@ref) can be used for constructing the tensor network with optimized contraction order for solving a maximal independent set problem.
 # For a vertex ``v\in V``, we define a boolean degree of freedom ``s_v\in\{0, 1\}``.
 # We defined the restriction on its neighbourhood ``N[v]``:
 # ```math
@@ -58,9 +59,16 @@ max_config = solve(problem, GraphPolynomial())[]
 
 # ### Configuration properties
 # ##### finding all maximal independent set
-max_edge_config = solve(problem, ConfigsAll())[]
+maximal_configs = solve(problem, ConfigsAll())[]
+
+imgs = ntuple(k->show_graph(graph;
+                            locs=locations, scale=0.25,
+                            vertex_colors=[iszero(maximal_configs[k][i]) ? "white" : "red"
+                            for i=1:nv(graph)]), length(maximal_configs));
+
+Compose.set_default_graphic_size(18cm, 12cm); Compose.compose(context(), ntuple(k->(context((mod1(k,5)-1)/5, ((k-1)÷5)/3, 1.2/5, 1.0/3), imgs[k]), 15)...)
 
 # 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.
-maximal_cliques = maximal_cliques(complement(graph))
+cliques = maximal_cliques(complement(graph))
 
 # For sparse graphs, the generic tensor network approach is usually much faster and memory efficient.
\ No newline at end of file
diff --git a/examples/PaintShop.jl b/examples/PaintShop.jl
index 4a9c90ba..389a3ddb 100644
--- a/examples/PaintShop.jl
+++ b/examples/PaintShop.jl
@@ -9,10 +9,96 @@
 #     * how to compute weighted graphs and handle open vertices.
 
 # ## Problme Definition
-# The [binary paint shop problem](http://m-hikari.com/ams/ams-2012/ams-93-96-2012/popovAMS93-96-2012-2.pdf).
+# The [binary paint shop problem](http://m-hikari.com/ams/ams-2012/ams-93-96-2012/popovAMS93-96-2012-2.pdf) is defined as follows:
+# we are given a ``2m`` length sequence containing ``m`` cars, where each car appears twice.
+# Each car need to be colored red in one occurrence, and blue in the other.
+# We need to choose which occurrence for each car to color with which color — the goal is to minimize the number of times we need to change the current color.
 
-# In the following, we are going to defined a binary paint shop problem for the following string
+# In the following, we use a character to represent a car,
+# and defined a binary paint shop problem as a string that each character appear exactly twice.
 
 using GraphTensorNetworks, Graphs
 
-sequence = "abaccb"
\ No newline at end of file
+sequence = collect("iadgbeadfcchghebif")
+
+# We can visualize this paint shop problem as a graph
+rot(a, b, θ) = cos(θ)*a + sin(θ)*b, cos(θ)*b - sin(θ)*a
+
+locations = [rot(0.0, 1.0, -0.25π - 1.5*π*(i-0.5)/length(sequence)) for i=1:length(sequence)]
+
+graph = path_graph(length(sequence))
+for i=1:length(sequence) 
+    j = findlast(==(sequence[i]), sequence)
+    i != j && add_edge!(graph, i, j)
+end
+
+show_graph(graph; locs=locations, texts=string.(sequence), edge_colors=[sequence[e.src] == sequence[e.dst] ? "blue" : "black" for e in edges(graph)])
+
+# Vertices connected by blue edges must have different colors,
+# and the goal becomes a min-cut problem defined on black edges.
+
+# ## Tensor network representation
+# Type [`PaintShop`](@ref) can be used for constructing the tensor network with optimized contraction order for solving a binary paint shop problem.
+# To obtain its tensor network representation, we associating car ``c_i`` (the ``i``-th character in our example) with a degree of freedom ``s_{c_i} \in \{0, 1\}``,
+# where we use ``0`` to denote the first appearance of a car is colored red and ``1`` to denote the first appearance of a car is colored blue.
+# For each black edges ``(i, i+1)``, we define an edge tensor labeld by ``(s_{c_i}, s_{c_{i+1}})`` as follows:
+# If both cars on this edge are their first or second appearance
+# ```math
+# B^{\rm parallel} = \begin{matrix}
+# x & 1 \\
+# 1 & x \\
+# \end{matrix},
+#
+# otherwise,
+# B^{\rm anti-parallel} = B^{\rm 10} = \begin{matrix}
+# 1 & x \\
+# x & 1 \\
+# \end{matrix}.
+# ```
+# It can be understood as, when both cars are their first appearance,
+# they tend to have the same configuration so that the color is not changed.
+# Otherwise, they tend to have different configuration to keep the color unchanged.
+
+# Let us contruct the problem instance as bellow.
+problem = PaintShop(sequence);
+
+# ### Counting properties
+# ##### maximal independence polynomial
+# The graph polynomial defined for the maximal independent set problem is
+# ```math
+# I_{\rm max}(G, x) = \sum_{k=0}^{\alpha(G)} b_k x^k,
+# ```
+# where ``b_k`` is the number of maximal independent sets of size ``k`` in graph ``G=(V, E)``.
+
+max_config = solve(problem, GraphPolynomial())[]
+
+# Since it only counts the maximal independent sets, the first several coefficients are 0.
+
+# ### Counting properties
+# ##### graph polynomial
+# The graph polynomial of the binary paint shop problem in our convension is defined as
+# ```math
+# D(G, x) = \sum_{k=0}^{\delta(G)} d_k x^k 
+# ```
+# where ``2d_k`` is the number of possible coloring with number of color changes ``2m-1-k``.
+paint_polynomial = solve(problem, GraphPolynomial())[]
+
+# ### Configuration properties
+# ##### finding best solutions
+best_configs = solve(problem, ConfigsMax())[]
+
+painting1 = paint_shop_coloring_from_config(best_configs.c.data[1]; initial=false)
+
+show_graph(graph; locs=locations, texts=string.(sequence), edge_colors=[sequence[e.src] == sequence[e.dst] ? "blue" : "black" for e in edges(graph)],
+                vertex_colors=[isone(c) ? "red" : "black" for c in painting1], vertex_text_color="white")
+
+#
+
+painting2 = paint_shop_coloring_from_config(best_configs.c.data[2]; initial=false)
+
+show_graph(graph; locs=locations, texts=string.(sequence), edge_colors=[sequence[e.src] == sequence[e.dst] ? "blue" : "black" for e in edges(graph)],
+                vertex_colors=[isone(c) ? "red" : "black" for c in painting2], vertex_text_color="white")
+
+# Since we have different choices of initial color, the number of best solution is 4.
+# The following function will check the solution and return you the number of color switchs
+num_paint_shop_color_switch(sequence, painting2)
\ No newline at end of file
diff --git a/src/GraphTensorNetworks.jl b/src/GraphTensorNetworks.jl
index 3257eae3..6f0bde90 100644
--- a/src/GraphTensorNetworks.jl
+++ b/src/GraphTensorNetworks.jl
@@ -15,6 +15,7 @@ export GreedyMethod, TreeSA, SABipartite, KaHyParBipartite, MergeVectors, MergeG
 export StaticBitVector, StaticElementVector, @bv_str
 export is_commutative_semiring
 export Max2Poly, TruncatedPoly, Polynomial, Tropical, CountingTropical, StaticElementVector, Mod, ConfigEnumerator, onehotv, ConfigSampler
+export CountingTropicalF64, CountingTropicalF32, TropicalF64, TropicalF32
 
 # Lower level APIs
 export AllConfigs, SingleConfig
@@ -25,11 +26,13 @@ export contractx, contractf, graph_polynomial, max_size, max_size_count
 # Graphs
 export random_regular_graph, diagonal_coupled_graph, is_independent_set
 export square_lattice_graph, unit_disk_graph, random_diagonal_coupled_graph, random_square_lattice_graph
+export line_graph
 
 # Tensor Networks (Graph problems)
 export GraphProblem, IndependentSet, MaximalIS, Matching, Coloring, optimize_code, set_packing, MaxCut, PaintShop, paintshop_from_pairs, UnWeighted
 export flavors, symbols, nflavor, get_weights
-export mis_compactify!, cut_assign, cut_size
+export mis_compactify!, cut_assign, cut_size, num_paint_shop_color_switch, paint_shop_coloring_from_config
+export is_good_vertex_coloring
 
 # Interfaces
 export solve, SizeMax, CountingAll, CountingMax, GraphPolynomial, SingleConfigMax, ConfigsAll, ConfigsMax
@@ -38,7 +41,7 @@ export solve, SizeMax, CountingAll, CountingMax, GraphPolynomial, SingleConfigMa
 export save_configs, load_configs
 
 # Visualization
-export show_graph
+export show_graph, spring_layout
 
 project_relative_path(xs...) = normpath(joinpath(dirname(dirname(pathof(@__MODULE__))), xs...))
 
diff --git a/src/arithematics.jl b/src/arithematics.jl
index 505fd4c6..026093dd 100644
--- a/src/arithematics.jl
+++ b/src/arithematics.jl
@@ -220,15 +220,18 @@ end
 Base.:(==)(x::ConfigSampler{N,S,C}, y::ConfigSampler{N,S,C}) where {N,S,C} = x.data == y.data
 
 function Base.:+(x::ConfigSampler{N,S,C}, y::ConfigSampler{N,S,C}) where {N,S,C}  # biased sampling: return `x`, maybe using random sampler is better.
-    return x
+    return randn() > 0.5 ? x : y
 end
 
 function Base.:*(x::ConfigSampler{L,S,C}, y::ConfigSampler{L,S,C}) where {L,S,C}
     ConfigSampler(x.data | y.data)
 end
 
-Base.zero(::Type{ConfigSampler{N,S,C}}) where {N,S,C} = ConfigSampler{N,S,C}(statictrues(StaticElementVector{N,S,C}))
-Base.one(::Type{ConfigSampler{N,S,C}}) where {N,S,C} = ConfigSampler{N,S,C}(staticfalses(StaticElementVector{N,S,C}))
+@generated function Base.zero(::Type{ConfigSampler{N,S,C}}) where {N,S,C}
+    ex = Expr(:call, :(StaticElementVector{$N,$S,$C}), Expr(:tuple, fill(typemax(UInt64), C)...))
+    :(ConfigSampler{N,S,C}($ex))
+end
+Base.one(::Type{ConfigSampler{N,S,C}}) where {N,S,C} = ConfigSampler{N,S,C}(zero(StaticElementVector{N,S,C}))
 Base.zero(::ConfigSampler{N,S,C}) where {N,S,C} = zero(ConfigSampler{N,S,C})
 Base.one(::ConfigSampler{N,S,C}) where {N,S,C} = one(ConfigSampler{N,S,C})
 
diff --git a/src/configurations.jl b/src/configurations.jl
index 49ba6e8c..b7b8c669 100644
--- a/src/configurations.jl
+++ b/src/configurations.jl
@@ -75,7 +75,7 @@ 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 {K,BT}
+function fx_solutions(gp::GraphProblem, ::Type{BT}, all::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)])
diff --git a/src/graph_polynomials.jl b/src/graph_polynomials.jl
index 085858f7..7d302982 100644
--- a/src/graph_polynomials.jl
+++ b/src/graph_polynomials.jl
@@ -96,6 +96,6 @@ function improved_counting(ys::AbstractArray...)
 end
 improved_counting(ys::Mod...) = Mods.CRT(ys...)
 
-for GP in [:IndependentSet, :Matching, :MaxCut, :MaximalIS, :PaintShop]
+for GP in [:IndependentSet, :Matching, :MaxCut, :MaximalIS, :PaintShop, :Coloring]
     @eval contractx(gp::$GP, x::T; usecuda=false) where T = contractf(l->Ref(x) .^ get_weights(gp, l), gp; usecuda=usecuda)
 end
\ No newline at end of file
diff --git a/src/graphs.jl b/src/graphs.jl
index 1e0d08d0..a5f18879 100644
--- a/src/graphs.jl
+++ b/src/graphs.jl
@@ -61,3 +61,22 @@ function unit_disk_graph(locs::AbstractVector, unit::Real)
     end
     return g
 end
+
+"""
+    line_graph(g::SimpleGraph)
+
+Returns the line graph of `g`.
+The line graph is generated by mapping an edge to a vertex and two edges sharing a common vertex will be connected.
+"""
+function line_graph(graph::SimpleGraph)
+    edge_list = collect(edges(graph))
+    linegraph = SimpleGraph(length(edge_list))
+    for (i, ei) in enumerate(edge_list)
+        for (j, ej) in enumerate(edge_list)
+            if ei.src ∈ (ej.src, ej.dst) || ei.dst ∈ (ej.src, ej.dst)
+                add_edge!(linegraph, i, j)
+            end
+        end
+    end
+    return linegraph
+end
diff --git a/src/networks/Coloring.jl b/src/networks/Coloring.jl
index edbff2fb..dc28295e 100644
--- a/src/networks/Coloring.jl
+++ b/src/networks/Coloring.jl
@@ -2,7 +2,7 @@
     Coloring{K,CT<:AbstractEinsum} <: GraphProblem
     Coloring{K}(graph; openvertices=(), optimizer=GreedyMethod(), simplifier=nothing)
 
-[Vertex Coloring](https://psychic-meme-f4d866f8.pages.github.io/dev/tutorials/Coloring/) problem.
+[Vertex Coloring](https://psychic-meme-f4d866f8.pages.github.io/dev/tutorials/Coloring.html) problem.
 `optimizer` and `simplifier` are for tensor network optimization, check `optimize_code` for details.
 """
 struct Coloring{K,CT<:AbstractEinsum} <: GraphProblem
@@ -42,3 +42,15 @@ end
 # coloring vertex tensor
 coloringv(vals::Vector{T}) where T = vals
 
+# utilities
+"""
+    is_good_vertex_coloring(graph::SimpleGraph, config)
+
+Returns true if the coloring specified by config is a valid one, i.e. does not violate the contraints of vertices of an edges having different colors.
+"""
+function is_good_vertex_coloring(graph::SimpleGraph, config)
+    for e in edges(graph)
+        config[e.src] == config[e.dst] && return false
+    end
+    return true
+end
\ No newline at end of file
diff --git a/src/networks/IndependentSet.jl b/src/networks/IndependentSet.jl
index 6b07717d..0876a5dd 100644
--- a/src/networks/IndependentSet.jl
+++ b/src/networks/IndependentSet.jl
@@ -3,7 +3,7 @@
     IndependentSet(graph; weights=UnWeighted(), openvertices=(),
                  optimizer=GreedyMethod(), simplifier=nothing)
 
-The [independent set problem](https://psychic-meme-f4d866f8.pages.github.io/dev/tutorials/IndependentSet/) in graph theory.
+The [independent set problem](https://psychic-meme-f4d866f8.pages.github.io/dev/tutorials/IndependentSet.html) in graph theory.
 In the constructor, `weights` are the weights of vertices.
 `openvertices` specifies labels for the output tensor.
 `optimizer` and `simplifier` are for tensor network optimization, check `optimize_code` for details.
diff --git a/src/networks/Matching.jl b/src/networks/Matching.jl
index 8bc63537..14527240 100644
--- a/src/networks/Matching.jl
+++ b/src/networks/Matching.jl
@@ -2,7 +2,7 @@
     Matching{CT<:AbstractEinsum} <: GraphProblem
     Matching(graph; openvertices=(), optimizer=GreedyMethod(), simplifier=nothing)
 
-[Vertex matching](https://mathworld.wolfram.com/Matching.html) problem.
+[Vertex matching](https://psychic-meme-f4d866f8.pages.github.io/dev/tutorials/Matching.html) problem.
 `optimizer` and `simplifier` are for tensor network optimization, check `optimize_code` for details.
 """
 struct Matching{CT<:AbstractEinsum} <: GraphProblem
diff --git a/src/networks/MaxCut.jl b/src/networks/MaxCut.jl
index ccb8ff94..6cced226 100644
--- a/src/networks/MaxCut.jl
+++ b/src/networks/MaxCut.jl
@@ -3,7 +3,7 @@
     MaxCut(graph; weights=UnWeighted(), openvertices=(),
                 optimizer=GreedyMethod(), simplifier=nothing)
 
-[Cut](https://psychic-meme-f4d866f8.pages.github.io/dev/tutorials/MaxCut/) problem (or spin glass problem).
+[Cut](https://psychic-meme-f4d866f8.pages.github.io/dev/tutorials/MaxCut.html) problem (or spin glass problem).
 In the constructor, `weights` are the weights of edges.
 `optimizer` and `simplifier` are for tensor network optimization, check `optimize_code` for details.
 """
diff --git a/src/networks/MaximalIS.jl b/src/networks/MaximalIS.jl
index ea9da6bd..4e15a89b 100644
--- a/src/networks/MaximalIS.jl
+++ b/src/networks/MaximalIS.jl
@@ -3,7 +3,7 @@
     MaximalIS(graph; weights=UnWeighted(), openvertices=(),
                  optimizer=GreedyMethod(), simplifier=nothing)
 
-[Maximal independent set](https://psychic-meme-f4d866f8.pages.github.io/dev/tutorials/MaximalIS/) problem. In the constructor, `weights` are the weights of vertices.
+[Maximal independent set](https://psychic-meme-f4d866f8.pages.github.io/dev/tutorials/MaximalIS.html) problem. In the constructor, `weights` are the weights of vertices.
 `optimizer` and `simplifier` are for tensor network optimization, check `optimize_code` for details.
 """
 struct MaximalIS{CT<:AbstractEinsum,WT<:Union{UnWeighted, Vector}} <: GraphProblem
diff --git a/src/networks/PaintShop.jl b/src/networks/PaintShop.jl
index 1a3ced7a..9aa841bc 100644
--- a/src/networks/PaintShop.jl
+++ b/src/networks/PaintShop.jl
@@ -3,7 +3,7 @@
     PaintShop(labels::AbstractVector; openvertices=(),
              optimizer=GreedyMethod(), simplifier=nothing)
 
-The [binary paint shop problem](http://m-hikari.com/ams/ams-2012/ams-93-96-2012/popovAMS93-96-2012-2.pdf).
+The [binary paint shop problem](https://psychic-meme-f4d866f8.pages.github.io/dev/tutorials/PaintShop.html).
 
 Example
 -----------------------------------------
@@ -54,7 +54,7 @@ end
 
 flavors(::Type{<:PaintShop}) = [0, 1]
 get_weights(::PaintShop, sym) = [0, 1]
-symbols(gp::PaintShop) = getixsv(gp.code)
+symbols(gp::PaintShop) = getixsv(gp.code)  # !!! may not be unique
 
 function generate_tensors(fx, c::PaintShop)
     ixs = getixsv(c.code)
@@ -68,3 +68,45 @@ function paintshop_bond_tensor(a::T, b::T, if1::Bool, if2::Bool) where T
     return m
 end
 
+"""
+    num_paint_shop_color_switch(labels::AbstractVector, coloring::AbstractVector)
+
+Check the validity of the `coloring` and returns the number of color switches.
+"""
+function num_paint_shop_color_switch(labels::AbstractVector, coloring::AbstractVector)
+    # check validity of solution
+    @assert length(unique(coloring)) == 2 && length(labels) == length(coloring)
+    unique_labels = unique(labels)
+    for l in unique_labels
+        locs = findall(==(l), labels)
+        @assert length(locs) == 2
+        c1, c2 = coloring[locs]
+        @assert c1 != c2
+    end
+    # counting color switch
+    return count(i->coloring[i] != coloring[i+1], 1:length(coloring)-1)
+end
+
+"""
+    paint_shop_coloring_from_config(config; initial=false)
+
+Return a valid painting from the paint shop configuration (given by the configuration solvers).
+The `config` is a sequence of 0 and 1, where 0 means the color changed, 1 mean color unchanged.
+"""
+function paint_shop_coloring_from_config(config; initial::Bool=false)
+    res = falses(length(config)+1)
+    res[1] = initial
+    @inbounds for i=2:length(res)
+        res[i] = res[i-1] ⊻ (1-config[i-1])
+    end
+    return res
+end
+
+function fx_solutions(gp::PaintShop, ::Type{BT}, all::Bool) where {BT}
+    syms = symbols(gp)
+    T = (all ? set_type : sampler_type)(BT, length(syms), nflavor(gp))
+    counter = Ref(0)
+    return function (l)
+        _onehotv.(Ref(T), (counter[]+=1; counter[]), flavors(gp), get_weights(gp, l))
+    end
+end
\ No newline at end of file
diff --git a/src/networks/networks.jl b/src/networks/networks.jl
index 0b079ffa..2ad81e58 100644
--- a/src/networks/networks.jl
+++ b/src/networks/networks.jl
@@ -83,7 +83,7 @@ julia> getixsv(gp.code)
  [8, 10]
 
 julia> gp.code(GraphTensorNetworks.generate_tensors(f, gp)...)
-0-dimensional Array{TropicalNumbers.TropicalF64, 0}:
+0-dimensional Array{TropicalF64, 0}:
 4.0ₜ
 ```
 """
diff --git a/src/visualize.jl b/src/visualize.jl
index e42c13b6..983d9ce5 100644
--- a/src/visualize.jl
+++ b/src/visualize.jl
@@ -52,19 +52,22 @@ function get_rescaler(locs::AbstractVector{<:Tuple}, pad)
     return Rescaler(promote(xmin, xmax, ymin, ymax, pad)...)
 end
 
-default_node_style(scale, stroke_color, fill_color) = compose(context(), Viznet.nodestyle(:default, r=0.15cm*scale), Compose.stroke(stroke_color), fill(fill_color), linewidth(0.3mm*scale))
+default_vertex_style(scale, stroke_color, fill_color) = compose(context(), Viznet.nodestyle(:default, r=0.15cm*scale), Compose.stroke(stroke_color), fill(fill_color), linewidth(0.3mm*scale))
 default_text_style(scale, color) = Viznet.textstyle(:default, fontsize(4pt*scale), fill(color))
-default_bond_style(scale, color) = Viznet.bondstyle(:default, Compose.stroke(color), linewidth(0.3mm*scale))
+default_edge_style(scale, color) = Viznet.bondstyle(:default, Compose.stroke(color), linewidth(0.3mm*scale))
 
 """
     show_graph(locations, edges;
-        colors=["black", "black", ...],
+        vertex_colors=["black", "black", ...],
+        edge_colors=["black", "black", ...],
         texts=["1", "2", ...],
         format=SVG,
-        bond_color="black",
+        edge_color="black",
+        kwargs...
         )
+    show_graph(graph::SimpleGraph; locs=spring_layout(graph), kwargs...)
 
-Plots vertices at `locations` with colors specified by `colors` and texts specified by `texts`.
+Plots vertices at `locations` with vertex colors specified by `vertex_colors` and texts specified by `texts`.
 You will need a `VSCode`, `Pluto` notebook or `Jupyter` notebook to show the image.
 If you want to write this image to the disk without displaying it in a frontend, please try
 
@@ -74,10 +77,25 @@ julia> open("test.png", "w") do f
        end
 ```
 
-The `format` keyword argument can also be `Compose.SVG` or `Compose.PDF`.
+The `format` keyword argument can be `Compose.SVG` or `Compose.PDF`.
+
+Other keyword arguments
+-------------------------------
+* line, vertex and text
+    * scale::Float64 = 1.0
+    * pad::Float64 = 1.5
+* vertex
+    * vertex_text_color::String = "black"
+    * vertex_stroke_color = "black"
+    * vertex_fill_color = "white"
+* edge
+    * edge_color::String = "black"
+* image size in `cm`
+    * image_size::Float64 = 12
 """
 function show_graph(locations, edges;
-        colors=nothing,
+        vertex_colors=nothing,
+        edge_colors=nothing,
         texts = nothing,
         format=SVG, io=nothing,
         kwargs...)
@@ -85,7 +103,7 @@ function show_graph(locations, edges;
         dx, dy = 12cm, 12cm
         img = Compose.compose(context())
     else
-        img, (dx, dy) = viz_atoms(locations, edges; colors=colors, texts=texts, config=GraphDisplayConfig(; config_plotting(locations)..., kwargs...))
+        img, (dx, dy) = viz_atoms(locations, edges; vertex_colors=vertex_colors, edge_colors=edge_colors, texts=texts, config=GraphDisplayConfig(; config_plotting(locations)..., kwargs...))
     end
     if io === nothing
         Compose.set_default_graphic_size(dx, dy)
@@ -94,7 +112,7 @@ function show_graph(locations, edges;
         return format(io, dx, dy)(img)
     end
 end
-function show_graph(graph::SimpleGraph; locs, kwargs...)
+function show_graph(graph::SimpleGraph; locs=spring_layout(graph), kwargs...)
     show_graph(locs, [(e.src, e.dst) for e in edges(graph)]; kwargs...)
 end
 
@@ -110,50 +128,133 @@ function fit_image(rescaler::Rescaler, image_size, imgs...)
 end
 
 # Returns a 2-tuple of (image::Context, size)
-function viz_atoms(locs, edges; colors, texts, config)
+function viz_atoms(locs, edges; vertex_colors, edge_colors, texts, config)
     rescaler = get_rescaler(locs, config.pad)
-    img = _viz_atoms(rescaler.(locs), edges, colors, texts, config, getscale(rescaler))
+    img = _viz_atoms(rescaler.(locs), edges, vertex_colors, edge_colors, texts, config, getscale(rescaler))
     return fit_image(rescaler, config.image_size, img)
 end
 
 Base.@kwdef struct GraphDisplayConfig
-    # line, node and text
+    # line, vertex and text
     scale::Float64 = 1.0
     pad::Float64 = 1.5
 
-    # node
-    node_text_color::String = "black"
-    node_stroke_color = "black"
-    node_fill_color = "white"
-    # bond
-    bond_color::String = "black"
+    # vertex
+    vertex_text_color::String = "black"
+    vertex_stroke_color = "black"
+    vertex_fill_color = "white"
+    # edge
+    edge_color::String = "black"
     # image size in cm
     image_size::Float64 = 12
 end
 
-function _viz_atoms(locs, edges, colors, texts, config, rescale)
+function _viz_atoms(locs, edges, vertex_colors, edge_colors, texts, config, rescale)
     rescale = rescale * config.image_size * config.scale * 1.6
-    if colors !== nothing
-        @assert length(locs) == length(colors)
-        node_styles = [default_node_style(rescale, config.node_stroke_color, color) for color in colors]
+    if vertex_colors !== nothing
+        @assert length(locs) == length(vertex_colors)
+        vertex_styles = [default_vertex_style(rescale, config.vertex_stroke_color, color) for color in vertex_colors]
+    else
+        vertex_styles = fill(default_vertex_style(rescale, config.vertex_stroke_color, config.vertex_fill_color), length(locs))
+    end
+    if edge_colors !== nothing
+        @assert length(edges) == length(edge_colors)
+        edge_styles = [default_edge_style(rescale, color) for color in edge_colors]
     else
-        node_styles = fill(default_node_style(rescale, config.node_stroke_color, config.node_fill_color), length(locs))
+        edge_styles = fill(default_edge_style(rescale, config.edge_color), length(edges))
     end
     if texts !== nothing
         @assert length(locs) == length(texts)
     end
-    edge_style = default_bond_style(rescale, config.bond_color)
-    text_style = default_text_style(rescale, config.node_text_color)
+    text_style = default_text_style(rescale, config.vertex_text_color)
     img1 = Viznet.canvas() do
-        for (i, node) in enumerate(locs)
-            node_styles[i] >> node
-            if config.node_text_color !== "transparent"
-                text_style >> (node, texts === nothing ? "$i" : texts[i])
+        for (i, vertex) in enumerate(locs)
+            vertex_styles[i] >> vertex
+            if config.vertex_text_color !== "transparent"
+                text_style >> (vertex, texts === nothing ? "$i" : texts[i])
             end
         end
-        for (i, j) in edges
-            edge_style >> (locs[i], locs[j])
+        for (k, (i, j)) in enumerate(edges)
+            edge_styles[k] >> (locs[i], locs[j])
         end
     end
     Compose.compose(context(), img1)
+end
+
+"""
+Spring layout for graph plotting, returns a vector of vertex locations.
+
+!!! note
+    This function is copied from [`GraphPlot.jl`](https://github.com/JuliaGraphs/GraphPlot.jl),
+    where you can find more information about his function.
+"""
+function spring_layout(g::AbstractGraph,
+                       locs_x=2*rand(nv(g)).-1.0,
+                       locs_y=2*rand(nv(g)).-1.0;
+                       C=2.0,
+                       MAXITER=100,
+                       INITTEMP=2.0)
+
+    nvg = nv(g)
+    adj_matrix = adjacency_matrix(g)
+
+    # The optimal distance bewteen vertices
+    k = C * sqrt(4.0 / nvg)
+    k² = k * k
+
+    # Store forces and apply at end of iteration all at once
+    force_x = zeros(nvg)
+    force_y = zeros(nvg)
+
+    # Iterate MAXITER times
+    @inbounds for iter = 1:MAXITER
+        # Calculate forces
+        for i = 1:nvg
+            force_vec_x = 0.0
+            force_vec_y = 0.0
+            for j = 1:nvg
+                i == j && continue
+                d_x = locs_x[j] - locs_x[i]
+                d_y = locs_y[j] - locs_y[i]
+                dist²  = (d_x * d_x) + (d_y * d_y)
+                dist = sqrt(dist²)
+
+                if !( iszero(adj_matrix[i,j]) && iszero(adj_matrix[j,i]) )
+                    # Attractive + repulsive force
+                    # F_d = dist² / k - k² / dist # original FR algorithm
+                    F_d = dist / k - k² / dist²
+                else
+                    # Just repulsive
+                    # F_d = -k² / dist  # original FR algorithm
+                    F_d = -k² / dist²
+                end
+                force_vec_x += F_d*d_x
+                force_vec_y += F_d*d_y
+            end
+            force_x[i] = force_vec_x
+            force_y[i] = force_vec_y
+        end
+        # Cool down
+        temp = INITTEMP / iter
+        # Now apply them, but limit to temperature
+        for i = 1:nvg
+            fx = force_x[i]
+            fy = force_y[i]
+            force_mag  = sqrt((fx * fx) + (fy * fy))
+            scale      = min(force_mag, temp) / force_mag
+            locs_x[i] += force_x[i] * scale
+            locs_y[i] += force_y[i] * scale
+        end
+    end
+
+    # Scale to unit square
+    min_x, max_x = minimum(locs_x), maximum(locs_x)
+    min_y, max_y = minimum(locs_y), maximum(locs_y)
+    function scaler(z, a, b)
+        2.0*((z - a)/(b - a)) - 1.0
+    end
+    map!(z -> scaler(z, min_x, max_x), locs_x, locs_x)
+    map!(z -> scaler(z, min_y, max_y), locs_y, locs_y)
+
+    return collect(zip(locs_x, locs_y))
 end
\ No newline at end of file
diff --git a/test/arithematics.jl b/test/arithematics.jl
index 15aab854..b0f60b79 100644
--- a/test/arithematics.jl
+++ b/test/arithematics.jl
@@ -25,6 +25,7 @@ using GraphTensorNetworks: StaticBitVector
 end
 
 @testset "arithematics" begin
+    Random.seed!(2)
     for (a, b, c) in [
                     (TropicalF64(2), TropicalF64(8), TropicalF64(9)),
                     (CountingTropicalF64(2, 8), CountingTropicalF64(8, 9), CountingTropicalF64(9, 2)),
diff --git a/test/configurations.jl b/test/configurations.jl
index 87dfc274..911e08ea 100644
--- a/test/configurations.jl
+++ b/test/configurations.jl
@@ -59,36 +59,4 @@ end
             @test all(x->count_ones(x)==(i-1), s.data)
         end
     end
-end
-
-@testset "set packing" begin
-    sets = [[1, 2, 5], [1, 3], [2, 4], [3, 6], [2, 3, 6]]  # each set is a vertex
-    gp = set_packing(sets; optimizer=GreedyMethod())
-    res = best_solutions(gp; all=true)[]
-    @test res.n == 2
-    @test BitVector(Bool[0,0,1,1,0]) ∈ res.c.data
-    @test BitVector(Bool[1,0,0,1,0]) ∈ res.c.data
-    @test BitVector(Bool[0,1,1,0,0]) ∈ res.c.data
-end
-
-@testset "enumerating - coloring" begin
-    g = SimpleGraph(5)
-    for (i,j) in [(1,2),(2,3),(3,4),(4,1),(1,5),(2,4)]
-        add_edge!(g, i, j)
-    end
-    code = Coloring{3}(g; optimizer=GreedyMethod())
-    res = solutions(code, CountingTropicalF64; all=true)[]
-    @test length(res.c.data) == 12
-    g = smallgraph(:petersen)
-    code = Coloring{3}(g; optimizer=GreedyMethod())
-    res = solutions(code, CountingTropicalF64; all=true)[]
-    @test length(res.c.data) == 120
-end
-
-@testset "enumerating - matching" begin
-    g = smallgraph(:petersen)
-    code = Matching(g; optimizer=GreedyMethod())
-    res = solutions(code, CountingTropicalF64; all=true)[]
-    @test res.n == 5
-    @test length(res.c.data) == 6
 end
\ No newline at end of file
diff --git a/test/graph_polynomials.jl b/test/graph_polynomials.jl
index ebaefb9a..8f87d9d2 100644
--- a/test/graph_polynomials.jl
+++ b/test/graph_polynomials.jl
@@ -35,21 +35,4 @@ end
     p6 = graph_polynomial(gp, Val(:polynomial))[]
     p7 = graph_polynomial(gp, Val(:finitefield))[]
     @test p6 ≈ p7
-end
-
-@testset "match polynomial" begin
-    g = SimpleGraph(7)
-    for (i,j) in [(1,2),(2,3),(3,4),(4,5),(5,6),(6,1),(1,7)]
-        add_edge!(g, i, j)
-    end
-    @test graph_polynomial(Matching(g), Val(:polynomial))[] == Polynomial([1,7,13,5])
-    g = smallgraph(:petersen)
-    @test graph_polynomial(Matching(g), Val(:polynomial))[].coeffs == [6, 90, 145, 75, 15, 1][end:-1:1]
-end
-
-@testset "paint shop" begin
-    labels = collect("abaccb")
-    pb = PaintShop(labels)
-    @test solve(pb, SizeMax())[] == Tropical(3.0)
-    @test StaticBitVector(Bool[0,1,1,0,1]) ∈ solve(pb, ConfigsMax())[].c.data
 end
\ No newline at end of file
diff --git a/test/networks/Coloring.jl b/test/networks/Coloring.jl
new file mode 100644
index 00000000..008b9da8
--- /dev/null
+++ b/test/networks/Coloring.jl
@@ -0,0 +1,20 @@
+using Test, GraphTensorNetworks, Graphs
+
+@testset "enumerating - coloring" begin
+    g = SimpleGraph(5)
+    for (i,j) in [(1,2),(2,3),(3,4),(4,1),(1,5),(2,4)]
+        add_edge!(g, i, j)
+    end
+    code = Coloring{3}(g; optimizer=GreedyMethod())
+    res = solutions(code, CountingTropical{Float64,Float64}; all=true)[]
+    @test length(res.c.data) == 12
+    g = smallgraph(:petersen)
+    code = Coloring{3}(g; optimizer=GreedyMethod())
+    res = solutions(code, CountingTropicalF64; all=true)[]
+    @test length(res.c.data) == 120
+
+    c = solve(code, SingleConfigMax())[]
+    @test c.c.data ∈ res.c.data
+    @test is_good_vertex_coloring(g, c.c.data)
+end
+
diff --git a/test/networks/IndependentSet.jl b/test/networks/IndependentSet.jl
index 9567c84f..658ded80 100644
--- a/test/networks/IndependentSet.jl
+++ b/test/networks/IndependentSet.jl
@@ -11,4 +11,14 @@ using GraphTensorNetworks, Test, Graphs
     @test count(!iszero, m) == 12
     mis_compactify!(m)
     @test count(!iszero, m) == 3
-end
\ No newline at end of file
+end
+
+@testset "set packing" begin
+    sets = [[1, 2, 5], [1, 3], [2, 4], [3, 6], [2, 3, 6]]  # each set is a vertex
+    gp = set_packing(sets; optimizer=GreedyMethod())
+    res = best_solutions(gp; all=true)[]
+    @test res.n == 2
+    @test BitVector(Bool[0,0,1,1,0]) ∈ res.c.data
+    @test BitVector(Bool[1,0,0,1,0]) ∈ res.c.data
+    @test BitVector(Bool[0,1,1,0,0]) ∈ res.c.data
+end
diff --git a/test/networks/Matching.jl b/test/networks/Matching.jl
new file mode 100644
index 00000000..7f2005e5
--- /dev/null
+++ b/test/networks/Matching.jl
@@ -0,0 +1,19 @@
+using Test, GraphTensorNetworks, Graphs
+
+@testset "enumerating - matching" begin
+    g = smallgraph(:petersen)
+    code = Matching(g; optimizer=GreedyMethod())
+    res = solutions(code, CountingTropicalF64; all=true)[]
+    @test res.n == 5
+    @test length(res.c.data) == 6
+end
+
+@testset "match polynomial" begin
+    g = SimpleGraph(7)
+    for (i,j) in [(1,2),(2,3),(3,4),(4,5),(5,6),(6,1),(1,7)]
+        add_edge!(g, i, j)
+    end
+    @test graph_polynomial(Matching(g), Val(:polynomial))[] == Polynomial([1,7,13,5])
+    g = smallgraph(:petersen)
+    @test graph_polynomial(Matching(g), Val(:polynomial))[].coeffs == [6, 90, 145, 75, 15, 1][end:-1:1]
+end
\ No newline at end of file
diff --git a/test/networks/PaintShop.jl b/test/networks/PaintShop.jl
new file mode 100644
index 00000000..d7226625
--- /dev/null
+++ b/test/networks/PaintShop.jl
@@ -0,0 +1,11 @@
+using GraphTensorNetworks, Test
+
+@testset "paint shop" begin
+    labels = collect("abaccb")
+    pb = PaintShop(labels)
+    @test solve(pb, SizeMax())[] == Tropical(3.0)
+    c = solve(pb, SingleConfigMax())[].c.data
+    coloring = paint_shop_coloring_from_config(c)
+    @test num_paint_shop_color_switch(labels, coloring) == 2
+    @test StaticBitVector(Bool[0,1,1,0,1]) ∈ solve(pb, ConfigsMax())[].c.data
+end
\ No newline at end of file
diff --git a/test/networks/networks.jl b/test/networks/networks.jl
index 8a93e38f..6edfae45 100644
--- a/test/networks/networks.jl
+++ b/test/networks/networks.jl
@@ -1,3 +1,6 @@
 include("IndependentSet.jl")
 include("MaximalIS.jl")
-include("MaxCut.jl")
\ No newline at end of file
+include("MaxCut.jl")
+include("PaintShop.jl")
+include("Coloring.jl")
+include("Matching.jl")
\ No newline at end of file