diff --git a/src/solvers/abstract_problem.jl b/src/solvers/abstract_problem.jl
index e3bde03b..6e629090 100644
--- a/src/solvers/abstract_problem.jl
+++ b/src/solvers/abstract_problem.jl
@@ -1,4 +1,4 @@
 
 abstract type AbstractProblem end
 
-set_truncation_info(P::AbstractProblem, args...; kws...) = P
+set_truncation_info!(P::AbstractProblem, args...; kws...) = P
diff --git a/src/solvers/adapters.jl b/src/solvers/adapters.jl
index 7c033d8e..e1c58d15 100644
--- a/src/solvers/adapters.jl
+++ b/src/solvers/adapters.jl
@@ -1,32 +1,46 @@
+"""
+  struct PauseAfterIncrement{S<:AbstractNetworkIterator}
 
-#
-# TupleRegionIterator
-#
-# Adapts outputs to be (region, region_kwargs) tuples
-#
-# More generic design? maybe just assuming RegionIterator
-# or its outputs implement some interface function that
-# generates each tuple?
-#
-
-mutable struct TupleRegionIterator{RegionIter}
-  region_iterator::RegionIter
+Iterator wrapper whos `compute!` function simply returns itself, doing nothing in the 
+process. This allows one to manually call a custom `compute!` or insert their own code it in
+the loop body in place of `compute!`.
+"""
+struct IncrementOnly{S<:AbstractNetworkIterator} <: AbstractNetworkIterator
+  parent::S
 end
 
-region_iterator(T::TupleRegionIterator) = T.region_iterator
+islaststep(adapter::IncrementOnly) = islaststep(adapter.parent)
+state(adapter::IncrementOnly) = state(adapter.parent)
+increment!(adapter::IncrementOnly) = increment!(adapter.parent)
+compute!(adapter::IncrementOnly) = adapter
 
-function Base.iterate(T::TupleRegionIterator, which=1)
-  state = iterate(region_iterator(T), which)
-  isnothing(state) && return nothing
-  (current_region, region_kwargs) = current_region_plan(region_iterator(T))
-  return (current_region, region_kwargs), last(state)
-end
+IncrementOnly(adapter::IncrementOnly) = adapter
 
 """
-  region_tuples(R::RegionIterator)
+  struct EachRegion{SweepIterator} <: AbstractNetworkIterator
 
-The `region_tuples` adapter converts a RegionIterator into an 
-iterator which outputs a tuple of the form (current_region, current_region_kwargs)
-at each step.
+Adapter that flattens each region iterator in the parent sweep iterator into a single
+iterator.
 """
-region_tuples(R::RegionIterator) = TupleRegionIterator(R)
+struct EachRegion{SI<:SweepIterator} <: AbstractNetworkIterator
+  parent::SI
+end
+
+# In keeping with Julia convention.
+eachregion(iter::SweepIterator) = EachRegion(iter)
+
+# Essential definitions
+function islaststep(adapter::EachRegion)
+  region_iter = region_iterator(adapter.parent)
+  return islaststep(adapter.parent) && islaststep(region_iter)
+end
+function increment!(adapter::EachRegion)
+  region_iter = region_iterator(adapter.parent)
+  islaststep(region_iter) ? increment!(adapter.parent) : increment!(region_iter)
+  return adapter
+end
+function compute!(adapter::EachRegion)
+  region_iter = region_iterator(adapter.parent)
+  compute!(region_iter)
+  return adapter
+end
diff --git a/src/solvers/applyexp.jl b/src/solvers/applyexp.jl
index 8bc070e3..3b91ff8e 100644
--- a/src/solvers/applyexp.jl
+++ b/src/solvers/applyexp.jl
@@ -1,5 +1,4 @@
 using Printf: @printf
-using Accessors: @set
 
 @kwdef mutable struct ApplyExpProblem{State} <: AbstractProblem
   operator
@@ -11,66 +10,69 @@ operator(A::ApplyExpProblem) = A.operator
 state(A::ApplyExpProblem) = A.state
 current_exponent(A::ApplyExpProblem) = A.current_exponent
 function current_time(A::ApplyExpProblem)
-  t = im*A.current_exponent
+  t = im * A.current_exponent
   return iszero(imag(t)) ? real(t) : t
 end
 
-set_operator(A::ApplyExpProblem, operator) = (@set A.operator = operator)
-set_state(A::ApplyExpProblem, state) = (@set A.state = state)
-set_current_exponent(A::ApplyExpProblem, exponent) = (@set A.current_exponent = exponent)
-
-function region_plan(A::ApplyExpProblem; nsites, time_step, sweep_kwargs...)
-  return applyexp_regions(state(A), time_step; nsites, sweep_kwargs...)
+# Rename region_plan
+function region_plan(A::ApplyExpProblem; nsites, exponent_step, sweep_kwargs...)
+  # The `exponent_step` kwarg for the `update!` function needs some pre-processing.
+  return applyexp_regions(state(A), exponent_step; nsites, sweep_kwargs...)
 end
 
