From c8459471e67c4163539808b8d432032bdde93838 Mon Sep 17 00:00:00 2001
From: Joey <joeytindall94@googlemail.com>
Date: Fri, 13 Sep 2024 10:19:49 -0400
Subject: [PATCH 01/14] Blah

---
 examples/test.jl | 21 +++++++++++++++++++++
 1 file changed, 21 insertions(+)
 create mode 100644 examples/test.jl

diff --git a/examples/test.jl b/examples/test.jl
new file mode 100644
index 00000000..6bf1bfcf
--- /dev/null
+++ b/examples/test.jl
@@ -0,0 +1,21 @@
+using ITensorNetworks: IndsNetwork, siteinds, ttn
+using ITensorNetworks.ModelHamiltonians: ising
+using ITensors: Index, OpSum, terms, sites
+using NamedGraphs.NamedGraphGenerators: named_grid
+using NamedGraphs.GraphsExtensions: rem_vertex
+
+function filter_terms(H, verts)
+    H_new = OpSum()
+    for term in terms(H)
+        if isempty(filter(v -> v ∈ verts, sites(term)))
+            H_new += term
+        end
+    end
+    return H_new
+end
+
+g = named_grid((8,1))
+s = siteinds("S=1/2", g)
+H = ising(s)
+H_mod = filter_terms(H, [(4,1)])
+ttno = ttn(H_mod, s)
\ No newline at end of file

From 6ff0cd572c947e9b1ed3642e690b43233277beb0 Mon Sep 17 00:00:00 2001
From: Joey <joeytindall94@googlemail.com>
Date: Thu, 17 Oct 2024 14:56:22 +0100
Subject: [PATCH 02/14] Bug fix in current ortho. Change test

---
 .../alternating_update/region_update.jl       | 45 ++++++++-----------
 .../test_solvers/test_dmrg.jl                 | 12 ++---
 2 files changed, 25 insertions(+), 32 deletions(-)

diff --git a/src/solvers/alternating_update/region_update.jl b/src/solvers/alternating_update/region_update.jl
index b92adc8c..97241c20 100644
--- a/src/solvers/alternating_update/region_update.jl
+++ b/src/solvers/alternating_update/region_update.jl
@@ -7,36 +7,27 @@ function current_ortho(sweep_plan, which_region_update)
   if !isa(region, AbstractEdge) && length(region) == 1
     return only(current_verts)
   end
