From 5b496a505ebe537ade0cc9ad3fe097e2e53276df Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Tue, 18 May 2021 15:52:48 -0400
Subject: [PATCH 01/12] Start optimizing swap gates in apply

---
 docs/src/MPSandMPO.md  |  1 +
 src/mps/abstractmps.jl | 51 +++++++++++++++++++++++++++++++-----------
 2 files changed, 39 insertions(+), 13 deletions(-)

diff --git a/docs/src/MPSandMPO.md b/docs/src/MPSandMPO.md
index 5031887e70..5b011bb9bd 100644
--- a/docs/src/MPSandMPO.md
+++ b/docs/src/MPSandMPO.md
@@ -33,6 +33,7 @@ MPO(::Type{<:Number}, ::Vector{<:Index}, ::String)
 
 ```@docs
 eltype(::ITensors.AbstractMPS)
+ITensors.promote_itensor_eltype(::ITensors.AbstractMPS)
 flux(::ITensors.AbstractMPS)
 hasqns(::ITensors.AbstractMPS)
 length(::ITensors.AbstractMPS)
diff --git a/src/mps/abstractmps.jl b/src/mps/abstractmps.jl
index ad0c1860f0..ecd8d2e9cf 100644
--- a/src/mps/abstractmps.jl
+++ b/src/mps/abstractmps.jl
@@ -42,7 +42,7 @@ end
 The element type of the MPS/MPO. Always returns `ITensor`.
 
 For the element type of the ITensors of the MPS/MPO,
-use `common_itensor_eltype`.
+use [`ITensors.promote_itensor_eltype`](@ref).
 """
 eltype(::AbstractMPS) = ITensor
 
@@ -1624,7 +1624,7 @@ function movesites(ψ::AbstractMPS, ns, ns′; kwargs...)
 end
 
 """
-    product(o::ITensor, ψ::Union{MPS, MPO}, [ns::Vector{Int}]; <keyword argument>)
+    product(o::ITensor, ψ::Union{MPS,MPO}, [ns::Vector{Int}]; <keyword argument>)
     apply([...])
 
 Get the product of the operator `o` with the MPS/MPO `ψ`,
@@ -1647,28 +1647,51 @@ function product(
   ns=findsites(ψ, o);
   move_sites_back::Bool=true,
   apply_dag::Bool=false,
+  next_gate::Union{ITensor,Nothing}=nothing,
   kwargs...,
 )
   N = length(ns)
   ns = sort(ns)
 
-  # TODO: make this smarter by minimizing
-  # distance to orthogonalization.
-  # For example, if ITensors.orthocenter(ψ) > ns[end],
-  # set to ns[end].
-  ψ = orthogonalize(ψ, ns[1])
+  next_ns = isnothing(next_gate) ? nothing : sort(findsites(ψ, next_gate))
+
+  @show ns
+  @show next_ns
+
+  if !isnothing(next_ns)
+    if first(next_ns) > last(ns)
+      @show first(next_ns) > last(ns)
+    end
+  end
+
   diff_ns = diff(ns)
   ns′ = ns
   if any(!=(1), diff_ns)
     ns′ = [ns[1] + n - 1 for n in 1:N]
     ψ = movesites(ψ, ns .=> ns′; kwargs...)
   end
-  ϕ = ψ[ns′[1]]
-  for n in 2:N
-    ϕ *= ψ[ns′[n]]
+  ns_range = first(ns′):last(ns′)
+
+  @show ns_range
+  @show ortho_lims(ψ)
+  # Minimize distance to orthogonalize
+  if first(ortho_lims(ψ)) > last(ns_range)
+    ψ = orthogonalize(ψ, last(ns_range))
+  elseif last(ortho_lims(ψ)) < first(ns_range)
+    ψ = orthogonalize(ψ, first(ns_range))
   end
+
+  @show ns_range
+  @show ortho_lims(ψ)
+
+  if first(ortho_lims(ψ)) < first(ns_range) || last(ortho_lims(ψ)) > last(ns_range)
+    error("Orthogonality center limits are $(ortho_lims(ψ)), but gate is being applied to $ns′. Must be within the range of where the gates are being applied.")
+  end
+
+  ϕ = prod(ψ[ns_range])
   ϕ = product(o, ϕ; apply_dag=apply_dag)
-  ψ[ns′[1]:ns′[end], kwargs...] = ϕ
+
+  ψ[ns_range, kwargs...] = ϕ
   if move_sites_back
     # Move the sites back to their original positions
     ψ = movesites(ψ, ns′ .=> ns; kwargs...)
@@ -1763,8 +1786,10 @@ function product(
   As::Vector{<:ITensor}, ψ::AbstractMPS; move_sites_back::Bool=true, kwargs...
 )
   Aψ = ψ
-  for A in As
-    Aψ = product(A, Aψ; move_sites_back=false, kwargs...)
+  for n in eachindex(As)
+    An = As[n]
+    next_gate = n == lastindex(As) ? nothing : As[n + 1]
+    Aψ = product(An, Aψ; move_sites_back=false, next_gate=next_gate, kwargs...)
   end
   if move_sites_back
     s = siteinds(Aψ)

From 8428f213754755ac512fe36279935ca3148fe7dd Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Tue, 18 May 2021 16:03:05 -0400
Subject: [PATCH 02/12] Update randomMPS syntax to use linkdims keyword
 argument in tests

---
 docs/src/MPSandMPO.md                 |  8 ++--
 docs/src/examples/DMRG.md             |  2 +-
 docs/src/getting_started/Tutorials.md | 14 +++----
 src/mps/abstractmps.jl                |  2 +-
 test/autompo.jl                       | 28 ++++++-------
 test/dmrg.jl                          | 18 ++++----
 test/mps.jl                           | 60 +++++++++++++--------------
 7 files changed, 66 insertions(+), 66 deletions(-)

diff --git a/docs/src/MPSandMPO.md b/docs/src/MPSandMPO.md
index d714612f11..cad2e747f9 100644
--- a/docs/src/MPSandMPO.md
+++ b/docs/src/MPSandMPO.md
@@ -15,10 +15,10 @@ MPS(::Type{<:Number}, ::Vector{<:Index})
 randomMPS(sites::Vector{<:Index})
 randomMPS(::Type{<:Number}, sites::Vector{<:Index})
 randomMPS(::Vector{<:Index}, ::Any)
-productMPS(::Vector{<:Index}, ::Any)
-productMPS(::Type{<:Number}, ::Vector{<:Index}, ::Any)
-productMPS(::Vector{<:IndexVal})
-productMPS(::Type{<:Number}, ::Vector{<:IndexVal})
+MPS(::Vector{<:Index}, ::Any)
+MPS(::Type{<:Number}, ::Vector{<:Index}, ::Any)
+MPS(::Vector{<:IndexVal})
+MPS(::Type{<:Number}, ::Vector{<:IndexVal})
 ```
 
 ## MPO Constructors
diff --git a/docs/src/examples/DMRG.md b/docs/src/examples/DMRG.md
index 90840f159d..ae5476e9d1 100644
--- a/docs/src/examples/DMRG.md
+++ b/docs/src/examples/DMRG.md
@@ -52,7 +52,7 @@ psi0 = randomMPS(sites,2)
 Here we have made a random MPS of bond dimension 2. We could have used a random product
 state instead, but choosing a slightly larger bond dimension can help DMRG avoid getting
 stuck in local minima. We could also set psi to some specific initial state using the 
-function `productMPS`, which is actually required if we were conserving QNs.
+function `MPS`, which is actually required if we were conserving QNs.
 
 Finally, we are ready to call DMRG:
 
diff --git a/docs/src/getting_started/Tutorials.md b/docs/src/getting_started/Tutorials.md
index db7efbe784..5f8b6aa015 100644
--- a/docs/src/getting_started/Tutorials.md
+++ b/docs/src/getting_started/Tutorials.md
@@ -216,7 +216,7 @@ by the lines
 
 ```julia
 state = [isodd(n) ? "Up" : "Dn" for n=1:N]
-psi0 = productMPS(sites,state)
+psi0 = MPS(sites,state)
 ```
 
 The first line of the new code above makes an array of strings which 
@@ -251,13 +251,13 @@ computing the quantum-number flux of `psi0`
     use the initial state:
     ```julia
     state = ["Up" for n=1:N]
-    psi0 = productMPS(sites,state)
+    psi0 = MPS(sites,state)
     ```
     Or to initialize this 10-site system to have a total "Sz" of +16
     in ITensor units (``S^z=8`` in physical units):
     ```julia
     state = ["Dn","Up","Up","Up","Up","Up","Up","Up","Up","Up"]
-    psi0 = productMPS(sites,state)
+    psi0 = MPS(sites,state)
     ```
     would work (as would any `state` with one "Dn" and nine "Up"'s
     in any order).
@@ -265,7 +265,7 @@ computing the quantum-number flux of `psi0`
     in ITensor units (``S^z=9`` in physical units) as
     ```julia
     state = ["Z0","Up","Up","Up","Up","Up","Up","Up","Up","Up"]
-    psi0 = productMPS(sites,state)
+    psi0 = MPS(sites,state)
     ```
     where "Z0" refers to the ``S^z=0`` state of a spin-one spin.
 
@@ -294,7 +294,7 @@ let
   H = MPO(ampo,sites)
 
   state = [isodd(n) ? "Up" : "Dn" for n=1:N]
-  psi0 = productMPS(sites,state)
+  psi0 = MPS(sites,state)
   @show flux(psi0)
 
   sweeps = Sweeps(5)
@@ -413,7 +413,7 @@ let
   end
 
   # Initialize psi to be a product state (alternating up and down)
-  psi = productMPS(s, n -> isodd(n) ? "Up" : "Dn")
+  psi = MPS(s, n -> isodd(n) ? "Up" : "Dn")
 
   c = div(N,2)
 
@@ -490,7 +490,7 @@ we will take to be near the center of the MPS. The details of this
 function are outside the scope of this tutorial, but are explained in the
 example code for measuring MPS.
 
-The line of code `psi = productMPS(s, n -> isodd(n) ? "Up" : "Dn")`
+The line of code `psi = MPS(s, n -> isodd(n) ? "Up" : "Dn")`
 initializes our MPS `psi` as a product state of alternating 
 up and down spins. We call `measure_Sz` before starting the
 time evolution.
