Skip to content

Commit

Permalink
Correct minor bugs and dependency from Quantics.jl
Browse files Browse the repository at this point in the history
  • Loading branch information
gianlucagrosso committed Jan 16, 2025
1 parent adb5abf commit 7942d31
Show file tree
Hide file tree
Showing 12 changed files with 220 additions and 89 deletions.
3 changes: 1 addition & 2 deletions Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,6 @@ julia = "1.6"
[extras]
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
Quantics = "87f76fb3-a40a-40c9-a63c-29fcfe7b7547"

[targets]
test = ["Test", "Random", "Quantics"]
test = ["Test", "Random"]
2 changes: 1 addition & 1 deletion src/PartitionedMPSs.jl
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@ import OrderedCollections: OrderedSet, OrderedDict
using EllipsisNotation
using LinearAlgebra: LinearAlgebra

import ITensors: ITensors, Index, ITensor, dim, inds, qr, commoninds, uniqueinds
import ITensors: ITensors, Index, ITensor, dim, inds, qr, commoninds, uniqueinds, hasplev
import ITensorMPS: ITensorMPS, AbstractMPS, MPS, MPO, siteinds, findsites
import ITensors.TagSets: hastag, hastags

Expand Down
6 changes: 3 additions & 3 deletions src/adaptivemul.jl
Original file line number Diff line number Diff line change
Expand Up @@ -74,9 +74,9 @@ function adaptivecontract(

result_blocks = SubDomainMPS[]
for (p, muls) in patches
prjmpss = [contract(m.a, m.b; alg, cutoff, maxdim, kwargs...) for m in muls]
#patches[p] = +(prjmpss...; alg="fit", cutoff, maxdim)
push!(result_blocks, +(prjmpss...; alg="fit", cutoff, maxdim))
subdmps = [contract(m.a, m.b; alg, cutoff, maxdim, kwargs...) for m in muls]
#patches[p] = +(subdmps...; alg="fit", cutoff, maxdim)
push!(result_blocks, +(subdmps...; alg="fit", cutoff, maxdim))
end

return PartitionedMPS(result_blocks)
Expand Down
6 changes: 3 additions & 3 deletions src/contract.jl
Original file line number Diff line number Diff line change
Expand Up @@ -72,12 +72,12 @@ function projcontract(
_, external_sites = _projector_after_contract(M1, M2)

if !_is_externalsites_compatible_with_projector(external_sites, proj)
error("The projector contains projection onto a site what is not a external sites.")
error("The projector contains projection onto a site that is not an external site.")
end

t1 = time_ns()
# t1 = time_ns()
r = contract(M1, M2; alg, cutoff, maxdim, kwargs...)
t2 = time_ns()
# t2 = time_ns()
#println("contract: $((t2 - t1)*1e-9) s")
return r
end
Expand Down
32 changes: 17 additions & 15 deletions src/partitionedmps.jl
Original file line number Diff line number Diff line change
Expand Up @@ -6,11 +6,11 @@ struct PartitionedMPS
data::OrderedDict{Projector,SubDomainMPS}

function PartitionedMPS(data::AbstractVector{SubDomainMPS})
sites_all = [siteinds(prjmps) for prjmps in data]
sites_all = [siteinds(subdmps) for subdmps in data]
for n in 2:length(data)
Set(sites_all[n]) == Set(sites_all[1]) || error("Sitedims mismatch")
end
isdisjoint([prjmps.projector for prjmps in data]) || error("Projectors are overlapping")
isdisjoint([subdmps.projector for subdmps in data]) || error("Projectors are overlapping")

dict_ = OrderedDict{Projector,SubDomainMPS}(
data[i].projector => data[i] for i in 1:length(data)
Expand Down Expand Up @@ -58,19 +58,19 @@ Base.length(obj::PartitionedMPS) = length(obj.data)
"""
Indexing for PartitionedMPS. This is deprecated and will be removed in the future.
"""
function Base.getindex(bmps::PartitionedMPS, i::Integer)::SubDomainMPS
@warn "Indexing for PartitionedMPS is deprecated. Use getindex(bmps, p::Projector) instead."
return first(Iterators.drop(values(bmps.data), i - 1))
function Base.getindex(partmps::PartitionedMPS, i::Integer)::SubDomainMPS
@warn "Indexing for PartitionedMPS is deprecated. Use getindex(partmps, p::Projector) instead."
return first(Iterators.drop(values(partmps.data), i - 1))
end

Base.getindex(obj::PartitionedMPS, p::Projector) = obj.data[p]

function Base.iterate(bmps::PartitionedMPS, state)
return iterate(bmps.data, state)
function Base.iterate(partmps::PartitionedMPS, state)
return iterate(partmps.data, state)
end

function Base.iterate(bmps::PartitionedMPS)
return iterate(bmps.data)
function Base.iterate(partmps::PartitionedMPS)
return iterate(partmps.data)
end

"""
Expand All @@ -92,11 +92,13 @@ Rearrange the site indices of the PartitionedMPS according to the given order.
If nessecary, tensors are fused or split to match the new order.
"""
function rearrange_siteinds(obj::PartitionedMPS, sites)
return PartitionedMPS([rearrange_siteinds(prjmps, sites) for prjmps in values(obj)])
return PartitionedMPS([rearrange_siteinds(subdmps, sites) for subdmps in values(obj)])
end

function prime::PartitionedMPS, args...; kwargs...)
return PartitionedMPS([prime(prjmps, args...; kwargs...) for prjmps in values.data)])
return PartitionedMPS([
prime(subdmps, args...; kwargs...) for subdmps in values.data)
])
end