-  if which_region_update == length(regions)
-    # look back by one should be sufficient, but may be brittle?
-    overlapping_vertex = only(
-      intersect(current_verts, support(regions[which_region_update - 1]))
-    )
-    return overlapping_vertex
-  else
-    # look forward
-    other_regions = filter(
-      x -> !(issetequal(x, current_verts)), support.(regions[(which_region_update + 1):end])
+  # look forward
+  other_regions = filter(
+    x -> !(issetequal(x, current_verts)), support.(regions[(which_region_update + 1):end])
+  )
+  # find the first region that has overlapping support with current region 
+  ind = findfirst(x -> !isempty(intersect(support(x), support(region))), other_regions)
+  if isnothing(ind)
+    # look backward
+    other_regions = reverse(
+      filter(
+        x -> !(issetequal(x, current_verts)), support.(regions[1:(which_region_update - 1)])
+      ),
     )
-    # find the first region that has overlapping support with current region 
     ind = findfirst(x -> !isempty(intersect(support(x), support(region))), other_regions)
-    if isnothing(ind)
-      # look backward
-      other_regions = reverse(
-        filter(
-          x -> !(issetequal(x, current_verts)),
-          support.(regions[1:(which_region_update - 1)]),
-        ),
-      )
-      ind = findfirst(x -> !isempty(intersect(support(x), support(region))), other_regions)
-    end
-    @assert !isnothing(ind)
-    future_verts = union(support(other_regions[ind]))
-    # return ortho_ceter as the vertex in current region that does not overlap with following one
-    overlapping_vertex = intersect(current_verts, future_verts)
-    nonoverlapping_vertex = only(setdiff(current_verts, overlapping_vertex))
-    return nonoverlapping_vertex
   end
+  @assert !isnothing(ind)
+  future_verts = union(support(other_regions[ind]))
+  # return ortho_ceter as the vertex in current region that does not overlap with following one
+  overlapping_vertex = intersect(current_verts, future_verts)
+  nonoverlapping_vertex = only(setdiff(current_verts, overlapping_vertex))
+  return nonoverlapping_vertex
 end
 
 function region_update(
diff --git a/test/test_treetensornetworks/test_solvers/test_dmrg.jl b/test/test_treetensornetworks/test_solvers/test_dmrg.jl
index cf8a1caf..004ec561 100644
--- a/test/test_treetensornetworks/test_solvers/test_dmrg.jl
+++ b/test/test_treetensornetworks/test_solvers/test_dmrg.jl
@@ -1,7 +1,7 @@
 @eval module $(gensym())
 using DataGraphs: edge_data, vertex_data
 using Dictionaries: Dictionary
-using Graphs: nv, vertices
+using Graphs: nv, vertices, uniform_tree
 using ITensorMPS: ITensorMPS
 using ITensorNetworks:
   ITensorNetworks,
@@ -19,6 +19,7 @@ using ITensorNetworks.ITensorsExtensions: replace_vertices
 using ITensorNetworks.ModelHamiltonians: ModelHamiltonians
 using ITensors: ITensors
 using KrylovKit: eigsolve
+using NamedGraphs: NamedGraph, rename_vertices
 using NamedGraphs.NamedGraphGenerators: named_comb_tree
 using Observers: observer
 using StableRNGs: StableRNG
@@ -313,11 +314,12 @@ end
   nsites = 2
   nsweeps = 10
 
-  c = named_comb_tree((3, 2))
-  s = siteinds("S=1/2", c)
-  os = ModelHamiltonians.heisenberg(c)
-  H = ttn(os, s)
   rng = StableRNG(1234)
+  g = NamedGraph(uniform_tree(10))
+  g = rename_vertices(v -> (v, 1), g)
+  s = siteinds("S=1/2", g)
+  os = ModelHamiltonians.heisenberg(g)
+  H = ttn(os, s)
   psi = random_ttn(rng, s; link_space=5)
   e, psi = dmrg(H, psi; nsweeps, maxdim, nsites)
 

From 4ca1ce1b6c71942d84e244d3a04fce6f24ee0db8 Mon Sep 17 00:00:00 2001
From: Joey <joeytindall94@googlemail.com>
Date: Thu, 17 Oct 2024 15:04:00 +0100
Subject: [PATCH 03/14] Remove file

---
 examples/test.jl | 21 ---------------------
 1 file changed, 21 deletions(-)
 delete mode 100644 examples/test.jl

diff --git a/examples/test.jl b/examples/test.jl
deleted file mode 100644
index 6bf1bfcf..00000000
--- a/examples/test.jl
+++ /dev/null
@@ -1,21 +0,0 @@
-using ITensorNetworks: IndsNetwork, siteinds, ttn
-using ITensorNetworks.ModelHamiltonians: ising
-using ITensors: Index, OpSum, terms, sites
-using NamedGraphs.NamedGraphGenerators: named_grid
-using NamedGraphs.GraphsExtensions: rem_vertex
-
-function filter_terms(H, verts)
-    H_new = OpSum()
-    for term in terms(H)
-        if isempty(filter(v -> v ∈ verts, sites(term)))
-            H_new += term
-        end
-    end
-    return H_new
-end
-
-g = named_grid((8,1))
-s = siteinds("S=1/2", g)
-H = ising(s)
-H_mod = filter_terms(H, [(4,1)])
-ttno = ttn(H_mod, s)
\ No newline at end of file

From 5296405f0588a9b304e1b8aad1ffbc7bc916d4dc Mon Sep 17 00:00:00 2001
From: Joey <joeytindall94@googlemail.com>
Date: Thu, 17 Oct 2024 15:05:50 +0100
Subject: [PATCH 04/14] Ortho fix and test

---
 Project.toml                                  | 10 ++---
 src/caches/beliefpropagationcache.jl          | 17 +++----
 .../alternating_update/alternating_update.jl  |  6 +--
 .../alternating_update/region_update.jl       | 45 ++++++++-----------
 src/solvers/defaults.jl                       | 18 +++++---
 test/Project.toml                             |  2 +-
 .../solvers.jl                                |  4 +-
 .../test_solvers/Project.toml                 |  3 +-
 .../test_solvers/test_dmrg.jl                 | 40 ++++++++++++++---
 .../test_solvers/test_tdvp_time_dependent.jl  |  4 +-
 10 files changed, 87 insertions(+), 62 deletions(-)

diff --git a/Project.toml b/Project.toml
index 64af22b6..fb5d1921 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "ITensorNetworks"
 uuid = "2919e153-833c-4bdc-8836-1ea460a35fc7"
 authors = ["Matthew Fishman <mfishman@flatironinstitute.org>, Joseph Tindall <jtindall@flatironinstitute.org> and contributors"]
-version = "0.11.15"
+version = "0.11.21"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
@@ -64,18 +64,18 @@ Graphs = "1.8"
 GraphsFlows = "0.1.1"
 ITensorMPS = "0.2.2"
 ITensors = "0.6.8"
-IsApprox = "0.1"
+IsApprox = "0.1, 1, 2"
 IterTools = "1.4.0"
-KrylovKit = "0.6, 0.7"
+KrylovKit = "0.6, 0.7, 0.8"
 MacroTools = "0.5"
 NDTensors = "0.3"
 NamedGraphs = "0.6.0"
-OMEinsumContractionOrders = "0.8.3"
+OMEinsumContractionOrders = "0.8.3, 0.9"
 Observers = "0.2.4"
 PackageExtensionCompat = "1"
 SerializedElementArrays = "0.1"
 SimpleTraits = "0.9"
-SparseArrayKit = "0.3"
+SparseArrayKit = "0.3, 0.4"
 SplitApplyCombine = "1.2"
 StaticArrays = "1.5.12"
 StructWalk = "0.2"
diff --git a/src/caches/beliefpropagationcache.jl b/src/caches/beliefpropagationcache.jl
index 3da4ca67..d5886ec7 100644
--- a/src/caches/beliefpropagationcache.jl
+++ b/src/caches/beliefpropagationcache.jl
@@ -16,11 +16,11 @@ using NDTensors: NDTensors
 
 default_message(elt, inds_e) = ITensor[denseblocks(delta(elt, i)) for i in inds_e]
 default_messages(ptn::PartitionedGraph) = Dictionary()
-function default_message_update(contract_list::Vector{ITensor}; kwargs...)
+function default_message_update(contract_list::Vector{ITensor}; normalize=true, kwargs...)
   sequence = optimal_contraction_sequence(contract_list)
   updated_messages = contract(contract_list; sequence, kwargs...)
   message_norm = norm(updated_messages)
-  if !iszero(message_norm)
+  if normalize && !iszero(message_norm)
     updated_messages /= message_norm
   end
   return ITensor[updated_messages]
@@ -106,7 +106,7 @@ end
 
 function message(bp_cache::BeliefPropagationCache, edge::PartitionEdge)
   mts = messages(bp_cache)
-  return get(mts, edge, default_message(bp_cache, edge))
+  return get(() -> default_message(bp_cache, edge), mts, edge)
 end
 function messages(bp_cache::BeliefPropagationCache, edges; kwargs...)
   return map(edge -> message(bp_cache, edge; kwargs...), edges)
@@ -152,15 +152,16 @@ end
 function environment(bp_cache::BeliefPropagationCache, verts::Vector)
   partition_verts = partitionvertices(bp_cache, verts)
   messages = environment(bp_cache, partition_verts)
-  central_tensors = ITensor[
-    tensornetwork(bp_cache)[v] for v in setdiff(vertices(bp_cache, partition_verts), verts)
-  ]
+  central_tensors = factors(bp_cache, setdiff(vertices(bp_cache, partition_verts), verts))
   return vcat(messages, central_tensors)
 end
 
+function factors(bp_cache::BeliefPropagationCache, verts::Vector)
+  return ITensor[tensornetwork(bp_cache)[v] for v in verts]
+end
+
 function factor(bp_cache::BeliefPropagationCache, vertex::PartitionVertex)
-  ptn = partitioned_tensornetwork(bp_cache)
-  return collect(eachtensor(subgraph(ptn, vertex)))
+  return factors(bp_cache, vertices(bp_cache, vertex))
 end
 
 """
diff --git a/src/solvers/alternating_update/alternating_update.jl b/src/solvers/alternating_update/alternating_update.jl
index 2cd5de71..750f3f36 100644
--- a/src/solvers/alternating_update/alternating_update.jl
+++ b/src/solvers/alternating_update/alternating_update.jl
@@ -9,8 +9,8 @@ function alternating_update(
   nsites, # define default for each level of solver implementation
   updater,  # this specifies the update performed locally
   outputlevel=default_outputlevel(),
-  region_printer=nothing,
-  sweep_printer=nothing,
+  region_printer=default_region_printer,
+  sweep_printer=default_sweep_printer,
   (sweep_observer!)=nothing,
   (region_observer!)=nothing,
   root_vertex=GraphsExtensions.default_root_vertex(init_state),
@@ -59,7 +59,7 @@ function alternating_update(
   (sweep_observer!)=nothing,
   sweep_printer=default_sweep_printer,#?
   (region_observer!)=nothing,
-  region_printer=nothing,
+  region_printer=default_region_printer,
 )
   state = copy(init_state)
   @assert !isnothing(sweep_plans)
diff --git a/src/solvers/alternating_update/region_update.jl b/src/solvers/alternating_update/region_update.jl
index b92adc8c..97241c20 100644
--- a/src/solvers/alternating_update/region_update.jl
+++ b/src/solvers/alternating_update/region_update.jl
@@ -7,36 +7,27 @@ function current_ortho(sweep_plan, which_region_update)
   if !isa(region, AbstractEdge) && length(region) == 1
     return only(current_verts)
   end
-  if which_region_update == length(regions)
-    # look back by one should be sufficient, but may be brittle?
-    overlapping_vertex = only(
-      intersect(current_verts, support(regions[which_region_update - 1]))
-    )
-    return overlapping_vertex
-  else
-    # look forward
-    other_regions = filter(
-      x -> !(issetequal(x, current_verts)), support.(regions[(which_region_update + 1):end])
+  # look forward
+  other_regions = filter(
+    x -> !(issetequal(x, current_verts)), support.(regions[(which_region_update + 1):end])
+  )
+  # find the first region that has overlapping support with current region 
+  ind = findfirst(x -> !isempty(intersect(support(x), support(region))), other_regions)
+  if isnothing(ind)
+    # look backward
+    other_regions = reverse(
+      filter(
+        x -> !(issetequal(x, current_verts)), support.(regions[1:(which_region_update - 1)])
+      ),
     )
-    # find the first region that has overlapping support with current region 
     ind = findfirst(x -> !isempty(intersect(support(x), support(region))), other_regions)
-    if isnothing(ind)
-      # look backward
-      other_regions = reverse(
-        filter(
-          x -> !(issetequal(x, current_verts)),
-          support.(regions[1:(which_region_update - 1)]),
-        ),
-      )
-      ind = findfirst(x -> !isempty(intersect(support(x), support(region))), other_regions)
-    end
-    @assert !isnothing(ind)
-    future_verts = union(support(other_regions[ind]))
-    # return ortho_ceter as the vertex in current region that does not overlap with following one
-    overlapping_vertex = intersect(current_verts, future_verts)
-    nonoverlapping_vertex = only(setdiff(current_verts, overlapping_vertex))
-    return nonoverlapping_vertex
   end
+  @assert !isnothing(ind)
+  future_verts = union(support(other_regions[ind]))
+  # return ortho_ceter as the vertex in current region that does not overlap with following one
+  overlapping_vertex = intersect(current_verts, future_verts)
+  nonoverlapping_vertex = only(setdiff(current_verts, overlapping_vertex))
+  return nonoverlapping_vertex
 end
 
 function region_update(
diff --git a/src/solvers/defaults.jl b/src/solvers/defaults.jl
index b5d315ff..09c2ae2f 100644
--- a/src/solvers/defaults.jl
+++ b/src/solvers/defaults.jl
@@ -1,4 +1,4 @@
-using Printf: @printf
+using Printf: @printf, @sprintf
 using ITensorMPS: maxlinkdim
 default_outputlevel() = 0
 default_nsites() = 2
@@ -7,10 +7,12 @@ default_extracter() = default_extracter
 default_inserter() = default_inserter
 default_checkdone() = (; kws...) -> false
 default_transform_operator() = nothing
+
+format(x) = @sprintf("%s", x)
+format(x::AbstractFloat) = @sprintf("%.1E", x)
+
 function default_region_printer(;
-  cutoff,
-  maxdim,
-  mindim,
+  inserter_kwargs,
   outputlevel,
   state,
   sweep_plan,
@@ -23,9 +25,11 @@ function default_region_printer(;
     region = first(sweep_plan[which_region_update])
     @printf("Sweep %d, region=%s \n", which_sweep, region)
     print("  Truncated using")
-    @printf(" cutoff=%.1E", cutoff)
-    @printf(" maxdim=%d", maxdim)
-    @printf(" mindim=%d", mindim)
+    for key in [:cutoff, :maxdim, :mindim]
+      if haskey(inserter_kwargs, key)
+        print(" ", key, "=", format(inserter_kwargs[key]))
+      end
+    end
     println()
     if spec != nothing
       @printf(
diff --git a/test/Project.toml b/test/Project.toml
index 70eb14d3..923e3bbe 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -21,7 +21,7 @@ NDTensors = "23ae76d9-e61a-49c4-8f12-3f1a16adf9cf"
 NamedGraphs = "678767b0-92e7-4007-89e4-4527a8725b19"
 OMEinsumContractionOrders = "6f22d1fd-8eed-4bb7-9776-e7d684900715"
 Observers = "338f10d5-c7f1-4033-a7d1-f9dec39bcaa0"
-OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SplitApplyCombine = "03a91e81-4c3e-53e1-a0a4-9c0c8f19dd66"
diff --git a/test/test_treetensornetworks/test_solvers/ITensorNetworksTestSolversUtils/solvers.jl b/test/test_treetensornetworks/test_solvers/ITensorNetworksTestSolversUtils/solvers.jl
index 82924f74..9b9568be 100644
--- a/test/test_treetensornetworks/test_solvers/ITensorNetworksTestSolversUtils/solvers.jl
+++ b/test/test_treetensornetworks/test_solvers/ITensorNetworksTestSolversUtils/solvers.jl
@@ -1,7 +1,7 @@
-using OrdinaryDiffEq: ODEProblem, Tsit5, solve
-using ITensors: ITensor
 using ITensorNetworks: TimeDependentSum, to_vec
+using ITensors: ITensor
 using KrylovKit: exponentiate
+using OrdinaryDiffEqTsit5: ODEProblem, Tsit5, solve
 
 function ode_solver(
   H::TimeDependentSum,
diff --git a/test/test_treetensornetworks/test_solvers/Project.toml b/test/test_treetensornetworks/test_solvers/Project.toml
index 77225041..e4716249 100644
--- a/test/test_treetensornetworks/test_solvers/Project.toml
+++ b/test/test_treetensornetworks/test_solvers/Project.toml
@@ -7,7 +7,8 @@ ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5"
 KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
 NamedGraphs = "678767b0-92e7-4007-89e4-4527a8725b19"
 Observers = "338f10d5-c7f1-4033-a7d1-f9dec39bcaa0"
-OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
+Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/test/test_treetensornetworks/test_solvers/test_dmrg.jl b/test/test_treetensornetworks/test_solvers/test_dmrg.jl
index b352d43c..004ec561 100644
--- a/test/test_treetensornetworks/test_solvers/test_dmrg.jl
+++ b/test/test_treetensornetworks/test_solvers/test_dmrg.jl
@@ -1,7 +1,7 @@
 @eval module $(gensym())
 using DataGraphs: edge_data, vertex_data
 using Dictionaries: Dictionary
-using Graphs: nv, vertices
+using Graphs: nv, vertices, uniform_tree
 using ITensorMPS: ITensorMPS
 using ITensorNetworks:
   ITensorNetworks,
@@ -19,9 +19,11 @@ using ITensorNetworks.ITensorsExtensions: replace_vertices
 using ITensorNetworks.ModelHamiltonians: ModelHamiltonians
 using ITensors: ITensors
 using KrylovKit: eigsolve
+using NamedGraphs: NamedGraph, rename_vertices
 using NamedGraphs.NamedGraphGenerators: named_comb_tree
 using Observers: observer
 using StableRNGs: StableRNG
+using Suppressor: @capture_out
 using Test: @test, @test_broken, @testset
 
 # This is needed since `eigen` is broken
@@ -76,6 +78,31 @@ ITensors.disable_auto_fermion()
   new_E = inner(psi', H, psi)
   @test new_E ≈ orig_E
   =#
+
+  #
+  # Test outputlevels are working
+  #
+  prev_output = ""
+  for outputlevel in 0:2
+    output = @capture_out begin
+      e, psi = dmrg(
+        H,
+        psi;
+        outputlevel,
+        nsweeps,
+        maxdim,
+        cutoff,
+        nsites,
+        updater_kwargs=(; krylovdim=3, maxiter=1),
+      )
+    end
+    if outputlevel == 0
+      @test length(output) == 0
+    else
+      @test length(output) > length(prev_output)
+    end
+    prev_output = output
+  end
 end
 
 @testset "Observers" begin
@@ -139,7 +166,7 @@ end
     nsweeps,
     maxdim,
     cutoff,
-    outputlevel=2,
+    outputlevel=0,
     transform_operator=ITensorNetworks.cache_operator_to_disk,
     transform_operator_kwargs=(; write_when_maxdim_exceeds=11),
   )
@@ -287,11 +314,12 @@ end
   nsites = 2
   nsweeps = 10
 
-  c = named_comb_tree((3, 2))
-  s = siteinds("S=1/2", c)
-  os = ModelHamiltonians.heisenberg(c)
-  H = ttn(os, s)
   rng = StableRNG(1234)
+  g = NamedGraph(uniform_tree(10))
+  g = rename_vertices(v -> (v, 1), g)
+  s = siteinds("S=1/2", g)
+  os = ModelHamiltonians.heisenberg(g)
+  H = ttn(os, s)
   psi = random_ttn(rng, s; link_space=5)
   e, psi = dmrg(H, psi; nsweeps, maxdim, nsites)
 
diff --git a/test/test_treetensornetworks/test_solvers/test_tdvp_time_dependent.jl b/test/test_treetensornetworks/test_solvers/test_tdvp_time_dependent.jl
index 17f1cc71..4101bc83 100644
--- a/test/test_treetensornetworks/test_solvers/test_tdvp_time_dependent.jl
+++ b/test/test_treetensornetworks/test_solvers/test_tdvp_time_dependent.jl
@@ -1,12 +1,12 @@
 @eval module $(gensym())
-using ITensors: contract
 using ITensorNetworks: ITensorNetworks, TimeDependentSum, ttn, mpo, mps, siteinds, tdvp
 using ITensorNetworks.ModelHamiltonians: ModelHamiltonians
-using OrdinaryDiffEq: Tsit5
+using ITensors: contract
 using KrylovKit: exponentiate
 using LinearAlgebra: norm
 using NamedGraphs: AbstractNamedEdge
 using NamedGraphs.NamedGraphGenerators: named_comb_tree
+using OrdinaryDiffEqTsit5: Tsit5
 using Test: @test, @test_broken, @testset
 
 include(

From d4608840a0a0600e9666ac18c99d9b9e49114b3c Mon Sep 17 00:00:00 2001
From: Joey <joeytindall94@googlemail.com>
Date: Fri, 18 Oct 2024 17:06:34 +0100
Subject: [PATCH 05/14] Improve orthogonalize. Simplify code

---
 src/abstractitensornetwork.jl                 | 49 ++++++++++++++-----
 .../alternating_update/region_update.jl       | 39 +--------------
 src/solvers/extract/extract.jl                | 13 ++---
 src/solvers/insert/insert.jl                  | 22 ++++-----
 .../abstracttreetensornetwork.jl              | 40 ++++++---------
 5 files changed, 71 insertions(+), 92 deletions(-)

diff --git a/src/abstractitensornetwork.jl b/src/abstractitensornetwork.jl
index 6f6ee164..653303aa 100644
--- a/src/abstractitensornetwork.jl
+++ b/src/abstractitensornetwork.jl
@@ -585,16 +585,25 @@ end
 
 # For ambiguity error; TODO: decide whether to use graph mutating methods when resulting graph is unchanged?
 function _orthogonalize_edge(tn::AbstractITensorNetwork, edge::AbstractEdge; kwargs...)
+  return _orthogonalize_edges(tn, [edge]; kwargs...)
+end
+
+# For ambiguity error; TODO: decide whether to use graph mutating methods when resulting graph is unchanged?
+function _orthogonalize_edges(
+  tn::AbstractITensorNetwork, edges::Vector{<:AbstractEdge}; kwargs...
+)
   # tn = factorize(tn, edge; kwargs...)
   # # TODO: Implement as `only(common_neighbors(tn, src(edge), dst(edge)))`
   # new_vertex = only(neighbors(tn, src(edge)) ∩ neighbors(tn, dst(edge)))
   # return contract(tn, new_vertex => dst(edge))
   tn = copy(tn)
-  left_inds = uniqueinds(tn, edge)
-  ltags = tags(tn, edge)
-  X, Y = factorize(tn[src(edge)], left_inds; tags=ltags, ortho="left", kwargs...)
-  tn[src(edge)] = X
-  tn[dst(edge)] *= Y
+  for edge in edges
+    left_inds = uniqueinds(tn, edge)
+    ltags = tags(tn, edge)
+    X, Y = factorize(tn[src(edge)], left_inds; tags=ltags, ortho="left", kwargs...)
+    tn[src(edge)] = X
+    tn[dst(edge)] *= Y
+  end
   return tn
 end
 
@@ -602,19 +611,35 @@ function ITensorMPS.orthogonalize(tn::AbstractITensorNetwork, edge::AbstractEdge
   return _orthogonalize_edge(tn, edge; kwargs...)
 end
 
+function ITensorMPS.orthogonalize(
+  tn::AbstractITensorNetwork, edges::Vector{<:AbstractEdge}; kwargs...
+)
+  return _orthogonalize_edges(tn, edges; kwargs...)
+end
+
 function ITensorMPS.orthogonalize(tn::AbstractITensorNetwork, edge::Pair; kwargs...)
   return orthogonalize(tn, edgetype(tn)(edge); kwargs...)
 end
 
-# Orthogonalize an ITensorNetwork towards a source vertex, treating
+function ITensorMPS.orthogonalize(
+  tn::AbstractITensorNetwork, edges::Vector{Pair}; kwargs...
+)
+  return orthogonalize(tn, edgetype(tn).(edges); kwargs...)
+end
+
+# Orthogonalize an ITensorNetwork towards a region, treating
 # the network as a tree spanned by a spanning tree.
 # TODO: Rename `tree_orthogonalize`.
-function ITensorMPS.orthogonalize(ψ::AbstractITensorNetwork, source_vertex)
-  spanning_tree_edges = post_order_dfs_edges(bfs_tree(ψ, source_vertex), source_vertex)
-  for e in spanning_tree_edges
-    ψ = orthogonalize(ψ, e)
-  end
-  return ψ
+function ITensorMPS.orthogonalize(ψ::AbstractITensorNetwork, region::Vector)
+  spanning_tree_edges = post_order_dfs_edges(bfs_tree(ψ, first(region)), first(region))
+  spanning_tree_edges = filter(
+    e -> !(src(e) ∈ region && dst(e) ∈ region), spanning_tree_edges
+  )
+  return orthogonalize(ψ, spanning_tree_edges)
+end
+
+function ITensorMPS.orthogonalize(ψ::AbstractITensorNetwork, region)
+  return orthogonalize(ψ, [region])
 end
 
 # TODO: decide whether to use graph mutating methods when resulting graph is unchanged?
diff --git a/src/solvers/alternating_update/region_update.jl b/src/solvers/alternating_update/region_update.jl
index 97241c20..6ad096b4 100644
--- a/src/solvers/alternating_update/region_update.jl
+++ b/src/solvers/alternating_update/region_update.jl
@@ -1,35 +1,3 @@
-#ToDo: generalize beyond 2-site
-#ToDo: remove concept of orthogonality center for generality
-function current_ortho(sweep_plan, which_region_update)
-  regions = first.(sweep_plan)
-  region = regions[which_region_update]
-  current_verts = support(region)
-  if !isa(region, AbstractEdge) && length(region) == 1
-    return only(current_verts)
-  end
-  # look forward
-  other_regions = filter(
-    x -> !(issetequal(x, current_verts)), support.(regions[(which_region_update + 1):end])
-  )
-  # find the first region that has overlapping support with current region 
-  ind = findfirst(x -> !isempty(intersect(support(x), support(region))), other_regions)
-  if isnothing(ind)
-    # look backward
-    other_regions = reverse(
-      filter(
-        x -> !(issetequal(x, current_verts)), support.(regions[1:(which_region_update - 1)])
-      ),
-    )
-    ind = findfirst(x -> !isempty(intersect(support(x), support(region))), other_regions)
-  end
-  @assert !isnothing(ind)
-  future_verts = union(support(other_regions[ind]))
-  # return ortho_ceter as the vertex in current region that does not overlap with following one
-  overlapping_vertex = intersect(current_verts, future_verts)
-  nonoverlapping_vertex = only(setdiff(current_verts, overlapping_vertex))
-  return nonoverlapping_vertex
-end
-
 function region_update(
   projected_operator,
   state;
@@ -55,14 +23,13 @@ function region_update(
 
   # ToDo: remove orthogonality center on vertex for generality
   # region carries same information 
-  ortho_vertex = current_ortho(sweep_plan, which_region_update)
   if !isnothing(transform_operator)
     projected_operator = transform_operator(
       state, projected_operator; outputlevel, transform_operator_kwargs...
     )
   end
   state, projected_operator, phi = extracter(
-    state, projected_operator, region, ortho_vertex; extracter_kwargs..., internal_kwargs
+    state, projected_operator, region; extracter_kwargs..., internal_kwargs
   )
   # create references, in case solver does (out-of-place) modify PH or state
   state! = Ref(state)
@@ -88,9 +55,7 @@ function region_update(
   #  drho = noise * noiseterm(PH, phi, ortho) # TODO: actually implement this for trees...
   # so noiseterm is a solver
   #end
-  state, spec = inserter(
-    state, phi, region, ortho_vertex; inserter_kwargs..., internal_kwargs
-  )
+  state, spec = inserter(state, phi, region; inserter_kwargs..., internal_kwargs)
   all_kwargs = (;
     which_region_update,
     sweep_plan,
diff --git a/src/solvers/extract/extract.jl b/src/solvers/extract/extract.jl
index feb57c2f..431394ea 100644
--- a/src/solvers/extract/extract.jl
+++ b/src/solvers/extract/extract.jl
@@ -7,18 +7,19 @@
 # insert_local_tensors takes that tensor and factorizes it back
 # apart and puts it back into the network.
 #
-function default_extracter(state, projected_operator, region, ortho; internal_kwargs)
-  state = orthogonalize(state, ortho)
+function default_extracter(state, projected_operator, region; internal_kwargs)
   if isa(region, AbstractEdge)
-    other_vertex = only(setdiff(support(region), [ortho]))
-    left_inds = uniqueinds(state[ortho], state[other_vertex])
+    vsrc, vdst = src(region), dst(region)
+    state = orthogonalize(state, vsrc)
+    left_inds = uniqueinds(state[vsrc], state[vdst])
     #ToDo: replace with call to factorize
     U, S, V = svd(
-      state[ortho], left_inds; lefttags=tags(state, region), righttags=tags(state, region)
+      state[vsrc], left_inds; lefttags=tags(state, region), righttags=tags(state, region)
     )
-    state[ortho] = U
+    state[vsrc] = U
     local_tensor = S * V
   else
+    state = orthogonalize(state, region)
     local_tensor = prod(state[v] for v in region)
   end
   projected_operator = position(projected_operator, state, region)
diff --git a/src/solvers/insert/insert.jl b/src/solvers/insert/insert.jl
index 11aed223..42b32023 100644
--- a/src/solvers/insert/insert.jl
+++ b/src/solvers/insert/insert.jl
@@ -6,8 +6,7 @@
 function default_inserter(
   state::AbstractTTN,
   phi::ITensor,
-  region,
-  ortho_vert;
+  region;
   normalize=false,
   maxdim=nothing,
   mindim=nothing,
@@ -16,16 +15,14 @@ function default_inserter(
 )
   state = copy(state)
   spec = nothing
-  other_vertex = setdiff(support(region), [ortho_vert])
-  if !isempty(other_vertex)
-    v = only(other_vertex)
-    e = edgetype(state)(ortho_vert, v)
-    indsTe = inds(state[ortho_vert])
+  if length(region) == 2
+    v = last(region)
+    e = edgetype(state)(first(region), last(region))
+    indsTe = inds(state[first(region)])
     L, phi, spec = factorize(phi, indsTe; tags=tags(state, e), maxdim, mindim, cutoff)
-    state[ortho_vert] = L
-
+    state[first(region)] = L
   else
-    v = ortho_vert
+    v = only(region)
   end
   state[v] = phi
   state = set_ortho_region(state, [v])
@@ -44,8 +41,7 @@ function default_inserter(
   normalize=false,
   internal_kwargs,
 )
-  v = only(setdiff(support(region), [ortho]))
-  state[v] *= phi
-  state = set_ortho_region(state, [v])
+  state[dst(region)] *= phi
+  state = set_ortho_region(state, [dst(region)])
   return state, nothing
 end
diff --git a/src/treetensornetworks/abstracttreetensornetwork.jl b/src/treetensornetworks/abstracttreetensornetwork.jl
index c8dccb1f..13d78197 100644
--- a/src/treetensornetworks/abstracttreetensornetwork.jl
+++ b/src/treetensornetworks/abstracttreetensornetwork.jl
@@ -1,6 +1,12 @@
 using Graphs: has_vertex
 using NamedGraphs.GraphsExtensions:
-  GraphsExtensions, edge_path, leaf_vertices, post_order_dfs_edges, post_order_dfs_vertices
+  GraphsExtensions,
+  edge_path,
+  leaf_vertices,
+  post_order_dfs_edges,
+  post_order_dfs_vertices,
+  a_star
+using NamedGraphs: namedgraph_a_star
 using IsApprox: IsApprox, Approx
 using ITensors: @Algorithm_str, directsum, hasinds, permute, plev
 using ITensorMPS: linkind, loginner, lognorm, orthogonalize
@@ -29,30 +35,16 @@ function set_ortho_region(tn::AbstractTTN, new_region)
   return error("Not implemented")
 end
 
-# 
-# Orthogonalization
-# 
-
-function ITensorMPS.orthogonalize(tn::AbstractTTN, ortho_center; kwargs...)
-  if isone(length(ortho_region(tn))) && ortho_center == only(ortho_region(tn))
-    return tn
-  end
-  # TODO: Rewrite this in a more general way.
-  if isone(length(ortho_region(tn)))
-    edge_list = edge_path(tn, only(ortho_region(tn)), ortho_center)
-  else
-    edge_list = post_order_dfs_edges(tn, ortho_center)
-  end
-  for e in edge_list
-    tn = orthogonalize(tn, e)
+function ITensorMPS.orthogonalize(ttn::AbstractTTN, region::Vector; kwargs...)
+  paths = [
+    namedgraph_a_star(underlying_graph(ttn), rp, r) for r in region for
+    rp in ortho_region(ttn)
+  ]
+  path = unique(reduce(vcat, paths))
+  if !isempty(path)
+    ttn = typeof(ttn)(orthogonalize(ITensorNetwork(ttn), path; kwargs...))
   end
-  return set_ortho_region(tn, typeof(ortho_region(tn))([ortho_center]))
-end
-
-# For ambiguity error
-
-function ITensorMPS.orthogonalize(tn::AbstractTTN, edge::AbstractEdge; kwargs...)
-  return typeof(tn)(orthogonalize(ITensorNetwork(tn), edge; kwargs...))
+  return set_ortho_region(ttn, region)
 end
 
 # 

From 32ebca76e9029b15249f65625204a4f5548e957a Mon Sep 17 00:00:00 2001
From: Joey <joeytindall94@googlemail.com>
Date: Fri, 18 Oct 2024 17:53:43 +0100
Subject: [PATCH 06/14] Post order dsf edges region

---
 src/abstractitensornetwork.jl                       | 5 +----
 src/edge_sequences.jl                               | 5 +++++
 src/treetensornetworks/abstracttreetensornetwork.jl | 9 +++------
 3 files changed, 9 insertions(+), 10 deletions(-)

diff --git a/src/abstractitensornetwork.jl b/src/abstractitensornetwork.jl
index 653303aa..244a24d7 100644
--- a/src/abstractitensornetwork.jl
+++ b/src/abstractitensornetwork.jl
@@ -631,10 +631,7 @@ end
 # the network as a tree spanned by a spanning tree.
 # TODO: Rename `tree_orthogonalize`.
 function ITensorMPS.orthogonalize(ψ::AbstractITensorNetwork, region::Vector)
-  spanning_tree_edges = post_order_dfs_edges(bfs_tree(ψ, first(region)), first(region))
-  spanning_tree_edges = filter(
-    e -> !(src(e) ∈ region && dst(e) ∈ region), spanning_tree_edges
-  )
+  spanning_tree_edges = post_order_dfs_edges_region(bfs_tree(ψ, first(region)), region)
   return orthogonalize(ψ, spanning_tree_edges)
 end
 
diff --git a/src/edge_sequences.jl b/src/edge_sequences.jl
index 9dab9fff..9b385c68 100644
--- a/src/edge_sequences.jl
+++ b/src/edge_sequences.jl
@@ -44,3 +44,8 @@ end
 @traitfn function edge_sequence(::Algorithm"parallel", g::::(!IsDirected))
   return [[e] for e in vcat(edges(g), reverse.(edges(g)))]
 end
+
+function post_order_dfs_edges_region(g::AbstractGraph, region)
+  es = post_order_dfs_edges(g, first(region))
+  return filter(e -> !((src(e) ∈ region) && (dst(e) ∈ region)), es)
+end
diff --git a/src/treetensornetworks/abstracttreetensornetwork.jl b/src/treetensornetworks/abstracttreetensornetwork.jl
index 13d78197..f75e9de2 100644
--- a/src/treetensornetworks/abstracttreetensornetwork.jl
+++ b/src/treetensornetworks/abstracttreetensornetwork.jl
@@ -10,7 +10,6 @@ using NamedGraphs: namedgraph_a_star
 using IsApprox: IsApprox, Approx
 using ITensors: @Algorithm_str, directsum, hasinds, permute, plev
 using ITensorMPS: linkind, loginner, lognorm, orthogonalize
-using TupleTools: TupleTools
 
 abstract type AbstractTreeTensorNetwork{V} <: AbstractITensorNetwork{V} end
 
@@ -36,11 +35,9 @@ function set_ortho_region(tn::AbstractTTN, new_region)
 end
 
 function ITensorMPS.orthogonalize(ttn::AbstractTTN, region::Vector; kwargs...)
-  paths = [
-    namedgraph_a_star(underlying_graph(ttn), rp, r) for r in region for
-    rp in ortho_region(ttn)
-  ]
-  path = unique(reduce(vcat, paths))
+  new_path = post_order_dfs_edges_region(ttn, region)
+  existing_path = post_order_dfs_edges_region(ttn, ortho_region(ttn))
+  path = setdiff(new_path, existing_path)
   if !isempty(path)
     ttn = typeof(ttn)(orthogonalize(ITensorNetwork(ttn), path; kwargs...))
   end

From 8a0b34e3aba7f85fae849b9252f29122f9c54183 Mon Sep 17 00:00:00 2001
From: Joey <joeytindall94@googlemail.com>
Date: Sun, 20 Oct 2024 17:52:00 +0100
Subject: [PATCH 07/14] Reverse edges in the reverse sweep of TDVP

---
 src/solvers/alternating_update/region_update.jl |  1 +
 src/solvers/insert/insert.jl                    |  3 +--
 src/solvers/sweep_plans/sweep_plans.jl          | 13 ++++++++-----
 3 files changed, 10 insertions(+), 7 deletions(-)

diff --git a/src/solvers/alternating_update/region_update.jl b/src/solvers/alternating_update/region_update.jl
index 6ad096b4..c741c82a 100644
--- a/src/solvers/alternating_update/region_update.jl
+++ b/src/solvers/alternating_update/region_update.jl
@@ -55,6 +55,7 @@ function region_update(
   #  drho = noise * noiseterm(PH, phi, ortho) # TODO: actually implement this for trees...
   # so noiseterm is a solver
   #end
+  #if isa(region, AbstractEdge) && 
   state, spec = inserter(state, phi, region; inserter_kwargs..., internal_kwargs)
   all_kwargs = (;
     which_region_update,
diff --git a/src/solvers/insert/insert.jl b/src/solvers/insert/insert.jl
index 42b32023..01fb35bd 100644
--- a/src/solvers/insert/insert.jl
+++ b/src/solvers/insert/insert.jl
@@ -33,8 +33,7 @@ end
 function default_inserter(
   state::AbstractTTN,
   phi::ITensor,
-  region::NamedEdge,
-  ortho;
+  region::NamedEdge;
   cutoff=nothing,
   maxdim=nothing,
   mindim=nothing,
diff --git a/src/solvers/sweep_plans/sweep_plans.jl b/src/solvers/sweep_plans/sweep_plans.jl
index 69221995..384cf8dc 100644
--- a/src/solvers/sweep_plans/sweep_plans.jl
+++ b/src/solvers/sweep_plans/sweep_plans.jl
@@ -13,10 +13,11 @@ end
 
 support(r) = r
 
-function reverse_region(edges, which_edge; nsites=1, region_kwargs=(;))
+function reverse_region(edges, which_edge; reverse_edge=false, nsites=1, region_kwargs=(;))
   current_edge = edges[which_edge]
   if nsites == 1
-    return [(current_edge, region_kwargs)]
+    !reverse_edge && return [(current_edge, region_kwargs)]
+    reverse_edge && return [(reverse(current_edge), region_kwargs)]
   elseif nsites == 2
     if last(edges) == current_edge
       return ()
@@ -62,6 +63,7 @@ function forward_sweep(
   dir::Base.ForwardOrdering,
   graph::AbstractGraph;
   root_vertex=GraphsExtensions.default_root_vertex(graph),
+  reverse_edges=false,
   region_kwargs,
   reverse_kwargs=region_kwargs,
   reverse_step=false,
@@ -71,12 +73,13 @@ function forward_sweep(
   regions = collect(
     flatten(map(i -> forward_region(edges, i; region_kwargs, kwargs...), eachindex(edges)))
   )
-
   if reverse_step
     reverse_regions = collect(
       flatten(
         map(
-          i -> reverse_region(edges, i; region_kwargs=reverse_kwargs, kwargs...),
+          i -> reverse_region(
+            edges, i; reverse_edge=reverse_edges, region_kwargs=reverse_kwargs, kwargs...
+          ),
           eachindex(edges),
         ),
       ),
@@ -90,7 +93,7 @@ end
 
 #ToDo: is there a better name for this? unidirectional_sweep? traversal?
 function forward_sweep(dir::Base.ReverseOrdering, args...; kwargs...)
-  return reverse(forward_sweep(Base.Forward, args...; kwargs...))
+  return reverse(forward_sweep(Base.Forward, args...; reverse_edges=true, kwargs...))
 end
 
 function default_sweep_plans(

From be0e7dcdf4a9924bf4def8ebfd96948c6d41bb5e Mon Sep 17 00:00:00 2001
From: Joey <joeytindall94@googlemail.com>
Date: Sun, 3 Nov 2024 15:18:43 +0000
Subject: [PATCH 08/14] Bug fix and renaming

---
 src/abstractitensornetwork.jl                 | 30 ++++++-------------
 .../abstracttreetensornetwork.jl              |  2 +-
 2 files changed, 10 insertions(+), 22 deletions(-)

diff --git a/src/abstractitensornetwork.jl b/src/abstractitensornetwork.jl
index 244a24d7..57133bd4 100644
--- a/src/abstractitensornetwork.jl
+++ b/src/abstractitensornetwork.jl
@@ -584,12 +584,16 @@ function LinearAlgebra.factorize(tn::AbstractITensorNetwork, edge::Pair; kwargs.
 end
 
 # For ambiguity error; TODO: decide whether to use graph mutating methods when resulting graph is unchanged?
-function _orthogonalize_edge(tn::AbstractITensorNetwork, edge::AbstractEdge; kwargs...)
-  return _orthogonalize_edges(tn, [edge]; kwargs...)
+function orthogonalize_edge(tn::AbstractITensorNetwork, edge::AbstractEdge; kwargs...)
+  return orthogonalize_path(tn, [edge]; kwargs...)
+end
+
+function orthogonalize_edge(tn::AbstractITensorNetwork, edge::Pair; kwargs...)
+  return orthogonalize_edge(tn, edgetype(tn)(edge); kwargs...)
 end
 
 # For ambiguity error; TODO: decide whether to use graph mutating methods when resulting graph is unchanged?
-function _orthogonalize_edges(
+function orthogonalize_path(
   tn::AbstractITensorNetwork, edges::Vector{<:AbstractEdge}; kwargs...
 )
   # tn = factorize(tn, edge; kwargs...)
@@ -607,24 +611,8 @@ function _orthogonalize_edges(
   return tn
 end
 
-function ITensorMPS.orthogonalize(tn::AbstractITensorNetwork, edge::AbstractEdge; kwargs...)
-  return _orthogonalize_edge(tn, edge; kwargs...)
-end
-
-function ITensorMPS.orthogonalize(
-  tn::AbstractITensorNetwork, edges::Vector{<:AbstractEdge}; kwargs...
-)
-  return _orthogonalize_edges(tn, edges; kwargs...)
-end
-
-function ITensorMPS.orthogonalize(tn::AbstractITensorNetwork, edge::Pair; kwargs...)
-  return orthogonalize(tn, edgetype(tn)(edge); kwargs...)
-end
-
-function ITensorMPS.orthogonalize(
-  tn::AbstractITensorNetwork, edges::Vector{Pair}; kwargs...
-)
-  return orthogonalize(tn, edgetype(tn).(edges); kwargs...)
+function orthogonalize_path(tn::AbstractITensorNetwork, edges::Vector{<:Pair}; kwargs...)
+  return orthogonalize_path(tn, edgetype(tn).(edges); kwargs...)
 end
 
 # Orthogonalize an ITensorNetwork towards a region, treating
diff --git a/src/treetensornetworks/abstracttreetensornetwork.jl b/src/treetensornetworks/abstracttreetensornetwork.jl
index f75e9de2..d0029334 100644
--- a/src/treetensornetworks/abstracttreetensornetwork.jl
+++ b/src/treetensornetworks/abstracttreetensornetwork.jl
@@ -39,7 +39,7 @@ function ITensorMPS.orthogonalize(ttn::AbstractTTN, region::Vector; kwargs...)
   existing_path = post_order_dfs_edges_region(ttn, ortho_region(ttn))
   path = setdiff(new_path, existing_path)
   if !isempty(path)
-    ttn = typeof(ttn)(orthogonalize(ITensorNetwork(ttn), path; kwargs...))
+    ttn = typeof(ttn)(orthogonalize_path(ITensorNetwork(ttn), path; kwargs...))
   end
   return set_ortho_region(ttn, region)
 end

From f0e8bf55b99d0785ace54a20ba9d8d4238c9396f Mon Sep 17 00:00:00 2001
From: Joey <joeytindall94@googlemail.com>
Date: Sun, 3 Nov 2024 18:46:19 -0500
Subject: [PATCH 09/14] Further renaming

---
 src/abstractitensornetwork.jl | 24 ++++++++++++++++++++----
 1 file changed, 20 insertions(+), 4 deletions(-)

diff --git a/src/abstractitensornetwork.jl b/src/abstractitensornetwork.jl
index ba35c0d2..09d1eab9 100644
--- a/src/abstractitensornetwork.jl
+++ b/src/abstractitensornetwork.jl
@@ -584,12 +584,12 @@ function LinearAlgebra.factorize(tn::AbstractITensorNetwork, edge::Pair; kwargs.
 end
 
 # For ambiguity error; TODO: decide whether to use graph mutating methods when resulting graph is unchanged?
-function orthogonalize_edge(tn::AbstractITensorNetwork, edge::AbstractEdge; kwargs...)
+function orthogonalize_path(tn::AbstractITensorNetwork, edge::AbstractEdge; kwargs...)
   return orthogonalize_path(tn, [edge]; kwargs...)
 end
 
-function orthogonalize_edge(tn::AbstractITensorNetwork, edge::Pair; kwargs...)
-  return orthogonalize_edge(tn, edgetype(tn)(edge); kwargs...)
+function orthogonalize_path(tn::AbstractITensorNetwork, edge::Pair; kwargs...)
+  return orthogonalize_path(tn, edgetype(tn)(edge); kwargs...)
 end
 
 # For ambiguity error; TODO: decide whether to use graph mutating methods when resulting graph is unchanged?
@@ -620,13 +620,29 @@ end
 # TODO: Rename `tree_orthogonalize`.
 function ITensorMPS.orthogonalize(ψ::AbstractITensorNetwork, region::Vector)
   spanning_tree_edges = post_order_dfs_edges_region(bfs_tree(ψ, first(region)), region)
-  return orthogonalize(ψ, spanning_tree_edges)
+  return orthogonalize_path(ψ, spanning_tree_edges)
 end
 
 function ITensorMPS.orthogonalize(ψ::AbstractITensorNetwork, region)
   return orthogonalize(ψ, [region])
 end
 
+function ITensorMPS.orthogonalize(ψ::AbstractITensorNetwork, edges::Vector{<:AbstractEdge})
+  return orthogonalize(ψ, unique(vcat([src(e) for e in edges], [dst(e) for e in edges])))
+end
+
+function ITensorMPS.orthogonalize(ψ::AbstractITensorNetwork, edges::Vector{<:Pair})
+  return orthogonalize(ψ, edgetype(ψ).(edges))
+end
+
+function ITensorMPS.orthogonalize(ψ::AbstractITensorNetwork, edge::AbstractEdge)
+  return orthogonalize(ψ, [edge])
+end
+
+function ITensorMPS.orthogonalize(ψ::AbstractITensorNetwork, edge::Pair)
+  return orthogonalize(ψ, edgetype(ψ)(edge))
+end
+
 # TODO: decide whether to use graph mutating methods when resulting graph is unchanged?
 function _truncate_edge(tn::AbstractITensorNetwork, edge::AbstractEdge; kwargs...)
   tn = copy(tn)

From f5d3aa42b7211d453beb7d18e47e7bea22d1b186 Mon Sep 17 00:00:00 2001
From: Joey <joeytindall94@googlemail.com>
Date: Thu, 7 Nov 2024 16:50:40 -0500
Subject: [PATCH 10/14] Clean up orthogonalize code

---
 src/abstractitensornetwork.jl                 | 31 +++++--------------
 src/edge_sequences.jl                         |  5 ---
 .../abstracttreetensornetwork.jl              | 12 ++++---
 3 files changed, 16 insertions(+), 32 deletions(-)

diff --git a/src/abstractitensornetwork.jl b/src/abstractitensornetwork.jl
index 09d1eab9..7f771c8e 100644
--- a/src/abstractitensornetwork.jl
+++ b/src/abstractitensornetwork.jl
@@ -40,7 +40,7 @@ using ITensorMPS: ITensorMPS, add, linkdim, linkinds, siteinds
 using .ITensorsExtensions: ITensorsExtensions, indtype, promote_indtype
 using LinearAlgebra: LinearAlgebra, factorize
 using MacroTools: @capture
-using NamedGraphs: NamedGraphs, NamedGraph, not_implemented
+using NamedGraphs: NamedGraphs, NamedGraph, not_implemented, steiner_tree
 using NamedGraphs.GraphsExtensions:
   ⊔, directed_graph, incident_edges, rename_vertices, vertextype
 using NDTensors: NDTensors, dim
@@ -617,30 +617,15 @@ end
 
 # Orthogonalize an ITensorNetwork towards a region, treating
 # the network as a tree spanned by a spanning tree.
-# TODO: Rename `tree_orthogonalize`.
-function ITensorMPS.orthogonalize(ψ::AbstractITensorNetwork, region::Vector)
-  spanning_tree_edges = post_order_dfs_edges_region(bfs_tree(ψ, first(region)), region)
-  return orthogonalize_path(ψ, spanning_tree_edges)
+function tree_orthogonalize(ψ::AbstractITensorNetwork, region::Vector)
+  region = collect(vertices(steiner_tree(underlying_graph(ψ), region)))
+  path = post_order_dfs_edges(bfs_tree(ψ, first(region)), first(region))
+  path = filter(e -> !((src(e) ∈ region) && (dst(e) ∈ region)), path)
+  return orthogonalize_path(ψ, path)
 end
 
-function ITensorMPS.orthogonalize(ψ::AbstractITensorNetwork, region)
-  return orthogonalize(ψ, [region])
-end
-
-function ITensorMPS.orthogonalize(ψ::AbstractITensorNetwork, edges::Vector{<:AbstractEdge})
-  return orthogonalize(ψ, unique(vcat([src(e) for e in edges], [dst(e) for e in edges])))
-end
-
-function ITensorMPS.orthogonalize(ψ::AbstractITensorNetwork, edges::Vector{<:Pair})
-  return orthogonalize(ψ, edgetype(ψ).(edges))
-end
-
-function ITensorMPS.orthogonalize(ψ::AbstractITensorNetwork, edge::AbstractEdge)
-  return orthogonalize(ψ, [edge])
-end
-
-function ITensorMPS.orthogonalize(ψ::AbstractITensorNetwork, edge::Pair)
-  return orthogonalize(ψ, edgetype(ψ)(edge))
+function tree_orthogonalize(ψ::AbstractITensorNetwork, region)
+  return tree_orthogonalize(ψ, [region])
 end
 
 # TODO: decide whether to use graph mutating methods when resulting graph is unchanged?
diff --git a/src/edge_sequences.jl b/src/edge_sequences.jl
index 9b385c68..9dab9fff 100644
--- a/src/edge_sequences.jl
+++ b/src/edge_sequences.jl
@@ -44,8 +44,3 @@ end
 @traitfn function edge_sequence(::Algorithm"parallel", g::::(!IsDirected))
   return [[e] for e in vcat(edges(g), reverse.(edges(g)))]
 end
-
-function post_order_dfs_edges_region(g::AbstractGraph, region)
-  es = post_order_dfs_edges(g, first(region))
-  return filter(e -> !((src(e) ∈ region) && (dst(e) ∈ region)), es)
-end
diff --git a/src/treetensornetworks/abstracttreetensornetwork.jl b/src/treetensornetworks/abstracttreetensornetwork.jl
index 0c424cd6..5f079299 100644
--- a/src/treetensornetworks/abstracttreetensornetwork.jl
+++ b/src/treetensornetworks/abstracttreetensornetwork.jl
@@ -6,7 +6,7 @@ using NamedGraphs.GraphsExtensions:
   post_order_dfs_edges,
   post_order_dfs_vertices,
   a_star
-using NamedGraphs: namedgraph_a_star
+using NamedGraphs: namedgraph_a_star, steiner_tree
 using IsApprox: IsApprox, Approx
 using ITensors: ITensors, @Algorithm_str, directsum, hasinds, permute, plev
 using ITensorMPS: ITensorMPS, linkind, loginner, lognorm, orthogonalize
@@ -36,15 +36,19 @@ function set_ortho_region(tn::AbstractTTN, new_region)
 end
 
 function ITensorMPS.orthogonalize(ttn::AbstractTTN, region::Vector; kwargs...)
-  new_path = post_order_dfs_edges_region(ttn, region)
-  existing_path = post_order_dfs_edges_region(ttn, ortho_region(ttn))
-  path = setdiff(new_path, existing_path)
+  st = steiner_tree(ttn, union(region, ortho_region(ttn)))
+  path = post_order_dfs_edges(st, first(region))
+  path = filter(e -> !((src(e) ∈ region) && (dst(e) ∈ region)), path)
   if !isempty(path)
     ttn = typeof(ttn)(orthogonalize_path(ITensorNetwork(ttn), path; kwargs...))
   end
   return set_ortho_region(ttn, region)
 end
 
+function ITensorMPS.orthogonalize(ttn::AbstractTTN, region; kwargs...)
+  return orthogonalize(ttn, [region]; kwargs...)
+end
+
 # 
 # Truncation
 # 

From fcf5b989a7089e8007805ddc4e83713be9397e3f Mon Sep 17 00:00:00 2001
From: Joey <joeytindall94@googlemail.com>
Date: Fri, 8 Nov 2024 08:48:51 -0500
Subject: [PATCH 11/14] Improve orthogonalize method efficiency

---
 src/abstractitensornetwork.jl                       | 6 ++++--
 src/solvers/extract/extract.jl                      | 3 ++-
 src/treetensornetworks/abstracttreetensornetwork.jl | 5 +++++
 test/test_itensornetwork.jl                         | 5 +++--
 4 files changed, 14 insertions(+), 5 deletions(-)

diff --git a/src/abstractitensornetwork.jl b/src/abstractitensornetwork.jl
index 7f771c8e..f7c148d0 100644
--- a/src/abstractitensornetwork.jl
+++ b/src/abstractitensornetwork.jl
@@ -7,6 +7,7 @@ using Graphs:
   add_edge!,
   add_vertex!,
   bfs_tree,
+  center,
   dst,
   edges,
   edgetype,
@@ -618,8 +619,9 @@ end
 # Orthogonalize an ITensorNetwork towards a region, treating
 # the network as a tree spanned by a spanning tree.
 function tree_orthogonalize(ψ::AbstractITensorNetwork, region::Vector)
-  region = collect(vertices(steiner_tree(underlying_graph(ψ), region)))
-  path = post_order_dfs_edges(bfs_tree(ψ, first(region)), first(region))
+  region_centre =
+    length(region) != 1 ? first(center(steiner_tree(ψ, region))) : only(region)
+  path = post_order_dfs_edges(bfs_tree(ψ, region_centre), region_centre)
   path = filter(e -> !((src(e) ∈ region) && (dst(e) ∈ region)), path)
   return orthogonalize_path(ψ, path)
 end
diff --git a/src/solvers/extract/extract.jl b/src/solvers/extract/extract.jl
index 431394ea..1013d1bd 100644
--- a/src/solvers/extract/extract.jl
+++ b/src/solvers/extract/extract.jl
@@ -7,12 +7,13 @@
 # insert_local_tensors takes that tensor and factorizes it back
 # apart and puts it back into the network.
 #
+
 function default_extracter(state, projected_operator, region; internal_kwargs)
   if isa(region, AbstractEdge)
+    # TODO: add functionality for orthogonalizing onto a bond so that can be called instead
     vsrc, vdst = src(region), dst(region)
     state = orthogonalize(state, vsrc)
     left_inds = uniqueinds(state[vsrc], state[vdst])
-    #ToDo: replace with call to factorize
     U, S, V = svd(
       state[vsrc], left_inds; lefttags=tags(state, region), righttags=tags(state, region)
     )
diff --git a/src/treetensornetworks/abstracttreetensornetwork.jl b/src/treetensornetworks/abstracttreetensornetwork.jl
index 5f079299..d3f608b0 100644
--- a/src/treetensornetworks/abstracttreetensornetwork.jl
+++ b/src/treetensornetworks/abstracttreetensornetwork.jl
@@ -36,6 +36,11 @@ function set_ortho_region(tn::AbstractTTN, new_region)
 end
 
 function ITensorMPS.orthogonalize(ttn::AbstractTTN, region::Vector; kwargs...)
+  return orthogonalize_ttn(ttn, region; kwargs...)
+end
+
+function orthogonalize_ttn(ttn::AbstractTTN, region::Vector; kwargs...)
+  issetequal(region, ortho_region(ttn)) && return ttn
   st = steiner_tree(ttn, union(region, ortho_region(ttn)))
   path = post_order_dfs_edges(st, first(region))
   path = filter(e -> !((src(e) ∈ region) && (dst(e) ∈ region)), path)
diff --git a/test/test_itensornetwork.jl b/test/test_itensornetwork.jl
index ba3caa01..53e2928f 100644
--- a/test/test_itensornetwork.jl
+++ b/test/test_itensornetwork.jl
@@ -51,6 +51,7 @@ using ITensorNetworks:
   orthogonalize,
   random_tensornetwork,
   siteinds,
+  tree_orthogonalize,
   ttn
 using LinearAlgebra: factorize
 using NamedGraphs: NamedEdge
@@ -287,13 +288,13 @@ const elts = (Float32, Float64, Complex{Float32}, Complex{Float64})
     @test nv(tn_ortho) == 5
     @test nv(tn) == 4
     @test Z ≈ Z̃
-    tn_ortho = orthogonalize(tn, 4 => 3)
+    tn_ortho = tree_orthogonalize(tn, [3, 4])
     Z̃ = norm_sqr(tn_ortho)
     @test nv(tn_ortho) == 4
     @test nv(tn) == 4
     @test Z ≈ Z̃
 
-    tn_ortho = orthogonalize(tn, 1)
+    tn_ortho = tree_orthogonalize(tn, 1)
     Z̃ = norm_sqr(tn_ortho)
     @test Z ≈ Z̃
     Z̃ = inner(tn_ortho, tn)

From f5098e121c16552012dc3d505340e07af88e2600 Mon Sep 17 00:00:00 2001
From: Joey <joeytindall94@googlemail.com>
Date: Fri, 8 Nov 2024 08:51:01 -0500
Subject: [PATCH 12/14] Fix function name

---
 src/treetensornetworks/abstracttreetensornetwork.jl | 4 ----
 1 file changed, 4 deletions(-)

diff --git a/src/treetensornetworks/abstracttreetensornetwork.jl b/src/treetensornetworks/abstracttreetensornetwork.jl
index d3f608b0..a08fe27d 100644
--- a/src/treetensornetworks/abstracttreetensornetwork.jl
+++ b/src/treetensornetworks/abstracttreetensornetwork.jl
@@ -36,10 +36,6 @@ function set_ortho_region(tn::AbstractTTN, new_region)
 end
 
 function ITensorMPS.orthogonalize(ttn::AbstractTTN, region::Vector; kwargs...)
-  return orthogonalize_ttn(ttn, region; kwargs...)
-end
-
-function orthogonalize_ttn(ttn::AbstractTTN, region::Vector; kwargs...)
   issetequal(region, ortho_region(ttn)) && return ttn
   st = steiner_tree(ttn, union(region, ortho_region(ttn)))
   path = post_order_dfs_edges(st, first(region))

From 629af711e5931b36d85d965777b372fa141c60b9 Mon Sep 17 00:00:00 2001
From: Joey <joeytindall94@googlemail.com>
Date: Fri, 8 Nov 2024 10:02:21 -0500
Subject: [PATCH 13/14] Wrap tree orthogonalize for ttns

---
 src/apply.jl                                        | 4 ++--
 src/tebd.jl                                         | 2 +-
 src/treetensornetworks/abstracttreetensornetwork.jl | 4 ++++
 3 files changed, 7 insertions(+), 3 deletions(-)

diff --git a/src/apply.jl b/src/apply.jl
index d38f04f9..6a55f45f 100644
--- a/src/apply.jl
+++ b/src/apply.jl
@@ -200,7 +200,7 @@ function ITensors.apply(
   v⃗ = neighbor_vertices(ψ, o)
   if length(v⃗) == 1
     if ortho
-      ψ = orthogonalize(ψ, v⃗[1])
+      ψ = tree_orthogonalize(ψ, v⃗[1])
     end
     oψᵥ = apply(o, ψ[v⃗[1]])
     if normalize
@@ -215,7 +215,7 @@ function ITensors.apply(
       error("Vertices where the gates are being applied must be neighbors for now.")
     end
     if ortho
-      ψ = orthogonalize(ψ, v⃗[1])
+      ψ = tree_orthogonalize(ψ, v⃗[1])
     end
     if variational_optimization_only || !is_product_env
       ψᵥ₁, ψᵥ₂ = full_update_bp(
diff --git a/src/tebd.jl b/src/tebd.jl
index edf5a188..d1d96017 100644
--- a/src/tebd.jl
+++ b/src/tebd.jl
@@ -23,7 +23,7 @@ function tebd(
     ψ = apply(u⃗, ψ; cutoff, maxdim, normalize=true, ortho, kwargs...)
     if ortho
       for v in vertices(ψ)
-        ψ = orthogonalize(ψ, v)
+        ψ = tree_orthogonalize(ψ, v)
       end
     end
   end
diff --git a/src/treetensornetworks/abstracttreetensornetwork.jl b/src/treetensornetworks/abstracttreetensornetwork.jl
index a08fe27d..4635ac8a 100644
--- a/src/treetensornetworks/abstracttreetensornetwork.jl
+++ b/src/treetensornetworks/abstracttreetensornetwork.jl
@@ -50,6 +50,10 @@ function ITensorMPS.orthogonalize(ttn::AbstractTTN, region; kwargs...)
   return orthogonalize(ttn, [region]; kwargs...)
 end
 
+function tree_orthogonalize(ttn::AbstractTTN, args...; kwargs...)
+  return orthogonalize(ttn, args...; kwargs...)
+end
+
 # 
 # Truncation
 # 

From 9a3ca1c88ef8ec0e9255822250f2470d5350e3a0 Mon Sep 17 00:00:00 2001
From: Joey <joeytindall94@googlemail.com>
Date: Fri, 8 Nov 2024 15:52:49 -0500
Subject: [PATCH 14/14] Improvements

---
 src/abstractitensornetwork.jl                 | 20 ++++++++--------
 src/solvers/sweep_plans/sweep_plans.jl        | 23 ++++++++-----------
 .../abstracttreetensornetwork.jl              |  2 +-
 3 files changed, 21 insertions(+), 24 deletions(-)

diff --git a/src/abstractitensornetwork.jl b/src/abstractitensornetwork.jl
index f7c148d0..fc0edce4 100644
--- a/src/abstractitensornetwork.jl
+++ b/src/abstractitensornetwork.jl
@@ -585,16 +585,16 @@ function LinearAlgebra.factorize(tn::AbstractITensorNetwork, edge::Pair; kwargs.
 end
 
 # For ambiguity error; TODO: decide whether to use graph mutating methods when resulting graph is unchanged?
-function orthogonalize_path(tn::AbstractITensorNetwork, edge::AbstractEdge; kwargs...)
-  return orthogonalize_path(tn, [edge]; kwargs...)
+function orthogonalize_walk(tn::AbstractITensorNetwork, edge::AbstractEdge; kwargs...)
+  return orthogonalize_walk(tn, [edge]; kwargs...)
 end
 
-function orthogonalize_path(tn::AbstractITensorNetwork, edge::Pair; kwargs...)
-  return orthogonalize_path(tn, edgetype(tn)(edge); kwargs...)
+function orthogonalize_walk(tn::AbstractITensorNetwork, edge::Pair; kwargs...)
+  return orthogonalize_walk(tn, edgetype(tn)(edge); kwargs...)
 end
 
 # For ambiguity error; TODO: decide whether to use graph mutating methods when resulting graph is unchanged?
-function orthogonalize_path(
+function orthogonalize_walk(
   tn::AbstractITensorNetwork, edges::Vector{<:AbstractEdge}; kwargs...
 )
   # tn = factorize(tn, edge; kwargs...)
@@ -612,18 +612,18 @@ function orthogonalize_path(
   return tn
 end
 
-function orthogonalize_path(tn::AbstractITensorNetwork, edges::Vector{<:Pair}; kwargs...)
-  return orthogonalize_path(tn, edgetype(tn).(edges); kwargs...)
+function orthogonalize_walk(tn::AbstractITensorNetwork, edges::Vector{<:Pair}; kwargs...)
+  return orthogonalize_walk(tn, edgetype(tn).(edges); kwargs...)
 end
 
 # Orthogonalize an ITensorNetwork towards a region, treating
 # the network as a tree spanned by a spanning tree.
 function tree_orthogonalize(ψ::AbstractITensorNetwork, region::Vector)
-  region_centre =
+  region_center =
     length(region) != 1 ? first(center(steiner_tree(ψ, region))) : only(region)
-  path = post_order_dfs_edges(bfs_tree(ψ, region_centre), region_centre)
+  path = post_order_dfs_edges(bfs_tree(ψ, region_center), region_center)
   path = filter(e -> !((src(e) ∈ region) && (dst(e) ∈ region)), path)
-  return orthogonalize_path(ψ, path)
+  return orthogonalize_walk(ψ, path)
 end
 
 function tree_orthogonalize(ψ::AbstractITensorNetwork, region)
diff --git a/src/solvers/sweep_plans/sweep_plans.jl b/src/solvers/sweep_plans/sweep_plans.jl
index 384cf8dc..52915e2b 100644
--- a/src/solvers/sweep_plans/sweep_plans.jl
+++ b/src/solvers/sweep_plans/sweep_plans.jl
@@ -70,20 +70,17 @@ function forward_sweep(
   kwargs...,
 )
   edges = post_order_dfs_edges(graph, root_vertex)
-  regions = collect(
-    flatten(map(i -> forward_region(edges, i; region_kwargs, kwargs...), eachindex(edges)))
-  )
+  regions = map(eachindex(edges)) do i
+    forward_region(edges, i; region_kwargs, kwargs...)
+  end
+  regions = collect(flatten(regions))
   if reverse_step
-    reverse_regions = collect(
-      flatten(
-        map(
-          i -> reverse_region(
-            edges, i; reverse_edge=reverse_edges, region_kwargs=reverse_kwargs, kwargs...
-          ),
-          eachindex(edges),
-        ),
-      ),
-    )
+    reverse_regions = map(eachindex(edges)) do i
+      reverse_region(
+        edges, i; reverse_edge=reverse_edges, region_kwargs=reverse_kwargs, kwargs...
+      )
+    end
+    reverse_regions = collect(flatten(reverse_regions))
     _check_reverse_sweeps(regions, reverse_regions, graph; kwargs...)
     regions = interleave(regions, reverse_regions)
   end
diff --git a/src/treetensornetworks/abstracttreetensornetwork.jl b/src/treetensornetworks/abstracttreetensornetwork.jl
index 4635ac8a..8815b33f 100644
--- a/src/treetensornetworks/abstracttreetensornetwork.jl
+++ b/src/treetensornetworks/abstracttreetensornetwork.jl
@@ -41,7 +41,7 @@ function ITensorMPS.orthogonalize(ttn::AbstractTTN, region::Vector; kwargs...)
   path = post_order_dfs_edges(st, first(region))
   path = filter(e -> !((src(e) ∈ region) && (dst(e) ∈ region)), path)
   if !isempty(path)
-    ttn = typeof(ttn)(orthogonalize_path(ITensorNetwork(ttn), path; kwargs...))
+    ttn = typeof(ttn)(orthogonalize_walk(ITensorNetwork(ttn), path; kwargs...))
   end
   return set_ortho_region(ttn, region)
 end