Fix cutoff bug fourier transform and format

Author: gianlucagrosso

Project.toml

@@ -9,8 +9,8 @@ FastMPOContractions = "f6e391d2-8ffa-4d7a-98cd-7e70024481ca"
 ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2"
 ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 QuanticsTCI = "b11687fd-3a1c-4c41-97d0-998ab401d50e"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SparseIR = "4fe2279e-80f0-4adb-8463-ee114ff56b7d"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 TCIITensorConversion = "9f0aa9f4-9415-4e6a-8795-331ebf40aa04"

ext/QuantiscPartitionedMPSsExt/QuantiscPartitionedMPSsExt.jl

@@ -8,7 +8,8 @@ using ITensorMPS: ITensorMPS, MPS, MPO, AbstractMPS
 using ITensorMPS: findsite, linkinds, linkind, findsites
 
 import Quantics: Quantics, _find_site_allplevs, combinesites, extractdiagonal, _asdiagonal
-import PartitionedMPSs: PartitionedMPSs, SubDomainMPS, PartitionedMPS, isprojectedat, project
+import PartitionedMPSs: PartitionedMPSs, SubDomainMPS, PartitionedMPS, isprojectedat,
+                        project
 
 function Quantics.makesitediagonal(subdmps::SubDomainMPS, site::Index)
     return _makesitediagonal(subdmps, site; baseplev=0)
@@ -18,7 +19,6 @@ function Quantics.makesitediagonal(subdmps::SubDomainMPS, sites::AbstractVector{
     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)
@@ -38,7 +38,7 @@ function Quantics.makesitediagonal(subdmps::SubDomainMPS, tag::String)
 end
 
 function _makesitediagonal(
-    subdmps::SubDomainMPS, sites::AbstractVector{Index{IndsT}}; baseplev=0
+        subdmps::SubDomainMPS, sites::AbstractVector{Index{IndsT}}; baseplev=0
 ) where {IndsT}
     M_ = deepcopy(MPO(collect(MPS(subdmps))))
     for site in sites
@@ -52,9 +52,8 @@ function _makesitediagonal(subdmps::SubDomainMPS, site::Index; baseplev=0)
     return _makesitediagonal(subdmps, [site]; baseplev=baseplev)
 end
 
-
 function Quantics.extractdiagonal(
-    subdmps::SubDomainMPS, sites::AbstractVector{Index{IndsT}}
+        subdmps::SubDomainMPS, sites::AbstractVector{Index{IndsT}}
 ) where {IndsT}
     tensors = collect(subdmps.data)
     for i in eachindex(tensors)
@@ -96,20 +95,17 @@ 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)
-    ])
+    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)
-    ])
+    return PartitionedMPS([_makesitediagonal(prjmps, site; baseplev=baseplev)
+                           for prjmps in values(obj)])
 end
 
 """
@@ -120,20 +116,19 @@ 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...,
+        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")
@@ -154,9 +149,11 @@ function Quantics.automul(
     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))
+    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))
+    M2 = Quantics.rearrange_siteinds(
+        M2, combinesites(sites_M2_diag, sites_shared, sites_col))
 
     M = PartitionedMPSs.contract(M1, M2; alg=alg, kwargs...)
 
@@ -177,7 +174,8 @@ function Quantics.automul(
             end
         end
     end
-    return PartitionedMPSs.truncate(Quantics.rearrange_siteinds(M, ressites); cutoff=cutoff, maxdim=maxdim)
+    return PartitionedMPSs.truncate(
+        Quantics.rearrange_siteinds(M, ressites); cutoff=cutoff, maxdim=maxdim)
 end
 
 function _findallsiteinds_by_tag(M::PartitionedMPS; tag=tag)

src/fouriertransform.jl

@@ -28,13 +28,12 @@ function _qft(sites; cutoff::Float64=1e-14, sign::Int=1)
     R > 1 || error("The number of bits must be greater than 1")
 
     sites_MPO = collect.(zip(prime.(sites), sites))
-    fouriertt = QuanticsTCI.quanticsfouriermpo(R; sign = Float64(sign), normalize = true)
+    fouriertt = QuanticsTCI.quanticsfouriermpo(R; sign=Float64(sign), normalize=true)
     M = MPO(fouriertt; sites=sites_MPO)
 
     return truncate(M; cutoff)
 end
 
-
 abstract type AbstractFT end
 
 struct FTCore
@@ -129,7 +128,7 @@ function fouriertransform(M::MPS;
 
     # Phase shift from origindst
     M_result = phase_rotation(M, sign * 2π * origindst / (2.0^length(sitepos));
-        targetsites=target_sites)
+        targetsites=target_sites, kwargs...)
 
     # Apply QFT
     M_result = _apply(MQ, M_result; kwargs...)
@@ -141,7 +140,7 @@ function fouriertransform(M::MPS;
 
     # Phase shift from originsrc
     M_result = phase_rotation(M_result, sign * 2π * originsrc / (2.0^length(sitepos));
-        targetsites=sitesdst)
+        targetsites=sitesdst, kwargs...)
 
     M_result *= exp(sign * im * 2π * originsrc * origindst / 2.0^length(sitepos))
 

src/mul.jl

@@ -79,7 +79,8 @@ function _asdiagonal(t, site::Index{T}; baseplev=0)::ITensor where {T<:Number}
         tensor[.., i, i] = rawdata[.., i]
     end
     return ITensor(
-        tensor, links..., ITensors.prime(site, baseplev + 1), ITensors.prime(site, baseplev)
+        tensor, links..., ITensors.prime(site, baseplev + 1), ITensors.prime(
+            site, baseplev)
     )
 end
 
@@ -162,4 +163,4 @@ function automul(M1::MPS, M2::MPS; tag_row::String="", tag_shared::String="",
     end
 
     return asMPS(M)
-end
+end

src/transformer.jl

@@ -116,9 +116,9 @@ end
 """
 Multiply by exp(i θ x), where x = (x_1, ..., x_R)_2.
 """
-function phase_rotation(M::MPS, θ::Float64; targetsites=nothing, tag="")::MPS
+function phase_rotation(M::MPS, θ::Float64; targetsites=nothing, tag="", kwargs...)::MPS
     transformer = phase_rotation_mpo(siteinds(M), θ; targetsites=targetsites, tag=tag)
-    _apply(transformer, M)
+    _apply(transformer, M; kwargs...)
 end
 
 """