-function update(
-  prob::ApplyExpProblem,
-  local_state,
-  region_iterator;
+function update!(
+  region_iter::RegionIterator{<:ApplyExpProblem},
+  local_state;
   nsites,
   exponent_step,
   solver=runge_kutta_solver,
-  outputlevel,
-  kws...,
 )
-  iszero(abs(exponent_step)) && return prob, local_state
+  prob = problem(region_iter)
+
+  if iszero(abs(exponent_step))
+    return region_iter, local_state
+  end
 
-  local_state, info = solver(
-    x->optimal_map(operator(prob), x), exponent_step, local_state; kws...
+  solver_kwargs = region_kwargs(solver, region_iter)
+
+  local_state, _ = solver(
+    x -> optimal_map(operator(prob), x), exponent_step, local_state; solver_kwargs...
   )
-  if nsites==1
-    curr_reg = current_region(region_iterator)
-    next_reg = next_region(region_iterator)
+  if nsites == 1
+    curr_reg = current_region(region_iter)
+    next_reg = next_region(region_iter)
     if !isnothing(next_reg) && next_reg != curr_reg
       next_edge = first(edge_sequence_between_regions(state(prob), curr_reg, next_reg))
       v1, v2 = src(next_edge), dst(next_edge)
       psi = copy(state(prob))
       psi[v1], R = qr(local_state, uniqueinds(local_state, psi[v2]))
-      shifted_operator = position(operator(prob), psi, NamedEdge(v1=>v2))
-      R_t, _ = solver(x->optimal_map(shifted_operator, x), -exponent_step, R; kws...)
-      local_state = psi[v1]*R_t
+      shifted_operator = position(operator(prob), psi, NamedEdge(v1 => v2))
+      R_t, _ = solver(
+        x -> optimal_map(shifted_operator, x), -exponent_step, R; solver_kwargs...
+      )
+      local_state = psi[v1] * R_t
     end
   end
 
-  prob = set_current_exponent(prob, current_exponent(prob)+exponent_step)
+  prob.current_exponent += exponent_step
 
-  return prob, local_state
+  return region_iter, local_state
 end
 
-function sweep_callback(
-  problem::ApplyExpProblem;
+function default_sweep_callback(
+  sweep_iterator::SweepIterator{<:ApplyExpProblem};
   exponent_description="exponent",
-  outputlevel,
-  sweep,
-  nsweeps,
+  outputlevel=0,
   process_time=identity,
-  kws...,
 )
   if outputlevel >= 1
+    the_problem = problem(sweep_iterator)
     @printf(
-      "  Current %s = %s, ", exponent_description, process_time(current_exponent(problem))
+      "  Current %s = %s, ",
+      exponent_description,
+      process_time(current_exponent(the_problem))
     )
-    @printf("maxlinkdim=%d", maxlinkdim(state(problem)))
+    @printf("maxlinkdim=%d", maxlinkdim(state(the_problem)))
     println()
     flush(stdout)
   end
@@ -79,19 +81,20 @@ end
 function applyexp(
   init_prob::AbstractProblem,
   exponents;
-  extract_kwargs=(;),
-  update_kwargs=(;),
-  insert_kwargs=(;),
-  outputlevel=0,
-  nsites=1,
+  sweep_callback=default_sweep_callback,
   order=4,
-  kws...,
+  nsites=2,
+  sweep_kwargs...,
 )
   exponent_steps = diff([zero(eltype(exponents)); exponents])
-  sweep_kws = (; outputlevel, extract_kwargs, insert_kwargs, nsites, order, update_kwargs)
-  kws_array = [(; sweep_kws..., time_step=t) for t in exponent_steps]
-  sweep_iter = sweep_iterator(init_prob, kws_array)
-  converged_prob = sweep_solve(sweep_iter; outputlevel, kws...)
+
+  kws_array = [
+    (; order, nsites, sweep_kwargs..., exponent_step) for exponent_step in exponent_steps
+  ]
+  sweep_iter = SweepIterator(init_prob, kws_array)
+
+  converged_prob = problem(sweep_solve!(sweep_callback, sweep_iter))
+
   return state(converged_prob)
 end
 
@@ -111,11 +114,10 @@ function time_evolve(
   time_points,
   init_state;
   process_time=process_real_times,
-  sweep_callback=(
-    a...; k...
-  )->sweep_callback(a...; exponent_description="time", process_time, k...),
-  kws...,
+  sweep_callback=iter ->
+    default_sweep_callback(iter; exponent_description="time", process_time),
+  sweep_kwargs...,
 )
-  exponents = [-im*t for t in time_points]
-  return applyexp(operator, exponents, init_state; sweep_callback, kws...)
+  exponents = [-im * t for t in time_points]
+  return applyexp(operator, exponents, init_state; sweep_callback, sweep_kwargs...)
 end
diff --git a/src/solvers/eigsolve.jl b/src/solvers/eigsolve.jl
index 6916406a..907e71f5 100644
--- a/src/solvers/eigsolve.jl
+++ b/src/solvers/eigsolve.jl
@@ -1,4 +1,3 @@
-using Accessors: @set
 using Printf: @printf
 using ITensors: truncerror
 
@@ -14,42 +13,43 @@ state(E::EigsolveProblem) = E.state
 operator(E::EigsolveProblem) = E.operator
 max_truncerror(E::EigsolveProblem) = E.max_truncerror
 
-set_operator(E::EigsolveProblem, operator) = (@set E.operator = operator)
-set_eigenvalue(E::EigsolveProblem, eigenvalue) = (@set E.eigenvalue = eigenvalue)
-set_state(E::EigsolveProblem, state) = (@set E.state = state)
-set_max_truncerror(E::EigsolveProblem, truncerror) = (@set E.max_truncerror = truncerror)
-
-function set_truncation_info(E::EigsolveProblem; spectrum=nothing)
+function set_truncation_info!(E::EigsolveProblem; spectrum=nothing)
   if !isnothing(spectrum)
-    E = set_max_truncerror(E, max(max_truncerror(E), truncerror(spectrum)))
+    E.max_truncerror = max(max_truncerror(E), truncerror(spectrum))
   end
   return E
 end
 
-function update(
-  prob::EigsolveProblem,
-  local_state,
-  region_iterator;
-  outputlevel,
+function update!(
+  region_iter::RegionIterator{<:EigsolveProblem},
+  local_state;
+  outputlevel=0,
   solver=eigsolve_solver,
-  kws...,
 )
-  eigval, local_state = solver(ψ->optimal_map(operator(prob), ψ), local_state; kws...)
-  prob = set_eigenvalue(prob, eigval)
+  prob = problem(region_iter)
+
+  eigval, local_state = solver(
+    ψ -> optimal_map(operator(prob), ψ), local_state; region_kwargs(solver, region_iter)...
+  )
+
+  prob.eigenvalue = eigval
+
   if outputlevel >= 2
-    @printf(
-      "  Region %s: energy = %.12f\n", current_region(region_iterator), eigenvalue(prob)
-    )
+    @printf("  Region %s: energy = %.12f\n", current_region(region_iter), eigenvalue(prob))
   end
-  return prob, local_state
+  return region_iter, local_state
 end
 
-function sweep_callback(problem::EigsolveProblem; outputlevel, sweep, nsweeps, kws...)
+function default_sweep_callback(
+  sweep_iterator::SweepIterator{<:EigsolveProblem}; outputlevel=0
+)
   if outputlevel >= 1
-    if nsweeps >= 10
-      @printf("After sweep %02d/%d ", sweep, nsweeps)
+    nsweeps = length(sweep_iterator)
+    current_sweep = sweep_iterator.which_sweep
+    if length(sweep_iterator) >= 10
+      @printf("After sweep %02d/%d ", current_sweep, nsweeps)
     else
-      @printf("After sweep %d/%d ", sweep, nsweeps)
+      @printf("After sweep %d/%d ", current_sweep, nsweeps)
     end
     @printf("eigenvalue=%.12f", eigenvalue(problem))
     @printf(" maxlinkdim=%d", maxlinkdim(state(problem)))
@@ -60,24 +60,22 @@ function sweep_callback(problem::EigsolveProblem; outputlevel, sweep, nsweeps, k
 end
 
 function eigsolve(
-  operator,
-  init_state;
-  nsweeps,
-  nsites=1,
-  outputlevel=0,
-  extract_kwargs=(;),
-  update_kwargs=(;),
-  insert_kwargs=(;),
-  kws...,
+  operator, init_state; nsweeps, nsites=1, outputlevel=0, factorize_kwargs, sweep_kwargs...
 )
   init_prob = EigsolveProblem(;
     state=align_indices(init_state), operator=ProjTTN(align_indices(operator))
   )
-  sweep_iter = sweep_iterator(
-    init_prob, nsweeps; nsites, outputlevel, extract_kwargs, update_kwargs, insert_kwargs
+  sweep_iter = SweepIterator(
+    init_prob,
+    nsweeps;
+    nsites,
+    outputlevel,
+    factorize_kwargs,
+    subspace_expand!_kwargs=(; eigen_kwargs=factorize_kwargs),
+    sweep_kwargs...,
   )
-  prob = sweep_solve(sweep_iter; outputlevel, kws...)
+  prob = problem(sweep_solve!(sweep_iter))
   return eigenvalue(prob), state(prob)
 end
 
-dmrg(args...; kws...) = eigsolve(args...; kws...)
+dmrg(operator, init_state; kwargs...) = eigsolve(operator, init_state; kwargs...)
diff --git a/src/solvers/extract.jl b/src/solvers/extract.jl
index 011058af..629b70f2 100644
--- a/src/solvers/extract.jl
+++ b/src/solvers/extract.jl
@@ -1,12 +1,16 @@
-function extract(problem, region_iterator; sweep, trunc=(;), kws...)
-  trunc = truncation_parameters(sweep; trunc...)
-  region = current_region(region_iterator)
-  psi = orthogonalize(state(problem), region)
+function extract!(region_iter::RegionIterator; subspace_algorithm="nothing")
+  prob = problem(region_iter)
+  region = current_region(region_iter)
+
+  psi = orthogonalize(state(prob), region)
   local_state = prod(psi[v] for v in region)
-  problem = set_state(problem, psi)
-  problem, local_state = subspace_expand(
-    problem, local_state, region_iterator; sweep, trunc, kws...
-  )
-  shifted_operator = position(operator(problem), state(problem), region)
-  return set_operator(problem, shifted_operator), local_state
+
+  prob.state = psi
+
+  _, local_state = subspace_expand!(region_iter, local_state; subspace_algorithm)
+  shifted_operator = position(operator(prob), state(prob), region)
+
+  prob.operator = shifted_operator
+
+  return region_iter, local_state
 end
diff --git a/src/solvers/fitting.jl b/src/solvers/fitting.jl
index e04df71d..c03852e6 100644
--- a/src/solvers/fitting.jl
+++ b/src/solvers/fitting.jl
@@ -1,9 +1,7 @@
-using Accessors: @set
 using Graphs: vertices
 using NamedGraphs: AbstractNamedGraph, NamedEdge
 using NamedGraphs.PartitionedGraphs: partitionedges
 using Printf: @printf
-using ConstructionBase: setproperties
 
 @kwdef mutable struct FittingProblem{State<:AbstractBeliefPropagationCache} <:
                       AbstractProblem
@@ -18,19 +16,18 @@ ket_graph(F::FittingProblem) = F.ket_graph
 overlap(F::FittingProblem) = F.overlap
 gauge_region(F::FittingProblem) = F.gauge_region
 
-set_state(F::FittingProblem, state) = (@set F.state = state)
-set_overlap(F::FittingProblem, overlap) = (@set F.overlap = overlap)
-
 function ket(F::FittingProblem)
   ket_vertices = vertices(ket_graph(F))
   return first(induced_subgraph(tensornetwork(state(F)), ket_vertices))
 end
 
-function extract(problem::FittingProblem, region_iterator; sweep, kws...)
-  region = current_region(region_iterator)
-  prev_region = gauge_region(problem)
-  tn = state(problem)
-  path = edge_sequence_between_regions(ket_graph(problem), prev_region, region)
+function extract!(region_iter::RegionIterator{<:FittingProblem})
+  prob = problem(region_iter)
+
+  region = current_region(region_iter)
+  prev_region = gauge_region(prob)
+  tn = state(prob)
+  path = edge_sequence_between_regions(ket_graph(prob), prev_region, region)
   tn = gauge_walk(Algorithm("orthogonalize"), tn, path)
   pe_path = partitionedges(partitioned_tensornetwork(tn), path)
   tn = update(
@@ -40,16 +37,28 @@ function extract(problem::FittingProblem, region_iterator; sweep, kws...)
   sequence = contraction_sequence(local_tensor; alg="optimal")
   local_tensor = dag(contract(local_tensor; sequence))
   #problem, local_tensor = subspace_expand(problem, local_tensor, region; sweep, kws...)
-  return setproperties(problem; state=tn, gauge_region=region), local_tensor
+
+  prob.state = tn
+  prob.gauge_region = region
+
+  return region_iter, local_tensor
 end
 
-function update(F::FittingProblem, local_tensor, region; outputlevel, kws...)
+function update!(
+  region_iter::RegionIterator{<:FittingProblem}, local_tensor; outputlevel=0
+)
+  F = problem(region_iter)
+
+  region = current_region(region_iter)
+
   n = (local_tensor * dag(local_tensor))[]
-  F = set_overlap(F, n / sqrt(n))
+  F.overlap = n / sqrt(n)
+
   if outputlevel >= 2
     @printf("  Region %s: squared overlap = %.12f\n", region, overlap(F))
   end
-  return F, local_tensor
+
+  return region_iter, local_tensor
 end
 
 function region_plan(F::FittingProblem; nsites, sweep_kwargs...)
@@ -62,11 +71,9 @@ function fit_tensornetwork(
   nsweeps=25,
   nsites=1,
   outputlevel=0,
-  extract_kwargs=(;),
-  update_kwargs=(;),
-  insert_kwargs=(;),
   normalize=true,
-  kws...,
+  factorize_kwargs,
+  extra_sweep_kwargs...,
 )
   bpc = BeliefPropagationCache(overlap_network, args...)
   ket_graph = first(
@@ -76,11 +83,15 @@ function fit_tensornetwork(
     ket_graph, state=bpc, gauge_region=collect(vertices(ket_graph))
   )
 
-  insert_kwargs = (; insert_kwargs..., normalize, set_orthogonal_region=false)
-  common_sweep_kwargs = (; nsites, outputlevel, update_kwargs, insert_kwargs)
-  kwargs_array = [(; common_sweep_kwargs..., sweep=s) for s in 1:nsweeps]
-  sweep_iter = sweep_iterator(init_prob, kwargs_array)
-  converged_prob = sweep_solve(sweep_iter; outputlevel, kws...)
+  insert!_kwargs = (; normalize, set_orthogonal_region=false)
+  update!_kwargs = (; outputlevel)
+
+  sweep_kwargs = (; nsites, outputlevel, update!_kwargs, insert!_kwargs, factorize_kwargs)
+  kwargs_array = [(; sweep_kwargs..., extra_sweep_kwargs..., sweep) for sweep in 1:nsweeps]
+
+  sweep_iter = SweepIterator(init_prob, kwargs_array)
+  converged_prob = problem(sweep_solve!(sweep_iter))
+
   return rename_vertices(inv_vertex_map(overlap_network), ket(converged_prob))
 end
 
@@ -98,14 +109,15 @@ end
 #end
 
 function ITensors.apply(
-  A::ITensorNetwork,
-  x::ITensorNetwork;
-  maxdim=default_maxdim(),
-  cutoff=default_cutoff(),
-  kwargs...,
+  A::AbstractITensorNetwork,
+  x::AbstractITensorNetwork;
+  maxdim=typemax(Int),
+  cutoff=0.0,
+  sweep_kwargs...,
 )
   init_state = ITensorNetwork(v -> inds -> delta(inds), siteinds(x); link_space=maxdim)
   overlap_network = inner_network(x, A, init_state)
-  insert_kwargs = (; trunc=(; cutoff, maxdim))
-  return fit_tensornetwork(overlap_network; insert_kwargs, kwargs...)
+  return fit_tensornetwork(
+    overlap_network; factorize_kwargs=(; maxdim, cutoff), sweep_kwargs...
+  )
 end
diff --git a/src/solvers/insert.jl b/src/solvers/insert.jl
index b3c60645..87ffaf6d 100644
--- a/src/solvers/insert.jl
+++ b/src/solvers/insert.jl
@@ -1,28 +1,23 @@
 using NamedGraphs: edgetype
 
-function insert(
-  problem,
-  local_tensor,
-  region_iterator;
-  normalize=false,
-  set_orthogonal_region=true,
-  sweep,
-  trunc=(;),
-  outputlevel=0,
-  kws...,
-)
-  trunc = truncation_parameters(sweep; trunc...)
-  region = current_region(region_iterator)
-  psi = copy(state(problem))
+function insert!(region_iter, local_tensor; normalize=false, set_orthogonal_region=true)
+  prob = problem(region_iter)
+
+  region = current_region(region_iter)
+  psi = copy(state(prob))
   if length(region) == 1
     C = local_tensor
   elseif length(region) == 2
     e = edgetype(psi)(first(region), last(region))
     indsTe = inds(psi[first(region)])
     tags = ITensors.tags(psi, e)
-    U, C, spectrum = factorize(local_tensor, indsTe; tags, trunc...)
+
+    U, C, spectrum = factorize(
+      local_tensor, indsTe; tags, region_kwargs(factorize, region_iter)...
+    )
+
     @preserve_graph psi[first(region)] = U
-    problem = set_truncation_info(problem; spectrum)
+    prob = set_truncation_info!(prob; spectrum)
   else
     error("Region of length $(length(region)) not currently supported")
   end
@@ -30,6 +25,8 @@ function insert(
   @preserve_graph psi[v] = C
   psi = set_orthogonal_region ? set_ortho_region(psi, [v]) : psi
   normalize && @preserve_graph psi[v] = psi[v] / norm(psi[v])
-  problem = set_state(problem, psi)
-  return problem
+
+  prob.state = psi
+
+  return region_iter
 end
diff --git a/src/solvers/iterators.jl b/src/solvers/iterators.jl
index 023818ba..16497f0e 100644
--- a/src/solvers/iterators.jl
+++ b/src/solvers/iterators.jl
@@ -1,100 +1,159 @@
+"""
+  abstract type AbstractNetworkIterator
+
+A stateful iterator with two states: `increment!` and `compute!`. Each iteration begins
+with a call to `increment!` before executing `compute!`, however the initial call to
+`iterate` skips the `increment!` call as it is assumed the iterator is initalized such that 
+this call is implict. Termination of the iterator is controlled by the function `done`.
+"""
+abstract type AbstractNetworkIterator end
+
+# We use greater than or equals here as we increment the state at the start of the iteration
+islaststep(iterator::AbstractNetworkIterator) = state(iterator) >= length(iterator)
+
+function Base.iterate(iterator::AbstractNetworkIterator, init=true)
+  islaststep(iterator) && return nothing
+  # We seperate increment! from step! and demand that any AbstractNetworkIterator *must*
+  # define a method for increment! This way we avoid cases where one may wish to nest
+  # calls to different step! methods accidentaly incrementing multiple times.
+  init || increment!(iterator)
+  rv = compute!(iterator)
+  return rv, false
+end
+
+function increment! end
+compute!(iterator::AbstractNetworkIterator) = iterator
+
+step!(iterator::AbstractNetworkIterator) = step!(identity, iterator)
+function step!(f, iterator::AbstractNetworkIterator)
+  compute!(iterator)
+  f(iterator)
+  increment!(iterator)
+  return iterator
+end
+
 #
-# SweepIterator
+# RegionIterator
 #
+"""
+  struct RegionIterator{Problem, RegionPlan} <: AbstractNetworkIterator
+"""
+mutable struct RegionIterator{Problem,RegionPlan} <: AbstractNetworkIterator
+  problem::Problem
+  region_plan::RegionPlan
+  which_region::Int
+  const which_sweep::Int
+  function RegionIterator(problem::P, region_plan::R, sweep::Int) where {P,R}
+    return new{P,R}(problem, region_plan, 1, sweep)
+  end
+end
 
-mutable struct SweepIterator
-  sweep_kws
-  region_iter
-  which_sweep::Int
+function RegionIterator(problem; sweep, sweep_kwargs...)
+  plan = region_plan(problem; sweep_kwargs...)
+  return RegionIterator(problem, plan, sweep)
 end
 
-problem(S::SweepIterator) = problem(S.region_iter)
+function new_region_iterator(iterator::RegionIterator; sweep_kwargs...)
+  return RegionIterator(iterator.problem; sweep_kwargs...)
+end
 
-Base.length(S::SweepIterator) = length(S.sweep_kws)
+state(region_iter::RegionIterator) = region_iter.which_region
+Base.length(region_iter::RegionIterator) = length(region_iter.region_plan)
 
-function Base.iterate(S::SweepIterator, which=nothing)
-  if isnothing(which)
-    sweep_kws_state = iterate(S.sweep_kws)
-  else
-    sweep_kws_state = iterate(S.sweep_kws, which)
-  end
-  isnothing(sweep_kws_state) && return nothing
-  current_sweep_kws, next = sweep_kws_state
+problem(region_iter::RegionIterator) = region_iter.problem
 
-  if !isnothing(which)
-    S.region_iter = region_iterator(
-      problem(S.region_iter); sweep=S.which_sweep, current_sweep_kws...
-    )
-  end
-  S.which_sweep += 1
-  return S.region_iter, next
+function current_region_plan(region_iter::RegionIterator)
+  return region_iter.region_plan[region_iter.which_region]
+end
+
+function current_region(region_iter::RegionIterator)
+  region, _ = current_region_plan(region_iter)
+  return region
+end
+
+function region_kwargs(region_iter::RegionIterator)
+  _, kwargs = current_region_plan(region_iter)
+  return kwargs
+end
+function region_kwargs(f::Function, iter::RegionIterator)
+  return get(region_kwargs(iter), Symbol(f, :_kwargs), (;))
 end
 
-function sweep_iterator(problem, sweep_kws)
-  region_iter = region_iterator(problem; sweep=1, first(sweep_kws)...)
-  return SweepIterator(sweep_kws, region_iter, 1)
+function prev_region(region_iter::RegionIterator)
+  state(region_iter) <= 1 && return nothing
+  prev, _ = region_iter.region_plan[region_iter.which_region - 1]
+  return prev
 end
 
-function sweep_iterator(problem, nsweeps::Integer; sweep_kws...)
-  return sweep_iterator(problem, Iterators.repeated(sweep_kws, nsweeps))
+function next_region(region_iter::RegionIterator)
+  islaststep(region_iter) && return nothing
+  next, _ = region_iter.region_plan[region_iter.which_region + 1]
+  return next
 end
 
 #
-# RegionIterator
+# Functions associated with RegionIterator
 #
+function increment!(region_iter::RegionIterator)
+  region_iter.which_region += 1
+  return region_iter
+end
 
-@kwdef mutable struct RegionIterator{Problem,RegionPlan}
-  problem::Problem
-  region_plan::RegionPlan
-  which_region::Int = 1
+function compute!(iter::RegionIterator)
+  _, local_state = extract!(iter; region_kwargs(extract!, iter)...)
+  _, local_state = update!(iter, local_state; region_kwargs(update!, iter)...)
+  insert!(iter, local_state; region_kwargs(insert!, iter)...)
+
+  return iter
 end
 
-problem(R::RegionIterator) = R.problem
-current_region_plan(R::RegionIterator) = R.region_plan[R.which_region]
-current_region(R::RegionIterator) = current_region_plan(R)[1]
-region_kwargs(R::RegionIterator) = current_region_plan(R)[2]
-function previous_region(R::RegionIterator)
-  R.which_region==1 ? nothing : R.region_plan[R.which_region - 1][1]
+region_plan(problem; sweep_kwargs...) = euler_sweep(state(problem); sweep_kwargs...)
+
+#
+# SweepIterator
+#
+
+mutable struct SweepIterator{Problem,Iter} <: AbstractNetworkIterator
+  region_iter::RegionIterator{Problem}
+  sweep_kwargs::Iterators.Stateful{Iter}
+  which_sweep::Int
+  function SweepIterator(problem::Prob, sweep_kwargs::Iter) where {Prob,Iter}
+    stateful_sweep_kwargs = Iterators.Stateful(sweep_kwargs)
+    first_kwargs, _ = Iterators.peel(stateful_sweep_kwargs)
+    region_iter = RegionIterator(problem; sweep=1, first_kwargs...)
+    return new{Prob,Iter}(region_iter, stateful_sweep_kwargs, 1)
+  end
 end
-function next_region(R::RegionIterator)
-  R.which_region==length(R.region_plan) ? nothing : R.region_plan[R.which_region + 1][1]
+
+islaststep(sweep_iter::SweepIterator) = isnothing(peek(sweep_iter.sweep_kwargs))
+
+region_iterator(sweep_iter::SweepIterator) = sweep_iter.region_iter
+problem(sweep_iter::SweepIterator) = problem(region_iterator(sweep_iter))
+
+state(sweep_iter::SweepIterator) = sweep_iter.which_sweep
+Base.length(sweep_iter::SweepIterator) = length(sweep_iter.sweep_kwargs)
+function increment!(sweep_iter::SweepIterator)
+  sweep_iter.which_sweep += 1
+  sweep_kwargs, _ = Iterators.peel(sweep_iter.sweep_kwargs)
+  update_region_iterator!(sweep_iter; sweep_kwargs...)
+  return sweep_iter
 end
-is_last_region(R::RegionIterator) = isnothing(next_region(R))
 
-function Base.iterate(R::RegionIterator, which=1)
-  R.which_region = which
-  region_plan_state = iterate(R.region_plan, which)
-  isnothing(region_plan_state) && return nothing
-  (current_region, region_kwargs), next = region_plan_state
-  R.problem = region_step(problem(R), R; region_kwargs...)
-  return R, next
+function update_region_iterator!(iterator::SweepIterator; kwargs...)
+  sweep = state(iterator)
+  iterator.region_iter = new_region_iterator(iterator.region_iter; sweep, kwargs...)
+  return iterator
 end
 
-#
-# Functions associated with RegionIterator
-#
+function compute!(sweep_iter::SweepIterator)
+  for _ in sweep_iter.region_iter
+    # TODO: Is it sensible to execute the default region callback function?
+  end
+end
 
-function region_iterator(problem; sweep_kwargs...)
-  return RegionIterator(; problem, region_plan=region_plan(problem; sweep_kwargs...))
-end
-
-function region_step(
-  problem,
-  region_iterator;
-  extract_kwargs=(;),
-  update_kwargs=(;),
-  insert_kwargs=(;),
-  sweep,
-  kws...,
-)
-  problem, local_state = extract(problem, region_iterator; extract_kwargs..., sweep, kws...)
-  problem, local_state = update(
-    problem, local_state, region_iterator; update_kwargs..., kws...
-  )
-  problem = insert(problem, local_state, region_iterator; sweep, insert_kwargs..., kws...)
-  return problem
-end
-
-function region_plan(problem; kws...)
-  return euler_sweep(state(problem); kws...)
+# More basic constructor where sweep_kwargs are constant throughout sweeps
+function SweepIterator(problem, nsweeps::Int; sweep_kwargs...)
+  # Initialize this to an empty RegionIterator
+  sweep_kwargs_iter = Iterators.repeated(sweep_kwargs, nsweeps)
+  return SweepIterator(problem, sweep_kwargs_iter)
 end
diff --git a/src/solvers/previous_interfaces/contract.jl b/src/solvers/previous_interfaces/contract.jl
deleted file mode 100644
index 00e5c4d6..00000000
--- a/src/solvers/previous_interfaces/contract.jl
+++ /dev/null
@@ -1,86 +0,0 @@
-using Graphs: nv, vertices
-using ITensors: ITensors, sim
-using ITensors.NDTensors: Algorithm, @Algorithm_str, contract
-using NamedGraphs: vertextype
-
-function sum_contract(
-  ::Algorithm"fit",
-  tns::Vector{<:Tuple{<:AbstractTTN,<:AbstractTTN}};
-  init,
-  nsites=2,
-  nsweeps=1,
-  cutoff=eps(),
-  updater=contract_updater,
-  kwargs...,
-)
-  tn1s = first.(tns)
-  tn2s = last.(tns)
-  ns = nv.(tn1s)
-  n = first(ns)
-  any(ns .!= nv.(tn2s)) && throw(
-    DimensionMismatch("Number of sites operator ($n) and state ($(nv(tn2))) do not match")
-  )
-  any(ns .!= n) &&
-    throw(DimensionMismatch("Number of sites in different operators ($n) do not match"))
-  # ToDo: Write test for single-vertex ttn, this implementation has not been tested.
-  if n == 1
-    res = 0
-    for (tn1, tn2) in zip(tn1s, tn2s)
-      v = only(vertices(tn2))
-      res += tn1[v] * tn2[v]
-    end
-    return typeof(tn2)([res])
-  end
-
-  # In case `tn1` and `tn2` have the same internal indices
-  operator = ProjOuterProdTTN{vertextype(first(tn1s))}[]
-  for (tn1, tn2) in zip(tn1s, tn2s)
-    tn1 = sim(linkinds, tn1)
-
-    # In case `init` and `tn2` have the same internal indices
-    init = sim(linkinds, init)
-    push!(operator, ProjOuterProdTTN(tn2, tn1))
-  end
-  operator = isone(length(operator)) ? only(operator) : ProjTTNSum(operator)
-  #ToDo: remove?
-  # Fix site and link inds of init
-  ## init = deepcopy(init)
-  ## init = sim(linkinds, init)
-  ## for v in vertices(tn2)
-  ##   replaceinds!(
-  ##     init[v], siteinds(init, v), uniqueinds(siteinds(tn1, v), siteinds(tn2, v))
-  ##   )
-  ## end
-
-  return alternating_update(operator, init; nsweeps, nsites, updater, cutoff, kwargs...)
-end
-
-function NDTensors.contract(
-  a::Algorithm"fit", tn1::AbstractTTN, tn2::AbstractTTN; kwargs...
-)
-  return sum_contract(a, [(tn1, tn2)]; kwargs...)
-end
-
-"""
-Overload of `ITensors.contract`.
-"""
-function NDTensors.contract(tn1::AbstractTTN, tn2::AbstractTTN; alg="fit", kwargs...)
-  return contract(Algorithm(alg), tn1, tn2; kwargs...)
-end
-
-"""
-Overload of `ITensors.apply`.
-"""
-function ITensors.apply(tn1::AbstractTTN, tn2::AbstractTTN; init, kwargs...)
-  init = init'
-  tn12 = contract(tn1, tn2; init, kwargs...)
-  return replaceprime(tn12, 1 => 0)
-end
-
-function sum_apply(
-  tns::Vector{<:Tuple{<:AbstractTTN,<:AbstractTTN}}; alg="fit", init, kwargs...
-)
-  init = init'
-  tn12 = sum_contract(Algorithm(alg), tns; init, kwargs...)
-  return replaceprime(tn12, 1 => 0)
-end
diff --git a/src/solvers/previous_interfaces/dmrg.jl b/src/solvers/previous_interfaces/dmrg.jl
deleted file mode 100644
index 1acbde35..00000000
--- a/src/solvers/previous_interfaces/dmrg.jl
+++ /dev/null
@@ -1,23 +0,0 @@
-using KrylovKit: KrylovKit
-
-function dmrg(
-  operator,
-  init_state;
-  nsweeps,
-  nsites=2,
-  updater=eigsolve_updater,
-  (region_observer!)=nothing,
-  kwargs...,
-)
-  eigvals_ref = Ref{Any}()
-  region_observer! = compose_observers(
-    region_observer!, ValuesObserver((; eigvals=eigvals_ref))
-  )
-  state = alternating_update(
-    operator, init_state; nsweeps, nsites, updater, region_observer!, kwargs...
-  )
-  eigval = only(eigvals_ref[])
-  return eigval, state
-end
-
-KrylovKit.eigsolve(H, init::AbstractTTN; kwargs...) = dmrg(H, init; kwargs...)
diff --git a/src/solvers/previous_interfaces/dmrg_x.jl b/src/solvers/previous_interfaces/dmrg_x.jl
deleted file mode 100644
index 7ab9d8cd..00000000
--- a/src/solvers/previous_interfaces/dmrg_x.jl
+++ /dev/null
@@ -1,19 +0,0 @@
-function dmrg_x(
-  operator,
-  init_state::AbstractTTN;
-  nsweeps,
-  nsites=2,
-  updater=dmrg_x_updater,
-  (region_observer!)=nothing,
-  kwargs...,
-)
-  eigvals_ref = Ref{Any}()
-  region_observer! = compose_observers(
-    region_observer!, ValuesObserver((; eigvals=eigvals_ref))
-  )
-  state = alternating_update(
-    operator, init_state; nsweeps, nsites, updater, region_observer!, kwargs...
-  )
-  eigval = only(eigvals_ref[])
-  return eigval, state
-end
diff --git a/src/solvers/previous_interfaces/linsolve.jl b/src/solvers/previous_interfaces/linsolve.jl
deleted file mode 100644
index acd93cef..00000000
--- a/src/solvers/previous_interfaces/linsolve.jl
+++ /dev/null
@@ -1,49 +0,0 @@
-using DocStringExtensions: TYPEDSIGNATURES
-using KrylovKit: KrylovKit
-
-"""
-$(TYPEDSIGNATURES)
-
-Compute a solution x to the linear system:
-
-(a₀ + a₁ * A)*x = b
-
-using starting guess x₀. Leaving a₀, a₁
-set to their default values solves the 
-system A*x = b.
-
-To adjust the balance between accuracy of solution
-and speed of the algorithm, it is recommed to first try
-adjusting the `solver_tol` keyword argument descibed below.
-
-Keyword arguments:
-  - `ishermitian::Bool=false` - should set to true if the MPO A is Hermitian
-  - `solver_krylovdim::Int=30` - max number of Krylov vectors to build on each solver iteration
-  - `solver_maxiter::Int=100` - max number outer iterations (restarts) to do in the solver step
-  - `solver_tol::Float64=1E-14` - tolerance or error goal of the solver
-
-Overload of `KrylovKit.linsolve`.
-"""
-function KrylovKit.linsolve(
-  A::AbstractTTN,
-  b::AbstractTTN,
-  x₀::AbstractTTN,
-  a₀::Number=0,
-  a₁::Number=1;
-  updater=linsolve_updater,
-  nsites=2,
-  nsweeps,  #it makes sense to require this to be defined
-  updater_kwargs=(;),
-  kwargs...,
-)
-  updater_kwargs = (; a₀, a₁, updater_kwargs...)
-  error("`linsolve` for TTN not yet implemented.")
-
-  # TODO: Define `itensornetwork_cache`
-  # TODO: Define `linsolve_cache`
-
-  P = linsolve_cache(itensornetwork_cache(x₀', A, x₀), itensornetwork_cache(x₀', b))
-  return alternating_update(
-    P, x₀; nsweeps, nsites, updater=linsolve_updater, updater_kwargs, kwargs...
-  )
-end
diff --git a/src/solvers/previous_interfaces/tdvp.jl b/src/solvers/previous_interfaces/tdvp.jl
deleted file mode 100644
index 7a58fe1b..00000000
--- a/src/solvers/previous_interfaces/tdvp.jl
+++ /dev/null
@@ -1,154 +0,0 @@
-using NamedGraphs.GraphsExtensions: GraphsExtensions
-
-#ToDo: Cleanup _compute_nsweeps, maybe restrict flexibility to simplify code
-function _compute_nsweeps(nsweeps::Int, t::Number, time_step::Number)
-  return error("Cannot specify both nsweeps and time_step in tdvp")
-end
-
-function _compute_nsweeps(nsweeps::Nothing, t::Number, time_step::Nothing)
-  return 1, [t]
-end
-
-function _compute_nsweeps(nsweeps::Nothing, t::Number, time_step::Number)
-  @assert isfinite(time_step) && abs(time_step) > 0.0
-  nsweeps = convert(Int, ceil(abs(t / time_step)))
-  if !(nsweeps * time_step ≈ t)
-    println("Time that will be reached = nsweeps * time_step = ", nsweeps * time_step)
-    println("Requested total time t = ", t)
-    error("Time step $time_step not commensurate with total time t=$t")
-  end
-  return nsweeps, extend_or_truncate(time_step, nsweeps)
-end
-
-function _compute_nsweeps(nsweeps::Int, t::Number, time_step::Nothing)
-  time_step = extend_or_truncate(t / nsweeps, nsweeps)
-  return nsweeps, time_step
-end
-
-function _compute_nsweeps(nsweeps, t::Number, time_step::Vector)
-  diff_time = t - sum(time_step)
-
-  isnothing(nsweeps)
-  if isnothing(nsweeps)
-    #extend_or_truncate time_step to reach final time t
-    last_time_step = last(time_step)
-    nsweepstopad = Int(ceil(abs(diff_time / last_time_step)))
-    if !(sum(time_step) + nsweepstopad * last_time_step ≈ t)
-      println(
-        "Time that will be reached = nsweeps * time_step = ",
-        sum(time_step) + nsweepstopad * last_time_step,
-      )
-      println("Requested total time t = ", t)
-      error("Time step $time_step not commensurate with total time t=$t")
-    end
-    time_step = extend_or_truncate(time_step, length(time_step) + nsweepstopad)
-    nsweeps = length(time_step)
-  else
-    nsweepstopad = nsweeps - length(time_step)
-    if abs(diff_time) < eps() && !iszero(nsweepstopad)
-      warn(
-        "A vector of timesteps that sums up to total time t=$t was supplied,
-   but its length (=$(length(time_step))) does not agree with supplied number of sweeps (=$(nsweeps)).",
-      )
-      return length(time_step), time_step
-    end
-    remaining_time_step = diff_time / nsweepstopad
-    append!(time_step, extend_or_truncate(remaining_time_step, nsweepstopad))
-  end
-  return nsweeps, time_step
-end
-
-function sub_time_steps(order)
-  if order == 1
-    return [1.0]
-  elseif order == 2
-    return [1 / 2, 1 / 2]
-  elseif order == 4
-    s = 1.0 / (2 - 2^(1 / 3))
-    return [s / 2, s / 2, (1 - 2 * s) / 2, (1 - 2 * s) / 2, s / 2, s / 2]
-  else
-    error("Trotter order of $order not supported")
-  end
-end
-
-"""
-    tdvp(operator::TTN, t::Number, init_state::TTN; kwargs...)
-
-Use the time dependent variational principle (TDVP) algorithm
-to approximately compute `exp(operator*t)*init_state` using an efficient algorithm based
-on alternating optimization of the state tensors and local Krylov
-exponentiation of operator. The time parameter `t` can be a real or complex number.
-                    
-Returns:
-* `state` - time-evolved state
-
-Optional keyword arguments:
-* `time_step::Number = t` - time step to use when evolving the state. Smaller time steps generally give more accurate results but can make the algorithm take more computational time to run.
-* `nsteps::Integer` - evolve by the requested total time `t` by performing `nsteps` of the TDVP algorithm. More steps can result in more accurate results but require more computational time to run. (Note that only one of the `time_step` or `nsteps` parameters can be provided, not both.)
-* `outputlevel::Int = 1` - larger outputlevel values resulting in printing more information and 0 means no output
-* `observer` - object implementing the Observer interface which can perform measurements and stop early
-* `write_when_maxdim_exceeds::Int` - when the allowed maxdim exceeds this value, begin saving tensors to disk to free memory in large calculations
-"""
-function tdvp(
-  operator,
-  t::Number,
-  init_state::AbstractTTN;
-  t_start=0.0,
-  time_step=nothing,
-  nsites=2,
-  nsweeps=nothing,
-  order::Integer=2,
-  outputlevel=default_outputlevel(),
-  region_printer=nothing,
-  sweep_printer=nothing,
-  (sweep_observer!)=nothing,
-  (region_observer!)=nothing,
-  root_vertex=GraphsExtensions.default_root_vertex(init_state),
-  reverse_step=true,
-  extracter_kwargs=(;),
-  extracter=default_extracter(), # ToDo: extracter could be inside extracter_kwargs, at the cost of having to extract it in region_update
-  updater_kwargs=(;),
-  updater=exponentiate_updater,
-  inserter_kwargs=(;),
-  inserter=default_inserter(),
-  transform_operator_kwargs=(;),
-  transform_operator=default_transform_operator(),
-  kwargs...,
-)
-  # move slurped kwargs into inserter
-  inserter_kwargs = (; inserter_kwargs..., kwargs...)
-  # process nsweeps and time_step
-  nsweeps, time_step = _compute_nsweeps(nsweeps, t, time_step)
-  t_evolved = t_start .+ cumsum(time_step)
-  sweep_plans = default_sweep_plans(
-    nsweeps,
-    init_state;
-    sweep_plan_func=tdvp_sweep_plan,
-    root_vertex,
-    reverse_step,
-    extracter,
-    extracter_kwargs,
-    updater,
-    updater_kwargs,
-    inserter,
-    inserter_kwargs,
-    transform_operator,
-    transform_operator_kwargs,
-    time_step,
-    order,
-    nsites,
-    t_evolved,
-  )
-
-  return alternating_update(
-    operator,
-    init_state,
-    sweep_plans;
-    outputlevel,
-    sweep_observer!,
-    region_observer!,
-    sweep_printer,
-    region_printer,
-  )
-  return state
-end
diff --git a/src/solvers/region_plans/dfs_plans.jl b/src/solvers/region_plans/dfs_plans.jl
index 074fa94a..9b44b980 100644
--- a/src/solvers/region_plans/dfs_plans.jl
+++ b/src/solvers/region_plans/dfs_plans.jl
@@ -3,14 +3,14 @@ using NamedGraphs.GraphsExtensions:
   default_root_vertex, post_order_dfs_edges, post_order_dfs_vertices
 
 function post_order_dfs_plan(
-  graph; nsites, root_vertex=default_root_vertex(graph), sweep_kwargs...
+  graph, sweep_kwargs; nsites, root_vertex=default_root_vertex(graph)
 )
   if nsites == 1
     vertices = post_order_dfs_vertices(graph, root_vertex)
-    fwd_sweep = [([v], sweep_kwargs) for v in vertices]
+    fwd_sweep = [[v] => sweep_kwargs for v in vertices]
   elseif nsites == 2
     edges = post_order_dfs_edges(graph, root_vertex)
-    fwd_sweep = [([src(e), dst(e)], sweep_kwargs) for e in edges]
+    fwd_sweep = [[src(e), dst(e)] => sweep_kwargs for e in edges]
   end
   return fwd_sweep
 end
diff --git a/src/solvers/region_plans/euler_plans.jl b/src/solvers/region_plans/euler_plans.jl
index cf661d0d..68548fdc 100644
--- a/src/solvers/region_plans/euler_plans.jl
+++ b/src/solvers/region_plans/euler_plans.jl
@@ -2,12 +2,14 @@ using Graphs: dst, src
 using NamedGraphs.GraphsExtensions: default_root_vertex
 
 function euler_sweep(graph; nsites, root_vertex=default_root_vertex(graph), sweep_kwargs...)
+  sweep_kwargs = (; nsites, root_vertex, sweep_kwargs...)
+
   if nsites == 1
     vertices = euler_tour_vertices(graph, root_vertex)
-    sweep = [([v], sweep_kwargs) for v in vertices]
+    sweep = [[v] => sweep_kwargs for v in vertices]
   elseif nsites == 2
     edges = euler_tour_edges(graph, root_vertex)
-    sweep = [([src(e), dst(e)], sweep_kwargs) for e in edges]
+    sweep = [[src(e), dst(e)] => sweep_kwargs for e in edges]
   end
   return sweep
 end
diff --git a/src/solvers/region_plans/tdvp_region_plans.jl b/src/solvers/region_plans/tdvp_region_plans.jl
index c03ad4eb..8ee90086 100644
--- a/src/solvers/region_plans/tdvp_region_plans.jl
+++ b/src/solvers/region_plans/tdvp_region_plans.jl
@@ -1,3 +1,5 @@
+using Accessors: @modify
+
 function applyexp_sub_steps(order)
   if order == 1
     return [1.0]
@@ -5,40 +7,54 @@ function applyexp_sub_steps(order)
     return [1 / 2, 1 / 2]
   elseif order == 4
     s = (2 - 2^(1 / 3))^(-1)
-    return [s/2, s/2, 1/2 - s, 1/2 - s, s/2, s/2]
+    return [s / 2, s / 2, 1 / 2 - s, 1 / 2 - s, s / 2, s / 2]
   else
     error("Applyexp order of $order not supported")
   end
 end
 
-function first_order_sweep(
-  graph, exponent_step, dir=Base.Forward; update_kwargs, nsites, kws...
-)
-  basic_fwd_sweep = post_order_dfs_plan(graph; nsites, kws...)
-  update_kwargs = (; nsites, exponent_step, update_kwargs...)
-  sweep = []
-  for (j, (region, region_kws)) in enumerate(basic_fwd_sweep)
-    push!(sweep, (region, (; nsites, update_kwargs, region_kws...)))
+function first_order_sweep(graph, sweep_kwargs; nsites)
+  basic_fwd_sweep = post_order_dfs_plan(graph, sweep_kwargs; nsites)
+  region_plan = []
+
+  for (j, (region, region_kwargs)) in enumerate(basic_fwd_sweep)
+    push!(region_plan, region => region_kwargs)
+
     if length(region) == 2 && j < length(basic_fwd_sweep)
-      rev_kwargs = (; update_kwargs..., exponent_step=(-update_kwargs.exponent_step))
-      push!(sweep, ([last(region)], (; update_kwargs=rev_kwargs, region_kws...)))
+      region_kwargs = @modify(-, region_kwargs.update!_kwargs.exponent_step)
+      push!(region_plan, [last(region)] => region_kwargs)
     end
   end
-  if dir==Base.Reverse
-    # Reverse regions as well as ordering of regions
-    sweep = [(reverse(reg_kws[1]), reg_kws[2]) for reg_kws in reverse(sweep)]
+
+  return region_plan
+end
+
+function reverse_regions(region_plan)
+  return map(reverse(region_plan)) do (region, kwargs)
+    return reverse(region) => kwargs
   end
-  return sweep
 end
 
-function applyexp_regions(graph, exponent_step; update_kwargs, order, nsites, kws...)
+# Generate the kwargs for each region.
+function applyexp_regions(
+  graph, raw_exponent_step; order, nsites, update!_kwargs=(; nsites), remaining_kwargs...
+)
   sweep_plan = []
+
   for (step, weight) in enumerate(applyexp_sub_steps(order))
-    dir = isodd(step) ? Base.Forward : Base.Reverse
-    append!(
-      sweep_plan,
-      first_order_sweep(graph, weight*exponent_step, dir; update_kwargs, nsites, kws...),
-    )
+    # Use this exponent step only if none provided
+    new_update!_kwargs = (; exponent_step=weight * raw_exponent_step, update!_kwargs...)
+
+    sweep_kwargs = (; remaining_kwargs..., update!_kwargs=new_update!_kwargs)
+
+    region_plan = first_order_sweep(graph, sweep_kwargs; nsites)
+
+    if iseven(step)
+      region_plan = reverse_regions(region_plan)
+    end
+
+    append!(sweep_plan, region_plan)
   end
+
   return sweep_plan
 end
diff --git a/src/solvers/subspace/densitymatrix.jl b/src/solvers/subspace/densitymatrix.jl
index 0c4cfa69..ae2ff507 100644
--- a/src/solvers/subspace/densitymatrix.jl
+++ b/src/solvers/subspace/densitymatrix.jl
@@ -1,70 +1,71 @@
 using NamedGraphs.GraphsExtensions: incident_edges
 using Printf: @printf
 
-function subspace_expand(
+function subspace_expand!(
   ::Backend"densitymatrix",
-  problem,
-  local_state::ITensor,
-  region_iterator;
-  expansion_factor,
-  max_expand,
+  region_iter,
+  local_state;
+  expansion_factor=1.5,
+  maxexpand=typemax(Int),
   north_pass=1,
-  trunc,
-  kws...,
+  eigen_kwargs=(;),
 )
-  region = current_region(region_iterator)
-  psi = copy(state(problem))
+  prob = problem(region_iter)
 
-  prev_vertex_set = setdiff(pos(operator(problem)), region)
-  (length(prev_vertex_set) != 1) && return problem, local_state
+  region = current_region(region_iter)
+  psi = copy(state(prob))
+
+  prev_vertex_set = setdiff(pos(operator(prob)), region)
+  (length(prev_vertex_set) != 1) && return region_iter, local_state
   prev_vertex = only(prev_vertex_set)
   A = psi[prev_vertex]
 
   next_vertices = filter(v -> (hascommoninds(psi[v], A)), region)
-  isempty(next_vertices) && return problem, local_state
+  isempty(next_vertices) && return region_iter, local_state
   next_vertex = only(next_vertices)
   C = psi[next_vertex]
 
   a = commonind(A, C)
-  isnothing(a) && return problem, local_state
+  isnothing(a) && return region_iter, local_state
   basis_size = prod(dim.(uniqueinds(A, C)))
 
   expanded_maxdim = compute_expansion(
-    dim(a), basis_size; expansion_factor, max_expand, trunc.maxdim
+    dim(a), basis_size; expansion_factor, maxexpand, eigen_kwargs.maxdim
   )
-  expanded_maxdim <= 0 && return problem, local_state
-  trunc = (; trunc..., maxdim=expanded_maxdim)
+  expanded_maxdim <= 0 && return region_iter, local_state
 
-  envs = environments(operator(problem))
-  H = operator(operator(problem))
+  envs = environments(operator(prob))
+  H = operator(operator(prob))
   sqrt_rho = A
-  for e in incident_edges(operator(problem))
+  for e in incident_edges(operator(prob))
     (src(e) ∈ region || dst(e) ∈ region) && continue
     sqrt_rho *= envs[e]
   end
   sqrt_rho *= H[prev_vertex]
 
-  conj_proj_A(T) = (T - prime(A)*(dag(prime(A))*T))
-  for pass in 1:north_pass
+  conj_proj_A(T) = (T - prime(A) * (dag(prime(A)) * T))
+  for _ in 1:north_pass
     sqrt_rho = conj_proj_A(sqrt_rho)
   end
   rho = sqrt_rho * dag(noprime(sqrt_rho))
-  D, U = eigen(rho; trunc..., ishermitian=true)
+  D, U = eigen(rho; eigen_kwargs..., ishermitian=true)
 
-  Uproj(T) = (T - prime(A, a)*(dag(prime(A, a))*T))
-  for pass in 1:north_pass
+  Uproj(T) = (T - prime(A, a) * (dag(prime(A, a)) * T))
+  for _ in 1:north_pass
     U = Uproj(U)
   end
-  if norm(dag(U)*A) > 1E-10
-    @printf("Warning: |U*A| = %.3E in subspace expansion\n", norm(dag(U)*A))
-    return problem, local_state
+  if norm(dag(U) * A) > 1E-10
+    @printf("Warning: |U*A| = %.3E in subspace expansion\n", norm(dag(U) * A))
+    return region_iter, local_state
   end
 
-  Ax, ax = directsum(A=>a, U=>commonind(U, D))
+  Ax, ax = directsum(A => a, U => commonind(U, D))
   expander = dag(Ax) * A
   psi[prev_vertex] = Ax
   psi[next_vertex] = expander * C
-  local_state = expander*local_state
+  local_state = expander * local_state
+
+  prob.state = psi
 
-  return set_state(problem, psi), local_state
+  return region_iter, local_state
 end
diff --git a/src/solvers/subspace/ortho_subspace.jl b/src/solvers/subspace/ortho_subspace.jl
index 7d7ca6c2..26465309 100644
--- a/src/solvers/subspace/ortho_subspace.jl
+++ b/src/solvers/subspace/ortho_subspace.jl
@@ -28,7 +28,7 @@ function subspace_expand!(
   max_expand=default_max_expand(),
   kws...,
 )
-  prev_region = previous_region(region_iterator)
+  prev_region = prev_region(region_iterator)
   region = current_region(region_iterator)
   if isnothing(prev_region) || isa(region, AbstractEdge)
     return local_tensor
diff --git a/src/solvers/subspace/subspace.jl b/src/solvers/subspace/subspace.jl
index 1c5dec87..d5388245 100644
--- a/src/solvers/subspace/subspace.jl
+++ b/src/solvers/subspace/subspace.jl
@@ -1,64 +1,41 @@
 using NDTensors: NDTensors
 using NDTensors.BackendSelection: Backend, @Backend_str
 
-default_expansion_factor() = 1.5
-default_max_expand() = typemax(Int)
+function subspace_expand!(region_iter, local_state; subspace_algorithm="nothing")
+  backend = Backend(subspace_algorithm)
 
-function subspace_expand(
-  problem,
-  local_state,
-  region_iterator;
-  expansion_factor=default_expansion_factor(),
-  max_expand=default_max_expand(),
-  subspace_algorithm=nothing,
-  sweep,
-  trunc,
-  kws...,
-)
-  expansion_factor = get_or_last(expansion_factor, sweep)
-  max_expand = get_or_last(max_expand, sweep)
-  return subspace_expand(
-    Backend(subspace_algorithm),
-    problem,
-    local_state,
-    region_iterator;
-    expansion_factor,
-    max_expand,
-    trunc,
-    kws...,
+  if backend isa Backend"nothing"
+    return region_iter, local_state
+  end
+
+  _, local_state = subspace_expand!(
+    backend, region_iter, local_state; region_kwargs(subspace_expand!, region_iter)...
   )
+
+  return region_iter, local_state
 end
 
-function subspace_expand(backend, problem, local_state, region_iterator; kws...)
-  error(
+function subspace_expand!(backend, region_iterator, local_state; kwargs...)
+  # We allow passing of any kwargs here is this method throws an error anyway
+  return error(
     "Subspace expansion (subspace_expand!) not defined for requested combination of subspace_algorithm and problem types",
   )
 end
 
-function subspace_expand(
-  backend::Backend{:nothing}, problem, local_state, region_iterator; kws...
-)
-  problem, local_state
-end
-
-function compute_expansion(
-  current_dim,
-  basis_size;
-  expansion_factor=default_expansion_factor(),
-  max_expand=default_max_expand(),
-  maxdim=default_maxdim(),
-)
+# Have these defaults set per backend in `subspace_expand!`
+function compute_expansion(current_dim, basis_size; expansion_factor, maxexpand, maxdim)
   # Note: expand_maxdim will be *added* to current bond dimension
   # Obtain expand_maxdim from expansion_factor
   expand_maxdim = ceil(Int, expansion_factor * current_dim)
   # Enforce max_expand keyword
-  expand_maxdim = min(max_expand, expand_maxdim)
+  expand_maxdim = min(maxexpand, expand_maxdim)
 
   # Restrict expand_maxdim below theoretical upper limit
-  expand_maxdim = min(basis_size-current_dim, expand_maxdim)
+  expand_maxdim = min(basis_size - current_dim, expand_maxdim)
   # Enforce total maxdim setting (e.g. used in insert step)
-  expand_maxdim = min(maxdim-current_dim, expand_maxdim)
+  expand_maxdim = min(maxdim - current_dim, expand_maxdim)
   # Ensure expand_maxdim is non-negative
   expand_maxdim = max(0, expand_maxdim)
+
   return expand_maxdim
 end
diff --git a/src/solvers/sweep_solve.jl b/src/solvers/sweep_solve.jl
index 3da97728..f67ce4de 100644
--- a/src/solvers/sweep_solve.jl
+++ b/src/solvers/sweep_solve.jl
@@ -1,40 +1,43 @@
 
-region_callback(problem; kws...) = nothing
+default_region_callback(sweep_iterator) = sweep_iterator
+default_sweep_callback(sweep_iterator) = sweep_iterator
 
-function sweep_callback(problem; outputlevel, sweep, nsweeps, kws...)
-  if outputlevel >= 1
-    println("Done with sweep $sweep/$nsweeps")
+# In this implementation the function `sweep_solve` is essentially just a wrapper around 
+# the iterate interface that allows one to pass callbacks.
+function sweep_solve!(
+  sweep_iterator;
+  sweep_callback=default_sweep_callback,
+  region_callback=default_region_callback,
+)
+  # Don't compute the region iteration automatically as we wish to insert a callback.
+  for _ in IncrementOnly(sweep_iterator)
+    for _ in region_iterator(sweep_iterator)
+      region_callback(sweep_iterator)
+    end
+    sweep_callback(sweep_iterator)
   end
+  return sweep_iterator
 end
 
-function sweep_solve(
-  sweep_iterator;
-  outputlevel=0,
-  region_callback=region_callback,
-  sweep_callback=sweep_callback,
-  kwargs...,
+# I suspect that `sweep_callback` is the more commonly used callback, so allow this to
+# be set using the `do` syntax.
+function sweep_solve!(
+  sweep_callback, sweep_iterator; region_callback=default_region_callback
 )
-  for (sweep, region_iter) in enumerate(sweep_iterator)
-    for (region, region_kwargs) in region_tuples(region_iter)
-      region_callback(
-        problem(region_iter);
-        nsweeps=length(sweep_iterator),
-        outputlevel,
-        region_iterator=region_iter,
-        region,
-        region_kwargs,
-        sweep,
-        kwargs...,
-      )
-    end
-    sweep_callback(
-      problem(region_iter);
-      nsweeps=length(sweep_iterator),
-      outputlevel,
-      region_iterator=region_iter,
-      sweep,
-      kwargs...,
-    )
+  return sweep_solve!(sweep_iterator; sweep_callback, region_callback)
+end
+
+function sweep_solve!(
+  each_region_iterator::EachRegion; region_callback=default_region_callback
+)
+  return sweep_solve!(region_callback, each_region_iterator)
+end
+function sweep_solve!(region_callback, each_region_iterator::EachRegion)
+  for _ in each_region_iterator
+    # I don't think it is obvious what object this particular callback should take,
+    # but for now be consistant and pass the parent sweep iterator.
+    sweep_iterator = each_region_iterator.parent
+    region_callback(sweep_iterator)
   end
-  return problem(sweep_iterator)
+  return each_region_iterator
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index 98b2d2b8..fb2673d0 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -24,7 +24,7 @@ end
 
 @time begin
   # tests in groups based on folder structure
-  for testgroup in filter(isdir, readdir(@__DIR__))
+  for testgroup in filter(f -> isdir(joinpath(@__DIR__, f)), readdir(@__DIR__))
     if GROUP == "ALL" || GROUP == uppercase(testgroup)
       groupdir = joinpath(@__DIR__, testgroup)
       for file in filter(istestfile, readdir(groupdir))
diff --git a/test/solvers/test_applyexp.jl b/test/solvers/test_applyexp.jl
index 69b83b7c..c3464c8d 100644
--- a/test/solvers/test_applyexp.jl
+++ b/test/solvers/test_applyexp.jl
@@ -12,110 +12,111 @@ function chain_plus_ancilla(; nchain)
   for j in 1:nchain
     add_vertex!(g, j)
   end
-  for j in 1:(nchain - 1)
-    add_edge!(g, j=>j+1)
+  for j in 1:(nchain-1)
+    add_edge!(g, j => j + 1)
   end
   # Add ancilla vertex near middle of chain
   add_vertex!(g, 0)
-  add_edge!(g, 0=>nchain÷2)
+  add_edge!(g, 0 => nchain ÷ 2)
   return g
 end
 
-@testset "Test Tree Time Evolution" begin
-  outputlevel = 0
+@testset "Time Evolution" begin
 
-  N = 10
-  g = chain_plus_ancilla(; nchain=N)
+  @testset "Test Tree Time Evolution" begin
+    outputlevel = 0
 
-  sites = siteinds("S=1/2", g)
+    N = 10
+    g = chain_plus_ancilla(; nchain=N)
 
-  # Make Heisenberg model Hamiltonian
-  h = OpSum()
-  for j in 1:(N - 1)
-    h += "Sz", j, "Sz", j+1
-    h += 1/2, "S+", j, "S-", j+1
-    h += 1/2, "S-", j, "S+", j+1
-  end
-  H = ttn(h, sites)
+    sites = siteinds("S=1/2", g)
 
-  # Make initial product state
-  state = Dict{Int,String}()
-  for (j, v) in enumerate(vertices(sites))
-    state[v] = iseven(j) ? "Up" : "Dn"
-  end
-  psi0 = ttn(state, sites)
+    # Make Heisenberg model Hamiltonian
+    h = OpSum()
+    for j in 1:(N-1)
+      h += "Sz", j, "Sz", j + 1
+      h += 1 / 2, "S+", j, "S-", j + 1
+      h += 1 / 2, "S-", j, "S+", j + 1
+    end
+    H = ttn(h, sites)
 
-  cutoff = 1E-10
-  maxdim = 100
-  nsweeps = 5
+    # Make initial product state
+    state = Dict{Int,String}()
+    for (j, v) in enumerate(vertices(sites))
+      state[v] = iseven(j) ? "Up" : "Dn"
+    end
+    psi0 = ttn(state, sites)
 
-  nsites = 2
-  trunc = (; cutoff, maxdim)
-  E, gs_psi = dmrg(H, psi0; insert_kwargs=(; trunc), nsites, nsweeps, outputlevel)
-  (outputlevel >= 1) && println("2-site DMRG energy = ", E)
+    cutoff = 1E-10
+    maxdim = 100
+    nsweeps = 5
 
-  insert_kwargs=(; trunc)
-  nsites = 1
-  tmax = 0.10
-  time_range = 0.0:0.02:tmax
-  psi1_t = time_evolve(H, time_range, gs_psi; insert_kwargs, nsites, outputlevel)
-  (outputlevel >= 1) && println("Done with $nsites-site TDVP")
+    nsites = 2
+    factorize_kwargs = (; cutoff, maxdim)
+    E, gs_psi = dmrg(H, psi0; factorize_kwargs, nsites, nsweeps, outputlevel)
+    (outputlevel >= 1) && println("2-site DMRG energy = ", E)
 
-  @test norm(psi1_t) > 0.999
+    nsites = 1
+    tmax = 0.10
+    time_range = 0.0:0.02:tmax
+    psi1_t = time_evolve(H, time_range, gs_psi; factorize_kwargs, nsites, outputlevel)
+    (outputlevel >= 1) && println("Done with $nsites-site TDVP")
 
-  nsites = 2
-  psi2_t = time_evolve(H, time_range, gs_psi; insert_kwargs, nsites, outputlevel)
-  (outputlevel >= 1) && println("Done with $nsites-site TDVP")
-  @test norm(psi2_t) > 0.999
+    @test norm(psi1_t) > 0.999
 
-  @test abs(inner(psi1_t, gs_psi)) > 0.99
-  @test abs(inner(psi1_t, psi2_t)) > 0.99
+    nsites = 2
+    psi2_t = time_evolve(H, time_range, gs_psi; factorize_kwargs, nsites, outputlevel)
+    (outputlevel >= 1) && println("Done with $nsites-site TDVP")
+    @test norm(psi2_t) > 0.999
 
-  # Test that accumulated phase angle is E*tmax
-  z = inner(psi1_t, gs_psi)
-  @test abs(atan(imag(z)/real(z)) - E*tmax) < 1E-4
-end
+    @test abs(inner(psi1_t, gs_psi)) > 0.99
+    @test abs(inner(psi1_t, psi2_t)) > 0.99
 
-@testset "Applyexp Time Point Handling" begin
-  N = 10
-  g = named_path_graph(N)
-  sites = siteinds("S=1/2", g)
-
-  # Make Heisenberg model Hamiltonian
-  h = OpSum()
-  for j in 1:(N - 1)
-    h += "Sz", j, "Sz", j+1
-    h += 1/2, "S+", j, "S-", j+1
-    h += 1/2, "S-", j, "S+", j+1
+    # Test that accumulated phase angle is E*tmax
+    z = inner(psi1_t, gs_psi)
+    @test atan(imag(z) / real(z)) ≈ E * tmax atol = 1E-4
   end
-  H = ttn(h, sites)
 
-  # Initial product state
-  state = Dict{Int,String}()
-  for (j, v) in enumerate(vertices(sites))
-    state[v] = iseven(j) ? "Up" : "Dn"
-  end
-  psi0 = ttn(state, sites)
-
-  nsites = 2
-  trunc = (; cutoff=1E-8, maxdim=100)
-  insert_kwargs=(; trunc)
-
-  # Test that all time points are reached and reported correctly
-  time_points = [0.0,0.1,0.25,0.32,0.4]
-  times = Real[]
-  function collect_times(problem; kws...)
-    push!(times, ITensorNetworks.current_time(problem))
-  end
-  time_evolve(H, time_points, psi0; insert_kwargs, nsites, sweep_callback=collect_times,outputlevel=1)
-  @test norm(times - time_points) < 10*eps(Float64)
-
-  # Test that all exponents are reached and reported correctly
-  exponent_points = [-0.0,-0.1,-0.25,-0.32,-0.4]
-  exponents = Real[]
-  function collect_exponents(problem; kws...)
-    push!(exponents, ITensorNetworks.current_exponent(problem))
+  @testset "Applyexp Time Point Handling" begin
+    N = 10
+    g = named_path_graph(N)
+    sites = siteinds("S=1/2", g)
+
+    # Make Heisenberg model Hamiltonian
+    h = OpSum()
+    for j in 1:(N-1)
+      h += "Sz", j, "Sz", j + 1
+      h += 1 / 2, "S+", j, "S-", j + 1
+      h += 1 / 2, "S-", j, "S+", j + 1
+    end
+    H = ttn(h, sites)
+
+    # Initial product state
+    state = Dict{Int,String}()
+    for (j, v) in enumerate(vertices(sites))
+      state[v] = iseven(j) ? "Up" : "Dn"
+    end
+    psi0 = ttn(state, sites)
+
+    nsites = 2
+    factorize_kwargs = (; cutoff=1E-8, maxdim=100)
+
+    # Test that all time points are reached and reported correctly
+    time_points = [0.0, 0.1, 0.25, 0.32, 0.4]
+    times = Real[]
+    function collect_times(sweep_iterator; kws...)
+      push!(times, ITensorNetworks.current_time(ITensorNetworks.problem(sweep_iterator)))
+    end
+    time_evolve(H, time_points, psi0; factorize_kwargs, nsites, sweep_callback=collect_times, outputlevel=1)
+    @test times ≈ time_points atol = 10 * eps(Float64)
+
+    # Test that all exponents are reached and reported correctly
+    exponent_points = [-0.0, -0.1, -0.25, -0.32, -0.4]
+    exponents = Real[]
+    function collect_exponents(sweep_iterator; kws...)
+      push!(exponents, ITensorNetworks.current_exponent(ITensorNetworks.problem(sweep_iterator)))
+    end
+    applyexp(H, exponent_points, psi0; factorize_kwargs, nsites, sweep_callback=collect_exponents, outputlevel=1)
+    @test exponents ≈ exponent_points atol = 10 * eps(Float64)
   end
-  applyexp(H, exponent_points, psi0; insert_kwargs, nsites, sweep_callback=collect_exponents,outputlevel=1)
-  @test norm(exponents - exponent_points) < 10*eps(Float64)
 end
diff --git a/test/solvers/test_eigsolve.jl b/test/solvers/test_eigsolve.jl
index 75194cae..5a29c2c4 100644
--- a/test/solvers/test_eigsolve.jl
+++ b/test/solvers/test_eigsolve.jl
@@ -20,8 +20,8 @@ include("utilities/tree_graphs.jl")
   for edge in edges(sites)
     i, j = src(edge), dst(edge)
     h += "Sz", i, "Sz", j
-    h += 1/2, "S+", i, "S-", j
-    h += 1/2, "S-", i, "S+", j
+    h += 1 / 2, "S+", i, "S-", j
+    h += 1 / 2, "S-", i, "S+", j
   end
   H = ttn(h, sites)
 
@@ -38,27 +38,26 @@ include("utilities/tree_graphs.jl")
 
   cutoff = 1E-5
   maxdim = 40
+
+  factorize_kwargs = (; cutoff, maxdim)
+
   nsweeps = 5
 
   #
   # Test 2-site DMRG without subspace expansion
   #
   nsites = 2
-  trunc = (; cutoff, maxdim)
-  insert_kwargs = (; trunc)
-  E, psi = dmrg(H, psi0; insert_kwargs, nsites, nsweeps, outputlevel)
+  E, psi = dmrg(H, psi0; factorize_kwargs, nsites, nsweeps, outputlevel)
   (outputlevel >= 1) && println("2-site DMRG energy = ", E)
-  @test abs(E-Ex) < 1E-5
+  @test E ≈ Ex atol = 1E-5
 
   #
   # Test 1-site DMRG with subspace expansion
   #
   nsites = 1
   nsweeps = 5
-  trunc = (; cutoff, maxdim)
-  extract_kwargs = (; trunc, subspace_algorithm="densitymatrix")
-  insert_kwargs = (; trunc)
-  E, psi = dmrg(H, psi0; extract_kwargs, insert_kwargs, nsites, nsweeps, outputlevel)
+  extract!_kwargs = (; subspace_algorithm="densitymatrix")
+  E, psi = dmrg(H, psi0; extract!_kwargs, factorize_kwargs, nsites, nsweeps, outputlevel)
   (outputlevel >= 1) && println("1-site+subspace DMRG energy = ", E)
-  @test abs(E-Ex) < 1E-5
+  @test E ≈ Ex atol = 1E-5
 end
diff --git a/test/solvers/test_iterators.jl b/test/solvers/test_iterators.jl
new file mode 100644
index 00000000..438f067e
--- /dev/null
+++ b/test/solvers/test_iterators.jl
@@ -0,0 +1,176 @@
+using Test: @test, @testset
+using ITensorNetworks: SweepIterator, islaststep, state, increment!, compute!, eachregion
+
+module TestIteratorUtils
+
+using ITensorNetworks
+
+struct TestProblem <: ITensorNetworks.AbstractProblem
+  data::Vector{Int}
+end
+ITensorNetworks.region_plan(::TestProblem) = [:a => (; val=1), :b => (; val=2)]
+function ITensorNetworks.compute!(iter::ITensorNetworks.RegionIterator{<:TestProblem})
+  kwargs = ITensorNetworks.region_kwargs(iter)
+  push!(ITensorNetworks.problem(iter).data, kwargs.val)
+  return iter
+end
+
+
+mutable struct TestIterator <: ITensorNetworks.AbstractNetworkIterator
+  state::Int
+  max::Int
+  output::Vector{Int}
+end
+
+ITensorNetworks.increment!(TI::TestIterator) = TI.state += 1
+Base.length(TI::TestIterator) = TI.max
+ITensorNetworks.state(TI::TestIterator) = TI.state
+function ITensorNetworks.compute!(TI::TestIterator)
+  push!(TI.output, ITensorNetworks.state(TI))
+  return TI
+end
+
+mutable struct SquareAdapter <: ITensorNetworks.AbstractNetworkIterator
+  parent::TestIterator
+end
+
+Base.length(SA::SquareAdapter) = length(SA.parent)
+ITensorNetworks.increment!(SA::SquareAdapter) = ITensorNetworks.increment!(SA.parent)
+ITensorNetworks.state(SA::SquareAdapter) = ITensorNetworks.state(SA.parent)
+function ITensorNetworks.compute!(SA::SquareAdapter)
+  ITensorNetworks.compute!(SA.parent)
+  return last(SA.parent.output)^2
+end
+
+end
+
+@testset "Iterators" begin
+
+  import .TestIteratorUtils
+
+  @testset "`AbstractNetworkIterator` Interface" begin
+    TI = TestIteratorUtils.TestIterator(1, 4, [])
+
+    @test !islaststep((TI))
+
+    # First iterator should compute only
+    rv, st = iterate(TI)
+    @test !islaststep((TI))
+    @test !st
+    @test rv === TI
+    @test length(TI.output) == 1
+    @test only(TI.output) == 1
+    @test state(TI) == 1
+    @test !st
+
+    rv, st = iterate(TI, st)
+    @test !islaststep((TI))
+    @test !st
+    @test length(TI.output) == 2
+    @test state(TI) == 2
+    @test TI.output == [1, 2]
+
+    increment!(TI)
+    @test !islaststep((TI))
+    @test state(TI) == 3
+    @test length(TI.output) == 2
+    @test TI.output == [1, 2]
+
+    compute!(TI)
+    @test !islaststep((TI))
+    @test state(TI) == 3
+    @test length(TI.output) == 3
+    @test TI.output == [1, 2, 3]
+
+    # Final Step
+    iterate(TI, false)
+    @test islaststep((TI))
+    @test state(TI) == 4
+    @test length(TI.output) == 4
+    @test TI.output == [1, 2, 3, 4]
+
+    @test iterate(TI, false) === nothing
+
+    TI = TestIteratorUtils.TestIterator(1, 5, [])
+
+    cb = []
+
+    for _ in TI
+      @test length(cb) == length(TI.output) - 1
+      @test cb == (TI.output)[1:end-1]
+      push!(cb, state(TI))
+      @test cb == TI.output
+    end
+
+    @test islaststep((TI))
+    @test length(TI.output) == 5
+    @test length(cb) == 5
+    @test cb == TI.output
+
+
+    TI = TestIteratorUtils.TestIterator(1, 5, [])
+  end
+
+  @testset "Adapters" begin
+    TI = TestIteratorUtils.TestIterator(1, 5, [])
+    SA = TestIteratorUtils.SquareAdapter(TI)
+
+    @testset "Generic" begin
+
+      i = 0
+      for rv in SA
+        i += 1
+        @test rv isa Int
+        @test rv == i^2
+        @test state(SA) == i
+      end
+
+      @test islaststep((SA))
+
+      TI = TestIteratorUtils.TestIterator(1, 5, [])
+      SA = TestIteratorUtils.SquareAdapter(TI)
+
+      SA_c = collect(SA)
+
+      @test SA_c isa Vector
+      @test length(SA_c) == 5
+      @test SA_c == [1, 4, 9, 16, 25]
+
+    end
+
+    @testset "EachRegion" begin
+      prob = TestIteratorUtils.TestProblem([])
+      prob_region = TestIteratorUtils.TestProblem([])
+
+      SI = SweepIterator(prob, 5)
+      SI_region = SweepIterator(prob_region, 5)
+
+      callback = []
+      callback_region = []
+
+      let i = 1
+        for _ in SI
+          push!(callback, i)
+          i += 1
+        end
+      end
+
+      @test length(callback) == 5
+
+      let i = 1
+        for _ in eachregion(SI_region)
+          push!(callback_region, i)
+          i += 1
+        end
+      end
+
+      @test length(callback_region) == 10
+
+      @test prob.data == prob_region.data
+
+      @test prob.data[1:2:end] == fill(1, 5)
+      @test prob.data[2:2:end] == fill(2, 5)
+
+    end
+  end
+end
diff --git a/test/test_ttn_contract.jl b/test/test_ttn_contract.jl
deleted file mode 100644
index 500826c3..00000000
--- a/test/test_ttn_contract.jl
+++ /dev/null
@@ -1,152 +0,0 @@
-@eval module $(gensym())
-using Graphs: vertices
-using ITensorNetworks:
-  ITensorNetworks,
-  OpSum,
-  ProjOuterProdTTN,
-  ProjTTNSum,
-  ttn,
-  apply,
-  contract,
-  delta,
-  dmrg,
-  inner,
-  mpo,
-  random_mps,
-  random_ttn,
-  siteinds
-using ITensorNetworks.ModelHamiltonians: ModelHamiltonians
-using ITensors: prime, replaceinds, replaceprime
-using LinearAlgebra: norm, normalize
-using NamedGraphs.NamedGraphGenerators: named_comb_tree
-using StableRNGs: StableRNG
-using Test: @test, @test_broken, @testset
-
-@testset "Contract MPO" begin
-  N = 20
-  s = siteinds("S=1/2", N)
-  rng = StableRNG(1234)
-  psi = random_mps(rng, s; link_space=8)
-  os = OpSum()
-  for j in 1:(N - 1)
-    os += 0.5, "S+", j, "S-", j + 1
-    os += 0.5, "S-", j, "S+", j + 1
-    os += "Sz", j, "Sz", j + 1
-  end
-  for j in 1:(N - 2)
-    os += 0.5, "S+", j, "S-", j + 2
-    os += 0.5, "S-", j, "S+", j + 2
-    os += "Sz", j, "Sz", j + 2
-  end
-  H = mpo(os, s)
-
-  # Test basic usage with default parameters
-  Hpsi = apply(H, psi; alg="fit", init=psi, nsweeps=1)
-  @test inner(psi, Hpsi) ≈ inner(psi', H, psi) rtol = 1e-5
-  # Test variational compression via DMRG
-  Hfit = ProjOuterProdTTN(psi', H)
-  e, Hpsi_via_dmrg = dmrg(Hfit, psi; updater_kwargs=(; which_eigval=:LR,), nsweeps=1)
-  @test abs(inner(Hpsi_via_dmrg, Hpsi / norm(Hpsi))) ≈ 1 rtol = 1e-4
-  # Test whether the interface works for ProjTTNSum with factors
-  Hfit = ProjTTNSum([ProjOuterProdTTN(psi', H), ProjOuterProdTTN(psi', H)], [-0.2, -0.8])
-  e, Hpsi_via_dmrg = dmrg(Hfit, psi; nsweeps=1, updater_kwargs=(; which_eigval=:SR,))
-  @test abs(inner(Hpsi_via_dmrg, Hpsi / norm(Hpsi))) ≈ 1 rtol = 1e-4
-
-  # Test basic usage for use with multiple ProjOuterProdTTN with default parameters
-  # BLAS.axpy-like test
-  os_id = OpSum()
-  os_id += -1, "Id", 1, "Id", 2
-  minus_identity = mpo(os_id, s)
-  os_id = OpSum()
-  os_id += +1, "Id", 1, "Id", 2
-  identity = mpo(os_id, s)
-  Hpsi = ITensorNetworks.sum_apply(
-    [(H, psi), (minus_identity, psi)]; alg="fit", init=psi, nsweeps=3
-  )
-  @test inner(psi, Hpsi) ≈ (inner(psi', H, psi) - norm(psi)^2) rtol = 1e-5
-  # Test the above via DMRG
-  # ToDo: Investigate why this is broken
-  Hfit = ProjTTNSum([ProjOuterProdTTN(psi', H), ProjOuterProdTTN(psi', identity)], [-1, 1])
-  e, Hpsi_normalized = dmrg(Hfit, psi; nsweeps=3, updater_kwargs=(; which_eigval=:SR))
-  @test_broken abs(inner(Hpsi, (Hpsi_normalized) / norm(Hpsi))) ≈ 1 rtol = 1e-5
-
-  #
-  # Change "top" indices of MPO to be a different set
-  #
-  t = siteinds("S=1/2", N)
-  psit = deepcopy(psi)
-
-  for j in 1:N
-    H[j] *= delta(s[j]', t[j])
-    psit[j] *= delta(s[j], t[j])
-  end
-  # Test with nsweeps=3
-  Hpsi = contract(H, psi; alg="fit", init=psit, nsweeps=3)
-  @test inner(psit, Hpsi) ≈ inner(psit, H, psi) rtol = 1e-5
-  # Test with less good initial guess MPS not equal to psi
-  psi_guess = truncate(psit; maxdim=2)
-  Hpsi = contract(H, psi; alg="fit", nsweeps=4, init=psi_guess)
-  @test inner(psit, Hpsi) ≈ inner(psit, H, psi) rtol = 1e-5
-
-  # Test with nsite=1
-  rng = StableRNG(1234)
-  Hpsi_guess = random_mps(rng, t; link_space=32)
-  Hpsi = contract(H, psi; alg="fit", init=Hpsi_guess, nsites=1, nsweeps=4)
-  @test inner(psit, Hpsi) ≈ inner(psit, H, psi) rtol = 1e-4
-end
-
-@testset "Contract TTN" begin
-  tooth_lengths = fill(4, 4)
-  root_vertex = (1, 4)
-  c = named_comb_tree(tooth_lengths)
-
-  s = siteinds("S=1/2", c)
-  rng = StableRNG(1234)
-  psi = normalize(random_ttn(rng, s; link_space=8))
-
-  os = ModelHamiltonians.heisenberg(c; J1=1, J2=1)
-  H = ttn(os, s)
-
-  # Test basic usage with default parameters
-  Hpsi = apply(H, psi; alg="fit", init=psi, nsweeps=1, cutoff=eps())
-  @test inner(psi, Hpsi) ≈ inner(psi', H, psi) rtol = 1e-5
-  # Test usage with non-default parameters
-  Hpsi = apply(
-    H, psi; alg="fit", init=psi, nsweeps=5, maxdim=[16, 32], cutoff=[1e-4, 1e-8, 1e-12]
-  )
-  @test inner(psi, Hpsi) ≈ inner(psi', H, psi) rtol = 1e-2
-
-  # Test basic usage for multiple ProjOuterProdTTN with default parameters
-  # BLAS.axpy-like test
-  os_id = OpSum()
-  os_id += -1, "Id", first(vertices(s)), "Id", first(vertices(s))
-  minus_identity = ttn(os_id, s)
-  Hpsi = ITensorNetworks.sum_apply(
-    [(H, psi), (minus_identity, psi)]; alg="fit", init=psi, nsweeps=1
-  )
-  @test inner(psi, Hpsi) ≈ (inner(psi', H, psi) - norm(psi)^2) rtol = 1e-5
-
-  #
-  # Change "top" indices of TTN to be a different set
-  #
-  t = siteinds("S=1/2", c)
-  psit = deepcopy(psi)
-  psit = replaceinds(psit, s => t)
-  H = replaceinds(H, prime(s; links=[]) => t)
-
-  # Test with nsweeps=2
-  Hpsi = contract(H, psi; alg="fit", init=psit, nsweeps=2)
-  @test inner(psit, Hpsi) ≈ inner(psit, H, psi) rtol = 1e-5
-
-  # Test with less good initial guess MPS not equal to psi
-  Hpsi_guess = truncate(psit; maxdim=2)
-  Hpsi = contract(H, psi; alg="fit", nsweeps=4, init=Hpsi_guess)
-  @test inner(psit, Hpsi) ≈ inner(psit, H, psi) rtol = 1e-5
-
-  # Test with nsite=1
-  rng = StableRNG(1234)
-  Hpsi_guess = random_ttn(rng, t; link_space=32)
-  Hpsi = contract(H, psi; alg="fit", nsites=1, nsweeps=10, init=Hpsi_guess)
-  @test inner(psit, Hpsi) ≈ inner(psit, H, psi) rtol = 1e-2
-end
-end
diff --git a/test/test_ttn_dmrg.jl b/test/test_ttn_dmrg.jl
deleted file mode 100644
index b8a8cdb8..00000000
--- a/test/test_ttn_dmrg.jl
+++ /dev/null
@@ -1,328 +0,0 @@
-@eval module $(gensym())
-using DataGraphs: edge_data, vertex_data
-using Dictionaries: Dictionary
-using Graphs: nv, vertices, uniform_tree
-using ITensorMPS: ITensorMPS
-using ITensorNetworks:
-  ITensorNetworks,
-  OpSum,
-  ttn,
-  apply,
-  dmrg,
-  inner,
-  mpo,
-  random_mps,
-  random_ttn,
-  linkdims,
-  siteinds
-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
-# if there are no QNs and auto-fermion
-# is enabled.
-ITensors.disable_auto_fermion()
-
-@testset "MPS DMRG" for nsites in [1, 2]
-  N = 10
-  cutoff = 1e-12
-
-  s = siteinds("S=1/2", N)
-
-  os = OpSum()
-  for j in 1:(N - 1)
-    os += 0.5, "S+", j, "S-", j + 1
-    os += 0.5, "S-", j, "S+", j + 1
-    os += "Sz", j, "Sz", j + 1
-  end
-
-  H = mpo(os, s)
-
-  rng = StableRNG(1234)
-  psi = random_mps(rng, s; link_space=20)
-
-  nsweeps = 10
-  maxdim = [10, 20, 40, 100]
-
-  # Compare to `ITensors.MPO` version of `dmrg`
-  H_mpo = ITensorMPS.MPO([H[v] for v in 1:nv(H)])
-  psi_mps = ITensorMPS.MPS([psi[v] for v in 1:nv(psi)])
-  e2, psi2 = ITensorMPS.dmrg(H_mpo, psi_mps; nsweeps, maxdim, outputlevel=0)
-
-  e, psi = dmrg(
-    H, psi; nsweeps, maxdim, cutoff, nsites, updater_kwargs=(; krylovdim=3, maxiter=1)
-  )
-  @test inner(psi', H, psi) ≈ e
-  @test inner(psi', H, psi) ≈ inner(psi2', H_mpo, psi2)
-
-  # Alias for `ITensorNetworks.dmrg`
-  e, psi = eigsolve(
-    H, psi; nsweeps, maxdim, cutoff, nsites, updater_kwargs=(; krylovdim=3, maxiter=1)
-  )
-  @test inner(psi', H, psi) ≈ e
-  @test inner(psi', H, psi) ≈ inner(psi2', H_mpo, psi2)
-
-  # Test custom sweep regions #BROKEN, ToDo: Make proper custom sweep regions for test
-  #=
-  orig_E = inner(psi', H, psi)
-  sweep_regions = [[1], [2], [3], [3], [2], [1]]
-  e, psi = dmrg(H, psi; nsweeps, maxdim, cutoff, sweep_regions)
-  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
-  N = 10
-  cutoff = 1e-12
-  s = siteinds("S=1/2", N)
-  os = OpSum()
-  for j in 1:(N - 1)
-    os += 0.5, "S+", j, "S-", j + 1
-    os += 0.5, "S-", j, "S+", j + 1
-    os += "Sz", j, "Sz", j + 1
-  end
-  H = mpo(os, s)
-  rng = StableRNG(1234)
-  psi = random_mps(rng, s; link_space=20)
-
-  nsweeps = 4
-  maxdim = [20, 40, 80, 80]
-  cutoff = [1e-10]
-
-  #
-  # Make observers
-  #
-  sweep(; which_sweep, kw...) = which_sweep
-  sweep_observer! = observer(sweep)
-
-  region(; which_region_update, sweep_plan, kw...) = first(sweep_plan[which_region_update])
-  energy(; eigvals, kw...) = eigvals[1]
-  region_observer! = observer(region, sweep, energy)
-
-  e, psi = dmrg(H, psi; nsweeps, maxdim, cutoff, sweep_observer!, region_observer!)
-
-  #
-  # Test out certain values
-  #
-  @test region_observer![9, :region] == [2, 1]
-  @test region_observer![30, :energy] < -4.25
-  @test region_observer![30, :energy] ≈ e rtol = 1e-6
-end
-
-@testset "Cache to Disk" begin
-  N = 10
-  cutoff = 1e-12
-  s = siteinds("S=1/2", N)
-  os = OpSum()
-  for j in 1:(N - 1)
-    os += 0.5, "S+", j, "S-", j + 1
-    os += 0.5, "S-", j, "S+", j + 1
-    os += "Sz", j, "Sz", j + 1
-  end
-  H = mpo(os, s)
-  rng = StableRNG(1234)
-  psi = random_mps(rng, s; link_space=10)
-
-  nsweeps = 4
-  maxdim = [10, 20, 40, 80]
-
-  @test_broken e, psi = dmrg(
-    H,
-    psi;
-    nsweeps,
-    maxdim,
-    cutoff,
-    outputlevel=0,
-    transform_operator=ITensorNetworks.cache_operator_to_disk,
-    transform_operator_kwargs=(; write_when_maxdim_exceeds=11),
-  )
-end
-
-@testset "Regression test: Arrays of Parameters" begin
-  N = 10
-  cutoff = 1e-12
-
-  s = siteinds("S=1/2", N)
-
-  os = OpSum()
-  for j in 1:(N - 1)
-    os += 0.5, "S+", j, "S-", j + 1
-    os += 0.5, "S-", j, "S+", j + 1
-    os += "Sz", j, "Sz", j + 1
-  end
-
-  H = mpo(os, s)
-
-  rng = StableRNG(1234)
-  psi = random_mps(rng, s; link_space=20)
-
-  # Choose nsweeps to be less than length of arrays
-  nsweeps = 5
-  maxdim = [200, 250, 400, 600, 800, 1200, 2000, 2400, 2600, 3000]
-  cutoff = [1e-10, 1e-10, 1e-12, 1e-12, 1e-12, 1e-12, 1e-14, 1e-14, 1e-14, 1e-14]
-
-  e, psi = dmrg(H, psi; nsweeps, maxdim, cutoff)
-end
-
-@testset "Tree DMRG" for nsites in [2]
-  cutoff = 1e-12
-
-  tooth_lengths = fill(2, 3)
-  c = named_comb_tree(tooth_lengths)
-
-  @testset "SVD approach" for use_qns in [false, true]
-    auto_fermion_enabled = ITensors.using_auto_fermion()
-    if use_qns  # test whether autofermion breaks things when using non-fermionic QNs
-      ITensors.enable_auto_fermion()
-    else        # when using no QNs, autofermion breaks # ToDo reference Issue in ITensors
-      ITensors.disable_auto_fermion()
-    end
-    s = siteinds("S=1/2", c; conserve_qns=use_qns)
-
-    os = ModelHamiltonians.heisenberg(c)
-
-    H = ttn(os, s)
-
-    # make init_state
-    d = Dict()
-    for (i, v) in enumerate(vertices(s))
-      d[v] = isodd(i) ? "Up" : "Dn"
-    end
-    states = v -> d[v]
-    psi = ttn(states, s)
-
-    # rng = StableRNG(1234)
-    # psi = random_ttn(rng, s; link_space=20) #FIXME: random_ttn broken for QN conserving case
-
-    nsweeps = 10
-    maxdim = [10, 20, 40, 100]
-    @show use_qns
-    e, psi = dmrg(
-      H, psi; nsweeps, maxdim, cutoff, nsites, updater_kwargs=(; krylovdim=3, maxiter=1)
-    )
-
-    # Compare to `ITensors.MPO` version of `dmrg`
-    linear_order = [4, 1, 2, 5, 3, 6]
-    vmap = Dictionary(collect(vertices(s))[linear_order], 1:length(linear_order))
-    sline = only.(collect(vertex_data(s)))[linear_order]
-    Hline = ITensorMPS.MPO(replace_vertices(v -> vmap[v], os), sline)
-    rng = StableRNG(1234)
-    psiline = ITensorMPS.random_mps(rng, sline, i -> isodd(i) ? "Up" : "Dn"; linkdims=20)
-    e2, psi2 = ITensorMPS.dmrg(Hline, psiline; nsweeps, maxdim, cutoff, outputlevel=0)
-
-    @test inner(psi', H, psi) ≈ ITensorMPS.inner(psi2', Hline, psi2) atol = 1e-5
-
-    if !auto_fermion_enabled
-      ITensors.disable_auto_fermion()
-    end
-  end
-end
-
-@testset "Tree DMRG for Fermions" for nsites in [2] #ToDo: change to [1,2] when random_ttn works with QNs
-  auto_fermion_enabled = ITensors.using_auto_fermion()
-  use_qns = true
-  cutoff = 1e-12
-  nsweeps = 10
-  maxdim = [10, 20, 40, 100]
-
-  # setup model
-  tooth_lengths = fill(2, 3)
-  c = named_comb_tree(tooth_lengths)
-  s = siteinds("Electron", c; conserve_qns=use_qns)
-  U = 2.0
-  t = 1.3
-  tp = 0.6
-  os = ModelHamiltonians.hubbard(c; U, t, tp)
-
-  # for conversion to ITensors.MPO
-  linear_order = [4, 1, 2, 5, 3, 6]
-  vmap = Dictionary(collect(vertices(s))[linear_order], 1:length(linear_order))
-  sline = only.(collect(vertex_data(s)))[linear_order]
-
-  # get MPS / MPO with JW string result
-  ITensors.disable_auto_fermion()
-  Hline = ITensorMPS.MPO(replace_vertices(v -> vmap[v], os), sline)
-  rng = StableRNG(1234)
-  psiline = ITensorMPS.random_mps(rng, sline, i -> isodd(i) ? "Up" : "Dn"; linkdims=20)
-  e_jw, psi_jw = ITensorMPS.dmrg(Hline, psiline; nsweeps, maxdim, cutoff, outputlevel=0)
-  ITensors.enable_auto_fermion()
-
-  # now get auto-fermion results 
-  H = ttn(os, s)
-  # make init_state
-  d = Dict()
-  for (i, v) in enumerate(vertices(s))
-    d[v] = isodd(i) ? "Up" : "Dn"
-  end
-  states = v -> d[v]
-  psi = ttn(states, s)
-  e, psi = dmrg(
-    H, psi; nsweeps, maxdim, cutoff, nsites, updater_kwargs=(; krylovdim=3, maxiter=1)
-  )
-
-  # Compare to `ITensors.MPO` version of `dmrg`
-  Hline = ITensorMPS.MPO(replace_vertices(v -> vmap[v], os), sline)
-  rng = StableRNG(1234)
-  psiline = ITensorMPS.random_mps(rng, sline, i -> isodd(i) ? "Up" : "Dn"; linkdims=20)
-  e2, psi2 = ITensorMPS.dmrg(Hline, psiline; nsweeps, maxdim, cutoff, outputlevel=0)
-
-  @test inner(psi', H, psi) ≈ ITensorMPS.inner(psi2', Hline, psi2) atol = 1e-5
-  @test e2 ≈ e_jw atol = 1e-5
-  @test inner(psi2', Hline, psi2) ≈ e_jw atol = 1e-5
-
-  if !auto_fermion_enabled
-    ITensors.disable_auto_fermion()
-  end
-end
-
-@testset "Regression test: tree truncation" begin
-  maxdim = 4
-  nsites = 2
-  nsweeps = 10
-
-  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)
-
-  @test all(edge_data(linkdims(psi)) .<= maxdim)
-end
-end
diff --git a/test/test_ttn_dmrg_x.jl b/test/test_ttn_dmrg_x.jl
deleted file mode 100644
index 4f2583ac..00000000
--- a/test/test_ttn_dmrg_x.jl
+++ /dev/null
@@ -1,70 +0,0 @@
-@eval module $(gensym())
-using Dictionaries: Dictionary
-using Graphs: nv, vertices
-using ITensorNetworks:
-  OpSum, ttn, apply, contract, dmrg_x, inner, mpo, mps, random_mps, siteinds
-using ITensorNetworks.ModelHamiltonians: ModelHamiltonians
-using ITensors: @disable_warn_order, array, dag, onehot, uniqueind
-using LinearAlgebra: eigen, normalize
-using NamedGraphs.NamedGraphGenerators: named_comb_tree
-using StableRNGs: StableRNG
-using Test: @test, @testset
-# TODO: Combine MPS and TTN tests.
-@testset "MPS DMRG-X" for conserve_qns in (false, true)
-  n = 10
-  s = siteinds("S=1/2", n; conserve_qns)
-  W = 12
-  # Random fields h ∈ [-W, W]
-  rng = StableRNG(1234)
-  h = W * (2 * rand(rng, n) .- 1)
-  H = mpo(ModelHamiltonians.heisenberg(n; h), s)
-  ψ = mps(v -> rand(rng, ["↑", "↓"]), s)
-  dmrg_x_kwargs = (nsweeps=20, normalize=true, maxdim=20, cutoff=1e-10, outputlevel=0)
-  e, ϕ = dmrg_x(H, ψ; nsites=2, dmrg_x_kwargs...)
-  @test inner(ϕ', H, ϕ) / inner(ϕ, ϕ) ≈ e
-  @test inner(ψ', H, ψ) / inner(ψ, ψ) ≈ inner(ϕ', H, ϕ) / inner(ϕ, ϕ) rtol = 1e-1
-  @test inner(H, ψ, H, ψ) ≉ inner(ψ', H, ψ)^2 rtol = 1e-7
-  @test inner(H, ϕ, H, ϕ) ≈ inner(ϕ', H, ϕ)^2 rtol = 1e-7
-  e, ϕ̃ = dmrg_x(H, ϕ; nsites=1, dmrg_x_kwargs...)
-  @test inner(ϕ̃', H, ϕ̃) / inner(ϕ̃, ϕ̃) ≈ e
-  @test inner(ψ', H, ψ) / inner(ψ, ψ) ≈ inner(ϕ̃', H, ϕ̃) / inner(ϕ̃, ϕ̃) rtol = 1e-1
-  @test inner(H, ϕ̃, H, ϕ̃) ≈ inner(ϕ̃', H, ϕ̃)^2 rtol = 1e-3
-  # Sometimes broken, sometimes not
-  # @test abs(loginner(ϕ̃, ϕ) / n) ≈ 0.0 atol = 1e-6
-end
-@testset "Tree DMRG-X" for conserve_qns in (false, true)
-  # TODO: Combine with tests above into a loop over graph structures.
-  tooth_lengths = fill(2, 3)
-  root_vertex = (3, 2)
-  c = named_comb_tree(tooth_lengths)
-  s = siteinds("S=1/2", c; conserve_qns)
-  W = 12
-  # Random fields h ∈ [-W, W]
-  rng = StableRNG(123)
-  h = Dictionary(vertices(c), W * (2 * rand(rng, nv(c)) .- 1))
-  H = ttn(ModelHamiltonians.heisenberg(c; h), s)
-  ψ = normalize(ttn(v -> rand(rng, ["↑", "↓"]), s))
-  dmrg_x_kwargs = (nsweeps=20, normalize=true, maxdim=20, cutoff=1e-10, outputlevel=0)
-  e, ϕ = dmrg_x(H, ψ; nsites=2, dmrg_x_kwargs...)
-  @test inner(ϕ', H, ϕ) / inner(ϕ, ϕ) ≈ e
-  @test inner(ψ', H, ψ) / inner(ψ, ψ) ≈ inner(ϕ', H, ϕ) / inner(ϕ, ϕ) rtol = 1e-1
-  @test inner(H, ψ, H, ψ) ≉ inner(ψ', H, ψ)^2 rtol = 1e-2
-  @test inner(H, ϕ, H, ϕ) ≈ inner(ϕ', H, ϕ)^2 rtol = 1e-7
-  e, ϕ̃ = dmrg_x(H, ϕ; nsites=1, dmrg_x_kwargs...)
-  @test inner(ϕ̃', H, ϕ̃) / inner(ϕ̃, ϕ̃) ≈ e
-  @test inner(ψ', H, ψ) / inner(ψ, ψ) ≈ inner(ϕ̃', H, ϕ̃) / inner(ϕ̃, ϕ̃) rtol = 1e-1
-  @test inner(H, ϕ̃, H, ϕ̃) ≈ inner(ϕ̃', H, ϕ̃)^2 rtol = 1e-6
-  # Sometimes broken, sometimes not
-  # @test abs(loginner(ϕ̃, ϕ) / nv(c)) ≈ 0.0 atol = 1e-8
-  # compare against ED
-  @disable_warn_order U0 = contract(ψ, root_vertex)
-  @disable_warn_order T = contract(H, root_vertex)
-  D, U = eigen(T; ishermitian=true)
-  u = uniqueind(U, T)
-  _, max_ind = findmax(abs, array(dag(U0) * U))
-  U_exact = U * dag(onehot(u => max_ind))
-  @disable_warn_order U_dmrgx = contract(ϕ, root_vertex)
-  @test inner(ϕ', H, ϕ) ≈ (dag(U_exact') * T * U_exact)[] atol = 1e-6
-  @test abs(inner(U_dmrgx, U_exact)) ≈ 1 atol = 1e-6
-end
-end
diff --git a/test/test_ttn_linsolve.jl b/test/test_ttn_linsolve.jl
deleted file mode 100644
index dab969ed..00000000
--- a/test/test_ttn_linsolve.jl
+++ /dev/null
@@ -1,48 +0,0 @@
-@eval module $(gensym())
-using ITensorNetworks: ITensorNetworks, OpSum, apply, dmrg, inner, mpo, random_mps, siteinds
-using KrylovKit: linsolve
-using StableRNGs: StableRNG
-using Test: @test, @test_broken, @testset
-
-@testset "Linsolve" begin
-  @testset "Linsolve Basics" begin
-    cutoff = 1E-11
-    maxdim = 8
-    nsweeps = 2
-
-    N = 8
-    # s = siteinds("S=1/2", N; conserve_qns=true)
-    s = siteinds("S=1/2", N; conserve_qns=false)
-
-    os = OpSum()
-    for j in 1:(N - 1)
-      os += 0.5, "S+", j, "S-", j + 1
-      os += 0.5, "S-", j, "S+", j + 1
-      os += "Sz", j, "Sz", j + 1
-    end
-    H = mpo(os, s)
-
-    #
-    # Test complex case
-    #
-
-    rng = StableRNG(1234)
-    ## TODO: Need to add support for `random_mps`/`random_tensornetwork` with state input.
-    ## states = [isodd(n) ? "Up" : "Dn" for n in 1:N]
-    ## x_c = random_mps(rng, states, s; link_space=4) + 0.1im * random_mps(rng, states, s; link_space=2)
-    x_c = random_mps(rng, s; link_space=4) + 0.1im * random_mps(rng, s; link_space=2)
-
-    b = apply(H, x_c; alg="fit", nsweeps=3, init=x_c) #cutoff is unsupported kwarg for apply/contract
-
-    ## TODO: Need to add support for `random_mps`/`random_tensornetwork` with state input.
-    ## x0 = random_mps(rng, states, s; link_space=10)
-    x0 = random_mps(rng, s; link_space=10)
-
-    x = @test_broken linsolve(
-      H, b, x0; cutoff, maxdim, nsweeps, updater_kwargs=(; tol=1E-6, ishermitian=true)
-    )
-
-    # @test norm(x - x_c) < 1E-3
-  end
-end
-end
diff --git a/test/test_ttn_tdvp.jl b/test/test_ttn_tdvp.jl
deleted file mode 100644
index f4426d21..00000000
--- a/test/test_ttn_tdvp.jl
+++ /dev/null
@@ -1,663 +0,0 @@
-@eval module $(gensym())
-using Graphs: dst, edges, src
-using ITensors: ITensor, contract, dag, inner, noprime, normalize, prime, scalar
-using ITensorNetworks:
-  ITensorNetworks,
-  OpSum,
-  ttn,
-  apply,
-  expect,
-  mpo,
-  mps,
-  op,
-  random_mps,
-  random_ttn,
-  siteinds,
-  tdvp
-using ITensorNetworks.ModelHamiltonians: ModelHamiltonians
-using LinearAlgebra: norm
-using NamedGraphs.NamedGraphGenerators: named_binary_tree, named_comb_tree
-using Observers: observer
-using StableRNGs: StableRNG
-using Test: @testset, @test
-@testset "MPS TDVP" begin
-  @testset "Basic TDVP" begin
-    N = 10
-    cutoff = 1e-12
-
-    s = siteinds("S=1/2", N)
-    os = OpSum()
-    for j in 1:(N - 1)
-      os += 0.5, "S+", j, "S-", j + 1
-      os += 0.5, "S-", j, "S+", j + 1
-      os += "Sz", j, "Sz", j + 1
-    end
-
-    H = mpo(os, s)
-
-    rng = StableRNG(1234)
-    ψ0 = random_mps(rng, s; link_space=10)
-
-    # Time evolve forward:
-    ψ1 = tdvp(H, -0.1im, ψ0; nsweeps=1, cutoff, nsites=1)
-    @test norm(ψ1) ≈ 1.0
-
-    ## Should lose fidelity:
-    #@test abs(inner(ψ0,ψ1)) < 0.9
-
-    # Average energy should be conserved:
-    @test real(inner(ψ1', H, ψ1)) ≈ inner(ψ0', H, ψ0)
-
-    # Time evolve backwards:
-    ψ2 = tdvp(
-      H,
-      +0.1im,
-      ψ1;
-      nsweeps=1,
-      cutoff,
-      updater_kwargs=(; krylovdim=20, maxiter=20, tol=1e-8),
-    )
-
-    @test norm(ψ2) ≈ 1.0
-
-    # Should rotate back to original state:
-    @test abs(inner(ψ0, ψ2)) > 0.99
-
-    # test different ways to specify time-step specifications
-    ψa = tdvp(H, -0.1im, ψ0; nsweeps=4, cutoff, nsites=1)
-    ψb = tdvp(H, -0.1im, ψ0; time_step=-0.025im, cutoff, nsites=1)
-    ψc = tdvp(
-      H, -0.1im, ψ0; time_step=[-0.02im, -0.03im, -0.015im, -0.035im], cutoff, nsites=1
-    )
-    ψd = tdvp(
-      H, -0.1im, ψ0; nsweeps=4, time_step=[-0.02im, -0.03im, -0.025im], cutoff, nsites=1
-    )
-    @test inner(ψa, ψb) ≈ 1.0 rtol = 1e-7
-    @test inner(ψa, ψc) ≈ 1.0 rtol = 1e-7
-    @test inner(ψa, ψd) ≈ 1.0 rtol = 1e-7
-  end
-
-  @testset "TDVP: Sum of Hamiltonians" begin
-    N = 10
-    cutoff = 1e-10
-
-    s = siteinds("S=1/2", N)
-
-    os1 = OpSum()
-    for j in 1:(N - 1)
-      os1 += 0.5, "S+", j, "S-", j + 1
-      os1 += 0.5, "S-", j, "S+", j + 1
-    end
-    os2 = OpSum()
-    for j in 1:(N - 1)
-      os2 += "Sz", j, "Sz", j + 1
-    end
-
-    H1 = mpo(os1, s)
-    H2 = mpo(os2, s)
-    Hs = [H1, H2]
-
-    rng = StableRNG(1234)
-    ψ0 = random_mps(rng, s; link_space=10)
-
-    ψ1 = tdvp(Hs, -0.1im, ψ0; nsweeps=1, cutoff, nsites=1)
-
-    @test norm(ψ1) ≈ 1.0
-
-    ## Should lose fidelity:
-    #@test abs(inner(ψ0,ψ1)) < 0.9
-
-    # Average energy should be conserved:
-    @test real(sum(H -> inner(ψ1', H, ψ1), Hs)) ≈ sum(H -> inner(ψ0', H, ψ0), Hs)
-
-    # Time evolve backwards:
-    ψ2 = tdvp(Hs, +0.1im, ψ1; nsweeps=1, cutoff)
-
-    @test norm(ψ2) ≈ 1.0
-
-    # Should rotate back to original state:
-    @test abs(inner(ψ0, ψ2)) > 0.99
-  end
-
-  @testset "Higher-Order TDVP" begin
-    N = 10
-    cutoff = 1e-12
-    order = 4
-
-    s = siteinds("S=1/2", N)
-
-    os = OpSum()
-    for j in 1:(N - 1)
-      os += 0.5, "S+", j, "S-", j + 1
-      os += 0.5, "S-", j, "S+", j + 1
-      os += "Sz", j, "Sz", j + 1
-    end
-
-    H = mpo(os, s)
-
-    rng = StableRNG(1234)
-    ψ0 = random_mps(rng, s; link_space=10)
-
-    # Time evolve forward:
-    ψ1 = tdvp(H, -0.1im, ψ0; time_step=-0.05im, order, cutoff, nsites=1)
-
-    @test norm(ψ1) ≈ 1.0
-
-    # Average energy should be conserved:
-    @test real(inner(ψ1', H, ψ1)) ≈ inner(ψ0', H, ψ0)
-
-    # Time evolve backwards:
-    ψ2 = tdvp(H, +0.1im, ψ1; time_step=+0.05im, order, cutoff)
-
-    @test norm(ψ2) ≈ 1.0
-
-    # Should rotate back to original state:
-    @test abs(inner(ψ0, ψ2)) > 0.99
-  end
-
-  @testset "Accuracy Test" begin
-    N = 4
-    tau = 0.1
-    ttotal = 1.0
-    cutoff = 1e-12
-
-    s = siteinds("S=1/2", N; conserve_qns=false)
-
-    os = OpSum()
-    for j in 1:(N - 1)
-      os += 0.5, "S+", j, "S-", j + 1
-      os += 0.5, "S-", j, "S+", j + 1
-      os += "Sz", j, "Sz", j + 1
-    end
-    H = mpo(os, s)
-    HM = contract(H)
-
-    Ut = exp(-im * tau * HM)
-
-    state = mps(n -> isodd(n) ? "Up" : "Dn", s)
-    psi2 = deepcopy(state)
-    psix = contract(state)
-
-    Sz_tdvp = Float64[]
-    Sz_tdvp2 = Float64[]
-    Sz_exact = Float64[]
-
-    c = div(N, 2)
-    Szc = op("Sz", s[c])
-
-    Nsteps = Int(ttotal / tau)
-    for step in 1:Nsteps
-      psix = noprime(Ut * psix)
-      psix /= norm(psix)
-
-      state = tdvp(
-        H,
-        -im * tau,
-        state;
-        cutoff,
-        normalize=false,
-        updater_kwargs=(; tol=1e-12, maxiter=500, krylovdim=25),
-      )
-      # TODO: What should `expect` output? Right now
-      # it outputs a dictionary.
-      push!(Sz_tdvp, real(expect("Sz", state; vertices=[c])[c]))
-
-      psi2 = tdvp(
-        H,
-        -im * tau,
-        psi2;
-        cutoff,
-        normalize=false,
-        updater_kwargs=(; tol=1e-12, maxiter=500, krylovdim=25),
-        updater=ITensorNetworks.exponentiate_updater,
-      )
-      # TODO: What should `expect` output? Right now
-      # it outputs a dictionary.
-      push!(Sz_tdvp2, real(expect("Sz", psi2; vertices=[c])[c]))
-
-      push!(Sz_exact, real(scalar(dag(prime(psix, s[c])) * Szc * psix)))
-      F = abs(scalar(dag(psix) * contract(state)))
-    end
-
-    @test norm(Sz_tdvp - Sz_exact) < 1e-5
-    @test norm(Sz_tdvp2 - Sz_exact) < 1e-5
-  end
-
-  @testset "TEBD Comparison" begin
-    N = 10
-    cutoff = 1e-12
-    tau = 0.1
-    ttotal = 1.0
-
-    s = siteinds("S=1/2", N; conserve_qns=true)
-
-    os = OpSum()
-    for j in 1:(N - 1)
-      os += 0.5, "S+", j, "S-", j + 1
-      os += 0.5, "S-", j, "S+", j + 1
-      os += "Sz", j, "Sz", j + 1
-    end
-
-    H = mpo(os, s)
-
-    gates = ITensor[]
-    for j in 1:(N - 1)
-      s1 = s[j]
-      s2 = s[j + 1]
-      hj =
-        op("Sz", s1) * op("Sz", s2) +
-        1 / 2 * op("S+", s1) * op("S-", s2) +
-        1 / 2 * op("S-", s1) * op("S+", s2)
-      Gj = exp(-1.0im * tau / 2 * hj)
-      push!(gates, Gj)
-    end
-    append!(gates, reverse(gates))
-
-    state = mps(n -> isodd(n) ? "Up" : "Dn", s)
-    phi = deepcopy(state)
-    c = div(N, 2)
-
-    #
-    # Evolve using TEBD
-    # 
-
-    Nsteps = convert(Int, ceil(abs(ttotal / tau)))
-    Sz1 = zeros(Nsteps)
-    En1 = zeros(Nsteps)
-    #Sz2 = zeros(Nsteps)
-    #En2 = zeros(Nsteps)
-
-    for step in 1:Nsteps
-      state = apply(gates, state; cutoff)
-
-      nsites = (step <= 3 ? 2 : 1)
-      phi = tdvp(
-        H,
-        -tau * im,
-        phi;
-        nsweeps=1,
-        cutoff,
-        nsites,
-        normalize=true,
-        updater_kwargs=(; krylovdim=15),
-      )
-
-      Sz1[step] = real(expect("Sz", state; vertices=[c])[c])
-      #Sz2[step] = real(expect("Sz", phi; vertices=[c])[c])
-      En1[step] = real(inner(state', H, state))
-      #En2[step] = real(inner(phi', H, phi))
-    end
-
-    #
-    # Evolve using TDVP
-    # 
-
-    phi = mps(n -> isodd(n) ? "Up" : "Dn", s)
-
-    obs = observer(
-      "Sz" => (; state) -> expect("Sz", state; vertices=[c])[c],
-      "En" => (; state) -> real(inner(state', H, state)),
-    )
-
-    phi = tdvp(
-      H,
-      -im * ttotal,
-      phi;
-      time_step=-im * tau,
-      cutoff,
-      normalize=false,
-      (sweep_observer!)=obs,
-      root_vertex=N, # defaults to 1, which breaks observer equality
-    )
-
-    Sz2 = obs.Sz
-    En2 = obs.En
-    @test norm(Sz1 - Sz2) < 1e-3
-    @test norm(En1 - En2) < 1e-3
-  end
-
-  @testset "Imaginary Time Evolution" for reverse_step in [true, false]
-    cutoff = 1e-12
-    tau = 1.0
-    ttotal = 10.0
-    N = 10
-    s = siteinds("S=1/2", N)
-
-    os = OpSum()
-    for j in 1:(N - 1)
-      os += 0.5, "S+", j, "S-", j + 1
-      os += 0.5, "S-", j, "S+", j + 1
-      os += "Sz", j, "Sz", j + 1
-    end
-
-    H = mpo(os, s)
-
-    rng = StableRNG(1234)
-    state = random_mps(rng, s; link_space=2)
-    en0 = inner(state', H, state)
-    nsites = [repeat([2], 10); repeat([1], 10)]
-    maxdim = 32
-    state = tdvp(
-      H,
-      -ttotal,
-      state;
-      time_step=(-tau),
-      maxdim,
-      cutoff,
-      nsites,
-      reverse_step,
-      normalize=true,
-      updater_kwargs=(; krylovdim=15),
-    )
-    en1 = inner(state', H, state)
-    @test en1 < en0
-  end
-
-  @testset "Observers" begin
-    N = 10
-    cutoff = 1e-12
-    tau = 0.1
-    ttotal = 1.0
-
-    s = siteinds("S=1/2", N; conserve_qns=true)
-
-    os = OpSum()
-    for j in 1:(N - 1)
-      os += 0.5, "S+", j, "S-", j + 1
-      os += 0.5, "S-", j, "S+", j + 1
-      os += "Sz", j, "Sz", j + 1
-    end
-    H = mpo(os, s)
-
-    c = div(N, 2)
-
-    #
-    # Using Observers.jl
-    # 
-
-    measure_sz(; state) = expect("Sz", state; vertices=[c])[c]
-    measure_en(; state) = real(inner(state', H, state))
-    sweep_obs = observer("Sz" => measure_sz, "En" => measure_en)
-
-    get_info(; info) = info
-    step_measure_sz(; state) = expect("Sz", state; vertices=[c])[c]
-    step_measure_en(; state) = real(inner(state', H, state))
-    region_obs = observer(
-      "Sz" => step_measure_sz, "En" => step_measure_en, "info" => get_info
-    )
-
-    state2 = mps(n -> isodd(n) ? "Up" : "Dn", s)
-    tdvp(
-      H,
-      -im * ttotal,
-      state2;
-      time_step=-im * tau,
-      cutoff,
-      normalize=false,
-      (sweep_observer!)=sweep_obs,
-      (region_observer!)=region_obs,
-      root_vertex=N, # defaults to 1, which breaks observer equality
-    )
-
-    Sz2 = sweep_obs.Sz
-    En2 = sweep_obs.En
-
-    Sz2_step = region_obs.Sz
-    En2_step = region_obs.En
-    infos = region_obs.info
-
-    #
-    # Could use ideas of other things to test here
-    #
-
-    @test all(x -> x.converged == 1, infos)
-  end
-end
-
-@testset "Tree TDVP" begin
-  @testset "Basic TDVP" for c in [named_comb_tree(fill(2, 3)), named_binary_tree(3)]
-    cutoff = 1e-12
-
-    tooth_lengths = fill(4, 4)
-    root_vertex = (1, 4)
-    c = named_comb_tree(tooth_lengths)
-    s = siteinds("S=1/2", c)
-
-    os = ModelHamiltonians.heisenberg(c)
-
-    H = ttn(os, s)
-
-    rng = StableRNG(1234)
-    ψ0 = normalize(random_ttn(rng, s))
-
-    # Time evolve forward:
-    ψ1 = tdvp(H, -0.1im, ψ0; root_vertex, nsweeps=1, cutoff, nsites=2)
-    @test norm(ψ1) ≈ 1.0
-
-    ## Should lose fidelity:
-    #@test abs(inner(ψ0,ψ1)) < 0.9
-
-    # Average energy should be conserved:
-    @test real(inner(ψ1', H, ψ1)) ≈ inner(ψ0', H, ψ0)
-
-    # Time evolve backwards:
-    ψ2 = tdvp(H, +0.1im, ψ1; nsweeps=1, cutoff)
-
-    @test norm(ψ2) ≈ 1.0
-
-    # Should rotate back to original state:
-    @test abs(inner(ψ0, ψ2)) > 0.99
-  end
-
-  @testset "TDVP: Sum of Hamiltonians" begin
-    cutoff = 1e-10
-
-    tooth_lengths = fill(2, 3)
-    c = named_comb_tree(tooth_lengths)
-    s = siteinds("S=1/2", c)
-
-    os1 = OpSum()
-    for e in edges(c)
-      os1 += 0.5, "S+", src(e), "S-", dst(e)
-      os1 += 0.5, "S-", src(e), "S+", dst(e)
-    end
-    os2 = OpSum()
-    for e in edges(c)
-      os2 += "Sz", src(e), "Sz", dst(e)
-    end
-
-    H1 = ttn(os1, s)
-    H2 = ttn(os2, s)
-    Hs = [H1, H2]
-
-    rng = StableRNG(1234)
-    ψ0 = normalize(random_ttn(rng, s; link_space=10))
-
-    ψ1 = tdvp(Hs, -0.1im, ψ0; nsweeps=1, cutoff, nsites=1)
-
-    @test norm(ψ1) ≈ 1.0
-
-    ## Should lose fidelity:
-    #@test abs(inner(ψ0,ψ1)) < 0.9
-
-    # Average energy should be conserved:
-    @test real(sum(H -> inner(ψ1', H, ψ1), Hs)) ≈ sum(H -> inner(ψ0', H, ψ0), Hs)
-
-    # Time evolve backwards:
-    ψ2 = tdvp(Hs, +0.1im, ψ1; nsweeps=1, cutoff)
-
-    @test norm(ψ2) ≈ 1.0
-
-    # Should rotate back to original state:
-    @test abs(inner(ψ0, ψ2)) > 0.99
-  end
-
-  @testset "Accuracy Test" begin
-    tau = 0.1
-    ttotal = 1.0
-    cutoff = 1e-12
-
-    tooth_lengths = fill(2, 3)
-    root_vertex = (3, 2)
-    c = named_comb_tree(tooth_lengths)
-    s = siteinds("S=1/2", c)
-
-    os = ModelHamiltonians.heisenberg(c)
-    H = ttn(os, s)
-    HM = contract(H)
-
-    Ut = exp(-im * tau * HM)
-
-    state = ttn(ComplexF64, v -> iseven(sum(isodd.(v))) ? "Up" : "Dn", s)
-    statex = contract(state)
-
-    Sz_tdvp = Float64[]
-    Sz_exact = Float64[]
-
-    c = (2, 1)
-    Szc = op("Sz", s[c])
-
-    Nsteps = Int(ttotal / tau)
-    for step in 1:Nsteps
-      statex = noprime(Ut * statex)
-      statex /= norm(statex)
-
-      state = tdvp(
-        H,
-        -im * tau,
-        state;
-        cutoff,
-        normalize=false,
-        updater_kwargs=(; tol=1e-12, maxiter=500, krylovdim=25),
-      )
-      push!(Sz_tdvp, real(expect("Sz", state; vertices=[c])[c]))
-      push!(Sz_exact, real(scalar(dag(prime(statex, s[c])) * Szc * statex)))
-      F = abs(scalar(dag(statex) * contract(state)))
-    end
-
-    @test norm(Sz_tdvp - Sz_exact) < 1e-5
-  end
-
-  # TODO: apply gates in ITensorNetworks
-
-  @testset "TEBD Comparison" begin
-    cutoff = 1e-12
-    maxdim = typemax(Int)
-    tau = 0.1
-    ttotal = 1.0
-
-    tooth_lengths = fill(2, 3)
-    c = named_comb_tree(tooth_lengths)
-    s = siteinds("S=1/2", c)
-
-    os = ModelHamiltonians.heisenberg(c)
-    H = ttn(os, s)
-
-    gates = ITensor[]
-    for e in edges(c)
-      s1 = s[src(e)]
-      s2 = s[dst(e)]
-      hj =
-        op("Sz", s1) * op("Sz", s2) +
-        1 / 2 * op("S+", s1) * op("S-", s2) +
-        1 / 2 * op("S-", s1) * op("S+", s2)
-      Gj = exp(-1.0im * tau / 2 * hj)
-      push!(gates, Gj)
-    end
-    append!(gates, reverse(gates))
-
-    state = ttn(v -> iseven(sum(isodd.(v))) ? "Up" : "Dn", s)
-    phi = copy(state)
-    c = (2, 1)
-
-    #
-    # Evolve using TEBD
-    # 
-
-    Nsteps = convert(Int, ceil(abs(ttotal / tau)))
-    Sz1 = zeros(Nsteps)
-    En1 = zeros(Nsteps)
-    Sz2 = zeros(Nsteps)
-    En2 = zeros(Nsteps)
-
-    for step in 1:Nsteps
-      state = apply(gates, state; cutoff, maxdim)
-
-      nsites = (step <= 3 ? 2 : 1)
-      phi = tdvp(
-        H,
-        -tau * im,
-        phi;
-        nsweeps=1,
-        cutoff,
-        nsites,
-        normalize=true,
-        updater_kwargs=(; krylovdim=15),
-      )
-
-      Sz1[step] = real(expect("Sz", state; vertices=[c])[c])
-      Sz2[step] = real(expect("Sz", phi; vertices=[c])[c])
-      En1[step] = real(inner(state', H, state))
-      En2[step] = real(inner(phi', H, phi))
-    end
-
-    #
-    # Evolve using TDVP
-    # 
-
-    phi = ttn(v -> iseven(sum(isodd.(v))) ? "Up" : "Dn", s)
-    obs = observer(
-      "Sz" => (; state) -> expect("Sz", state; vertices=[c])[c],
-      "En" => (; state) -> real(inner(state', H, state)),
-    )
-    phi = tdvp(
-      H,
-      -im * ttotal,
-      phi;
-      time_step=-im * tau,
-      cutoff,
-      normalize=false,
-      (sweep_observer!)=obs,
-      root_vertex=(3, 2),
-    )
-
-    @test norm(Sz1 - Sz2) < 5e-3
-    @test norm(En1 - En2) < 5e-3
-    @test abs.(last(Sz1) - last(obs.Sz)) .< 5e-3
-    @test abs.(last(Sz2) - last(obs.Sz)) .< 5e-3
-  end
-
-  @testset "Imaginary Time Evolution" for reverse_step in [true, false]
-    cutoff = 1e-12
-    tau = 1.0
-    ttotal = 50.0
-
-    tooth_lengths = fill(2, 3)
-    c = named_comb_tree(tooth_lengths)
-    s = siteinds("S=1/2", c)
-
-    os = ModelHamiltonians.heisenberg(c)
-    H = ttn(os, s)
-
-    rng = StableRNG(1234)
-    state = normalize(random_ttn(rng, s; link_space=2))
-
-    trange = 0.0:tau:ttotal
-    for (step, t) in enumerate(trange)
-      nsites = (step <= 10 ? 2 : 1)
-      state = tdvp(
-        H,
-        -tau,
-        state;
-        cutoff,
-        nsites,
-        reverse_step,
-        normalize=true,
-        updater_kwargs=(; krylovdim=15),
-      )
-    end
-
-    @test inner(state', H, state) < -2.47
-  end
-end
-end
diff --git a/test/test_ttn_tdvp_time_dependent.jl b/test/test_ttn_tdvp_time_dependent.jl
deleted file mode 100644
index 4101bc83..00000000
--- a/test/test_ttn_tdvp_time_dependent.jl
+++ /dev/null
@@ -1,236 +0,0 @@
-@eval module $(gensym())
-using ITensorNetworks: ITensorNetworks, TimeDependentSum, ttn, mpo, mps, siteinds, tdvp
-using ITensorNetworks.ModelHamiltonians: ModelHamiltonians
-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(
-  joinpath(
-    @__DIR__, "ITensorNetworksTestSolversUtils", "ITensorNetworksTestSolversUtils.jl"
-  ),
-)
-
-using .ITensorNetworksTestSolversUtils:
-  ITensorNetworksTestSolversUtils, krylov_solver, ode_solver
-
-# Functions need to be defined in global scope (outside
-# of the @testset macro)
-
-ω₁ = 0.1
-ω₂ = 0.2
-
-ode_alg = Tsit5()
-ode_kwargs = (; reltol=1e-8, abstol=1e-8)
-
-ω⃗ = [ω₁, ω₂]
-f⃗ = [t -> cos(ω * t) for ω in ω⃗]
-ode_updater_kwargs = (; f=[f⃗], solver_alg=ode_alg, ode_kwargs)
-
-function ode_updater(
-  init;
-  state!,
-  projected_operator!,
-  outputlevel,
-  which_sweep,
-  sweep_plan,
-  which_region_update,
-  internal_kwargs,
-  ode_kwargs,
-  solver_alg,
-  f,
-)
-  region = first(sweep_plan[which_region_update])
-  (; time_step, t) = internal_kwargs
-  t = isa(region, AbstractNamedEdge) ? t : t + time_step
-
-  H⃗₀ = projected_operator![]
-  result, info = ode_solver(
-    -im * TimeDependentSum(f, H⃗₀),
-    time_step,
-    init;
-    current_time=t,
-    solver_alg,
-    ode_kwargs...,
-  )
-  return result, (; info)
-end
-
-function tdvp_ode_solver(H⃗₀, ψ₀; time_step, kwargs...)
-  psi_t, info = ode_solver(
-    -im * TimeDependentSum(f⃗, H⃗₀), time_step, ψ₀; solver_alg=ode_alg, ode_kwargs...
-  )
-  return psi_t, (; info)
-end
-
-krylov_kwargs = (; tol=1e-8, krylovdim=15, eager=true)
-krylov_updater_kwargs = (; f=[f⃗], krylov_kwargs)
-
-function ITensorNetworksTestSolversUtils.krylov_solver(
-  H⃗₀, ψ₀; time_step, ishermitian=false, issymmetric=false, kwargs...
-)
-  psi_t, info = krylov_solver(
-    -im * TimeDependentSum(f⃗, H⃗₀),
-    time_step,
-    ψ₀;
-    krylov_kwargs...,
-    ishermitian,
-    issymmetric,
-  )
-  return psi_t, (; info)
-end
-
-function krylov_updater(
-  init;
-  state!,
-  projected_operator!,
-  outputlevel,
-  which_sweep,
-  sweep_plan,
-  which_region_update,
-  internal_kwargs,
-  ishermitian=false,
-  issymmetric=false,
-  f,
-  krylov_kwargs,
-)
-  (; time_step, t) = internal_kwargs
-  H⃗₀ = projected_operator![]
-  region = first(sweep_plan[which_region_update])
-  t = isa(region, AbstractNamedEdge) ? t : t + time_step
-
-  result, info = krylov_solver(
-    -im * TimeDependentSum(f, H⃗₀),
-    time_step,
-    init;
-    current_time=t,
-    krylov_kwargs...,
-    ishermitian,
-    issymmetric,
-  )
-  return result, (; info)
-end
-
-@testset "MPS: Time dependent Hamiltonian" begin
-  n = 4
-  J₁ = 1.0
-  J₂ = 0.1
-
-  time_step = 0.1
-  time_total = 1.0
-
-  nsites = 2
-  maxdim = 100
-  cutoff = 1e-8
-
-  s = siteinds("S=1/2", n)
-  ℋ₁₀ = ModelHamiltonians.heisenberg(n; J1=J₁, J2=0.0)
-  ℋ₂₀ = ModelHamiltonians.heisenberg(n; J1=0.0, J2=J₂)
-  ℋ⃗₀ = [ℋ₁₀, ℋ₂₀]
-  H⃗₀ = [mpo(ℋ₀, s) for ℋ₀ in ℋ⃗₀]
-
-  ψ₀ = complex(mps(j -> isodd(j) ? "↑" : "↓", s))
-
-  ψₜ_ode = tdvp(
-    H⃗₀,
-    time_total,
-    ψ₀;
-    time_step,
-    maxdim,
-    cutoff,
-    nsites,
-    updater=ode_updater,
-    updater_kwargs=ode_updater_kwargs,
-  )
-
-  ψₜ_krylov = tdvp(
-    H⃗₀,
-    time_total,
-    ψ₀;
-    time_step,
-    cutoff,
-    nsites,
-    updater=krylov_updater,
-    updater_kwargs=krylov_updater_kwargs,
-  )
-
-  ψₜ_full, _ = tdvp_ode_solver(contract.(H⃗₀), contract(ψ₀); time_step=time_total)
-
-  @test norm(ψ₀) ≈ 1
-  @test norm(ψₜ_ode) ≈ 1
-  @test norm(ψₜ_krylov) ≈ 1
-  @test norm(ψₜ_full) ≈ 1
-
-  ode_err = norm(contract(ψₜ_ode) - ψₜ_full)
-  krylov_err = norm(contract(ψₜ_krylov) - ψₜ_full)
-  #ToDo: Investigate why Krylov gives better result than ODE solver
-  @test_broken krylov_err > ode_err
-  @test ode_err < 1e-2
-  @test krylov_err < 1e-2
-end
-
-@testset "TTN: Time dependent Hamiltonian" begin
-  tooth_lengths = fill(2, 3)
-  root_vertex = (3, 2)
-  c = named_comb_tree(tooth_lengths)
-  s = siteinds("S=1/2", c)
-
-  J₁ = 1.0
-  J₂ = 0.1
-
-  time_step = 0.1
-  time_total = 1.0
-
-  nsites = 2
-  maxdim = 100
-  cutoff = 1e-8
-
-  s = siteinds("S=1/2", c)
-  ℋ₁₀ = ModelHamiltonians.heisenberg(c; J1=J₁, J2=0.0)
-  ℋ₂₀ = ModelHamiltonians.heisenberg(c; J1=0.0, J2=J₂)
-  ℋ⃗₀ = [ℋ₁₀, ℋ₂₀]
-  H⃗₀ = [ttn(ℋ₀, s) for ℋ₀ in ℋ⃗₀]
-
-  ψ₀ = ttn(ComplexF64, v -> iseven(sum(isodd.(v))) ? "↑" : "↓", s)
-
-  ψₜ_ode = tdvp(
-    H⃗₀,
-    time_total,
-    ψ₀;
-    time_step,
-    maxdim,
-    cutoff,
-    nsites,
-    updater=ode_updater,
-    updater_kwargs=ode_updater_kwargs,
-  )
-
-  ψₜ_krylov = tdvp(
-    H⃗₀,
-    time_total,
-    ψ₀;
-    time_step,
-    cutoff,
-    nsites,
-    updater=krylov_updater,
-    updater_kwargs=krylov_updater_kwargs,
-  )
-  ψₜ_full, _ = tdvp_ode_solver(contract.(H⃗₀), contract(ψ₀); time_step=time_total)
-
-  @test norm(ψ₀) ≈ 1
-  @test norm(ψₜ_ode) ≈ 1
-  @test norm(ψₜ_krylov) ≈ 1
-  @test norm(ψₜ_full) ≈ 1
-
-  ode_err = norm(contract(ψₜ_ode) - ψₜ_full)
-  krylov_err = norm(contract(ψₜ_krylov) - ψₜ_full)
-  #ToDo: Investigate why Krylov gives better result than ODE solver
-  @test_broken krylov_err > ode_err
-  @test ode_err < 1e-2
-  @test krylov_err < 1e-2
-end
-end