Remove support for MPS/MPO, Sweeps, AbstractObserver (#50)

Author: mtfishman

Project.toml

@@ -4,6 +4,7 @@ authors = ["Matthew Fishman <mfishman@flatironinstitute.org> and contributors"]
 version = "0.2.2"
 
 [deps]
+AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
 Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
 DataGraphs = "b5a273c3-7e6c-41f6-98bd-8d7f1525a36a"
 Dictionaries = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4"

README.md

@@ -273,8 +273,8 @@ julia> @visualize Z;
 
 julia> contraction_sequence(Z)
 2-element Vector{Vector}:
- [(1, 1), (1, 2)]
- Any[(2, 1), Any[(2, 2), [(3, 1), (3, 2)]]]
+ Key{Tuple{Int64, Int64}}[Key((1, 1)), Key((1, 2))]
+ Any[Key((2, 1)), Any[Key((2, 2)), Key{Tuple{Int64, Int64}}[Key((3, 1)), Key((3, 2))]]]
 
 julia> Z̃ = contract(Z, (1, 1) => (2, 1));
 

examples/treetensornetworks/solvers/04_tdvp_observers.jl

@@ -36,7 +36,7 @@ end
 
 function measure_sz(; psi, bond, half_sweep)
   if bond == 1 && half_sweep == 2
-    return expect(psi, "Sz"; sites=N ÷ 2)
+    return expect(psi, "Sz"; vertices=[N ÷ 2])
   end
   return nothing
 end

examples/treetensornetworks/solvers/05_tdvp_nonuniform_timesteps.jl

@@ -37,11 +37,11 @@ println("\nResults")
 println("=======")
 print("step = ", 0)
 print(", time = ", zero(ComplexF64))
-print(", ⟨Sᶻ⟩ = ", round(expect(psi0, "Sz"; sites=N ÷ 2); digits=3))
+print(", ⟨Sᶻ⟩ = ", round(expect(psi0, "Sz"; vertices=[N ÷ 2]); digits=3))
 println()
 for n in 1:length(times)
   print("step = ", n)
   print(", time = ", round(times[n]; digits=3))
-  print(", ⟨Sᶻ⟩ = ", round(expect(psis[n], "Sz"; sites=N ÷ 2); digits=3))
+  print(", ⟨Sᶻ⟩ = ", round(expect(psis[n], "Sz"; vertices=[N ÷ 2]); digits=3))
   println()
 end

src/ITensorNetworks.jl

@@ -1,5 +1,6 @@
 module ITensorNetworks
 
+using AbstractTrees
 using Compat
 using DataGraphs
 using Dictionaries
@@ -10,6 +11,7 @@ using IsApprox
 using ITensors
 using ITensors.ContractionSequenceOptimization
 using ITensors.ITensorVisualizationCore
+using ITensors.LazyApply
 using IterTools
 using KrylovKit: KrylovKit
 using NamedGraphs
@@ -56,6 +58,7 @@ function iterate(::AbstractDataGraph)
   )
 end
 
+include("observers.jl")
 include("utils.jl")
 include("visualize.jl")
 include("graphs.jl")
@@ -65,7 +68,7 @@ include("lattices.jl")
 include("abstractindsnetwork.jl")
 include("indextags.jl")
 include("indsnetwork.jl")
-include("opsum.jl") # Requires IndsNetwork
+include("opsum.jl")
 include("sitetype.jl")
 include("abstractitensornetwork.jl")
 include("contraction_sequences.jl")
@@ -86,14 +89,6 @@ include(joinpath("treetensornetworks", "projttns", "abstractprojttn.jl"))
 include(joinpath("treetensornetworks", "projttns", "projttn.jl"))
 include(joinpath("treetensornetworks", "projttns", "projttnsum.jl"))
 include(joinpath("treetensornetworks", "projttns", "projttn_apply.jl"))
-# Compatibility of ITensors.MPS/MPO with tree sweeping routines
-include(joinpath("treetensornetworks", "solvers", "tree_patch.jl"))
-# Compatibility of ITensor observer and Observers
-# TODO: Delete this
-include(joinpath("treetensornetworks", "solvers", "update_observer.jl"))
-# Utilities for making it easier
-# to define solvers (like ODE solvers)
-# for TDVP
 include(joinpath("treetensornetworks", "solvers", "solver_utils.jl"))
 include(joinpath("treetensornetworks", "solvers", "applyexp.jl"))
 include(joinpath("treetensornetworks", "solvers", "tdvporder.jl"))