"""
Expand Down Expand Up @@ -234,21 +236,21 @@ end
function ITensorMPS.MPO(
obj::PartitionedMPS; cutoff=default_cutoff(), maxdim=default_maxdim()
)::MPO
return MPO(collect(MPS(obj; cutoff=cutoff, maxdim=maxdim, kwargs...)))
return MPO(collect(MPS(obj; cutoff=cutoff, maxdim=maxdim)))
end

"""
Make the PartitionedMPS diagonal for a given site index `s` by introducing a dummy index `s'`.
"""
function makesitediagonal(obj::PartitionedMPS, site)
return PartitionedMPS([
_makesitediagonal(prjmps, site; baseplev=baseplev) for prjmps in values(obj)
_makesitediagonal(subdmps, site; baseplev=baseplev) for subdmps in values(obj)
])
end

function _makesitediagonal(obj::PartitionedMPS, site; baseplev=0)
return PartitionedMPS([
_makesitediagonal(prjmps, site; baseplev=baseplev) for prjmps in values(obj)
_makesitediagonal(subdmps, site; baseplev=baseplev) for subdmps in values(obj)
])
end

Expand All @@ -257,7 +259,7 @@ Extract diagonal of the PartitionedMPS for `s`, `s'`, ... for a given site index
where `s` must have a prime level of 0.
"""
function extractdiagonal(obj::PartitionedMPS, site)
return PartitionedMPS([extractdiagonal(prjmps, site) for prjmps in values(obj)])
return PartitionedMPS([extractdiagonal(subdmps, site) for subdmps in values(obj)])
end

function dist(a::PartitionedMPS, b::PartitionedMPS)
Expand Down
46 changes: 24 additions & 22 deletions src/patching.jl
Original file line number Diff line number Diff line change
Expand Up @@ -5,35 +5,35 @@ If the bond dimension of the result reaches `maxdim`,
perform patching recursively to reduce the bond dimension.
"""
function _add_patching(
prjmpss::AbstractVector{SubDomainMPS};
subdmpss::AbstractVector{SubDomainMPS};
cutoff=0.0,
maxdim=typemax(Int),
alg="fit",
patchorder=Index[],
)::Vector{SubDomainMPS}
if length(unique([prjmps.projector for prjmps in prjmpss])) != 1
if length(unique([sudmps.projector for sudmps in subdmpss])) != 1
error("All SubDomainMPS objects must have the same projector.")
end

# First perform addition upto given maxdim
# TODO: Early termination if the bond dimension reaches maxdim
sum_approx = _add(prjmpss...; alg, cutoff, maxdim)
sum_approx = _add(subdmpss...; alg, cutoff, maxdim)

# If the bond dimension is less than maxdim, return the result
maxbonddim(sum_approx) < maxdim && return [sum_approx]

@assert maxbonddim(sum_approx) == maxdim
# @assert maxbonddim(sum_approx) == maxdim

nextprjidx = _next_projindex(prjmpss[1].projector, patchorder)
nextprjidx = _next_projindex(subdmpss[1].projector, patchorder)

nextprjidx === nothing && return PartitionedMPS(sum_approx)
nextprjidx === nothing && return [sum_approx]

blocks = SubDomainMPS[]
for prjval in 1:ITensors.dim(nextprjidx)
prj_ = prjmpss[1].projector & Projector(nextprjidx => prjval)
prj_ = subdmpss[1].projector & Projector(nextprjidx => prjval)
blocks =
blocks _add_patching(
[project(prjmps, prj_) for prjmps in prjmpss];
[project(sudmps, prj_) for sudmps in subdmpss];
cutoff,
maxdim,
alg,
Expand All @@ -60,13 +60,15 @@ end
Add multiple PartitionedMPS objects.
"""
function add_patching(
bmpss::AbstractVector{PartitionedMPS};
partmps::AbstractVector{PartitionedMPS};
cutoff=0.0,
maxdim=typemax(Int),
alg="fit",
patchorder=Index[],
)::PartitionedMPS
result = _add_patching(union(values(x) for x in bmpss); cutoff, maxdim, alg, patchorder)
result = _add_patching(
union(values(x) for x in partmps); cutoff, maxdim, alg, patchorder
)
return PartitionedMPS(result)
end

