Remove in-place operations.

src/ITensorNetworks.jl

@@ -55,6 +55,7 @@ include(joinpath("treetensornetwork", "opsum_to_ttno.jl"))
 include(joinpath("treetensornetwork", "abstractprojttno.jl"))
 include(joinpath("treetensornetwork", "projttno.jl"))
 include(joinpath("treetensornetwork", "projttnosum.jl"))
+include(joinpath("treetensornetwork", "projttno_apply.jl"))
 include("utility.jl")
 include("specialitensornetworks.jl")
 include("renameitensornetwork.jl")

src/abstractitensornetwork.jl

@@ -342,11 +342,11 @@ function tags(tn::AbstractITensorNetwork, edge)
   return commontags(is)
 end
 
-function svd!(tn::AbstractITensorNetwork, edge::Pair; kwargs...)
-  return svd!(tn, edgetype(tn)(edge))
+function svd(tn::AbstractITensorNetwork, edge::Pair; kwargs...)
+  return svd(tn, edgetype(tn)(edge))
 end
 
-function svd!(
+function svd(
   tn::AbstractITensorNetwork,
   edge::AbstractEdge;
   U_vertex=src(edge),
@@ -356,10 +356,11 @@ function svd!(
   v_tags=tags(tn, edge),
   kwargs...,
 )
+  tn = copy(tn)
   left_inds = uniqueinds(tn, edge)
   U, S, V = svd(tn[src(edge)], left_inds; lefttags=u_tags, right_tags=v_tags, kwargs...)
 
-  rem_vertex!(tn, src(edge)) # TODO: avoid this if we can?
+  rem_vertex!(tn, src(edge))
   add_vertex!(tn, U_vertex)
   tn[U_vertex] = U
 
@@ -372,22 +373,19 @@ function svd!(
   return tn
 end
 
-function svd(tn::AbstractITensorNetwork, edge; kwargs...)
-  return svd!(copy(tn), edge; kwargs...)
-end
-
-function qr!(
+function qr(
   tn::AbstractITensorNetwork,
   edge::AbstractEdge;
   Q_vertex=src(edge),
   R_vertex=(edge, "R"),
   tags=tags(tn, edge),
   kwargs...,
 )
+  tn = copy(tn)
   left_inds = uniqueinds(tn, edge)
   Q, R = factorize(tn[src(edge)], left_inds; tags, kwargs...)
 
-  rem_vertex!(tn, src(edge)) # TODO: avoid this if we can?
+  rem_vertex!(tn, src(edge))
   add_vertex!(tn, Q_vertex)
   tn[Q_vertex] = Q
 
@@ -397,11 +395,7 @@ function qr!(
   return tn
 end
 
-function qr(tn::AbstractITensorNetwork, edge; kwargs...)
-  return qr!(copy(tn), edge; kwargs...)
-end
-
-function factorize!(
+function factorize(
   tn::AbstractITensorNetwork,
   edge::AbstractEdge;
   X_vertex=src(edge),
@@ -417,7 +411,7 @@ function factorize!(
   left_inds = uniqueinds(tn, edge)
   X, Y = factorize(tn[src(edge)], left_inds; tags, kwargs...)
 
-  rem_vertex!(tn, src(edge)) # TODO: avoid this if we can?
+  rem_vertex!(tn, src(edge))
   add_vertex!(tn, X_vertex)
   add_vertex!(tn, Y_vertex)
 
@@ -436,16 +430,13 @@ function factorize!(
   return tn
 end
 
-function factorize(tn::AbstractITensorNetwork, edge::AbstractEdge; kwargs...)
-  return factorize!(copy(tn), edge; kwargs...)
-end
-
 # For ambiguity error; TODO: decide whether to use graph mutating methods when resulting graph is unchanged?
-function _orthogonalize_edge!(tn::AbstractITensorNetwork, edge::AbstractEdge; kwargs...)
-  # factorize!(tn, edge; kwargs...)
+function _orthogonalize_edge(tn::AbstractITensorNetwork, edge::AbstractEdge; kwargs...)
+  # tn = factorize(tn, edge; kwargs...)
   # # TODO: Implement as `only(common_neighbors(tn, src(edge), dst(edge)))`
   # new_vertex = only(neighbors(tn, src(edge)) ∩ neighbors(tn, dst(edge)))
   # return contract(tn, new_vertex => dst(edge))
+  tn = copy(tn)
   left_inds = uniqueinds(tn, edge)
   ltags = tags(tn, edge)
   X, Y = factorize(tn[src(edge)], left_inds; tags=ltags, ortho="left", kwargs...)
@@ -454,31 +445,28 @@ function _orthogonalize_edge!(tn::AbstractITensorNetwork, edge::AbstractEdge; kw
   return tn
 end
 
-function orthogonalize!(tn::AbstractITensorNetwork, edge::AbstractEdge; kwargs...)
-  return _orthogonalize_edge!(tn, edge; kwargs...)
+function orthogonalize(tn::AbstractITensorNetwork, edge::AbstractEdge; kwargs...)
+  return _orthogonalize_edge(tn, edge; kwargs...)
 end
 
-function orthogonalize!(tn::AbstractITensorNetwork, edge::Pair; kwargs...)
-  return orthogonalize!(tn, edgetype(tn)(edge); kwargs...)
+function orthogonalize(tn::AbstractITensorNetwork, edge::Pair; kwargs...)
+  return orthogonalize(tn, edgetype(tn)(edge); kwargs...)
 end
 
 # Orthogonalize an ITensorNetwork towards a source vertex, treating
 # the network as a tree spanned by a spanning tree.
 # TODO: Rename `tree_orthogonalize`.
-function orthogonalize!(ψ::AbstractITensorNetwork, source_vertex)
+function orthogonalize(ψ::AbstractITensorNetwork, source_vertex)
   spanning_tree_edges = post_order_dfs_edges(bfs_tree(ψ, source_vertex), source_vertex)
   for e in spanning_tree_edges
-    orthogonalize!(ψ, e)
+    ψ = orthogonalize(ψ, e)
   end
   return ψ
 end
 
-function orthogonalize(tn::AbstractITensorNetwork, args...; kwargs...)
-  return orthogonalize!(copy(tn), args...; kwargs...)
-end
-
 # TODO: decide whether to use graph mutating methods when resulting graph is unchanged?
-function _truncate_edge!(tn::AbstractITensorNetwork, edge::AbstractEdge; kwargs...)
+function _truncate_edge(tn::AbstractITensorNetwork, edge::AbstractEdge; kwargs...)
+  tn = copy(tn)
   left_inds = uniqueinds(tn, edge)
   ltags = tags(tn, edge)
   U, S, V = svd(tn[src(edge)], left_inds; lefttags=ltags, ortho="left", kwargs...)
@@ -487,16 +475,12 @@ function _truncate_edge!(tn::AbstractITensorNetwork, edge::AbstractEdge; kwargs.
   return tn
 end
 
-function truncate!(tn::AbstractITensorNetwork, edge::AbstractEdge; kwargs...)
-  return _truncate_edge!(tn, edge; kwargs...)
-end
-
-function truncate!(tn::AbstractITensorNetwork, edge::Pair; kwargs...)
-  return truncate!(tn, edgetype(tn)(edge); kwargs...)
+function truncate(tn::AbstractITensorNetwork, edge::AbstractEdge; kwargs...)
+  return _truncate_edge(tn, edge; kwargs...)
 end
 
-function truncate(tn::AbstractITensorNetwork, edge; kwargs...)
-  return truncate!(copy(tn), edge; kwargs...)
+function truncate(tn::AbstractITensorNetwork, edge::Pair; kwargs...)
+  return truncate(tn, edgetype(tn)(edge); kwargs...)
 end
 
 function Base.:*(c::Number, ψ::AbstractITensorNetwork)

src/exports.jl

@@ -82,6 +82,7 @@ export AbstractITensorNetwork,
   TTNO,
   ProjTTNO,
   ProjTTNOSum,
+  ProjTTNOApply,
   set_nsite,
   position,
   finite_state_machine,

src/imports.jl

@@ -34,15 +34,13 @@ import ITensors:
   # contraction
   contract,
   orthogonalize,
-  orthogonalize!,
   isortho,
   inner,
   loginner,
   norm,
   lognorm,
   expect,
   # truncation
-  truncate!,
   truncate,
   replacebond!,
   replacebond,

src/treetensornetwork/abstractprojttno.jl

@@ -131,7 +131,7 @@ end
 function position(
   P::AbstractProjTTNO{V}, psi::TTNS{V}, pos::Union{Vector{<:V},NamedEdge{V}}
 ) where {V}
-  # shift position; TODO: update for immutable struct
+  # shift position
   P = shift_position(P, pos)
   # invalidate environments corresponding to internal edges
   for e in internal_edges(P)

src/treetensornetwork/projttno_apply.jl

@@ -0,0 +1,55 @@
+struct ProjTTNOApply{V} <: AbstractProjTTNO{V}
+  pos::Union{Vector{<:V},NamedEdge{V}}
+  psi0::TTNS{V}
+  H::TTNO{V}
+  environments::Dictionary{NamedEdge{V},ITensor}
+end
+
+function ProjTTNOApply(psi0::TTNS, H::TTNO)
+  return ProjTTNOApply(vertextype(H)[], psi0, H, Dictionary{edgetype(H),ITensor}())
+end
+
+function copy(P::ProjTTNOApply)
+  return ProjTTNOApply(P.pos, copy(P.psi0), copy(P.H), copy(P.environments))
+end
+
+function set_nsite(P::ProjTTNOApply, nsite)
+  return P
+end
+
+function shift_position(P::ProjTTNOApply, pos)
+  return ProjTTNOApply(pos, P.psi0, P.H, P.environments)
+end
+
+function make_environment!(P::ProjTTNOApply{V}, psi::TTNS{V}, e::NamedEdge{V})::ITensor where {V}
+  # invalidate environment for opposite edge direction if necessary
+  reverse(e) ∈ incident_edges(P) || unset!(P.environments, reverse(e))
+  # do nothing if valid environment already present
+  if haskey(P.environments, e)
+    env = environment(P, e)
+  else
+    if is_leaf(underlying_graph(P), src(e))
+      # leaves are easy
+      env = P.psi0[src(e)] * P.H[src(e)] * dag(psi[src(e)])
+    else
+      # construct by contracting neighbors
+      neighbor_envs = ITensor[]
+      for n in setdiff(neighbors(underlying_graph(P), src(e)), [dst(e)])
+        push!(neighbor_envs, make_environment!(P, psi, edgetype(P)(n, src(e))))
+      end
+      # manually heuristic for contraction order: two environments, site tensors, then
+      # other environments
+      frst, scnd, rst = _separate_first_two(neighbor_envs)
+      itensor_map = vcat(P.psi0[src(e)], frst, scnd, P.H[src(e)], dag(psi[src(e)]), rst)
+      # TODO: actually use optimal contraction sequence here
+      env = reduce(*, itensor_map)
+    end
+    # cache
+    set!(P.environments, e, env)
+  end
+  @assert(
+    hascommoninds(environment(P, e), psi[src(e)]),
+    "Something went wrong, probably re-orthogonalized this edge in the same direction twice!"
+  )
+  return env
+end

test/test_ttns.jl

@@ -30,7 +30,6 @@ using Random
     # dense array constructor from Vector{Index} and NamedDimGraph
     @disable_warn_order s4 = TTNS(AS, sites_s, c; vertex_order, cutoff)
     # see if this actually worked
-    @show ortho_center(s1)
     root_vertex = only(ortho_center(s1))
     @disable_warn_order begin
       S1 = contract(s1, root_vertex)