Merge remote-tracking branch 'upstream/main' into TensorNetworkFitting

Author: JoeyT1994

Project.toml

@@ -1,7 +1,7 @@
 name = "ITensorNetworks"
 uuid = "2919e153-833c-4bdc-8836-1ea460a35fc7"
 authors = ["Matthew Fishman <[email protected]>, Joseph Tindall <[email protected]> and contributors"]
-version = "0.12.2"
+version = "0.13.2"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
@@ -13,7 +13,6 @@ Dictionaries = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
-ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2"
 ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5"
 IsApprox = "28f27b66-4bd8-47e7-9110-e2746eb8bed7"
 IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
@@ -26,7 +25,6 @@ Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SerializedElementArrays = "d3ce8812-9567-47e9-a7b5-65a6d70a3065"
 SimpleTraits = "699a6c99-e7fa-54fc-8d76-47d257e15c1d"
-SparseArrayKit = "a9a3c162-d163-4c15-8926-b8794fbefed2"
 SplitApplyCombine = "03a91e81-4c3e-53e1-a0a4-9c0c8f19dd66"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 StructWalk = "31cdf514-beb7-4750-89db-dda9d2eb8d3d"
@@ -63,24 +61,22 @@ DocStringExtensions = "0.9"
 EinExprs = "0.6.4"
 Graphs = "1.8"
 GraphsFlows = "0.1.1"
-ITensorMPS = "0.3"
-ITensors = "0.7, 0.8"
+ITensors = "0.7, 0.8, 0.9"
 IsApprox = "0.1, 1, 2"
 IterTools = "1.4.0"
-KrylovKit = "0.6, 0.7, 0.8"
+KrylovKit = "0.6, 0.7, 0.8, 0.9"
 MacroTools = "0.5"
 NDTensors = "0.3, 0.4"
 NamedGraphs = "0.6.0"
 OMEinsumContractionOrders = "0.8.3, 0.9"
 Observers = "0.2.4"
 SerializedElementArrays = "0.1"
 SimpleTraits = "0.9"
-SparseArrayKit = "0.3, 0.4"
 SplitApplyCombine = "1.2"
 StaticArrays = "1.5.12"
 StructWalk = "0.2"
 Suppressor = "0.2"
-TensorOperations = "5.1.4"
+TensorOperations = "5.2.0"
 TimerOutputs = "0.5.22"
 TupleTools = "1.4"
 julia = "1.10"

docs/Project.toml

@@ -5,5 +5,5 @@ Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 
 [compat]
 Documenter = "1.10.0"
-ITensorNetworks = "0.12.0"
+ITensorNetworks = "0.13.0"
 Literate = "2.20.1"

examples/Project.toml

@@ -2,4 +2,4 @@
 ITensorNetworks = "2919e153-833c-4bdc-8836-1ea460a35fc7"
 
 [compat]
-ITensorNetworks = "0.12.0"
+ITensorNetworks = "0.13.2"

ext/ITensorNetworksTensorOperationsExt/ITensorNetworksTensorOperationsExt.jl

@@ -7,7 +7,8 @@ using TensorOperations: TensorOperations, optimaltree
 
 function ITensorNetworks.contraction_sequence(::Algorithm"optimal", tn::ITensorList)
   network = collect.(inds.(tn))
-  inds_to_dims = Dict(i => dim(i) for i in unique(reduce(vcat, network)))
+  #Converting dims to Float64 to minimize overflow issues
+  inds_to_dims = Dict(i => Float64(dim(i)) for i in unique(reduce(vcat, network)))
   seq, _ = optimaltree(network, inds_to_dims)
   return seq
 end

src/ITensorNetworks.jl

@@ -23,7 +23,6 @@ include("contract_approx/partition.jl")
 include("contract_approx/binary_tree_partition.jl")
 include("contract.jl")
 include("specialitensornetworks.jl")
-include("boundarymps.jl")
 include("partitioneditensornetwork.jl")
 include("edge_sequences.jl")
 include("caches/abstractbeliefpropagationcache.jl")

src/abstractitensornetwork.jl

@@ -38,7 +38,6 @@ using ITensors:
   settags,
   sim,
   swaptags
