Make compatible with ITensors 0.9 (#221)

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.13.1"
+version = "0.13.2"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
@@ -62,7 +62,7 @@ DocStringExtensions = "0.9"
 EinExprs = "0.6.4"
 Graphs = "1.8"
 GraphsFlows = "0.1.1"
-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, 0.9"
@@ -78,7 +78,7 @@ 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"

examples/Project.toml

@@ -2,4 +2,4 @@
 ITensorNetworks = "2919e153-833c-4bdc-8836-1ea460a35fc7"
 
 [compat]
-ITensorNetworks = "0.13.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/apply.jl

@@ -22,7 +22,6 @@ using ITensors:
   replaceinds,
   unioninds,
   uniqueinds
-using ITensors.ContractionSequenceOptimization: optimal_contraction_sequence
 using KrylovKit: linsolve
 using LinearAlgebra: eigen, norm, svd
 using NamedGraphs: NamedEdge, has_edge
@@ -406,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(
     [
@@ -417,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)
@@ -446,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/caches/abstractbeliefpropagationcache.jl

@@ -16,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)

src/caches/beliefpropagationcache.jl

@@ -93,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

src/expect.jl

@@ -3,16 +3,18 @@ using ITensors: Op, op, contract, which_op
 
 default_expect_alg() = "bp"
 
-function expect(
-  ψIψ::AbstractFormNetwork, op::Op; contract_kwargs=(; sequence="automatic"), kwargs...
-)
+function expect(ψIψ::AbstractFormNetwork, op::Op; kwargs...)
   v = only(op.sites)
   ψIψ_v = ψIψ[operator_vertex(ψIψ, v)]
   s = commonind(ψIψ[ket_vertex(ψIψ, v)], ψIψ_v)
   operator = ITensors.op(op.which_op, s)
   ∂ψIψ_∂v = environment(ψIψ, operator_vertices(ψIψ, [v]); kwargs...)
-  numerator = contract(vcat(∂ψIψ_∂v, operator); contract_kwargs...)[]
-  denominator = contract(vcat(∂ψIψ_∂v, ψIψ_v); contract_kwargs...)[]
+  numerator_ts = vcat(∂ψIψ_∂v, operator)
+  denominator_ts = vcat(∂ψIψ_∂v, ψIψ_v)
+  numerator_seq = contraction_sequence(numerator_ts; alg="optimal")
+  denominator_seq = contraction_sequence(denominator_ts; alg="optimal")
+  numerator = contract(numerator_ts; sequence=numerator_seq)[]
+  denominator = contract(denominator_ts; sequence=denominator_seq)[]
 
   return numerator / denominator
 end

test/test_gauging.jl

@@ -14,6 +14,7 @@ using ITensors.NDTensors: vector
 using LinearAlgebra: diag
 using NamedGraphs.NamedGraphGenerators: named_grid
 using StableRNGs: StableRNG
+using TensorOperations: TensorOperations
 using Test: @test, @testset
 
 @testset "gauging" begin