diff --git a/examples/belief_propagation/bpexample.jl b/examples/belief_propagation/bpexample.jl
index 6fc96ed1..a0af8818 100644
--- a/examples/belief_propagation/bpexample.jl
+++ b/examples/belief_propagation/bpexample.jl
@@ -4,6 +4,7 @@ using Metis
 using ITensorNetworks
 using Random
 using SplitApplyCombine
+using NamedGraphs
 
 using ITensorNetworks:
   belief_propagation,
@@ -34,7 +35,8 @@ function main()
     ψψ; subgraph_vertices=collect(values(group(v -> v[1], vertices(ψψ))))
   )
 
-  mts = belief_propagation(ψψ, mts; contract_kwargs=(; alg="exact"))
+  mts = belief_propagation(ψψ, mts; contract_kwargs=(; alg="exact"), niters=20)
+
   numerator_network = approx_network_region(
     ψψ, mts, [(v, 1)]; verts_tn=ITensorNetwork([apply(op("Sz", s[v]), ψ[v])])
   )
@@ -52,7 +54,7 @@ function main()
   )
   Zpp = partition(ψψ; subgraph_vertices=nested_graph_leaf_vertices(Zp))
   mts = message_tensors(Zpp)
-  mts = belief_propagation(ψψ, mts; contract_kwargs=(; alg="exact"))
+  mts = belief_propagation(ψψ, mts; contract_kwargs=(; alg="exact"), niters=20)
   numerator_network = approx_network_region(
     ψψ, mts, [(v, 1)]; verts_tn=ITensorNetwork([apply(op("Sz", s[v]), ψ[v])])
   )
diff --git a/examples/belief_propagation/bpsequences.jl b/examples/belief_propagation/bpsequences.jl
new file mode 100644
index 00000000..456e96c7
--- /dev/null
+++ b/examples/belief_propagation/bpsequences.jl
@@ -0,0 +1,81 @@
+using Compat
+using ITensors
+using Metis
+using ITensorNetworks
+using Random
+using SplitApplyCombine
+using Graphs
+using NamedGraphs
+
+using ITensorNetworks:
+  belief_propagation,
+  approx_network_region,
+  contract_inner,
+  message_tensors,
+  nested_graph_leaf_vertices,
+  edge_sequence
+
+function main()
+  g_labels = [
+    "Comb Tree",
+    "100 Site Random Regular Graph z = 3",
+    "6x6 Square Grid",
+    "4x4 Hexagonal Lattice",
+  ]
+  gs = [
+    named_comb_tree((6, 6)),
+    NamedGraph(Graphs.random_regular_graph(100, 3)),
+    named_grid((6, 6)),
+    NamedGraphs.hexagonal_lattice_graph(4, 4),
+  ]
+  χs = [4, 4, 2, 3]
+
+  for (i, g) in enumerate(gs)
+    Random.seed!(5467)
+    g_label = g_labels[i]
+    χ = χs[i]
+    s = siteinds("S=1/2", g)
+    ψ = randomITensorNetwork(s; link_space=χ)
+    ψψ = ψ ⊗ prime(dag(ψ); sites=[])
+
+    #Initial message tensors for BP
+    mts_init = message_tensors(
+      ψψ; subgraph_vertices=collect(values(group(v -> v[1], vertices(ψψ))))
+    )
+
+    println("\nFirst testing out a $g_label. Random network with bond dim $χ")
+
+    #Now test out various sequences
+    print("Parallel updates (sequence is irrelevant): ")
+    belief_propagation(
+      ψψ,
+      mts_init;
+      contract_kwargs=(; alg="exact"),
+      target_precision=1e-10,
+      niters=100,
+      edges=edge_sequence(mts_init; alg="parallel"),
+      verbose=true,
+    )
+    print("Sequential updates (sequence is default edge list of the message tensors): ")
+    belief_propagation(
+      ψψ,
+      mts_init;
+      contract_kwargs=(; alg="exact"),
+      target_precision=1e-10,
+      niters=100,
+      edges=[e for e in edges(mts_init)],
+      verbose=true,
+    )
+    print("Sequential updates (sequence is our custom sequence finder): ")
+    belief_propagation(
+      ψψ,
+      mts_init;
+      contract_kwargs=(; alg="exact"),
+      target_precision=1e-10,
+      niters=100,
+      verbose=true,
+    )
+  end
+end
+
+main()
diff --git a/examples/dynamics/heavy_hex_ising_real_tebd.jl b/examples/dynamics/heavy_hex_ising_real_tebd.jl
index 4ce55b20..e8dcc863 100644
--- a/examples/dynamics/heavy_hex_ising_real_tebd.jl
+++ b/examples/dynamics/heavy_hex_ising_real_tebd.jl
@@ -24,10 +24,10 @@ end
 function ibm_processor_graph(n::Int64, m::Int64)
   g = heavy_hex_lattice_graph(n, m)
   dims = maximum(vertices(hexagonal_lattice_graph(n, m)))
