Merge pull request #253 from LHerviou/Temp_ImprovingInfiniteMPOMatrix

Temp improving infinite mpo matrix

src/ITensors.jl

@@ -183,16 +183,12 @@ end
 # Handle orthogonality center correctly
 Base.getindex(ψ::MPS, r::UnitRange{Int}) = MPS([ψ[n] for n in r])
 
-# Base.zero(T::ITensor) = false * T
-#
-# function ITensors.NDTensors.datatype(
-#   ::ITensors.EmptyStorage{Float64,ITensors.NDTensors.BlockSparse{Float64,Vector{Float64},N}}
-# ) where {N}
-#   return Vector{Float64}
-# end
-#
-# function ITensors.NDTensors.datatype(T::ITensors.NDTensors.TensorStorage{Float64})
-#   return typeof(data(T))
-# end
+#TODO Remove if everything is working nicely
+#Was still crashing on my laptop after updating ITensors
+Base.fill!(::NDTensors.NoData, ::Any) = NDTensors.NoData()
 
-Base.fill!(::NDTensors.NoData, ::Any) = 0
+function ITensors.NDTensors.contraction_output(
+  A::NDTensors.EmptyTensor, B::NDTensors.DiagBlockSparseTensor, label
+)
+  return NDTensors.EmptyTensor(eltype(B), label)
+end

src/infinitempo.jl

@@ -85,3 +85,36 @@ function InfiniteMPO(Hm::InfiniteMPOMatrix)
   end
   return InfiniteMPO(Hi)
 end
+
+"""
+	fuse_legs!(Hcl::InfiniteMPO, L, R)
+
+    Fuse the non-site legs of the infiniteMPO Hcl and the corresponding left L and right R environments.
+		Essentially the inverse of splitblocks. It becomes useful for the very dense MPOs once get after compression sometimes.
+		Hcl is modified on site, L and R are not.
+		Input: Hcl the infinite MPO
+		       L   the left environment (an ITensor)
+					 R   the right environment (an ITensor)
+		Output: the tuple (newL, newR) the updated environments
+"""
+function fuse_legs!(Hcl::InfiniteMPO, L, R)
+  N = nsites(Hcl)
+  for j in 1:(N - 1)
+    right_link = only(commoninds(Hcl[j], Hcl[j + 1]))
+    comb = combiner(right_link; tags=tags(right_link), dir=dir(right_link))
+    comb_ind = combinedind(comb)
+    Hcl.data.data[j] = Hcl[j] * comb
+    Hcl.data.data[j + 1] = dag(comb) * Hcl[j + 1]
+  end
+  right_link = only(commoninds(Hcl[N], Hcl[N + 1]))
+  right_link2 = translatecell(fermion_momentum_translator, right_link, -1)
+  comb = combiner(right_link; tags=tags(right_link), dir=dir(right_link))
+  comb_ind = combinedind(comb)
+  comb2 = translatecell(fermion_momentum_translator, comb, -1)
+  comb2_ind = combinedind(comb2)
+  Hcl.data.data[N] = Hcl[N] * comb
+  Hcl.data.data[1] = dag(comb2) * Hcl[1]
+  L = L * comb2
+  R = dag(comb) * R
+  return L, R
+end

src/infinitempomatrix.jl

@@ -114,7 +114,15 @@ function convert_itensor_to_itensormatrix(tensor; kwargs...)
   end
 end
 
-"Build the projectors on the three parts of the itensor used to split a MPO into an InfiniteMPOMatrix"
+"""
+    build_three_projectors_from_index(is::Index; kwargs...)
+
+  Build the projectors on the three parts of the itensor used to split a MPO into an InfiniteMPOMatrix
+  More precisely, create projectors on the first dimension, the 2:end-1 and the last dimension of the index
+  Input: is the Index to split
+  Output: the triplet of projectors (first, middle, last)
+  Optional arguments: tags: if we want to change the tags of the index. Else default to tags(is)
+"""
 function build_three_projectors_from_index(is::Index; kwargs...)
   old_dim = dim(is)
   new_tags = get(kwargs, :tags, tags(is))
@@ -148,31 +156,50 @@ function build_three_projectors_from_index(is::Index; kwargs...)
   return top, middle, bottom
 end
 
-function convert_itensor_33matrix(tensor; leftdir=ITensors.In, kwargs...)
+"""
+    convert_itensor_33matrix(tensor; leftdir=ITensors.In, kwargs...)
+
+  Converts a 4-legged tensor (coming from an (infinite) MPO) with two site indices and a left and a right leg into a 3 x 3 matrix of ITensor.
+  We assume the normal form for MPO (before full compression) where the top left and bottom right corners are identity matrices.
+  The goal is to write the tensor in the form
+         1      M_12   M_13
+         M_21   M_22   M_23
+         M_31   M_32   1
+  such that we can then easily compress it. Note that for most of our tensors, the upper triangular part will be 0.
+  Input: tensor the four leg tensors
+  Output: the 3x3 matrix of tensors
+  Optional arguments: leftdir: the direction of the left_index of the matrix to select it properly. Default to Itensors.In
+                      leftindex: instead one can directly provide the left index. Default to nothing
+                      siteindex: instead of relying on the tag, one can provide the two site indices. Default to tag filtering.
+                      left_tags: if we want to change the tags of the left indices. Default to tags(leftindex).
+                      right_tags: if we want to change the tags of the right indices. Default to tags(rightindex).
+"""
+function convert_itensor_33matrix(tensor; kwargs...)
   @assert order(tensor) == 4