-using ITensorMPS: ITensorMPS, add, linkdim, linkinds, siteinds
 using .ITensorsExtensions: ITensorsExtensions, indtype, promote_indtype
 using LinearAlgebra: LinearAlgebra, factorize
 using MacroTools: @capture
@@ -255,7 +254,7 @@ indsnetwork(tn::AbstractITensorNetwork) = IndsNetwork(tn)
 
 # TODO: Output a `VertexDataGraph`? Unfortunately
 # `IndsNetwork` doesn't allow iterating over vertex data.
-function ITensorMPS.siteinds(tn::AbstractITensorNetwork)
+function siteinds(tn::AbstractITensorNetwork)
   is = IndsNetwork(underlying_graph(tn))
   for v in vertices(tn)
     is[v] = uniqueinds(tn, v)
@@ -268,7 +267,7 @@ function flatten_siteinds(tn::AbstractITensorNetwork)
   return identity.(flatten(map(v -> siteinds(tn, v), vertices(tn))))
 end
 
-function ITensorMPS.linkinds(tn::AbstractITensorNetwork)
+function linkinds(tn::AbstractITensorNetwork)
   is = IndsNetwork(underlying_graph(tn))
   for e in edges(tn)
     is[e] = commoninds(tn, e)
@@ -302,7 +301,11 @@ function ITensors.uniqueinds(tn::AbstractITensorNetwork, edge::Pair)
   return uniqueinds(tn, edgetype(tn)(edge))
 end
 
-function ITensors.siteinds(tn::AbstractITensorNetwork, vertex)
+function siteinds(tn::AbstractITensorNetwork, vertex)
+  return uniqueinds(tn, vertex)
+end
+# Fix ambiguity error with IndsNetwork constructor.
+function siteinds(tn::AbstractITensorNetwork, vertex::Int)
   return uniqueinds(tn, vertex)
 end
 
@@ -311,7 +314,7 @@ function ITensors.commoninds(tn::AbstractITensorNetwork, edge)
   return commoninds(tn[src(e)], tn[dst(e)])
 end
 
-function ITensorMPS.linkinds(tn::AbstractITensorNetwork, edge)
+function linkinds(tn::AbstractITensorNetwork, edge)
   return commoninds(tn, edge)
 end
 
@@ -807,24 +810,24 @@ end
 # Link dimensions
 # 
 
-function ITensorMPS.maxlinkdim(tn::AbstractITensorNetwork)
+function maxlinkdim(tn::AbstractITensorNetwork)
   md = 1
   for e in edges(tn)
     md = max(md, linkdim(tn, e))
   end
   return md
 end
 
-function ITensorMPS.linkdim(tn::AbstractITensorNetwork, edge::Pair)
+function linkdim(tn::AbstractITensorNetwork, edge::Pair)
   return linkdim(tn, edgetype(tn)(edge))
 end
 
-function ITensorMPS.linkdim(tn::AbstractITensorNetwork{V}, edge::AbstractEdge{V}) where {V}
+function linkdim(tn::AbstractITensorNetwork{V}, edge::AbstractEdge{V}) where {V}
   ls = linkinds(tn, edge)
   return prod([isnothing(l) ? 1 : dim(l) for l in ls])
 end
 
-function ITensorMPS.linkdims(tn::AbstractITensorNetwork{V}) where {V}
+function linkdims(tn::AbstractITensorNetwork{V}) where {V}
   ld = DataGraph{V}(
     copy(underlying_graph(tn)); vertex_data_eltype=Nothing, edge_data_eltype=Int
   )
@@ -882,7 +885,7 @@ is_multi_edge(tn::AbstractITensorNetwork, e) = length(linkinds(tn, e)) > 1
 is_multi_edge(tn::AbstractITensorNetwork) = Base.Fix1(is_multi_edge, tn)
 
 """Add two itensornetworks together by growing the bond dimension. The network structures need to be have the same vertex names, same site index on each vertex """
-function ITensorMPS.add(tn1::AbstractITensorNetwork, tn2::AbstractITensorNetwork)
+function add(tn1::AbstractITensorNetwork, tn2::AbstractITensorNetwork)
   @assert issetequal(vertices(tn1), vertices(tn2))
 
   tn1 = combine_linkinds(tn1; edges=filter(is_multi_edge(tn1), edges(tn1)))

