From 830fdc403abac76aae61538a461e243121d3e220 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Wed, 31 Jan 2024 18:54:52 -0500 Subject: [PATCH 01/27] Generalize ProjTTNSum to vector of concrete subtype of AbstractProjTTN instead of ProjTTN. --- src/treetensornetworks/projttns/projttnsum.jl | 21 ++++++++++++------- 1 file changed, 14 insertions(+), 7 deletions(-) diff --git a/src/treetensornetworks/projttns/projttnsum.jl b/src/treetensornetworks/projttns/projttnsum.jl index 7ce816f9..994912d1 100644 --- a/src/treetensornetworks/projttns/projttnsum.jl +++ b/src/treetensornetworks/projttns/projttnsum.jl @@ -1,18 +1,24 @@ """ ProjTTNSum """ -struct ProjTTNSum{V} - pm::Vector{ProjTTN{V}} - function ProjTTNSum(pm::Vector{ProjTTN{V}}) where {V} - return new{V}(pm) +struct ProjTTNSum{T<:AbstractProjTTN} + pm::Vector{T} where T + function ProjTTNSum(pm::Vector{T}) where {T} + return new{T}(pm) end end +#ToDo: define accessor functions copy(P::ProjTTNSum) = ProjTTNSum(copy.(P.pm)) +# The following constructors don't generalize well, maybe require to pass AbstractProjTTN instead of TTN and remove these? ProjTTNSum(ttnos::Vector{<:TTN}) = ProjTTNSum([ProjTTN(M) for M in ttnos]) +function ProjTTNSum(init_states::Vector{<:TTN},ttnos::Vector{<:TTN}) + return ProjTTNSum([ProjTTNApply(state,operator) for (state,operator) in zip(init_states,ttnos)]) +end -ProjTTNSum(Ms::TTN{V}...) where {V} = ProjTTNSum([Ms...]) +# This constructor can't differentiate between different concrete AbstractProjTTN, remove? +ProjTTNSum(Ms::TTN{V}...) where {V} = ProjTTNSum([Ms...]) on_edge(P::ProjTTNSum) = on_edge(P.pm[1]) @@ -52,8 +58,9 @@ end Base.size(P::ProjTTNSum) = size(P.pm[1]) +#ToDo remove parametrization? function position( - P::ProjTTNSum{V}, psi::TTN{V}, pos::Union{Vector{<:V},NamedEdge{V}} -) where {V} + P::ProjTTNSum, psi::TTN, pos +) return ProjTTNSum(map(M -> position(M, psi, pos), P.pm)) end From 4983fa45e8684c51685f30920961a6fca8fd68d5 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Wed, 31 Jan 2024 19:03:40 -0500 Subject: [PATCH 02/27] Format. --- src/treetensornetworks/projttns/projttnsum.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/treetensornetworks/projttns/projttnsum.jl b/src/treetensornetworks/projttns/projttnsum.jl index 994912d1..dc3d9fb0 100644 --- a/src/treetensornetworks/projttns/projttnsum.jl +++ b/src/treetensornetworks/projttns/projttnsum.jl @@ -2,7 +2,7 @@ ProjTTNSum """ struct ProjTTNSum{T<:AbstractProjTTN} - pm::Vector{T} where T + pm::Vector{T} where {T} function ProjTTNSum(pm::Vector{T}) where {T} return new{T}(pm) end @@ -13,12 +13,14 @@ copy(P::ProjTTNSum) = ProjTTNSum(copy.(P.pm)) # The following constructors don't generalize well, maybe require to pass AbstractProjTTN instead of TTN and remove these? ProjTTNSum(ttnos::Vector{<:TTN}) = ProjTTNSum([ProjTTN(M) for M in ttnos]) -function ProjTTNSum(init_states::Vector{<:TTN},ttnos::Vector{<:TTN}) - return ProjTTNSum([ProjTTNApply(state,operator) for (state,operator) in zip(init_states,ttnos)]) +function ProjTTNSum(init_states::Vector{<:TTN}, ttnos::Vector{<:TTN}) + return ProjTTNSum([ + ProjTTNApply(state, operator) for (state, operator) in zip(init_states, ttnos) + ]) end # This constructor can't differentiate between different concrete AbstractProjTTN, remove? -ProjTTNSum(Ms::TTN{V}...) where {V} = ProjTTNSum([Ms...]) +ProjTTNSum(Ms::TTN{V}...) where {V} = ProjTTNSum([Ms...]) on_edge(P::ProjTTNSum) = on_edge(P.pm[1]) @@ -59,8 +61,6 @@ end Base.size(P::ProjTTNSum) = size(P.pm[1]) #ToDo remove parametrization? -function position( - P::ProjTTNSum, psi::TTN, pos -) +function position(P::ProjTTNSum, psi::TTN, pos) return ProjTTNSum(map(M -> position(M, psi, pos), P.pm)) end From c86bf86782e3983a4c4c803c74edfad68de9f30e Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Thu, 1 Feb 2024 11:07:53 -0500 Subject: [PATCH 03/27] Implement missing methods for generalized ProjTTNSum and add test . --- src/solvers/contract.jl | 9 +-- .../projttns/abstractprojttn.jl | 8 +++ src/treetensornetworks/projttns/projttn.jl | 8 +++ .../projttns/projttn_apply.jl | 31 ++++++++++ src/treetensornetworks/projttns/projttnsum.jl | 8 +++ src/treetensornetworks/solvers/contract.jl | 58 ++++++++++++++----- .../test_solvers/test_contract.jl | 11 +++- 7 files changed, 111 insertions(+), 22 deletions(-) diff --git a/src/solvers/contract.jl b/src/solvers/contract.jl index ed55a729..dd41f844 100644 --- a/src/solvers/contract.jl +++ b/src/solvers/contract.jl @@ -9,11 +9,6 @@ function contract_updater( region_kwargs, updater_kwargs, ) - v = ITensor(true) - projected_operator = projected_operator![] - for j in sites(projected_operator) - v *= init_state(projected_operator)[j] - end - vp = contract(projected_operator, v) - return vp, (;) + P=projected_operator![] + return contract(P), (;) end diff --git a/src/treetensornetworks/projttns/abstractprojttn.jl b/src/treetensornetworks/projttns/abstractprojttn.jl index e3b0e140..f23b8dc2 100644 --- a/src/treetensornetworks/projttns/abstractprojttn.jl +++ b/src/treetensornetworks/projttns/abstractprojttn.jl @@ -88,6 +88,10 @@ function contract(P::AbstractProjTTN, v::ITensor)::ITensor return Hv end +function contract(P::AbstractProjTTN)::ITensor + error("Not implemented") +end + function product(P::AbstractProjTTN, v::ITensor)::ITensor Pv = contract(P, v) if order(Pv) != order(v) @@ -105,6 +109,10 @@ function product(P::AbstractProjTTN, v::ITensor)::ITensor return noprime(Pv) end +function product(P::AbstractProjTTN)::ITensor + error("Not implemented") +end + (P::AbstractProjTTN)(v::ITensor) = product(P, v) function Base.eltype(P::AbstractProjTTN)::Type diff --git a/src/treetensornetworks/projttns/projttn.jl b/src/treetensornetworks/projttns/projttn.jl index 725650b2..805ef960 100644 --- a/src/treetensornetworks/projttns/projttn.jl +++ b/src/treetensornetworks/projttns/projttn.jl @@ -68,3 +68,11 @@ function make_environment(P::ProjTTN, state::AbstractTTN, e::AbstractEdge) ) return P end + +function contract(P::ProjTTN)::ITensor + error("Not implemented") +end + +function product(P::ProjTTN)::ITensor + error("Not implemented") +end \ No newline at end of file diff --git a/src/treetensornetworks/projttns/projttn_apply.jl b/src/treetensornetworks/projttns/projttn_apply.jl index d4f691b7..d52855fe 100644 --- a/src/treetensornetworks/projttns/projttn_apply.jl +++ b/src/treetensornetworks/projttns/projttn_apply.jl @@ -71,3 +71,34 @@ function make_environment(P::ProjTTNApply, state::AbstractTTN, e::AbstractEdge) ) return P end + +function contract(P::ProjTTNApply)::ITensor + environments = ITensor[environment(P, edge) for edge in incident_edges(P)] + # manual heuristic for contraction order fixing: for each site in ProjTTN, apply up to + # two environments, then TTN tensor, then other environments + if on_edge(P) + itensor_map = environments + else + itensor_map = Union{ITensor,OneITensor}[] # TODO: will a Hamiltonian TTN tensor ever be a OneITensor? + for s in sites(P) + site_envs = filter(hascommoninds(operator(P)[s]), environments) + frst, scnd, rst = _separate_first_two(site_envs) + site_tensors = vcat(frst, scnd, operator(P)[s], rst) + append!(itensor_map, site_tensors) + end + end + Hv = ITensor(true) + for j in sites(P) + Hv *= init_state(P)[j] + end + # TODO: actually use optimal contraction sequence here + for it in itensor_map + Hv *= it + end + return Hv +end + +#ToDo: Is this a good idea? +function product(P::ProjTTNApply)::ITensor + return noprime(contract(P::ProjTTNApply)) +end \ No newline at end of file diff --git a/src/treetensornetworks/projttns/projttnsum.jl b/src/treetensornetworks/projttns/projttnsum.jl index dc3d9fb0..7b3242aa 100644 --- a/src/treetensornetworks/projttns/projttnsum.jl +++ b/src/treetensornetworks/projttns/projttnsum.jl @@ -48,6 +48,14 @@ function product(P::ProjTTNSum, v::ITensor)::ITensor return Pv end +function contract(P::ProjTTNSum)::ITensor + Pv = contract(P.pm[1]) + for n in 2:length(P.pm) + Pv += contract(P.pm[n]) + end + return Pv +end + function Base.eltype(P::ProjTTNSum) elT = eltype(P.pm[1]) for n in 2:length(P.pm) diff --git a/src/treetensornetworks/solvers/contract.jl b/src/treetensornetworks/solvers/contract.jl index 398138e7..4c7fb1fb 100644 --- a/src/treetensornetworks/solvers/contract.jl +++ b/src/treetensornetworks/solvers/contract.jl @@ -1,30 +1,45 @@ function contract( ::Algorithm"fit", - tn1::AbstractTTN, - tn2::AbstractTTN; - init=random_ttn(flatten_external_indsnetwork(tn1, tn2); link_space=trivial_space(tn1)), + tn1s::Vector{<:AbstractTTN}, + tn2s::Vector{<:AbstractTTN}; + #ToDo: this default probably doesn't work with QNs? + init=random_ttn(flatten_external_indsnetwork(first(tn1s), first(tn2s)); link_space=trivial_space(first(tn1s))), ##not generic yet nsweeps=1, nsites=2, # used to be default of call to default_sweep_regions updater_kwargs=(;), kwargs..., -) - n = nv(tn1) - n != nv(tn2) && throw( +) + ns = nv.(tn1s) + n = first(ns) + any(ns .!= nv.(tn2s)) && throw( DimensionMismatch("Number of sites operator ($n) and state ($(nv(tn2))) do not match") ) + any(ns .!= n) && throw( + DimensionMismatch("Number of sites in different operators ($n) do not match") + ) + # ToDo: Write test for single-vertex TTN, this implementation has not been tested. if n == 1 - v = only(vertices(tn2)) - return typeof(tn2)([tn1[v] * tn2[v]]) + res = 0 + for (tn1,tn2) in zip(tn1s,tn2s) + v = only(vertices(tn2)) + res += tn1[v] * tn2[v] + end + typeof([res]) + return end # check_hascommoninds(siteinds, tn1, tn2) # In case `tn1` and `tn2` have the same internal indices - tn1 = sim(linkinds, tn1) - - # In case `init` and `tn2` have the same internal indices - init = sim(linkinds, init) + PHs=ProjTTNApply{vertextype(first(tn1s))}[] + for (tn1,tn2) in zip(tn1s,tn2s) + tn1 = sim(linkinds, tn1) + # In case `init` and `tn2` have the same internal indices + init = sim(linkinds, init) + push!(PHs,ProjTTNApply(tn2, tn1)) + end + PH = isone(length(PHs)==1) ? only(PHs) : ProjTTNSum(PHs) # Fix site and link inds of init ## init = deepcopy(init) ## init = sim(linkinds, init) @@ -33,8 +48,6 @@ function contract( ## init[v], siteinds(init, v), uniqueinds(siteinds(tn1, v), siteinds(tn2, v)) ## ) ## end - - PH = ProjTTNApply(tn2, tn1) sweep_plan = default_sweep_regions(nsites, init; kwargs...) psi = alternating_update( contract_updater, PH, init; nsweeps, sweep_plan, updater_kwargs, kwargs... @@ -43,6 +56,12 @@ function contract( return psi end +function contract( + a::Algorithm"fit", tn1::AbstractTTN, tn2::AbstractTTN; kwargs...) + return contract(a, [tn1,], [tn2,]; kwargs...) +end + + """ Overload of `ITensors.contract`. """ @@ -50,6 +69,12 @@ function contract(tn1::AbstractTTN, tn2::AbstractTTN; alg="fit", kwargs...) return contract(Algorithm(alg), tn1, tn2; kwargs...) end +function contract(tn1s::Vector{<:AbstractTTN}, tn2s::Vector{<:AbstractTTN}; alg="fit", kwargs...) + return contract(Algorithm(alg), tn1s, tn2s; kwargs...) +end + + + """ Overload of `ITensors.apply`. """ @@ -57,3 +82,8 @@ function apply(tn1::AbstractTTN, tn2::AbstractTTN; kwargs...) tn12 = contract(tn1, tn2; kwargs...) return replaceprime(tn12, 1 => 0) end + +function apply(tn1s::Vector{<:AbstractTTN}, tn2s::Vector{<:AbstractTTN}; alg="fit", kwargs...) + tn12 = contract(tn1s, tn2s; alg="fit",kwargs...) + return replaceprime(tn12, 1 => 0) +end \ No newline at end of file diff --git a/test/test_treetensornetworks/test_solvers/test_contract.jl b/test/test_treetensornetworks/test_solvers/test_contract.jl index 810db6fb..3b6ca493 100644 --- a/test/test_treetensornetworks/test_solvers/test_contract.jl +++ b/test/test_treetensornetworks/test_solvers/test_contract.jl @@ -20,11 +20,19 @@ using Test os += "Sz", j, "Sz", j + 2 end H = mpo(os, s) - + # Test basic usage with default parameters Hpsi = apply(H, psi; alg="fit", init=psi') @test inner(psi, Hpsi) ≈ inner(psi', H, psi) atol = 1E-5 + # Test basic usage for use with multiple ProjTTNApply with default parameters + # BLAS.axpy-like test + os_id = OpSum() + os_id += -1, "Id", 1, "Id", 2 + minus_identity = mpo(os_id,s) + Hpsi = apply([H,minus_identity], [psi,copy(psi)]; alg="fit", init=psi') + @test inner(psi, Hpsi) ≈ (inner(psi', H, psi)-norm(psi)^2) atol = 1E-5 + # # Change "top" indices of MPO to be a different set # @@ -50,6 +58,7 @@ using Test @test inner(psit, Hpsi) ≈ inner(psit, H, psi) atol = 1E-4 end + @testset "Contract TTN" begin tooth_lengths = fill(2, 3) root_vertex = (3, 2) From 4a75a4dbd63c0921fe56f1a528287e34ca945a64 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Thu, 1 Feb 2024 11:27:39 -0500 Subject: [PATCH 04/27] Rename field and add its accessor for ProjTTNSum. --- src/imports.jl | 4 ++ src/treetensornetworks/projttns/projttnsum.jl | 60 +++++++++---------- 2 files changed, 31 insertions(+), 33 deletions(-) diff --git a/src/imports.jl b/src/imports.jl index 71406878..8d6a52a8 100644 --- a/src/imports.jl +++ b/src/imports.jl @@ -96,6 +96,10 @@ import ITensors: scalartype, #adding add + +import ITensors.LazyApply: + # extracting terms from a sum + terms #Algorithm Algorithm diff --git a/src/treetensornetworks/projttns/projttnsum.jl b/src/treetensornetworks/projttns/projttnsum.jl index 7b3242aa..bfae1f59 100644 --- a/src/treetensornetworks/projttns/projttnsum.jl +++ b/src/treetensornetworks/projttns/projttnsum.jl @@ -2,73 +2,67 @@ ProjTTNSum """ struct ProjTTNSum{T<:AbstractProjTTN} - pm::Vector{T} where {T} - function ProjTTNSum(pm::Vector{T}) where {T} - return new{T}(pm) + terms::Vector{T} + + function ProjTTNSum(terms::Vector{<:AbstractProjTTN}) + return new{eltype(terms)}(terms) end end -#ToDo: define accessor functions -copy(P::ProjTTNSum) = ProjTTNSum(copy.(P.pm)) +terms(P::ProjTTNSum) = P.terms -# The following constructors don't generalize well, maybe require to pass AbstractProjTTN instead of TTN and remove these? -ProjTTNSum(ttnos::Vector{<:TTN}) = ProjTTNSum([ProjTTN(M) for M in ttnos]) -function ProjTTNSum(init_states::Vector{<:TTN}, ttnos::Vector{<:TTN}) - return ProjTTNSum([ - ProjTTNApply(state, operator) for (state, operator) in zip(init_states, ttnos) - ]) -end +copy(P::ProjTTNSum) = ProjTTNSum(copy.(terms(P))) -# This constructor can't differentiate between different concrete AbstractProjTTN, remove? -ProjTTNSum(Ms::TTN{V}...) where {V} = ProjTTNSum([Ms...]) +#ToDo: Should we get rid of this one as well? +ProjTTNSum(ttnos::Vector{<:TTN}) = ProjTTNSum([ProjTTN(M) for M in ttnos]) -on_edge(P::ProjTTNSum) = on_edge(P.pm[1]) +on_edge(P::ProjTTNSum) = on_edge(terms(P)[1]) -nsite(P::ProjTTNSum) = nsite(P.pm[1]) +nsite(P::ProjTTNSum) = nsite(terms(P)[1]) function set_nsite(Ps::ProjTTNSum, nsite) - return ProjTTNSum(map(M -> set_nsite(M, nsite), Ps.pm)) + return ProjTTNSum(map(M -> set_nsite(M, nsite), Ps.terms)) end -underlying_graph(P::ProjTTNSum) = underlying_graph(P.pm[1]) +underlying_graph(P::ProjTTNSum) = underlying_graph(terms(P)[1]) -Base.length(P::ProjTTNSum) = length(P.pm[1]) +Base.length(P::ProjTTNSum) = length(terms(P)[1]) -sites(P::ProjTTNSum) = sites(P.pm[1]) +sites(P::ProjTTNSum) = sites(terms(P)[1]) -incident_edges(P::ProjTTNSum) = incident_edges(P.pm[1]) +incident_edges(P::ProjTTNSum) = incident_edges(terms(P)[1]) -internal_edges(P::ProjTTNSum) = internal_edges(P.pm[1]) +internal_edges(P::ProjTTNSum) = internal_edges(terms(P)[1]) function product(P::ProjTTNSum, v::ITensor)::ITensor - Pv = product(P.pm[1], v) - for n in 2:length(P.pm) - Pv += product(P.pm[n], v) + Pv = product(terms(P)[1], v) + for n in 2:length(terms(P)) + Pv += product(terms(P)[n], v) end return Pv end function contract(P::ProjTTNSum)::ITensor - Pv = contract(P.pm[1]) - for n in 2:length(P.pm) - Pv += contract(P.pm[n]) + Pv = contract(terms(P)[1]) + for n in 2:length(terms(P)) + Pv += contract(terms(P)[n]) end return Pv end function Base.eltype(P::ProjTTNSum) - elT = eltype(P.pm[1]) - for n in 2:length(P.pm) - elT = promote_type(elT, eltype(P.pm[n])) + elT = eltype(terms(P)[1]) + for n in 2:length(terms(P)) + elT = promote_type(elT, eltype(terms(P)[n])) end return elT end (P::ProjTTNSum)(v::ITensor) = product(P, v) -Base.size(P::ProjTTNSum) = size(P.pm[1]) +Base.size(P::ProjTTNSum) = size(terms(P)[1]) #ToDo remove parametrization? function position(P::ProjTTNSum, psi::TTN, pos) - return ProjTTNSum(map(M -> position(M, psi, pos), P.pm)) + return ProjTTNSum(map(M -> position(M, psi, pos), terms(P))) end From afa6982d2360e57027664b2a10a58015e8dd4ca1 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Thu, 1 Feb 2024 11:28:10 -0500 Subject: [PATCH 05/27] Add TTN test. --- .../test_treetensornetworks/test_solvers/test_contract.jl | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/test/test_treetensornetworks/test_solvers/test_contract.jl b/test/test_treetensornetworks/test_solvers/test_contract.jl index 3b6ca493..8c893530 100644 --- a/test/test_treetensornetworks/test_solvers/test_contract.jl +++ b/test/test_treetensornetworks/test_solvers/test_contract.jl @@ -73,6 +73,14 @@ end # Test basic usage with default parameters Hpsi = apply(H, psi; alg="fit") @test inner(psi, Hpsi) ≈ inner(psi', H, psi) atol = 1E-5 + + # Test basic usage for multiple ProjTTNApply with default parameters + # BLAS.axpy-like test + os_id = OpSum() + os_id += -1, "Id", vertices(s)[1], "Id", vertices(s)[1] + minus_identity = TTN(os_id,s) + Hpsi = apply([H,minus_identity], [psi,copy(psi)]; alg="fit", init=psi') + @test inner(psi, Hpsi) ≈ (inner(psi', H, psi)-norm(psi)^2) atol = 1E-5 # # Change "top" indices of TTN to be a different set From 788d788001380320fd151c51d27fce23800646ba Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Thu, 1 Feb 2024 11:35:22 -0500 Subject: [PATCH 06/27] Format. --- src/imports.jl | 2 +- src/solvers/contract.jl | 2 +- .../projttns/abstractprojttn.jl | 4 +- src/treetensornetworks/projttns/projttn.jl | 6 +-- .../projttns/projttn_apply.jl | 2 +- src/treetensornetworks/projttns/projttnsum.jl | 4 +- src/treetensornetworks/solvers/contract.jl | 44 ++++++++++--------- .../test_solvers/test_contract.jl | 17 ++++--- 8 files changed, 41 insertions(+), 40 deletions(-) diff --git a/src/imports.jl b/src/imports.jl index 8d6a52a8..e7071450 100644 --- a/src/imports.jl +++ b/src/imports.jl @@ -96,7 +96,7 @@ import ITensors: scalartype, #adding add - + import ITensors.LazyApply: # extracting terms from a sum terms diff --git a/src/solvers/contract.jl b/src/solvers/contract.jl index dd41f844..065ed73e 100644 --- a/src/solvers/contract.jl +++ b/src/solvers/contract.jl @@ -9,6 +9,6 @@ function contract_updater( region_kwargs, updater_kwargs, ) - P=projected_operator![] + P = projected_operator![] return contract(P), (;) end diff --git a/src/treetensornetworks/projttns/abstractprojttn.jl b/src/treetensornetworks/projttns/abstractprojttn.jl index f23b8dc2..7978d180 100644 --- a/src/treetensornetworks/projttns/abstractprojttn.jl +++ b/src/treetensornetworks/projttns/abstractprojttn.jl @@ -89,7 +89,7 @@ function contract(P::AbstractProjTTN, v::ITensor)::ITensor end function contract(P::AbstractProjTTN)::ITensor - error("Not implemented") + return error("Not implemented") end function product(P::AbstractProjTTN, v::ITensor)::ITensor @@ -110,7 +110,7 @@ function product(P::AbstractProjTTN, v::ITensor)::ITensor end function product(P::AbstractProjTTN)::ITensor - error("Not implemented") + return error("Not implemented") end (P::AbstractProjTTN)(v::ITensor) = product(P, v) diff --git a/src/treetensornetworks/projttns/projttn.jl b/src/treetensornetworks/projttns/projttn.jl index 805ef960..549c6f3c 100644 --- a/src/treetensornetworks/projttns/projttn.jl +++ b/src/treetensornetworks/projttns/projttn.jl @@ -70,9 +70,9 @@ function make_environment(P::ProjTTN, state::AbstractTTN, e::AbstractEdge) end function contract(P::ProjTTN)::ITensor - error("Not implemented") + return error("Not implemented") end function product(P::ProjTTN)::ITensor - error("Not implemented") -end \ No newline at end of file + return error("Not implemented") +end diff --git a/src/treetensornetworks/projttns/projttn_apply.jl b/src/treetensornetworks/projttns/projttn_apply.jl index d52855fe..7029b510 100644 --- a/src/treetensornetworks/projttns/projttn_apply.jl +++ b/src/treetensornetworks/projttns/projttn_apply.jl @@ -101,4 +101,4 @@ end #ToDo: Is this a good idea? function product(P::ProjTTNApply)::ITensor return noprime(contract(P::ProjTTNApply)) -end \ No newline at end of file +end diff --git a/src/treetensornetworks/projttns/projttnsum.jl b/src/treetensornetworks/projttns/projttnsum.jl index bfae1f59..6cda6244 100644 --- a/src/treetensornetworks/projttns/projttnsum.jl +++ b/src/treetensornetworks/projttns/projttnsum.jl @@ -3,8 +3,8 @@ ProjTTNSum """ struct ProjTTNSum{T<:AbstractProjTTN} terms::Vector{T} - - function ProjTTNSum(terms::Vector{<:AbstractProjTTN}) + + function ProjTTNSum(terms::Vector{<:AbstractProjTTN}) return new{eltype(terms)}(terms) end end diff --git a/src/treetensornetworks/solvers/contract.jl b/src/treetensornetworks/solvers/contract.jl index 4c7fb1fb..a4f10b04 100644 --- a/src/treetensornetworks/solvers/contract.jl +++ b/src/treetensornetworks/solvers/contract.jl @@ -3,43 +3,45 @@ function contract( tn1s::Vector{<:AbstractTTN}, tn2s::Vector{<:AbstractTTN}; #ToDo: this default probably doesn't work with QNs? - init=random_ttn(flatten_external_indsnetwork(first(tn1s), first(tn2s)); link_space=trivial_space(first(tn1s))), ##not generic yet + init=random_ttn( + flatten_external_indsnetwork(first(tn1s), first(tn2s)); + link_space=trivial_space(first(tn1s)), + ), ##not generic yet nsweeps=1, nsites=2, # used to be default of call to default_sweep_regions updater_kwargs=(;), kwargs..., -) +) ns = nv.(tn1s) n = first(ns) any(ns .!= nv.(tn2s)) && throw( DimensionMismatch("Number of sites operator ($n) and state ($(nv(tn2))) do not match") ) - any(ns .!= n) && throw( - DimensionMismatch("Number of sites in different operators ($n) do not match") - ) + any(ns .!= n) && + throw(DimensionMismatch("Number of sites in different operators ($n) do not match")) # ToDo: Write test for single-vertex TTN, this implementation has not been tested. if n == 1 res = 0 - for (tn1,tn2) in zip(tn1s,tn2s) + for (tn1, tn2) in zip(tn1s, tn2s) v = only(vertices(tn2)) res += tn1[v] * tn2[v] end typeof([res]) - return + return nothing end # check_hascommoninds(siteinds, tn1, tn2) # In case `tn1` and `tn2` have the same internal indices - PHs=ProjTTNApply{vertextype(first(tn1s))}[] - for (tn1,tn2) in zip(tn1s,tn2s) + PHs = ProjTTNApply{vertextype(first(tn1s))}[] + for (tn1, tn2) in zip(tn1s, tn2s) tn1 = sim(linkinds, tn1) # In case `init` and `tn2` have the same internal indices init = sim(linkinds, init) - push!(PHs,ProjTTNApply(tn2, tn1)) + push!(PHs, ProjTTNApply(tn2, tn1)) end - PH = isone(length(PHs)==1) ? only(PHs) : ProjTTNSum(PHs) + PH = isone(length(PHs) == 1) ? only(PHs) : ProjTTNSum(PHs) # Fix site and link inds of init ## init = deepcopy(init) ## init = sim(linkinds, init) @@ -56,12 +58,10 @@ function contract( return psi end -function contract( - a::Algorithm"fit", tn1::AbstractTTN, tn2::AbstractTTN; kwargs...) - return contract(a, [tn1,], [tn2,]; kwargs...) +function contract(a::Algorithm"fit", tn1::AbstractTTN, tn2::AbstractTTN; kwargs...) + return contract(a, [tn1], [tn2]; kwargs...) end - """ Overload of `ITensors.contract`. """ @@ -69,12 +69,12 @@ function contract(tn1::AbstractTTN, tn2::AbstractTTN; alg="fit", kwargs...) return contract(Algorithm(alg), tn1, tn2; kwargs...) end -function contract(tn1s::Vector{<:AbstractTTN}, tn2s::Vector{<:AbstractTTN}; alg="fit", kwargs...) +function contract( + tn1s::Vector{<:AbstractTTN}, tn2s::Vector{<:AbstractTTN}; alg="fit", kwargs... +) return contract(Algorithm(alg), tn1s, tn2s; kwargs...) end - - """ Overload of `ITensors.apply`. """ @@ -83,7 +83,9 @@ function apply(tn1::AbstractTTN, tn2::AbstractTTN; kwargs...) return replaceprime(tn12, 1 => 0) end -function apply(tn1s::Vector{<:AbstractTTN}, tn2s::Vector{<:AbstractTTN}; alg="fit", kwargs...) - tn12 = contract(tn1s, tn2s; alg="fit",kwargs...) +function apply( + tn1s::Vector{<:AbstractTTN}, tn2s::Vector{<:AbstractTTN}; alg="fit", kwargs... +) + tn12 = contract(tn1s, tn2s; alg="fit", kwargs...) return replaceprime(tn12, 1 => 0) -end \ No newline at end of file +end diff --git a/test/test_treetensornetworks/test_solvers/test_contract.jl b/test/test_treetensornetworks/test_solvers/test_contract.jl index 8c893530..0d37e401 100644 --- a/test/test_treetensornetworks/test_solvers/test_contract.jl +++ b/test/test_treetensornetworks/test_solvers/test_contract.jl @@ -20,7 +20,7 @@ using Test os += "Sz", j, "Sz", j + 2 end H = mpo(os, s) - + # Test basic usage with default parameters Hpsi = apply(H, psi; alg="fit", init=psi') @test inner(psi, Hpsi) ≈ inner(psi', H, psi) atol = 1E-5 @@ -29,9 +29,9 @@ using Test # BLAS.axpy-like test os_id = OpSum() os_id += -1, "Id", 1, "Id", 2 - minus_identity = mpo(os_id,s) - Hpsi = apply([H,minus_identity], [psi,copy(psi)]; alg="fit", init=psi') - @test inner(psi, Hpsi) ≈ (inner(psi', H, psi)-norm(psi)^2) atol = 1E-5 + minus_identity = mpo(os_id, s) + Hpsi = apply([H, minus_identity], [psi, copy(psi)]; alg="fit", init=psi') + @test inner(psi, Hpsi) ≈ (inner(psi', H, psi) - norm(psi)^2) atol = 1E-5 # # Change "top" indices of MPO to be a different set @@ -58,7 +58,6 @@ using Test @test inner(psit, Hpsi) ≈ inner(psit, H, psi) atol = 1E-4 end - @testset "Contract TTN" begin tooth_lengths = fill(2, 3) root_vertex = (3, 2) @@ -73,14 +72,14 @@ end # Test basic usage with default parameters Hpsi = apply(H, psi; alg="fit") @test inner(psi, Hpsi) ≈ inner(psi', H, psi) atol = 1E-5 - + # Test basic usage for multiple ProjTTNApply with default parameters # BLAS.axpy-like test os_id = OpSum() os_id += -1, "Id", vertices(s)[1], "Id", vertices(s)[1] - minus_identity = TTN(os_id,s) - Hpsi = apply([H,minus_identity], [psi,copy(psi)]; alg="fit", init=psi') - @test inner(psi, Hpsi) ≈ (inner(psi', H, psi)-norm(psi)^2) atol = 1E-5 + minus_identity = TTN(os_id, s) + Hpsi = apply([H, minus_identity], [psi, copy(psi)]; alg="fit", init=psi') + @test inner(psi, Hpsi) ≈ (inner(psi', H, psi) - norm(psi)^2) atol = 1E-5 # # Change "top" indices of TTN to be a different set From 616ccaeea85d0980361a9c9aad601199175c19d4 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Thu, 1 Feb 2024 11:42:19 -0500 Subject: [PATCH 07/27] Deleted some comments and unnecessary type restrictions. --- src/treetensornetworks/projttns/abstractprojttn.jl | 4 ++-- src/treetensornetworks/projttns/projttnsum.jl | 1 - src/treetensornetworks/solvers/contract.jl | 2 +- 3 files changed, 3 insertions(+), 4 deletions(-) diff --git a/src/treetensornetworks/projttns/abstractprojttn.jl b/src/treetensornetworks/projttns/abstractprojttn.jl index 7978d180..45808685 100644 --- a/src/treetensornetworks/projttns/abstractprojttn.jl +++ b/src/treetensornetworks/projttns/abstractprojttn.jl @@ -88,7 +88,7 @@ function contract(P::AbstractProjTTN, v::ITensor)::ITensor return Hv end -function contract(P::AbstractProjTTN)::ITensor +function contract(P::AbstractProjTTN) return error("Not implemented") end @@ -109,7 +109,7 @@ function product(P::AbstractProjTTN, v::ITensor)::ITensor return noprime(Pv) end -function product(P::AbstractProjTTN)::ITensor +function product(P::AbstractProjTTN) return error("Not implemented") end diff --git a/src/treetensornetworks/projttns/projttnsum.jl b/src/treetensornetworks/projttns/projttnsum.jl index 6cda6244..4906e34a 100644 --- a/src/treetensornetworks/projttns/projttnsum.jl +++ b/src/treetensornetworks/projttns/projttnsum.jl @@ -13,7 +13,6 @@ terms(P::ProjTTNSum) = P.terms copy(P::ProjTTNSum) = ProjTTNSum(copy.(terms(P))) -#ToDo: Should we get rid of this one as well? ProjTTNSum(ttnos::Vector{<:TTN}) = ProjTTNSum([ProjTTN(M) for M in ttnos]) on_edge(P::ProjTTNSum) = on_edge(terms(P)[1]) diff --git a/src/treetensornetworks/solvers/contract.jl b/src/treetensornetworks/solvers/contract.jl index a4f10b04..1a0c184a 100644 --- a/src/treetensornetworks/solvers/contract.jl +++ b/src/treetensornetworks/solvers/contract.jl @@ -6,7 +6,7 @@ function contract( init=random_ttn( flatten_external_indsnetwork(first(tn1s), first(tn2s)); link_space=trivial_space(first(tn1s)), - ), ##not generic yet + ), nsweeps=1, nsites=2, # used to be default of call to default_sweep_regions updater_kwargs=(;), From 6fdc0e636024739c6382f6451d03a3df2806eaf4 Mon Sep 17 00:00:00 2001 From: b-kloss Date: Thu, 1 Feb 2024 11:53:32 -0500 Subject: [PATCH 08/27] Replace access by field by accessor function. Co-authored-by: Matt Fishman --- src/treetensornetworks/projttns/projttnsum.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/treetensornetworks/projttns/projttnsum.jl b/src/treetensornetworks/projttns/projttnsum.jl index 4906e34a..fef61f3e 100644 --- a/src/treetensornetworks/projttns/projttnsum.jl +++ b/src/treetensornetworks/projttns/projttnsum.jl @@ -20,7 +20,7 @@ on_edge(P::ProjTTNSum) = on_edge(terms(P)[1]) nsite(P::ProjTTNSum) = nsite(terms(P)[1]) function set_nsite(Ps::ProjTTNSum, nsite) - return ProjTTNSum(map(M -> set_nsite(M, nsite), Ps.terms)) + return ProjTTNSum(map(p -> set_nsite(p, nsite), terms(Ps))) end underlying_graph(P::ProjTTNSum) = underlying_graph(terms(P)[1]) From 120204952f7d947c4a8d3431b53dcdf0201ff90b Mon Sep 17 00:00:00 2001 From: b-kloss Date: Thu, 1 Feb 2024 12:13:55 -0500 Subject: [PATCH 09/27] Rename ttnos to operators in constructor. Co-authored-by: Matt Fishman --- src/treetensornetworks/projttns/projttnsum.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/treetensornetworks/projttns/projttnsum.jl b/src/treetensornetworks/projttns/projttnsum.jl index fef61f3e..87496252 100644 --- a/src/treetensornetworks/projttns/projttnsum.jl +++ b/src/treetensornetworks/projttns/projttnsum.jl @@ -13,7 +13,7 @@ terms(P::ProjTTNSum) = P.terms copy(P::ProjTTNSum) = ProjTTNSum(copy.(terms(P))) -ProjTTNSum(ttnos::Vector{<:TTN}) = ProjTTNSum([ProjTTN(M) for M in ttnos]) +ProjTTNSum(operators::Vector{<:AbstractTTN}) = ProjTTNSum(ProjTTN.(operators)) on_edge(P::ProjTTNSum) = on_edge(terms(P)[1]) From aa464bf3d30fb4aac19a97734d833358d8fd2c07 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Thu, 1 Feb 2024 13:27:58 -0500 Subject: [PATCH 10/27] refactor contract --- src/solvers/contract.jl | 2 +- .../projttns/abstractprojttn.jl | 34 ++++--------------- src/treetensornetworks/projttns/projttn.jl | 22 ++++++++---- .../projttns/projttn_apply.jl | 24 +++++-------- src/treetensornetworks/projttns/projttnsum.jl | 6 ++-- 5 files changed, 36 insertions(+), 52 deletions(-) diff --git a/src/solvers/contract.jl b/src/solvers/contract.jl index 065ed73e..2b8d6209 100644 --- a/src/solvers/contract.jl +++ b/src/solvers/contract.jl @@ -10,5 +10,5 @@ function contract_updater( updater_kwargs, ) P = projected_operator![] - return contract(P), (;) + return contract(P, ITensor(true)), (;) end diff --git a/src/treetensornetworks/projttns/abstractprojttn.jl b/src/treetensornetworks/projttns/abstractprojttn.jl index 45808685..af3314ff 100644 --- a/src/treetensornetworks/projttns/abstractprojttn.jl +++ b/src/treetensornetworks/projttns/abstractprojttn.jl @@ -65,31 +65,15 @@ function _separate_first_two(V::Vector) return frst, scnd, rst end +projected_operator_tensors(P::AbstractProjTTN) = error("Not implemented.") + +#contract(P::AbstractProjTTN, v::ITensor)::ITensor = foldl(*,projected_operator_tensors(P);init=v) function contract(P::AbstractProjTTN, v::ITensor)::ITensor - environments = ITensor[environment(P, edge) for edge in incident_edges(P)] - # manual heuristic for contraction order fixing: for each site in ProjTTN, apply up to - # two environments, then TTN tensor, then other environments - if on_edge(P) - itensor_map = environments - else - itensor_map = Union{ITensor,OneITensor}[] # TODO: will a Hamiltonian TTN tensor ever be a OneITensor? - for s in sites(P) - site_envs = filter(hascommoninds(operator(P)[s]), environments) - frst, scnd, rst = _separate_first_two(site_envs) - site_tensors = vcat(frst, scnd, operator(P)[s], rst) - append!(itensor_map, site_tensors) - end + itensor_map = projected_operator_tensors(P) + for t in itensor_map + v *= t end - # TODO: actually use optimal contraction sequence here - Hv = v - for it in itensor_map - Hv *= it - end - return Hv -end - -function contract(P::AbstractProjTTN) - return error("Not implemented") + return v end function product(P::AbstractProjTTN, v::ITensor)::ITensor @@ -109,10 +93,6 @@ function product(P::AbstractProjTTN, v::ITensor)::ITensor return noprime(Pv) end -function product(P::AbstractProjTTN) - return error("Not implemented") -end - (P::AbstractProjTTN)(v::ITensor) = product(P, v) function Base.eltype(P::AbstractProjTTN)::Type diff --git a/src/treetensornetworks/projttns/projttn.jl b/src/treetensornetworks/projttns/projttn.jl index 549c6f3c..4d9636ab 100644 --- a/src/treetensornetworks/projttns/projttn.jl +++ b/src/treetensornetworks/projttns/projttn.jl @@ -69,10 +69,20 @@ function make_environment(P::ProjTTN, state::AbstractTTN, e::AbstractEdge) return P end -function contract(P::ProjTTN)::ITensor - return error("Not implemented") -end - -function product(P::ProjTTN)::ITensor - return error("Not implemented") +function projected_operator_tensors(P::ProjTTN) + environments = ITensor[environment(P, edge) for edge in incident_edges(P)] + # manual heuristic for contraction order fixing: for each site in ProjTTN, apply up to + # two environments, then TTN tensor, then other environments + if on_edge(P) + itensor_map = environments + else + itensor_map = Union{ITensor,OneITensor}[] # TODO: will a Hamiltonian TTN tensor ever be a OneITensor? + for s in sites(P) + site_envs = filter(hascommoninds(operator(P)[s]), environments) + frst, scnd, rst = _separate_first_two(site_envs) + site_tensors = vcat(frst, scnd, operator(P)[s], rst) + append!(itensor_map, site_tensors) + end + end + return itensor_map end diff --git a/src/treetensornetworks/projttns/projttn_apply.jl b/src/treetensornetworks/projttns/projttn_apply.jl index 7029b510..2072db9f 100644 --- a/src/treetensornetworks/projttns/projttn_apply.jl +++ b/src/treetensornetworks/projttns/projttn_apply.jl @@ -72,14 +72,17 @@ function make_environment(P::ProjTTNApply, state::AbstractTTN, e::AbstractEdge) return P end -function contract(P::ProjTTNApply)::ITensor +function projected_operator_tensors(P::ProjTTNApply) environments = ITensor[environment(P, edge) for edge in incident_edges(P)] # manual heuristic for contraction order fixing: for each site in ProjTTN, apply up to # two environments, then TTN tensor, then other environments + itensor_map = Union{ITensor,OneITensor}[] # TODO: will a Hamiltonian TTN tensor ever be a OneITensor? + for j in sites(P) + push!(itensor_map, init_state(P)[j]) + end if on_edge(P) - itensor_map = environments + append!(itensor_map, environments) else - itensor_map = Union{ITensor,OneITensor}[] # TODO: will a Hamiltonian TTN tensor ever be a OneITensor? for s in sites(P) site_envs = filter(hascommoninds(operator(P)[s]), environments) frst, scnd, rst = _separate_first_two(site_envs) @@ -87,18 +90,9 @@ function contract(P::ProjTTNApply)::ITensor append!(itensor_map, site_tensors) end end - Hv = ITensor(true) - for j in sites(P) - Hv *= init_state(P)[j] - end - # TODO: actually use optimal contraction sequence here - for it in itensor_map - Hv *= it - end - return Hv + return itensor_map end -#ToDo: Is this a good idea? -function product(P::ProjTTNApply)::ITensor - return noprime(contract(P::ProjTTNApply)) +function product(P::ProjTTNApply, v::ITensor)::ITensor + return noprime(contract(P, v)) end diff --git a/src/treetensornetworks/projttns/projttnsum.jl b/src/treetensornetworks/projttns/projttnsum.jl index 4906e34a..b724f5b0 100644 --- a/src/treetensornetworks/projttns/projttnsum.jl +++ b/src/treetensornetworks/projttns/projttnsum.jl @@ -41,10 +41,10 @@ function product(P::ProjTTNSum, v::ITensor)::ITensor return Pv end -function contract(P::ProjTTNSum)::ITensor - Pv = contract(terms(P)[1]) +function contract(P::ProjTTNSum, v::ITensor)::ITensor + Pv = contract(terms(P)[1], v) for n in 2:length(terms(P)) - Pv += contract(terms(P)[n]) + Pv += contract(terms(P)[n], v) end return Pv end From b6a8f32fdf924d69446e1ea4b2e3979eeaa9f217 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Thu, 1 Feb 2024 16:15:55 -0500 Subject: [PATCH 11/27] Fix priming behaviours of apply and remove default for init. --- src/treetensornetworks/projttns/projttnsum.jl | 25 ++------- src/treetensornetworks/solvers/contract.jl | 53 ++++++++++++------- .../test_solvers/test_contract.jl | 19 ++++--- 3 files changed, 49 insertions(+), 48 deletions(-) diff --git a/src/treetensornetworks/projttns/projttnsum.jl b/src/treetensornetworks/projttns/projttnsum.jl index 8a75f97e..60c01477 100644 --- a/src/treetensornetworks/projttns/projttnsum.jl +++ b/src/treetensornetworks/projttns/projttnsum.jl @@ -33,35 +33,18 @@ incident_edges(P::ProjTTNSum) = incident_edges(terms(P)[1]) internal_edges(P::ProjTTNSum) = internal_edges(terms(P)[1]) -function product(P::ProjTTNSum, v::ITensor)::ITensor - Pv = product(terms(P)[1], v) - for n in 2:length(terms(P)) - Pv += product(terms(P)[n], v) - end - return Pv -end +product(P::ProjTTNSum, v::ITensor) = noprime(contract(P, v)) -function contract(P::ProjTTNSum, v::ITensor)::ITensor - Pv = contract(terms(P)[1], v) - for n in 2:length(terms(P)) - Pv += contract(terms(P)[n], v) - end - return Pv -end +contract(P::ProjTTNSum, v::ITensor)::ITensor = sum(p -> contract(p, v), terms(P)) function Base.eltype(P::ProjTTNSum) - elT = eltype(terms(P)[1]) - for n in 2:length(terms(P)) - elT = promote_type(elT, eltype(terms(P)[n])) - end - return elT + return mapreduce(eltype, promote_type, terms(P)) end (P::ProjTTNSum)(v::ITensor) = product(P, v) Base.size(P::ProjTTNSum) = size(terms(P)[1]) -#ToDo remove parametrization? -function position(P::ProjTTNSum, psi::TTN, pos) +function position(P::ProjTTNSum, psi::AbstractTTN, pos) return ProjTTNSum(map(M -> position(M, psi, pos), terms(P))) end diff --git a/src/treetensornetworks/solvers/contract.jl b/src/treetensornetworks/solvers/contract.jl index 1a0c184a..4bf945b2 100644 --- a/src/treetensornetworks/solvers/contract.jl +++ b/src/treetensornetworks/solvers/contract.jl @@ -1,17 +1,18 @@ function contract( ::Algorithm"fit", - tn1s::Vector{<:AbstractTTN}, - tn2s::Vector{<:AbstractTTN}; - #ToDo: this default probably doesn't work with QNs? - init=random_ttn( - flatten_external_indsnetwork(first(tn1s), first(tn2s)); - link_space=trivial_space(first(tn1s)), - ), + tns::Vector{<:Tuple{<:AbstractTTN,<:AbstractTTN}}; + init, + #init=random_ttn( + # flatten_external_indsnetwork(first(tn1s), first(tn2s)); + # link_space=trivial_space(first(tn1s)), + #), nsweeps=1, nsites=2, # used to be default of call to default_sweep_regions updater_kwargs=(;), kwargs..., ) + tn1s = first.(tns) + tn2s = last.(tns) ns = nv.(tn1s) n = first(ns) any(ns .!= nv.(tn2s)) && throw( @@ -59,7 +60,7 @@ function contract( end function contract(a::Algorithm"fit", tn1::AbstractTTN, tn2::AbstractTTN; kwargs...) - return contract(a, [tn1], [tn2]; kwargs...) + return contract(a, [(tn1, tn2)]; kwargs...) end """ @@ -69,23 +70,37 @@ function contract(tn1::AbstractTTN, tn2::AbstractTTN; alg="fit", kwargs...) return contract(Algorithm(alg), tn1, tn2; kwargs...) end -function contract( - tn1s::Vector{<:AbstractTTN}, tn2s::Vector{<:AbstractTTN}; alg="fit", kwargs... -) - return contract(Algorithm(alg), tn1s, tn2s; kwargs...) -end - """ Overload of `ITensors.apply`. """ -function apply(tn1::AbstractTTN, tn2::AbstractTTN; kwargs...) - tn12 = contract(tn1, tn2; kwargs...) +function apply(tn1::AbstractTTN, tn2::AbstractTTN; init, kwargs...) + #plin=plev(first(externalinds(init))) + #outputindsnet= + #plout=plev(outputindsnet[first(vertices(outputindsnet))]) + plev_inc = plev_diff(flatten_external_indsnetwork(tn1, tn2), external_indsnetwork(init)) + init = prime(init, plev_inc) + tn12 = contract(tn1, tn2; init, kwargs...) return replaceprime(tn12, 1 => 0) end -function apply( - tn1s::Vector{<:AbstractTTN}, tn2s::Vector{<:AbstractTTN}; alg="fit", kwargs... +function sum_apply( + tns::Vector{<:Tuple{<:AbstractTTN,<:AbstractTTN}}; alg="fit", init, kwargs... ) - tn12 = contract(tn1s, tn2s; alg="fit", kwargs...) + #plin=plev(first(externalinds(init))) + #outputindsnet = flatten_external_indsnetwork(first(tns), first(tn2s)) + #plout=plev(outputindsnet[first(vertices(outputindsnet))]) + plev_inc = plev_diff( + flatten_external_indsnetwork(first(first(tns)), last(first(tns))), + external_indsnetwork(init), + ) + init = prime(init, plev_inc) + alg != "fit" && error("sum_apply not implemented for other algorithms than fit.") + tn12 = contract(Algorithm(alg), tns; init, kwargs...) return replaceprime(tn12, 1 => 0) end + +function plev_diff(a::IndsNetwork, b::IndsNetwork) + pla = plev(only(a[first(vertices(a))])) + plb = plev(only(b[first(vertices(b))])) + return pla - plb +end diff --git a/test/test_treetensornetworks/test_solvers/test_contract.jl b/test/test_treetensornetworks/test_solvers/test_contract.jl index 0d37e401..96b7f338 100644 --- a/test/test_treetensornetworks/test_solvers/test_contract.jl +++ b/test/test_treetensornetworks/test_solvers/test_contract.jl @@ -22,7 +22,7 @@ using Test H = mpo(os, s) # Test basic usage with default parameters - Hpsi = apply(H, psi; alg="fit", init=psi') + Hpsi = apply(H, psi; alg="fit", init=psi) @test inner(psi, Hpsi) ≈ inner(psi', H, psi) atol = 1E-5 # Test basic usage for use with multiple ProjTTNApply with default parameters @@ -30,7 +30,9 @@ using Test os_id = OpSum() os_id += -1, "Id", 1, "Id", 2 minus_identity = mpo(os_id, s) - Hpsi = apply([H, minus_identity], [psi, copy(psi)]; alg="fit", init=psi') + Hpsi = ITensorNetworks.sum_apply( + [(H, psi), (minus_identity, copy(psi))]; alg="fit", init=psi + ) @test inner(psi, Hpsi) ≈ (inner(psi', H, psi) - norm(psi)^2) atol = 1E-5 # @@ -38,15 +40,14 @@ using Test # t = siteinds("S=1/2", N) psit = deepcopy(psi) + for j in 1:N H[j] *= delta(s[j]', t[j]) psit[j] *= delta(s[j], t[j]) end - # Test with nsweeps=3 - Hpsi = apply(H, psi; alg="fit", nsweeps=3) + Hpsi = apply(H, psi; alg="fit", init=psit, nsweeps=3) @test inner(psit, Hpsi) ≈ inner(psit, H, psi) atol = 1E-5 - # Test with less good initial guess MPS not equal to psi psi_guess = truncate(psit; maxdim=2) Hpsi = apply(H, psi; alg="fit", nsweeps=4, init=psi_guess) @@ -70,7 +71,7 @@ end H = TTN(os, s) # Test basic usage with default parameters - Hpsi = apply(H, psi; alg="fit") + Hpsi = apply(H, psi; alg="fit", init=psi) @test inner(psi, Hpsi) ≈ inner(psi', H, psi) atol = 1E-5 # Test basic usage for multiple ProjTTNApply with default parameters @@ -78,7 +79,9 @@ end os_id = OpSum() os_id += -1, "Id", vertices(s)[1], "Id", vertices(s)[1] minus_identity = TTN(os_id, s) - Hpsi = apply([H, minus_identity], [psi, copy(psi)]; alg="fit", init=psi') + Hpsi = ITensorNetworks.sum_apply( + [(H, psi), (minus_identity, copy(psi))]; alg="fit", init=psi + ) @test inner(psi, Hpsi) ≈ (inner(psi', H, psi) - norm(psi)^2) atol = 1E-5 # @@ -90,7 +93,7 @@ end H = replaceinds(H, prime(s; links=[]) => t) # Test with nsweeps=2 - Hpsi = apply(H, psi; alg="fit", nsweeps=2) + Hpsi = apply(H, psi; alg="fit", init=psit, nsweeps=2) @test inner(psit, Hpsi) ≈ inner(psit, H, psi) atol = 1E-5 # Test with less good initial guess MPS not equal to psi From 80ccddc2aa31463da70c022ce48cca7f1a1b327e Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Thu, 1 Feb 2024 16:24:07 -0500 Subject: [PATCH 12/27] Remove prime-level logic from apply, and use contract when primelevel is not increased by operator. --- src/treetensornetworks/solvers/contract.jl | 21 ++----------------- .../test_solvers/test_contract.jl | 12 +++++------ 2 files changed, 8 insertions(+), 25 deletions(-) diff --git a/src/treetensornetworks/solvers/contract.jl b/src/treetensornetworks/solvers/contract.jl index 4bf945b2..8ce97aec 100644 --- a/src/treetensornetworks/solvers/contract.jl +++ b/src/treetensornetworks/solvers/contract.jl @@ -74,11 +74,7 @@ end Overload of `ITensors.apply`. """ function apply(tn1::AbstractTTN, tn2::AbstractTTN; init, kwargs...) - #plin=plev(first(externalinds(init))) - #outputindsnet= - #plout=plev(outputindsnet[first(vertices(outputindsnet))]) - plev_inc = plev_diff(flatten_external_indsnetwork(tn1, tn2), external_indsnetwork(init)) - init = prime(init, plev_inc) + init = init' tn12 = contract(tn1, tn2; init, kwargs...) return replaceprime(tn12, 1 => 0) end @@ -86,21 +82,8 @@ end function sum_apply( tns::Vector{<:Tuple{<:AbstractTTN,<:AbstractTTN}}; alg="fit", init, kwargs... ) - #plin=plev(first(externalinds(init))) - #outputindsnet = flatten_external_indsnetwork(first(tns), first(tn2s)) - #plout=plev(outputindsnet[first(vertices(outputindsnet))]) - plev_inc = plev_diff( - flatten_external_indsnetwork(first(first(tns)), last(first(tns))), - external_indsnetwork(init), - ) - init = prime(init, plev_inc) + init = init' alg != "fit" && error("sum_apply not implemented for other algorithms than fit.") tn12 = contract(Algorithm(alg), tns; init, kwargs...) return replaceprime(tn12, 1 => 0) end - -function plev_diff(a::IndsNetwork, b::IndsNetwork) - pla = plev(only(a[first(vertices(a))])) - plb = plev(only(b[first(vertices(b))])) - return pla - plb -end diff --git a/test/test_treetensornetworks/test_solvers/test_contract.jl b/test/test_treetensornetworks/test_solvers/test_contract.jl index 96b7f338..85982b69 100644 --- a/test/test_treetensornetworks/test_solvers/test_contract.jl +++ b/test/test_treetensornetworks/test_solvers/test_contract.jl @@ -46,16 +46,16 @@ using Test psit[j] *= delta(s[j], t[j]) end # Test with nsweeps=3 - Hpsi = apply(H, psi; alg="fit", init=psit, nsweeps=3) + Hpsi = contract(H, psi; alg="fit", init=psit, nsweeps=3) @test inner(psit, Hpsi) ≈ inner(psit, H, psi) atol = 1E-5 # Test with less good initial guess MPS not equal to psi psi_guess = truncate(psit; maxdim=2) - Hpsi = apply(H, psi; alg="fit", nsweeps=4, init=psi_guess) + Hpsi = contract(H, psi; alg="fit", nsweeps=4, init=psi_guess) @test inner(psit, Hpsi) ≈ inner(psit, H, psi) atol = 1E-5 # Test with nsite=1 Hpsi_guess = random_mps(t; internal_inds_space=32) - Hpsi = apply(H, psi; alg="fit", init=Hpsi_guess, nsites=1, nsweeps=4) + Hpsi = contract(H, psi; alg="fit", init=Hpsi_guess, nsites=1, nsweeps=4) @test inner(psit, Hpsi) ≈ inner(psit, H, psi) atol = 1E-4 end @@ -93,17 +93,17 @@ end H = replaceinds(H, prime(s; links=[]) => t) # Test with nsweeps=2 - Hpsi = apply(H, psi; alg="fit", init=psit, nsweeps=2) + Hpsi = contract(H, psi; alg="fit", init=psit, nsweeps=2) @test inner(psit, Hpsi) ≈ inner(psit, H, psi) atol = 1E-5 # Test with less good initial guess MPS not equal to psi Hpsi_guess = truncate(psit; maxdim=2) - Hpsi = apply(H, psi; alg="fit", nsweeps=4, init=Hpsi_guess) + Hpsi = contract(H, psi; alg="fit", nsweeps=4, init=Hpsi_guess) @test inner(psit, Hpsi) ≈ inner(psit, H, psi) atol = 1E-5 # Test with nsite=1 Hpsi_guess = random_ttn(t; link_space=4) - Hpsi = apply(H, psi; alg="fit", nsites=1, nsweeps=4, init=Hpsi_guess) + Hpsi = contract(H, psi; alg="fit", nsites=1, nsweeps=4, init=Hpsi_guess) @test inner(psit, Hpsi) ≈ inner(psit, H, psi) atol = 1E-4 end From bd69886cfb96674b4fdc724e7bd72cc13629ee60 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Thu, 1 Feb 2024 16:30:41 -0500 Subject: [PATCH 13/27] Add the primelevel logic back in as check in apply. --- src/treetensornetworks/solvers/contract.jl | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/src/treetensornetworks/solvers/contract.jl b/src/treetensornetworks/solvers/contract.jl index 8ce97aec..3799eda7 100644 --- a/src/treetensornetworks/solvers/contract.jl +++ b/src/treetensornetworks/solvers/contract.jl @@ -74,6 +74,11 @@ end Overload of `ITensors.apply`. """ function apply(tn1::AbstractTTN, tn2::AbstractTTN; init, kwargs...) + if !isone(plev_diff(flatten_external_indsnetwork(tn1, tn2), external_indsnetwork(init))) + error( + "Initial guess `init` needs to primelevel one less than the contraction tn1 and tn2." + ) + end init = init' tn12 = contract(tn1, tn2; init, kwargs...) return replaceprime(tn12, 1 => 0) @@ -82,8 +87,25 @@ end function sum_apply( tns::Vector{<:Tuple{<:AbstractTTN,<:AbstractTTN}}; alg="fit", init, kwargs... ) + if !isone( + plev_diff( + flatten_external_indsnetwork(first(first(tns)), last(first(tns))), + external_indsnetwork(init), + ), + ) + error( + "Initial guess `init` needs to primelevel one less than the contraction tn1 and tn2." + ) + end + init = init' alg != "fit" && error("sum_apply not implemented for other algorithms than fit.") tn12 = contract(Algorithm(alg), tns; init, kwargs...) return replaceprime(tn12, 1 => 0) end + +function plev_diff(a::IndsNetwork, b::IndsNetwork) + pla = plev(only(a[first(vertices(a))])) + plb = plev(only(b[first(vertices(b))])) + return pla - plb +end From 0cec308a18a3b6600e741bead5e859a145a885c9 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Thu, 1 Feb 2024 19:55:17 -0500 Subject: [PATCH 14/27] Rename ProjTTNApply to ProjOuterProdTTN and redefine its functionality, making it formally behave similarly to ProjTTN, with the old functionality accessible via contract_ket. --- src/exports.jl | 2 +- .../projttns/projttn_apply.jl | 54 ++++++++++++------- src/treetensornetworks/projttns/projttnsum.jl | 11 ++-- src/treetensornetworks/solvers/contract.jl | 20 +++---- .../test_solvers/test_contract.jl | 14 ++--- 5 files changed, 55 insertions(+), 46 deletions(-) diff --git a/src/exports.jl b/src/exports.jl index 97194466..e76c3ea9 100644 --- a/src/exports.jl +++ b/src/exports.jl @@ -76,7 +76,7 @@ export AbstractITensorNetwork, random_ttn, ProjTTN, ProjTTNSum, - ProjTTNApply, + ProjOuterProdTTN, set_nsite, position, finite_state_machine, diff --git a/src/treetensornetworks/projttns/projttn_apply.jl b/src/treetensornetworks/projttns/projttn_apply.jl index 2072db9f..a32bc5fc 100644 --- a/src/treetensornetworks/projttns/projttn_apply.jl +++ b/src/treetensornetworks/projttns/projttn_apply.jl @@ -1,45 +1,47 @@ -struct ProjTTNApply{V} <: AbstractProjTTN{V} +struct ProjOuterProdTTN{V} <: AbstractProjTTN{V} pos::Union{Vector{<:V},NamedEdge{V}} init_state::TTN{V} operator::TTN{V} environments::Dictionary{NamedEdge{V},ITensor} end -environments(p::ProjTTNApply) = p.environments -operator(p::ProjTTNApply) = p.operator -underlying_graph(p::ProjTTNApply) = underlying_graph(operator(p)) -pos(p::ProjTTNApply) = p.pos -init_state(p::ProjTTNApply) = p.init_state +environments(p::ProjOuterProdTTN) = p.environments +operator(p::ProjOuterProdTTN) = p.operator +underlying_graph(p::ProjOuterProdTTN) = underlying_graph(operator(p)) +pos(p::ProjOuterProdTTN) = p.pos +init_state(p::ProjOuterProdTTN) = p.init_state -function ProjTTNApply(init_state::AbstractTTN, operator::AbstractTTN) - return ProjTTNApply( +function ProjOuterProdTTN(init_state::AbstractTTN, operator::AbstractTTN) + return ProjOuterProdTTN( vertextype(operator)[], init_state, operator, Dictionary{edgetype(operator),ITensor}() ) end -function copy(P::ProjTTNApply) - return ProjTTNApply(P.pos, copy(init_state(P)), copy(operator(P)), copy(environments(P))) +function copy(P::ProjOuterProdTTN) + return ProjOuterProdTTN( + P.pos, copy(init_state(P)), copy(operator(P)), copy(environments(P)) + ) end -function set_nsite(P::ProjTTNApply, nsite) +function set_nsite(P::ProjOuterProdTTN, nsite) return P end -function shift_position(P::ProjTTNApply, pos) - return ProjTTNApply(pos, init_state(P), operator(P), environments(P)) +function shift_position(P::ProjOuterProdTTN, pos) + return ProjOuterProdTTN(pos, init_state(P), operator(P), environments(P)) end -function set_environments(p::ProjTTNApply, environments) - return ProjTTNApply(pos(p), init_state(p), operator(p), environments) +function set_environments(p::ProjOuterProdTTN, environments) + return ProjOuterProdTTN(pos(p), init_state(p), operator(p), environments) end -set_environment(p::ProjTTNApply, edge, env) = set_environment!(copy(p), edge, env) -function set_environment!(p::ProjTTNApply, edge, env) +set_environment(p::ProjOuterProdTTN, edge, env) = set_environment!(copy(p), edge, env) +function set_environment!(p::ProjOuterProdTTN, edge, env) set!(environments(p), edge, env) return p end -function make_environment(P::ProjTTNApply, state::AbstractTTN, e::AbstractEdge) +function make_environment(P::ProjOuterProdTTN, state::AbstractTTN, e::AbstractEdge) # invalidate environment for opposite edge direction if necessary reverse(e) ∈ incident_edges(P) || (P = invalidate_environment(P, reverse(e))) # do nothing if valid environment already present @@ -72,7 +74,7 @@ function make_environment(P::ProjTTNApply, state::AbstractTTN, e::AbstractEdge) return P end -function projected_operator_tensors(P::ProjTTNApply) +function projected_operator_tensors(P::ProjOuterProdTTN) environments = ITensor[environment(P, edge) for edge in incident_edges(P)] # manual heuristic for contraction order fixing: for each site in ProjTTN, apply up to # two environments, then TTN tensor, then other environments @@ -93,6 +95,18 @@ function projected_operator_tensors(P::ProjTTNApply) return itensor_map end -function product(P::ProjTTNApply, v::ITensor)::ITensor +function contract_ket(P::ProjOuterProdTTN, v::ITensor) + itensor_map = projected_operator_tensors(P) + for t in itensor_map + v *= t + end + return v +end + +function contract(P::ProjOuterProdTTN, x::ITensor) + return conj(contract_ket(P, x)) * contract_ket(P, ITensor(true)) +end + +function product(P::ProjOuterProdTTN, v::ITensor)::ITensor return noprime(contract(P, v)) end diff --git a/src/treetensornetworks/projttns/projttnsum.jl b/src/treetensornetworks/projttns/projttnsum.jl index 60c01477..03502521 100644 --- a/src/treetensornetworks/projttns/projttnsum.jl +++ b/src/treetensornetworks/projttns/projttnsum.jl @@ -1,11 +1,10 @@ """ ProjTTNSum """ -struct ProjTTNSum{T<:AbstractProjTTN} - terms::Vector{T} - - function ProjTTNSum(terms::Vector{<:AbstractProjTTN}) - return new{eltype(terms)}(terms) +struct ProjTTNSum{V} <: AbstractProjTTN{V} + terms::Vector{AbstractProjTTN{V}} + function ProjTTNSum(terms::Vector{<:AbstractProjTTN{V}}) where {V} + return new{V}(terms) end end @@ -37,6 +36,8 @@ product(P::ProjTTNSum, v::ITensor) = noprime(contract(P, v)) contract(P::ProjTTNSum, v::ITensor)::ITensor = sum(p -> contract(p, v), terms(P)) +contract_ket(P::ProjTTNSum, v::ITensor)::ITensor = sum(p -> contract_ket(p, v), terms(P)) + function Base.eltype(P::ProjTTNSum) return mapreduce(eltype, promote_type, terms(P)) end diff --git a/src/treetensornetworks/solvers/contract.jl b/src/treetensornetworks/solvers/contract.jl index 3799eda7..90e8c40a 100644 --- a/src/treetensornetworks/solvers/contract.jl +++ b/src/treetensornetworks/solvers/contract.jl @@ -1,12 +1,8 @@ -function contract( +function sum_contract( ::Algorithm"fit", tns::Vector{<:Tuple{<:AbstractTTN,<:AbstractTTN}}; init, - #init=random_ttn( - # flatten_external_indsnetwork(first(tn1s), first(tn2s)); - # link_space=trivial_space(first(tn1s)), - #), - nsweeps=1, + nsweeps, nsites=2, # used to be default of call to default_sweep_regions updater_kwargs=(;), kwargs..., @@ -27,20 +23,19 @@ function contract( v = only(vertices(tn2)) res += tn1[v] * tn2[v] end - typeof([res]) - return nothing + return typeof(tn2)([res]) end # check_hascommoninds(siteinds, tn1, tn2) # In case `tn1` and `tn2` have the same internal indices - PHs = ProjTTNApply{vertextype(first(tn1s))}[] + PHs = ProjOuterProdTTN{vertextype(first(tn1s))}[] for (tn1, tn2) in zip(tn1s, tn2s) tn1 = sim(linkinds, tn1) # In case `init` and `tn2` have the same internal indices init = sim(linkinds, init) - push!(PHs, ProjTTNApply(tn2, tn1)) + push!(PHs, ProjOuterProdTTN(tn2, tn1)) end PH = isone(length(PHs) == 1) ? only(PHs) : ProjTTNSum(PHs) # Fix site and link inds of init @@ -60,7 +55,7 @@ function contract( end function contract(a::Algorithm"fit", tn1::AbstractTTN, tn2::AbstractTTN; kwargs...) - return contract(a, [(tn1, tn2)]; kwargs...) + return sum_contract(a, [(tn1, tn2)]; kwargs...) end """ @@ -99,8 +94,7 @@ function sum_apply( end init = init' - alg != "fit" && error("sum_apply not implemented for other algorithms than fit.") - tn12 = contract(Algorithm(alg), tns; init, kwargs...) + tn12 = sum_contract(Algorithm(alg), tns; init, kwargs...) return replaceprime(tn12, 1 => 0) end diff --git a/test/test_treetensornetworks/test_solvers/test_contract.jl b/test/test_treetensornetworks/test_solvers/test_contract.jl index 85982b69..46eadf11 100644 --- a/test/test_treetensornetworks/test_solvers/test_contract.jl +++ b/test/test_treetensornetworks/test_solvers/test_contract.jl @@ -22,16 +22,16 @@ using Test H = mpo(os, s) # Test basic usage with default parameters - Hpsi = apply(H, psi; alg="fit", init=psi) + Hpsi = apply(H, psi; alg="fit", init=psi, nsweeps=1) @test inner(psi, Hpsi) ≈ inner(psi', H, psi) atol = 1E-5 - # Test basic usage for use with multiple ProjTTNApply with default parameters + # 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 + [(H, psi), (minus_identity, copy(psi))]; alg="fit", init=psi, nsweeps=1 ) @test inner(psi, Hpsi) ≈ (inner(psi', H, psi) - norm(psi)^2) atol = 1E-5 @@ -71,16 +71,16 @@ end H = TTN(os, s) # Test basic usage with default parameters - Hpsi = apply(H, psi; alg="fit", init=psi) + Hpsi = apply(H, psi; alg="fit", init=psi, nsweeps=1) @test inner(psi, Hpsi) ≈ inner(psi', H, psi) atol = 1E-5 - # Test basic usage for multiple ProjTTNApply with default parameters + # Test basic usage for multiple ProjOuterProdTTN with default parameters # BLAS.axpy-like test os_id = OpSum() 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 + [(H, psi), (minus_identity, copy(psi))]; alg="fit", init=psi, nsweeps=1 ) @test inner(psi, Hpsi) ≈ (inner(psi', H, psi) - norm(psi)^2) atol = 1E-5 @@ -121,7 +121,7 @@ end t2 = TreeTensorNetwork([M2[v] for v in eachindex(M2)]) # Test with good initial guess - @test contract(t1, t2; alg="fit", init=t12_ref) ≈ t12_ref rtol = 1e-7 + @test contract(t1, t2; alg="fit", init=t12_ref, nsweeps=1) ≈ t12_ref rtol = 1e-7 end nothing From ac96ce53900a374fa92201f022f0ba13f92fc9b7 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Thu, 1 Feb 2024 19:59:17 -0500 Subject: [PATCH 15/27] Rename init_state to internal_state for ProjOuterProdTTN. --- src/solvers/contract.jl | 2 +- .../projttns/projttn_apply.jl | 20 +++++++++---------- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/src/solvers/contract.jl b/src/solvers/contract.jl index 2b8d6209..6a659f48 100644 --- a/src/solvers/contract.jl +++ b/src/solvers/contract.jl @@ -10,5 +10,5 @@ function contract_updater( updater_kwargs, ) P = projected_operator![] - return contract(P, ITensor(true)), (;) + return contract_ket(P, ITensor(true)), (;) end diff --git a/src/treetensornetworks/projttns/projttn_apply.jl b/src/treetensornetworks/projttns/projttn_apply.jl index a32bc5fc..a52f54e9 100644 --- a/src/treetensornetworks/projttns/projttn_apply.jl +++ b/src/treetensornetworks/projttns/projttn_apply.jl @@ -1,6 +1,6 @@ struct ProjOuterProdTTN{V} <: AbstractProjTTN{V} pos::Union{Vector{<:V},NamedEdge{V}} - init_state::TTN{V} + internal_state::TTN{V} operator::TTN{V} environments::Dictionary{NamedEdge{V},ITensor} end @@ -9,17 +9,17 @@ environments(p::ProjOuterProdTTN) = p.environments operator(p::ProjOuterProdTTN) = p.operator underlying_graph(p::ProjOuterProdTTN) = underlying_graph(operator(p)) pos(p::ProjOuterProdTTN) = p.pos -init_state(p::ProjOuterProdTTN) = p.init_state +internal_state(p::ProjOuterProdTTN) = p.internal_state -function ProjOuterProdTTN(init_state::AbstractTTN, operator::AbstractTTN) +function ProjOuterProdTTN(internal_state::AbstractTTN, operator::AbstractTTN) return ProjOuterProdTTN( - vertextype(operator)[], init_state, operator, Dictionary{edgetype(operator),ITensor}() + vertextype(operator)[], internal_state, operator, Dictionary{edgetype(operator),ITensor}() ) end function copy(P::ProjOuterProdTTN) return ProjOuterProdTTN( - P.pos, copy(init_state(P)), copy(operator(P)), copy(environments(P)) + P.pos, copy(internal_state(P)), copy(operator(P)), copy(environments(P)) ) end @@ -28,11 +28,11 @@ function set_nsite(P::ProjOuterProdTTN, nsite) end function shift_position(P::ProjOuterProdTTN, pos) - return ProjOuterProdTTN(pos, init_state(P), operator(P), environments(P)) + return ProjOuterProdTTN(pos, internal_state(P), operator(P), environments(P)) end function set_environments(p::ProjOuterProdTTN, environments) - return ProjOuterProdTTN(pos(p), init_state(p), operator(p), environments) + return ProjOuterProdTTN(pos(p), internal_state(p), operator(p), environments) end set_environment(p::ProjOuterProdTTN, edge, env) = set_environment!(copy(p), edge, env) @@ -48,7 +48,7 @@ function make_environment(P::ProjOuterProdTTN, state::AbstractTTN, e::AbstractEd if !haskey(environments(P), e) if is_leaf(underlying_graph(P), src(e)) # leaves are easy - env = init_state(P)[src(e)] * operator(P)[src(e)] * dag(state[src(e)]) + env = internal_state(P)[src(e)] * operator(P)[src(e)] * dag(state[src(e)]) else # construct by contracting neighbors neighbor_envs = ITensor[] @@ -60,7 +60,7 @@ function make_environment(P::ProjOuterProdTTN, state::AbstractTTN, e::AbstractEd # other environments frst, scnd, rst = _separate_first_two(neighbor_envs) itensor_map = vcat( - init_state(P)[src(e)], frst, scnd, operator(P)[src(e)], dag(state[src(e)]), rst + internal_state(P)[src(e)], frst, scnd, operator(P)[src(e)], dag(state[src(e)]), rst ) # no prime here in comparison to the same routine for Projttn # TODO: actually use optimal contraction sequence here env = reduce(*, itensor_map) @@ -80,7 +80,7 @@ function projected_operator_tensors(P::ProjOuterProdTTN) # two environments, then TTN tensor, then other environments itensor_map = Union{ITensor,OneITensor}[] # TODO: will a Hamiltonian TTN tensor ever be a OneITensor? for j in sites(P) - push!(itensor_map, init_state(P)[j]) + push!(itensor_map, internal_state(P)[j]) end if on_edge(P) append!(itensor_map, environments) From 1d5f722f054b3f354c28edeb50c8b6dbe7106f4f Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Thu, 1 Feb 2024 20:00:43 -0500 Subject: [PATCH 16/27] Rename file. --- .../projttns/{projttn_apply.jl => projouterprodttn.jl} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename src/treetensornetworks/projttns/{projttn_apply.jl => projouterprodttn.jl} (100%) diff --git a/src/treetensornetworks/projttns/projttn_apply.jl b/src/treetensornetworks/projttns/projouterprodttn.jl similarity index 100% rename from src/treetensornetworks/projttns/projttn_apply.jl rename to src/treetensornetworks/projttns/projouterprodttn.jl From 572932afd064089f6116576c7383498e81413b9a Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Thu, 1 Feb 2024 20:10:07 -0500 Subject: [PATCH 17/27] Remove remaining type parametrization on output. --- src/treetensornetworks/projttns/abstractprojttn.jl | 2 +- src/treetensornetworks/projttns/projouterprodttn.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/treetensornetworks/projttns/abstractprojttn.jl b/src/treetensornetworks/projttns/abstractprojttn.jl index af3314ff..9ead5e6c 100644 --- a/src/treetensornetworks/projttns/abstractprojttn.jl +++ b/src/treetensornetworks/projttns/abstractprojttn.jl @@ -76,7 +76,7 @@ function contract(P::AbstractProjTTN, v::ITensor)::ITensor return v end -function product(P::AbstractProjTTN, v::ITensor)::ITensor +function product(P::AbstractProjTTN, v::ITensor) Pv = contract(P, v) if order(Pv) != order(v) error( diff --git a/src/treetensornetworks/projttns/projouterprodttn.jl b/src/treetensornetworks/projttns/projouterprodttn.jl index a52f54e9..6e4e20b9 100644 --- a/src/treetensornetworks/projttns/projouterprodttn.jl +++ b/src/treetensornetworks/projttns/projouterprodttn.jl @@ -107,6 +107,6 @@ function contract(P::ProjOuterProdTTN, x::ITensor) return conj(contract_ket(P, x)) * contract_ket(P, ITensor(true)) end -function product(P::ProjOuterProdTTN, v::ITensor)::ITensor +function product(P::ProjOuterProdTTN, v::ITensor) return noprime(contract(P, v)) end From e43b01611bdbfc31797255912a1956de2a41834c Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Thu, 1 Feb 2024 20:12:25 -0500 Subject: [PATCH 18/27] Remove more return type parametrizations. --- src/treetensornetworks/projttns/abstractprojttn.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/treetensornetworks/projttns/abstractprojttn.jl b/src/treetensornetworks/projttns/abstractprojttn.jl index 9ead5e6c..68aaabc1 100644 --- a/src/treetensornetworks/projttns/abstractprojttn.jl +++ b/src/treetensornetworks/projttns/abstractprojttn.jl @@ -67,8 +67,7 @@ end projected_operator_tensors(P::AbstractProjTTN) = error("Not implemented.") -#contract(P::AbstractProjTTN, v::ITensor)::ITensor = foldl(*,projected_operator_tensors(P);init=v) -function contract(P::AbstractProjTTN, v::ITensor)::ITensor +function contract(P::AbstractProjTTN, v::ITensor) itensor_map = projected_operator_tensors(P) for t in itensor_map v *= t From a9d75d3ce3acac97cfe71825c2542bcc611cff3a Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Thu, 1 Feb 2024 20:15:33 -0500 Subject: [PATCH 19/27] Format and removed another return type parametrization. --- src/treetensornetworks/projttns/projouterprodttn.jl | 5 ++++- src/treetensornetworks/projttns/projttnsum.jl | 4 ++-- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/src/treetensornetworks/projttns/projouterprodttn.jl b/src/treetensornetworks/projttns/projouterprodttn.jl index 6e4e20b9..278a66cf 100644 --- a/src/treetensornetworks/projttns/projouterprodttn.jl +++ b/src/treetensornetworks/projttns/projouterprodttn.jl @@ -13,7 +13,10 @@ internal_state(p::ProjOuterProdTTN) = p.internal_state function ProjOuterProdTTN(internal_state::AbstractTTN, operator::AbstractTTN) return ProjOuterProdTTN( - vertextype(operator)[], internal_state, operator, Dictionary{edgetype(operator),ITensor}() + vertextype(operator)[], + internal_state, + operator, + Dictionary{edgetype(operator),ITensor}(), ) end diff --git a/src/treetensornetworks/projttns/projttnsum.jl b/src/treetensornetworks/projttns/projttnsum.jl index 03502521..59c1a5ca 100644 --- a/src/treetensornetworks/projttns/projttnsum.jl +++ b/src/treetensornetworks/projttns/projttnsum.jl @@ -34,9 +34,9 @@ internal_edges(P::ProjTTNSum) = internal_edges(terms(P)[1]) product(P::ProjTTNSum, v::ITensor) = noprime(contract(P, v)) -contract(P::ProjTTNSum, v::ITensor)::ITensor = sum(p -> contract(p, v), terms(P)) +contract(P::ProjTTNSum, v::ITensor) = sum(p -> contract(p, v), terms(P)) -contract_ket(P::ProjTTNSum, v::ITensor)::ITensor = sum(p -> contract_ket(p, v), terms(P)) +contract_ket(P::ProjTTNSum, v::ITensor) = sum(p -> contract_ket(p, v), terms(P)) function Base.eltype(P::ProjTTNSum) return mapreduce(eltype, promote_type, terms(P)) From 7acb0d0ec07cc2dfb1fc6e713193dc185df82306 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Thu, 1 Feb 2024 20:24:38 -0500 Subject: [PATCH 20/27] Change filename in includes. --- src/ITensorNetworks.jl | 2 +- src/treetensornetworks/projttns/abstractprojttn.jl | 6 +----- 2 files changed, 2 insertions(+), 6 deletions(-) diff --git a/src/ITensorNetworks.jl b/src/ITensorNetworks.jl index 24fea311..d4c993a2 100644 --- a/src/ITensorNetworks.jl +++ b/src/ITensorNetworks.jl @@ -116,7 +116,7 @@ include(joinpath("treetensornetworks", "opsum_to_ttn.jl")) include(joinpath("treetensornetworks", "projttns", "abstractprojttn.jl")) include(joinpath("treetensornetworks", "projttns", "projttn.jl")) include(joinpath("treetensornetworks", "projttns", "projttnsum.jl")) -include(joinpath("treetensornetworks", "projttns", "projttn_apply.jl")) +include(joinpath("treetensornetworks", "projttns", "projouterprodttn.jl")) include(joinpath("treetensornetworks", "solvers", "solver_utils.jl")) include(joinpath("treetensornetworks", "solvers", "update_step.jl")) include(joinpath("treetensornetworks", "solvers", "alternating_update.jl")) diff --git a/src/treetensornetworks/projttns/abstractprojttn.jl b/src/treetensornetworks/projttns/abstractprojttn.jl index 68aaabc1..2739d45e 100644 --- a/src/treetensornetworks/projttns/abstractprojttn.jl +++ b/src/treetensornetworks/projttns/abstractprojttn.jl @@ -68,11 +68,7 @@ end projected_operator_tensors(P::AbstractProjTTN) = error("Not implemented.") function contract(P::AbstractProjTTN, v::ITensor) - itensor_map = projected_operator_tensors(P) - for t in itensor_map - v *= t - end - return v + return foldl(*, projected_operator_tensors(P); init=v) end function product(P::AbstractProjTTN, v::ITensor) From a90bb0372ec74e65ec9c4783a8cb93cb73447940 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Thu, 1 Feb 2024 20:29:04 -0500 Subject: [PATCH 21/27] Access position in copy via accessor function. --- src/treetensornetworks/projttns/projouterprodttn.jl | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/src/treetensornetworks/projttns/projouterprodttn.jl b/src/treetensornetworks/projttns/projouterprodttn.jl index 278a66cf..3052e0bd 100644 --- a/src/treetensornetworks/projttns/projouterprodttn.jl +++ b/src/treetensornetworks/projttns/projouterprodttn.jl @@ -22,7 +22,7 @@ end function copy(P::ProjOuterProdTTN) return ProjOuterProdTTN( - P.pos, copy(internal_state(P)), copy(operator(P)), copy(environments(P)) + pos(P), copy(internal_state(P)), copy(operator(P)), copy(environments(P)) ) end @@ -109,7 +109,3 @@ end function contract(P::ProjOuterProdTTN, x::ITensor) return conj(contract_ket(P, x)) * contract_ket(P, ITensor(true)) end - -function product(P::ProjOuterProdTTN, v::ITensor) - return noprime(contract(P, v)) -end From 88351d75544c016febb313657ae0d3af7698698f Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Thu, 1 Feb 2024 20:31:30 -0500 Subject: [PATCH 22/27] Access vertextype by function instead of via where in type definition. --- src/treetensornetworks/projttns/projttnsum.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/treetensornetworks/projttns/projttnsum.jl b/src/treetensornetworks/projttns/projttnsum.jl index 59c1a5ca..cd1cdf72 100644 --- a/src/treetensornetworks/projttns/projttnsum.jl +++ b/src/treetensornetworks/projttns/projttnsum.jl @@ -3,8 +3,8 @@ ProjTTNSum """ struct ProjTTNSum{V} <: AbstractProjTTN{V} terms::Vector{AbstractProjTTN{V}} - function ProjTTNSum(terms::Vector{<:AbstractProjTTN{V}}) where {V} - return new{V}(terms) + function ProjTTNSum(terms::Vector{<:AbstractProjTTN}) + return new{eltype(terms)}(terms) end end From 59cbf791a180d2c529aee6faa7930e8d04560652 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Fri, 2 Feb 2024 06:57:21 -0500 Subject: [PATCH 23/27] Fix DMRG variational apply test and type signature in type definition. --- src/treetensornetworks/projttns/projouterprodttn.jl | 4 +++- src/treetensornetworks/projttns/projttnsum.jl | 6 +++--- test/test_treetensornetworks/test_solvers/test_contract.jl | 4 ++++ 3 files changed, 10 insertions(+), 4 deletions(-) diff --git a/src/treetensornetworks/projttns/projouterprodttn.jl b/src/treetensornetworks/projttns/projouterprodttn.jl index 3052e0bd..27f532e6 100644 --- a/src/treetensornetworks/projttns/projouterprodttn.jl +++ b/src/treetensornetworks/projttns/projouterprodttn.jl @@ -106,6 +106,8 @@ function contract_ket(P::ProjOuterProdTTN, v::ITensor) return v end +# ToDo: debug this, not yet working +# probably function contract(P::ProjOuterProdTTN, x::ITensor) - return conj(contract_ket(P, x)) * contract_ket(P, ITensor(true)) + return conj(contract_ket(P, dag(x))) * contract_ket(P, ITensor(true)) end diff --git a/src/treetensornetworks/projttns/projttnsum.jl b/src/treetensornetworks/projttns/projttnsum.jl index cd1cdf72..31eed8f5 100644 --- a/src/treetensornetworks/projttns/projttnsum.jl +++ b/src/treetensornetworks/projttns/projttnsum.jl @@ -2,9 +2,9 @@ ProjTTNSum """ struct ProjTTNSum{V} <: AbstractProjTTN{V} - terms::Vector{AbstractProjTTN{V}} - function ProjTTNSum(terms::Vector{<:AbstractProjTTN}) - return new{eltype(terms)}(terms) + terms::Vector{T} where {T<:AbstractProjTTN{V}} + function ProjTTNSum(terms::Vector{T}) where {V,T<:AbstractProjTTN{V}} + return new{V}(terms) end end diff --git a/test/test_treetensornetworks/test_solvers/test_contract.jl b/test/test_treetensornetworks/test_solvers/test_contract.jl index 46eadf11..4be8a574 100644 --- a/test/test_treetensornetworks/test_solvers/test_contract.jl +++ b/test/test_treetensornetworks/test_solvers/test_contract.jl @@ -24,6 +24,10 @@ using Test # Test basic usage with default parameters Hpsi = apply(H, psi; alg="fit", init=psi, nsweeps=1) @test inner(psi, Hpsi) ≈ inner(psi', H, psi) atol = 1E-5 + # Test variational compression via DMRG + 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 basic usage for use with multiple ProjOuterProdTTN with default parameters # BLAS.axpy-like test From 559c94805a9f4e667ad8b87795a854a443b97e7b Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Fri, 2 Feb 2024 08:05:20 -0500 Subject: [PATCH 24/27] Fix test. --- test/test_treetensornetworks/test_solvers/test_linsolve.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_treetensornetworks/test_solvers/test_linsolve.jl b/test/test_treetensornetworks/test_solvers/test_linsolve.jl index 72291642..168471bb 100644 --- a/test/test_treetensornetworks/test_solvers/test_linsolve.jl +++ b/test/test_treetensornetworks/test_solvers/test_linsolve.jl @@ -41,7 +41,7 @@ using Random x_c = random_mps(s; states, internal_inds_space=4) + 0.1im * random_mps(s; states, internal_inds_space=2) - b = apply(H, x_c; alg="fit", nsweeps=3) #cutoff is unsupported kwarg for apply/contract + b = apply(H, x_c; alg="fit", nsweeps=3, init=x_c) #cutoff is unsupported kwarg for apply/contract x0 = random_mps(s; states, internal_inds_space=10) x = @test_broken linsolve( From cb8a78fb30ad5076e91e251ca6fd9c1f2ca1ef9b Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Fri, 2 Feb 2024 08:10:15 -0500 Subject: [PATCH 25/27] Fix TimeDependentSum. --- src/treetensornetworks/solvers/solver_utils.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/treetensornetworks/solvers/solver_utils.jl b/src/treetensornetworks/solvers/solver_utils.jl index 61b228bb..552ba5aa 100644 --- a/src/treetensornetworks/solvers/solver_utils.jl +++ b/src/treetensornetworks/solvers/solver_utils.jl @@ -28,7 +28,7 @@ struct TimeDependentSum{S,T} f::Vector{S} H0::T end -TimeDependentSum(f::Vector, H0::ProjTTNSum) = TimeDependentSum(f, H0.pm) +TimeDependentSum(f::Vector, H0::ProjTTNSum) = TimeDependentSum(f, terms(H0)) Base.length(H::TimeDependentSum) = length(H.f) function Base.:*(c::Number, H::TimeDependentSum) From e42ef5a384f8d95823ed6fa601138f93864328be Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Fri, 2 Feb 2024 10:10:40 -0500 Subject: [PATCH 26/27] Add prefactors to ProjTTNSum, optimize ProjOuterProdTTN --- src/solvers/contract.jl | 2 +- .../projttns/abstractprojttn.jl | 3 +++ .../projttns/projouterprodttn.jl | 6 ++--- src/treetensornetworks/projttns/projttnsum.jl | 27 ++++++++++++++----- .../test_solvers/test_contract.jl | 13 +++++++-- 5 files changed, 38 insertions(+), 13 deletions(-) diff --git a/src/solvers/contract.jl b/src/solvers/contract.jl index 6a659f48..cf5cddd2 100644 --- a/src/solvers/contract.jl +++ b/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 diff --git a/src/treetensornetworks/projttns/abstractprojttn.jl b/src/treetensornetworks/projttns/abstractprojttn.jl index 2739d45e..7ac54106 100644 --- a/src/treetensornetworks/projttns/abstractprojttn.jl +++ b/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) diff --git a/src/treetensornetworks/projttns/projouterprodttn.jl b/src/treetensornetworks/projttns/projouterprodttn.jl index 27f532e6..d507202e 100644 --- a/src/treetensornetworks/projttns/projouterprodttn.jl +++ b/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 diff --git a/src/treetensornetworks/projttns/projttnsum.jl b/src/treetensornetworks/projttns/projttnsum.jl index 31eed8f5..bbe93759 100644 --- a/src/treetensornetworks/projttns/projttnsum.jl +++ b/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)) diff --git a/test/test_treetensornetworks/test_solvers/test_contract.jl b/test/test_treetensornetworks/test_solvers/test_contract.jl index 4be8a574..42e0f068 100644 --- a/test/test_treetensornetworks/test_solvers/test_contract.jl +++ b/test/test_treetensornetworks/test_solvers/test_contract.jl @@ -28,6 +28,10 @@ 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 @@ -35,10 +39,15 @@ using Test 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 From 1ce884db39daa2e5271220571b0bf2bed8ea8c88 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Fri, 2 Feb 2024 11:48:18 -0500 Subject: [PATCH 27/27] Fix factors and edit tests. --- .../projttns/abstractprojttn.jl | 2 +- src/treetensornetworks/projttns/projttnsum.jl | 25 ++++++++++++------- .../test_solvers/test_contract.jl | 21 ++++++++++------ 3 files changed, 30 insertions(+), 18 deletions(-) diff --git a/src/treetensornetworks/projttns/abstractprojttn.jl b/src/treetensornetworks/projttns/abstractprojttn.jl index 7ac54106..fc8eeb7e 100644 --- a/src/treetensornetworks/projttns/abstractprojttn.jl +++ b/src/treetensornetworks/projttns/abstractprojttn.jl @@ -102,7 +102,7 @@ function Base.eltype(P::AbstractProjTTN)::Type end vertextype(::Type{<:AbstractProjTTN{V}}) where {V} = V -vertextype(p::AbstractProjTTN) = typeof(p) +vertextype(p::AbstractProjTTN) = vertextype(typeof(p)) function Base.size(P::AbstractProjTTN)::Tuple{Int,Int} d = 1 diff --git a/src/treetensornetworks/projttns/projttnsum.jl b/src/treetensornetworks/projttns/projttnsum.jl index bbe93759..b731b989 100644 --- a/src/treetensornetworks/projttns/projttnsum.jl +++ b/src/treetensornetworks/projttns/projttnsum.jl @@ -12,13 +12,13 @@ end terms(P::ProjTTNSum) = P.terms factors(P::ProjTTNSum) = P.factors -copy(P::ProjTTNSum) = ProjTTNSum(copy.(terms(P))) +copy(P::ProjTTNSum) = ProjTTNSum(copy.(terms(P)), copy(factors(P))) function ProjTTNSum(operators::Vector{<:AbstractProjTTN}) - return ProjTTNSum(operators, fill(1, length(operators))) + return ProjTTNSum(operators, fill(one(Bool), length(operators))) end function ProjTTNSum(operators::Vector{<:AbstractTTN}) - return ProjTTNSum(ProjTTN.(operators), fill(1, length(operators))) + return ProjTTNSum(ProjTTN.(operators)) end on_edge(P::ProjTTNSum) = on_edge(terms(P)[1]) @@ -26,7 +26,7 @@ on_edge(P::ProjTTNSum) = on_edge(terms(P)[1]) nsite(P::ProjTTNSum) = nsite(terms(P)[1]) function set_nsite(Ps::ProjTTNSum, nsite) - return ProjTTNSum(map(p -> set_nsite(p, nsite), terms(Ps))) + return ProjTTNSum(map(p -> set_nsite(p, nsite), terms(Ps)), factors(Ps)) end underlying_graph(P::ProjTTNSum) = underlying_graph(terms(P)[1]) @@ -41,15 +41,19 @@ internal_edges(P::ProjTTNSum) = internal_edges(terms(P)[1]) product(P::ProjTTNSum, v::ITensor) = noprime(contract(P, v)) -contract(P::ProjTTNSum, v::ITensor) = - mapreduce(+, zip(factors(P), terms(P))) do (f, p) +function contract(P::ProjTTNSum, v::ITensor) + res = mapreduce(+, zip(factors(P), terms(P))) do (f, p) f * contract(p, v) end + return res +end -contract_ket(P::ProjTTNSum, v::ITensor) = - mapreduce(+, zip(factors(P), terms(P))) do (f, p) +function contract_ket(P::ProjTTNSum, v::ITensor) + res = mapreduce(+, zip(factors(P), terms(P))) do (f, p) f * contract_ket(p, v) end + return res +end function Base.eltype(P::ProjTTNSum) return mapreduce(eltype, promote_type, terms(P)) @@ -60,5 +64,8 @@ end Base.size(P::ProjTTNSum) = size(terms(P)[1]) function position(P::ProjTTNSum, psi::AbstractTTN, pos) - return ProjTTNSum(map(M -> position(M, psi, pos), terms(P))) + theterms = map(M -> position(M, psi, pos), terms(P)) + #@show typeof(theterms) + #@show factors(P) + return ProjTTNSum(theterms, factors(P)) end diff --git a/test/test_treetensornetworks/test_solvers/test_contract.jl b/test/test_treetensornetworks/test_solvers/test_contract.jl index 42e0f068..49f79e57 100644 --- a/test/test_treetensornetworks/test_solvers/test_contract.jl +++ b/test/test_treetensornetworks/test_solvers/test_contract.jl @@ -29,8 +29,8 @@ using Test 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,)) + Hfit = ProjTTNSum([ProjOuterProdTTN(psi', H), ProjOuterProdTTN(psi', H)], [-0.2, -0.8]) + Hpsi_via_dmrg = dmrg(Hfit, psi; nsweeps=1, updater_kwargs=(; which_eigval=:SR,)) @test abs(inner(Hpsi_via_dmrg, Hpsi / norm(Hpsi))) ≈ 1 atol = 1E-4 # Test basic usage for use with multiple ProjOuterProdTTN with default parameters @@ -38,15 +38,20 @@ using Test os_id = OpSum() os_id += -1, "Id", 1, "Id", 2 minus_identity = mpo(os_id, s) + os_id = OpSum() + os_id += +1, "Id", 1, "Id", 2 + identity = mpo(os_id, s) Hpsi = ITensorNetworks.sum_apply( - [(H, psi), (minus_identity, psi)]; alg="fit", init=psi, nsweeps=1 + [(H, psi), (minus_identity, psi)]; alg="fit", init=psi, nsweeps=3 ) @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 + # Test the above via DMRG + # ToDo: Investigate why this is broken + Hfit = ProjTTNSum([ProjOuterProdTTN(psi', H), ProjOuterProdTTN(psi', identity)], [-1, 1]) + Hpsi_normalized = ITensorNetworks.dmrg( + Hfit, psi; nsweeps=3, updater_kwargs=(; which_eigval=:SR) + ) + @test_broken abs(inner(Hpsi, (Hpsi_normalized) / norm(Hpsi))) ≈ 1 atol = 1E-5 # # Change "top" indices of MPO to be a different set