From da7636e5cd0c5f83fd4dd7b277cb5b0eecfabde9 Mon Sep 17 00:00:00 2001 From: b-kloss Date: Fri, 2 Feb 2024 13:51:58 -0500 Subject: [PATCH] Generalize `ProjTTNSum`, replace `ProjTTNApply` with `ProjOuterProdTTN` (#132) --- src/ITensorNetworks.jl | 2 +- src/exports.jl | 2 +- src/imports.jl | 4 + src/solvers/contract.jl | 9 +- .../projttns/abstractprojttn.jl | 30 ++--- .../projttns/projouterprodttn.jl | 113 ++++++++++++++++++ src/treetensornetworks/projttns/projttn.jl | 18 +++ .../projttns/projttn_apply.jl | 73 ----------- src/treetensornetworks/projttns/projttnsum.jl | 72 ++++++----- src/treetensornetworks/solvers/contract.jl | 80 ++++++++++--- .../solvers/solver_utils.jl | 2 +- .../test_solvers/test_contract.jl | 59 +++++++-- .../test_solvers/test_linsolve.jl | 2 +- 13 files changed, 302 insertions(+), 164 deletions(-) create mode 100644 src/treetensornetworks/projttns/projouterprodttn.jl delete mode 100644 src/treetensornetworks/projttns/projttn_apply.jl 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/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/imports.jl b/src/imports.jl index 71406878..e7071450 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/solvers/contract.jl b/src/solvers/contract.jl index ed55a729..cf5cddd2 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_ket(P, ITensor(one(Bool))), (;) end diff --git a/src/treetensornetworks/projttns/abstractprojttn.jl b/src/treetensornetworks/projttns/abstractprojttn.jl index e3b0e140..fc8eeb7e 100644 --- a/src/treetensornetworks/projttns/abstractprojttn.jl +++ b/src/treetensornetworks/projttns/abstractprojttn.jl @@ -65,30 +65,13 @@ function _separate_first_two(V::Vector) return frst, scnd, rst end -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 - end - # TODO: actually use optimal contraction sequence here - Hv = v - for it in itensor_map - Hv *= it - end - return Hv +projected_operator_tensors(P::AbstractProjTTN) = error("Not implemented.") + +function contract(P::AbstractProjTTN, v::ITensor) + return foldl(*, projected_operator_tensors(P); init=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( @@ -118,6 +101,9 @@ function Base.eltype(P::AbstractProjTTN)::Type return ElType end +vertextype(::Type{<:AbstractProjTTN{V}}) where {V} = V +vertextype(p::AbstractProjTTN) = vertextype(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 new file mode 100644 index 00000000..d507202e --- /dev/null +++ b/src/treetensornetworks/projttns/projouterprodttn.jl @@ -0,0 +1,113 @@ +struct ProjOuterProdTTN{V} <: AbstractProjTTN{V} + pos::Union{Vector{<:V},NamedEdge{V}} + internal_state::TTN{V} + operator::TTN{V} + environments::Dictionary{NamedEdge{V},ITensor} +end + +environments(p::ProjOuterProdTTN) = p.environments +operator(p::ProjOuterProdTTN) = p.operator +underlying_graph(p::ProjOuterProdTTN) = underlying_graph(operator(p)) +pos(p::ProjOuterProdTTN) = p.pos +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}(), + ) +end + +function copy(P::ProjOuterProdTTN) + return ProjOuterProdTTN( + pos(P), copy(internal_state(P)), copy(operator(P)), copy(environments(P)) + ) +end + +function set_nsite(P::ProjOuterProdTTN, nsite) + return P +end + +function shift_position(P::ProjOuterProdTTN, pos) + return ProjOuterProdTTN(pos, internal_state(P), operator(P), environments(P)) +end + +function set_environments(p::ProjOuterProdTTN, environments) + return ProjOuterProdTTN(pos(p), internal_state(p), operator(p), environments) +end + +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::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 + if !haskey(environments(P), e) + if is_leaf(underlying_graph(P), src(e)) + # leaves are easy + env = internal_state(P)[src(e)] * operator(P)[src(e)] * dag(state[src(e)]) + else + # construct by contracting neighbors + neighbor_envs = ITensor[] + for n in setdiff(neighbors(underlying_graph(P), src(e)), [dst(e)]) + P = make_environment(P, state, edgetype(P)(n, src(e))) + push!(neighbor_envs, environment(P, edgetype(P)(n, src(e)))) + end + # manually heuristic for contraction order: two environments, site tensors, then + # other environments + frst, scnd, rst = _separate_first_two(neighbor_envs) + itensor_map = vcat( + 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) + end + P = set_environment(P, e, env) + end + @assert( + hascommoninds(environment(P, e), state[src(e)]), + "Something went wrong, probably re-orthogonalized this edge in the same direction twice!" + ) + return P +end + +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 + itensor_map = Union{ITensor,OneITensor}[] # TODO: will a Hamiltonian TTN tensor ever be a OneITensor? + for j in sites(P) + push!(itensor_map, internal_state(P)[j]) + end + if on_edge(P) + append!(itensor_map, environments) + else + 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 + +function contract_ket(P::ProjOuterProdTTN, v::ITensor) + itensor_map = projected_operator_tensors(P) + for t in itensor_map + v *= t + end + return v +end + +# ToDo: verify conjugation etc. with complex AbstractTTN +function contract(P::ProjOuterProdTTN, x::ITensor) + ket = contract_ket(P, ITensor(one(Bool))) + return (dag(ket) * x) * ket +end diff --git a/src/treetensornetworks/projttns/projttn.jl b/src/treetensornetworks/projttns/projttn.jl index 725650b2..4d9636ab 100644 --- a/src/treetensornetworks/projttns/projttn.jl +++ b/src/treetensornetworks/projttns/projttn.jl @@ -68,3 +68,21 @@ function make_environment(P::ProjTTN, state::AbstractTTN, e::AbstractEdge) ) return P end + +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 deleted file mode 100644 index d4f691b7..00000000 --- a/src/treetensornetworks/projttns/projttn_apply.jl +++ /dev/null @@ -1,73 +0,0 @@ -struct ProjTTNApply{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 - -function ProjTTNApply(init_state::AbstractTTN, operator::AbstractTTN) - return ProjTTNApply( - 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))) -end - -function set_nsite(P::ProjTTNApply, nsite) - return P -end - -function shift_position(P::ProjTTNApply, pos) - return ProjTTNApply(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) -end - -set_environment(p::ProjTTNApply, edge, env) = set_environment!(copy(p), edge, env) -function set_environment!(p::ProjTTNApply, edge, env) - set!(environments(p), edge, env) - return p -end - -function make_environment(P::ProjTTNApply, 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 - 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)]) - else - # construct by contracting neighbors - neighbor_envs = ITensor[] - for n in setdiff(neighbors(underlying_graph(P), src(e)), [dst(e)]) - P = make_environment(P, state, edgetype(P)(n, src(e))) - push!(neighbor_envs, environment(P, edgetype(P)(n, src(e)))) - end - # manually heuristic for contraction order: two environments, site tensors, then - # 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 - ) # no prime here in comparison to the same routine for Projttn - # TODO: actually use optimal contraction sequence here - env = reduce(*, itensor_map) - end - P = set_environment(P, e, env) - end - @assert( - hascommoninds(environment(P, e), state[src(e)]), - "Something went wrong, probably re-orthogonalized this edge in the same direction twice!" - ) - return P -end diff --git a/src/treetensornetworks/projttns/projttnsum.jl b/src/treetensornetworks/projttns/projttnsum.jl index 7ce816f9..b731b989 100644 --- a/src/treetensornetworks/projttns/projttnsum.jl +++ b/src/treetensornetworks/projttns/projttnsum.jl @@ -1,59 +1,71 @@ """ ProjTTNSum """ -struct ProjTTNSum{V} - pm::Vector{ProjTTN{V}} - function ProjTTNSum(pm::Vector{ProjTTN{V}}) where {V} - return new{V}(pm) +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 -copy(P::ProjTTNSum) = ProjTTNSum(copy.(P.pm)) +terms(P::ProjTTNSum) = P.terms +factors(P::ProjTTNSum) = P.factors -ProjTTNSum(ttnos::Vector{<:TTN}) = ProjTTNSum([ProjTTN(M) for M in ttnos]) +copy(P::ProjTTNSum) = ProjTTNSum(copy.(terms(P)), copy(factors(P))) -ProjTTNSum(Ms::TTN{V}...) where {V} = ProjTTNSum([Ms...]) +function ProjTTNSum(operators::Vector{<:AbstractProjTTN}) + return ProjTTNSum(operators, fill(one(Bool), length(operators))) +end +function ProjTTNSum(operators::Vector{<:AbstractTTN}) + return ProjTTNSum(ProjTTN.(operators)) +end -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(p -> set_nsite(p, nsite), terms(Ps)), factors(Ps)) 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) +product(P::ProjTTNSum, v::ITensor) = noprime(contract(P, v)) + +function contract(P::ProjTTNSum, v::ITensor) + res = mapreduce(+, zip(factors(P), terms(P))) do (f, p) + f * contract(p, v) end - return Pv + return res 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])) +function contract_ket(P::ProjTTNSum, v::ITensor) + res = mapreduce(+, zip(factors(P), terms(P))) do (f, p) + f * contract_ket(p, v) end - return elT + return res +end + +function Base.eltype(P::ProjTTNSum) + return mapreduce(eltype, promote_type, terms(P)) end (P::ProjTTNSum)(v::ITensor) = product(P, v) -Base.size(P::ProjTTNSum) = size(P.pm[1]) +Base.size(P::ProjTTNSum) = size(terms(P)[1]) -function position( - P::ProjTTNSum{V}, psi::TTN{V}, pos::Union{Vector{<:V},NamedEdge{V}} -) where {V} - return ProjTTNSum(map(M -> position(M, psi, pos), P.pm)) +function position(P::ProjTTNSum, psi::AbstractTTN, pos) + theterms = map(M -> position(M, psi, pos), terms(P)) + #@show typeof(theterms) + #@show factors(P) + return ProjTTNSum(theterms, factors(P)) end diff --git a/src/treetensornetworks/solvers/contract.jl b/src/treetensornetworks/solvers/contract.jl index 398138e7..90e8c40a 100644 --- a/src/treetensornetworks/solvers/contract.jl +++ b/src/treetensornetworks/solvers/contract.jl @@ -1,30 +1,43 @@ -function contract( +function sum_contract( ::Algorithm"fit", - tn1::AbstractTTN, - tn2::AbstractTTN; - init=random_ttn(flatten_external_indsnetwork(tn1, tn2); link_space=trivial_space(tn1)), - nsweeps=1, + tns::Vector{<:Tuple{<:AbstractTTN,<:AbstractTTN}}; + init, + nsweeps, nsites=2, # used to be default of call to default_sweep_regions updater_kwargs=(;), kwargs..., ) - n = nv(tn1) - n != nv(tn2) && throw( + tn1s = first.(tns) + tn2s = last.(tns) + 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 + return typeof(tn2)([res]) 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 = 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, ProjOuterProdTTN(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 +46,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 +54,10 @@ function contract( return psi end +function contract(a::Algorithm"fit", tn1::AbstractTTN, tn2::AbstractTTN; kwargs...) + return sum_contract(a, [(tn1, tn2)]; kwargs...) +end + """ Overload of `ITensors.contract`. """ @@ -53,7 +68,38 @@ 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...) + 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) 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' + tn12 = sum_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/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) diff --git a/test/test_treetensornetworks/test_solvers/test_contract.jl b/test/test_treetensornetworks/test_solvers/test_contract.jl index 810db6fb..49f79e57 100644 --- a/test/test_treetensornetworks/test_solvers/test_contract.jl +++ b/test/test_treetensornetworks/test_solvers/test_contract.jl @@ -22,31 +22,58 @@ 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 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 whether the interface works for ProjTTNSum with factors + 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 + # BLAS.axpy-like 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=3 + ) + @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 # 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 = 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 @@ -62,9 +89,19 @@ 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, nsweeps=1) @test inner(psi, Hpsi) ≈ inner(psi', H, psi) atol = 1E-5 + # 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, 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 TTN to be a different set # @@ -74,17 +111,17 @@ end H = replaceinds(H, prime(s; links=[]) => t) # Test with nsweeps=2 - Hpsi = apply(H, psi; alg="fit", 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 @@ -102,7 +139,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 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(