-  left_ind = get(kwargs, :leftindex, nothing)
+  leftind = get(kwargs, :leftindex, nothing)
+  leftdir = get(kwargs, :leftdir, ITensors.In)
   #Identify the different indices
-  sit = filterinds(inds(tensor); tags="Site")
+  sit = get(kwargs, :siteindex, filterinds(inds(tensor); tags="Site"))
   local_sit = dag(only(filterinds(sit; plev=0)))
   #A bit roundabout as filterinds does not accept dir
-  if isnothing(left_ind)
+  if isnothing(leftind)
     temp = uniqueinds(tensor, sit)
     if dir(temp[1]) == leftdir
-      left_ind = temp[1]
+      leftind = temp[1]
       right_ind = temp[2]
     else
-      left_ind = temp[2]
+      leftind = temp[2]
       right_ind = temp[1]
     end
   else
-    right_ind = only(uniqueinds(tensor, sit, left_ind))
+    right_ind = only(uniqueinds(tensor, sit, leftind))
   end
-  left_dim = dim(left_ind)
+  left_dim = dim(leftind)
   right_dim = dim(right_ind)
   #Build the local projectors.
-  left_tags = get(kwargs, :left_tags, tags(left_ind))
+  left_tags = get(kwargs, :left_tags, tags(leftind))
   top_left, middle_left, bottom_left = build_three_projectors_from_index(
-    left_ind; tags=left_tags
+    leftind; tags=left_tags
   )
   right_tags = get(kwargs, :righ_tags, tags(right_ind))
   top_right, middle_right, bottom_right = build_three_projectors_from_index(
@@ -188,12 +215,27 @@ function convert_itensor_33matrix(tensor; leftdir=ITensors.In, kwargs...)
   return matrix
 end
 
+"""
+    convert_itensor_3vector(tensor; leftdir=ITensors.In, kwargs...)
+
+  Converts a 3-legged tensor (the extremity of a MPO) with two site indices and one leg into a 3 Vector of ITensor.
+  We assume the normal form for MPO (before full compression) where the top left and bottom right corners are identity matrices in the bulk.
+
+  Input: tensor the three leg tensors
+  Output: the 3x1 or 1x3 vector of tensors
+  Optional arguments: leftdir: the direction of the left_index of the matrix to decide whether we get a 3x1 or 1x3 vector. Default to Itensors.In
+                      first: Alternative way to specify the shape, where it is the first tensor of the MPO. Default to false
+                      last:  Alternative way to specify the shape, where it is the last tensor of the MPO. Default to false
+                      siteindex: instead of relying on the tag, one can provide the two site indices.
+                      left_tags: if we want to change the tags of the left indices. Else default to tags(leftindex)
+                      right_tags: if we want to change the tags of the right indices. Else default to tags(rightindex)
+"""
 function convert_itensor_3vector(
   tensor; leftdir=ITensors.In, first=false, last=false, kwargs...
 )
   @assert order(tensor) == 3
   #Identify the different indices
-  sit = filterinds(inds(tensor); tags="Site")
+  sit = get(kwargs, :siteindex, filterinds(inds(tensor); tags="Site"))
   local_sit = dag(only(filterinds(sit; plev=0)))
   #A bit roundabout as filterinds does not accept dir
   old_ind = only(uniqueinds(tensor, sit))
@@ -211,24 +253,3 @@ function convert_itensor_3vector(
   end
   return vector
 end
-
-function apply_tensor(A::Array{ITensor,N}, B::ITensor...) where {N}
-  new_A = copy(A)
-  for x in 1:length(new_A)
-    new_A[x] = contract(new_A[x], B...)
-  end
-  return new_A
-end
-
-# import ITensors._add
-# function ITensors._add(
-#   A::T1, B::T2
-# ) where {T1<:ITensors.NDTensors.EmptyTensor,T2<:ITensors.NDTensors.EmptyTensor}
-#   ndims(A) != ndims(B) && throw(
-#     ITensors.DimensionMismatch("cannot add ITensors with different numbers of indices")
-#   )
-#   indices = inds(A)
-#   length(commoninds(indices, B)) != length(inds(B)) &&
-#     error("cannot add ITensors with different indices")
-#   return ITensor(promote_type(eltype(A), eltype(B)), indices)
-# end

src/models/models.jl

@@ -166,48 +166,46 @@ function InfiniteMPOMatrix(model::Model, s::CelledVector, translator::Function;
     sp = prime(s[x])
     sd = dag(s[x])
     left_inds = [
-      only(uniqueinds(mpos.data.data[x][j, 1], mpos.data.data[x][1, 1])) for
-      j in 2:(size(mpos.data.data[x], 1) - 1)
+      only(uniqueinds(mpos[x][j, 1], mpos[x][1, 1])) for j in 2:(size(mpos[x], 1) - 1)
     ]
     right_inds = [
-      only(uniqueinds(mpos.data.data[x][end, j], mpos.data.data[x][1, 1])) for
-      j in 2:(size(mpos.data.data[x], 2) - 1)
+      only(uniqueinds(mpos[x][end, j], mpos[x][1, 1])) for j in 2:(size(mpos[x], 2) - 1)
     ]
     if x == N
       new_right_inds = [
-        dag(only(uniqueinds(mpos.data.data[1][j, 1], mpos.data.data[1][1, 1]))) for
-        j in 2:(size(mpos.data.data[1], 1) - 1)
+        dag(only(uniqueinds(mpos[1][j, 1], mpos[1][1, 1]))) for
+        j in 2:(size(mpos[1], 1) - 1)
       ]
       for j in 1:length(new_right_inds)
         new_right_inds[j] = translatecell(translator, new_right_inds[j], 1)
       end
     else
       new_right_inds = [
-        dag(only(uniqueinds(mpos.data.data[x + 1][j, 1], mpos.data.data[x + 1][1, 1]))) for
-        j in 2:(size(mpos.data.data[x], 2) - 1)
+        dag(only(uniqueinds(mpos[x + 1][j, 1], mpos[x + 1][1, 1]))) for
+        j in 2:(size(mpos[x], 2) - 1)
       ]
     end
-    for j in 2:(size(mpos.data.data[x], 1) - 1)
-      for k in 2:(size(mpos.data.data[x], 2) - 1)
-        if isempty(mpos.data.data[x][j, k])
+    for j in 2:(size(mpos[x], 1) - 1)
+      for k in 2:(size(mpos[x], 2) - 1)
+        if isempty(mpos[x][j, k])
           mpos.data.data[x][j, k] = ITensor(left_inds[j - 1], sd, sp, new_right_inds[k - 1])
         else
           replaceinds!(mpos.data.data[x][j, k], right_inds[k - 1] => new_right_inds[k - 1])
         end
       end
     end
-    for j in [1, size(mpos.data.data[x], 1)]
-      for k in 2:(size(mpos.data.data[x], 2) - 1)
-        if isempty(mpos.data.data[x][j, k])
+    for j in [1, size(mpos[x], 1)]
+      for k in 2:(size(mpos[x], 2) - 1)
+        if isempty(mpos[x][j, k])
           mpos.data.data[x][j, k] = ITensor(sd, sp, new_right_inds[k - 1])
         else
           replaceinds!(mpos.data.data[x][j, k], right_inds[k - 1] => new_right_inds[k - 1])
         end
       end
     end
-    for j in 2:(size(mpos.data.data[x], 1) - 1)
-      for k in [1, size(mpos.data.data[x], 2)]
-        if isempty(mpos.data.data[x][j, k])
+    for j in 2:(size(mpos[x], 1) - 1)
+      for k in [1, size(mpos[x], 2)]
+        if isempty(mpos[x][j, k])
           mpos.data.data[x][j, k] = ITensor(left_inds[j - 1], sd, sp)
         end
       end

src/vumps_mpo.jl

@@ -1,9 +1,3 @@
-function ITensors.NDTensors.contraction_output(
-  A::NDTensors.EmptyTensor, B::NDTensors.DiagBlockSparseTensor, label
-)
-  return ITensor(eltype(B), label)
-end
-
 # Struct for use in linear system solver
 struct AOᴸ
   ψ::InfiniteCanonicalMPS

test/test_iMPOConversions.jl

@@ -86,12 +86,12 @@ function ITensors.expect(ψ::InfiniteCanonicalMPS, h::InfiniteMPOMatrix)
   Ncell = nsites(h)
   L, R = generate_edges(h)
   l = commoninds(ψ.AL[0], ψ.AL[1])
-  L = ITensorInfiniteMPS.apply_tensor(L, δ(l, dag(prime(l))))
+  L = map(a -> contract(a, δ(l, dag(prime(l)))), L)
   r = commoninds(ψ.AR[Ncell + 1], ψ.AR[Ncell])
-  R = ITensorInfiniteMPS.apply_tensor(R, δ(r, dag(prime(r))))
-  L = ITensorInfiniteMPS.apply_tensor(L, ψ.C[0], dag(prime(ψ.C[0])))
+  R = map(a -> contract(a, δ(r, dag(prime(r)))), R)
+  L = map(a -> contract(a, ψ.C[0], dag(prime(ψ.C[0]))), L)
   for j in 1:nsites(ψ)
-    temp = ITensorInfiniteMPS.apply_tensor(h[j], ψ.AR[j], dag(prime(ψ.AR[j])))
+    temp = map(a -> contract(a, ψ.AR[j], dag(prime(ψ.AR[j]))), h[j])
     L = L * temp
   end
   return (L * R)[1][1]#ITensorInfiniteMPS.scalar_product(L, R)[1]