Correct minor bugs and dependency from Quantics.jl

Project.toml

@@ -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"]

src/PartitionedMPSs.jl

@@ -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
 

src/adaptivemul.jl

@@ -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)

src/contract.jl

@@ -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

src/partitionedmps.jl

@@ -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)
@@ -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
 
 """
@@ -92,11 +92,11 @@ 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
 
 """
@@ -242,13 +242,13 @@ Make the PartitionedMPS diagonal for a given site index `s` by introducing a dum
 """
 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
 
@@ -257,7 +257,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)

src/patching.jl

@@ -26,7 +26,7 @@ function _add_patching(
 
     nextprjidx = _next_projindex(prjmpss[1].projector, patchorder)
 
-    nextprjidx === nothing && return PartitionedMPS(sum_approx)
+    nextprjidx === nothing && return [sum_approx]
 
     blocks = SubDomainMPS[]
     for prjval in 1:ITensors.dim(nextprjidx)
@@ -60,13 +60,13 @@ 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
 

src/subdomainmps.jl

@@ -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(
@@ -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
@@ -252,7 +254,8 @@ 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
@@ -261,9 +264,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)
@@ -275,7 +278,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')
@@ -284,9 +287,9 @@ function extractdiagonal(
     return SubDomainMPS(MPS(tensors), projector)
 end
 
-function extractdiagonal(SubDomainMPS::SubDomainMPS, tag::String)::SubDomainMPS
+function extractdiagonal(obj::SubDomainMPS, tag::String)::SubDomainMPS
     targetsites = findallsiteinds_by_tag(
-        unique(ITensors.noprime.(PartitionedMPSs._allsites(SubDomainMPS))); tag=tag
+        unique(ITensors.noprime.(_allsites(obj))); tag=tag
     )
-    return extractdiagonal(SubDomainMPS, targetsites)
+    return extractdiagonal(obj, targetsites)
 end

src/util.jl

@@ -105,3 +105,123 @@ function rearrange_siteinds(M::AbstractMPS, sites::Vector{Vector{Index{T}}})::MP
     tensors[end] *= t
     return MPS(tensors)
 end
+
+# A valid tag should not contain "=".
+_valid_tag(tag::String)::Bool = !occursin("=", tag)
+
+"""
+Find sites with the given tag
+
+For tag = `x`, if `sites` contains an Index object with `x`, the function returns a vector containing only its positon.
+
+If not, the function seach for all Index objects with tags `x=1`, `x=2`, ..., and return their positions.
+
+If no Index object is found, an empty vector will be returned.
+"""
+function findallsites_by_tag(
+    sites::Vector{Index{T}}; tag::String="x", maxnsites::Int=1000
+)::Vector{Int} where {T}
+    _valid_tag(tag) || error("Invalid tag: $tag")
+
+    # 1) Check if there is an Index with exactly `tag`
+    idx = findall(hastags(tag), sites)
+    if !isempty(idx)
+        if length(idx) > 1
+            error("Found more than one site index with tag $(tag)!")
+        end
+        return idx
+    end
+
+    # 2) If not found, search for tag=1, tag=2, ...
+    result = Int[]
+    for n in 1:maxnsites
+        tag_ = tag * "=$n"
+        idx = findall(hastags(tag_), sites)
+        if length(idx) == 0
+            break
+        elseif length(idx) > 1
+            error("Found more than one site indices with $(tag_)!")
+        end
+        push!(result, idx[1])
+    end
+    return result
+end
+
+function findallsiteinds_by_tag(
+    sites::AbstractVector{Index{T}}; tag::String="x", maxnsites::Int=1000
+) where {T}
+    _valid_tag(tag) || error("Invalid tag: $tag")
+    positions = findallsites_by_tag(sites; tag=tag, maxnsites=maxnsites)
+    return [sites[p] for p in positions]
+end
+
+function findallsites_by_tag(
+    sites::Vector{Vector{Index{T}}}; tag::String="x", maxnsites::Int=1000
+)::Vector{NTuple{2,Int}} where {T}
+    _valid_tag(tag) || error("Invalid tag: $tag")
+
+    sites_dict = Dict{Index{T},NTuple{2,Int}}()
+    for i in 1:length(sites)
+        for j in 1:length(sites[i])
+            sites_dict[sites[i][j]] = (i, j)
+        end
+    end
+
+    sitesflatten = collect(Iterators.flatten(sites))
+
+    idx_exact = findall(i -> hastags(i, tag) && hasplev(i, 0), sitesflatten)
+    if !isempty(idx_exact)
+        if length(idx_exact) > 1
+            error("Found more than one site index with tag '$tag'!")
+        end
+        # Return a single position
+        return [sites_dict[sitesflatten[only(idx_exact)]]]
+    end
+
+    result = NTuple{2,Int}[]
+    for n in 1:maxnsites
+        tag_ = tag * "=$n"
+        idx = findall(i -> hastags(i, tag_) && hasplev(i, 0), sitesflatten)
+        if length(idx) == 0
+            break
+        elseif length(idx) > 1
+            error("Found more than one site indices with $(tag_)!")
+        end
+
+        push!(result, sites_dict[sitesflatten[only(idx)]])
+    end
+    return result
+end
+
+function findallsiteinds_by_tag(
+    sites::Vector{Vector{Index{T}}}; tag::String="x", maxnsites::Int=1000
+)::Vector{Index{T}} where {T}
+    _valid_tag(tag) || error("Invalid tag: $tag")
+    positions = findallsites_by_tag(sites; tag=tag, maxnsites=maxnsites)
+    return [sites[i][j] for (i, j) in positions]
+end
+
+function makesitediagonal(M::AbstractMPS, tag::String)::MPS
+    M_ = deepcopy(MPO(collect(M)))
+    sites = siteinds(M_)
+
+    target_positions = findallsites_by_tag(siteinds(M_); tag=tag)
+
+    for t in eachindex(target_positions)
+        i, j = target_positions[t]
+        M_[i] = _asdiagonal(M_[i], sites[i][j])
+    end
+
+    return MPS(collect(M_))
+end
+
+function _extract_diagonal(t, site::Index{T}, site2::Index{T}) where {T<:Number}
+    dim(site) == dim(site2) || error("Dimension mismatch")
+    restinds = uniqueinds(inds(t), site, site2)
+    newdata = zeros(eltype(t), dim.(restinds)..., dim(site))
+    olddata = Array(t, restinds..., site, site2)
+    for i in 1:dim(site)
+        newdata[.., i] = olddata[.., i, i]
+    end
+    return ITensor(newdata, restinds..., site)
+end

test/partitionedmps_tests.jl

@@ -30,8 +30,8 @@ import PartitionedMPSs: PartitionedMPSs, Projector, project, SubDomainMPS, Parti
         @test length([(k, v) for (k, v) in PartitionedMPS(prjΨ1)]) == 1
 
         Ψreconst = PartitionedMPS(prjΨ1) + PartitionedMPS(prjΨ2)
-        @test Ψreconst[1] ≈ prjΨ1
-        @test Ψreconst[2] ≈ prjΨ2
+        @test Ψreconst[Projector(sitesx[1] => 1)] ≈ prjΨ1
+        @test Ψreconst[Projector(sitesx[1] => 2)] ≈ prjΨ2
         @test MPS(Ψreconst) ≈ Ψ
         @test ITensors.norm(Ψreconst) ≈ ITensors.norm(MPS(Ψreconst))