-  v1, v2 = (1, dims[2]), (dims[1], 1)
+  v1, v2 = (dims[1], 1), (1, dims[2])
   add_vertices!(g, [v1, v2])
-  add_edge!(g, v1 => v1 .- (0, 1))
-  add_edge!(g, v2 => v2 .+ (0, 1))
+  add_edge!(g, v1 => v1 .- (1, 0))
+  add_edge!(g, v2 => v2 .+ (1, 0))
 
   return g
 end
diff --git a/examples/gauging/gauging_itns.jl b/examples/gauging/gauging_itns.jl
index 6c30cfd4..7463d54a 100644
--- a/examples/gauging/gauging_itns.jl
+++ b/examples/gauging/gauging_itns.jl
@@ -4,7 +4,6 @@ using Metis
 using ITensorNetworks
 using Random
 using SplitApplyCombine
-using ProfileView
 
 using ITensorNetworks:
   message_tensors,
@@ -18,7 +17,8 @@ using ITensorNetworks:
   vidal_to_symmetric_gauge,
   initialize_bond_tensors,
   vidal_itn_isometries,
-  norm_network
+  norm_network,
+  edge_sequence
 
 using NamedGraphs
 using NamedGraphs: add_edges!, rem_vertex!, hexagonal_lattice_graph
@@ -45,7 +45,10 @@ end
 
 """Bring an ITN into the Vidal gauge, various methods possible. Result is timed"""
 function benchmark_state_gauging(
-  ψ::ITensorNetwork; mode="BeliefPropagation", no_iterations=50
+  ψ::ITensorNetwork;
+  mode="belief_propagation",
+  no_iterations=50,
+  BP_update_order::String="sequential",
 )
   s = siteinds(ψ)
 
