Belief Propagation

examples/belief_propagation/bpexample.jl

@@ -0,0 +1,50 @@
+using ITensors
+using ITensorNetworks
+using ITensorNetworks:
+  delta_network,
+  formsubgraphs,
+  partition,
+  flatten_thicken_bonds,
+  flattened_inner_network,
+  construct_initial_mts,
+  update_all_mts,
+  get_single_site_expec,
+  iterate_single_site_expec,
+  contract_boundary_mps
+using KaHyPar
+using Compat
+
+#nxn GRID
+n = 4
+g = named_grid((n, n))
+s = siteinds("S=1/2", g)
+chi = 3
+
+#Random Tensor Network, Flatten it too
+psi = randomITensorNetwork(s; link_space=chi)
+psiflat, combiners = flatten_thicken_bonds(deepcopy(psi))
+
+v = (1, 1)
+
+#Apply Sz to site v and flatten that
+psiflatO = flatten_thicken_bonds(psi; ops=["Sz"], vops=[v], s=s, combiners=combiners)
+
+#Get the value of sz on v via exact contraction
+actual_sz = ITensors.contract(psiflatO)[1] / ITensors.contract(psiflat)[1]
+
+println("Actual value of Sz on site " * string(v) * " is " * string(actual_sz))
+
+nsites = 1
+println("First " * string(nsites) * " sites form a subgraph")
+dg_subgraphs = formsubgraphs(g, Int(n * n / nsites))
+mts = construct_initial_mts(psiflat, dg_subgraphs; init=(I...) -> allequal(I) ? 1 : 0)
+niters = 5
+
+iterate_single_site_expec(psiflat, psiflatO, mts, dg_subgraphs, niters, v)
+
+nsites = 4
+println("Now " * string(nsites) * " sites form a subgraph")
+dg_subgraphs = formsubgraphs(g, Int(n * n / nsites))
+mts = construct_initial_mts(psiflat, dg_subgraphs; init=(I...) -> allequal(I) ? 1 : 0)
+
+iterate_single_site_expec(psiflat, psiflatO, mts, dg_subgraphs, niters, v)

src/ITensorNetworks.jl

@@ -125,6 +125,11 @@ include("expect.jl")
 include("models.jl")
 include("tebd.jl")
 include("itensornetwork.jl")
+include("specialitensornetworks.jl")
+include("renameitensornetwork.jl")
+include("boundarymps.jl")
+include("subgraphs.jl")
+include("beliefpropagation.jl")
 include(joinpath("treetensornetwork", "treetensornetwork.jl"))
 
 include("exports.jl")

src/abstractitensornetwork.jl

@@ -450,6 +450,80 @@ end
 # TODO: rename `sqnorm` to match https://github.com/JuliaStats/Distances.jl?
 norm2(ψ::AbstractITensorNetwork; sequence) = contract_inner(ψ, ψ; sequence)
 
+#Run over all edges of an ITENSORNETWORK and combine any duplicated indices into one
+#Can take pre-existing combiners as an argument
+function combine_link_indices(psi::AbstractITensorNetwork; combiners=nothing)
+  if (combiners == nothing)
+    outcombiners = Dict{NamedDimEdge{Tuple},ITensor}()
+  end
+  es = edges(psi)
+  new_psi = deepcopy(psi)
+
+  for e in es
+    v1, v2 = src(e), dst(e)
+    if (combiners == nothing)
+      is = commoninds(psi, e)
+      C = combiner(is; tags=tags(is[1]))
+    else
+      C = combiners[e]
+    end
+    new_psi[v1] = new_psi[v1] * C
+    new_psi[v2] = new_psi[v2] * C
+    if (combiners == nothing)
+      outcombiners[e] = C
+    end
+  end
+
+  if (combiners == nothing)
+    return new_psi, outcombiners
+  else
+    return new_psi
+  end
+end
+
+#Run over all edges of an ITENSORNETWORK and uncombine any indices which were previously combined
+#Can take pre-existing combiners as an argument
+function uncombine_link_indices(
+  psi::AbstractITensorNetwork, combiners=Dict{NamedDimEdge{Tuple},ITensor}()
+)
+  es = edges(psi)
+  new_psi = deepcopy(psi)
+
+  for e in es
+    v1, v2 = src(e), dst(e)
+    if (haskey(combiners, e))
+      C = combiners[e]
+      new_psi[v1] = new_psi[v1] * C
+      new_psi[v2] = new_psi[v2] * C
+    end
+  end
+
+  return new_psi
+end
+
+#GIVEN AN ITENSOR NETWORK WITH DANGLING (PHYSICAL) BONDS, CONTRACT IT ONTO ITS CONJUGATE AND COMBINE ALL
+#DOUBLE INDICES. RETURN THE COMBINERS USED ON THE LINK INDICES
+#CAN APPLY LOCAL OPS ALONG THE PHYSICAL INDICES OF THE ITENSORNETWORK BEFORE CONTRACTION (ops is a list of strings, vops is a list of vertices, s is the indsnetwork)
+#CAN ALSO PASS COMBINERS FROM BEFORE TO ENSURE THE INDICES ARE CONSISTENT BETWEEN DIFFERENT APPLICATIONS OF THIS FUNCTION
+function flatten_thicken_bonds(
+  psi::AbstractITensorNetwork; ops=nothing, vops=nothing, s=nothing, combiners=nothing
+)
+  psiin = deepcopy(psi)
+  outcombiners = Dict{NamedDimEdge{Tuple},ITensor}()
+
+  if (ops != nothing)
+    for (op, v) in zip(ops, vops)
+      O = ITensor(Op(op, v), s)
+      psiin[v] = noprime!(O * psiin[v])
+    end
+  end
+
+  flatpsi = flattened_inner_network(psiin, psi)
+  flatpsi = rename_vertices_itn(flatpsi, Dictionary(vertices(flatpsi), vertices(psi)))
+
+  return combine_link_indices(flatpsi; combiners=combiners)
+end
+
 #
 # Printing
 #