Revert

src/beliefpropagationdmrg/graphsextensions.jl

@@ -0,0 +1,132 @@
+using Graphs: AbstractGraph
+using NamedGraphs.GraphsExtensions: a_star, add_edge!, dfs_tree, post_order_dfs_edges
+using NamedGraphs.GraphsExtensions:
+  default_root_vertex, degree, degrees, neighbors, edgetype, rem_edge!, vertextype
+
+#Given a graph, traverse it from start vertex to end vertex, covering each edge exactly once.
+#Complexity is O(length(edges(g)))
+function eulerian_path(g::AbstractGraph, start_vertex, end_vertex)
+  #Conditions on g for the required path to exist
+  if start_vertex != end_vertex
+    @assert isodd(degree(g, start_vertex) % 2)
+    @assert isodd(degree(g, end_vertex) % 2)
+    @assert all(
+      x -> iseven(x), degrees(g, setdiff(vertices(g), [start_vertex, end_vertex]))
+    )
+  else
+    @assert all(x -> iseven(x), degrees(g))
+  end
+
+  path = vertextype(g)[]
+  stack = vertextype(g)[]
+  current_vertex = end_vertex
+  g_modified = copy(g)
+  while !isempty(stack) || !iszero(degree(g_modified, current_vertex))
+    if iszero(degree(g_modified, current_vertex))
+      push!(path, current_vertex)
+      last_vertex = pop!(stack)
+      current_vertex = last_vertex
+    else
+      push!(stack, current_vertex)
+      vn = first(neighbors(g_modified, current_vertex))
+      rem_edge!(g_modified, edgetype(g_modified)(current_vertex, vn))
+      current_vertex = vn
+    end
+  end
+
+  push!(path, current_vertex)
+
+  return edgetype(g_modified)[
+    edgetype(g_modified)(path[i], path[i + 1]) for i in 1:(length(path) - 1)
+  ]
+end
+
+function eulerian_cycle(g::AbstractGraph, start_vertex)
+  return eulerian_path(g, start_vertex, start_vertex)
+end
+
+function make_all_degrees_even(g::AbstractGraph)
+  g_modified = copy(g)
+  vertices_odd_degree = collect(filter(v -> isodd(degree(g, v)), vertices(g_modified)))
+  while !isempty(vertices_odd_degree)
+    vertex_pairs = [
+      (vertices_odd_degree[i], vertices_odd_degree[j]) for
+      i in 1:length(vertices_odd_degree) for j in (i + 1):length(vertices_odd_degree)
+    ]
+    vertex_pair = first(sort(vertex_pairs; by=vp -> length(a_star(g, vp...))))
+    add_edge!(g_modified, edgetype(g_modified)(vertex_pair...))
+    vertices_odd_degree = filter(
+      v -> v != first(vertex_pair) && v != last(vertex_pair), vertices_odd_degree
+    )
+  end
+  return g_modified
+end
+
+#Given a graph, traverse it in a cycle from start_vertex and try to minimise the number of edges traversed more than once
+function _eulerian_cycle(g::AbstractGraph, start_vertex; add_additional_traversal=false)
+  g_modified = make_all_degrees_even(g::AbstractGraph)
+  path = eulerian_cycle(g_modified, start_vertex)
+
+  !add_additional_traversal && return path
+
+  modified_path = edgetype(g_modified)[]
+
+  for e in path
+    if src(e) ∉ neighbors(g, dst(e))
+      inner_path = a_star(g, src(e), dst(e))
+      append!(modified_path, inner_path)
+    else
+      push!(modified_path, e)
+    end
+  end
+
+  return modified_path
+end
+
+function _bp_region_plan(
+  g::AbstractGraph,
+  start_vertex=default_root_vertex(g);
+  nsites::Int=1,
+  add_additional_traversal=false,
+)
+  path = _eulerian_cycle(g, start_vertex; add_additional_traversal)
+  if nsites == 1
+    regions = [[v] for v in vcat(src.(path))]
+    @assert all(v -> only(v) ∈ vertices(g), regions)
+    return vcat(regions, reverse(regions))
+  else
+    regions = [e for e in path]
+    regions = filter(e -> e ∈ edges(g) || reverse(e) ∈ edges(g), regions)
+    regions = [[src(e), dst(e)] for e in regions]
+    return vcat(regions, reverse(regions))
+  end
+end
+
+function path_to_path(path)
+  verts = []
+  for e in path
+    if isempty(verts) || src(e) ≠ last(verts)
+      push!(verts, src(e))
+      push!(verts, dst(e))
+    end
+  end
+  return [[v] for v in verts]
+end
+
+function bp_region_plan(
+  g::AbstractGraph,
+  start_vertex=default_root_vertex(g);
+  nsites::Int=1,
+  add_additional_traversal=false,
+)
+  path = post_order_dfs_edges(g, start_vertex)
+  if nsites == 1
+    regions = path_to_path(path)
+    @assert all(v -> only(v) ∈ vertices(g), regions)
+    return vcat(regions, reverse(regions))
+  else
+    regions = filter(e -> e ∈ edges(g) || reverse(e) ∈ edges(g), path)
+    @assert all(e -> e ∈ regions || reverse(e) ∈ regions, edges(g))
+    return regions
+  end
+end