Expand All @@ -77,29 +79,29 @@ Do patching recursively to reduce the bond dimension.
If the bond dimension of a SubDomainMPS exceeds `maxdim`, perform patching.
"""
function adaptive_patching(
prjmps::SubDomainMPS, patchorder; cutoff=0.0, maxdim=typemax(Int)
subdmps::SubDomainMPS, patchorder; cutoff=0.0, maxdim=typemax(Int)
)::Vector{SubDomainMPS}
if maxbonddim(prjmps) <= maxdim
return [prjmps]
if maxbonddim(subdmps) <= maxdim
return [subdmps]
end

# If the bond dimension exceeds maxdim, perform patching
refined_prjmpss = SubDomainMPS[]
nextprjidx = _next_projindex(prjmps.projector, patchorder)
refined_subdmpss = SubDomainMPS[]
nextprjidx = _next_projindex(subdmps.projector, patchorder)
if nextprjidx === nothing
return [prjmps]
return [subdmps]
end

for prjval in 1:ITensors.dim(nextprjidx)
prj_ = prjmps.projector & Projector(nextprjidx => prjval)
prjmps_ = truncate(project(prjmps, prj_); cutoff, maxdim)
if maxbonddim(prjmps_) <= maxdim
push!(refined_prjmpss, prjmps_)
prj_ = subdmps.projector & Projector(nextprjidx => prjval)
subdmps_ = truncate(project(subdmps, prj_); cutoff, maxdim)
if maxbonddim(subdmps_) <= maxdim
push!(refined_subdmpss, subdmps_)
else
append!(refined_prjmpss, adaptive_patching(prjmps_, patchorder; cutoff, maxdim))
append!(refined_subdmpss, adaptive_patching(subdmps_, patchorder; cutoff, maxdim))
end
end
return refined_prjmpss
return refined_subdmpss
end

"""
Expand Down
53 changes: 28 additions & 25 deletions src/subdomainmps.jl
Original file line number Diff line number Diff line change
Expand Up @@ -150,9 +150,11 @@ function _add(ψ::AbstractMPS...; alg="fit", cutoff=1e-15, maxdim=typemax(Int),
return +(ITensors.Algorithm(alg), ψ...)
elseif alg == "densitymatrix"
if cutoff < 1e-15
@warn "Cutoff is very small, it may suffer from numerical round errors. The densitymatrix algorithm squares the singular values of the reduce density matrix. Please consider increasing it or using fit algorithm."
@warn "Cutoff is very small, it may suffer from numerical round errors.
The densitymatrix algorithm squares the singular values of the reduce density matrix.
Please consider increasing it or using fit algorithm."
end
return +(ITensors.Algorithm"densitymatrix"(), ψ...; cutoff, maxdim, kwargs...)
return +(ITensors.Algorithm("densitymatrix"), ψ...; cutoff, maxdim, kwargs...)
elseif alg == "fit"
function f(x, y)
return ITensors.truncate(
Expand Down Expand Up @@ -211,39 +213,39 @@ function LinearAlgebra.norm(M::SubDomainMPS)
end

function _makesitediagonal(
SubDomainMPS::SubDomainMPS, sites::AbstractVector{Index{IndsT}}; baseplev=0
obj::SubDomainMPS, sites::AbstractVector{Index{IndsT}}; baseplev=0
) where {IndsT}
M_ = deepcopy(MPO(collect(MPS(SubDomainMPS))))
M_ = deepcopy(MPO(collect(MPS(obj))))
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_, SubDomainMPS.projector)
return project(M_, obj.projector)
end

function _makesitediagonal(SubDomainMPS::SubDomainMPS, site::Index; baseplev=0)
return _makesitediagonal(SubDomainMPS, [site]; baseplev=baseplev)
function _makesitediagonal(obj::SubDomainMPS, site::Index; baseplev=0)
return _makesitediagonal(obj, [site]; baseplev=baseplev)
end

function makesitediagonal(SubDomainMPS::SubDomainMPS, site::Index)
return _makesitediagonal(SubDomainMPS, site; baseplev=0)
function makesitediagonal(obj::SubDomainMPS, site::Index)
return _makesitediagonal(obj, site; baseplev=0)
end

function makesitediagonal(SubDomainMPS::SubDomainMPS, sites::AbstractVector{Index})
return _makesitediagonal(SubDomainMPS, sites; baseplev=0)
function makesitediagonal(obj::SubDomainMPS, sites::AbstractVector{Index})
return _makesitediagonal(obj, sites; baseplev=0)
end

function makesitediagonal(SubDomainMPS::SubDomainMPS, tag::String)
mps_diagonal = makesitediagonal(MPS(SubDomainMPS), tag)
function makesitediagonal(obj::SubDomainMPS, tag::String)
mps_diagonal = makesitediagonal(MPS(obj), tag)
SubDomainMPS_diagonal = SubDomainMPS(mps_diagonal)

target_sites = findallsiteinds_by_tag(
unique(ITensors.noprime.(Iterators.flatten(siteinds(SubDomainMPS)))); tag=tag
unique(ITensors.noprime.(Iterators.flatten(siteinds(obj)))); tag=tag
)

newproj = deepcopy(SubDomainMPS.projector)
newproj = deepcopy(obj.projector)
for s in target_sites
if isprojectedat(SubDomainMPS.projector, s)
if isprojectedat(obj.projector, s)
newproj[ITensors.prime(s)] = newproj[s]
end
end
Expand All @@ -252,7 +254,10 @@ function makesitediagonal(SubDomainMPS::SubDomainMPS, tag::String)
end

# FIXME: may be type unstable
function _find_site_allplevs(tensor::ITensor, site::Index; maxplev=10)
# Gianluca: FIXED (?)
function _find_site_allplevs(
tensor::ITensor, site::Index{T}; maxplev=10
)::Vector{Index{T}} where {T}
ITensors.plev(site) == 0 || error("Site index must be unprimed.")
return [
ITensors.prime(site, plev) for
Expand All @@ -261,9 +266,9 @@ function _find_site_allplevs(tensor::ITensor, site::Index; maxplev=10)
end

function extractdiagonal(
SubDomainMPS::SubDomainMPS, sites::AbstractVector{Index{IndsT}}
obj::SubDomainMPS, sites::AbstractVector{Index{IndsT}}
) where {IndsT}
tensors = collect(SubDomainMPS.data)
tensors = collect(obj.data)
for i in eachindex(tensors)
for site in intersect(sites, ITensors.inds(tensors[i]))
sitewithallplevs = _find_site_allplevs(tensors[i], site)
Expand All @@ -275,7 +280,7 @@ function extractdiagonal(
end
end

projector = deepcopy(SubDomainMPS.projector)
projector = deepcopy(obj.projector)
for site in sites
if site' in keys(projector.data)
delete!(projector.data, site')
Expand All @@ -284,9 +289,7 @@ function extractdiagonal(
return SubDomainMPS(MPS(tensors), projector)
end

function extractdiagonal(SubDomainMPS::SubDomainMPS, tag::String)::SubDomainMPS
targetsites = findallsiteinds_by_tag(
unique(ITensors.noprime.(PartitionedMPSs._allsites(SubDomainMPS))); tag=tag
)
return extractdiagonal(SubDomainMPS, targetsites)
function extractdiagonal(obj::SubDomainMPS, tag::String)::SubDomainMPS
targetsites = findallsiteinds_by_tag(unique(ITensors.noprime.(_allsites(obj))); tag=tag)
return extractdiagonal(obj, targetsites)
end
Loading

0 comments on commit 7942d31

Please sign in to comment.