@@ -65,12 +68,19 @@ function benchmark_state_gauging(
   for i in 1:no_iterations
     println("On Iteration " * string(i))
 
-    if mode == "BeliefPropagation"
-      times_iters[i] = @elapsed mts, _ = belief_propagation_iteration(
-        ψψ, mts; contract_kwargs=(; alg="exact")
-      )
+    if mode == "belief_propagation"
+      if BP_update_order != "parallel"
+        times_iters[i] = @elapsed mts, _ = belief_propagation_iteration(
+          ψψ, mts; contract_kwargs=(; alg="exact")
+        )
+      else
+        times_iters[i] = @elapsed mts, _ = belief_propagation_iteration(
+          ψψ, mts; contract_kwargs=(; alg="exact"), edges=edge_sequence(mts; alg="parallel")
+        )
+      end
+
       times_gauging[i] = @elapsed ψ, bond_tensors = vidal_gauge(ψinit, mts)
-    elseif mode == "Eager"
+    elseif mode == "eager"
       times_iters[i] = @elapsed ψ, bond_tensors, mts = eager_gauging(ψ, bond_tensors, mts)
     else
       times_iters[i] = @elapsed begin
@@ -82,7 +92,7 @@ function benchmark_state_gauging(
 
     C[i] = vidal_itn_canonicalness(ψ, bond_tensors)
   end
-
+  @show times_iters, time
   simulation_times = cumsum(times_iters) + times_gauging
 
   return simulation_times, C
@@ -94,23 +104,41 @@ s = siteinds("S=1/2", g)
 ψ = randomITensorNetwork(s; link_space=χ)
 no_iterations = 30
 
-BPG_simulation_times, BPG_Cs = benchmark_state_gauging(ψ; no_iterations)
-Eager_simulation_times, Eager_Cs = benchmark_state_gauging(ψ; mode="Eager", no_iterations)
+BPG_simulation_times, BPG_Cs = benchmark_state_gauging(
+  ψ; no_iterations, BP_update_order="parallel"
+)
+BPG_sequential_simulation_times, BPG_sequential_Cs = benchmark_state_gauging(
+  ψ; no_iterations
+)
+Eager_simulation_times, Eager_Cs = benchmark_state_gauging(ψ; mode="eager", no_iterations)
 SU_simulation_times, SU_Cs = benchmark_state_gauging(ψ; mode="SU", no_iterations)
 
-epsilon = 1e-6
+epsilon = 1e-10
+
 println(
-  "Time for BPG to reach C < epsilon was " *
+  "Time for BPG (with parallel updates) to reach C < epsilon was " *
   string(BPG_simulation_times[findfirst(x -> x < 0, BPG_Cs .- epsilon)]) *
-  " seconds",
+  " seconds. No iters was " *
+  string(findfirst(x -> x < 0, BPG_Cs .- epsilon)),
+)
+println(
+  "Time for BPG (with sequential updates) to reach C < epsilon was " *
+  string(
+    BPG_sequential_simulation_times[findfirst(x -> x < 0, BPG_sequential_Cs .- epsilon)]
+  ) *
+  " seconds. No iters was " *
+  string(findfirst(x -> x < 0, BPG_sequential_Cs .- epsilon)),
 )
+
 println(
-  "Time for Eager to reach C < epsilon was " *
+  "Time for Eager Gauging to reach C < epsilon was " *
   string(Eager_simulation_times[findfirst(x -> x < 0, Eager_Cs .- epsilon)]) *
-  " seconds",
+  " seconds. No iters was " *
+  string(findfirst(x -> x < 0, Eager_Cs .- epsilon)),
 )
 println(
-  "Time for SU to reach C < epsilon was " *
+  "Time for SU Gauging (with sequential updates) to reach C < epsilon was " *
   string(SU_simulation_times[findfirst(x -> x < 0, SU_Cs .- epsilon)]) *
-  " seconds",
+  " seconds. No iters was " *
+  string(findfirst(x -> x < 0, SU_Cs .- epsilon)),
 )
diff --git a/src/ITensorNetworks.jl b/src/ITensorNetworks.jl
index 55f261ed..5b660bc9 100644
--- a/src/ITensorNetworks.jl
+++ b/src/ITensorNetworks.jl
@@ -97,6 +97,7 @@ include("specialitensornetworks.jl")
 include("renameitensornetwork.jl")
 include("boundarymps.jl")
 include(joinpath("beliefpropagation", "beliefpropagation.jl"))
+include(joinpath("beliefpropagation", "beliefpropagation_schedule.jl"))
 include(joinpath("beliefpropagation", "sqrt_beliefpropagation.jl"))
 include("contraction_tree_to_graph.jl")
 include("gauging.jl")
diff --git a/src/beliefpropagation/beliefpropagation.jl b/src/beliefpropagation/beliefpropagation.jl
index dc23d4a5..8a1f17b4 100644
--- a/src/beliefpropagation/beliefpropagation.jl
+++ b/src/beliefpropagation/beliefpropagation.jl
@@ -22,12 +22,14 @@ function message_tensors_skeleton(subgraphs::DataGraph)
 end
 
 function message_tensors(
-  subgraphs::DataGraph; itensor_constructor=inds_e -> dense(delta(inds_e))
+  subgraphs::DataGraph;
+  itensor_constructor=inds_e -> ITensor[dense(delta(i)) for i in inds_e],
 )
   mts = message_tensors_skeleton(subgraphs)
   for e in edges(subgraphs)
     inds_e = commoninds(subgraphs[src(e)], subgraphs[dst(e)])
-    mts[e] = ITensorNetwork(map(itensor_constructor, inds_e))
+    itensors = itensor_constructor(inds_e)
+    mts[e] = ITensorNetwork(itensors)
     mts[reverse(e)] = dag(mts[e])
   end
   return mts
@@ -74,24 +76,24 @@ function update_message_tensor(
 end
 
 """
-Do an update of all message tensors for a given ITensornetwork and its partition into sub graphs
+Do a sequential update of message tensors on `edges` for a given ITensornetwork and its partition into sub graphs
 """
 function belief_propagation_iteration(
   tn::ITensorNetwork,
-  mts::DataGraph;
+  mts::DataGraph,
+  edges::Vector{<:AbstractEdge};
   contract_kwargs=(; alg="density_matrix", output_structure=path_graph_structure, maxdim=1),
   compute_norm=false,
 )
   new_mts = copy(mts)
   c = 0
-  es = edges(mts)
-  for e in es
+  for e in edges
     environment_tensornetworks = ITensorNetwork[
-      mts[e_in] for e_in in setdiff(boundary_edges(mts, [src(e)]; dir=:in), [reverse(e)])
+      new_mts[e_in] for
+      e_in in setdiff(boundary_edges(new_mts, [src(e)]; dir=:in), [reverse(e)])
     ]
-
     new_mts[src(e) => dst(e)] = update_message_tensor(
-      tn, mts[src(e)], environment_tensornetworks; contract_kwargs
+      tn, new_mts[src(e)], environment_tensornetworks; contract_kwargs
     )
 
     if compute_norm
@@ -102,25 +104,64 @@ function belief_propagation_iteration(
       c += 0.5 * norm(denseblocks(LHS) - denseblocks(RHS))
     end
   end
-  return new_mts, c / (length(es))
+  return new_mts, c / (length(edges))
+end
+
+"""
+Do parallel updates between groups of edges of all message tensors for a given ITensornetwork and its partition into sub graphs.
+Currently we send the full message tensor data struct to belief_propagation_iteration for each subgraph. But really we only need the
+mts relevant to that subgraph.
+"""
+function belief_propagation_iteration(
+  tn::ITensorNetwork,
+  mts::DataGraph,
+  edge_groups::Vector{<:Vector{<:AbstractEdge}};
+  contract_kwargs=(; alg="density_matrix", output_structure=path_graph_structure, maxdim=1),
+  compute_norm=false,
+)
+  new_mts = copy(mts)
+  c = 0
+  for edges in edge_groups
+    updated_mts, ct = belief_propagation_iteration(
+      tn, mts, edges; contract_kwargs, compute_norm
+    )
+    for e in edges
+      new_mts[e] = updated_mts[e]
+    end
+    c += ct
+  end
+  return new_mts, c / (length(edge_groups))
+end
+
+function belief_propagation_iteration(
+  tn::ITensorNetwork,
+  mts::DataGraph;
+  contract_kwargs=(; alg="density_matrix", output_structure=path_graph_structure, maxdim=1),
+  compute_norm=false,
+  edges=edge_sequence(mts),
+)
+  return belief_propagation_iteration(tn, mts, edges; contract_kwargs, compute_norm)
 end
 
 function belief_propagation(
   tn::ITensorNetwork,
   mts::DataGraph;
   contract_kwargs=(; alg="density_matrix", output_structure=path_graph_structure, maxdim=1),
-  niters=20,
-  target_precision::Union{Float64,Nothing}=nothing,
+  niters=default_bp_niters(mts),
+  target_precision=nothing,
+  edges=edge_sequence(mts),
+  verbose=false,
 )
-  compute_norm = target_precision == nothing ? false : true
+  compute_norm = !isnothing(target_precision)
+  if isnothing(niters)
+    error("You need to specify a number of iterations for BP!")
+  end
   for i in 1:niters
-    mts, c = belief_propagation_iteration(tn, mts; contract_kwargs, compute_norm)
+    mts, c = belief_propagation_iteration(tn, mts, edges; contract_kwargs, compute_norm)
     if compute_norm && c <= target_precision
-      println(
-        "Belief Propagation finished. Reached a canonicalness of " *
-        string(c) *
-        " after $i iterations. ",
-      )
+      if verbose
+        println("BP converged to desired precision after $i iterations.")
+      end
       break
     end
   end
@@ -133,11 +174,12 @@ function belief_propagation(
   nvertices_per_partition=nothing,
   npartitions=nothing,
   subgraph_vertices=nothing,
-  niters=20,
-  target_precision::Union{Float64,Nothing}=nothing,
+  niters=default_bp_niters(mts),
+  target_precision=nothing,
+  verbose=false,
 )
   mts = message_tensors(tn; nvertices_per_partition, npartitions, subgraph_vertices)
-  return belief_propagation(tn, mts; contract_kwargs, niters, target_precision)
+  return belief_propagation(tn, mts; contract_kwargs, niters, target_precision, verbose)
 end
 
 """
diff --git a/src/beliefpropagation/beliefpropagation_schedule.jl b/src/beliefpropagation/beliefpropagation_schedule.jl
new file mode 100644
index 00000000..cbc5f206
--- /dev/null
+++ b/src/beliefpropagation/beliefpropagation_schedule.jl
@@ -0,0 +1,40 @@
+default_edge_sequence_alg() = "forest_cover"
+
+@traitfn default_bp_niters(g::::(!IsDirected)) = is_tree(g) ? 1 : nothing
+@traitfn function default_bp_niters(g::::IsDirected)
+  return default_bp_niters(undirected_graph(underlying_graph(g)))
+end
+
+@traitfn function edge_sequence(
+  g::::(!IsDirected); alg=default_edge_sequence_alg(), kwargs...
+)
+  return edge_sequence(Algorithm(alg), g; kwargs...)
+end
+
+@traitfn function edge_sequence(g::::IsDirected; alg=default_edge_sequence_alg(), kwargs...)
+  return edge_sequence(Algorithm(alg), undirected_graph(underlying_graph(g)); kwargs...)
+end
+
+@traitfn function edge_sequence(alg::Algorithm, g::::IsDirected; kwargs...)
+  return edge_sequence(alg, undirected_graph(underlying_graph(g)); kwargs...)
+end
+
+@traitfn function edge_sequence(
+  ::Algorithm"forest_cover", g::::(!IsDirected); root_vertex=NamedGraphs.default_root_vertex
+)
+  forests = NamedGraphs.forest_cover(g)
+  edges = edgetype(g)[]
+  for forest in forests
+    trees = [forest[vs] for vs in connected_components(forest)]
+    for tree in trees
+      tree_edges = post_order_dfs_edges(tree, root_vertex(tree))
+      push!(edges, vcat(tree_edges, reverse(reverse.(tree_edges)))...)
+    end
+  end
+
+  return edges
+end
+
+@traitfn function edge_sequence(::Algorithm"parallel", g::::(!IsDirected))
+  return [[e] for e in vcat(edges(g), reverse.(edges(g)))]
+end
diff --git a/src/beliefpropagation/sqrt_beliefpropagation.jl b/src/beliefpropagation/sqrt_beliefpropagation.jl
index 15662367..20ff3d8e 100644
--- a/src/beliefpropagation/sqrt_beliefpropagation.jl
+++ b/src/beliefpropagation/sqrt_beliefpropagation.jl
@@ -1,44 +1,19 @@
 # using ITensors: scalartype
 # using ITensorNetworks: find_subgraph, map_diag, sqrt_diag, boundary_edges
 
-function sqrt_belief_propagation(
-  tn::ITensorNetwork,
-  mts::DataGraph;
-  niters=20,
-  # target_precision::Union{Float64,Nothing}=nothing,
-)
-  # compute_norm = target_precision == nothing ? false : true
-  sqrt_mts = sqrt_message_tensors(tn, mts)
-  for i in 1:niters
-    sqrt_mts, c = sqrt_belief_propagation_iteration(tn, sqrt_mts) #; compute_norm)
-    # if compute_norm && c <= target_precision
-    #   println(
-    #     "Belief Propagation finished. Reached a canonicalness of " *
-    #     string(c) *
-    #     " after $i iterations. ",
-    #   )
-    #   break
-    # end
-  end
-  return sqr_message_tensors(sqrt_mts)
-end
-
 function sqrt_belief_propagation_iteration(
-  tn::ITensorNetwork,
-  sqrt_mts::DataGraph;
-  # compute_norm=false,
+  tn::ITensorNetwork, sqrt_mts::DataGraph, edges::Vector{<:AbstractEdge}
 )
   new_sqrt_mts = copy(sqrt_mts)
   c = 0.0
-  es = edges(sqrt_mts)
-  for e in es
+  for e in edges
     environment_tensornetworks = ITensorNetwork[
-      sqrt_mts[e_in] for
-      e_in in setdiff(boundary_edges(sqrt_mts, [src(e)]; dir=:in), [reverse(e)])
+      new_sqrt_mts[e_in] for
+      e_in in setdiff(boundary_edges(new_sqrt_mts, [src(e)]; dir=:in), [reverse(e)])
     ]
 
     new_sqrt_mts[src(e) => dst(e)] = update_sqrt_message_tensor(
-      tn, sqrt_mts[src(e)], environment_tensornetworks;
+      tn, new_sqrt_mts[src(e)], environment_tensornetworks;
     )
 
     # if compute_norm
@@ -49,7 +24,54 @@ function sqrt_belief_propagation_iteration(
     #   c += 0.5 * norm(LHS - RHS)
     # end
   end
-  return new_sqrt_mts, c / (length(es))
+  return new_sqrt_mts, c / (length(edges))
+end
+
+function sqrt_belief_propagation_iteration(
+  tn::ITensorNetwork, sqrt_mts::DataGraph, edges::Vector{<:Vector{<:AbstractEdge}}
+)
+  new_sqrt_mts = copy(sqrt_mts)
+  c = 0.0
+  for e_group in edges
+    updated_sqrt_mts, ct = sqrt_belief_propagation_iteration(tn, sqr_mts, e_group)
+    for e in e_group
+      new_sqrt_mts[e] = updated_sqrt_mts[e]
+    end
+    c += ct
+  end
+  return new_sqrt_mts, c / (length(edges))
+end
+
+function sqrt_belief_propagation_iteration(
+  tn::ITensorNetwork, sqrt_mts::DataGraph; edges=edge_sequence(mts)
+)
+  return sqrt_belief_propagation_iteration(tn, sqrt_mts, edges)
+end
+
+function sqrt_belief_propagation(
+  tn::ITensorNetwork,
+  mts::DataGraph;
+  niters=default_bp_niters(mts),
+  edges=edge_sequence(mts),
+  # target_precision::Union{Float64,Nothing}=nothing,
+)
+  # compute_norm = target_precision == nothing ? false : true
+  sqrt_mts = sqrt_message_tensors(tn, mts)
+  if isnothing(niters)
+    error("You need to specify a number of iterations for BP!")
+  end
+  for i in 1:niters
+    sqrt_mts, c = sqrt_belief_propagation_iteration(tn, sqrt_mts, edges) #; compute_norm)
+    # if compute_norm && c <= target_precision
+    #   println(
+    #     "Belief Propagation finished. Reached a canonicalness of " *
+    #     string(c) *
+    #     " after $i iterations. ",
+    #   )
+    #   break
+    # end
+  end
+  return sqr_message_tensors(sqrt_mts)
 end
 
 function update_sqrt_message_tensor(
diff --git a/src/gauging.jl b/src/gauging.jl
index 9b321a6e..00347046 100644
--- a/src/gauging.jl
+++ b/src/gauging.jl
@@ -17,11 +17,12 @@ function vidal_gauge(
   bond_tensors::DataGraph;
   eigen_message_tensor_cutoff=10 * eps(real(scalartype(ψ))),
   regularization=10 * eps(real(scalartype(ψ))),
+  edges=NamedGraphs.edges(ψ),
   svd_kwargs...,
 )
   ψ_vidal = copy(ψ)
 
-  for e in edges(ψ_vidal)
+  for e in edges
     vsrc, vdst = src(e), dst(e)
     ψvsrc, ψvdst = ψ_vidal[vsrc], ψ_vidal[vdst]
 
@@ -79,11 +80,12 @@ function vidal_gauge(
   mts::DataGraph;
   eigen_message_tensor_cutoff=10 * eps(real(scalartype(ψ))),
   regularization=10 * eps(real(scalartype(ψ))),
+  edges=NamedGraphs.edges(ψ),
   svd_kwargs...,
 )
   bond_tensors = initialize_bond_tensors(ψ)
   return vidal_gauge(
-    ψ, mts, bond_tensors; eigen_message_tensor_cutoff, regularization, svd_kwargs...
+    ψ, mts, bond_tensors; eigen_message_tensor_cutoff, regularization, edges, svd_kwargs...
   )
 end
 
@@ -94,6 +96,7 @@ function vidal_gauge(
   regularization=10 * eps(real(scalartype(ψ))),
   niters=30,
   target_canonicalness::Union{Nothing,Float64}=nothing,
+  verbose=false,
   svd_kwargs...,
 )
   ψψ = norm_network(ψ)
@@ -101,7 +104,12 @@ function vidal_gauge(
   mts = message_tensors(Z)
 
   mts = belief_propagation(
-    ψψ, mts; contract_kwargs=(; alg="exact"), niters, target_precision=target_canonicalness
+    ψψ,
+    mts;
+    contract_kwargs=(; alg="exact"),
+    niters,
+    target_precision=target_canonicalness,
+    verbose,
   )
   return vidal_gauge(ψ, mts; eigen_message_tensor_cutoff, regularization, svd_kwargs...)
 end
@@ -171,10 +179,14 @@ function symmetric_to_vidal_gauge(
 end
 
 """Function to measure the 'isometries' of a state in the Vidal Gauge"""
-function vidal_itn_isometries(ψ::ITensorNetwork, bond_tensors::DataGraph)
+function vidal_itn_isometries(
+  ψ::ITensorNetwork,
+  bond_tensors::DataGraph;
+  edges=vcat(NamedGraphs.edges(ψ), reverse.(NamedGraphs.edges(ψ))),
+)
   isometries = DataGraph{vertextype(ψ),ITensor,ITensor}(directed_graph(underlying_graph(ψ)))
 
-  for e in vcat(edges(ψ), reverse.(edges(ψ)))
+  for e in edges
     vsrc, vdst = src(e), dst(e)
     ψv = copy(ψ[vsrc])
     for vn in setdiff(neighbors(ψ, vsrc), [vdst])
diff --git a/src/imports.jl b/src/imports.jl
index 3186de4d..71406878 100644
--- a/src/imports.jl
+++ b/src/imports.jl
@@ -96,6 +96,8 @@ import ITensors:
   scalartype,
   #adding
   add
+#Algorithm
+Algorithm
 
 using ITensors.ContractionSequenceOptimization: deepmap
 
diff --git a/test/test_belief_propagation.jl b/test/test_belief_propagation.jl
index 6bc1eb01..d3747598 100644
--- a/test/test_belief_propagation.jl
+++ b/test/test_belief_propagation.jl
@@ -22,7 +22,7 @@ ITensors.disable_warn_order()
 
 @testset "belief_propagation" begin
 
-  #FIRST TEST SINGLE SITE ON AN MPS, SHOULD BE EXACT
+  #First test on an MPS, should be exact
   g_dims = (1, 6)
   g = named_grid(g_dims)
   s = siteinds("S=1/2", g)
@@ -40,18 +40,44 @@ ITensors.disable_warn_order()
 
   Z = partition(ψψ; subgraph_vertices=collect(values(group(v -> v[1], vertices(ψψ)))))
   mts = message_tensors(Z)
-  mts = belief_propagation(ψψ, mts; contract_kwargs=(; alg="exact"), target_precision=1e-8)
+  mts = belief_propagation(ψψ, mts; contract_kwargs=(; alg="exact"))
+
+  numerator_network = approx_network_region(
+    ψψ, mts, [(v, 1)]; verts_tn=ITensorNetwork(ITensor[apply(op("Sz", s[v]), ψ[v])])
+  )
+  denominator_network = approx_network_region(ψψ, mts, [(v, 1)])
+  bp_sz = ITensors.contract(numerator_network)[] / ITensors.contract(denominator_network)[]
+
+  @test abs.(bp_sz - exact_sz) <= 1e-14
+
+  #Now test on a tree, should also be exact
+  g = named_comb_tree((6, 6))
+  s = siteinds("S=1/2", g)
+  χ = 2
+  Random.seed!(1564)
+  ψ = randomITensorNetwork(s; link_space=χ)
+
+  ψψ = ψ ⊗ prime(dag(ψ); sites=[])
+
+  v = (1, 3)
+
+  Oψ = copy(ψ)
+  Oψ[v] = apply(op("Sz", s[v]), ψ[v])
+  exact_sz = contract_inner(Oψ, ψ) / contract_inner(ψ, ψ)
+
+  Z = partition(ψψ; subgraph_vertices=collect(values(group(v -> v[1], vertices(ψψ)))))
+  mts = message_tensors(Z)
+  mts = belief_propagation(ψψ, mts; contract_kwargs=(; alg="exact"))
 
   numerator_network = approx_network_region(
     ψψ, mts, [(v, 1)]; verts_tn=ITensorNetwork(ITensor[apply(op("Sz", s[v]), ψ[v])])
   )
   denominator_network = approx_network_region(ψψ, mts, [(v, 1)])
-  bp_sz = contract(numerator_network)[] / contract(denominator_network)[]
+  bp_sz = ITensors.contract(numerator_network)[] / ITensors.contract(denominator_network)[]
 
   @test abs.(bp_sz - exact_sz) <= 1e-14
 
-  #NOW TEST TWO_SITE_EXPEC TAKING ON THE PARTITION FUNCTION OF THE RECTANGULAR ISING. SHOULD BE REASONABLE AND 
-  #INDEPENDENT OF INIT CONDITIONS, FOR SMALL BETA
+  #Now test two-site expec taking on the partition function of the Ising model. Not exact, but close
   g_dims = (3, 4)
   g = named_grid(g_dims)
   s = IndsNetwork(g; link_space=2)
@@ -68,17 +94,18 @@ ITensors.disable_warn_order()
   nsites = 2
   Z = partition(ψψ; nvertices_per_partition=nsites)
   mts = message_tensors(Z)
-  mts = belief_propagation(ψψ, mts; contract_kwargs=(; alg="exact"))
+  mts = belief_propagation(ψψ, mts; contract_kwargs=(; alg="exact"), niters=20)
   numerator_network = approx_network_region(
     ψψ, mts, vs; verts_tn=ITensorNetwork(ITensor[ψOψ[v] for v in vs])
   )
 
   denominator_network = approx_network_region(ψψ, mts, vs)
-  bp_szsz = contract(numerator_network)[] / contract(denominator_network)[]
+  bp_szsz =
+    ITensors.contract(numerator_network)[] / ITensors.contract(denominator_network)[]
 
   @test abs.(bp_szsz - actual_szsz) <= 0.05
 
-  #TEST FORMING OF A TWO SITE RDM. JUST CHECK THAT IS HAS THE CORRECT SIZE, TRACE AND IS PSD
+  #Test forming a two-site RDM. Check it has the correct size, trace 1 and is PSD
   g_dims = (4, 4)
   g = named_grid(g_dims)
   s = siteinds("S=1/2", g)
@@ -93,10 +120,10 @@ ITensors.disable_warn_order()
   )
   Zpp = partition(ψψ; subgraph_vertices=nested_graph_leaf_vertices(Zp))
   mts = message_tensors(Zpp)
-  mts = belief_propagation(ψψ, mts; contract_kwargs=(; alg="exact"))
+  mts = belief_propagation(ψψ, mts; contract_kwargs=(; alg="exact"), niters=20)
 
   ψψsplit = split_index(ψψ, NamedEdge.([(v, 1) => (v, 2) for v in vs]))
-  rdm = contract(
+  rdm = ITensors.contract(
     approx_network_region(
       ψψ,
       mts,
@@ -112,7 +139,7 @@ ITensors.disable_warn_order()
   @test all(>=(0), real(eigs)) && all(==(0), imag(eigs))
   @test size(rdm) == (2^length(vs), 2^length(vs))
 
-  #TEST MORE ADVANCED BP WITH ITENSORNETWORK MESSAGE TENSORS. IN THIS CASE IT SHOULD BE LIKE BOUNDARY MPS
+  #Test more advanced block BP with MPS message tensors on a grid 
   g_dims = (5, 4)
   g = named_grid(g_dims)
   s = siteinds("S=1/2", g)
@@ -143,7 +170,7 @@ ITensors.disable_warn_order()
   numerator_network = approx_network_region(ψψ, mts, [v]; verts_tn=ITensorNetwork(ψOψ[v]))
 
   denominator_network = approx_network_region(ψψ, mts, [v])
-  bp_sz = contract(numerator_network)[] / contract(denominator_network)[]
+  bp_sz = ITensors.contract(numerator_network)[] / ITensors.contract(denominator_network)[]
 
   exact_sz =
     contract_boundary_mps(ψOψ; cutoff=1e-16) / contract_boundary_mps(ψψ; cutoff=1e-16)
diff --git a/test/test_examples/test_examples.jl b/test/test_examples/test_examples.jl
index d4532451..8afc3be4 100644
--- a/test/test_examples/test_examples.jl
+++ b/test/test_examples/test_examples.jl
@@ -14,7 +14,7 @@ using Test
     "peps.jl",
     "steiner_tree.jl",
     joinpath("belief_propagation", "bpexample.jl"),
-    joinpath("belief_propagation", "sqrt_bp.jl"),
+    joinpath("belief_propagation", "bpsequences.jl"),
     joinpath("dynamics", "2d_ising_imag_tebd.jl"),
     joinpath("dynamics", "heavy_hex_ising_real_tebd.jl"),
     joinpath("treetensornetworks", "comb_tree.jl"),
diff --git a/test/test_gauging.jl b/test/test_gauging.jl
index b8d9fb5f..5d5e6c77 100644
--- a/test/test_gauging.jl
+++ b/test/test_gauging.jl
@@ -21,7 +21,7 @@ using SplitApplyCombine
   dims = (n, n)
   g = named_grid(dims)
   s = siteinds("S=1/2", g)
-  χ = 3
+  χ = 6
 
   Random.seed!(5467)
   ψ = randomITensorNetwork(s; link_space=χ)
diff --git a/test/test_treetensornetworks/test_solvers/test_dmrg.jl b/test/test_treetensornetworks/test_solvers/test_dmrg.jl
index 61c2febc..939b4fee 100644
--- a/test/test_treetensornetworks/test_solvers/test_dmrg.jl
+++ b/test/test_treetensornetworks/test_solvers/test_dmrg.jl
@@ -28,7 +28,7 @@ using Observers
   # Compare to `ITensors.MPO` version of `dmrg`
   H_mpo = MPO([H[v] for v in 1:nv(H)])
   psi_mps = MPS([psi[v] for v in 1:nv(psi)])
-  e2, psi2 = dmrg(H_mpo, psi_mps; nsweeps, maxdim, normalize=false, outputlevel=0)
+  e2, psi2 = dmrg(H_mpo, psi_mps; nsweeps, maxdim, outputlevel=0)
 
   psi = dmrg(H, psi; nsweeps, maxdim, cutoff, nsite, solver_krylovdim=3, solver_maxiter=1)
   @test inner(psi', H, psi) ≈ inner(psi2', H_mpo, psi2)
@@ -134,7 +134,7 @@ end
   sline = only.(collect(vertex_data(s)))[linear_order]
   Hline = MPO(relabel_sites(os, vmap), sline)
   psiline = randomMPS(sline; linkdims=20)
-  e2, psi2 = dmrg(Hline, psiline, sweeps; normalize=false, outputlevel=0)
+  e2, psi2 = dmrg(Hline, psiline, sweeps; outputlevel=0)
 
   @test inner(psi', H, psi) ≈ inner(psi2', Hline, psi2) atol = 1e-5
 end