From e4b9ff86dd232a595253f88a3d337f3d9db6afcf Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Sergio=20S=C3=A1nchez=20Ram=C3=ADrez?=
 <sergio.sanchez.ramirez+git@bsc.es>
Date: Wed, 17 Apr 2024 17:12:19 +0200
Subject: [PATCH] Fix truncation of 0s on `evolve!`

---
 src/Ansatz/Chain.jl | 24 ++++++++++--------------
 1 file changed, 10 insertions(+), 14 deletions(-)

diff --git a/src/Ansatz/Chain.jl b/src/Ansatz/Chain.jl
index ce7566a..ded195b 100644
--- a/src/Ansatz/Chain.jl
+++ b/src/Ansatz/Chain.jl
@@ -318,7 +318,7 @@ Truncate the dimension of the virtual `bond`` of the [`Chain`](@ref) Tensor Netw
   - Either `threshold` or `maxdim` must be provided. If both are provided, `maxdim` is used.
   - The bond must contain the Schmidt coefficients, i.e. a site canonization must be performed before calling `truncate!`.
 """
-function truncate!(qtn::Chain, bond; threshold::Union{Nothing,Real} = nothing, maxdim::Union{Nothing,Int} = nothing)
+function truncate!(qtn::Chain, bond; threshold::Real = 1e-16, maxdim::Union{Nothing,Int} = nothing)
     # TODO replace for select(:between)
     vind = rightindex(qtn, bond[1])
     if vind != leftindex(qtn, bond[2])
@@ -332,12 +332,15 @@ function truncate!(qtn::Chain, bond; threshold::Union{Nothing,Real} = nothing, m
     tensor = TensorNetwork(qtn)[vind]
     spectrum = parent(tensor)
 
-    extent = if !isnothing(maxdim)
-        1:min(size(TensorNetwork(qtn), vind), maxdim)
-    elseif !isnothing(threshold)
-        findall(>(threshold) ∘ abs, spectrum)
+    extent = if isnothing(maxdim)
+        1:size(size(TensorNetwork(qtn), vind))
     else
-        throw(ArgumentError("Either `threshold` or `maxdim` must be provided"))
+        1:min(size(TensorNetwork(qtn), vind), maxdim)
+    end |> collect
+
+    # remove 0s from spectrum
+    filter!(extent) do i
+        abs(spectrum[i]) > threshold
     end
 
     slice!(TensorNetwork(qtn), vind, extent)
@@ -461,14 +464,7 @@ end
 
 Applies a local operator `gate` to the [`Chain`](@ref) tensor network.
 """
-function evolve!(
-    qtn::Chain,
-    gate::Dense;
-    threshold = nothing,
-    maxdim = nothing,
-    iscanonical = false,
-    renormalize = false,
-)
+function evolve!(qtn::Chain, gate::Dense; threshold = 1e-16, maxdim = nothing, iscanonical = false, renormalize = false)
     # check gate is a valid operator
     if !(socket(gate) isa Operator)
         throw(ArgumentError("Gate must be an operator, but got $(socket(gate))"))