From 706b7aa6a6b5ec576206a83f2f1d53850e4e13a1 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Fri, 5 Apr 2024 14:05:26 -0400
Subject: [PATCH] Rewrite tensor network constructors

---
 src/abstractitensornetwork.jl                 |  18 +-
 src/itensornetwork.jl                         | 325 ++++++++-------
 src/itensors.jl                               |   3 -
 src/solvers/insert/insert.jl                  |   9 +-
 src/specialitensornetworks.jl                 |  10 +-
 .../abstracttreetensornetwork.jl              | 240 ++++++-----
 src/treetensornetworks/ttn.jl                 | 377 +++++++++---------
 test/test_additensornetworks.jl               |   4 +-
 test/test_itensornetwork.jl                   |   8 +-
 test/test_ttno.jl                             |   4 +-
 test/test_ttns.jl                             |   4 +-
 11 files changed, 536 insertions(+), 466 deletions(-)

diff --git a/src/abstractitensornetwork.jl b/src/abstractitensornetwork.jl
index 2a6ed37b..756ec055 100644
--- a/src/abstractitensornetwork.jl
+++ b/src/abstractitensornetwork.jl
@@ -116,7 +116,8 @@ end
 # TODO: broadcasting
 
 function Base.union(tn1::AbstractITensorNetwork, tn2::AbstractITensorNetwork; kwargs...)
-  tn = ITensorNetwork(union(data_graph(tn1), data_graph(tn2)); kwargs...)
+  # TODO: Use a different constructor call here?
+  tn = _ITensorNetwork(union(data_graph(tn1), data_graph(tn2)); kwargs...)
   # Add any new edges that are introduced during the union
   for v1 in vertices(tn1)
     for v2 in vertices(tn2)
