[WIP] Add approximation tensor network contraction subpackage

Project.toml

@@ -4,23 +4,31 @@ authors = ["Matthew Fishman <mfishman@flatironinstitute.org> and contributors"]
 version = "0.1.0"
 
 [deps]
+Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 DataGraphs = "b5a273c3-7e6c-41f6-98bd-8d7f1525a36a"
 Dictionaries = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4"
+GraphRecipes = "bd48cda9-67a9-57be-86fa-5b3c104eda73"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
+GraphsFlows = "06909019-6f44-4949-96fc-b9d9aaa02889"
 ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MultiDimDictionaries = "87ff4268-a46e-478f-b30a-76b83dd64e3c"
 NamedGraphs = "678767b0-92e7-4007-89e4-4527a8725b19"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
+SimpleWeightedGraphs = "47aef6b3-ad0c-573a-a1e2-d07658019622"
 Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
+SweepContractor = "75a5deae-e917-4509-af32-a989148c8d5f"
+TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 
 [compat]
 Dictionaries = "0.3.15"
 Graphs = "1.6.0"
 ITensors = "0.3"
 Requires = "1.3.0"
 Suppressor = "0.2.1"
-julia = "1.6"
+julia = "1.6, 1.8"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

examples/peps/utils.jl

@@ -10,7 +10,9 @@ function cartesian_to_linear(dims::Tuple)
 end
 
 NamedGraphs.NamedDimGraph(vertices::Vector) = NamedDimGraph(tuple.(vertices))
-NamedGraphs.NamedDimGraph(vertices::Vector{<:Tuple}) = NamedDimGraph(Graph(length(vertices)); vertices)
+function NamedGraphs.NamedDimGraph(vertices::Vector{<:Tuple})
+  return NamedDimGraph(Graph(length(vertices)); vertices)
+end
 
 function rename_vertices(e::AbstractEdge, name_map::Dictionary)
   return typeof(e)(name_map[src(e)], name_map[dst(e)])
@@ -92,7 +94,9 @@ function ITensors.apply(o::ITensor, ψ::ITensorNetwork; cutoff, maxdim, normaliz
       error("Vertices where the gates are being applied must be neighbors for now.")
     end
     oψᵥ = apply(o, ψ[v⃗[1]] * ψ[v⃗[2]])
-    ψᵥ₁, ψᵥ₂ = factorize(oψᵥ, inds(ψ[v⃗[1]]); cutoff, maxdim, tags=ITensorNetworks.edge_tag(e))
+    ψᵥ₁, ψᵥ₂ = factorize(
+      oψᵥ, inds(ψ[v⃗[1]]); cutoff, maxdim, tags=ITensorNetworks.edge_tag(e)
+    )
     if normalize
       ψᵥ₁ ./= norm(ψᵥ₁)
       ψᵥ₂ ./= norm(ψᵥ₂)
@@ -107,7 +111,9 @@ function ITensors.apply(o::ITensor, ψ::ITensorNetwork; cutoff, maxdim, normaliz
   return ψ
 end
 
-function ITensors.apply(o⃗::Vector{ITensor}, ψ::ITensorNetwork; cutoff, maxdim, normalize=false)
+function ITensors.apply(
+  o⃗::Vector{ITensor}, ψ::ITensorNetwork; cutoff, maxdim, normalize=false
+)
   o⃗ψ = ψ
   for oᵢ in o⃗
     o⃗ψ = apply(oᵢ, o⃗ψ; cutoff, maxdim, normalize)

src/ApproximateTNContraction/ApproximateTNContraction.jl

@@ -0,0 +1,21 @@
+@reexport module ApproximateTNContraction
+
+using ITensors
+
+using ITensors: data, contract
+
+using TimerOutputs
+
+const timer = TimerOutput()
+
+include("ITensors.jl")
+include("orthogonal_tensor.jl")
+include("networks/lattices.jl")
+include("networks/inds_network.jl")
+include("networks/itensor_network.jl")
+include("networks/3d_classical_ising.jl")
+include("models/models.jl")
+include("interfaces/sweep_contractor.jl")
+include("contract/contract.jl")
+
+end

src/ApproximateTNContraction/ITensors.jl

@@ -0,0 +1,49 @@
+#
+# ITensors.jl extensions
+#
+
+# Generalize siteind to n-dimensional lattice
+function ITensors.siteind(st::SiteType, N1::Integer, N2::Integer, Ns::Integer...; kwargs...)
+  s = siteind(st; kwargs...)
+  if !isnothing(s)
+    ts = "n1=$N1,n2=$N2"
+    for i in eachindex(Ns)
+      ts *= ",n$(i + 2)=$(Ns[i])"
+    end
+    return addtags(s, ts)
+  end
+  return isnothing(s) && error(space_error_message(st))
+end
+
+# Generalize siteinds to n-dimensional lattice
+function ITensors.siteinds(
+  str::AbstractString, N1::Integer, N2::Integer, Ns::Integer...; kwargs...
+)
+  st = SiteType(str)
+  return [siteind(st, ns...) for ns in Base.product(1:N1, 1:N2, UnitRange.(1, Ns)...)]
+end
+
+# Get the promoted type of the Index objects in a collection
+# of Index (Tuple, Vector, ITensor, etc.)
+indtype(i::Index) = typeof(i)
+indtype(T::Type{<:Index}) = T
+indtype(is::Tuple{Vararg{<:Index}}) = eltype(is)
+indtype(is::Vector{<:Index}) = eltype(is)
+indtype(A::ITensor...) = indtype(inds.(A))
+
+indtype(tn1, tn2) = promote_type(indtype(tn1), indtype(tn2))
+indtype(tn) = mapreduce(indtype, promote_type, tn)
+
+#
+# MPS functionality extensions
+#
+
+Base.keytype(m::MPS) = keytype(data(m))
+
+# A version of indexing which returns an empty order-0 ITensor
+# when out of bounds
+get_itensor(x::MPS, n::Int) = n in 1:length(x) ? x[n] : ITensor()
+
+# Reverse the site ordering of an MPS.
+# XXX: also reverse the orthogonality limits.
+Base.reverse(x::MPS) = MPS(reverse(x.data))