src/apply.jl

@@ -22,8 +22,6 @@ using ITensors:
   replaceinds,
   unioninds,
   uniqueinds
-using ITensors.ContractionSequenceOptimization: optimal_contraction_sequence
-using ITensorMPS: siteinds
 using KrylovKit: linsolve
 using LinearAlgebra: eigen, norm, svd
 using NamedGraphs: NamedEdge, has_edge
@@ -407,7 +405,8 @@ function fidelity(
     ],
     envs,
   )
-  term1 = ITensors.contract(term1_tns; sequence=optimal_contraction_sequence(term1_tns))
+  sequence = contraction_sequence(term1_tns; alg="optimal")
+  term1 = ITensors.contract(term1_tns; sequence)
 
   term2_tns = vcat(
     [
@@ -418,9 +417,11 @@ function fidelity(
     ],
     envs,
   )
-  term2 = ITensors.contract(term2_tns; sequence=optimal_contraction_sequence(term2_tns))
+  sequence = contraction_sequence(term2_tns; alg="optimal")
+  term2 = ITensors.contract(term2_tns; sequence)
   term3_tns = vcat([p_prev, q_prev, prime(dag(p_cur)), prime(dag(q_cur)), gate], envs)
-  term3 = ITensors.contract(term3_tns; sequence=optimal_contraction_sequence(term3_tns))
+  sequence = contraction_sequence(term3_tns; alg="optimal")
+  term3 = ITensors.contract(term3_tns; sequence)
 
   f = term3[] / sqrt(term1[] * term2[])
   return f * conj(f)
@@ -447,61 +448,32 @@ function optimise_p_q(
   qs_ind = setdiff(inds(q_cur), collect(Iterators.flatten(inds.(vcat(envs, p_cur)))))
   ps_ind = setdiff(inds(p_cur), collect(Iterators.flatten(inds.(vcat(envs, q_cur)))))
 
-  opt_b_seq = optimal_contraction_sequence(vcat(ITensor[p, q, o, dag(prime(q_cur))], envs))
-  opt_b_tilde_seq = optimal_contraction_sequence(
-    vcat(ITensor[p, q, o, dag(prime(p_cur))], envs)
-  )
-  opt_M_seq = optimal_contraction_sequence(
-    vcat(ITensor[q_cur, replaceinds(prime(dag(q_cur)), prime(qs_ind), qs_ind), p_cur], envs)
-  )
-  opt_M_tilde_seq = optimal_contraction_sequence(
-    vcat(ITensor[p_cur, replaceinds(prime(dag(p_cur)), prime(ps_ind), ps_ind), q_cur], envs)
-  )
-
-  function b(
-    p::ITensor,
-    q::ITensor,
-    o::ITensor,
-    envs::Vector{ITensor},
-    r::ITensor;
-    opt_sequence=nothing,
-  )
-    return noprime(
-      ITensors.contract(vcat(ITensor[p, q, o, dag(prime(r))], envs); sequence=opt_sequence)
-    )
+  function b(p::ITensor, q::ITensor, o::ITensor, envs::Vector{ITensor}, r::ITensor)
+    ts = vcat(ITensor[p, q, o, dag(prime(r))], envs)
+    sequence = contraction_sequence(ts; alg="optimal")
+    return noprime(ITensors.contract(ts; sequence))
   end
 
-  function M_p(
-    envs::Vector{ITensor},
-    p_q_tensor::ITensor,
-    s_ind,
-    apply_tensor::ITensor;
-    opt_sequence=nothing,
-  )
-    return noprime(
-      ITensors.contract(
-        vcat(
-          ITensor[
-            p_q_tensor,
-            replaceinds(prime(dag(p_q_tensor)), prime(s_ind), s_ind),
-            apply_tensor,
-          ],
-          envs,
-        );
-        sequence=opt_sequence,
-      ),
+  function M_p(envs::Vector{ITensor}, p_q_tensor::ITensor, s_ind, apply_tensor::ITensor)
+    ts = vcat(
+      ITensor[
+        p_q_tensor, replaceinds(prime(dag(p_q_tensor)), prime(s_ind), s_ind), apply_tensor
+      ],
+      envs,
     )
+    sequence = contraction_sequence(ts; alg="optimal")
+    return noprime(ITensors.contract(ts; sequence))
   end
   for i in 1:nfullupdatesweeps
-    b_vec = b(p, q, o, envs, q_cur; opt_sequence=opt_b_seq)
-    M_p_partial = partial(M_p, envs, q_cur, qs_ind; opt_sequence=opt_M_seq)
+    b_vec = b(p, q, o, envs, q_cur)
+    M_p_partial = partial(M_p, envs, q_cur, qs_ind)
 
     p_cur, info = linsolve(
       M_p_partial, b_vec, p_cur; isposdef=envisposdef, ishermitian=false
     )
 
-    b_tilde_vec = b(p, q, o, envs, p_cur; opt_sequence=opt_b_tilde_seq)
-    M_p_tilde_partial = partial(M_p, envs, p_cur, ps_ind; opt_sequence=opt_M_tilde_seq)
+    b_tilde_vec = b(p, q, o, envs, p_cur)
+    M_p_tilde_partial = partial(M_p, envs, p_cur, ps_ind)
 
     q_cur, info = linsolve(
       M_p_tilde_partial, b_tilde_vec, q_cur; isposdef=envisposdef, ishermitian=false

src/boundarymps.jl

@@ -2,7 +2,6 @@ using Graphs: IsDirected
 using SplitApplyCombine: group
 using LinearAlgebra: diag, dot
 using ITensors: dir
-using ITensorMPS: ITensorMPS
 using NamedGraphs.PartitionedGraphs:
   PartitionedGraphs,
   PartitionedGraph,
@@ -17,7 +16,7 @@ using NDTensors: NDTensors
 abstract type AbstractBeliefPropagationCache end
 
 function default_message_update(contract_list::Vector{ITensor}; normalize=true, kwargs...)
-  sequence = optimal_contraction_sequence(contract_list)
+  sequence = contraction_sequence(contract_list; alg="optimal")
   updated_messages = contract(contract_list; sequence, kwargs...)
   message_norm = norm(updated_messages)
   if normalize && !iszero(message_norm)
@@ -140,7 +139,7 @@ for f in [
   :(PartitionedGraphs.partitionvertices),
   :(PartitionedGraphs.vertices),
   :(PartitionedGraphs.boundary_partitionedges),
-  :(ITensorMPS.linkinds),
+  :(linkinds),
 ]
   @eval begin
     function $f(bpc::AbstractBeliefPropagationCache, args...; kwargs...)

src/caches/abstractbeliefpropagationcache.jl

@@ -2,7 +2,6 @@ using Graphs: IsDirected
 using SplitApplyCombine: group
 using LinearAlgebra: diag, dot
 using ITensors: dir
-using ITensorMPS: ITensorMPS
 using NamedGraphs.PartitionedGraphs:
   PartitionedGraphs,
   PartitionedGraph,
@@ -94,22 +93,16 @@ function environment(bpc::BeliefPropagationCache, verts::Vector; kwargs...)
   return vcat(messages, central_tensors)
 end
 
-function region_scalar(
-  bp_cache::BeliefPropagationCache,
-  pv::PartitionVertex;
-  contract_kwargs=(; sequence="automatic"),
-)
+function region_scalar(bp_cache::BeliefPropagationCache, pv::PartitionVertex)
   incoming_mts = incoming_messages(bp_cache, [pv])
   local_state = factors(bp_cache, pv)
-  return contract(vcat(incoming_mts, local_state); contract_kwargs...)[]
+  ts = vcat(incoming_mts, local_state)
+  sequence = contraction_sequence(ts; alg="optimal")
+  return contract(ts; sequence)[]
 end
 
-function region_scalar(
-  bp_cache::BeliefPropagationCache,
-  pe::PartitionEdge;
-  contract_kwargs=(; sequence="automatic"),
-)
-  return contract(
-    vcat(message(bp_cache, pe), message(bp_cache, reverse(pe))); contract_kwargs...
-  )[]
+function region_scalar(bp_cache::BeliefPropagationCache, pe::PartitionEdge)
+  ts = vcat(message(bp_cache, pe), message(bp_cache, reverse(pe)))
+  sequence = contraction_sequence(ts; alg="optimal")
+  return contract(ts; sequence)[]
 end