@@ -103,10 +98,7 @@ include(joinpath("treetensornetworks", "solvers", "tdvp_generic.jl"))
 include(joinpath("treetensornetworks", "solvers", "tdvp.jl"))
 include(joinpath("treetensornetworks", "solvers", "dmrg.jl"))
 include(joinpath("treetensornetworks", "solvers", "dmrg_x.jl"))
-include(joinpath("treetensornetworks", "solvers", "projmpo_apply.jl"))
-include(joinpath("treetensornetworks", "solvers", "contract_operator_state.jl"))
-include(joinpath("treetensornetworks", "solvers", "projmps2.jl"))
-include(joinpath("treetensornetworks", "solvers", "projmpo_mps2.jl"))
+include(joinpath("treetensornetworks", "solvers", "contract.jl"))
 include(joinpath("treetensornetworks", "solvers", "linsolve.jl"))
 include(joinpath("treetensornetworks", "solvers", "tree_sweeping.jl"))
 

src/abstractitensornetwork.jl

@@ -118,6 +118,10 @@ function LinearAlgebra.promote_leaf_eltypes(tn::AbstractITensorNetwork)
   return LinearAlgebra.promote_leaf_eltypes(itensors(tn))
 end
 
+function trivial_space(tn::AbstractITensorNetwork)
+  return trivial_space(tn[first(vertices(tn))])
+end
+
 function ITensors.promote_itensor_eltype(tn::AbstractITensorNetwork)
   return LinearAlgebra.promote_leaf_eltypes(tn)
 end
@@ -133,13 +137,18 @@ function ITensors.convert_leaf_eltype(eltype::Type, tn::AbstractITensorNetwork)
   return tn
 end
 
-# TODO: mimic ITensors.AbstractMPS implementation using map
+# TODO: Mimic ITensors.AbstractMPS implementation using map
+# TODO: Implement using `adapt`
 function NDTensors.convert_scalartype(eltype::Type{<:Number}, tn::AbstractITensorNetwork)
   tn = copy(tn)
   vertex_data(tn) .= ITensors.adapt.(Ref(eltype), vertex_data(tn))
   return tn
 end
 
+function Base.complex(tn::AbstractITensorNetwork)
+  return NDTensors.convert_scalartype(complex(LinearAlgebra.promote_leaf_eltypes(tn)), tn)
+end
+
 #
 # Conversion to Graphs
 #
@@ -168,22 +177,48 @@ function IndsNetwork(tn::AbstractITensorNetwork)
   return is
 end
 
-function siteinds(tn::AbstractITensorNetwork)
+# Alias
+indsnetwork(tn::AbstractITensorNetwork) = IndsNetwork(tn)
+
+function external_indsnetwork(tn::AbstractITensorNetwork)
   is = IndsNetwork(underlying_graph(tn))
   for v in vertices(tn)
     is[v] = uniqueinds(tn, v)
   end
   return is
 end
 
-function linkinds(tn::AbstractITensorNetwork)
+# For backwards compatibility
+# TODO: Delete this
+siteinds(tn::AbstractITensorNetwork) = external_indsnetwork(tn)
+
+# External indsnetwork of the flattened network, with vertices
+# mapped back to `tn1`.
+function flatten_external_indsnetwork(
+  tn1::AbstractITensorNetwork,
+  tn2::AbstractITensorNetwork,
+)
+  is = external_indsnetwork(sim(tn1; sites=[]) ⊗ tn2)
+  flattened_is = IndsNetwork(underlying_graph(tn1))
+  for v in vertices(flattened_is)
+    # setindex_preserve_graph!(flattened_is, unioninds(is[v, 1], is[v, 2]), v)
+    flattened_is[v] = unioninds(is[v, 1], is[v, 2])
+  end
+  return flattened_is
+end
+
+function internal_indsnetwork(tn::AbstractITensorNetwork)
   is = IndsNetwork(underlying_graph(tn))
   for e in edges(tn)
     is[e] = commoninds(tn, e)
   end
   return is
 end
 
+# For backwards compatibility
+# TODO: Delete this
+linkinds(tn::AbstractITensorNetwork) = internal_indsnetwork(tn)
+
 #
 # Index access
 #
@@ -313,6 +348,22 @@ end
 # TODO: how to define this lazily?
 #norm(tn::AbstractITensorNetwork) = sqrt(inner(tn, tn))
 
