From 5a257c8ae75552a8f76b181b2046b2dd0e94e8e4 Mon Sep 17 00:00:00 2001
From: Ben Corbett <corbett5@llnl.gov>
Date: Thu, 26 Sep 2024 18:22:24 -0600
Subject: [PATCH 01/10] Added option for MPO-MPS zipup.

---
 src/lib/ITensorMPS/src/mpo.jl | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/lib/ITensorMPS/src/mpo.jl b/src/lib/ITensorMPS/src/mpo.jl
index 647ca1bad6..12a41a5f9b 100644
--- a/src/lib/ITensorMPS/src/mpo.jl
+++ b/src/lib/ITensorMPS/src/mpo.jl
@@ -822,7 +822,7 @@ end
 function ITensors.contract(
   ::Algorithm"zipup",
   A::MPO,
-  B::MPO;
+  B::AbstractMPS;
   cutoff=1e-14,
   maxdim=maxlinkdim(A) * maxlinkdim(B),
   mindim=1,
@@ -836,14 +836,15 @@ function ITensors.contract(
   N = length(A)
   N != length(B) &&
     throw(DimensionMismatch("lengths of MPOs A ($N) and B ($(length(B))) do not match"))
+  ResultType = typeof(B)
   # Special case for a single site
-  N == 1 && return MPO([A[1] * B[1]])
+  N == 1 && return ResultType([A[1] * B[1]])
   A = orthogonalize(A, 1)
   B = orthogonalize(B, 1)
   A = sim(linkinds, A)
   sA = siteinds(uniqueinds, A, B)
   sB = siteinds(uniqueinds, B, A)
-  C = MPO(N)
+  C = ResultType(N)
   lCᵢ = Index[]
   R = ITensor(true)
   for i in 1:(N - 2)
@@ -874,7 +875,6 @@ function ITensors.contract(
     mindim,
     kwargs...,
   )
-  truncate!(C; kwargs...)
   return C
 end
 

From ec126bab4cc37cd71062990ca8d639167f32138a Mon Sep 17 00:00:00 2001
From: Ben Corbett <corbett5@llnl.gov>
Date: Tue, 8 Oct 2024 15:31:33 -0600
Subject: [PATCH 02/10] Removed ResultType.

---
 src/lib/ITensorMPS/src/mpo.jl | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/src/lib/ITensorMPS/src/mpo.jl b/src/lib/ITensorMPS/src/mpo.jl
index 12a41a5f9b..637f79470b 100644
--- a/src/lib/ITensorMPS/src/mpo.jl
+++ b/src/lib/ITensorMPS/src/mpo.jl
@@ -836,15 +836,14 @@ function ITensors.contract(
   N = length(A)
   N != length(B) &&
     throw(DimensionMismatch("lengths of MPOs A ($N) and B ($(length(B))) do not match"))
-  ResultType = typeof(B)
   # Special case for a single site
-  N == 1 && return ResultType([A[1] * B[1]])
+  N == 1 && return typeof(B)([A[1] * B[1]])
   A = orthogonalize(A, 1)
   B = orthogonalize(B, 1)
   A = sim(linkinds, A)
   sA = siteinds(uniqueinds, A, B)
   sB = siteinds(uniqueinds, B, A)
-  C = ResultType(N)
+  C = typeof(B)(N)
   lCᵢ = Index[]
   R = ITensor(true)
   for i in 1:(N - 2)

From 70404628211240dd5c17ef353b393aff450ac757 Mon Sep 17 00:00:00 2001
From: Ben Corbett <corbett5@llnl.gov>
Date: Wed, 9 Oct 2024 11:09:02 -0600
Subject: [PATCH 03/10] Added truncate_kwargs parameter

---
 src/lib/ITensorMPS/src/mpo.jl | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/lib/ITensorMPS/src/mpo.jl b/src/lib/ITensorMPS/src/mpo.jl
index 637f79470b..2ceef25972 100644
--- a/src/lib/ITensorMPS/src/mpo.jl
+++ b/src/lib/ITensorMPS/src/mpo.jl
@@ -826,6 +826,7 @@ function ITensors.contract(
   cutoff=1e-14,
   maxdim=maxlinkdim(A) * maxlinkdim(B),
   mindim=1,
+  truncate_kwargs=(;),
   kwargs...,
 )
   if hassameinds(siteinds, A, B)
@@ -874,6 +875,7 @@ function ITensors.contract(
     mindim,
     kwargs...,
   )
+  truncate!(C, truncate_kwargs...)
   return C
 end
 

From 74f452a60bfcb5aa8465ae87d672f6805dcec07a Mon Sep 17 00:00:00 2001
From: Ben Corbett <corbett5@llnl.gov>
Date: Wed, 9 Oct 2024 22:31:27 -0600
Subject: [PATCH 04/10] Changed truncate_kwargs default value.

---
 src/lib/ITensorMPS/src/mpo.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/lib/ITensorMPS/src/mpo.jl b/src/lib/ITensorMPS/src/mpo.jl
index 2ceef25972..fea89e39e8 100644
--- a/src/lib/ITensorMPS/src/mpo.jl
+++ b/src/lib/ITensorMPS/src/mpo.jl
@@ -826,7 +826,7 @@ function ITensors.contract(
   cutoff=1e-14,
   maxdim=maxlinkdim(A) * maxlinkdim(B),
   mindim=1,
-  truncate_kwargs=(;),
+  truncate_kwargs=(; cutoff, maxdim, mindim),
   kwargs...,
 )
   if hassameinds(siteinds, A, B)

From 07f7534785a4b657bfb90acc2e4e26674c793867 Mon Sep 17 00:00:00 2001
From: Ben Corbett <corbett5@llnl.gov>
Date: Thu, 10 Oct 2024 00:25:54 -0600
Subject: [PATCH 05/10] Fixed mistake with truncate kwargs.

---
 src/lib/ITensorMPS/src/mpo.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/lib/ITensorMPS/src/mpo.jl b/src/lib/ITensorMPS/src/mpo.jl
index fea89e39e8..362b1364eb 100644
--- a/src/lib/ITensorMPS/src/mpo.jl
+++ b/src/lib/ITensorMPS/src/mpo.jl
@@ -875,7 +875,7 @@ function ITensors.contract(
     mindim,
     kwargs...,
   )
-  truncate!(C, truncate_kwargs...)
+  truncate!(C; truncate_kwargs...)
   return C
 end
 

From a332ac986d8a666a0ded1e2c34c1ae11f9e344dd Mon Sep 17 00:00:00 2001
From: Ben Corbett <corbett5@llnl.gov>
Date: Tue, 15 Oct 2024 15:15:16 -0600
Subject: [PATCH 06/10] Updated docs and added tests.

---
 src/lib/ITensorMPS/src/mpo.jl            | 27 +++++++++++++++++-------
 src/lib/ITensorMPS/test/base/test_mpo.jl | 22 +++++++++----------
 2 files changed, 30 insertions(+), 19 deletions(-)

diff --git a/src/lib/ITensorMPS/src/mpo.jl b/src/lib/ITensorMPS/src/mpo.jl
index 362b1364eb..e3e1cf88db 100644
--- a/src/lib/ITensorMPS/src/mpo.jl
+++ b/src/lib/ITensorMPS/src/mpo.jl
@@ -662,11 +662,16 @@ Choose the method with the `method` keyword, for example
 - `mindim::Int=1`: the minimal bond dimension of the resulting MPS.
 - `normalize::Bool=false`: whether or not to normalize the resulting MPS.
 - `method::String="densitymatrix"`: the algorithm to use for the contraction.
-   Currently the options are "densitymatrix", where the network formed by the
-   MPO and MPS is squared and contracted down to a density matrix which is
-   diagonalized iteratively at each site, and "naive", where the MPO and MPS
-   tensor are contracted exactly at each site and then a truncation of the
-   resulting MPS is performed.
+    - "densitymatrix": The network formed by the MPO and MPS is squared and contracted down to
+    a density matrix which is diagonalized iteratively at each site.
+    - "naive": The MPO and MPS tensor are contracted exactly at each site and then a truncation 
+    of the resulting MPS is performed.
+    - "zipup": The MPO and MPS tensors are contracted then truncated at each site without enforcing
+    the appropriate orthogonal gauge. Once this sweep is complete a call to `truncate!` occurs.
+    Because the initial truncation is not locally optimal it is recommended to use a loose
+    `cutoff` and `maxdim` and then pass the desired truncation parameters to the locally optimal
+    `truncate!` sweep via the additional keyword argument `truncate_kwargs`.
+    Suggested use is `contract(A, ψ; method="zipup", cutoff=cutoff / 10, maxdim=2 * maxdim, truncate_kwargs=(; cutoff, maxdim))`.
 
 See also [`apply`](@ref).
 """
@@ -946,14 +951,20 @@ C = apply(A, B; alg="naive", truncate=false)
    in general you should set a `cutoff` value.
 - `maxdim::Int=maxlinkdim(A) * maxlinkdim(B))`: the maximal bond dimension of the results MPS.
 - `mindim::Int=1`: the minimal bond dimension of the resulting MPS.
-- `alg="zipup"`: Either `"zipup"` or `"naive"`. `"zipup"` contracts pairs of
-   site tensors and truncates with SVDs in a sweep across the sites, while `"naive"`
-   first contracts pairs of tensor exactly and then truncates at the end if `truncate=true`.
 - `truncate=true`: Enable or disable truncation. If `truncate=false`, ignore
    other truncation parameters like `cutoff` and `maxdim`. This is most relevant
    for the `"naive"` version, if you just want to contract the tensors pairwise
    exactly. This can be useful if you are contracting MPOs that have diverging
    norms, such as MPOs originating from sums of local operators.
+- `alg="zipup"`: the algorithm to use for the contraction. Supported algorithms are
+    - "naive": The MPO tensors are contracted exactly at each site and then a truncation 
+    of the resulting MPO is performed.
+    - "zipup": The MPO and MPS tensors are contracted then truncated at each site without enforcing
+    the appropriate orthogonal gauge. Once this sweep is complete a call to `truncate!` occurs.
+    Because the initial truncation is not locally optimal it is recommended to use a loose
+    `cutoff` and `maxdim` and then pass the desired truncation parameters to the locally optimal
+    `truncate!` sweep via the additional keyword argument `truncate_kwargs`.
+    Suggested use is `contract(A, ψ; method="zipup", cutoff=cutoff / 10, maxdim=2 * maxdim, truncate_kwargs=(; cutoff, maxdim))`.
 
 See also [`apply`](@ref) for details about the arguments available.
 """
diff --git a/src/lib/ITensorMPS/test/base/test_mpo.jl b/src/lib/ITensorMPS/test/base/test_mpo.jl
index 32746a752b..005b95a205 100644
--- a/src/lib/ITensorMPS/test/base/test_mpo.jl
+++ b/src/lib/ITensorMPS/test/base/test_mpo.jl
@@ -240,14 +240,14 @@ end
     @test_throws DimensionMismatch error_contract(phi, K, badpsi)
   end
 
-  @testset "contract" begin
+  @testset "contract(::MPO, ::MPS)" for method in ["densitymatrix", "naive", "zipup"]
     phi = random_mps(sites)
     K = random_mpo(sites)
     @test maxlinkdim(K) == 1
     psi = random_mps(sites)
-    psi_out = contract(K, psi; maxdim=1)
+    psi_out = contract(K, psi; method, maxdim=1)
     @test inner(phi', psi_out) ≈ inner(phi', K, psi)
-    psi_out = contract(psi, K; maxdim=1)
+    psi_out = contract(psi, K; method, maxdim=1)
     @test inner(phi', psi_out) ≈ inner(phi', K, psi)
     psi_out = psi * K
     @test inner(phi', psi_out) ≈ inner(phi', K, psi)
@@ -255,7 +255,7 @@ end
 
     badsites = [Index(2, "Site") for n in 1:(N + 1)]
     badpsi = random_mps(badsites)
-    @test_throws DimensionMismatch contract(K, badpsi)
+    @test_throws DimensionMismatch contract(K, badpsi; method)
 
     # make bigger random MPO...
     for link_dim in 2:5
@@ -287,7 +287,7 @@ end
       orthogonalize!(psi, 1; maxdim=link_dim)
       orthogonalize!(K, 1; maxdim=link_dim)
       orthogonalize!(phi, 1; normalize=true, maxdim=link_dim)
-      psi_out = contract(deepcopy(K), deepcopy(psi); maxdim=10 * link_dim, cutoff=0.0)
+      psi_out = contract(deepcopy(K), deepcopy(psi); method, maxdim=10 * link_dim, cutoff=0.0)
       @test inner(phi', psi_out) ≈ inner(phi', K, psi)
     end
   end
@@ -354,14 +354,14 @@ end
     @test maxlinkdim(H) ≤ maxlinkdim(H₁) + maxlinkdim(H₂)
   end
 
-  @testset "contract(::MPO, ::MPO)" begin
+  @testset "contract(::MPO, ::MPO)" for alg in ["naive", "zipup"]
     psi = random_mps(sites)
     K = random_mpo(sites)
     L = random_mpo(sites)
     @test maxlinkdim(K) == 1
     @test maxlinkdim(L) == 1
-    KL = contract(prime(K), L; maxdim=1)
-    psi_kl_out = contract(prime(K), contract(L, psi; maxdim=1); maxdim=1)
+    KL = contract(prime(K), L; alg, maxdim=1)
+    psi_kl_out = contract(prime(K), contract(L, psi; alg, maxdim=1); alg, maxdim=1)
     @test inner(psi'', KL, psi) ≈ inner(psi'', psi_kl_out) atol = 5e-3
 
     # where both K and L have differently labelled sites
@@ -373,15 +373,15 @@ end
       replaceind!(K[ii], sites[ii]', othersitesk[ii])
       replaceind!(L[ii], sites[ii]', othersitesl[ii])
     end
-    KL = contract(K, L; maxdim=1)
+    KL = contract(K, L; alg, maxdim=1)
     psik = random_mps(othersitesk)
     psil = random_mps(othersitesl)
-    psi_kl_out = contract(K, contract(L, psil; maxdim=1); maxdim=1)
+    psi_kl_out = contract(K, contract(L, psil; alg, maxdim=1); alg, maxdim=1)
     @test inner(psik, KL, psil) ≈ inner(psik, psi_kl_out) atol = 5e-3
 
     badsites = [Index(2, "Site") for n in 1:(N + 1)]
     badL = random_mpo(badsites)
-    @test_throws DimensionMismatch contract(K, badL)
+    @test_throws DimensionMismatch contract(K, badL; alg)
   end
 
   @testset "*(::MPO, ::MPO)" begin

From e2979a12f8f2ca7d5e189473bac26de1a3037a90 Mon Sep 17 00:00:00 2001
From: Ben Corbett <corbett5@llnl.gov>
Date: Thu, 17 Oct 2024 11:54:10 -0600
Subject: [PATCH 07/10] Format changes.

---
 src/lib/ITensorMPS/src/abstractmps.jl    | 2 +-
 src/lib/ITensorMPS/test/base/test_mpo.jl | 4 +++-
 2 files changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/lib/ITensorMPS/src/abstractmps.jl b/src/lib/ITensorMPS/src/abstractmps.jl
index 629075f183..577415214d 100644
--- a/src/lib/ITensorMPS/src/abstractmps.jl
+++ b/src/lib/ITensorMPS/src/abstractmps.jl
@@ -1082,7 +1082,7 @@ function deprecate_make_inds_match!(
             inds(ψ[$n]) = $(inds(M2[n]))
 
         Make sure the site indices of your MPO/MPS match. You may need to prime
-        one of the MPS, such as `dot(ϕ', ψ)`."""
+        one of the MPS, such as `dot(ϕ', ψ)`.""",
       )
     end
     make_inds_match = false
diff --git a/src/lib/ITensorMPS/test/base/test_mpo.jl b/src/lib/ITensorMPS/test/base/test_mpo.jl
index 005b95a205..73e9ac0a73 100644
--- a/src/lib/ITensorMPS/test/base/test_mpo.jl
+++ b/src/lib/ITensorMPS/test/base/test_mpo.jl
@@ -287,7 +287,9 @@ end
       orthogonalize!(psi, 1; maxdim=link_dim)
       orthogonalize!(K, 1; maxdim=link_dim)
       orthogonalize!(phi, 1; normalize=true, maxdim=link_dim)
-      psi_out = contract(deepcopy(K), deepcopy(psi); method, maxdim=10 * link_dim, cutoff=0.0)
+      psi_out = contract(
+        deepcopy(K), deepcopy(psi); method, maxdim=10 * link_dim, cutoff=0.0
+      )
       @test inner(phi', psi_out) ≈ inner(phi', K, psi)
     end
   end

From 54387b92b3a4e2f8c1c62cdf977f258bac785d95 Mon Sep 17 00:00:00 2001
From: Ben Corbett <corbett5@llnl.gov>
Date: Thu, 17 Oct 2024 11:56:54 -0600
Subject: [PATCH 08/10] Manual format changes.

---
 src/lib/ITensorMPS/src/abstractmps.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/lib/ITensorMPS/src/abstractmps.jl b/src/lib/ITensorMPS/src/abstractmps.jl
index 577415214d..629075f183 100644
--- a/src/lib/ITensorMPS/src/abstractmps.jl
+++ b/src/lib/ITensorMPS/src/abstractmps.jl
@@ -1082,7 +1082,7 @@ function deprecate_make_inds_match!(
             inds(ψ[$n]) = $(inds(M2[n]))
 
         Make sure the site indices of your MPO/MPS match. You may need to prime
-        one of the MPS, such as `dot(ϕ', ψ)`.""",
+        one of the MPS, such as `dot(ϕ', ψ)`."""
       )
     end
     make_inds_match = false

From 08c6cd6379e9047c0992acf3f40ebab42ec34ae6 Mon Sep 17 00:00:00 2001
From: Ben Corbett <corbett5@llnl.gov>
Date: Wed, 23 Oct 2024 17:42:40 -0600
Subject: [PATCH 09/10] Combined documentation and added reference.

---
 src/lib/ITensorMPS/src/mpo.jl | 36 +++++++++++++++++++++++------------
 1 file changed, 24 insertions(+), 12 deletions(-)

diff --git a/src/lib/ITensorMPS/src/mpo.jl b/src/lib/ITensorMPS/src/mpo.jl
index e3e1cf88db..d882a4c15c 100644
--- a/src/lib/ITensorMPS/src/mpo.jl
+++ b/src/lib/ITensorMPS/src/mpo.jl
@@ -629,6 +629,24 @@ function ITensors.contract(A::MPO, ψ::MPS; alg=nothing, method=alg, kwargs...)
   return contract(Algorithm(alg), A, ψ; kwargs...)
 end
 
+function zipup_docstring(isMPOMPO::Bool)::String
+  rhsTypeString = isMPOMPO ? "MPO" : "MPS"
+  rhsString = isMPOMPO ? "B" : "ψ"
+  return """    - "zipup": The MPO and $rhsTypeString tensors are contracted then truncated at each site without enforcing
+  the appropriate orthogonal gauge. Once this sweep is complete a call to `truncate!` occurs.
+  Because the initial truncation is not locally optimal it is recommended to use a loose
+  `cutoff` and `maxdim` and then pass the desired truncation parameters to the locally optimal
+  `truncate!` sweep via the additional keyword argument `truncate_kwargs`.
+  A set of parameters suggested in [^Paeckel2019] is
+  `contract(A, $rhsString; method="zipup", cutoff=cutoff / 10, maxdim=2 * maxdim, truncate_kwargs=(; cutoff, maxdim))`."""
+end
+
+function Paeckel2019_citation_docstring()::String
+  return """
+  [^Paeckel2019]: Time-evolution methods for matrix-product states. Sebastian Paeckel et al. [arXiv:1901.05824](https://arxiv.org/abs/1901.05824)
+  """
+end
+
 contract_mpo_mps_doc = """
     contract(ψ::MPS, A::MPO; kwargs...) -> MPS
     *(::MPS, ::MPO; kwargs...) -> MPS
@@ -666,14 +684,11 @@ Choose the method with the `method` keyword, for example
     a density matrix which is diagonalized iteratively at each site.
     - "naive": The MPO and MPS tensor are contracted exactly at each site and then a truncation 
     of the resulting MPS is performed.
-    - "zipup": The MPO and MPS tensors are contracted then truncated at each site without enforcing
-    the appropriate orthogonal gauge. Once this sweep is complete a call to `truncate!` occurs.
-    Because the initial truncation is not locally optimal it is recommended to use a loose
-    `cutoff` and `maxdim` and then pass the desired truncation parameters to the locally optimal
-    `truncate!` sweep via the additional keyword argument `truncate_kwargs`.
-    Suggested use is `contract(A, ψ; method="zipup", cutoff=cutoff / 10, maxdim=2 * maxdim, truncate_kwargs=(; cutoff, maxdim))`.
+    $(zipup_docstring(false))
 
 See also [`apply`](@ref).
+
+$(Paeckel2019_citation_docstring())
 """
 
 @doc """
@@ -959,14 +974,11 @@ C = apply(A, B; alg="naive", truncate=false)
 - `alg="zipup"`: the algorithm to use for the contraction. Supported algorithms are
     - "naive": The MPO tensors are contracted exactly at each site and then a truncation 
     of the resulting MPO is performed.
-    - "zipup": The MPO and MPS tensors are contracted then truncated at each site without enforcing
-    the appropriate orthogonal gauge. Once this sweep is complete a call to `truncate!` occurs.
-    Because the initial truncation is not locally optimal it is recommended to use a loose
-    `cutoff` and `maxdim` and then pass the desired truncation parameters to the locally optimal
-    `truncate!` sweep via the additional keyword argument `truncate_kwargs`.
-    Suggested use is `contract(A, ψ; method="zipup", cutoff=cutoff / 10, maxdim=2 * maxdim, truncate_kwargs=(; cutoff, maxdim))`.
+    $(zipup_docstring(false))
 
 See also [`apply`](@ref) for details about the arguments available.
+
+$(Paeckel2019_citation_docstring())
 """
 
 @doc """

From f1f6cc2b034555597e30183673063c8432e0a1b6 Mon Sep 17 00:00:00 2001
From: Ben Corbett <corbett5@llnl.gov>
Date: Wed, 23 Oct 2024 18:15:20 -0600
Subject: [PATCH 10/10] Fixed typo in docs.

---
 src/lib/ITensorMPS/src/mpo.jl | 5 +----
 1 file changed, 1 insertion(+), 4 deletions(-)

diff --git a/src/lib/ITensorMPS/src/mpo.jl b/src/lib/ITensorMPS/src/mpo.jl
index d882a4c15c..301ec12091 100644
--- a/src/lib/ITensorMPS/src/mpo.jl
+++ b/src/lib/ITensorMPS/src/mpo.jl
@@ -971,10 +971,7 @@ C = apply(A, B; alg="naive", truncate=false)
    for the `"naive"` version, if you just want to contract the tensors pairwise
    exactly. This can be useful if you are contracting MPOs that have diverging
    norms, such as MPOs originating from sums of local operators.
-- `alg="zipup"`: the algorithm to use for the contraction. Supported algorithms are
-    - "naive": The MPO tensors are contracted exactly at each site and then a truncation 
-    of the resulting MPO is performed.
-    $(zipup_docstring(false))
+  $(zipup_docstring(true))
 
 See also [`apply`](@ref) for details about the arguments available.