From 367872e31fb357e51af94fab5db3af1136e2e112 Mon Sep 17 00:00:00 2001 From: Hiroshi Shinaoka Date: Tue, 26 Nov 2024 11:47:04 -0500 Subject: [PATCH] Weakdep on PartitionedMPSs.jl --- Project.toml | 11 +- .../QuantiscPartitionedMPSsExt.jl | 187 ++++++++++++++++++ src/Quantics.jl | 1 - src/mul.jl | 64 ------ src/util.jl | 107 ---------- test/mul_tests.jl | 8 +- test/util_tests.jl | 10 +- 7 files changed, 204 insertions(+), 184 deletions(-) create mode 100644 ext/QuantiscPartitionedMPSsExt/QuantiscPartitionedMPSsExt.jl diff --git a/Project.toml b/Project.toml index 3b86143..d53eaa8 100644 --- a/Project.toml +++ b/Project.toml @@ -9,18 +9,23 @@ FastMPOContractions = "f6e391d2-8ffa-4d7a-98cd-7e70024481ca" ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2" ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -ProjMPSs = "3cac4759-b6c4-45be-a543-65afe6a1360a" QuanticsTCI = "b11687fd-3a1c-4c41-97d0-998ab401d50e" SparseIR = "4fe2279e-80f0-4adb-8463-ee114ff56b7d" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" TCIITensorConversion = "9f0aa9f4-9415-4e6a-8795-331ebf40aa04" +[weakdeps] +PartitionedMPSs = "17ce1de9-5131-45e3-8a48-9723b6e2dc23" + +[extensions] +QuantiscPartitionedMPSsExt = "PartitionedMPSs" + [compat] EllipsisNotation = "1" FastMPOContractions = "^0.2.2" ITensorMPS = "0.3.1" ITensors = "0.7" -ProjMPSs = "0.4.1" +PartitionedMPSs = "0.5.2" QuanticsTCI = "0.7.0" SparseIR = "^0.96, 0.97, 1" StaticArrays = "1" @@ -34,4 +39,4 @@ ReTestItems = "817f1d60-ba6b-4fd5-9520-3cf149f6a823" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Test", "Random", "ReTestItems", "Aqua"] +test = ["Test", "Random", "ReTestItems", "Aqua", "PartitionedMPSs"] diff --git a/ext/QuantiscPartitionedMPSsExt/QuantiscPartitionedMPSsExt.jl b/ext/QuantiscPartitionedMPSsExt/QuantiscPartitionedMPSsExt.jl new file mode 100644 index 0000000..8d3cbd6 --- /dev/null +++ b/ext/QuantiscPartitionedMPSsExt/QuantiscPartitionedMPSsExt.jl @@ -0,0 +1,187 @@ +module QuantiscPartitionedMPSsExt + +using ITensors +import ITensors +using ITensors.SiteTypes: siteinds +import ITensors.NDTensors: Tensor, BlockSparseTensor, blockview +using ITensorMPS: ITensorMPS, MPS, MPO, AbstractMPS +using ITensorMPS: findsite, linkinds, linkind, findsites + +import Quantics: Quantics, _find_site_allplevs, combinesites, extractdiagonal +import PartitionedMPSs: PartitionedMPSs, SubDomainMPS, PartitionedMPS, isprojectedat, project + +function Quantics.makesitediagonal(subdmps::SubDomainMPS, site::Index) + return _makesitediagonal(subdmps, site; baseplev=0) +end + +function Quantics.makesitediagonal(subdmps::SubDomainMPS, sites::AbstractVector{Index}) + return _makesitediagonal(subdmps, sites; baseplev=0) +end + + +function Quantics.makesitediagonal(subdmps::SubDomainMPS, tag::String) + mps_diagonal = Quantics.makesitediagonal(MPS(subdmps), tag) + subdmps_diagonal = SubDomainMPS(mps_diagonal) + + target_sites = Quantics.findallsiteinds_by_tag( + unique(ITensors.noprime.(Iterators.flatten(siteinds(subdmps)))); tag=tag + ) + + newproj = deepcopy(subdmps.projector) + for s in target_sites + if isprojectedat(subdmps.projector, s) + newproj[ITensors.prime(s)] = newproj[s] + end + end + + return project(subdmps_diagonal, newproj) +end + +function _makesitediagonal( + subdmps::SubDomainMPS, sites::AbstractVector{Index{IndsT}}; baseplev=0 +) where {IndsT} + M_ = deepcopy(MPO(collect(MPS(subdmps)))) + for site in sites + target_site::Int = only(findsites(M_, site)) + M_[target_site] = _asdiagonal(M_[target_site], site; baseplev=baseplev) + end + return project(M_, subdmps.projector) +end + +function _makesitediagonal(subdmps::SubDomainMPS, site::Index; baseplev=0) + return _makesitediagonal(subdmps, [site]; baseplev=baseplev) +end + + +function Quantics.extractdiagonal( + subdmps::SubDomainMPS, sites::AbstractVector{Index{IndsT}} +) where {IndsT} + tensors = collect(subdmps.data) + for i in eachindex(tensors) + for site in intersect(sites, ITensors.inds(tensors[i])) + sitewithallplevs = _find_site_allplevs(tensors[i], site) + tensors[i] = if length(sitewithallplevs) > 1 + tensors[i] = Quantics._extract_diagonal(tensors[i], sitewithallplevs...) + else + tensors[i] + end + end + end + + projector = deepcopy(subdmps.projector) + for site in sites + if site' in keys(projector.data) + delete!(projector.data, site') + end + end + return SubDomainMPS(MPS(tensors), projector) +end + +function Quantics.extractdiagonal(subdmps::SubDomainMPS, site::Index{IndsT}) where {IndsT} + return Quantics.extractdiagonal(subdmps, [site]) +end + +function Quantics.extractdiagonal(subdmps::SubDomainMPS, tag::String)::SubDomainMPS + targetsites = Quantics.findallsiteinds_by_tag( + unique(ITensors.noprime.(PartitionedMPSs._allsites(subdmps))); tag=tag + ) + return Quantics.extractdiagonal(subdmps, targetsites) +end + +function Quantics.rearrange_siteinds(subdmps::SubDomainMPS, sites) + return PartitionedMPSs.rearrange_siteinds(subdmps, sites) +end + +function Quantics.rearrange_siteinds(partmps::PartitionedMPS, sites) + return PartitionedMPSs.rearrange_siteinds(partmps, sites) +end + + +""" +Make the PartitionedMPS diagonal for a given site index `s` by introducing a dummy index `s'`. +""" +function Quantics.makesitediagonal(obj::PartitionedMPS, site) + return PartitionedMPS([ + _makesitediagonal(prjmps, site; baseplev=baseplev) for prjmps in values(obj) + ]) +end + +function _makesitediagonal(obj::PartitionedMPS, site; baseplev=0) + return PartitionedMPS([ + _makesitediagonal(prjmps, site; baseplev=baseplev) for prjmps in values(obj) + ]) +end + +""" +Extract diagonal of the PartitionedMPS for `s`, `s'`, ... for a given site index `s`, +where `s` must have a prime level of 0. +""" +function Quantics.extractdiagonal(obj::PartitionedMPS, site) + return PartitionedMPS([extractdiagonal(prjmps, site) for prjmps in values(obj)]) +end + + +""" +By default, elementwise multiplication will be performed. +""" +function Quantics.automul( + M1::PartitionedMPS, + M2::PartitionedMPS; + tag_row::String="", + tag_shared::String="", + tag_col::String="", + alg="naive", + maxdim=typemax(Int), + cutoff=1e-25, + kwargs..., +) + all(length.(siteinds(M1)) .== 1) || error("M1 should have only 1 site index per site") + all(length.(siteinds(M2)) .== 1) || error("M2 should have only 1 site index per site") + + sites_row = _findallsiteinds_by_tag(M1; tag=tag_row) + sites_shared = _findallsiteinds_by_tag(M1; tag=tag_shared) + sites_col = _findallsiteinds_by_tag(M2; tag=tag_col) + sites_matmul = Set(Iterators.flatten([sites_row, sites_shared, sites_col])) + + sites1 = only.(siteinds(M1)) + sites1_ewmul = setdiff(only.(siteinds(M1)), sites_matmul) + sites2_ewmul = setdiff(only.(siteinds(M2)), sites_matmul) + sites2_ewmul == sites1_ewmul || error("Invalid sites for elementwise multiplication") + + M1 = _makesitediagonal(M1, sites1_ewmul; baseplev=1) + M2 = _makesitediagonal(M2, sites2_ewmul; baseplev=0) + + sites_M1_diag = [collect(x) for x in siteinds(M1)] + sites_M2_diag = [collect(x) for x in siteinds(M2)] + + M1 = Quantics.rearrange_siteinds(M1, combinesites(sites_M1_diag, sites_row, sites_shared)) + + M2 = Quantics.rearrange_siteinds(M2, combinesites(sites_M2_diag, sites_shared, sites_col)) + + M = PartitionedMPSs.contract(M1, M2; alg=alg, kwargs...) + + M = extractdiagonal(M, sites1_ewmul) + + ressites = Vector{eltype(siteinds(M1)[1])}[] + for s in siteinds(M) + s_ = unique(ITensors.noprime.(s)) + if length(s_) == 1 + push!(ressites, s_) + else + if s_[1] ∈ sites1 + push!(ressites, [s_[1]]) + push!(ressites, [s_[2]]) + else + push!(ressites, [s_[2]]) + push!(ressites, [s_[1]]) + end + end + end + return PartitionedMPSs.truncate(Quantics.rearrange_siteinds(M, ressites); cutoff=cutoff, maxdim=maxdim) +end + +function _findallsiteinds_by_tag(M::PartitionedMPS; tag=tag) + return Quantics.findallsiteinds_by_tag(only.(siteinds(M)); tag=tag) +end + +end diff --git a/src/Quantics.jl b/src/Quantics.jl index c683f0e..bfdda89 100644 --- a/src/Quantics.jl +++ b/src/Quantics.jl @@ -6,7 +6,6 @@ using ITensors.SiteTypes: siteinds import ITensors.NDTensors: Tensor, BlockSparseTensor, blockview using ITensorMPS: ITensorMPS, MPS, MPO, AbstractMPS using ITensorMPS: findsite, linkinds, linkind, findsites -import ProjMPSs: ProjMPSs, ProjMPS, BlockedMPS, isprojectedat, project import SparseIR: Fermionic, Bosonic, Statistics import LinearAlgebra: I diff --git a/src/mul.jl b/src/mul.jl index aeea43a..b875a9e 100644 --- a/src/mul.jl +++ b/src/mul.jl @@ -162,68 +162,4 @@ function automul(M1::MPS, M2::MPS; tag_row::String="", tag_shared::String="", end return asMPS(M) -end - - -""" -By default, elementwise multiplication will be performed. -""" -function automul( - M1::BlockedMPS, - M2::BlockedMPS; - tag_row::String="", - tag_shared::String="", - tag_col::String="", - alg="naive", - maxdim=typemax(Int), - cutoff=1e-25, - kwargs..., -) - all(length.(siteinds(M1)) .== 1) || error("M1 should have only 1 site index per site") - all(length.(siteinds(M2)) .== 1) || error("M2 should have only 1 site index per site") - - sites_row = _findallsiteinds_by_tag(M1; tag=tag_row) - sites_shared = _findallsiteinds_by_tag(M1; tag=tag_shared) - sites_col = _findallsiteinds_by_tag(M2; tag=tag_col) - sites_matmul = Set(Iterators.flatten([sites_row, sites_shared, sites_col])) - - sites1 = only.(siteinds(M1)) - sites1_ewmul = setdiff(only.(siteinds(M1)), sites_matmul) - sites2_ewmul = setdiff(only.(siteinds(M2)), sites_matmul) - sites2_ewmul == sites1_ewmul || error("Invalid sites for elementwise multiplication") - - M1 = _makesitediagonal(M1, sites1_ewmul; baseplev=1) - M2 = _makesitediagonal(M2, sites2_ewmul; baseplev=0) - - sites_M1_diag = [collect(x) for x in siteinds(M1)] - sites_M2_diag = [collect(x) for x in siteinds(M2)] - - M1 = rearrange_siteinds(M1, combinesites(sites_M1_diag, sites_row, sites_shared)) - - M2 = rearrange_siteinds(M2, combinesites(sites_M2_diag, sites_shared, sites_col)) - - M = ProjMPSs.contract(M1, M2; alg=alg, kwargs...) - - M = extractdiagonal(M, sites1_ewmul) - - ressites = Vector{eltype(siteinds(M1)[1])}[] - for s in siteinds(M) - s_ = unique(ITensors.noprime.(s)) - if length(s_) == 1 - push!(ressites, s_) - else - if s_[1] ∈ sites1 - push!(ressites, [s_[1]]) - push!(ressites, [s_[2]]) - else - push!(ressites, [s_[2]]) - push!(ressites, [s_[1]]) - end - end - end - return ProjMPSs.truncate(rearrange_siteinds(M, ressites); cutoff=cutoff, maxdim=maxdim) -end - -function _findallsiteinds_by_tag(M::BlockedMPS; tag=tag) - return Quantics.findallsiteinds_by_tag(only.(siteinds(M)); tag=tag) end \ No newline at end of file diff --git a/src/util.jl b/src/util.jl index 4b79d54..550d35a 100644 --- a/src/util.jl +++ b/src/util.jl @@ -533,46 +533,6 @@ function makesitediagonal(M::AbstractMPS, tag::String)::MPS return MPS(collect(M_)) end -function _makesitediagonal( - projmps::ProjMPS, sites::AbstractVector{Index{IndsT}}; baseplev=0 -) where {IndsT} - M_ = deepcopy(MPO(collect(MPS(projmps)))) - for site in sites - target_site::Int = only(findsites(M_, site)) - M_[target_site] = _asdiagonal(M_[target_site], site; baseplev=baseplev) - end - return project(M_, projmps.projector) -end - -function _makesitediagonal(projmps::ProjMPS, site::Index; baseplev=0) - return _makesitediagonal(projmps, [site]; baseplev=baseplev) -end - -function makesitediagonal(projmps::ProjMPS, site::Index) - return _makesitediagonal(projmps, site; baseplev=0) -end - -function makesitediagonal(projmps::ProjMPS, sites::AbstractVector{Index}) - return _makesitediagonal(projmps, sites; baseplev=0) -end - -function makesitediagonal(projmps::ProjMPS, tag::String) - mps_diagonal = Quantics.makesitediagonal(MPS(projmps), tag) - projmps_diagonal = ProjMPS(mps_diagonal) - - target_sites = Quantics.findallsiteinds_by_tag( - unique(ITensors.noprime.(Iterators.flatten(siteinds(projmps)))); tag=tag - ) - - newproj = deepcopy(projmps.projector) - for s in target_sites - if isprojectedat(projmps.projector, s) - newproj[ITensors.prime(s)] = newproj[s] - end - end - - return project(projmps_diagonal, newproj) -end # FIXME: may be type unstable function _find_site_allplevs(tensor::ITensor, site::Index; maxplev=10) @@ -583,48 +543,6 @@ function _find_site_allplevs(tensor::ITensor, site::Index; maxplev=10) ] end -function extractdiagonal( - projmps::ProjMPS, sites::AbstractVector{Index{IndsT}} -) where {IndsT} - tensors = collect(projmps.data) - for i in eachindex(tensors) - for site in intersect(sites, ITensors.inds(tensors[i])) - sitewithallplevs = _find_site_allplevs(tensors[i], site) - tensors[i] = if length(sitewithallplevs) > 1 - tensors[i] = Quantics._extract_diagonal(tensors[i], sitewithallplevs...) - else - tensors[i] - end - end - end - - projector = deepcopy(projmps.projector) - for site in sites - if site' in keys(projector.data) - delete!(projector.data, site') - end - end - return ProjMPS(MPS(tensors), projector) -end - -function extractdiagonal(projmps::ProjMPS, site::Index{IndsT}) where {IndsT} - return Quantics.extractdiagonal(projmps, [site]) -end - -function extractdiagonal(projmps::ProjMPS, tag::String)::ProjMPS - targetsites = Quantics.findallsiteinds_by_tag( - unique(ITensors.noprime.(ProjMPSs._allsites(projmps))); tag=tag - ) - return extractdiagonal(projmps, targetsites) -end - -function rearrange_siteinds(projmps::ProjMPS, sites) - return ProjMPSs.rearrange_siteinds(projmps, sites) -end - -function rearrange_siteinds(bmps::BlockedMPS, sites) - return ProjMPSs.rearrange_siteinds(bmps, sites) -end """ @@ -677,28 +595,3 @@ function _apply(A::MPO, Ψ::MPS; alg::String="fit", cutoff::Real=1e-25, kwargs.. AΨ = noprime.(FastMPOContractions.contract_mpo_mpo(A, asMPO(Ψ); alg, kwargs...)) MPS(collect(AΨ)) end - - - -""" -Make the BlockedMPS diagonal for a given site index `s` by introducing a dummy index `s'`. -""" -function makesitediagonal(obj::BlockedMPS, site) - return BlockedMPS([ - _makesitediagonal(prjmps, site; baseplev=baseplev) for prjmps in values(obj) - ]) -end - -function _makesitediagonal(obj::BlockedMPS, site; baseplev=0) - return BlockedMPS([ - _makesitediagonal(prjmps, site; baseplev=baseplev) for prjmps in values(obj) - ]) -end - -""" -Extract diagonal of the BlockedMPS for `s`, `s'`, ... for a given site index `s`, -where `s` must have a prime level of 0. -""" -function extractdiagonal(obj::BlockedMPS, site) - return BlockedMPS([extractdiagonal(prjmps, site) for prjmps in values(obj)]) -end \ No newline at end of file diff --git a/test/mul_tests.jl b/test/mul_tests.jl index 69d3dd5..6b21f09 100644 --- a/test/mul_tests.jl +++ b/test/mul_tests.jl @@ -120,11 +120,11 @@ end @testitem "mul_tests.jl/batchedmatmul" begin using Test + import PartitionedMPSs: PartitionedMPSs, SubDomainMPS, PartitionedMPS, isprojectedat, project, Projector import Quantics using ITensors using ITensors.SiteTypes: siteinds using ITensorMPS: random_mps, MPS - import ProjMPSs: ProjMPSs, ProjMPS, project, BlockedMPS, Projector """ Reconstruct 3D matrix @@ -167,7 +167,7 @@ end @test ab_arr ≈ ab_arr_reconst end - @testset "BlockedMPS" begin + @testset "PartitionedMPS" begin @testset "batchedmatmul" for T in [Float64] """ C(x, z, k) = sum_y A(x, y, k) * B(y, z, k) @@ -194,14 +194,14 @@ end ab_arr[:, :, k] .= a_arr[:, :, k] * b_arr[:, :, k] end - a_ = BlockedMPS([ + a_ = PartitionedMPS([ project(a, Projector(Dict(sx[1] => 1, sy[1] => 1))), project(a, Projector(Dict(sx[1] => 1, sy[1] => 2))), project(a, Projector(Dict(sx[1] => 2, sy[1] => 1))), project(a, Projector(Dict(sx[1] => 2, sy[1] => 2))) ]) - b_ = BlockedMPS([ + b_ = PartitionedMPS([ project(b, Projector(Dict(sy[1] => 1, sz[1] => 1))), project(b, Projector(Dict(sy[1] => 1, sz[1] => 2))), project(b, Projector(Dict(sy[1] => 2, sz[1] => 1))), diff --git a/test/util_tests.jl b/test/util_tests.jl index 387d805..df45d56 100644 --- a/test/util_tests.jl +++ b/test/util_tests.jl @@ -1,10 +1,10 @@ @testitem "util.jl" begin using Test + import PartitionedMPSs: PartitionedMPSs, SubDomainMPS, PartitionedMPS, project, isprojectedat import Quantics using ITensors using ITensors.SiteTypes: siteinds using ITensorMPS: randomMPS, randomMPO, random_mps, MPO, MPS - import ProjMPSs: ProjMPSs, ProjMPS, BlockedMPS, project, isprojectedat include("_util.jl") @@ -251,7 +251,7 @@ end end - @testset "ProjMPS" begin + @testset "SubDomainMPS" begin @testset "rearrange_siteinds" begin N = 3 sitesx = [Index(2, "x=$n") for n in 1:N] @@ -261,7 +261,7 @@ Ψ = MPS(collect(_random_mpo(sites))) - prjΨ = ProjMPS(Ψ) + prjΨ = SubDomainMPS(Ψ) prjΨ1 = project(prjΨ, Dict(sitesx[1] => 1)) sitesxy = collect(collect.(zip(sitesx, sitesy))) @@ -273,7 +273,7 @@ prjΨ1_rearranged = Quantics.rearrange_siteinds(prjΨ1, sites_rearranged) @test reduce(*, MPS(prjΨ1)) ≈ reduce(*, MPS(prjΨ1_rearranged)) - @test ProjMPSs.siteinds(prjΨ1_rearranged) == sites_rearranged + @test PartitionedMPSs.siteinds(prjΨ1_rearranged) == sites_rearranged end @testset "makesitediagonal and extractdiagonal" begin @@ -288,7 +288,7 @@ Ψ = MPS(collect(_random_mpo(sites))) - prjΨ = ProjMPS(Ψ) + prjΨ = SubDomainMPS(Ψ) prjΨ1 = project(prjΨ, Dict(sitesx[1] => 1)) prjΨ1_diagonalz = Quantics.makesitediagonal(prjΨ1, "y")