revert broken TaylorCluster optimizations

src/MPSKit.jl

@@ -10,7 +10,6 @@ using RecipesBase
 using VectorInterface
 using Accessors
 using HalfIntegers
-import TupleTools as TT
 
 using LinearAlgebra: diag, Diagonal
 using LinearAlgebra: LinearAlgebra

src/algorithms/timestep/timeevmpo.jl

@@ -56,20 +56,23 @@ function make_time_mpo(H::MPOHamiltonian, dt::Number, alg::TaylorCluster)
     if alg.extension
         H_next = H_n * H′
         linds_next = LinearIndices(ntuple(i -> V, N + 1))
-        cinds_next = CartesianIndices(linds_next)
         for (i, slice) in enumerate(parent(H_n))
-            for Inext in nonzero_keys(H_next[i])
-                aₑ = cinds_next[Inext[1]]
-                bₑ = cinds_next[Inext[4]]
-                for c in findall(==(1), aₑ.I), d in findall(==(V), bₑ.I)
-                    a = TT.deleteat(Tuple(aₑ), c)
-                    b = TT.deleteat(Tuple(bₑ), d)
-                    if all(>(1), b) && !(all(in((1, V), a)) && any(==(V), a))
-                        n1 = count(==(1), a) + 1
-                        n3 = count(==(V), b) + 1
-                        factor = τ * factorial(N) / (factorial(N + 1) * n1 * n3)
-                        slice[linds[a...], 1, 1, linds[b...]] += factor * H_next[i][Inext]
-                    end
+            for a in cinds, b in cinds
+                all(>(1), b.I) || continue
+                all(in((1, V)), a.I) && any(==(V), a.I) && continue
+
+                n1 = count(==(1), a.I) + 1
+                n3 = count(==(V), b.I) + 1
+                factor = τ * factorial(N) / (factorial(N + 1) * n1 * n3)
+
+                for c in 1:(N + 1), d in 1:(N + 1)
+                    aₑ = insert!([a.I...], c, 1)
+                    bₑ = insert!([b.I...], d, V)
+
+                    # TODO: use VectorInterface for memory efficiency
+                    slice[linds[a], 1, 1, linds[b]] += factor *
+                                                       H_next[i][linds_next[aₑ...], 1, 1,
+                                                                 linds_next[bₑ...]]
                 end
             end
         end
@@ -83,46 +86,39 @@ function make_time_mpo(H::MPOHamiltonian, dt::Number, alg::TaylorCluster)
         mpo = InfiniteMPO(parent(H_n))
     end
     for slice in parent(mpo)
-        for (I, v) in nonzero_pairs(slice)
-            a = cinds[I[1]]
-            b = cinds[I[4]]
-            if all(in((1, V)), a.I) && any(!=(1), a.I)
-                delete!(slice, I)
-            elseif any(!=(1), b.I) && all(in((1, V)), b.I)
-                exponent = count(==(V), b.I)
-                factor = τ^exponent * factorial(N - exponent) / factorial(N)
-                Idst = CartesianIndex(Base.setindex(I.I, 1, 4))
-                slice[Idst] += factor * v
-                delete!(slice, I)
+        for b in cinds[2:end]
+            all(in((1, V)), b.I) || continue
+
+            b_lin = linds[b]
+            a = count(==(V), b.I)
+            factor = τ^a * factorial(N - a) / factorial(N)
+            slice[:, 1, 1, 1] = slice[:, 1, 1, 1] + factor * slice[:, 1, 1, b_lin]
+            for I in nonzero_keys(slice)
+                (I[1] == b_lin || I[4] == b_lin) && delete!(slice, I)
             end
         end
     end
 
     # Remove equivalent rows and columns: Algorithm 2
     for slice in parent(mpo)
-        for I in nonzero_keys(slice)
-            a = cinds[I[1]]
-            a_sort = CartesianIndex(sort(collect(a.I); by=(!=(1)))...)
-            n1_a = count(==(1), a.I)
-            n3_a = count(==(V), a.I)
-            if n3_a <= n1_a && a_sort != a
-                Idst = CartesianIndex(Base.setindex(I.I, linds[a_sort], 1))
-                slice[Idst] += slice[I]
-                delete!(slice, I)
-            elseif a != a_sort
-                delete!(slice, I)
-            end
+        for c in cinds
+            c_lin = linds[c]
+            s_c = CartesianIndex(sort(collect(c.I); by=(!=(1)))...)
+            s_r = CartesianIndex(sort(collect(c.I); by=(!=(V)))...)
+
+            n1 = count(==(1), c.I)
+            n3 = count(==(V), c.I)
 
-            b = cinds[I[4]]
-            b_sort = CartesianIndex(sort(collect(b.I); by=(!=(V)))...)
-            n1_b = count(==(1), b.I)
-            n3_b = count(==(V), b.I)
-            if n3_b > n1_b && b_sort != b
-                Idst = CartesianIndex(Base.setindex(I.I, linds[b_sort], 4))
-                slice[Idst] += slice[I]
-                delete!(slice, I)
-            elseif b != b_sort
-                delete!(slice, I)
+            if n3 <= n1 && s_c != c
+                slice[linds[s_c], 1, 1, :] += slice[c_lin, 1, 1, :]
+                for I in nonzero_keys(slice)
+                    (I[1] == c_lin || I[4] == c_lin) && delete!(slice, I)
+                end
+            elseif n3 > n1 && s_r != c
+                slice[:, 1, 1, linds[s_r]] += slice[:, 1, 1, c_lin]
+                for I in nonzero_keys(slice)
+                    (I[1] == c_lin || I[4] == c_lin) && delete!(slice, I)
+                end
             end
         end
     end
@@ -137,14 +133,10 @@ function make_time_mpo(H::MPOHamiltonian, dt::Number, alg::TaylorCluster)
                 b = CartesianIndex(replace(a.I, V => 1))
                 b_lin = linds[b]
                 factor = τ^n1 * factorial(N - n1) / factorial(N)
+                slice[:, 1, 1, b_lin] += factor * slice[:, 1, 1, a_lin]
+
                 for I in nonzero_keys(slice)
-                    if I[4] == a_lin
-                        Idst = CartesianIndex(Base.setindex(I.I, b_lin, 4))
-                        slice[Idst] += factor * slice[I]
-                        delete!(slice, I)
-                    elseif I[1] == a_lin
-                        delete!(slice, I)
-                    end
+                    (I[1] == a_lin || I[4] == a_lin) && delete!(slice, I)
                 end
             end
         end