src/beliefpropagationdmrg/main.jl

@@ -1,6 +1,5 @@
 using ITensorNetworks: random_tensornetwork
 using NamedGraphs.NamedGraphGenerators: named_comb_tree
-using NPZ
 using ITensors: expect
 using Random
 
@@ -12,16 +11,14 @@ Random.seed!(5634)
 
 function main()
   g = named_grid((24, 1); periodic=true)
-  g = renamer(g)
-  save = true
   L = length(vertices(g))
-  h, hl, J = 0.6, 0.2, 1.0
+  h, hlongitudinal, J = 0.6, 0.2, 1.0
   s = siteinds("S=1/2", g)
   χ = 2
   ψ0 = random_tensornetwork(s; link_space=χ)
 
-  H = ising(s; h, hl, J1=J)
-  tnos = opsum_to_tno(s, H)
+  H = ising(s; h, hl=hlongitudinal, J1=J)
+  tnos = opsum_to_tnos(s, H)
 
   energy_calc_fun =
     (tn, bp_cache) -> sum(expect(tn, H; alg="bp", (cache!)=Ref(bp_cache))) / L

src/beliefpropagationdmrg/tensornetworkoperators.jl

@@ -19,6 +19,7 @@ using Dictionaries
 
 using LinearAlgebra: eigvals
 
+#Convert the forests in forest_cover to trees by adding "dummy" edges and keep track of them
 function connect_forests(s::IndsNetwork)
   g = underlying_graph(s)
   fs = forest_cover(g)
@@ -39,9 +40,8 @@ function connect_forests(s::IndsNetwork)
   return fs_connected
 end
 
-function opsum_to_tno(
-  s::IndsNetwork, H::OpSum; cutoff::Float64=1e-14, insert_dummy_inds=false
-)
+#Convert an Opsum to a vector of tensor network operators whose sum matches opsum
+function opsum_to_tnos(s::IndsNetwork, H::OpSum; cutoff::Float64=1e-14)
   s_fs_connected = connect_forests(s)
   no_forests = length(s_fs_connected)
   tnos = ITensorNetwork[]