diff --git a/src/mps/abstractmps.jl b/src/mps/abstractmps.jl
index cf15835c31..f4bacac421 100644
--- a/src/mps/abstractmps.jl
+++ b/src/mps/abstractmps.jl
@@ -1774,7 +1774,7 @@ function product(
   @show ortho_lims(ψ)
 
   if first(ortho_lims(ψ)) < first(ns_range) || last(ortho_lims(ψ)) > last(ns_range)
-    error("Orthogonality center limits are $(ortho_lims(ψ)), but gate is being applied to $ns′. Must be within the range of where the gates are being applied.")
+    error("Orthogonality center limits are $(ortho_lims(ψ)), but gate is being applied to sites $ns_range. Orthogonality center must be within the range of where the gates are being applied.")
   end
 
   ϕ = prod(ψ[ns_range])
diff --git a/test/autompo.jl b/test/autompo.jl
index ff65919aea..ff6901f806 100644
--- a/test/autompo.jl
+++ b/test/autompo.jl
@@ -723,15 +723,15 @@ end
     a3 += "Cdag", 1, "N", 2, "C", 3
     M3 = MPO(a3, s)
 
-    p011 = productMPS(s, [1, 2, 2, 1, 1])
-    p110 = productMPS(s, [2, 2, 1, 1, 1])
+    p011 = MPS(s, [1, 2, 2, 1, 1])
+    p110 = MPS(s, [2, 2, 1, 1, 1])
 
     @test inner(p110, M1, p011) ≈ -1.0
     @test inner(p110, M2, p011) ≈ -1.0
     @test inner(p110, M3, p011) ≈ -1.0
 
-    p001 = productMPS(s, [1, 1, 2, 1, 1])
-    p100 = productMPS(s, [2, 1, 1, 1, 1])
+    p001 = MPS(s, [1, 1, 2, 1, 1])
+    p100 = MPS(s, [2, 1, 1, 1, 1])
 
     @test inner(p100, M1, p001) ≈ +1.0
     @test inner(p100, M2, p001) ≈ +1.0
@@ -752,14 +752,14 @@ end
     a2 += -1, "Cdn", 3, "Cdagdn", 1
     M2 = MPO(a2, s)
 
-    p0uu = productMPS(s, [1, 2, 2, 1, 1])
-    puu0 = productMPS(s, [2, 2, 1, 1, 1])
-    p0ud = productMPS(s, [1, 2, 3, 1, 1])
-    pdu0 = productMPS(s, [3, 2, 1, 1, 1])
-    p00u = productMPS(s, [1, 1, 2, 1, 1])
-    pu00 = productMPS(s, [2, 1, 1, 1, 1])
-    p00d = productMPS(s, [1, 1, 3, 1, 1])
-    pd00 = productMPS(s, [3, 1, 1, 1, 1])
+    p0uu = MPS(s, [1, 2, 2, 1, 1])
+    puu0 = MPS(s, [2, 2, 1, 1, 1])
+    p0ud = MPS(s, [1, 2, 3, 1, 1])
+    pdu0 = MPS(s, [3, 2, 1, 1, 1])
+    p00u = MPS(s, [1, 1, 2, 1, 1])
+    pu00 = MPS(s, [2, 1, 1, 1, 1])
+    p00d = MPS(s, [1, 1, 3, 1, 1])
+    pd00 = MPS(s, [3, 1, 1, 1, 1])
 
     @test inner(puu0, M1, p0uu) ≈ -1.0
     @test inner(pdu0, M2, p0ud) ≈ -1.0
@@ -778,8 +778,8 @@ end
         ampo += -1im, "S-", i, "S+", i + 1
       end
       H = MPO(ampo, sites)
-      psiud = productMPS(sites, [1, 2, 1, 2])
-      psidu = productMPS(sites, [2, 1, 1, 2])
+      psiud = MPS(sites, [1, 2, 1, 2])
+      psidu = MPS(sites, [2, 1, 1, 2])
       @test inner(psiud, H, psidu) ≈ +1im
       @test inner(psidu, H, psiud) ≈ -1im
     end
diff --git a/test/dmrg.jl b/test/dmrg.jl
index 25c9185898..ca3f6971a1 100644
--- a/test/dmrg.jl
+++ b/test/dmrg.jl
@@ -40,7 +40,7 @@ using ITensors, Test, Random
     H = MPO(ampo, sites)
 
     state = [isodd(n) ? "Up" : "Dn" for n in 1:N]
-    psi = randomMPS(sites, state, 4)
+    psi = randomMPS(sites, state; linkdims=4)
 
     sweeps = Sweeps(3)
     @test length(sweeps) == 3
@@ -67,7 +67,7 @@ using ITensors, Test, Random
     H = MPO(ampo, sites)
 
     state = [isodd(n) ? "Up" : "Dn" for n in 1:N]
-    psi = randomMPS(sites, state, 4)
+    psi = randomMPS(sites, state; linkdims=4)
 
     sweeps = Sweeps(3)
     @test length(sweeps) == 3
@@ -94,7 +94,7 @@ using ITensors, Test, Random
     H = MPO(ampo, sites)
 
     state = [isodd(n) ? "Up" : "Dn" for n in 1:N]
-    psi = randomMPS(sites, state, 4)
+    psi = randomMPS(sites, state; linkdims=4)
     PH = ProjMPO(H)
 
     n = 4
@@ -254,7 +254,7 @@ using ITensors, Test, Random
     end
     H = MPO(ampo, sites)
 
-    psi0i = randomMPS(sites, 10)
+    psi0i = randomMPS(sites; linkdims=10)
 
     sweeps = Sweeps(4)
     maxdim!(sweeps, 10, 20, 100, 100)
@@ -264,7 +264,7 @@ using ITensors, Test, Random
     energy0, psi0 = dmrg(H, psi0i, sweeps; outputlevel=0)
     @test energy0 < -11.5
 
-    psi1i = randomMPS(sites, 10)
+    psi1i = randomMPS(sites; linkdims=10)
     energy1, psi1 = dmrg(H, [psi0], psi1i, sweeps; outputlevel=0, weight=weight)
 
     @test energy1 > energy0
@@ -285,7 +285,7 @@ using ITensors, Test, Random
     state[3] = 2
     state[5] = 2
     state[7] = 2
-    psi0 = productMPS(s, state)
+    psi0 = MPS(s, state)
 
     ampo = OpSum()
     for j in 1:(N - 1)
@@ -331,7 +331,7 @@ using ITensors, Test, Random
     maxdim!(sweeps, 50, 100, 200, 400, 800, 800)
     cutoff!(sweeps, 1E-10)
     state = ["Up", "Dn", "Dn", "Up", "Emp", "Up", "Up", "Emp", "Dn", "Dn"]
-    psi0 = randomMPS(sites, state, 10)
+    psi0 = randomMPS(sites, state; linkdims=10)
     energy, psi = dmrg(H, psi0, sweeps; outputlevel=0)
     @test (-8.02 < energy < -8.01)
   end
@@ -352,7 +352,7 @@ using ITensors, Test, Random
     maxdim!(sweeps, 10)
     cutoff!(sweeps, 1E-11)
 
-    psi0 = randomMPS(sites, 4)
+    psi0 = randomMPS(sites; linkdims=4)
 
     # Test that input works with wrong ortho center:
     orthogonalize!(psi0, 5)
@@ -395,7 +395,7 @@ using ITensors, Test, Random
     H = MPO(ampo, sites)
     Hsplit = splitblocks(linkinds, H)
     state = [isodd(n) ? "↑" : "↓" for n in 1:N]
-    ψ0 = productMPS(sites, state)
+    ψ0 = MPS(sites, state)
     using_threaded_blocksparse = ITensors.enable_threaded_blocksparse()
     energy, _ = dmrg(H, ψ0, sweeps; outputlevel=outputlevel)
     energy_split, _ = dmrg(Hsplit, ψ0, sweeps; outputlevel=outputlevel)
diff --git a/test/mps.jl b/test/mps.jl
index a1b0cc0bb2..b05f0292d1 100644
--- a/test/mps.jl
+++ b/test/mps.jl
@@ -60,7 +60,7 @@ include("util.jl")
     @test psi ⋅ psi ≈ *(dag(psi)..., psi...)[]
   end
 
-  @testset "productMPS" begin
+  @testset "MPS" begin
     @testset "vector of string input" begin
       sites = siteinds("S=1/2", N)
       state = fill("", N)
@@ -72,13 +72,13 @@ include("util.jl")
         sign = isodd(j) ? +1.0 : -1.0
         @test (psi[j] * op(sites, "Sz", j) * dag(prime(psi[j], "Site")))[] ≈ sign / 2
       end
-      psi = productMPS(sites, state)
+      psi = MPS(sites, state)
       for j in 1:N
         sign = isodd(j) ? +1.0 : -1.0
         @test (psi[j] * op(sites, "Sz", j) * dag(prime(psi[j], "Site")))[] ≈ sign / 2
       end
       @test_throws DimensionMismatch MPS(sites, fill("", N - 1))
-      @test_throws DimensionMismatch productMPS(sites, fill("", N - 1))
+      @test_throws DimensionMismatch MPS(sites, fill("", N - 1))
     end
 
     @testset "String input" begin
@@ -88,7 +88,7 @@ include("util.jl")
         sign = -1.0
         @test (psi[j] * op(sites, "Sz", j) * dag(prime(psi[j], "Site")))[] ≈ sign / 2
       end
-      psi = productMPS(sites, "Dn")
+      psi = MPS(sites, "Dn")
       for j in 1:N
         sign = -1.0
         @test (psi[j] * op(sites, "Sz", j) * dag(prime(psi[j], "Site")))[] ≈ sign / 2
@@ -102,7 +102,7 @@ include("util.jl")
         sign = -1.0
         @test (psi[j] * op(sites, "Sz", j) * dag(prime(psi[j], "Site")))[] ≈ sign / 2
       end
-      psi = productMPS(sites, 2)
+      psi = MPS(sites, 2)
       for j in 1:N
         sign = -1.0
         @test (psi[j] * op(sites, "Sz", j) * dag(prime(psi[j], "Site")))[] ≈ sign / 2
@@ -120,7 +120,7 @@ include("util.jl")
         sign = isodd(j) ? +1.0 : -1.0
         @test (psi[j] * op(sites, "Sz", j) * dag(prime(psi[j], "Site")))[] ≈ sign / 2
       end
-      psi = productMPS(sites, state)
+      psi = MPS(sites, state)
       for j in 1:N
         sign = isodd(j) ? +1.0 : -1.0
         @test (psi[j] * op(sites, "Sz", j) * dag(prime(psi[j], "Site")))[] ≈ sign / 2
@@ -139,7 +139,7 @@ include("util.jl")
         sign = isodd(j) ? +1.0 : -1.0
         @test (psi[j] * op(sites, "Sz", j) * dag(prime(psi[j], "Site")))[] ≈ sign / 2
       end
-      psi = productMPS(ivals)
+      psi = MPS(ivals)
       for j in 1:N
         sign = isodd(j) ? +1.0 : -1.0
         @test (psi[j] * op(sites, "Sz", j) * dag(prime(psi[j], "Site")))[] ≈ sign / 2
@@ -151,7 +151,7 @@ include("util.jl")
         for j in 1:N
           @test eltype(psi[j]) == ComplexF64
         end
-        psi = productMPS(ComplexF64, sites, fill(1, N))
+        psi = MPS(ComplexF64, sites, fill(1, N))
         for j in 1:N
           @test eltype(psi[j]) == ComplexF64
         end
@@ -165,13 +165,13 @@ include("util.jl")
       psi = MPS([site], [1])
       @test psi[1][1] ≈ 1.0
       @test psi[1][2] ≈ 0.0
-      psi = productMPS([site], [1])
+      psi = MPS([site], [1])
       @test psi[1][1] ≈ 1.0
       @test psi[1][2] ≈ 0.0
       psi = MPS([site], [2])
       @test psi[1][1] ≈ 0.0
       @test psi[1][2] ≈ 1.0
-      psi = productMPS([site], [2])
+      psi = MPS([site], [2])
       @test psi[1][1] ≈ 0.0
       @test psi[1][2] ≈ 1.0
     end
@@ -260,7 +260,7 @@ include("util.jl")
   end
 
   @testset "norm MPS" begin
-    psi = randomMPS(sites, 10)
+    psi = randomMPS(sites; linkdims=10)
     psidag = sim(linkinds, dag(psi))
     psi² = ITensor(1)
     for j in 1:N
@@ -269,11 +269,11 @@ include("util.jl")
     @test psi²[] ≈ psi ⋅ psi
     @test sqrt(psi²[]) ≈ norm(psi)
 
-    psi = randomMPS(sites, 10)
+    psi = randomMPS(sites; linkdims=10)
     psi .*= 1:N
     @test norm(psi) ≈ factorial(N)
 
-    psi = randomMPS(sites, 10)
+    psi = randomMPS(sites; linkdims=10)
     for j in 1:N
       psi[j] .*= j
     end
@@ -285,7 +285,7 @@ include("util.jl")
   end
 
   @testset "lognorm MPS" begin
-    psi = randomMPS(sites, 10)
+    psi = randomMPS(sites; linkdims=10)
     for j in 1:N
       psi[j] .*= j
     end
@@ -351,9 +351,9 @@ include("util.jl")
     s = siteinds("S=1/2", N; conserve_qns=conserve_qns)
     state = n -> isodd(n) ? "↑" : "↓"
 
-    ψ₁ = randomMPS(s, state, 4)
-    ψ₂ = randomMPS(s, state, 4)
-    ψ₃ = randomMPS(s, state, 4)
+    ψ₁ = randomMPS(s, state; linkdims=4)
+    ψ₂ = randomMPS(s, state; linkdims=4)
+    ψ₃ = randomMPS(s, state; linkdims=4)
 
     ψ = ψ₁ + ψ₂
 
@@ -535,7 +535,7 @@ end
   @testset "sample! method" begin
     N = 10
     sites = [Index(3, "Site,n=$n") for n in 1:N]
-    psi = randomMPS(sites, 3)
+    psi = randomMPS(sites; linkdims=3)
     nrm2 = inner(psi, psi)
     psi[1] *= (1.0 / sqrt(nrm2))
 
@@ -566,7 +566,7 @@ end
     N = 20
     chi = 8
     sites = siteinds(2, N)
-    M = randomMPS(sites, chi)
+    M = randomMPS(sites; linkdims=chi)
 
     @test ITensors.leftlim(M) == 0
     @test ITensors.rightlim(M) == 2
@@ -586,7 +586,7 @@ end
     end
 
     # Complex case
-    Mc = randomMPS(sites, chi)
+    Mc = randomMPS(sites; linkdims=chi)
     @test inner(Mc, Mc) ≈ 1.0 + 0.0im
   end
 
@@ -597,7 +597,7 @@ end
 
     # Make flux-zero random MPS
     state = [isodd(n) ? 1 : 2 for n in 1:N]
-    M = randomMPS(sites, state, chi)
+    M = randomMPS(sites, state; linkdims=chi)
     @test flux(M) == QN("Sz", 0)
 
     @test ITensors.leftlim(M) == 0
@@ -610,10 +610,10 @@ end
 
     # Test making random MPS with different flux
     state[1] = 2
-    M = randomMPS(sites, state, chi)
+    M = randomMPS(sites, state; linkdims=chi)
     @test flux(M) == QN("Sz", -2)
     state[3] = 2
-    M = randomMPS(sites, state, chi)
+    M = randomMPS(sites, state; linkdims=chi)
     @test flux(M) == QN("Sz", -4)
   end
 
@@ -623,7 +623,7 @@ end
 
     # Non-fermionic case - spin system
     s = siteinds("S=1/2", N; conserve_qns=true)
-    psi = randomMPS(s, n -> isodd(n) ? "Up" : "Dn", m)
+    psi = randomMPS(s, n -> isodd(n) ? "Up" : "Dn"; linkdims=m)
     Cpm = correlation_matrix(psi, "S+", "S-")
     # Check using OpSum:
     for i in 1:N, j in i:N
@@ -639,7 +639,7 @@ end
 
     # With start_site, end_site arguments:
     s = siteinds("S=1/2", N)
-    psi = randomMPS(ComplexF64, s, m)
+    psi = randomMPS(ComplexF64, s; linkdims=m)
     ss, es = 3, 6
     Nb = es - ss + 1
     Cpm = correlation_matrix(psi, "S+", "S-"; site_range=ss:es)
@@ -654,7 +654,7 @@ end
 
     # Fermionic case
     s = siteinds("Electron", N)
-    psi = randomMPS(s, m)
+    psi = randomMPS(s; linkdims=m)
     Cuu = correlation_matrix(psi, "Cdagup", "Cup")
     # Check using OpSum:
     for i in 1:N, j in i:N
@@ -874,7 +874,7 @@ end
     @test ITensors.orthocenter(ψ) == N
     @test maxlinkdim(ψ) == 1
 
-    ψ0 = randomMPS(s, 2)
+    ψ0 = randomMPS(s; linkdims=2)
     A = prod(ψ0)
     ψ = MPS(A, s; cutoff=1e-15, orthocenter=2)
     @test prod(ψ) ≈ A
@@ -914,7 +914,7 @@ end
   @testset "Set range of MPS tensors" begin
     N = 5
     s = siteinds("S=1/2", N)
-    ψ0 = randomMPS(s, 3)
+    ψ0 = randomMPS(s; linkdims=3)
 
     ψ = orthogonalize(ψ0, 2)
     A = prod(ITensors.data(ψ)[2:(N - 1)])
@@ -1545,7 +1545,7 @@ end
   @testset "dense conversion of MPS" begin
     N = 4
     s = siteinds("S=1/2", N; conserve_qns=true)
-    QM = randomMPS(s, ["Up", "Dn", "Up", "Dn"], 4)
+    QM = randomMPS(s, ["Up", "Dn", "Up", "Dn"]; linkdims=4)
     qsz1 = scalar(QM[1] * op("Sz", s[1]) * dag(prime(QM[1], "Site")))
 
     M = dense(QM)
@@ -1579,7 +1579,7 @@ end
       a .+= "Sz", j, "Sz", j + 1
     end
     H = MPO(a, s)
-    ψ = randomMPS(s, n -> isodd(n) ? "↑" : "↓", 10)
+    ψ = randomMPS(s, n -> isodd(n) ? "↑" : "↓"; linkdims=10)
     # Create MPO/MPS with pairs of sites merged
     H2 = MPO([H[b] * H[b + 1] for b in 1:2:N])
     ψ2 = MPS([ψ[b] * ψ[b + 1] for b in 1:2:N])

From 925fd40c78a5bd4103804d3d1cf3a9b88b70da0b Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Tue, 18 May 2021 16:36:14 -0400
Subject: [PATCH 03/12] Use location of next gate to choose a better
 orthogonality center

---
 src/mps/abstractmps.jl | 42 +++++++++++++++++++++++++++++++++---------
 1 file changed, 33 insertions(+), 9 deletions(-)

diff --git a/src/mps/abstractmps.jl b/src/mps/abstractmps.jl
index f4bacac421..d3767b738a 100644
--- a/src/mps/abstractmps.jl
+++ b/src/mps/abstractmps.jl
@@ -1437,8 +1437,8 @@ end
 # TODO: add a version that determines the sites
 # from common site indices of ψ and A
 """
-    setindex!(ψ::Union{MPS, MPO}, A::ITensor, r::UnitRange{Int};
-              orthocenter::Int = last(r), perm = nothing, kwargs...)
+    setindex!(ψ::Union{MPS,MPO}, A::ITensor, r::UnitRange{Int};
+              orthocenter::Int=last(r), perm=nothing, kwargs...)
     replacesites!([...])
     replacesites([...])
 
@@ -1739,6 +1739,10 @@ function product(
   next_gate::Union{ITensor,Nothing}=nothing,
   kwargs...,
 )
+
+  println("\n###################################")
+  println("product(::ITensor, ::$(typeof(ψ)))\n")
+
   N = length(ns)
   ns = sort(ns)
 
@@ -1747,12 +1751,6 @@ function product(
   @show ns
   @show next_ns
 
-  if !isnothing(next_ns)
-    if first(next_ns) > last(ns)
-      @show first(next_ns) > last(ns)
-    end
-  end
-
   diff_ns = diff(ns)
   ns′ = ns
   if any(!=(1), diff_ns)
@@ -1780,7 +1778,33 @@ function product(
   ϕ = prod(ψ[ns_range])
   ϕ = product(o, ϕ; apply_dag=apply_dag)
 
-  ψ[ns_range, kwargs...] = ϕ
+  println("\n Determine next desired orthocenter")
+
+  @show ns_range
+  @show next_ns
+
+  # Anticipate where the next orthogonality
+  # center should be based on the position
+  # of the next gate being applied
+  orthocenter = if !isnothing(next_ns)
+    if first(next_ns) > last(ns_range)
+      @show first(next_ns) > last(ns_range)
+      last(ns_range)
+    elseif last(next_ns) < first(ns_range)
+      first(ns_range)
+    else
+      last(next_ns)
+    end
+  else
+    last(ns_range)
+  end
+
+  @show orthocenter
+
+  ψ[ns_range, orthocenter=orthocenter, kwargs...] = ϕ
+
+  @show ortho_lims(ψ)
+
   if move_sites_back
     # Move the sites back to their original positions
     ψ = movesites(ψ, ns′ .=> ns; kwargs...)

From 53d51489ab39657dd8f42cdee90a57f8bbd8bf93 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Wed, 19 May 2021 17:43:42 -0400
Subject: [PATCH 04/12] Optimize ITensor -> MPS constructor by sweep outwards
 towards the desired orthogonality center

---
 src/mps/abstractmps.jl | 68 +++++++++++++-----------------------------
 1 file changed, 21 insertions(+), 47 deletions(-)

diff --git a/src/mps/abstractmps.jl b/src/mps/abstractmps.jl
index d3767b738a..b34b4fb4e6 100644
--- a/src/mps/abstractmps.jl
+++ b/src/mps/abstractmps.jl
@@ -1490,7 +1490,7 @@ function setindex!(
   # For MPO case, restrict to 0 prime level
   #sites = filter(hasplev(0), sites)
 
-  if !isnothing(perm)
+  if !isnothing(perm) # TODO: also check if it is a trivial permutation like [1,2,3]
     sites0 = sites
     sites = sites0[[perm...]]
     # Check if the site indices
@@ -1526,12 +1526,8 @@ function setindex!(
       end
     end
   end
-
   ψA = MPST(A, sites; leftinds=lind, orthocenter=orthocenter - first(r) + 1, kwargs...)
-  #@assert prod(ψA) ≈ A
-
   ψ[firstsite:lastsite] = ψA
-
   return ψ
 end
 
@@ -1570,30 +1566,34 @@ function (::Type{MPST})(
     @assert hasinds(A, s)
   end
   @assert isnothing(leftinds) || hasinds(A, leftinds)
-
   @assert 1 ≤ orthocenter ≤ N
 
+  l = isnothing(leftinds) ? Index[] : leftinds
+  # Need to flatten with reduce(vcat, sites)
+  r = uniqueinds(A, reduce(vcat, sites), l)
   ψ = Vector{ITensor}(undef, N)
   Ã = A
-  l = leftinds
-  # TODO: To minimize work, loop from
-  # 1:orthocenter and reverse(orthocenter:N)
-  # so the orthogonality center is set correctly.
-  for n in 1:(N - 1)
-    Lis = IndexSet(sites[n])
-    if !isnothing(l)
-      Lis = unioninds(Lis, l)
-    end
+
+  # Sweep left to orthocenter
+  for n in 1:(orthocenter - 1)
+    Lis = unioninds(sites[n], l)
     L, R = factorize(Ã, Lis; kwargs..., tags="Link,n=$n", ortho="left")
-    l = commonind(L, R)
+    l = commoninds(L, R)
     ψ[n] = L
     Ã = R
   end
-  ψ[N] = Ã
-  M = MPST(ψ)
-  setleftlim!(M, N - 1)
-  setrightlim!(M, N + 1)
-  orthogonalize!(M, orthocenter)
+
+  # Sweep right to orthocenter
+  for n in reverse((orthocenter + 1):N)
+    Ris = unioninds(sites[n], r)
+    R, L = factorize(Ã, Ris; kwargs..., tags="Link,n=$n", ortho="left")
+    r = commonind(R, L)
+    ψ[n] = R
+    Ã = L
+  end
+  ψ[orthocenter] = Ã
+  M = MPST(ψ; ortho_lims=orthocenter:orthocenter)
+  @assert ortho_lims(M) == orthocenter:orthocenter
   return M
 end
 
@@ -1739,18 +1739,10 @@ function product(
   next_gate::Union{ITensor,Nothing}=nothing,
   kwargs...,
 )
-
-  println("\n###################################")
-  println("product(::ITensor, ::$(typeof(ψ)))\n")
-
   N = length(ns)
   ns = sort(ns)
 
   next_ns = isnothing(next_gate) ? nothing : sort(findsites(ψ, next_gate))
-
-  @show ns
-  @show next_ns
-
   diff_ns = diff(ns)
   ns′ = ns
   if any(!=(1), diff_ns)
@@ -1758,9 +1750,6 @@ function product(
     ψ = movesites(ψ, ns .=> ns′; kwargs...)
   end
   ns_range = first(ns′):last(ns′)
-
-  @show ns_range
-  @show ortho_lims(ψ)
   # Minimize distance to orthogonalize
   if first(ortho_lims(ψ)) > last(ns_range)
     ψ = orthogonalize(ψ, last(ns_range))
@@ -1768,9 +1757,6 @@ function product(
     ψ = orthogonalize(ψ, first(ns_range))
   end
 
-  @show ns_range
-  @show ortho_lims(ψ)
-
   if first(ortho_lims(ψ)) < first(ns_range) || last(ortho_lims(ψ)) > last(ns_range)
     error("Orthogonality center limits are $(ortho_lims(ψ)), but gate is being applied to sites $ns_range. Orthogonality center must be within the range of where the gates are being applied.")
   end
@@ -1778,17 +1764,11 @@ function product(
   ϕ = prod(ψ[ns_range])
   ϕ = product(o, ϕ; apply_dag=apply_dag)
 
-  println("\n Determine next desired orthocenter")
-
-  @show ns_range
-  @show next_ns
-
   # Anticipate where the next orthogonality
   # center should be based on the position
   # of the next gate being applied
   orthocenter = if !isnothing(next_ns)
     if first(next_ns) > last(ns_range)
-      @show first(next_ns) > last(ns_range)
       last(ns_range)
     elseif last(next_ns) < first(ns_range)
       first(ns_range)
@@ -1798,13 +1778,7 @@ function product(
   else
     last(ns_range)
   end
-
-  @show orthocenter
-
   ψ[ns_range, orthocenter=orthocenter, kwargs...] = ϕ
-
-  @show ortho_lims(ψ)
-
   if move_sites_back
     # Move the sites back to their original positions
     ψ = movesites(ψ, ns′ .=> ns; kwargs...)

From 177afb16da8b945650e66bbdedc930996b564542 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Thu, 20 May 2021 16:25:28 -0400
Subject: [PATCH 05/12] Add splitblocks example to benchmark. Start creating
 logic for minimizing the number of swapgates

---
 benchmark/bench_dmrg.jl |  6 ++++
 src/mps/abstractmps.jl  | 64 +++++++++++++++++++++++++++++++++++++++--
 2 files changed, 67 insertions(+), 3 deletions(-)

diff --git a/benchmark/bench_dmrg.jl b/benchmark/bench_dmrg.jl
index 3732634b82..7b94bbc32a 100644
--- a/benchmark/bench_dmrg.jl
+++ b/benchmark/bench_dmrg.jl
@@ -45,6 +45,12 @@ let N = 100
   dmrg(H, psi, sweeps; outputlevel=0)
 
   suite["1d_S=1_heisenberg_qn"] = @benchmarkable dmrg($H, $psi, $sweeps; outputlevel=0)
+
+  # Precompile
+  H_splitblocks = splitblocks(linkinds, H)
+  dmrg(H_splitblocks, psi, sweeps; outputlevel=0)
+
+  suite["1d_S=1_heisenberg_qn_splitblocks"] = @benchmarkable dmrg($H_splitblocks, $psi, $sweeps; outputlevel=0)
 end
 
 #let Ny = 6,
diff --git a/src/mps/abstractmps.jl b/src/mps/abstractmps.jl
index b34b4fb4e6..04a6666bfe 100644
--- a/src/mps/abstractmps.jl
+++ b/src/mps/abstractmps.jl
@@ -1712,6 +1712,53 @@ function movesites(ψ::AbstractMPS, ns, ns′; kwargs...)
   return ψ
 end
 
+eachfront(itr, n = 1) = Iterators.take(itr, length(itr) - n)
+
+function areconsecutive(v)
+  for n in eachfront(eachindex(v))
+    if v[n] + 1 ≠ v[n + 1]
+      return false
+    end
+  end
+  return true
+end
+
+function minimal_swap_range(ns, ::Nothing)
+  new_start = first(ns)
+  return new_start:(new_start + length(ns) - 1)
+end
+
+# Compute the contiguous range of sites that
+# involve the minimal number of swaps to apply 
+# the current gate and the next gate
+function minimal_swap_range(ns, next_ns)
+  N = length(ns)
+  new_start = if first(next_ns) > last(ns)
+    # Next gate is completely to the right of the current gate
+    @show first(next_ns) > last(ns)
+    last(ns) - N + 1
+  elseif last(next_ns) < first(ns)
+    # Next gate is completely to the left of the current gate
+    @show last(next_ns) < first(ns)
+    first(ns)
+  elseif first(ns) < first(next_ns) && last(ns) < last(next_ns) && last(ns) ≥ first(next_ns)
+    @show first(ns) < first(next_ns) && last(ns) < last(next_ns) && last(ns) ≥ first(next_ns)
+    first(next_ns) - N + 1
+  elseif first(ns) ≤ first(next_ns) && last(ns) ≥ last(next_ns)
+    @show first(ns) ≤ first(next_ns) && last(ns) ≥ last(next_ns)
+    first(next_ns)
+  elseif first(ns) < first(next_ns) && last(ns) < last(next_ns)
+    @show first(ns) < first(next_ns) && last(ns) < last(next_ns)
+    first(ns)
+  elseif first(ns) ≥ first(next_ns) && last(ns) ≥ last(next_ns)
+    @show first(ns) ≥ first(next_ns) && last(ns) ≥ last(next_ns)
+    last(ns) - N + 1
+  else
+    error("The current gate is being applied to the sites $ns and the next gate is being applied to the sites $next_ns. Could not determine where to move the current sites to make the contiguous.")
+  end
+  return new_start:(new_start + N - 1)
+end
+
 """
     product(o::ITensor, ψ::Union{MPS,MPO}, [ns::Vector{Int}]; <keyword argument>)
     apply([...])
@@ -1743,10 +1790,21 @@ function product(
   ns = sort(ns)
 
   next_ns = isnothing(next_gate) ? nothing : sort(findsites(ψ, next_gate))
-  diff_ns = diff(ns)
   ns′ = ns
-  if any(!=(1), diff_ns)
-    ns′ = [ns[1] + n - 1 for n in 1:N]
+
+  println("\n##############################")
+  println("Inside product(::ITensor, ::$(typeof(ψ)))")
+  @show ns
+  @show next_ns
+  @show areconsecutive(ns)
+
+  if !areconsecutive(ns) #any(!=(1), diff_ns)
+    # TODO: change the new position depending on next_ns
+    ns′ = minimal_swap_range(ns, next_ns)
+
+    #ns′ = [ns[1] + n - 1 for n in 1:N]
+    @show ns′
+
     ψ = movesites(ψ, ns .=> ns′; kwargs...)
   end
   ns_range = first(ns′):last(ns′)

From 863cb720598b3432c4a17ddd8ce14d859717a732 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Fri, 21 May 2021 11:07:44 -0400
Subject: [PATCH 06/12] Add logic for minimizing the number of swaps

---
 src/mps/abstractmps.jl | 56 ++++++++++++++++++++++--------------------
 test/apply.jl          | 31 +++++++++++++++++++++++
 test/runtests.jl       |  1 +
 3 files changed, 61 insertions(+), 27 deletions(-)
 create mode 100644 test/apply.jl

diff --git a/src/mps/abstractmps.jl b/src/mps/abstractmps.jl
index 04a6666bfe..f8fd46597a 100644
--- a/src/mps/abstractmps.jl
+++ b/src/mps/abstractmps.jl
@@ -1731,30 +1731,32 @@ end
 # Compute the contiguous range of sites that
 # involve the minimal number of swaps to apply 
 # the current gate and the next gate
-function minimal_swap_range(ns, next_ns)
+function minimal_swap_range(ns, ns_next)
   N = length(ns)
-  new_start = if first(next_ns) > last(ns)
-    # Next gate is completely to the right of the current gate
-    @show first(next_ns) > last(ns)
+  N_next = length(ns_next)
+  new_start = if last(ns) ≤ first(ns_next)
+    # Next gate is to the right of the current gate
+    # Move right towards the current gate
     last(ns) - N + 1
-  elseif last(next_ns) < first(ns)
-    # Next gate is completely to the left of the current gate
-    @show last(next_ns) < first(ns)
+  elseif first(ns) ≥ last(ns_next)
+    # Next gate is to the left of the current gate
+    # Move left towards the current gate
     first(ns)
-  elseif first(ns) < first(next_ns) && last(ns) < last(next_ns) && last(ns) ≥ first(next_ns)
-    @show first(ns) < first(next_ns) && last(ns) < last(next_ns) && last(ns) ≥ first(next_ns)
-    first(next_ns) - N + 1
-  elseif first(ns) ≤ first(next_ns) && last(ns) ≥ last(next_ns)
-    @show first(ns) ≤ first(next_ns) && last(ns) ≥ last(next_ns)
-    first(next_ns)
-  elseif first(ns) < first(next_ns) && last(ns) < last(next_ns)
-    @show first(ns) < first(next_ns) && last(ns) < last(next_ns)
+  elseif first(ns) ≤ first(ns_next) && last(ns) ≥ last(ns_next)
+    # Next gate is completely within the current gate,
+    # move inwards to the start of the current gate
+    first(ns_next)
+  elseif first(ns) ≥ first(ns_next) && last(ns) ≤ last(ns_next)
+    # Current gate is completely within the next gate
     first(ns)
-  elseif first(ns) ≥ first(next_ns) && last(ns) ≥ last(next_ns)
-    @show first(ns) ≥ first(next_ns) && last(ns) ≥ last(next_ns)
-    last(ns) - N + 1
+  elseif first(ns) ≤ first(ns_next) && last(ns) ≤ last(ns_next)
+    # Next gate is overlapping to the right of the current gate
+    last(ns_next) - N_next - 1
+  elseif first(ns) ≥ first(ns_next) && last(ns) ≥ last(ns_next)
+    # Next gate is overlapping to the left of the current gate
+    first(ns_next) + N_next
   else
-    error("The current gate is being applied to the sites $ns and the next gate is being applied to the sites $next_ns. Could not determine where to move the current sites to make the contiguous.")
+    error("The current gate is being applied to the sites $ns and the next gate is being applied to the sites $ns_next. Could not determine where to move the current sites to make the contiguous.")
   end
   return new_start:(new_start + N - 1)
 end
@@ -1789,18 +1791,18 @@ function product(
   N = length(ns)
   ns = sort(ns)
 
-  next_ns = isnothing(next_gate) ? nothing : sort(findsites(ψ, next_gate))
+  ns_next = isnothing(next_gate) ? nothing : sort(findsites(ψ, next_gate))
   ns′ = ns
 
   println("\n##############################")
   println("Inside product(::ITensor, ::$(typeof(ψ)))")
   @show ns
-  @show next_ns
+  @show ns_next
   @show areconsecutive(ns)
 
   if !areconsecutive(ns) #any(!=(1), diff_ns)
-    # TODO: change the new position depending on next_ns
-    ns′ = minimal_swap_range(ns, next_ns)
+    # TODO: change the new position depending on ns_next
+    ns′ = minimal_swap_range(ns, ns_next)
 
     #ns′ = [ns[1] + n - 1 for n in 1:N]
     @show ns′
@@ -1825,13 +1827,13 @@ function product(
   # Anticipate where the next orthogonality
   # center should be based on the position
   # of the next gate being applied
-  orthocenter = if !isnothing(next_ns)
-    if first(next_ns) > last(ns_range)
+  orthocenter = if !isnothing(ns_next)
+    if first(ns_next) > last(ns_range)
       last(ns_range)
-    elseif last(next_ns) < first(ns_range)
+    elseif last(ns_next) < first(ns_range)
       first(ns_range)
     else
-      last(next_ns)
+      last(ns_next)
     end
   else
     last(ns_range)
diff --git a/test/apply.jl b/test/apply.jl
new file mode 100644
index 0000000000..7883edcc35
--- /dev/null
+++ b/test/apply.jl
@@ -0,0 +1,31 @@
+using ITensors
+using Test
+
+@testset "ITensors.minimal_swap_range" begin
+  @test ITensors.minimal_swap_range([1, 2], [1, 2]) == 1:2
+  @test ITensors.minimal_swap_range([1, 2], [3, 4]) == 1:2
+  @test ITensors.minimal_swap_range([1, 3], [5, 6]) == 2:3
+  @test ITensors.minimal_swap_range([1, 4], [5, 6]) == 3:4
+  @test ITensors.minimal_swap_range([1, 5], [5, 6]) == 4:5
+  @test ITensors.minimal_swap_range([1, 6], [5, 6]) == 5:6
+  @test ITensors.minimal_swap_range([1, 6], [4, 5, 6]) == 4:5
+  @test ITensors.minimal_swap_range([1, 3, 5], [7, 8]) == 3:5
+  @test ITensors.minimal_swap_range([1, 5], [2, 3]) == 2:3
+  @test ITensors.minimal_swap_range([1, 5], [2, 4]) == 2:3
+  @test ITensors.minimal_swap_range([1, 3, 5], [2, 4]) == 2:4
+  @test ITensors.minimal_swap_range([1, 3, 5], [5, 8]) == 3:5
+  @test ITensors.minimal_swap_range([3, 4], [1, 2]) == 3:4
+  @test ITensors.minimal_swap_range([4, 6], [1, 2]) == 4:5
+  @test ITensors.minimal_swap_range([4, 6], [1, 3]) == 4:5
+  @test ITensors.minimal_swap_range([4, 6], [1, 4]) == 4:5
+  @test ITensors.minimal_swap_range([2, 6], [1, 6]) == 2:3
+  @test ITensors.minimal_swap_range([2, 6], [1, 7]) == 2:3
+  @test ITensors.minimal_swap_range([2, 6], [1, 5]) == 3:4
+  @test ITensors.minimal_swap_range([2, 5, 6], [1, 5]) == 3:5
+  @test ITensors.minimal_swap_range([1, 5], [2, 6]) == 3:4
+  @test ITensors.minimal_swap_range([4, 8], [3, 8]) == 4:5
+  @test ITensors.minimal_swap_range([4, 6, 8], [3, 5]) == 5:7
+  @test ITensors.minimal_swap_range([4, 6, 8], [3, 4]) == 4:6
+  @test ITensors.minimal_swap_range([4, 6, 8], [3, 8]) == 4:6
+end
+
diff --git a/test/runtests.jl b/test/runtests.jl
index acc03f4c1f..f2324ff862 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -32,6 +32,7 @@ using Test
     "lattices.jl",
     "mps.jl",
     "mpo.jl",
+    "apply.jl",
     "sweeps.jl",
     "sweepnext.jl",
     "autompo.jl",

From 1d98be7b7522efa997e8b965972fcf3cfe7fed9d Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Fri, 21 May 2021 15:52:18 -0400
Subject: [PATCH 07/12] Fix bug in choosing new orthogonality center. Add tests
 for picking contiguous gate location based on next get location. Add counting
 number of swaps

---
 src/mps/abstractmps.jl |  76 ++++++++++----------
 test/apply.jl          | 158 ++++++++++++++++++++++++++++++++++-------
 2 files changed, 170 insertions(+), 64 deletions(-)

diff --git a/src/mps/abstractmps.jl b/src/mps/abstractmps.jl
index f8fd46597a..da9f3734fd 100644
--- a/src/mps/abstractmps.jl
+++ b/src/mps/abstractmps.jl
@@ -1666,13 +1666,19 @@ function _movesite(ns::Vector{Int}, n1n2::Pair{Int,Int})
   return ns
 end
 
-function _movesites(ψ::AbstractMPS, ns::Vector{Int}, ns′::Vector{Int}; kwargs...)
+function _movesites(ψ::AbstractMPS, ns::Vector{Int}, ns′::Vector{Int}; (nswaps!)=nothing, kwargs...)
   ψ = copy(ψ)
   N = length(ns)
   @assert N == length(ns′)
+  nswaps = 0
   for i in 1:N
-    ψ = movesite(ψ, ns[i] => ns′[i]; kwargs...)
-    ns = _movesite(ns, ns[i] => ns′[i])
+    n1, n2 = ns[i], ns′[i]
+    nswaps += n1 > n2 ? length(n2:(n1 - 1)) : length(n1:(n2 - 1))
+    ψ = movesite(ψ, n1 => n2; kwargs...)
+    ns = _movesite(ns, n1 => n2)
+  end
+  if !isnothing(nswaps!)
+    push!(nswaps!, nswaps)
   end
   return ψ, ns
 end
@@ -1680,7 +1686,7 @@ end
 # TODO: make a permutesites(::MPS/MPO, perm)
 # function that takes a permutation of the sites
 # p(1:N) for N sites
-function movesites(ψ::AbstractMPS, nsns′::Vector{Pair{Int,Int}}; kwargs...)
+function movesites(ψ::AbstractMPS, nsns′::Vector{Pair{Int,Int}}; (nswaps!)=nothing, kwargs...)
   ns = first.(nsns′)
   ns′ = last.(nsns′)
   ψ = copy(ψ)
@@ -1691,26 +1697,26 @@ function movesites(ψ::AbstractMPS, nsns′::Vector{Pair{Int,Int}}; kwargs...)
   ns′ = ns′[p]
   ns = collect(ns)
   while ns ≠ ns′
-    ψ, ns = _movesites(ψ, ns, ns′; kwargs...)
+    ψ, ns = _movesites(ψ, ns, ns′; (nswaps!)=nswaps!, kwargs...)
   end
   return ψ
 end
 
 # TODO: call the Vector{Pair{Int, Int}} version
-function movesites(ψ::AbstractMPS, ns, ns′; kwargs...)
-  ψ = copy(ψ)
-  N = length(ns)
-  @assert N == length(ns′)
-  p = sortperm(ns′)
-  ns = ns[p]
-  ns′ = ns′[p]
-  ns = collect(ns)
-  for i in 1:N
-    ψ = movesite(ψ, ns[i] => ns′[i]; kwargs...)
-    ns = _movesite(ns, ns[i] => ns′[i])
-  end
-  return ψ
-end
+## function movesites(ψ::AbstractMPS, ns, ns′; kwargs...)
+##   ψ = copy(ψ)
+##   N = length(ns)
+##   @assert N == length(ns′)
+##   p = sortperm(ns′)
+##   ns = ns[p]
+##   ns′ = ns′[p]
+##   ns = collect(ns)
+##   for i in 1:N
+##     ψ = movesite(ψ, ns[i] => ns′[i]; kwargs...)
+##     ns = _movesite(ns, ns[i] => ns′[i])
+##   end
+##   return ψ
+## end
 
 eachfront(itr, n = 1) = Iterators.take(itr, length(itr) - n)
 
@@ -1751,10 +1757,10 @@ function minimal_swap_range(ns, ns_next)
     first(ns)
   elseif first(ns) ≤ first(ns_next) && last(ns) ≤ last(ns_next)
     # Next gate is overlapping to the right of the current gate
-    last(ns_next) - N_next - 1
+    first(ns_next) - 1
   elseif first(ns) ≥ first(ns_next) && last(ns) ≥ last(ns_next)
     # Next gate is overlapping to the left of the current gate
-    first(ns_next) + N_next
+    last(ns_next) + 2 - N
   else
     error("The current gate is being applied to the sites $ns and the next gate is being applied to the sites $ns_next. Could not determine where to move the current sites to make the contiguous.")
   end
@@ -1786,6 +1792,7 @@ function product(
   move_sites_back::Bool=true,
   apply_dag::Bool=false,
   next_gate::Union{ITensor,Nothing}=nothing,
+  (nswaps!)=nothing,
   kwargs...,
 )
   N = length(ns)
@@ -1794,20 +1801,13 @@ function product(
   ns_next = isnothing(next_gate) ? nothing : sort(findsites(ψ, next_gate))
   ns′ = ns
 
-  println("\n##############################")
-  println("Inside product(::ITensor, ::$(typeof(ψ)))")
-  @show ns
-  @show ns_next
-  @show areconsecutive(ns)
-
   if !areconsecutive(ns) #any(!=(1), diff_ns)
-    # TODO: change the new position depending on ns_next
     ns′ = minimal_swap_range(ns, ns_next)
-
-    #ns′ = [ns[1] + n - 1 for n in 1:N]
-    @show ns′
-
-    ψ = movesites(ψ, ns .=> ns′; kwargs...)
+    ψ = movesites(ψ, ns .=> ns′; (nswaps!)=nswaps!, kwargs...)
+  else
+    if !isnothing(nswaps!)
+      push!(nswaps!, 0)
+    end
   end
   ns_range = first(ns′):last(ns′)
   # Minimize distance to orthogonalize
@@ -1815,10 +1815,10 @@ function product(
     ψ = orthogonalize(ψ, last(ns_range))
   elseif last(ortho_lims(ψ)) < first(ns_range)
     ψ = orthogonalize(ψ, first(ns_range))
-  end
-
-  if first(ortho_lims(ψ)) < first(ns_range) || last(ortho_lims(ψ)) > last(ns_range)
-    error("Orthogonality center limits are $(ortho_lims(ψ)), but gate is being applied to sites $ns_range. Orthogonality center must be within the range of where the gates are being applied.")
+  elseif first(ortho_lims(ψ)) < first(ns_range)
+    ψ = orthogonalize(ψ, first(ns_range))
+  elseif last(ortho_lims(ψ)) > last(ns_range)
+    ψ = orthogonalize(ψ, last(ns_range))
   end
 
   ϕ = prod(ψ[ns_range])
@@ -1833,7 +1833,7 @@ function product(
     elseif last(ns_next) < first(ns_range)
       first(ns_range)
     else
-      last(ns_next)
+      last(ns_range)
     end
   else
     last(ns_range)
diff --git a/test/apply.jl b/test/apply.jl
index 7883edcc35..f95ffd69f8 100644
--- a/test/apply.jl
+++ b/test/apply.jl
@@ -1,31 +1,137 @@
 using ITensors
 using Test
 
-@testset "ITensors.minimal_swap_range" begin
-  @test ITensors.minimal_swap_range([1, 2], [1, 2]) == 1:2
-  @test ITensors.minimal_swap_range([1, 2], [3, 4]) == 1:2
-  @test ITensors.minimal_swap_range([1, 3], [5, 6]) == 2:3
-  @test ITensors.minimal_swap_range([1, 4], [5, 6]) == 3:4
-  @test ITensors.minimal_swap_range([1, 5], [5, 6]) == 4:5
-  @test ITensors.minimal_swap_range([1, 6], [5, 6]) == 5:6
-  @test ITensors.minimal_swap_range([1, 6], [4, 5, 6]) == 4:5
-  @test ITensors.minimal_swap_range([1, 3, 5], [7, 8]) == 3:5
-  @test ITensors.minimal_swap_range([1, 5], [2, 3]) == 2:3
-  @test ITensors.minimal_swap_range([1, 5], [2, 4]) == 2:3
-  @test ITensors.minimal_swap_range([1, 3, 5], [2, 4]) == 2:4
-  @test ITensors.minimal_swap_range([1, 3, 5], [5, 8]) == 3:5
-  @test ITensors.minimal_swap_range([3, 4], [1, 2]) == 3:4
-  @test ITensors.minimal_swap_range([4, 6], [1, 2]) == 4:5
-  @test ITensors.minimal_swap_range([4, 6], [1, 3]) == 4:5
-  @test ITensors.minimal_swap_range([4, 6], [1, 4]) == 4:5
-  @test ITensors.minimal_swap_range([2, 6], [1, 6]) == 2:3
-  @test ITensors.minimal_swap_range([2, 6], [1, 7]) == 2:3
-  @test ITensors.minimal_swap_range([2, 6], [1, 5]) == 3:4
-  @test ITensors.minimal_swap_range([2, 5, 6], [1, 5]) == 3:5
-  @test ITensors.minimal_swap_range([1, 5], [2, 6]) == 3:4
-  @test ITensors.minimal_swap_range([4, 8], [3, 8]) == 4:5
-  @test ITensors.minimal_swap_range([4, 6, 8], [3, 5]) == 5:7
-  @test ITensors.minimal_swap_range([4, 6, 8], [3, 4]) == 4:6
-  @test ITensors.minimal_swap_range([4, 6, 8], [3, 8]) == 4:6
+@testset "apply" begin
+
+  @testset "ITensors.minimal_swap_range" begin
+    @test ITensors.minimal_swap_range([1, 2], [1, 2]) == 1:2
+    @test ITensors.minimal_swap_range([1, 2], [3, 4]) == 1:2
+    @test ITensors.minimal_swap_range([1, 3], [5, 6]) == 2:3
+    @test ITensors.minimal_swap_range([1, 4], [5, 6]) == 3:4
+    @test ITensors.minimal_swap_range([3, 6], [1, 4]) == 4:5
+    @test ITensors.minimal_swap_range([1, 5], [5, 6]) == 4:5
+    @test ITensors.minimal_swap_range([1, 6], [5, 6]) == 5:6
+    @test ITensors.minimal_swap_range([1, 6], [4, 5, 6]) == 4:5
+    @test ITensors.minimal_swap_range([1, 3, 5], [7, 8]) == 3:5
+    @test ITensors.minimal_swap_range([1, 5], [2, 3]) == 2:3
+    @test ITensors.minimal_swap_range([1, 5], [2, 4]) == 2:3
+    @test ITensors.minimal_swap_range([1, 3, 5], [2, 4]) == 2:4
+    @test ITensors.minimal_swap_range([1, 3, 5], [5, 8]) == 3:5
+    @test ITensors.minimal_swap_range([3, 4], [1, 2]) == 3:4
+    @test ITensors.minimal_swap_range([4, 6], [1, 2]) == 4:5
+    @test ITensors.minimal_swap_range([4, 6], [1, 3]) == 4:5
+    @test ITensors.minimal_swap_range([4, 6], [1, 4]) == 4:5
+    @test ITensors.minimal_swap_range([2, 6], [1, 6]) == 2:3
+    @test ITensors.minimal_swap_range([2, 6], [1, 7]) == 2:3
+    @test ITensors.minimal_swap_range([2, 6], [1, 5]) == 5:6
+    @test ITensors.minimal_swap_range([1, 5], [2, 6]) == 1:2
+    @test ITensors.minimal_swap_range([1, 3, 5], [2, 6]) == 1:3
+    @test ITensors.minimal_swap_range([2, 5, 6], [1, 5]) == 4:6
+    @test ITensors.minimal_swap_range([1, 5], [2, 6]) == 1:2
+    @test ITensors.minimal_swap_range([4, 8], [3, 8]) == 4:5
+    @test ITensors.minimal_swap_range([4, 6, 8], [3, 5]) == 4:6
+    @test ITensors.minimal_swap_range([4, 6, 8], [3, 4]) == 4:6
+    @test ITensors.minimal_swap_range([4, 6, 8], [3, 8]) == 4:6
+  end
+
+  @testset "apply swap" begin
+    N = 6
+    s = siteinds("Qubit", N)
+
+    Xop = [("X", n) for n in 1:N]
+    CXop = [("CX", i, j) for i in 1:N, j in 1:N]
+    CCXop = [("CCX", i, j, k) for i in 1:N, j in 1:N, k in 1:N]
+
+    X = op.(Xop, (s,))
+    CX = op.(CXop, (s,))
+    CCX = op.(CCXop, (s,))
+
+    ψ0 = MPS(s, "0")
+    orthogonalize!(ψ0, 1)
+    @test ortho_lims(ψ0) == 1:1
+
+    move_sites_back = false
+    kwargs = (; move_sites_back=move_sites_back)
+
+    nswaps = Int[]
+    ψ = apply(X[1], ψ0; (nswaps!)=nswaps, kwargs...)
+    @test ortho_lims(ψ) == 1:1
+    @test nswaps == [0]
+
+    nswaps = Int[]
+    ψ = apply(X[3], ψ0; (nswaps!)=nswaps, kwargs...)
+    @test ortho_lims(ψ) == 3:3
+    @test nswaps == [0]
+
+    nswaps = Int[]
+    ψ = apply(CX[1, 2], ψ0; (nswaps!)=nswaps, kwargs...)
+    @test ortho_lims(ψ) == 2:2
+    @test nswaps == [0]
+
+    nswaps = Int[]
+    ψ = apply(CX[2, 3], ψ0; (nswaps!)=nswaps, kwargs...)
+    @test ortho_lims(ψ) == 3:3
+    @test nswaps == [0]
+
+    nswaps = Int[]
+    ψ = apply(CX[3, 4], ψ0; (nswaps!)=nswaps, kwargs...)
+    @test ortho_lims(ψ) == 4:4
+    @test nswaps == [0]
+
+    nswaps = Int[]
+    ψ = apply([CX[3, 4], X[4]], ψ0; (nswaps!)=nswaps, kwargs...)
+    @test ortho_lims(ψ) == 4:4
+    @test nswaps == [0, 0]
+
+    nswaps = Int[]
+    ψ = apply(CX[1, 3], ψ0; (nswaps!)=nswaps, kwargs...)
+    @test ortho_lims(ψ) == 2:2
+    @test nswaps == [1]
+
+    nswaps = Int[]
+    ψ = apply(CX[1, 4], ψ0; (nswaps!)=nswaps, kwargs...)
+    @test ortho_lims(ψ) == 2:2
+    @test nswaps == [2]
+
+    nswaps = Int[]
+    ψ = apply([CCX[2, 3, 4], X[2]], ψ0; (nswaps!)=nswaps, kwargs...)
+    @test ortho_lims(ψ) == 2:2
+    @test nswaps == [0, 0]
+
+    nswaps = Int[]
+    ψ = apply([CX[3, 5], CX[1, 3]], ψ0; (nswaps!)=nswaps, kwargs...)
+
+    nswaps = Int[]
+    ψ = apply(CX[1, 4], ψ0; next_gate=CX[2, 3], (nswaps!)=nswaps, kwargs...)
+    @test nswaps == [2]
+
+    nswaps = Int[]
+    ψ = apply([CX[3, 6], CX[1, 3]], ψ0; (nswaps!)=nswaps, kwargs...)
+    @test nswaps == [2, 1]
+
+    @test ITensors.minimal_swap_range([3, 6], [1, 4]) == 4:5
+    nswaps = Int[]
+    ψ = apply([CX[3, 6], CX[1, 4]], ψ0; (nswaps!)=nswaps, kwargs...)
+    @test nswaps == [2, 1]
+
+    # XXX broken
+    nswaps = Int[]
+    ψ = apply([CX[1, 3], CX[3, 5]], ψ0; next_gate=CX[3, 5], (nswaps!)=nswaps, kwargs...)
+    @test nswaps == [1, 1]
+
+    # XXX broken
+    nswaps = Int[]
+    ψ = apply([CX[1, 3], CX[3, 5]], ψ0; next_gate=CX[3, 5], (nswaps!)=nswaps, kwargs...)
+    @test nswaps == [1, 1]
+
+  end
+
+  @testset "Debugging" begin
+    N = 4
+    s = siteinds("Qubit", N)
+    X = [op("X", s, n) for n in 1:N]
+    ψ = MPS(s, "0")
+    @test prod(product(X[1], ψ)) ≈ prod(MPS(s, n -> n == 1 ? "1" : "0"))
+  end
 end
 

From 6f51b75c9aa0cdcc1279efdf331acbee25a97bda Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Fri, 21 May 2021 16:17:33 -0400
Subject: [PATCH 08/12] Add keyword argument for choosing swapping strategy

---
 src/mps/abstractmps.jl | 11 +++++++++--
 1 file changed, 9 insertions(+), 2 deletions(-)

diff --git a/src/mps/abstractmps.jl b/src/mps/abstractmps.jl
index da9f3734fd..9ced4eac5e 100644
--- a/src/mps/abstractmps.jl
+++ b/src/mps/abstractmps.jl
@@ -1792,6 +1792,7 @@ function product(
   move_sites_back::Bool=true,
   apply_dag::Bool=false,
   next_gate::Union{ITensor,Nothing}=nothing,
+  swap_alg="minimize",
   (nswaps!)=nothing,
   kwargs...,
 )
@@ -1801,8 +1802,14 @@ function product(
   ns_next = isnothing(next_gate) ? nothing : sort(findsites(ψ, next_gate))
   ns′ = ns
 
-  if !areconsecutive(ns) #any(!=(1), diff_ns)
-    ns′ = minimal_swap_range(ns, ns_next)
+  if !areconsecutive(ns)
+    ns′ = if swap_alg == "minimize"
+      minimal_swap_range(ns, ns_next)
+    elseif swap_alg == "left"
+      ns′ = first(ns):(first(ns) + N - 1)
+    elseif swap_alg == "right"
+      ns′ = (last(ns) - N + 1):last(ns)
+    end
     ψ = movesites(ψ, ns .=> ns′; (nswaps!)=nswaps!, kwargs...)
   else
     if !isnothing(nswaps!)

From 504e618f491cb6e0e1c3f3798a180c331095d418 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Fri, 21 May 2021 17:58:32 -0400
Subject: [PATCH 09/12] Add customization of algorithm for deciding new gate
 position with consecutive sites

---
 src/mps/abstractmps.jl | 65 ++++++++++++++++++++++++++++++++++++------
 1 file changed, 56 insertions(+), 9 deletions(-)

diff --git a/src/mps/abstractmps.jl b/src/mps/abstractmps.jl
index 9ced4eac5e..e89eb37464 100644
--- a/src/mps/abstractmps.jl
+++ b/src/mps/abstractmps.jl
@@ -1729,15 +1729,56 @@ function areconsecutive(v)
   return true
 end
 
-function minimal_swap_range(ns, ::Nothing)
+abstract type SwapAlgorithm end
+struct SwapMinimal <: SwapAlgorithm end
+struct SwapMinimalSimple <: SwapAlgorithm end
+struct SwapLeft <: SwapAlgorithm end
+struct SwapRight <: SwapAlgorithm end
+struct SwapMiddle <: SwapAlgorithm end
+
+function consecutive_range(::SwapLeft, ns, ns_next)
+  N = length(ns)
   new_start = first(ns)
-  return new_start:(new_start + length(ns) - 1)
+  return new_start:(new_start + N - 1)
+end
+
+function consecutive_range(::SwapRight, ns, ns_next)
+  N = length(ns)
+  new_stop = last(ns)
+  return (new_stop - N + 1):new_stop
+end
+
+function consecutive_range(::SwapMiddle, ns, ns_next)
+  N = length(ns)
+  new_start = (first(ns) + last(ns)) ÷ 2 - (length(ns) - 1) ÷ 2
+  return new_start:(new_start + N - 1)
 end
 
-# Compute the contiguous range of sites that
-# involve the minimal number of swaps to apply 
-# the current gate and the next gate
-function minimal_swap_range(ns, ns_next)
+function consecutive_range(::SwapMinimalSimple, ns, ns_next)
+  N = length(ns)
+  N_next = length(ns_next)
+  new_start = if last(ns) ≤ first(ns_next)
+    # Next gate is to the right of the current gate
+    # Move right towards the current gate
+    last(ns) - N + 1
+  elseif first(ns) ≥ last(ns_next)
+    # Next gate is to the left of the current gate
+    # Move left towards the current gate
+    first(ns)
+  else
+    # If neither is true, just swap left
+    first(ns)
+  end
+  return new_start:(new_start + N - 1)
+end
+
+consecutive_range(::SwapMinimalSimple, ns, ns_next::Nothing) = consecutive_range(SwapRight(), ns, ns_next)
+
+# Compute a contiguous range of sites 
+# based on the current range of sites the gate
+# is being applied to as well as the range of
+# sites the next gate will be applied to
+function consecutive_range(::SwapMinimal, ns, ns_next)
   N = length(ns)
   N_next = length(ns_next)
   new_start = if last(ns) ≤ first(ns_next)
@@ -1767,6 +1808,8 @@ function minimal_swap_range(ns, ns_next)
   return new_start:(new_start + N - 1)
 end
 
+consecutive_range(::SwapMinimal, ns, ns_next::Nothing) = consecutive_range(SwapRight(), ns, ns_next)
+
 """
     product(o::ITensor, ψ::Union{MPS,MPO}, [ns::Vector{Int}]; <keyword argument>)
     apply([...])
@@ -1804,11 +1847,15 @@ function product(
 
   if !areconsecutive(ns)
     ns′ = if swap_alg == "minimize"
-      minimal_swap_range(ns, ns_next)
+      consecutive_range(SwapMinimal(), ns, ns_next)
+    elseif swap_alg == "simple_minimize"
+      consecutive_range(SwapMinimalSimple(), ns, ns_next)
     elseif swap_alg == "left"
-      ns′ = first(ns):(first(ns) + N - 1)
+      consecutive_range(SwapLeft(), ns, ns_next)
     elseif swap_alg == "right"
-      ns′ = (last(ns) - N + 1):last(ns)
+      consecutive_range(SwapRight(), ns, ns_next)
+    elseif swap_alg == "middle"
+      consecutive_range(SwapMiddle(), ns, ns_next)
     end
     ψ = movesites(ψ, ns .=> ns′; (nswaps!)=nswaps!, kwargs...)
   else

From 0b6b955c8a63c06e0436077a7624f6b9c881f9f5 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Tue, 25 May 2021 12:02:25 -0400
Subject: [PATCH 10/12] Extend MPO function to allow dangling indices. Allow
 for integer tags. Add dmrg with splitblocks benchmark.

---
 src/ITensors.jl |  4 ++++
 src/mps/mpo.jl  | 37 +++++++++++++++++++++++++++----------
 src/tagset.jl   |  6 +-----
 3 files changed, 32 insertions(+), 15 deletions(-)

diff --git a/src/ITensors.jl b/src/ITensors.jl
index 2bef5b5c52..199c239444 100644
--- a/src/ITensors.jl
+++ b/src/ITensors.jl
@@ -7,6 +7,10 @@ The ITensor library also includes composable and extensible algorithms for optim
 """
 module ITensors
 
+#if isdefined(Base, :Experimental) && isdefined(Base.Experimental, Symbol("@compiler_options"))
+#  @eval Base.Experimental.@compiler_options compile=min optimize=1 infer=false
+#end
+
 #####################################
 # External packages
 #
diff --git a/src/mps/mpo.jl b/src/mps/mpo.jl
index 0ab0cce1ca..1e6bdb243c 100644
--- a/src/mps/mpo.jl
+++ b/src/mps/mpo.jl
@@ -457,8 +457,15 @@ function _contract_densitymatrix(A::MPO, ψ::MPS; kwargs...)::MPS
     error("In `contract(A::MPO, x::MPS)`, `A` and `x` must share a set of site indices")
 
   # In case A and ψ have the same link indices
+  # For example, ψ could be A projected onto a product state,
+  # like `|ψ⟩ = A|0⟩`.
   A = sim(linkinds, A)
 
+  # Dangling indices of ψ. For now, only support dangling
+  # indices on the edges of the system.
+  sψ = siteinds(uniqueinds, ψ, A)
+  any(n -> n ≠ firstindex(sψ) && n ≠ lastindex(sψ) && !isempty(sψ[n]), eachindex(sψ)) && error("In `contract(A::MPO, x::MPS)`, `x` can only have dangling indices on the first and last tensor.")
+
   ψ_c = dag(ψ)
   A_c = dag(A)
 
@@ -467,13 +474,17 @@ function _contract_densitymatrix(A::MPO, ψ::MPS; kwargs...)::MPS
   sim!(linkinds, ψ_c)
   sim!(siteinds, commoninds, A_c, ψ_c)
 
+  # Prime any dangling indices on the ends of ψ
+  prime!(siteinds, uniqueinds, ψ_c, A_c)
+
   # A version helpful for making the density matrix
   simA_c = sim(siteinds, uniqueinds, A_c, ψ_c)
 
   # Store the left environment tensors
   E = Vector{ITensor}(undef, n - 1)
-
-  E[1] = ψ[1] * A[1] * A_c[1] * ψ_c[1]
+  sψ1 = sψ[1]
+  sψ1_c = sψ1'
+  E[1] = ψ[1] * A[1] * A_c[1] * replaceinds(ψ_c[1], sψ1_c => sψ1)
   for j in 2:(n - 1)
     E[j] = E[j - 1] * ψ[j] * A[j] * A_c[j] * ψ_c[j]
   end
@@ -481,11 +492,15 @@ function _contract_densitymatrix(A::MPO, ψ::MPS; kwargs...)::MPS
   simR_c = ψ_c[n] * simA_c[n]
   ρ = E[n - 1] * R * simR_c
   l = linkind(ψ, n - 1)
-  ts = isnothing(l) ? "" : tags(l)
-  Lis = siteinds(uniqueinds, A, ψ, n)
-  Ris = siteinds(uniqueinds, simA_c, ψ_c, n)
-  F = eigen(ρ, Lis, Ris; ishermitian=true, tags=ts, kwargs...)
+  tsₙ = isnothing(l) ? "" : tags(l)
+  plₙ = isnothing(l) ? 0 : plev(l)
+  sψn = sψ[n]
+  sψn_c = sψn'
+  Lis = unioninds(siteinds(uniqueinds, A, ψ, n), sψn)
+  Ris = unioninds(siteinds(uniqueinds, simA_c, ψ_c, n), sψn_c)
+  F = eigen(ρ, Lis, Ris; ishermitian=true, tags=tsₙ, plev=plₙ, kwargs...)
   D, U, Ut = F.D, F.V, F.Vt
+  Ut = setprime(Ut, plₙ; inds=commoninds(Ut, D))
   l_renorm, r_renorm = F.l, F.r
   ψ_out[n] = Ut
   R = R * dag(Ut) * ψ[n - 1] * A[n - 1]
@@ -495,11 +510,13 @@ function _contract_densitymatrix(A::MPO, ψ::MPS; kwargs...)::MPS
     s̃ = siteinds(uniqueinds, simA_c, ψ_c, j)
     ρ = E[j - 1] * R * simR_c
     l = linkind(ψ, j - 1)
-    ts = isnothing(l) ? "" : tags(l)
-    Lis = IndexSet(s..., l_renorm)
-    Ris = IndexSet(s̃..., r_renorm)
-    F = eigen(ρ, Lis, Ris; ishermitian=true, tags=ts, kwargs...)
+    tsⱼ = isnothing(l) ? "" : tags(l)
+    plⱼ = isnothing(l) ? 0 : plev(l)
+    Lis = [s..., l_renorm]
+    Ris = [s̃..., r_renorm]
+    F = eigen(ρ, Lis, Ris; ishermitian=true, tags=tsⱼ, plev=plⱼ, kwargs...)
     D, U, Ut = F.D, F.V, F.Vt
+    Ut = setprime(Ut, plₙ; inds=commoninds(Ut, D))
     l_renorm, r_renorm = F.l, F.r
     ψ_out[j] = Ut
     R = R * dag(Ut) * ψ[j - 1] * A[j - 1]
diff --git a/src/tagset.jl b/src/tagset.jl
index cd14725dbd..6618df38ce 100644
--- a/src/tagset.jl
+++ b/src/tagset.jl
@@ -56,11 +56,7 @@ end
 function _addtag!(ts::MTagSetStorage, ntags::Int, tag::IntTag)
   t = Tag(tag)
   if !isnull(t)
-    if isint(t)
-      error("Cannot use a bare integer as a tag.")
-    else
-      ntags = _addtag_ordered!(ts, ntags, tag)
-    end
+    ntags = _addtag_ordered!(ts, ntags, tag)
   end
   return ntags
 end

From 189a0a872507c0d75d268cbb1e52ae707aac50c3 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Tue, 25 May 2021 13:30:26 -0400
Subject: [PATCH 11/12] Fix bug in replaceinds for empty indices being replaced
 which fixes a new bug in contract MPO-MPS. Fix link tags from  to  for
 ITensor -> MPS/MPO constructor. Make MPO-MPS contraction return and MPS with
 prime level 0

---
 src/indexset.jl        | 33 +++++++++++++++++++--------------
 src/mps/abstractmps.jl |  6 +++---
 src/mps/mpo.jl         |  5 +++--
 3 files changed, 25 insertions(+), 19 deletions(-)

diff --git a/src/indexset.jl b/src/indexset.jl
index 23072dc4d7..215bb75d54 100644
--- a/src/indexset.jl
+++ b/src/indexset.jl
@@ -456,28 +456,33 @@ function swaptags(is::Indices, tags1, tags2, args...; kwargs...)
   return swaptags(fmatch(args...; kwargs...), is, tags1, tags2)
 end
 
-function replaceinds(is::Indices, rep_inds::Pair{<:Index,<:Index}...)
-  return replaceinds(is, zip(rep_inds...)...)
-end
-
-function replaceinds(is::Indices, rep_inds::Vector{<:Pair{<:Index,<:Index}})
-  return replaceinds(is, rep_inds...)
-end
-
-function replaceinds(is::Indices, rep_inds::Tuple{Vararg{Pair{<:Index,<:Index}}})
-  return replaceinds(is, rep_inds...)
-end
+# Check that the QNs are all the same
+hassameflux(i1::Index, i2::Index) = (dim(i1) == dim(i2))
 
+# For notation:
+# `replaceinds(is, [i, j] => [k, l])`
+# `replaceinds(is, (i, j) => (k, l))`
 function replaceinds(is::Indices, rep_inds::Pair)
-  return replaceinds(is, Tuple(first(rep_inds)) .=> Tuple(last(rep_inds)))
+  return replaceinds(is, first(rep_inds), last(rep_inds))
 end
 
-# Check that the QNs are all the same
-hassameflux(i1::Index, i2::Index) = (dim(i1) == dim(i2))
+# For notation:
+# `replaceinds(is, [i => k, j => l])`
+# `replaceinds(is, (i => k, j => l))`
+function replaceinds(is::Indices, inds1_inds2)
+  @assert all(x -> x isa Pair, inds1_inds2)
+  return replaceinds(is, first.(inds1_inds2), last.(inds1_inds2))
+end
 
+# For notation:
+# `replaceinds(is, [i, j], [k, l])`
+# `replaceinds(is, (i, j), (k, l))`
 function replaceinds(is::Indices, inds1, inds2)
+  @assert length(inds1) == length(inds2)
+  isempty(inds1) && return is
   is1 = (inds1)
   poss = indexin(is1, is)
+  # TODO: don't convert to Tuple in the case of a Vector
   is_tuple = Tuple(is)
   for (j, pos) in enumerate(poss)
     isnothing(pos) && continue
diff --git a/src/mps/abstractmps.jl b/src/mps/abstractmps.jl
index e89eb37464..0565a0c793 100644
--- a/src/mps/abstractmps.jl
+++ b/src/mps/abstractmps.jl
@@ -813,7 +813,7 @@ function hassamenuminds(::typeof(siteinds), M1::AbstractMPS, M2::AbstractMPS)
 end
 
 for fname in
-    (:sim, :prime, :setprime, :noprime, :addtags, :removetags, :replacetags, :settags)
+    (:sim, :prime, :setprime, :noprime, :replaceprime, :swapprime, :addtags, :removetags, :replacetags, :settags)
   fname! = Symbol(fname, :!)
   @eval begin
     """
@@ -1577,7 +1577,7 @@ function (::Type{MPST})(
   # Sweep left to orthocenter
   for n in 1:(orthocenter - 1)
     Lis = unioninds(sites[n], l)
-    L, R = factorize(Ã, Lis; kwargs..., tags="Link,n=$n", ortho="left")
+    L, R = factorize(Ã, Lis; kwargs..., tags="Link,l=$n", ortho="left")
     l = commoninds(L, R)
     ψ[n] = L
     Ã = R
@@ -1586,7 +1586,7 @@ function (::Type{MPST})(
   # Sweep right to orthocenter
   for n in reverse((orthocenter + 1):N)
     Ris = unioninds(sites[n], r)
-    R, L = factorize(Ã, Ris; kwargs..., tags="Link,n=$n", ortho="left")
+    R, L = factorize(Ã, Ris; kwargs..., tags="Link,l=$n", ortho="left")
     r = commonind(R, L)
     ψ[n] = R
     Ã = L
diff --git a/src/mps/mpo.jl b/src/mps/mpo.jl
index 1e6bdb243c..0f12311b8b 100644
--- a/src/mps/mpo.jl
+++ b/src/mps/mpo.jl
@@ -500,7 +500,6 @@ function _contract_densitymatrix(A::MPO, ψ::MPS; kwargs...)::MPS
   Ris = unioninds(siteinds(uniqueinds, simA_c, ψ_c, n), sψn_c)
   F = eigen(ρ, Lis, Ris; ishermitian=true, tags=tsₙ, plev=plₙ, kwargs...)
   D, U, Ut = F.D, F.V, F.Vt
-  Ut = setprime(Ut, plₙ; inds=commoninds(Ut, D))
   l_renorm, r_renorm = F.l, F.r
   ψ_out[n] = Ut
   R = R * dag(Ut) * ψ[n - 1] * A[n - 1]
@@ -516,7 +515,6 @@ function _contract_densitymatrix(A::MPO, ψ::MPS; kwargs...)::MPS
     Ris = [s̃..., r_renorm]
     F = eigen(ρ, Lis, Ris; ishermitian=true, tags=tsⱼ, plev=plⱼ, kwargs...)
     D, U, Ut = F.D, F.V, F.Vt
-    Ut = setprime(Ut, plₙ; inds=commoninds(Ut, D))
     l_renorm, r_renorm = F.l, F.r
     ψ_out[j] = Ut
     R = R * dag(Ut) * ψ[j - 1] * A[j - 1]
@@ -528,6 +526,9 @@ function _contract_densitymatrix(A::MPO, ψ::MPS; kwargs...)::MPS
   ψ_out[1] = R
   setleftlim!(ψ_out, 0)
   setrightlim!(ψ_out, 2)
+  # TODO: set the prime levels to those
+  # of the link indices of ψ
+  setprime!(linkinds, ψ_out, 0)
   return ψ_out
 end
 

From 30ea283e085d97b75b7dbc3d6da9ff922fccd990 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Tue, 25 May 2021 15:59:28 -0400
Subject: [PATCH 12/12] Fix some tests

---
 src/indexset.jl | 12 ++++++++++
 test/apply.jl   | 62 ++++++++++++++++++++++++-------------------------
 2 files changed, 43 insertions(+), 31 deletions(-)

diff --git a/src/indexset.jl b/src/indexset.jl
index 215bb75d54..6a721b4b99 100644
--- a/src/indexset.jl
+++ b/src/indexset.jl
@@ -466,6 +466,18 @@ function replaceinds(is::Indices, rep_inds::Pair)
   return replaceinds(is, first(rep_inds), last(rep_inds))
 end
 
+# For notation:
+# `replaceinds(is, i => k, j => l)`
+function replaceinds(is::Indices, rep_inds::Pair{<:Index,<:Index}...)
+  return replaceinds(is, first.(rep_inds), last.(rep_inds))
+end
+
+# For notation:
+# `replaceinds(is, i, k)`
+function replaceinds(is::Indices, inds1::Index, inds2::Index)
+  return replaceinds(is, (inds1,), (inds2,))
+end
+
 # For notation:
 # `replaceinds(is, [i => k, j => l])`
 # `replaceinds(is, (i => k, j => l))`
diff --git a/test/apply.jl b/test/apply.jl
index f95ffd69f8..a91119bc07 100644
--- a/test/apply.jl
+++ b/test/apply.jl
@@ -4,34 +4,34 @@ using Test
 @testset "apply" begin
 
   @testset "ITensors.minimal_swap_range" begin
-    @test ITensors.minimal_swap_range([1, 2], [1, 2]) == 1:2
-    @test ITensors.minimal_swap_range([1, 2], [3, 4]) == 1:2
-    @test ITensors.minimal_swap_range([1, 3], [5, 6]) == 2:3
-    @test ITensors.minimal_swap_range([1, 4], [5, 6]) == 3:4
-    @test ITensors.minimal_swap_range([3, 6], [1, 4]) == 4:5
-    @test ITensors.minimal_swap_range([1, 5], [5, 6]) == 4:5
-    @test ITensors.minimal_swap_range([1, 6], [5, 6]) == 5:6
-    @test ITensors.minimal_swap_range([1, 6], [4, 5, 6]) == 4:5
-    @test ITensors.minimal_swap_range([1, 3, 5], [7, 8]) == 3:5
-    @test ITensors.minimal_swap_range([1, 5], [2, 3]) == 2:3
-    @test ITensors.minimal_swap_range([1, 5], [2, 4]) == 2:3
-    @test ITensors.minimal_swap_range([1, 3, 5], [2, 4]) == 2:4
-    @test ITensors.minimal_swap_range([1, 3, 5], [5, 8]) == 3:5
-    @test ITensors.minimal_swap_range([3, 4], [1, 2]) == 3:4
-    @test ITensors.minimal_swap_range([4, 6], [1, 2]) == 4:5
-    @test ITensors.minimal_swap_range([4, 6], [1, 3]) == 4:5
-    @test ITensors.minimal_swap_range([4, 6], [1, 4]) == 4:5
-    @test ITensors.minimal_swap_range([2, 6], [1, 6]) == 2:3
-    @test ITensors.minimal_swap_range([2, 6], [1, 7]) == 2:3
-    @test ITensors.minimal_swap_range([2, 6], [1, 5]) == 5:6
-    @test ITensors.minimal_swap_range([1, 5], [2, 6]) == 1:2
-    @test ITensors.minimal_swap_range([1, 3, 5], [2, 6]) == 1:3
-    @test ITensors.minimal_swap_range([2, 5, 6], [1, 5]) == 4:6
-    @test ITensors.minimal_swap_range([1, 5], [2, 6]) == 1:2
-    @test ITensors.minimal_swap_range([4, 8], [3, 8]) == 4:5
-    @test ITensors.minimal_swap_range([4, 6, 8], [3, 5]) == 4:6
-    @test ITensors.minimal_swap_range([4, 6, 8], [3, 4]) == 4:6
-    @test ITensors.minimal_swap_range([4, 6, 8], [3, 8]) == 4:6
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [1, 2], [1, 2]) == 1:2
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [1, 2], [3, 4]) == 1:2
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [1, 3], [5, 6]) == 2:3
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [1, 4], [5, 6]) == 3:4
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [3, 6], [1, 4]) == 4:5
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [1, 5], [5, 6]) == 4:5
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [1, 6], [5, 6]) == 5:6
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [1, 6], [4, 5, 6]) == 4:5
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [1, 3, 5], [7, 8]) == 3:5
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [1, 5], [2, 3]) == 2:3
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [1, 5], [2, 4]) == 2:3
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [1, 3, 5], [2, 4]) == 2:4
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [1, 3, 5], [5, 8]) == 3:5
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [3, 4], [1, 2]) == 3:4
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [4, 6], [1, 2]) == 4:5
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [4, 6], [1, 3]) == 4:5
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [4, 6], [1, 4]) == 4:5
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [2, 6], [1, 6]) == 2:3
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [2, 6], [1, 7]) == 2:3
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [2, 6], [1, 5]) == 5:6
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [1, 5], [2, 6]) == 1:2
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [1, 3, 5], [2, 6]) == 1:3
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [2, 5, 6], [1, 5]) == 4:6
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [1, 5], [2, 6]) == 1:2
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [4, 8], [3, 8]) == 4:5
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [4, 6, 8], [3, 5]) == 4:6
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [4, 6, 8], [3, 4]) == 4:6
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [4, 6, 8], [3, 8]) == 4:6
   end
 
   @testset "apply swap" begin
@@ -85,12 +85,12 @@ using Test
 
     nswaps = Int[]
     ψ = apply(CX[1, 3], ψ0; (nswaps!)=nswaps, kwargs...)
-    @test ortho_lims(ψ) == 2:2
+    @test ortho_lims(ψ) == 3:3
     @test nswaps == [1]
 
     nswaps = Int[]
     ψ = apply(CX[1, 4], ψ0; (nswaps!)=nswaps, kwargs...)
-    @test ortho_lims(ψ) == 2:2
+    @test ortho_lims(ψ) == 4:4
     @test nswaps == [2]
 
     nswaps = Int[]
@@ -109,7 +109,7 @@ using Test
     ψ = apply([CX[3, 6], CX[1, 3]], ψ0; (nswaps!)=nswaps, kwargs...)
     @test nswaps == [2, 1]
 
-    @test ITensors.minimal_swap_range([3, 6], [1, 4]) == 4:5
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [3, 6], [1, 4]) == 4:5
     nswaps = Int[]
     ψ = apply([CX[3, 6], CX[1, 4]], ψ0; (nswaps!)=nswaps, kwargs...)
     @test nswaps == [2, 1]