@@ -129,7 +130,8 @@ function Base.union(tn1::AbstractITensorNetwork, tn2::AbstractITensorNetwork; kw
 end
 
 function NamedGraphs.rename_vertices(f::Function, tn::AbstractITensorNetwork)
-  return ITensorNetwork(rename_vertices(f, data_graph(tn)))
+  # TODO: Use a different constructor call here?
+  return _ITensorNetwork(rename_vertices(f, data_graph(tn)))
 end
 
 #
@@ -736,7 +738,8 @@ function norm_network(tn::AbstractITensorNetwork)
     setindex_preserve_graph!(tndag, dag(tndag[v]), v)
   end
   tnket = rename_vertices(v -> (v, 2), data_graph(prime(tndag; sites=[])))
-  tntn = ITensorNetwork(union(tnbra, tnket))
+  # TODO: Use a different constructor here?
+  tntn = _ITensorNetwork(union(tnbra, tnket))
   for v in vertices(tn)
     if !isempty(commoninds(tntn[(v, 1)], tntn[(v, 2)]))
       add_edge!(tntn, (v, 1) => (v, 2))
@@ -809,6 +812,9 @@ end
 
 Base.show(io::IO, graph::AbstractITensorNetwork) = show(io, MIME"text/plain"(), graph)
 
+# TODO: Move to an `ITensorNetworksVisualizationInterfaceExt`
+# package extension (and define a `VisualizationInterface` package
+# based on `ITensorVisualizationCore`.).
 function ITensorVisualizationCore.visualize(
   tn::AbstractITensorNetwork,
   args...;
@@ -865,6 +871,7 @@ function site_combiners(tn::AbstractITensorNetwork{V}) where {V}
   return Cs
 end
 
+# TODO: Combine with `insert_links`.
 function insert_missing_internal_inds(
   tn::AbstractITensorNetwork, edges; internal_inds_space=trivial_space(tn)
 )
@@ -880,12 +887,17 @@ function insert_missing_internal_inds(
   return tn
 end
 
+# TODO: Combine with `insert_links`.
 function insert_missing_internal_inds(
   tn::AbstractITensorNetwork; internal_inds_space=trivial_space(tn)
 )
   return insert_internal_inds(tn, edges(tn); internal_inds_space)
 end
 
+# TODO: What to output? Could be an `IndsNetwork`. Or maybe
+# that would be a different function `commonindsnetwork`.
+# Even in that case, this could output a `Dictionary`
+# from the edges to the common inds on that edge.
 function ITensors.commoninds(tn1::AbstractITensorNetwork, tn2::AbstractITensorNetwork)
   inds = Index[]
   for v1 in vertices(tn1)
diff --git a/src/itensornetwork.jl b/src/itensornetwork.jl
index 5e84138e..1b8ac152 100644
--- a/src/itensornetwork.jl
+++ b/src/itensornetwork.jl
@@ -10,8 +10,8 @@ struct Private end
 """
 struct ITensorNetwork{V} <: AbstractITensorNetwork{V}
   data_graph::DataGraph{V,ITensor,ITensor,NamedGraph{V},NamedEdge{V}}
-  function ITensorNetwork{V}(::Private, data_graph::DataGraph) where {V}
-    return new{V}(data_graph)
+  global function _ITensorNetwork(data_graph::DataGraph)
+    return new{vertextype(data_graph)}(data_graph)
   end
 end
 
@@ -26,108 +26,124 @@ function DataGraphs.underlying_graph_type(TN::Type{<:ITensorNetwork})
   return fieldtype(data_graph_type(TN), :underlying_graph)
 end
 
-function ITensorNetwork{V}(data_graph::DataGraph{V}) where {V}
-  return ITensorNetwork{V}(Private(), copy(data_graph))
-end
-function ITensorNetwork{V}(data_graph::DataGraph) where {V}
-  return ITensorNetwork{V}(Private(), DataGraph{V}(data_graph))
-end
-
-ITensorNetwork(data_graph::DataGraph) = ITensorNetwork{vertextype(data_graph)}(data_graph)
+## function ITensorNetwork{V}(data_graph::DataGraph{V}) where {V}
+##   return ITensorNetwork{V}(Private(), copy(data_graph))
+## end
+## function ITensorNetwork{V}(data_graph::DataGraph) where {V}
+##   return ITensorNetwork{V}(Private(), DataGraph{V}(data_graph))
+## end
+## 
+## ITensorNetwork(data_graph::DataGraph) = ITensorNetwork{vertextype(data_graph)}(data_graph)
 
+# Versions taking vertex types.
 function ITensorNetwork{V}() where {V}
-  return ITensorNetwork{V}(data_graph_type(ITensorNetwork{V})())
+  # TODO: Is there a better way to write this?
+  return _ITensorNetwork(data_graph_type(ITensorNetwork{V})())
+end
+function ITensorNetwork{V}(tn::ITensorNetwork) where {V}
+  # TODO: Is there a better way to write this?
+  return _ITensorNetwork(DataGraph{V}(data_graph(tn)))
+end
+function ITensorNetwork{V}(g::NamedGraph) where {V}
+  # TODO: Is there a better way to write this?
+  return ITensorNetwork(NamedGraph{V}(g))
 end
 
 ITensorNetwork() = ITensorNetwork{Any}()
 
 # Conversion
 ITensorNetwork(tn::ITensorNetwork) = copy(tn)
-ITensorNetwork{V}(tn::ITensorNetwork{V}) where {V} = copy(tn)
-function ITensorNetwork{V}(tn::AbstractITensorNetwork) where {V}
-  return ITensorNetwork{V}(Private(), DataGraph{V}(data_graph(tn)))
-end
-ITensorNetwork(tn::AbstractITensorNetwork) = ITensorNetwork{vertextype(tn)}(tn)
+## ITensorNetwork{V}(tn::ITensorNetwork{V}) where {V} = copy(tn)
+## function ITensorNetwork{V}(tn::AbstractITensorNetwork) where {V}
+##   return ITensorNetwork{V}(Private(), DataGraph{V}(data_graph(tn)))
+## end
+## ITensorNetwork(tn::AbstractITensorNetwork) = ITensorNetwork{vertextype(tn)}(tn)
 
 NamedGraphs.convert_vertextype(::Type{V}, tn::ITensorNetwork{V}) where {V} = tn
 NamedGraphs.convert_vertextype(V::Type, tn::ITensorNetwork) = ITensorNetwork{V}(tn)
 
-Base.copy(tn::ITensorNetwork) = ITensorNetwork(copy(data_graph(tn)))
+Base.copy(tn::ITensorNetwork) = _ITensorNetwork(copy(data_graph(tn)))
 
 #
 # Construction from collections of ITensors
 #
 
-ITensorNetwork(vs::Vector, ts::Vector{ITensor}) = ITensorNetwork(Dictionary(vs, ts))
-
-ITensorNetwork(ts::Vector{<:Pair{<:Any,ITensor}}) = ITensorNetwork(dictionary(ts))
+## function ITensorNetwork(ts::ITensorCollection)
+##   return ITensorNetwork{keytype(ts)}(ts)
+## end
 
-function ITensorNetwork(ts::ITensorCollection)
-  return ITensorNetwork{keytype(ts)}(ts)
-end
-
-function ITensorNetwork{V}(ts::ITensorCollection) where {V}
-  g = NamedGraph{V}(collect(eachindex(ts)))
+function itensors_to_itensornetwork(ts)
+  g = NamedGraph(collect(eachindex(ts)))
   tn = ITensorNetwork(g)
   for v in vertices(g)
     tn[v] = ts[v]
   end
   return tn
 end
-
+function ITensorNetwork(ts::AbstractVector{ITensor})
+  return itensors_to_itensornetwork(ts)
+end
+function ITensorNetwork(ts::AbstractDictionary{<:Any,ITensor})
+  return itensors_to_itensornetwork(ts)
+end
+function ITensorNetwork(ts::AbstractDict{<:Any,ITensor})
+  return itensors_to_itensornetwork(ts)
+end
+function ITensorNetwork(vs::AbstractVector, ts::AbstractVector{ITensor})
+  return itensors_to_itensornetwork(Dictionary(vs, ts))
+end
+function ITensorNetwork(ts::AbstractVector{<:Pair{<:Any,ITensor}})
+  return itensors_to_itensornetwork(dictionary(ts))
+end
+# TODO: Decide what this should do, maybe it should factorize?
 function ITensorNetwork(t::ITensor)
-  ts = ITensor[t]
-  return ITensorNetwork{keytype(ts)}(ts)
+  return itensors_to_itensornetwork([t])
 end
 
 #
 # Construction from underyling named graph
 #
 
-function ITensorNetwork{V}(
-  eltype::Type, undef::UndefInitializer, graph::AbstractNamedGraph; kwargs...
-) where {V}
-  return ITensorNetwork{V}(eltype, undef, IndsNetwork{V}(graph; kwargs...))
-end
-
-function ITensorNetwork{V}(eltype::Type, graph::AbstractNamedGraph; kwargs...) where {V}
-  return ITensorNetwork{V}(eltype, IndsNetwork{V}(graph; kwargs...))
-end
-
-function ITensorNetwork{V}(
-  undef::UndefInitializer, graph::AbstractNamedGraph; kwargs...
-) where {V}
-  return ITensorNetwork{V}(undef, IndsNetwork{V}(graph; kwargs...))
-end
-
-function ITensorNetwork{V}(graph::AbstractNamedGraph; kwargs...) where {V}
-  return ITensorNetwork{V}(IndsNetwork{V}(graph; kwargs...))
-end
-
 function ITensorNetwork(
   eltype::Type, undef::UndefInitializer, graph::AbstractNamedGraph; kwargs...
 )
-  return ITensorNetwork{vertextype(graph)}(eltype, undef, graph; kwargs...)
-end
-
-function ITensorNetwork(eltype::Type, graph::AbstractNamedGraph; kwargs...)
-  return ITensorNetwork{vertextype(graph)}(eltype, graph; kwargs...)
-end
-
-function ITensorNetwork(undef::UndefInitializer, graph::AbstractNamedGraph; kwargs...)
-  return ITensorNetwork{vertextype(graph)}(undef, graph; kwargs...)
-end
-
-function ITensorNetwork(graph::AbstractNamedGraph; kwargs...)
-  return ITensorNetwork{vertextype(graph)}(graph; kwargs...)
+  return ITensorNetwork(eltype, undef, IndsNetwork(graph; kwargs...))
 end
 
 function ITensorNetwork(
-  itensor_constructor::Function, underlying_graph::AbstractNamedGraph; kwargs...
+  f, graph::AbstractNamedGraph; kwargs...
 )
-  return ITensorNetwork(itensor_constructor, IndsNetwork(underlying_graph; kwargs...))
+  return ITensorNetwork(f, IndsNetwork(graph; kwargs...))
 end
 
+function ITensorNetwork(graph::AbstractNamedGraph; kwargs...)
+  return ITensorNetwork(IndsNetwork(graph; kwargs...))
+end
+
+## function ITensorNetwork(
+##   eltype::Type, undef::UndefInitializer, graph::AbstractNamedGraph; kwargs...
+## )
+##   return ITensorNetwork(eltype, undef, graph; kwargs...)
+## end
+## 
+## function ITensorNetwork(eltype::Type, graph::AbstractNamedGraph; kwargs...)
+##   return ITensorNetwork{vertextype(graph)}(eltype, graph; kwargs...)
+## end
+## 
+## function ITensorNetwork(undef::UndefInitializer, graph::AbstractNamedGraph; kwargs...)
+##   return ITensorNetwork{vertextype(graph)}(undef, graph; kwargs...)
+## end
+## 
+## function ITensorNetwork(graph::AbstractNamedGraph; kwargs...)
+##   return ITensorNetwork{vertextype(graph)}(graph; kwargs...)
+## end
+
+## function ITensorNetwork(
+##   itensor_constructor::Function, underlying_graph::AbstractNamedGraph; kwargs...
+## )
+##   return ITensorNetwork(itensor_constructor, IndsNetwork(underlying_graph; kwargs...))
+## end
+
 #
 # Construction from underyling simple graph
 #
@@ -135,69 +151,99 @@ end
 function ITensorNetwork(
   eltype::Type, undef::UndefInitializer, graph::AbstractSimpleGraph; kwargs...
 )
-  return ITensorNetwork(eltype, undef, NamedGraph(graph); kwargs...)
+  return ITensorNetwork(eltype, undef, IndsNetwork(graph; kwargs...))
 end
 
-function ITensorNetwork(eltype::Type, graph::AbstractSimpleGraph; kwargs...)
-  return ITensorNetwork(eltype, NamedGraph(graph); kwargs...)
+function ITensorNetwork(f, graph::AbstractSimpleGraph; kwargs...)
+  return ITensorNetwork(f, IndsNetwork(graph); kwargs...)
 end
 
-function ITensorNetwork(undef::UndefInitializer, graph::AbstractSimpleGraph; kwargs...)
-  return ITensorNetwork(undef, NamedGraph(graph); kwargs...)
-end
+## function ITensorNetwork(eltype::Type, graph::AbstractSimpleGraph; kwargs...)
+##   return ITensorNetwork(eltype, IndsNetwork(graph); kwargs...)
+## end
+## 
+## function ITensorNetwork(undef::UndefInitializer, graph::AbstractSimpleGraph; kwargs...)
+##   return ITensorNetwork(undef, IndsNetwork(graph); kwargs...)
+## end
 
 function ITensorNetwork(graph::AbstractSimpleGraph; kwargs...)
-  return ITensorNetwork(NamedGraph(graph); kwargs...)
+  return ITensorNetwork(IndsNetwork(graph); kwargs...)
 end
 
-function ITensorNetwork(
-  itensor_constructor::Function, underlying_graph::AbstractSimpleGraph; kwargs...
-)
-  return ITensorNetwork(itensor_constructor, NamedGraph(underlying_graph); kwargs...)
-end
+## function ITensorNetwork(
+##   itensor_constructor::Function, underlying_graph::AbstractSimpleGraph; kwargs...
+## )
+##   return ITensorNetwork(itensor_constructor, NamedGraph(underlying_graph); kwargs...)
+## end
 
 #
 # Construction from IndsNetwork
 #
 
-function ITensorNetwork{V}(
+function ITensorNetwork(
   eltype::Type, undef::UndefInitializer, inds_network::IndsNetwork; kwargs...
-) where {V}
-  return ITensorNetwork{V}(inds_network; kwargs...) do v, inds...
-    return ITensor(eltype, undef, inds...)
+)
+  return ITensorNetwork(inds_network; kwargs...) do v
+    return (inds...) -> ITensor(eltype, undef, inds...)
   end
 end
 
-function ITensorNetwork{V}(eltype::Type, inds_network::IndsNetwork; kwargs...) where {V}
-  return ITensorNetwork{V}(inds_network; kwargs...) do v, inds...
-    return ITensor(eltype, inds...)
+function ITensorNetwork(eltype::Type, inds_network::IndsNetwork; kwargs...)
+  return ITensorNetwork(inds_network; kwargs...) do v
+    return (inds...) -> ITensor(eltype, inds...)
   end
 end
 
-function ITensorNetwork{V}(
+function ITensorNetwork(
   undef::UndefInitializer, inds_network::IndsNetwork; kwargs...
-) where {V}
-  return ITensorNetwork{V}(inds_network; kwargs...) do v, inds...
-    return ITensor(undef, inds...)
+)
+  return ITensorNetwork(inds_network; kwargs...) do v
+    return (inds...) -> ITensor(undef, inds...)
   end
 end
 
-function ITensorNetwork{V}(inds_network::IndsNetwork; kwargs...) where {V}
-  return ITensorNetwork{V}(inds_network; kwargs...) do v, inds...
-    return ITensor(inds...)
+function ITensorNetwork(inds_network::IndsNetwork; kwargs...)
+  return ITensorNetwork(inds_network; kwargs...) do v
+    return (inds...) -> ITensor(inds...)
   end
 end
 
-function ITensorNetwork{V}(
+# TODO: Handle `eltype` and `undef` through `generic_state`.
+function generic_state(f, inds...)
+  return f(inds...)
+end
+function generic_state(a::AbstractArray, inds...)
+  return itensor(a, inds...)
+end
+function generic_state(s::AbstractString, inds...)
+  return state(s, inds...)
+end
+
+# TODO: This is repeated from `ModelHamiltonians`, put into a
+# single location (such as a `MakeCallable` submodule).
+to_callable(value::Type) = value
+to_callable(value::Function) = value
+to_callable(value::AbstractDict) = Base.Fix1(getindex, value)
+to_callable(value::AbstractDictionary) = Base.Fix1(getindex, value)
+to_callable(value::AbstractArray{<:Any,N}) where {N} = Base.Fix1(getindex, value) ∘ CartesianIndex
+to_callable(value) = Returns(value)
+
+function ITensorNetwork(
+  value, inds_network::IndsNetwork; kwargs...
+)
+  return ITensorNetwork(to_callable(value), inds_network; kwargs...)
+end
+
+function ITensorNetwork(
   itensor_constructor::Function, inds_network::IndsNetwork; link_space=1, kwargs...
-) where {V}
+)
   # Graphs.jl uses `zero` to create a graph of the same type
   # without any vertices or edges.
   inds_network_merge = typeof(inds_network)(
     underlying_graph(inds_network); link_space, kwargs...
   )
   inds_network = union(inds_network_merge, inds_network)
-  tn = ITensorNetwork{V}()
+  tn = ITensorNetwork{vertextype(inds_network)}()
   for v in vertices(inds_network)
     add_vertex!(tn, v)
   end
@@ -210,54 +256,59 @@ function ITensorNetwork{V}(
       get(inds_network, edgetype(inds_network)(v, nv), indtype(inds_network)[]) for
       nv in neighbors(inds_network, v)
     ]
-    setindex_preserve_graph!(tn, itensor_constructor(v, siteinds, linkinds...), v)
+    tensor_v = generic_state(itensor_constructor(v), siteinds...)
+    # TODO: Use `insert_missing_linkinds` instead.
+    tensor_v = contract(tensor_v, onehot.(vcat(linkinds...) .=> 1)...)
+    setindex_preserve_graph!(tn, tensor_v, v)
   end
   return tn
 end
 
-function ITensorNetwork(inds_network::IndsNetwork; kwargs...)
-  return ITensorNetwork{vertextype(inds_network)}(inds_network; kwargs...)
-end
-
-function ITensorNetwork(
-  eltype::Type, undef::UndefInitializer, inds_network::IndsNetwork; kwargs...
-)
-  return ITensorNetwork{vertextype(inds_network)}(eltype, undef, inds_network; kwargs...)
-end
-
-function ITensorNetwork(eltype::Type, inds_network::IndsNetwork; kwargs...)
-  return ITensorNetwork{vertextype(inds_network)}(eltype, inds_network; kwargs...)
-end
-
-function ITensorNetwork(undef::UndefInitializer, inds_network::IndsNetwork; kwargs...)
-  return ITensorNetwork{vertextype(inds_network)}(undef, inds_network; kwargs...)
-end
-
-function ITensorNetwork(itensor_constructor::Function, inds_network::IndsNetwork; kwargs...)
-  return ITensorNetwork{vertextype(inds_network)}(
-    itensor_constructor, inds_network; kwargs...
-  )
-end
-
-# TODO: Deprecate in favor of version above? Or use keyword argument?
-# This can be handled by `ITensorNetwork((v, inds...) -> state(inds...), inds_network)`
-function ITensorNetwork(eltype::Type, is::IndsNetwork, initstate::Function)
-  ψ = ITensorNetwork(eltype, is)
-  for v in vertices(ψ)
-    ψ[v] = convert_eltype(eltype, state(initstate(v), only(is[v])))
-  end
-  ψ = insert_links(ψ, edges(is))
-  return ψ
-end
-
-function ITensorNetwork(eltype::Type, is::IndsNetwork, initstate::Union{String,Integer})
-  return ITensorNetwork(eltype, is, v -> initstate)
-end
-
-function ITensorNetwork(is::IndsNetwork, initstate::Union{String,Integer,Function})
-  return ITensorNetwork(Number, is, initstate)
-end
-
+## function ITensorNetwork(inds_network::IndsNetwork; kwargs...)
+##   return ITensorNetwork{vertextype(inds_network)}(inds_network; kwargs...)
+## end
+## 
+## function ITensorNetwork(
+##   eltype::Type, undef::UndefInitializer, inds_network::IndsNetwork; kwargs...
+## )
+##   return ITensorNetwork{vertextype(inds_network)}(eltype, undef, inds_network; kwargs...)
+## end
+## 
+## function ITensorNetwork(eltype::Type, inds_network::IndsNetwork; kwargs...)
+##   return ITensorNetwork{vertextype(inds_network)}(eltype, inds_network; kwargs...)
+## end
+## 
+## function ITensorNetwork(undef::UndefInitializer, inds_network::IndsNetwork; kwargs...)
+##   return ITensorNetwork{vertextype(inds_network)}(undef, inds_network; kwargs...)
+## end
+
+## function ITensorNetwork(itensor_constructor::Function, inds_network::IndsNetwork; kwargs...)
+##   return ITensorNetwork{vertextype(inds_network)}(
+##     itensor_constructor, inds_network; kwargs...
+##   )
+## end
+
+## # TODO: Deprecate in favor of version above? Or use keyword argument?
+## # This can be handled by `ITensorNetwork((v, inds...) -> state(inds...), inds_network)`
+## function ITensorNetwork(eltype::Type, is::IndsNetwork, initstate::Function)
+##   ψ = ITensorNetwork(eltype, is)
+##   for v in vertices(ψ)
+##     ψ[v] = convert_eltype(eltype, state(initstate(v), only(is[v])))
+##   end
+##   ψ = insert_links(ψ, edges(is))
+##   return ψ
+## end
+
+## function ITensorNetwork(eltype::Type, is::IndsNetwork, initstate::Union{String,Integer})
+##   return ITensorNetwork(eltype, is, v -> initstate)
+## end
+## 
+## function ITensorNetwork(is::IndsNetwork, initstate::Union{String,Integer,Function})
+##   return ITensorNetwork(Number, is, initstate)
+## end
+
+# TODO: Remove this in favor of `insert_missing_internal_inds`
+# or call it a different name, such as `factorize_edges`.
 function insert_links(ψ::ITensorNetwork, edges::Vector=edges(ψ); cutoff=1e-15)
   for e in edges
     # Define this to work?
diff --git a/src/itensors.jl b/src/itensors.jl
index f49bdb59..a479fbaa 100644
--- a/src/itensors.jl
+++ b/src/itensors.jl
@@ -21,9 +21,6 @@ function tensor_sum(A::ITensor, B::ITensor)
   return A + B
 end
 
-# TODO: Replace with a trait?
-const ITensorCollection = Union{Vector{ITensor},Dictionary{<:Any,ITensor}}
-
 # Patch for contraction sequences with `Key`
 # leaf values.
 # TODO: Move patch to `ITensors.jl`.
diff --git a/src/solvers/insert/insert.jl b/src/solvers/insert/insert.jl
index e17ff39c..a0250c95 100644
--- a/src/solvers/insert/insert.jl
+++ b/src/solvers/insert/insert.jl
@@ -2,7 +2,7 @@
 # are essentially inverse operations, adapted for different kinds of 
 # algorithms and networks.
 
-# sort of 2-site replacebond!; TODO: use dense TTN constructor instead
+# TODO: use dense TTN constructor to make this more general.
 function default_inserter(
   state::AbstractTTN,
   phi::ITensor,
@@ -27,9 +27,8 @@ function default_inserter(
     v = ortho_vert
   end
   state[v] = phi
-  state = set_ortho_center(state, [v])
-  @assert isortho(state) && only(ortho_center(state)) == v
-  normalize && (state[v] ./= norm(state[v]))
+  state = set_ortho_region(state, [v])
+  normalize && (state[v] /= norm(state[v]))
   return state, spec
 end
 
@@ -46,6 +45,6 @@ function default_inserter(
 )
   v = only(setdiff(support(region), [ortho]))
   state[v] *= phi
-  state = set_ortho_center(state, [v])
+  state = set_ortho_region(state, [v])
   return state, nothing
 end
diff --git a/src/specialitensornetworks.jl b/src/specialitensornetworks.jl
index 806c397f..9f41c551 100644
--- a/src/specialitensornetworks.jl
+++ b/src/specialitensornetworks.jl
@@ -8,7 +8,9 @@ RETURN A TENSOR NETWORK WITH COPY TENSORS ON EACH VERTEX.
 Note that passing a link_space will mean the indices of the resulting network don't match those of the input indsnetwork
 """
 function delta_network(eltype::Type, s::IndsNetwork; link_space=nothing)
-  return ITensorNetwork((v, inds...) -> delta(eltype, inds...), s; link_space)
+  return ITensorNetwork(s; link_space) do v
+    return (inds...) -> delta(eltype, inds...)
+  end
 end
 
 function delta_network(s::IndsNetwork; link_space=nothing)
@@ -27,8 +29,10 @@ end
 Build an ITensor network on a graph specified by the inds network s. Bond_dim is given by link_space and entries are randomised (normal distribution, mean 0 std 1)
 """
 function random_tensornetwork(eltype::Type, s::IndsNetwork; link_space=nothing)
-  return ITensorNetwork(s; link_space) do v, inds...
-    itensor(randn(eltype, dim(inds)...), inds...)
+  return ITensorNetwork(s; link_space) do v
+    # TODO: This only makes a random product state right now!
+    # The rest is padded with zeros, if the bond domension is greater than 1.
+    return (inds...) -> itensor(randn(eltype, dim(inds)...), inds...)
   end
 end
 
diff --git a/src/treetensornetworks/abstracttreetensornetwork.jl b/src/treetensornetworks/abstracttreetensornetwork.jl
index 5bd3045c..08beb04c 100644
--- a/src/treetensornetworks/abstracttreetensornetwork.jl
+++ b/src/treetensornetworks/abstracttreetensornetwork.jl
@@ -1,8 +1,8 @@
 using Graphs: has_vertex
 using NamedGraphs: edge_path, leaf_vertices, post_order_dfs_edges, post_order_dfs_vertices
 using IsApprox: IsApprox, Approx
-using ITensors: directsum, hasinds, permute, plev
-using ITensors.ITensorMPS: isortho, linkind, loginner, lognorm, orthogonalize
+using ITensors: @Algorithm_str, directsum, hasinds, permute, plev
+using ITensors.ITensorMPS: linkind, loginner, lognorm, orthogonalize
 using TupleTools: TupleTools
 
 abstract type AbstractTreeTensorNetwork{V} <: AbstractITensorNetwork{V} end
@@ -17,8 +17,8 @@ end
 # Field access
 # 
 
-ITensorNetwork(ψ::AbstractTTN) = ψ.itensor_network
-ortho_center(ψ::AbstractTTN) = ψ.ortho_center
+ITensorNetwork(tn::AbstractTTN) = error("Not implemented")
+ortho_region(tn::AbstractTTN) = error("Not implemented")
 
 function default_root_vertex(gs::AbstractGraph...)
   # @assert all(is_tree.(gs))
@@ -29,29 +29,26 @@ end
 # Orthogonality center
 # 
 
-ITensorMPS.isortho(ψ::AbstractTTN) = isone(length(ortho_center(ψ)))
-
-function set_ortho_center(ψ::AbstractTTN{V}, new_center::Vector{<:V}) where {V}
-  return typeof(ψ)(itensor_network(ψ), new_center)
+function set_ortho_region(tn::AbstractTTN, new_region)
+  return error("Not implemented")
 end
 
-reset_ortho_center(ψ::AbstractTTN) = set_ortho_center(ψ, vertices(ψ))
-
 # 
 # Orthogonalization
 # 
 
-function ITensorMPS.orthogonalize(ψ::AbstractTTN{V}, root_vertex::V; kwargs...) where {V}
-  (isortho(ψ) && only(ortho_center(ψ)) == root_vertex) && return ψ
-  if isortho(ψ)
-    edge_list = edge_path(ψ, only(ortho_center(ψ)), root_vertex)
+function ITensorMPS.orthogonalize(tn::AbstractTTN, new_ortho_region; kwargs...)
+  new_ortho_region ∈ ortho_region(tn) && return tn
+  # TODO: Rewrite this in a more general way.
+  if isone(length(ortho_region(tn)))
+    edge_list = edge_path(tn, only(ortho_region(tn)), root_vertex)
   else
-    edge_list = post_order_dfs_edges(ψ, root_vertex)
+    edge_list = post_order_dfs_edges(tn, root_vertex)
   end
   for e in edge_list
-    ψ = orthogonalize(ψ, e)
+    tn = orthogonalize(tn, e)
   end
-  return set_ortho_center(ψ, [root_vertex])
+  return set_ortho_region(tn, [root_vertex])
 end
 
 # For ambiguity error
@@ -64,14 +61,14 @@ end
 # Truncation
 # 
 
-function Base.truncate(ψ::AbstractTTN; root_vertex=default_root_vertex(ψ), kwargs...)
-  for e in post_order_dfs_edges(ψ, root_vertex)
+function Base.truncate(tn::AbstractTTN; root_vertex=default_root_vertex(tn), kwargs...)
+  for e in post_order_dfs_edges(tn, root_vertex)
     # always orthogonalize towards source first to make truncations controlled
-    ψ = orthogonalize(ψ, src(e))
-    ψ = truncate(ψ, e; kwargs...)
-    ψ = set_ortho_center(ψ, [dst(e)])
+    tn = orthogonalize(tn, src(e))
+    tn = truncate(tn, e; kwargs...)
+    tn = set_ortho_region(tn, [dst(e)])
   end
-  return ψ
+  return tn
 end
 
 # For ambiguity error
@@ -85,126 +82,126 @@ end
 
 # TODO: decide on contraction order: reverse dfs vertices or forward dfs edges?
 function NDTensors.contract(
-  ψ::AbstractTTN{V}, root_vertex::V=default_root_vertex(ψ); kwargs...
-) where {V}
-  ψ = copy(ψ)
+  tn::AbstractTTN, root_vertex=default_root_vertex(tn); kwargs...
+)
+  tn = copy(tn)
   # reverse post order vertices
-  traversal_order = reverse(post_order_dfs_vertices(ψ, root_vertex))
-  return contract(ITensorNetwork(ψ); sequence=traversal_order, kwargs...)
+  traversal_order = reverse(post_order_dfs_vertices(tn, root_vertex))
+  return contract(ITensorNetwork(tn); sequence=traversal_order, kwargs...)
   # # forward post order edges
-  # ψ = copy(ψ)
-  # for e in post_order_dfs_edges(ψ, root_vertex)
-  #   ψ = contract(ψ, e)
+  # tn = copy(tn)
+  # for e in post_order_dfs_edges(tn, root_vertex)
+  #   tn = contract(tn, e)
   # end
-  # return ψ[root_vertex]
+  # return tn[root_vertex]
 end
 
 function ITensors.inner(
-  ϕ::AbstractTTN, ψ::AbstractTTN; root_vertex=default_root_vertex(ϕ, ψ)
+  x::AbstractTTN, y::AbstractTTN; root_vertex=default_root_vertex(x, y)
 )
-  ϕᴴ = sim(dag(ϕ); sites=[])
-  ψ = sim(ψ; sites=[])
-  ϕψ = ϕᴴ ⊗ ψ
+  xᴴ = sim(dag(x); sites=[])
+  y = sim(y; sites=[])
+  xy = xᴴ ⊗ y
   # TODO: find the largest tensor and use it as
   # the `root_vertex`.
-  for e in post_order_dfs_edges(ψ, root_vertex)
-    if has_vertex(ϕψ, (src(e), 2))
-      ϕψ = contract(ϕψ, (src(e), 2) => (src(e), 1))
+  for e in post_order_dfs_edges(y, root_vertex)
+    if has_vertex(xy, (src(e), 2))
+      xy = contract(xy, (src(e), 2) => (src(e), 1))
     end
-    ϕψ = contract(ϕψ, (src(e), 1) => (dst(e), 1))
-    if has_vertex(ϕψ, (dst(e), 2))
-      ϕψ = contract(ϕψ, (dst(e), 2) => (dst(e), 1))
+    xy = contract(xy, (src(e), 1) => (dst(e), 1))
+    if has_vertex(xy, (dst(e), 2))
+      xy = contract(xy, (dst(e), 2) => (dst(e), 1))
     end
   end
-  return ϕψ[root_vertex, 1][]
+  return xy[root_vertex, 1][]
 end
 
-function LinearAlgebra.norm(ψ::AbstractTTN)
-  if isortho(ψ)
-    return norm(ψ[only(ortho_center(ψ))])
+function LinearAlgebra.norm(tn::AbstractTTN)
+  if isone(length(ortho_region(tn)))
+    return norm(tn[only(ortho_region(tn))])
   end
-  return √(abs(real(inner(ψ, ψ))))
+  return √(abs(real(inner(tn, tn))))
 end
 
 # 
 # Utility
 # 
 
-function LinearAlgebra.normalize!(ψ::AbstractTTN)
-  c = ortho_center(ψ)
-  lognorm_ψ = lognorm(ψ)
-  if lognorm_ψ == -Inf
-    return ψ
+function LinearAlgebra.normalize!(tn::AbstractTTN)
+  c = ortho_region(tn)
+  lognorm_tn = lognorm(tn)
+  if lognorm_tn == -Inf
+    return tn
   end
-  z = exp(lognorm_ψ / length(c))
+  z = exp(lognorm_tn / length(c))
   for v in c
-    ψ[v] ./= z
+    tn[v] ./= z
   end
-  return ψ
+  return tn
 end
 
-function LinearAlgebra.normalize(ψ::AbstractTTN)
-  return normalize!(copy(ψ))
+function LinearAlgebra.normalize(tn::AbstractTTN)
+  return normalize!(copy(tn))
 end
 
-function _apply_to_orthocenter!(f, ψ::AbstractTTN, x)
-  v = first(ortho_center(ψ))
-  ψ[v] = f(ψ[v], x)
-  return ψ
+function _apply_to_ortho_region!(f, tn::AbstractTTN, x)
+  v = first(ortho_region(tn))
+  tn[v] = f(tn[v], x)
+  return tn
 end
 
-function _apply_to_orthocenter(f, ψ::AbstractTTN, x)
-  return _apply_to_orthocenter!(f, copy(ψ), x)
+function _apply_to_ortho_region(f, tn::AbstractTTN, x)
+  return _apply_to_ortho_region!(f, copy(tn), x)
 end
 
-Base.:*(ψ::AbstractTTN, α::Number) = _apply_to_orthocenter(*, ψ, α)
+Base.:*(tn::AbstractTTN, α::Number) = _apply_to_ortho_region(*, tn, α)
 
-Base.:*(α::Number, ψ::AbstractTTN) = ψ * α
+Base.:*(α::Number, tn::AbstractTTN) = tn * α
 
-Base.:/(ψ::AbstractTTN, α::Number) = _apply_to_orthocenter(/, ψ, α)
+Base.:/(tn::AbstractTTN, α::Number) = _apply_to_ortho_region(/, tn, α)
 
-Base.:-(ψ::AbstractTTN) = -1 * ψ
+Base.:-(tn::AbstractTTN) = -1 * tn
 
-function LinearAlgebra.rmul!(ψ::AbstractTTN, α::Number)
-  return _apply_to_orthocenter!(*, ψ, α)
+function LinearAlgebra.rmul!(tn::AbstractTTN, α::Number)
+  return _apply_to_ortho_region!(*, tn, α)
 end
 
-function ITensorMPS.lognorm(ψ::AbstractTTN)
-  if isortho(ψ)
-    return log(norm(ψ[only(ortho_center(ψ))]))
+function ITensorMPS.lognorm(tn::AbstractTTN)
+  if isone(length(ortho_region(tn)))
+    return log(norm(tn[only(ortho_region(tn))]))
   end
-  lognorm2_ψ = loginner(ψ, ψ)
-  rtol = eps(real(scalartype(ψ))) * 10
+  lognorm2_tn = loginner(tn, tn)
+  rtol = eps(real(scalartype(tn))) * 10
   atol = rtol
-  if !IsApprox.isreal(lognorm2_ψ, Approx(; rtol=rtol, atol=atol))
+  if !IsApprox.isreal(lognorm2_tn, Approx(; rtol=rtol, atol=atol))
     @warn "log(norm²) is $lognorm2_T, which is not real up to a relative tolerance of $rtol and an absolute tolerance of $atol. Taking the real part, which may not be accurate."
   end
-  return 0.5 * real(lognorm2_ψ)
+  return 0.5 * real(lognorm2_tn)
 end
 
-function logdot(ψ1::TTNT, ψ2::TTNT; kwargs...) where {TTNT<:AbstractTTN}
-  return loginner(ψ1, ψ2; kwargs...)
+function logdot(tn1::AbstractTTN, tn2::AbstractTTN; kwargs...)
+  return loginner(tn1, tn2; kwargs...)
 end
 
 # TODO: stick with this traversal or find optimal contraction sequence?
 function ITensorMPS.loginner(
-  ψ1::TTNT, ψ2::TTNT; root_vertex=default_root_vertex(ψ1, ψ2)
-)::Number where {TTNT<:AbstractTTN}
-  N = nv(ψ1)
-  if nv(ψ2) != N
-    throw(DimensionMismatch("inner: mismatched number of vertices $N and $(nv(ψ2))"))
+  tn1::AbstractTTN, tn2::AbstractTTN; root_vertex=default_root_vertex(tn1, tn2)
+)
+  N = nv(tn1)
+  if nv(tn2) != N
+    throw(DimensionMismatch("inner: mismatched number of vertices $N and $(nv(tn2))"))
   end
-  ψ1dag = sim(dag(ψ1); sites=[])
-  traversal_order = reverse(post_order_dfs_vertices(ψ1, root_vertex))
+  tn1dag = sim(dag(tn1); sites=[])
+  traversal_order = reverse(post_order_dfs_vertices(tn1, root_vertex))
 
-  O = ψ1dag[root_vertex] * ψ2[root_vertex]
+  O = tn1dag[root_vertex] * tn2[root_vertex]
 
   normO = norm(O)
   log_inner_tot = log(normO)
   O ./= normO
 
   for v in traversal_order[2:end]
-    O = (O * ψ1dag[v]) * ψ2[v]
+    O = (O * tn1dag[v]) * tn2[v]
     normO = norm(O)
     log_inner_tot += log(normO)
     O ./= normO
@@ -216,10 +213,10 @@ function ITensorMPS.loginner(
   return log_inner_tot
 end
 
-function _add_maxlinkdims(ψs::AbstractTTN...)
-  maxdims = Dictionary{edgetype(ψs[1]),Int}()
-  for e in edges(ψs[1])
-    maxdims[e] = sum(ψ -> linkdim(ψ, e), ψs)
+function _add_maxlinkdims(tns::AbstractTTN...)
+  maxdims = Dictionary{edgetype(tns[1]),Int}()
+  for e in edges(tns[1])
+    maxdims[e] = sum(tn -> linkdim(tn, e), tns)
     maxdims[reverse(e)] = maxdims[e]
   end
   return maxdims
@@ -227,57 +224,55 @@ end
 
 # TODO: actually implement this?
 function Base.:+(
-  ::ITensors.Algorithm"densitymatrix",
-  ψs::TTNT...;
+  ::Algorithm"densitymatrix",
+  tns::AbstractTTN...;
   cutoff=1e-15,
-  root_vertex=default_root_vertex(ψs...),
+  root_vertex=default_root_vertex(tns...),
   kwargs...,
-) where {TTNT<:AbstractTTN}
+)
   return error("Not implemented (yet) for trees.")
 end
 
 function Base.:+(
-  ::ITensors.Algorithm"directsum", ψs::TTNT...; root_vertex=default_root_vertex(ψs...)
-) where {TTNT<:AbstractTTN}
-  @assert all(ψ -> nv(first(ψs)) == nv(ψ), ψs)
+  ::Algorithm"directsum", tns::AbstractTTN...; root_vertex=default_root_vertex(tns...)
+)
+  @assert all(tn -> nv(first(tns)) == nv(tn), tns)
 
   # Output state
-  ϕ = ttn(siteinds(ψs[1]))
+  tn = ttn(siteinds(tns[1]))
 
-  vs = post_order_dfs_vertices(ϕ, root_vertex)
-  es = post_order_dfs_edges(ϕ, root_vertex)
-  link_space = Dict{edgetype(ϕ),Index}()
+  vs = post_order_dfs_vertices(tn, root_vertex)
+  es = post_order_dfs_edges(tn, root_vertex)
+  link_space = Dict{edgetype(tn),Index}()
 
   for v in reverse(vs)
     edges = filter(e -> dst(e) == v || src(e) == v, es)
     dims_in = findall(e -> dst(e) == v, edges)
     dim_out = findfirst(e -> src(e) == v, edges)
-
-    ls = [Tuple(only(linkinds(ψ, e)) for e in edges) for ψ in ψs]
-    ϕv, lv = directsum((ψs[i][v] => ls[i] for i in 1:length(ψs))...; tags=tags.(first(ls)))
+    ls = [Tuple(only(linkinds(tn, e)) for e in edges) for tn in tns]
+    tnv, lv = directsum((tns[i][v] => ls[i] for i in 1:length(tns))...; tags=tags.(first(ls)))
     for din in dims_in
       link_space[edges[din]] = lv[din]
     end
     if !isnothing(dim_out)
-      ϕv = replaceind(ϕv, lv[dim_out] => dag(link_space[edges[dim_out]]))
+      tnv = replaceind(tnv, lv[dim_out] => dag(link_space[edges[dim_out]]))
     end
-
-    ϕ[v] = ϕv
+    tn[v] = tnv
   end
-  return convert(TTNT, ϕ)
+  return tn
 end
 
 # TODO: switch default algorithm once more are implemented
-function Base.:+(ψs::AbstractTTN...; alg=ITensors.Algorithm"directsum"(), kwargs...)
-  return +(ITensors.Algorithm(alg), ψs...; kwargs...)
+function Base.:+(tns::AbstractTTN...; alg=Algorithm"directsum"(), kwargs...)
+  return +(Algorithm(alg), tns...; kwargs...)
 end
 
-Base.:+(ψ::AbstractTTN) = ψ
+Base.:+(tn::AbstractTTN) = tn
 
-ITensors.add(ψs::AbstractTTN...; kwargs...) = +(ψs...; kwargs...)
+ITensors.add(tns::AbstractTTN...; kwargs...) = +(tns...; kwargs...)
 
-function Base.:-(ψ1::AbstractTTN, ψ2::AbstractTTN; kwargs...)
-  return +(ψ1, -ψ2; kwargs...)
+function Base.:-(tn1::AbstractTTN, tn2::AbstractTTN; kwargs...)
+  return +(tn1, -tn2; kwargs...)
 end
 
 function ITensors.add(tn1::AbstractTTN, tn2::AbstractTTN; kwargs...)
@@ -286,18 +281,17 @@ end
 
 # TODO: Delete this
 function ITensors.permute(
-  ψ::AbstractTTN, ::Tuple{typeof(linkind),typeof(siteinds),typeof(linkind)}
+  tn::AbstractTTN, ::Tuple{typeof(linkind),typeof(siteinds),typeof(linkind)}
 )
-  ψ̃ = copy(ψ)
-  for v in vertices(ψ)
-    ls = [only(linkinds(ψ, n => v)) for n in neighbors(ψ, v)] # TODO: won't work for multiple indices per link...
-    ss = TupleTools.sort(Tuple(siteinds(ψ, v)); by=plev)
+  tn = copy(tn)
+  for v in vertices(tn)
+    ls = [only(linkinds(tn, n => v)) for n in neighbors(tn, v)] # TODO: won't work for multiple indices per link...
+    ss = TupleTools.sort(Tuple(siteinds(tn, v)); by=plev)
     setindex_preserve_graph!(
-      ψ̃, permute(ψ[v], filter(!isnothing, (ls[1], ss..., ls[2:end]...))), v
+      tn, permute(tn[v], filter(!isnothing, (ls[1], ss..., ls[2:end]...))), v
     )
   end
-  ψ̃ = set_ortho_center(ψ̃, ortho_center(ψ))
-  return ψ̃
+  return set_ortho_region(tn, ortho_region(tn))
 end
 
 function Base.isapprox(
diff --git a/src/treetensornetworks/ttn.jl b/src/treetensornetworks/ttn.jl
index 63afbc18..5abe8185 100644
--- a/src/treetensornetworks/ttn.jl
+++ b/src/treetensornetworks/ttn.jl
@@ -1,223 +1,236 @@
 using ITensors.ITensorMPS: randomMPS, replacebond!
 using ITensors.NDTensors: truncate!
 using LinearAlgebra: normalize
-using NamedGraphs: named_path_graph
+using NamedGraphs: named_path_graph, vertextype
 using Random: randn!
 
+
 """
     TreeTensorNetwork{V} <: AbstractTreeTensorNetwork{V}
-
-# Fields
-
-- itensor_network::ITensorNetwork{V}
-- ortho_lims::Vector{V}: A vector of vertices defining the orthogonality limits.
-
 """
 struct TreeTensorNetwork{V} <: AbstractTreeTensorNetwork{V}
-  itensor_network::ITensorNetwork{V}
-  ortho_center::Vector{V}
-  function TreeTensorNetwork{V}(
-    itensor_network::ITensorNetwork, ortho_center::Vector=vertices(itensor_network)
-  ) where {V}
-    @assert is_tree(itensor_network)
-    return new{V}(itensor_network, ortho_center)
+  tensornetwork::ITensorNetwork{V}
+  ortho_region::Vector{V}
+  global function _TreeTensorNetwork(
+    tensornetwork::ITensorNetwork, ortho_region
+  )
+    @assert is_tree(tensornetwork)
+    return new{vertextype(tensornetwork)}(tensornetwork, ortho_region)
+  end
+  global function _TreeTensorNetwork(
+    tensornetwork::ITensorNetwork
+  )
+    return _TreeTensorNetwork(tensornetwork, vertices(tensornetwork))
   end
 end
 
-const TTN = TreeTensorNetwork
-
-function data_graph_type(G::Type{<:TTN})
-  return data_graph_type(fieldtype(G, :itensor_network))
-end
+TreeTensorNetwork(tn::ITensorNetwork; ortho_region=vertices(tn)) = _TreeTensorNetwork(tn, ortho_region)
 
-function Base.copy(ψ::TTN)
-  return ttn(copy(ψ.itensor_network), copy(ψ.ortho_center))
-end
+const TTN = TreeTensorNetwork
 
 # Field access
-itensor_network(ψ::TTN) = getfield(ψ, :itensor_network)
+ITensorNetwork(tn::TTN) = getfield(tn, :tensornetwork)
+ortho_region(tn::TTN) = getfield(tn, :ortho_region)
 
 # Required for `AbstractITensorNetwork` interface
-data_graph(ψ::TTN) = data_graph(itensor_network(ψ))
-
-# 
-# Constructor
-# 
-
-ttn(tn::ITensorNetwork, args...) = TTN{vertextype(tn)}(tn, args...)
+data_graph(tn::TTN) = data_graph(ITensorNetwork(tn))
 
-# catch-all for default ElType
-function ttn(g::AbstractGraph, args...; kwargs...)
-  return ttn(Float64, g, args...; kwargs...)
+function data_graph_type(G::Type{<:TTN})
+  return data_graph_type(fieldtype(G, :tensornetwork))
 end
 
-function ttn(eltype::Type{<:Number}, graph::AbstractGraph, args...; kwargs...)
-  itensor_network = ITensorNetwork(eltype, graph; kwargs...)
-  return ttn(itensor_network, args...)
+function Base.copy(tn::TTN)
+  return _TreeTensorNetwork(copy(ITensorNetwork(tn)), copy(ortho_region(tn)))
 end
 
-# construct from given state (map)
-function ttn(::Type{ElT}, is::AbstractIndsNetwork, initstate, args...) where {ElT<:Number}
-  itensor_network = ITensorNetwork(ElT, is, initstate)
-  return ttn(itensor_network, args...)
-end
+# 
+# Constructor
+# 
 
-# Constructor from a collection of ITensors.
-# TODO: Support other collections like `Dictionary`,
-# interface for custom vertex names.
-function ttn(ts::ITensorCollection)
-  return ttn(ITensorNetwork(ts))
+function set_ortho_region(tn::TTN, ortho_region)
+  return ttn(ITensorNetwork(tn); ortho_region)
 end
 
-# TODO: Implement `random_circuit_ttn` for non-trivial
-# bond dimensions and correlations.
-# TODO: Implement random_ttn for QN-Index
-function random_ttn(args...; kwargs...)
-  T = ttn(args...; kwargs...)
-  randn!.(vertex_data(T))
-  normalize!.(vertex_data(T))
-  return T
+function ttn(args...; ortho_region=nothing)
+  tn = ITensorNetwork(args...)
+  if isnothing(ortho_region)
+    ortho_region = vertices(tn)
+  end
+  return _TreeTensorNetwork(tn, ortho_region)
 end
 
-function random_mps(
-  external_inds::Vector{<:Index};
-  states=nothing,
-  internal_inds_space=trivial_space(external_inds),
-)
-  # TODO: Implement in terms of general
-  # `random_ttn` constructor
-  tn_mps = if isnothing(states)
-    randomMPS(external_inds; linkdims=internal_inds_space)
-  else
-    randomMPS(external_inds, states; linkdims=internal_inds_space)
+# Construct from dense ITensor, using IndsNetwork of site indices.
+function ttn(a::ITensor, is::IndsNetwork; ortho_region=[default_root_vertex(is)], kwargs...)
+  for v in vertices(is)
+    @assert hasinds(a, is[v])
   end
-  return ttn([tn_mps[v] for v in eachindex(tn_mps)])
-end
+  @assert ortho_region ∈ vertices(is)
+  tn = ITensorNetwork(is)
+  ortho_center = first(ortho_region)
+  for e in post_order_dfs_edges(tn, ortho_center)
+    left_inds = uniqueinds(is, e)
+    a_l, a_r = factorize(a, left_inds; tags=edge_tag(e), ortho="left", kwargs...)
+    tn[src(e)] = a_l
+    is[e] = commoninds(a_l, a_r)
+    a = a_r
+  end
+  tn[ortho_center] = a
+  ttn_a = ttn(tn)
+  return orthogonalize(ttn_a, ortho_center)
+end
+
+## ttn(tn::ITensorNetwork, args...) = TTN{vertextype(tn)}(tn, args...)
+
+## # catch-all for default ElType
+## function ttn(g::AbstractGraph, args...; kwargs...)
+##   return ttn(Float64, g, args...; kwargs...)
+## end
+
+## function ttn(eltype::Type{<:Number}, graph::AbstractGraph, args...; kwargs...)
+##   itensornetwork = ITensorNetwork(eltype, graph; kwargs...)
+##   return ttn(itensornetwork, args...)
+## end
+
+## # construct from given state (map)
+## function ttn(::Type{ElT}, is::AbstractIndsNetwork, initstate, args...) where {ElT<:Number}
+##   itensornetwork = ITensorNetwork(ElT, is, initstate)
+##   return ttn(itensornetwork, args...)
+## end
+
+## # Constructor from a collection of ITensors.
+## # TODO: Support other collections like `Dictionary`,
+## # interface for custom vertex names.
+## function ttn(ts::ITensorCollection)
+##   return ttn(ITensorNetwork(ts))
+## end
+
+## # TODO: Implement `random_circuit_ttn` for non-trivial
+## # bond dimensions and correlations.
+## # TODO: Implement random_ttn for QN-Index
+## function random_ttn(args...; kwargs...)
+##   T = ttn(args...; kwargs...)
+##   randn!.(vertex_data(T))
+##   normalize!.(vertex_data(T))
+##   return T
+## end
+
+## function random_mps(
+##   external_inds::Vector{<:Index};
+##   states=nothing,
+##   internal_inds_space=trivial_space(external_inds),
+## )
+##   # TODO: Implement in terms of general
+##   # `random_ttn` constructor
+##   tn_mps = if isnothing(states)
+##     randomMPS(external_inds; linkdims=internal_inds_space)
+##   else
+##     randomMPS(external_inds, states; linkdims=internal_inds_space)
+##   end
+##   return ttn([tn_mps[v] for v in eachindex(tn_mps)])
+## end
 
 #
 # Construction from operator (map)
 #
 
-function ttn(
-  ::Type{ElT},
-  sites_map::Pair{<:AbstractIndsNetwork,<:AbstractIndsNetwork},
-  ops::Dictionary;
-  kwargs...,
-) where {ElT<:Number}
-  s = first(sites_map) # TODO: Use the sites_map
-  N = nv(sites)
-  os = Prod{Op}()
-  for v in vertices(sites)
-    os *= Op(ops[v], v)
-  end
-  T = ttn(ElT, os, sites; kwargs...)
-  # see https://github.com/ITensor/ITensors.jl/issues/526
-  lognormT = lognorm(T)
-  T /= exp(lognormT / N) # TODO: fix broadcasting for in-place assignment
-  truncate!(T; cutoff=1e-15)
-  T *= exp(lognormT / N)
-  return T
-end
-
-function ttn(
-  ::Type{ElT},
-  sites_map::Pair{<:AbstractIndsNetwork,<:AbstractIndsNetwork},
-  fops::Function;
-  kwargs...,
-) where {ElT<:Number}
-  sites = first(sites_map) # TODO: Use the sites_map
-  ops = Dictionary(vertices(sites), map(v -> fops(v), vertices(sites)))
-  return ttn(ElT, sites, ops; kwargs...)
-end
-
-function ttn(
-  ::Type{ElT},
-  sites_map::Pair{<:AbstractIndsNetwork,<:AbstractIndsNetwork},
-  op::String;
-  kwargs...,
-) where {ElT<:Number}
-  sites = first(sites_map) # TODO: Use the sites_map
-  ops = Dictionary(vertices(sites), fill(op, nv(sites)))
-  return ttn(ElT, sites, ops; kwargs...)
-end
-
-# construct from dense ITensor, using IndsNetwork of site indices
-function ttn(A::ITensor, is::IndsNetwork; ortho_center=default_root_vertex(is), kwargs...)
-  for v in vertices(is)
-    @assert hasinds(A, is[v])
-  end
-  @assert ortho_center ∈ vertices(is)
-  ψ = ITensorNetwork(is)
-  Ã = A
-  for e in post_order_dfs_edges(ψ, ortho_center)
-    left_inds = uniqueinds(is, e)
-    L, R = factorize(Ã, left_inds; tags=edge_tag(e), ortho="left", kwargs...)
-    l = commonind(L, R)
-    ψ[src(e)] = L
-    is[e] = [l]
-    Ã = R
-  end
-  ψ[ortho_center] = Ã
-  T = ttn(ψ)
-  T = orthogonalize(T, ortho_center)
-  return T
-end
+## function ttn(
+##   ::Type{ElT},
+##   sites_map::Pair{<:AbstractIndsNetwork,<:AbstractIndsNetwork},
+##   ops::Dictionary;
+##   kwargs...,
+## ) where {ElT<:Number}
+##   s = first(sites_map) # TODO: Use the sites_map
+##   N = nv(sites)
+##   os = Prod{Op}()
+##   for v in vertices(sites)
+##     os *= Op(ops[v], v)
+##   end
+##   T = ttn(ElT, os, sites; kwargs...)
+##   # see https://github.com/ITensor/ITensors.jl/issues/526
+##   lognormT = lognorm(T)
+##   T /= exp(lognormT / N) # TODO: fix broadcasting for in-place assignment
+##   truncate!(T; cutoff=1e-15)
+##   T *= exp(lognormT / N)
+##   return T
+## end
+
+## function ttn(
+##   ::Type{ElT},
+##   sites_map::Pair{<:AbstractIndsNetwork,<:AbstractIndsNetwork},
+##   fops::Function;
+##   kwargs...,
+## ) where {ElT<:Number}
+##   sites = first(sites_map) # TODO: Use the sites_map
+##   ops = Dictionary(vertices(sites), map(v -> fops(v), vertices(sites)))
+##   return ttn(ElT, sites, ops; kwargs...)
+## end
+
+## function ttn(
+##   ::Type{ElT},
+##   sites_map::Pair{<:AbstractIndsNetwork,<:AbstractIndsNetwork},
+##   op::String;
+##   kwargs...,
+## ) where {ElT<:Number}
+##   sites = first(sites_map) # TODO: Use the sites_map
+##   ops = Dictionary(vertices(sites), fill(op, nv(sites)))
+##   return ttn(ElT, sites, ops; kwargs...)
+## end
 
 # Special constructors
 
-function mps(external_inds::Vector{<:Index}; states)
-  return mps(map(i -> [i], external_inds); states)
-end
-
-function mps(external_inds::Vector{<:Vector{<:Index}}; states)
-  g = named_path_graph(length(external_inds))
-  tn = ITensorNetwork(g)
-  for v in vertices(tn)
-    tn[v] = state(only(external_inds[v]), states(v))
-  end
-  tn = insert_missing_internal_inds(
-    tn, edges(g); internal_inds_space=trivial_space(indtype(external_inds))
-  )
-  return ttn(tn)
-end
+## function mps(external_inds::Vector{<:Index}; states)
+##   return mps(map(i -> [i], external_inds); states)
+## end
+
+## function mps(external_inds::Vector{<:Vector{<:Index}}; states)
+##   g = named_path_graph(length(external_inds))
+##   tn = ITensorNetwork(g)
+##   for v in vertices(tn)
+##     tn[v] = state(only(external_inds[v]), states(v))
+##   end
+##   tn = insert_missing_internal_inds(
+##     tn, edges(g); internal_inds_space=trivial_space(indtype(external_inds))
+##   )
+##   return ttn(tn)
+## end
 
 # 
 # Utility
 # 
 
-function ITensorMPS.replacebond!(T::TTN, edge::AbstractEdge, phi::ITensor; kwargs...)
-  ortho::String = get(kwargs, :ortho, "left")
-  swapsites::Bool = get(kwargs, :swapsites, false)
-  which_decomp::Union{String,Nothing} = get(kwargs, :which_decomp, nothing)
-  normalize::Bool = get(kwargs, :normalize, false)
-
-  indsTe = inds(T[src(edge)])
-  if swapsites
-    sb = siteinds(M, src(edge))
-    sbp1 = siteinds(M, dst(edge))
-    indsTe = replaceinds(indsTe, sb, sbp1)
-  end
-
-  L, R, spec = factorize(
-    phi, indsTe; which_decomp=which_decomp, tags=tags(T, edge), kwargs...
-  )
-
-  T[src(edge)] = L
-  T[dst(edge)] = R
-  if ortho == "left"
-    normalize && (T[dst(edge)] ./= norm(T[dst(edge)]))
-    isortho(T) && (T = set_ortho_center(T, [dst(edge)]))
-  elseif ortho == "right"
-    normalize && (T[src(edge)] ./= norm(T[src(edge)]))
-    isortho(T) && (T = set_ortho_center(T, [src(edge)]))
-  end
-  return spec
-end
-
-function ITensorMPS.replacebond!(T::TTN, edge::Pair, phi::ITensor; kwargs...)
-  return replacebond!(T, edgetype(T)(edge), phi; kwargs...)
-end
-
-function ITensorMPS.replacebond(T0::TTN, args...; kwargs...)
-  return replacebond!(copy(T0), args...; kwargs...)
-end
+## function ITensorMPS.replacebond!(T::TTN, edge::AbstractEdge, phi::ITensor; kwargs...)
+##   ortho::String = get(kwargs, :ortho, "left")
+##   swapsites::Bool = get(kwargs, :swapsites, false)
+##   which_decomp::Union{String,Nothing} = get(kwargs, :which_decomp, nothing)
+##   normalize::Bool = get(kwargs, :normalize, false)
+## 
+##   indsTe = inds(T[src(edge)])
+##   if swapsites
+##     sb = siteinds(M, src(edge))
+##     sbp1 = siteinds(M, dst(edge))
+##     indsTe = replaceinds(indsTe, sb, sbp1)
+##   end
+## 
+##   L, R, spec = factorize(
+##     phi, indsTe; which_decomp=which_decomp, tags=tags(T, edge), kwargs...
+##   )
+## 
+##   T[src(edge)] = L
+##   T[dst(edge)] = R
+##   if ortho == "left"
+##     normalize && (T[dst(edge)] ./= norm(T[dst(edge)]))
+##     isone(length(ortho_region(T))) && (T = set_ortho_region(T, [dst(edge)]))
+##   elseif ortho == "right"
+##     normalize && (T[src(edge)] ./= norm(T[src(edge)]))
+##     isone(length(ortho_region(T))) && (T = set_ortho_region(T, [src(edge)]))
+##   end
+##   return spec
+## end
+## 
+## function ITensorMPS.replacebond!(T::TTN, edge::Pair, phi::ITensor; kwargs...)
+##   return replacebond!(T, edgetype(T)(edge), phi; kwargs...)
+## end
+## 
+## function ITensorMPS.replacebond(T0::TTN, args...; kwargs...)
+##   return replacebond!(copy(T0), args...; kwargs...)
+## end
diff --git a/test/test_additensornetworks.jl b/test/test_additensornetworks.jl
index fdc497d6..e2ab2d89 100644
--- a/test/test_additensornetworks.jl
+++ b/test/test_additensornetworks.jl
@@ -10,8 +10,8 @@ using Test: @test, @testset
   Random.seed!(5623)
   g = named_grid((2, 3))
   s = siteinds("S=1/2", g)
-  ψ1 = ITensorNetwork(s, v -> "↑")
-  ψ2 = ITensorNetwork(s, v -> "↓")
+  ψ1 = ITensorNetwork(v -> "↑", s)
+  ψ2 = ITensorNetwork(v -> "↓", s)
 
   ψ_GHZ = ψ1 + ψ2
 
diff --git a/test/test_itensornetwork.jl b/test/test_itensornetwork.jl
index c9d59c74..bafa6361 100644
--- a/test/test_itensornetwork.jl
+++ b/test/test_itensornetwork.jl
@@ -126,7 +126,7 @@ using Test: @test, @test_broken, @testset
     dims = (2, 2)
     g = named_grid(dims)
     s = siteinds("S=1/2", g)
-    ψ = ITensorNetwork(s, v -> "↑")
+    ψ = ITensorNetwork(v -> "↑", s)
     tn = inner_network(ψ, ψ)
     tn_2 = contract(tn, ((1, 2), 2) => ((1, 2), 1))
     @test !has_vertex(tn_2, ((1, 2), 2))
@@ -137,7 +137,7 @@ using Test: @test, @test_broken, @testset
     dims = (2, 2)
     g = named_grid(dims)
     s = siteinds("S=1/2", g)
-    ψ = ITensorNetwork(s, v -> "↑")
+    ψ = ITensorNetwork(v -> "↑", s)
     rem_vertex!(ψ, (1, 2))
     tn = inner_network(ψ, ψ)
     @test !has_vertex(tn, ((1, 2), 1))
@@ -159,8 +159,8 @@ using Test: @test, @test_broken, @testset
       siteinds("S=1/2", named_grid((3, 3))),
     )
 
-    ψ = ITensorNetwork(g; link_space) do v, inds...
-      return itensor(randn(elt, dims(inds)...), inds...)
+    ψ = ITensorNetwork(g; link_space) do v
+      return (inds...) -> itensor(randn(elt, dims(inds)...), inds...)
     end
     @test eltype(ψ[first(vertices(ψ))]) == elt
     ψ = ITensorNetwork(g; link_space) do v, inds...
diff --git a/test/test_ttno.jl b/test/test_ttno.jl
index b83a192f..79c25175 100644
--- a/test/test_ttno.jl
+++ b/test/test_ttno.jl
@@ -1,6 +1,6 @@
 @eval module $(gensym())
 using Graphs: vertices
-using ITensorNetworks: ttn, contract, ortho_center, siteinds, union_all_inds
+using ITensorNetworks: ttn, contract, ortho_region, siteinds, union_all_inds
 using ITensors: @disable_warn_order, prime, randomITensor
 using LinearAlgebra: norm
 using NamedGraphs: named_comb_tree
@@ -27,7 +27,7 @@ using Test: @test, @testset
     O = randomITensor(sites_o...)
     # dense TTN constructor from IndsNetwork
     @disable_warn_order o1 = ttn(O, is_isp; cutoff)
-    root_vertex = only(ortho_center(o1))
+    root_vertex = only(ortho_region(o1))
     @disable_warn_order begin
       O1 = contract(o1, root_vertex)
     end
diff --git a/test/test_ttns.jl b/test/test_ttns.jl
index 24f5eeae..c9ae7344 100644
--- a/test/test_ttns.jl
+++ b/test/test_ttns.jl
@@ -1,7 +1,7 @@
 @eval module $(gensym())
 using DataGraphs: vertex_data
 using Graphs: vertices
-using ITensorNetworks: ttn, contract, ortho_center, siteinds
+using ITensorNetworks: ttn, contract, ortho_region, siteinds
 using ITensors: @disable_warn_order, randomITensor
 using LinearAlgebra: norm
 using NamedGraphs: named_comb_tree
@@ -25,7 +25,7 @@ using Test: @test, @testset
     S = randomITensor(vertex_data(is)...)
     # dense TTN constructor from IndsNetwork
     @disable_warn_order s1 = ttn(S, is; cutoff)
-    root_vertex = only(ortho_center(s1))
+    root_vertex = only(ortho_region(s1))
     @disable_warn_order begin
       S1 = contract(s1, root_vertex)
     end