Add prefactors to ProjTTNSum, optimize ProjOuterProdTTN

src/solvers/contract.jl

@@ -10,5 +10,5 @@ function contract_updater(
   updater_kwargs,
 )
   P = projected_operator![]
-  return contract_ket(P, ITensor(true)), (;)
+  return contract_ket(P, ITensor(one(Bool))), (;)
 end

src/treetensornetworks/projttns/abstractprojttn.jl

@@ -101,6 +101,9 @@ function Base.eltype(P::AbstractProjTTN)::Type
   return ElType
 end
 
+vertextype(::Type{<:AbstractProjTTN{V}}) where {V} = V
+vertextype(p::AbstractProjTTN) = typeof(p)
+
 function Base.size(P::AbstractProjTTN)::Tuple{Int,Int}
   d = 1
   for e in incident_edges(P)

src/treetensornetworks/projttns/projouterprodttn.jl

@@ -106,8 +106,8 @@ function contract_ket(P::ProjOuterProdTTN, v::ITensor)
   return v
 end
 
-# ToDo: debug this, not yet working
-# probably 
+# ToDo: verify conjugation etc. with complex AbstractTTN
 function contract(P::ProjOuterProdTTN, x::ITensor)
-  return conj(contract_ket(P, dag(x))) * contract_ket(P, ITensor(true))
+  ket = contract_ket(P, ITensor(one(Bool)))
+  return (dag(ket) * x) * ket
 end

src/treetensornetworks/projttns/projttnsum.jl

@@ -1,18 +1,25 @@
 """
 ProjTTNSum
 """
-struct ProjTTNSum{V} <: AbstractProjTTN{V}
-  terms::Vector{T} where {T<:AbstractProjTTN{V}}
-  function ProjTTNSum(terms::Vector{T}) where {V,T<:AbstractProjTTN{V}}
-    return new{V}(terms)
+struct ProjTTNSum{V,T<:AbstractProjTTN{V},Z<:Number} <: AbstractProjTTN{V}
+  terms::Vector{T}
+  factors::Vector{Z}
+  function ProjTTNSum(terms::Vector{<:AbstractProjTTN}, factors::Vector{<:Number})
+    return new{vertextype(eltype(terms)),eltype(terms),eltype(factors)}(terms, factors)
   end
 end
 
 terms(P::ProjTTNSum) = P.terms
+factors(P::ProjTTNSum) = P.factors
 
 copy(P::ProjTTNSum) = ProjTTNSum(copy.(terms(P)))
 
-ProjTTNSum(operators::Vector{<:AbstractTTN}) = ProjTTNSum(ProjTTN.(operators))
+function ProjTTNSum(operators::Vector{<:AbstractProjTTN})
+  return ProjTTNSum(operators, fill(1, length(operators)))
+end
+function ProjTTNSum(operators::Vector{<:AbstractTTN})
+  return ProjTTNSum(ProjTTN.(operators), fill(1, length(operators)))
+end
 
 on_edge(P::ProjTTNSum) = on_edge(terms(P)[1])
 
@@ -34,9 +41,15 @@ internal_edges(P::ProjTTNSum) = internal_edges(terms(P)[1])
 
 product(P::ProjTTNSum, v::ITensor) = noprime(contract(P, v))
 
-contract(P::ProjTTNSum, v::ITensor) = sum(p -> contract(p, v), terms(P))
+contract(P::ProjTTNSum, v::ITensor) =
+  mapreduce(+, zip(factors(P), terms(P))) do (f, p)
+    f * contract(p, v)
+  end
 
-contract_ket(P::ProjTTNSum, v::ITensor) = sum(p -> contract_ket(p, v), terms(P))
+contract_ket(P::ProjTTNSum, v::ITensor) =
+  mapreduce(+, zip(factors(P), terms(P))) do (f, p)
+    f * contract_ket(p, v)
+  end
 
 function Base.eltype(P::ProjTTNSum)
   return mapreduce(eltype, promote_type, terms(P))

test/test_treetensornetworks/test_solvers/test_contract.jl

@@ -28,17 +28,26 @@ using Test
   Hfit = ProjOuterProdTTN(psi', H)
   Hpsi_via_dmrg = dmrg(Hfit, psi; updater_kwargs=(; which_eigval=:LR,), nsweeps=1)
   @test abs(inner(Hpsi_via_dmrg, Hpsi / norm(Hpsi))) ≈ 1 atol = 1E-4
+  # Test whether the interface works for ProjTTNSum with factors
+  Hfit = ProjTTNSum([ProjOuterProdTTN(psi', H), ProjOuterProdTTN(psi', H)], [0.5, 0.5])
+  Hpsi_via_dmrg = dmrg(Hfit, psi; nsweeps=1, updater_kwargs=(; which_eigval=:LR,))
+  @test abs(inner(Hpsi_via_dmrg, Hpsi / norm(Hpsi))) ≈ 1 atol = 1E-4
 
   # Test basic usage for use with multiple ProjOuterProdTTN with default parameters
   # BLAS.axpy-like test
   os_id = OpSum()
   os_id += -1, "Id", 1, "Id", 2
   minus_identity = mpo(os_id, s)
   Hpsi = ITensorNetworks.sum_apply(
-    [(H, psi), (minus_identity, copy(psi))]; alg="fit", init=psi, nsweeps=1
+    [(H, psi), (minus_identity, psi)]; alg="fit", init=psi, nsweeps=1
   )
   @test inner(psi, Hpsi) ≈ (inner(psi', H, psi) - norm(psi)^2) atol = 1E-5
 
+  #Hpsi = ITensorNetworks.dmrg(
+  #  [(H, psi), (minus_identity, psi)]; alg="fit", init=psi, nsweeps=1
+  #)
+  #@test inner(psi, Hpsi) ≈ (inner(psi', H, psi) - norm(psi)^2) atol = 1E-5
+
   #
   # Change "top" indices of MPO to be a different set
   #
@@ -84,7 +93,7 @@ end
   os_id += -1, "Id", vertices(s)[1], "Id", vertices(s)[1]
   minus_identity = TTN(os_id, s)
   Hpsi = ITensorNetworks.sum_apply(
-    [(H, psi), (minus_identity, copy(psi))]; alg="fit", init=psi, nsweeps=1
+    [(H, psi), (minus_identity, psi)]; alg="fit", init=psi, nsweeps=1
   )
   @test inner(psi, Hpsi) ≈ (inner(psi', H, psi) - norm(psi)^2) atol = 1E-5