diff --git a/Project.toml b/Project.toml index ff8fc4a9..cf796d01 100644 --- a/Project.toml +++ b/Project.toml @@ -2,7 +2,7 @@ name = "TensorNetworkQuantumSimulator" uuid = "4de3b72a-362e-43dd-83ff-3f381eda9f9c" license = "MIT" version = "0.4.1" -authors = ["JoeyT1994 ", "MSRudolph ", "and contributors"] +authors = ["JoeyT1994 ", "and contributors"] description = "A Julia package for quantum simulation with tensor networks of near-arbitrary topology." [deps] @@ -20,7 +20,6 @@ SimpleGraphConverter = "205b04f2-f585-4877-a239-566270b3f673" SplitApplyCombine = "03a91e81-4c3e-53e1-a0a4-9c0c8f19dd66" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" -TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2" TypeParameterAccessors = "7e5a90cf-f82e-492e-a09b-e3e26432c138" [compat] @@ -31,14 +30,13 @@ ITensors = "0.9" KrylovKit = "0.10.2" LinearAlgebra = "1.11.0" NamedGraphs = "0.12.1" -OMEinsumContractionOrders = "1.2" +OMEinsumContractionOrders = "1.3" Revise = "3.8.0" SimpleGraphAlgorithms = "0.6.0" SimpleGraphConverter = "0.1.0" SplitApplyCombine = "1.2.3" Statistics = "1.11.1" StatsBase = "0.34.4" -TensorOperations = "5.2" TypeParameterAccessors = "0.3.10, 0.4" julia = "1.10" diff --git a/src/TensorNetworkQuantumSimulator.jl b/src/TensorNetworkQuantumSimulator.jl index a2273a07..7cf72a81 100644 --- a/src/TensorNetworkQuantumSimulator.jl +++ b/src/TensorNetworkQuantumSimulator.jl @@ -97,6 +97,7 @@ export Treewidth, ExactTreewidth, HyperND, + ExhaustiveSearch, norm_sqr, map_virtualinds, map_virtualinds!, diff --git a/src/contraction_sequences.jl b/src/contraction_sequences.jl index 83455e12..d4ee7b6e 100644 --- a/src/contraction_sequences.jl +++ b/src/contraction_sequences.jl @@ -1,35 +1,20 @@ using ITensors: Index, ITensor, @Algorithm_str, inds, noncommoninds, dim -using TensorOperations: TensorOperations, optimaltree -using OMEinsumContractionOrders: OMEinsumContractionOrders, optimize_code, EinCode, NestedEinsum, TreeSA, GreedyMethod, SABipartite, Treewidth, ExactTreewidth, HyperND +using OMEinsumContractionOrders: OMEinsumContractionOrders, optimize_code, EinCode, NestedEinsum, TreeSA, GreedyMethod, SABipartite, Treewidth, ExactTreewidth, HyperND, ExhaustiveSearch -function prune_trivial_tensors(tensors::Vector{<:ITensor}) - pruned_tensors = copy(tensors) - for (i, t) in enumerate(pruned_tensors) - if all(d -> d == 1, dim.(inds(tensors[i]))) - pruned_tensors[i] = adapt_like(t, ITensor(1)) - end - end - return pruned_tensors -end - -function contraction_sequence(::Algorithm"optimal", tensors::Vector{<:ITensor}; prune_tensors = false) - #Needed because tensor operations bugs on trivial tensors - if prune_tensors - ITensors.disable_warn_order() - tensors = prune_trivial_tensors(tensors) - end - network = collect.(inds.(tensors)) - #Converting dims to Float64 to minimize overflow issues - inds_to_dims = Dict(i => Float64(dim(i)) for i in unique(Iterators.flatten(network))) - seq, _ = optimaltree(network, inds_to_dims) - seq = typeof(seq) <: Int ? [seq] : seq - return seq +# The exact "optimal" contraction order (Pfeifer 2014 netcon) is provided by +# OMEinsumContractionOrders' `ExhaustiveSearch` optimizer, which ported this routine from +# TensorOperations. It handles trivial 1-/2-tensor inputs directly, so the previous +# trivial-tensor pruning and scalar-`Int` workarounds are no longer needed. +function contraction_sequence(::Algorithm"optimal", tensors::Vector{<:ITensor}) + return contraction_sequence(Algorithm("omeinsum"), tensors; optimizer = ExhaustiveSearch()) end function contraction_sequence(::Algorithm"omeinsum", tensors::Vector{<:ITensor}; optimizer = TreeSA()) code, size_dict = to_eincode(tensors) optcode = optimize_code(code, size_dict, optimizer) - return to_contraction_sequence(optcode) + seq = to_contraction_sequence(optcode) + #A single-tensor network optimizes to a lone leaf; wrap it so a Vector is always returned. + return seq isa Integer ? [seq] : seq end function contraction_sequence(tensors::Vector{<:ITensor}; alg = "optimal", kwargs...) @@ -40,7 +25,9 @@ end function to_eincode(tensors::Vector{<:ITensor}) ixs = map(t -> collect(inds(t)), tensors) LT = eltype(eltype(ixs)) - iy = collect(LT, reduce(noncommoninds, tensors)) + #`reduce` over a single tensor returns the tensor itself, not its indices; a one-tensor + #network is trivial and its open indices are all of that tensor's indices. + iy = length(tensors) == 1 ? collect(LT, inds(only(tensors))) : collect(LT, reduce(noncommoninds, tensors)) size_dict = Dict{LT, Int}(i => dim(i) for ix in ixs for i in ix) return EinCode(ixs, iy), size_dict end diff --git a/src/expect.jl b/src/expect.jl index 003f9f72..cd81eaf4 100644 --- a/src/expect.jl +++ b/src/expect.jl @@ -71,7 +71,7 @@ function expect( function contract_region(op_string_f) tensors = norm_factors(network(cache), steiner_vs; op_strings = op_string_f) append!(tensors, incoming_ms) - seq = contraction_sequence(tensors; alg = "optimal", prune_tensors = true) + seq = contraction_sequence(tensors; alg = "optimal") return scalar(contract(tensors; sequence = seq)) end diff --git a/src/imports.jl b/src/imports.jl index f1b43234..febbbf68 100644 --- a/src/imports.jl +++ b/src/imports.jl @@ -44,8 +44,6 @@ using NamedGraphs.GraphsExtensions: using NamedGraphs.NamedGraphGenerators: named_grid, named_hexagonal_lattice_graph, named_comb_tree, named_path_graph -using TensorOperations - using ITensors: ITensors using ITensors: Index, ITensor, hasqns, noncommonind, combinedind, combiner, replaceinds, sim, onehot, delta, plev, dense, unioninds, uniqueinds, commonind, commoninds, replaceind, datatype, inds, dag, noprime, factorize_svd, prime, hascommoninds, itensor, map_diag!, @Algorithm_str, scalar, @OpName_str, @SiteType_str, denseblocks, tags, op, apply, contract, inner diff --git a/src/rdm.jl b/src/rdm.jl index bf942e00..6e6c6bd6 100644 --- a/src/rdm.jl +++ b/src/rdm.jl @@ -63,7 +63,7 @@ function reduced_density_matrix( #TODO: If there are a lot of tensors here, (more than 100 say), we need to think about defining a custom sequence as optimal may be too slow ρ_tensors = norm_factors(network(cache), steiner_vs; op_strings = op_string_f) append!(ρ_tensors, incoming_ms) - seq = contraction_sequence(ρ_tensors; alg = "optimal", prune_tensors = true) + seq = contraction_sequence(ρ_tensors; alg = "optimal") ρ = contract(ρ_tensors; sequence = seq) if normalize diff --git a/test/test_contraction_sequences.jl b/test/test_contraction_sequences.jl index 99e50b28..347c2169 100644 --- a/test/test_contraction_sequences.jl +++ b/test/test_contraction_sequences.jl @@ -37,7 +37,7 @@ collect_leaves!(acc, x) = (for y in x; collect_leaves!(acc, y); end; acc) tn = random_tensornetwork(Float64, g; bond_dimension = 2) tensors = [tn[v] for v in vertices(tn)] n = length(tensors) - for optimizer in (GreedyMethod(), TreeSA()) + for optimizer in (GreedyMethod(), TreeSA(), ExhaustiveSearch()) seq = TNQS.contraction_sequence(tensors; alg = "omeinsum", optimizer) @test sort(collect_leaves!(Int[], seq)) == collect(1:n) # every tensor exactly once @test seq isa AbstractVector && any(x -> x isa AbstractVector, seq) # nested tree, not a flat list @@ -46,7 +46,7 @@ collect_leaves!(acc, x) = (for y in x; collect_leaves!(acc, y); end; acc) # --- the sequence the backend returns is a *correct* contraction: executing it gives the # same scalar as the independent `optimal` backend. ref = scalar(ITensors.contract(tensors; sequence = TNQS.contraction_sequence(tensors; alg = "optimal"))) - for optimizer in (GreedyMethod(), TreeSA()) + for optimizer in (GreedyMethod(), TreeSA(), ExhaustiveSearch()) seq = TNQS.contraction_sequence(tensors; alg = "omeinsum", optimizer) @test scalar(ITensors.contract(tensors; sequence = seq)) ≈ ref end