+function isapprox(
+  x::AbstractITensorNetwork,
+  y::AbstractITensorNetwork;
+  atol::Real=0,
+  rtol::Real=Base.rtoldefault(
+    LinearAlgebra.promote_leaf_eltypes(x), LinearAlgebra.promote_leaf_eltypes(y), atol
+  ),
+)
+  error("Not implemented")
+  d = norm(x - y)
+  if !isfinite(d)
+    error("In `isapprox(x::AbstractITensorNetwork, y::AbstractITensorNetwork)`, `norm(x - y)` is not finite")
+  end
+  return d <= max(atol, rtol * max(norm(x), norm(y)))
+end
+
 function contract(tn::AbstractITensorNetwork; sequence=vertices(tn), kwargs...)
   sequence_linear_index = deepmap(v -> vertex_to_parent_vertex(tn, v), sequence)
   return contract(Vector{ITensor}(tn); sequence=sequence_linear_index, kwargs...)
@@ -739,6 +790,23 @@ function site_combiners(tn::AbstractITensorNetwork{V}) where {V}
   return Cs
 end
 
+function insert_missing_internal_inds(tn::AbstractITensorNetwork, edges; internal_inds_space=trivial_space(tn))
+  tn = copy(tn)
+  for e in edges
+    if !hascommoninds(tn[src(e)], tn[dst(e)])
+      iₑ = Index(internal_inds_space, edge_tag(e))
+      X = onehot(iₑ => 1)
+      tn[src(e)] *= X
+      tn[dst(e)] *= dag(X)
+    end
+  end
+  return tn
+end
+
+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: should this make sure that internal indices
 ## # don't clash?
 ## function hvncat(

src/apply.jl

@@ -47,7 +47,7 @@ function ITensors.apply(
   o⃗::Vector{ITensor},
   ψ::AbstractITensorNetwork;
   cutoff,
-  maxdim,
+  maxdim=typemax(Int),
   normalize=false,
   ortho=false,
 )

src/contraction_sequences.jl

@@ -4,8 +4,8 @@ end
 
 function contraction_sequence(tn::AbstractITensorNetwork; kwargs...)
   seq_linear_index = contraction_sequence(Vector{ITensor}(tn); kwargs...)
-  # TODO: use Functors.fmap
-  return deepmap(n -> vertices(tn)[n], seq_linear_index)
+  # TODO: Use Functors.fmap?
+  return deepmap(n -> Key(vertices(tn)[n]), seq_linear_index)
 end
 
 function contraction_sequence(::Algorithm"optimal", tn::Vector{ITensor})

src/expect.jl

@@ -5,17 +5,17 @@ function expect(
   maxdim=nothing,
   ortho=false,
   sequence=nothing,
-  sites=vertices(ψ),
+  vertices=vertices(ψ),
 )
   s = siteinds(ψ)
   ElT = promote_itensor_eltype(ψ)
-  # ElT = ishermitian(ITensors.op(op, s[sites[1]])) ? real(ElT) : ElT
-  res = Dictionary(sites, Vector{ElT}(undef, length(sites)))
+  # ElT = ishermitian(ITensors.op(op, s[vertices[1]])) ? real(ElT) : ElT
+  res = Dictionary(vertices, Vector{ElT}(undef, length(vertices)))
   if isnothing(sequence)
     sequence = contraction_sequence(inner_network(ψ, ψ; flatten=true))
   end
   normψ² = norm_sqr(ψ; sequence)
-  for v in sites
+  for v in vertices
     O = ITensor(Op(op, v), s)
     Oψ = apply(O, ψ; cutoff, maxdim, ortho)
     res[v] = contract_inner(ψ, Oψ; sequence) / normψ²

src/exports.jl

@@ -23,7 +23,8 @@ export grid,
   vertex_path
 
 # NamedGraphs
-export named_binary_tree,
+export Key,
+  named_binary_tree,
   named_grid,
   is_tree,
   parent_vertex,
@@ -45,7 +46,6 @@ export IndsNetwork, union_all_inds
 # ITensorNetworks: itensornetwork.jl
 export AbstractITensorNetwork,
   ITensorNetwork,
-  randomITensorNetwork,
   ⊗,
   itensors,
   reverse_bfs_edges,
@@ -59,13 +59,17 @@ export AbstractITensorNetwork,
   combine_linkinds,
   subgraphs,
   reverse_bfs_edges,
+  randomITensorNetwork,
+  random_mps,
   # treetensornetwork
   default_root_vertex,
+  mpo,
+  mps,
   ortho_center,
   set_ortho_center,
   TreeTensorNetwork,
   TTN,
-  randomTTN,
+  random_ttn,
   ProjTTN,
   ProjTTNSum,
   ProjTTNApply,