JoeyT1994 · JoeyT1994 · Jul 13, 2026 · Jul 8, 2026 · Jul 10, 2026
diff --git 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 <jtindall@flatironinstitute.org>", "MSRudolph <manuel.rudolph@web.de>", "and contributors"]
+authors = ["JoeyT1994 <jtindall@flatironinstitute.org>", "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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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