src/util.jl

@@ -364,11 +364,10 @@ function combinesites(M::MPO, site1::Index, site2::Index)
     return MPO(tensors)
 end
 
-
 function combinesites(
-    sites::Vector{Vector{Index{IndsT}}},
-    site1::AbstractVector{Index{IndsT}},
-    site2::AbstractVector{Index{IndsT}},
+        sites::Vector{Vector{Index{IndsT}}},
+        site1::AbstractVector{Index{IndsT}},
+        site2::AbstractVector{Index{IndsT}}
 ) where {IndsT}
     length(site1) == length(site2) || error("Length mismatch")
     for (s1, s2) in zip(site1, site2)
@@ -378,7 +377,7 @@ function combinesites(
 end
 
 function combinesites(
-    sites::Vector{Vector{Index{IndsT}}}, site1::Index, site2::Index
+        sites::Vector{Vector{Index{IndsT}}}, site1::Index, site2::Index
 ) where {IndsT}
     sites = deepcopy(sites)
     p1 = findfirst(x -> x[1] == site1, sites)
@@ -395,8 +394,6 @@ function combinesites(
     return sites
 end
 
-
-
 function directprod(::Type{T}, sites, indices) where {T}
     length(sites) == length(indices) || error("Length mismatch between sites and indices")
     any(0 .== indices) && error("indices must be 1-based")
@@ -533,18 +530,14 @@ function makesitediagonal(M::AbstractMPS, tag::String)::MPS
     return MPS(collect(M_))
 end
 
-
 # FIXME: may be type unstable
 function _find_site_allplevs(tensor::ITensor, site::Index; maxplev=10)
     ITensors.plev(site) == 0 || error("Site index must be unprimed.")
-    return [
-        ITensors.prime(site, plev) for
-        plev in 0:maxplev if ITensors.prime(site, plev) ∈ ITensors.inds(tensor)
-    ]
+    return [ITensors.prime(site, plev)
+            for
+            plev in 0:maxplev if ITensors.prime(site, plev) ∈ ITensors.inds(tensor)]
 end
 
-
-
 """
 Extract diagonal components
 """
@@ -581,7 +574,7 @@ function _apply(A::MPO, Ψ::MPO; alg::String="fit", cutoff::Real=1e-25, kwargs..
         @warn "cutoff is too small for densitymatrix algorithm. Use fit algorithm instead."
     end
     AΨ = replaceprime(
-        FastMPOContractions.contract_mpo_mpo(A', asMPO(Ψ); alg, kwargs...), 2 => 1)
+        FastMPOContractions.contract_mpo_mpo(A', asMPO(Ψ); alg, cutoff, kwargs...), 2 => 1)
     MPO(collect(AΨ))
 end
 
@@ -592,6 +585,6 @@ function _apply(A::MPO, Ψ::MPS; alg::String="fit", cutoff::Real=1e-25, kwargs..
     if alg == "densitymatrix" && cutoff <= 1e-10
         @warn "cutoff is too small for densitymatrix algorithm. Use fit algorithm instead."
     end
-    AΨ = noprime.(FastMPOContractions.contract_mpo_mpo(A, asMPO(Ψ); alg, kwargs...))
+    AΨ = noprime.(FastMPOContractions.contract_mpo_mpo(A, asMPO(Ψ); alg, cutoff, kwargs...))
     MPS(collect(AΨ))
 end

test/_util.jl

@@ -3,7 +3,7 @@ import ITensorMPS: random_mps
 using Random
 
 function _random_mpo(
-    rng::AbstractRNG, sites::AbstractVector{<:AbstractVector{Index{T}}}; linkdims::Int=1
+        rng::AbstractRNG, sites::AbstractVector{<:AbstractVector{Index{T}}}; linkdims::Int=1
 ) where {T}
     sites_ = collect(Iterators.flatten(sites))
     Ψ = random_mps(rng, sites_; linkdims)
@@ -17,7 +17,7 @@ function _random_mpo(
 end
 
 function _random_mpo(
-    sites::AbstractVector{<:AbstractVector{Index{T}}}; linkdims::Int=1
+        sites::AbstractVector{<:AbstractVector{Index{T}}}; linkdims::Int=1
 ) where {T}
     return _random_mpo(Random.default_rng(), sites; linkdims)
 end

test/fouriertransform_tests.jl

@@ -29,11 +29,11 @@
 
     @testset "qft_mpo" for sign in [1, -1], nbit in [2, 3]
         N = 2^nbit
-    
+
         sites = siteinds("Qubit", nbit)
         M = Quantics._qft(sites; sign=sign)
         M_ref = _qft_ref(sites; sign=sign)
-    
+
         @test M ≈ M_ref
     end
 end

test/mul_tests.jl

@@ -120,7 +120,8 @@ end
 
 @testitem "mul_tests.jl/batchedmatmul" begin
     using Test
-    import PartitionedMPSs: PartitionedMPSs, SubDomainMPS, PartitionedMPS, isprojectedat, project, Projector
+    import PartitionedMPSs: PartitionedMPSs, SubDomainMPS, PartitionedMPS, isprojectedat,
+                            project, Projector
     import Quantics
     using ITensors
     using ITensors.SiteTypes: siteinds

test/util_tests.jl

@@ -1,6 +1,7 @@
 @testitem "util.jl" begin
     using Test
-    import PartitionedMPSs: PartitionedMPSs, SubDomainMPS, PartitionedMPS, project, isprojectedat
+    import PartitionedMPSs: PartitionedMPSs, SubDomainMPS, PartitionedMPS, project,
+                            isprojectedat
     import Quantics
     using ITensors
     using ITensors.SiteTypes: siteinds