@@ -60,9 +60,6 @@ function opsum_to_tno(
     end
     tno = truncate(ttn(new_opsum, s_f); cutoff)
     tno = ITensorNetwork(tno)
-    if insert_dummy_inds
-      tno = insert_linkinds(tno, setdiff(edges(s), edges(tno)))
-    end
     push!(tnos, tno)
   end
 

src/beliefpropagationdmrg/utils.jl

@@ -0,0 +1,202 @@
+using ITensors: siteinds, Op, prime, OpSum, apply
+using Graphs: AbstractGraph, SimpleGraph, edges, vertices, is_tree, connected_components
+using NamedGraphs: NamedGraph, NamedEdge, NamedGraphs, rename_vertices
+using NamedGraphs.NamedGraphGenerators: named_grid, named_hexagonal_lattice_graph
+using NamedGraphs.GraphsExtensions:
+  decorate_graph_edges,
+  forest_cover,
+  add_edges,
+  rem_edges,
+  add_vertices,
+  rem_vertices,
+  disjoint_union,
+  subgraph,
+  src,
+  dst
+using NamedGraphs.PartitionedGraphs: PartitionVertex, partitionedge, unpartitioned_graph
+using ITensorNetworks:
+  BeliefPropagationCache,
+  AbstractITensorNetwork,
+  AbstractFormNetwork,
+  IndsNetwork,
+  ITensorNetwork,
+  insert_linkinds,
+  ttn,
+  union_all_inds,
+  neighbor_vertices,
+  environment,
+  messages,
+  update_factor,
+  message,
+  partitioned_tensornetwork,
+  bra_vertex,
+  ket_vertex,
+  operator_vertex,
+  default_cache_update_kwargs,
+  dual_index_map,
+  region_scalar,
+  renormalize_messages,
+  scalar_factors_quotient,
+  norm_sqr_network
+using DataGraphs: underlying_graph
+using ITensorNetworks.ModelHamiltonians: heisenberg
+using ITensors:
+  ITensor,
+  noprime,
+  dag,
+  noncommonind,
+  commonind,
+  replaceind,
+  dim,
+  noncommoninds,
+  delta,
+  replaceinds,
+  Trotter,
+  apply
+using ITensors.NDTensors: denseblocks
+using SplitApplyCombine: group
+
+using ITensors: siteinds, Op, prime, OpSum, apply
+using Graphs: AbstractGraph, SimpleGraph, edges, vertices, is_tree, connected_components
+using NamedGraphs: NamedGraph, NamedEdge, NamedGraphs, rename_vertices
+using NamedGraphs.NamedGraphGenerators: named_grid, named_hexagonal_lattice_graph
+using NamedGraphs.GraphsExtensions:
+  decorate_graph_edges,
+  forest_cover,
+  add_edges,
+  rem_edges,
+  add_vertices,
+  rem_vertices,
+  disjoint_union,
+  subgraph,
+  src,
+  dst
+using NamedGraphs.PartitionedGraphs: PartitionVertex, partitionedge, unpartitioned_graph
+using ITensorNetworks:
+  BeliefPropagationCache,
+  AbstractITensorNetwork,
+  AbstractFormNetwork,
+  IndsNetwork,
+  ITensorNetwork,
+  insert_linkinds,
+  ttn,
+  union_all_inds,
+  neighbor_vertices,
+  environment,
+  messages,
+  update_factor,
+  message,
+  partitioned_tensornetwork,
+  bra_vertex,
+  ket_vertex,
+  operator_vertex,
+  default_cache_update_kwargs,
+  dual_index_map
+using DataGraphs: underlying_graph
+using ITensorNetworks.ModelHamiltonians: heisenberg
+using ITensors:
+  ITensor,
+  noprime,
+  dag,
+  noncommonind,
+  commonind,
+  replaceind,
+  dim,
+  noncommoninds,
+  delta,
+  replaceinds
+using ITensors.NDTensors: denseblocks
+using Dictionaries: set!
+
+function BP_apply(
+  o::ITensor, ψ::AbstractITensorNetwork, bpc::BeliefPropagationCache; apply_kwargs...
+)
+  bpc = copy(bpc)
+  ψ = copy(ψ)
+  vs = neighbor_vertices(ψ, o)
+  envs = environment(bpc, PartitionVertex.(vs))
+  singular_values! = Ref(ITensor())
+  ψ = noprime(apply(o, ψ; envs, singular_values!, normalize=true, apply_kwargs...))
+  ψdag = prime(dag(ψ); sites=[])
+  if length(vs) == 2
+    v1, v2 = vs
+    pe = partitionedge(bpc, (v1, "bra") => (v2, "bra"))
+    mts = messages(bpc)
+    ind1, ind2 = noncommonind(singular_values![], ψ[v1]),
+    commonind(singular_values![], ψ[v1])
+    singular_values![] = denseblocks(replaceind(singular_values![], ind1, ind2'))
+    set!(mts, pe, ITensor[singular_values![]])
+    set!(mts, reverse(pe), ITensor[singular_values![]])
+  end
+  for v in vs
+    bpc = update_factor(bpc, (v, "ket"), ψ[v])
+    bpc = update_factor(bpc, (v, "bra"), ψdag[v])
+  end
+  return ψ, bpc
+end
+
+function get_local_term(bpc::BeliefPropagationCache, v)
+  qf = tensornetwork(bpc)
+  return qf[(v, "ket")] * qf[(v, "operator")] * qf[(v, "bra")]
+end
+
+function exact_energy(g::AbstractGraph, bpc::BeliefPropagationCache)
+  tn = ITensorNetwork(g)
+  for v in vertices(g)
+    tn[v] = get_local_term(bpc, v)
+  end
+  degree_two_sites = filter(v -> degree(tn, v) == 2, vertices(tn))
+  while !isempty(degree_two_sites)
+    v = first(degree_two_sites)
+    vn = first(neighbors(g, v))
+    tn = contract(tn, NamedEdge(v => vn); merged_vertex=vn)
+    degree_two_sites = filter(v -> degree(tn, v) == 2, vertices(tn))
+  end
+  return ITensors.contract(ITensor[tn[v] for v in vertices(tn)]; sequence="automatic")[]
+end
+
+function renamer(g)
+  vertex_rename = Dictionary()
+  for (i, v) in enumerate(vertices(g))
+    set!(vertex_rename, v, (i,))
+  end
+  return rename_vertices(v -> vertex_rename[v], g)
+end
+
+function imaginary_time_evo(
+  s::IndsNetwork,
+  ψ::ITensorNetwork,
+  model::Function,
+  dbetas::Vector{<:Tuple};
+  model_params,
+  bp_update_kwargs=(; maxiter=10, tol=1e-10),
+  apply_kwargs=(; cutoff=1e-12, maxdim=10),
+)
+  ψ = copy(ψ)
+  g = underlying_graph(ψ)
+
+  ℋ = model(g; model_params...)
+  ψψ = norm_sqr_network(ψ)
+  bpc = BeliefPropagationCache(ψψ, group(v -> v[1], vertices(ψψ)))
+  bpc = update(bpc; bp_update_kwargs...)
+  L = length(vertices(ψ))
+  println("Starting Imaginary Time Evolution")
+  β = 0
+  for (i, period) in enumerate(dbetas)
+    nbetas, dβ = first(period), last(period)
+    println("Entering evolution period $i , β = $β, dβ = $dβ")
+    U = exp(-dβ * ℋ; alg=Trotter{1}())
+    gates = Vector{ITensor}(U, s)
+    for i in 1:nbetas
+      for gate in gates
+        ψ, bpc = BP_apply(gate, ψ, bpc; apply_kwargs...)
+      end
+      β += dβ
+      bpc = update(bpc; bp_update_kwargs...)
+    end
+    e = sum(expect(ψ, ℋ; alg="bp"))
+    println("Energy is $(e / L)")
+  end
+
+  return ψ
+end