diff --git a/Project.toml b/Project.toml index 4d5ec8ca..1f89bd42 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ITensorNetworks" uuid = "2919e153-833c-4bdc-8836-1ea460a35fc7" authors = ["Matthew Fishman , Joseph Tindall and contributors"] -version = "0.13.1" +version = "0.13.2" [deps] AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c" @@ -25,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" @@ -62,7 +61,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" @@ -73,12 +72,11 @@ 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" diff --git a/examples/Project.toml b/examples/Project.toml index 7c06c098..f8109c92 100644 --- a/examples/Project.toml +++ b/examples/Project.toml @@ -2,4 +2,4 @@ ITensorNetworks = "2919e153-833c-4bdc-8836-1ea460a35fc7" [compat] -ITensorNetworks = "0.13.0" +ITensorNetworks = "0.13.2" diff --git a/ext/ITensorNetworksTensorOperationsExt/ITensorNetworksTensorOperationsExt.jl b/ext/ITensorNetworksTensorOperationsExt/ITensorNetworksTensorOperationsExt.jl index 8e39b4ee..242b0201 100644 --- a/ext/ITensorNetworksTensorOperationsExt/ITensorNetworksTensorOperationsExt.jl +++ b/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 diff --git a/src/apply.jl b/src/apply.jl index b4946210..5250c806 100644 --- a/src/apply.jl +++ b/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 diff --git a/src/caches/abstractbeliefpropagationcache.jl b/src/caches/abstractbeliefpropagationcache.jl index 63b7e43d..55e1e9fe 100644 --- a/src/caches/abstractbeliefpropagationcache.jl +++ b/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) diff --git a/src/caches/beliefpropagationcache.jl b/src/caches/beliefpropagationcache.jl index 1e820107..2c0a23a5 100644 --- a/src/caches/beliefpropagationcache.jl +++ b/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 diff --git a/src/expect.jl b/src/expect.jl index ac9752cf..c1e92efe 100644 --- a/src/expect.jl +++ b/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 diff --git a/test/test_gauging.jl b/test/test_gauging.jl index 67b1b5dd..b56f1b93 100644 --- a/test/test_gauging.jl +++ b/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