Merge pull request #25 from tensor4all/6-use-qft-mpo-from-tensorcross…

…interpolationjl

Use QuanticsTCI to generate a QFT MPO

Project.toml

@@ -10,18 +10,22 @@ 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"
 
 [compat]
 EllipsisNotation = "1"
 FastMPOContractions = "^0.2.2"
 ITensorMPS = "0.3.1"
 ITensors = "0.7"
+ProjMPSs = "0.4.1"
+QuanticsTCI = "0.7.0"
 SparseIR = "^0.96, 0.97, 1"
 StaticArrays = "1"
+TCIITensorConversion = "0.2.0"
 julia = "1"
-ProjMPSs = "0.4.1"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"

src/Quantics.jl

@@ -4,14 +4,16 @@ using ITensors
 import ITensors
 using ITensors.SiteTypes: siteinds
 import ITensors.NDTensors: Tensor, BlockSparseTensor, blockview
-using ITensorMPS: MPS, MPO, AbstractMPS
+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
 using StaticArrays
 
+import QuanticsTCI
+import TCIITensorConversion
 import FastMPOContractions
 
 using EllipsisNotation

src/fouriertransform.jl

@@ -20,124 +20,20 @@ We denote the input and output MPS's by ``X`` and ``Y``, respectively.
 
 """
 function _qft(sites; cutoff::Float64=1e-14, sign::Int=1)
-    if any([!hastags(inds(s), "Qubit") for s in sites])
+    if !all(dim.(sites) .== 2)
         error("All siteinds for qft must has Qubit tag")
     end
-    M = _qft_wo_norm(sites; cutoff=cutoff, sign=sign)
-    M *= 2.0^(-0.5 * length(sites))
 
-    # Quick hack: In the Markus's note,
-    # the digits are ordered oppositely from the present convention.
-    M = MPO([M[n] for n in length(M):-1:1])
-    _replace_mpo_siteinds!(M, reverse(sites), sites)
+    R = length(sites)
+    R > 1 || error("The number of bits must be greater than 1")
 
-    return M
-end
-
-function _assign!(M::MPO, n::Int, arr; autoreshape=false)
-    if autoreshape
-        arr = reshape(arr, map(dim, inds(M[n]))...)
-    end
-    M[n] = ITensor(arr, inds(M[n])...)
-    return nothing
-end
-
-"""
-For length(sites) == 1
-The resultant MPO is NOT renormalized.
-"""
-function _qft_nsite1_wo_norm(sites; sign::Int=1)
-    length(sites) == 1 || error("num sites > 1")
-    _exp(x, k) = exp(sign * im * π * (x - 1) * (k - 1))
-
-    arr = zeros(ComplexF64, 2, 2)
-    for out in 1:2, in in 1:2
-        arr[out, in] = _exp(out, in)
-    end
-
-    M = Quantics._zero_mpo(sites)
-    _assign!(M, 1, arr)
-
-    return M
-end
-
-function _qft_wo_norm(sites; cutoff::Float64=1e-14, sign::Int=1)
-    N = length(sites)
-    if N == 1
-        return _qft_nsite1_wo_norm(sites; sign=sign)
-    end
-
-    M_prev = _qft_wo_norm(sites[2:end]; cutoff=cutoff, sign=sign)
-    M_top = _qft_toplayer(sites; sign=sign)
+    sites_MPO = collect.(zip(prime.(sites), sites))
+    fouriertt = QuanticsTCI.quanticsfouriermpo(R; sign = Float64(sign), normalize = true)
+    M = MPO(fouriertt; sites=sites_MPO)
 
-    M = _contract(M_top, M_prev)
-    ITensors.truncate!(M; cutoff=cutoff)
-
-    return M
+    return truncate(M; cutoff)
 end
 
-function _qft_toplayer(sites; sign::Int=1)
-    N = length(sites)
-    N > 1 || error("N must be greater than 1")
-
-    tensors = []
-
-    # site = 1
-    arr = zeros(ComplexF64, 2, 2, 2)
-    for x in 1:2, k in 1:2
-        # arr: (out, in, link)
-        arr[x, k, k] = exp(sign * im * π * (x - 1) * (k - 1))
-    end
-    push!(tensors, arr)
-
-    for n in 2:N
-        ϕ = π * 0.5^(n - 1)
-        _exp(x, k) = exp(sign * im * ϕ * (x - 1) * (k - 1))
-        # Right most tensor
-        if n == N
-            # arr: (link, out, in)
-            arr = zeros(ComplexF64, 2, 2, 2)
-            for x in 1:2, k in 1:2
-                arr[k, x, x] = _exp(x, k)
-            end
-            push!(tensors, arr)
-        else
-            # arr: (link_left, out, in, link_right)
-            arr = zeros(ComplexF64, 2, 2, 2, 2)
-            for x in 1:2, k in 1:2
-                arr[k, x, x, k] = _exp(x, k)
-            end
-            push!(tensors, arr)
-        end
-    end
-
-    M = Quantics._zero_mpo(sites; linkdims=fill(2, N - 1))
-    for n in 1:N
-        _assign!(M, n, tensors[n])
-    end
-
-    return M
-end
-
-function _contract(M_top, M_prev)
-    length(M_top) == length(M_prev) + 1 || error("Length mismatch")
-    N = length(M_top)
-    M_top = ITensors.replaceprime(M_top, 1 => 2; tags="Qubit")
-    M_top = ITensors.replaceprime(M_top, 0 => 1; tags="Qubit")
-    M_top_ = ITensors.data(M_top)
-    M_prev_ = ITensors.data(M_prev)
-
-    M_data = [M_top_[1]]
-    for n in 1:(N - 1)
-        push!(M_data, M_top_[n + 1] * M_prev_[n])
-    end
-
-    M = MPO(M_data)
-    M = ITensors.replaceprime(M, 1 => 0; tags="Qubit")
-    M = ITensors.replaceprime(M, 2 => 1; tags="Qubit")
-
-    return M
-end
 
 abstract type AbstractFT end
 

test/fouriertransform_tests.jl

@@ -26,15 +26,15 @@
         M = MPO(trans_t, sites; cutoff=cutoff)
         return M
     end
-    @testset "qft_mpo" for sign in [1, -1], nbit in [1, 2, 3]
-        N = 2^nbit
 
+    @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 Array(reduce(*, M), vcat(sites, sites')) ≈
-              Array(reduce(*, M_ref), vcat(sites, sites'))
+
+        @test M ≈ M_ref
     end
 end