move some utilities in YaoExtensions to Yao.EasyBuild submodule. (#315)

* integrate utilities in YaoExtensions.

* not exporting EasyBuild

* clean up EasyBuild

* fsim_circuit

* rm KMod and mathgate, add hadamard test

* add phase estimation

* rename

* clean up

* add docstrings

src/EasyBuild/block_extension/FSimGate.jl

@@ -0,0 +1,51 @@
+export FSimGate, fsim_block
+
+"""
+    FSimGate{T<:Number} <: PrimitiveBlock{2}
+
+The two parameter `FSim` gate.
+
+References
+-------------------------
+* Arute, Frank, et al. "Quantum supremacy using a programmable superconducting processor." Nature 574.7779 (2019): 505-510.
+"""
+mutable struct FSimGate{T<:Number} <: PrimitiveBlock{2}
+    theta::T
+    phi::T
+end
+
+function Base.:(==)(fs1::FSimGate, fs2::FSimGate)
+    return fs1.theta == fs2.theta && fs1.phi == fs2.phi
+end
+
+function YaoAPI.mat(::Type{T}, fs::FSimGate) where T
+    θ, ϕ = fs.theta, fs.phi
+    T[1 0          0          0;
+     0 cos(θ)     -im*sin(θ) 0;
+     0 -im*sin(θ) cos(θ)     0;
+     0 0          0          exp(-im*ϕ)]
+end
+
+YaoAPI.iparams_eltype(::FSimGate{T}) where T = T
+YaoAPI.getiparams(fs::FSimGate{T}) where T = (fs.theta, fs.phi)
+function YaoAPI.setiparams!(fs::FSimGate{T}, θ, ϕ) where T
+    fs.theta = θ
+    fs.phi = ϕ
+    return fs
+end
+
+YaoBlocks.@dumpload_fallback FSimGate FSimGate
+YaoBlocks.Optimise.to_basictypes(fs::FSimGate) = fsim_block(fs.theta, fs.phi)
+
+"""
+    fsim_block(θ::Real, ϕ::Real)
+
+The circuit representation of FSim gate.
+"""
+function fsim_block(θ::Real, ϕ::Real)
+    if θ ≈ π/2
+        return cphase(2,2,1,-ϕ)*SWAP*rot(kron(Z,Z), -π/2)*put(2,1=>phase(-π/4))
+    else
+        return cphase(2,2,1,-ϕ)*rot(SWAP,2*θ)*rot(kron(Z,Z), -θ)*put(2,1=>phase(θ/2))
+    end
+end

src/EasyBuild/block_extension/blocks.jl

@@ -0,0 +1,3 @@
+import .YaoBlocks: _apply!
+include("shortcuts.jl")
+include("FSimGate.jl")

src/EasyBuild/block_extension/shortcuts.jl

@@ -0,0 +1,17 @@
+export ISWAP, SqrtX, SqrtY, SqrtW, singlet_block
+export ISWAPGate, SqrtXGate, SqrtYGate, SqrtWGate, CPhaseGate
+
+const CPhaseGate{N, T} = ControlBlock{N,<:ShiftGate{T},<:Any}
+
+@const_gate ISWAP = PermMatrix([1,3,2,4], [1,1.0im,1.0im,1])
+@const_gate SqrtX = [0.5+0.5im 0.5-0.5im; 0.5-0.5im 0.5+0.5im]
+@const_gate SqrtY = [0.5+0.5im -0.5-0.5im; 0.5+0.5im 0.5+0.5im]
+# √W is a non-Clifford gate
+@const_gate SqrtW = mat(rot((X+Y)/sqrt(2), π/2))
+
+"""
+    singlet_block(θ::Real, ϕ::Real)
+
+The circuit block for initialzing a singlet state.
+"""
+singlet_block() = chain(put(2, 1=>chain(X, H)), control(2, -1, 2=>X))

src/EasyBuild/easybuild.jl

@@ -0,0 +1,13 @@
+module EasyBuild
+using YaoBlocks, YaoBlocks.LuxurySparse, YaoBlocks.YaoAPI, YaoBlocks.YaoArrayRegister
+using YaoBlocks.LinearAlgebra
+include("block_extension/blocks.jl")
+include("general_U4.jl")
+include("phaseestimation.jl")
+include("qft_circuit.jl")
+include("hamiltonians.jl")
+include("variational_circuit.jl")
+include("supremacy_circuit.jl")
+include("google53.jl")
+include("hadamardtest.jl")
+end

src/EasyBuild/general_U4.jl

@@ -0,0 +1,47 @@
+"""
+Impliments PRA 69.062321.
+"""
+
+export general_U2, general_U4
+
+"""
+    general_U2(θ1, θ2, θ3; ϕ=nothing)
+
+A general single qubits gate: ``e^(iϕ)R_z(θ_3)R_y(θ_2)R_z(θ_1)``.
+Leave `ϕ` as `nothing` to fix the global phase.
+"""
+function general_U2(θ1, θ2, θ3; ϕ=nothing)
+    gate = Rz(θ3) * Ry(θ2) * Rz(θ1)
+    if ϕ !== nothing
+        push!(gate, phase(ϕ))
+    end
+    return gate
+end
+
+"""
+    general_U4([params...]) -> AbstractBlock
+
+A general two qubits gate decomposed to (CNOT, Ry, Rz), parameters default to 0.
+
+!!!note
+
+    Although the name is U(4), This is actually a SU(4) gate up to a phase, the phase `det(dispatch!(general_U4(), :random))` is fixed to -1.
+"""
+general_U4() = general_U4(zeros(15))
+function general_U4(params)
+    if length(params) != 15
+        throw(ArgumentError("The number of parameters must be 15, got $(length(params))"))
+    end
+    return chain(2, [
+        put(2, 1=>general_U2(params[1:3]...)),
+        put(2, 2=>general_U2(params[4:6]...)),
+        cnot(2, 2, 1),
+        put(2, 1=>Rz(params[7])),
+        put(2, 2=>Ry(params[8])),
+        cnot(2, 1, 2),
+        put(2, 2=>Ry(params[9])),
+        cnot(2, 2, 1),
+        put(2, 1=>general_U2(params[10:12]...)),
+        put(2, 2=>general_U2(params[13:15]...))
+    ])
+end

src/EasyBuild/google53.jl

@@ -0,0 +1,104 @@
+export Lattice53, rand_google53
+
+entangler_google53(nbits::Int, i::Int, j::Int) = put(nbits, (i,j)=>FSimGate(π/2, π/6))
+
+struct Lattice53
+    labels::Matrix{Int}
+end
+
+function Lattice53(;nbits::Int=53)
+    config = ones(Bool, 5, 12)
+    config[end,2:2:end] .= false
+    config[1, 7] = false
+    labels = zeros(Int, 5, 12)
+    k = 0
+    for (i,c) in enumerate(config)
+        if c
+            k += 1
+            labels[i] = k
+            k>=nbits && break
+        end
+    end
+    return Lattice53(labels)
+end
+
+nbits(lattice::Lattice53) = maximum(lattice.labels)
+
+function Base.getindex(lattice::Lattice53, i, j)
+    1<=i<=size(lattice.labels, 1) && 1<=j<=size(lattice.labels, 2) ? lattice.labels[i,j] : 0
+end
+upperleft(lattice::Lattice53,i,j) = lattice[i-j%2,j-1]
+lowerleft(lattice::Lattice53,i,j) = lattice[i+(j-1)%2,j-1]
+upperright(lattice::Lattice53,i,j) = lattice[i-j%2,j+1]
+lowerright(lattice::Lattice53,i,j) = lattice[i+(j-1)%2,j+1]
+
+function pattern53(lattice::Lattice53, chr::Char)
+    res = Tuple{Int,Int}[]
+    # i0, di, j0, dj and direction
+    di = 1 + (chr>'D')
+    dj = 2 - (chr>'D')
+    j0 = 1 + min(dj-1, mod(chr-'A',2))
+    direction = 'C'<=chr<='F' ? lowerright : upperright
+    for j=j0:dj:12
+        i0 = chr>'D' ? mod((chr-'D') + (j-(chr>='G'))÷2, 2) : 1
+        for i = i0:di:5
+            src = lattice[i, j]
+            dest = direction(lattice, i, j)
+            src!=0 && dest !=0 && push!(res, (src, dest))
+        end
+    end
+    return res
+end
+
+function print_lattice53(lattice, pattern)
+    for i_=1:10
+        i = (i_+1)÷2
+        for j=1:12
+            if i_%2 == j%2 && lattice[i,j]!=0
+                print(" ∘  ")
+            else
+                print("    ")
+            end
+        end
+        println()
+        for j=1:12
+            if i_%2 == j%2 && lattice[i,j]!=0
+                hasll = (lowerleft(lattice, i, j), lattice[i,j]) in pattern
+                haslr = (lattice[i,j], lowerright(lattice, i, j)) in pattern
+                print(hasll ? "/ " : "  ")
+                print(haslr ? " \\" : "  ")
+            else
+                print("    ")
+            end
+        end
+        println()
+    end
+end
+
+"""
+    rand_google53(depth::Int; nbits=53) -> AbstactBlock
+
+Google supremacy circuit with 53 qubits, also know as the Sycamore quantum supremacy circuits.
+
+References
+-------------------------
+* Arute, Frank, et al. "Quantum supremacy using a programmable superconducting processor." Nature 574.7779 (2019): 505-510.
+"""
+function rand_google53(depth::Int; nbits::Int=53)
+    c = chain(nbits)
+    lattice = Lattice53(nbits=nbits)
+    k = 0
+    for pattern in Iterators.cycle(['A', 'B', 'C', 'D', 'C', 'D', 'A', 'B'])
+        push!(c, rand_google53_layer(lattice, pattern))
+        k += 1
+        k>=depth && break
+    end
+    return c
+end
+
+function rand_google53_layer(lattice, pattern)
+    nbit = nbits(lattice)
+    chain(nbit, chain(nbit, [put(nbit, i=>rand([SqrtW, SqrtX, SqrtY])) for i=1:nbit]),
+        chain(nbit, [entangler_google53(nbit,i,j) for (i,j) in pattern53(lattice, pattern)])
+        )
+end

src/EasyBuild/hadamardtest.jl

@@ -0,0 +1,44 @@
+export hadamard_test, hadamard_test_circuit, swap_test_circuit
+
+"""
+    hadamard_test_circuit(U::AbstractBlock, ϕ::Real)
+
+The Hadamard test circuit.
+
+References
+-----------------------
+* [Wiki](https://en.wikipedia.org/wiki/Hadamard_test_(quantum_computation))
+"""
+function hadamard_test_circuit(U::AbstractBlock{N}, ϕ::Real) where N
+    chain(N+1, put(N+1, 1=>H),
+        put(N+1, 1=>Rz(ϕ)),
+        control(N+1, 1, 2:N+1=>U),  # get matrix first, very inefficient
+        put(N+1, 1=>H)
+        )
+end
+
+function hadamard_test(U::AbstractBlock{N}, reg::AbstractRegister, ϕ::Real) where N
+    c = hadamard_test_circuit(U, ϕ::Real)
+    reg = join(reg, zero_state(1))
+    expect(put(N+1, 1=>Z), reg |> c)
+end
+
+"""
+    swap_test_circuit(nbit::Int, nstate::Int, ϕ::Real)
+
+The swap test circuit for computing the overlap between multiple density matrices.
+The `nbit` and `nstate` specifies the number of qubit in each state and how many state we want to compare.
+
+References
+-----------------------
+* Ekert, Artur K., et al. "Direct estimations of linear and nonlinear functionals of a quantum state." Physical review letters 88.21 (2002): 217901.
+"""
+function swap_test_circuit(nbit::Int, nstate::Int, ϕ::Real)
+    N = nstate*nbit + 1
+
+    chain(N, put(N, 1=>H),
+        put(N, 1=>Rz(ϕ)),
+        chain(N, [chain(N, [control(N, 1, (i+(k*nbit-nbit)+1, i+k*nbit+1)=>SWAP) for i=1:nbit]) for k=1:nstate-1]),  # get matrix first, very inefficient
+        put(N, 1=>H)
+        )
+end

src/EasyBuild/hamiltonians.jl

@@ -0,0 +1,31 @@
+export heisenberg, transverse_ising
+
+"""
+    heisenberg(nbit::Int; periodic::Bool=true)
+
+1D Heisenberg hamiltonian defined as ``\\sum_{i=1}^{n} X_{i}X_{i+1} + Y_{i}Y_{i+1} + Z_{i}Z_{i+1}``, where ``n`` is specified by `nbit`.
+`periodic` means the boundary condition is periodic.
+
+References
+----------------------
+* de Oliveira, Mário J. "Ground-state properties of the spin-1/2 antiferromagnetic Heisenberg chain obtained by use of a Monte Carlo method." Physical Review B 48.9 (1993): 6141-6143.
+"""
+function heisenberg(nbit::Int; periodic::Bool=true)
+    map(1:(periodic ? nbit : nbit-1)) do i
+        j=i%nbit+1
+        repeat(nbit,X,(i,j)) + repeat(nbit, Y, (i,j)) + repeat(nbit, Z, (i,j))
+    end |> sum
+end
+
+"""
+    transverse_ising(nbit::Int, h::Number; periodic::Bool=true)
+
+1D transverse Ising hamiltonian defined as ``\\sum_{i=1}^{n} hX_{i} + Z_{i}Z_{i+1}``, where ``n`` is specified by `nbit`.
+`periodic` means the boundary condition is periodic.
+"""
+function transverse_ising(nbit::Int, h::Number; periodic::Bool=true)
+    ising_term = map(1:(periodic ? nbit : nbit-1)) do i
+        repeat(nbit,Z,(i,i%nbit+1))
+    end |> sum
+    ising_term + h*sum(map(i->put(nbit,i=>X), 1:nbit))
+end

src/EasyBuild/phaseestimation.jl

@@ -0,0 +1,53 @@
+export phase_estimation_circuit, phase_estimation_analysis
+
+"""
+    phase_estimation_circuit(unitarygate::GeneralMatrixBlock, n_reg, n_b) -> ChainBlock
+
+Phase estimation circuit. Input arguments are
+
+* `unitarygate`: the input unitary matrix.
+* `n_reg`: the number of bits to store phases,
+* `n_b`: the number of bits to store vector.
+
+References
+----------------------
+[Wiki](https://en.wikipedia.org/wiki/Quantum_phase_estimation_algorithm)
+"""
+function phase_estimation_circuit(unitarygate::GeneralMatrixBlock, n_reg::Int, n_b::Int)
+    nbit = n_b + n_reg
+    # Apply Hadamard Gate.
+    hs = repeat(nbit, H, 1:n_reg)
+
+    # Construct a control circuit.
+    control_circuit = chain(nbit)
+    for i = 1:n_reg
+        push!(control_circuit, control(nbit, (i,), (n_reg+1:nbit...,)=>unitarygate))
+        if i != n_reg
+            unitarygate = matblock(mat(unitarygate) * mat(unitarygate))
+        end
+    end
+
+    # Inverse QFT Block.
+    iqft = subroutine(nbit, qft_circuit(n_reg)',[1:n_reg...,])
+    chain(hs, control_circuit, iqft)
+end
+
+"""
+    phase_estimation_analysis(eigenvectors::Matrix, reg::ArrayReg) -> Tuple
+
+Analyse phase estimation result using state projection.
+It returns a tuple of (most probable configuration, the overlap matrix, the relative probability for this configuration)
+`eigenvectors` is the eigen vectors of the unitary gate matrix, while `reg` is the result of phase estimation.
+"""
+function phase_estimation_analysis(eigenvectors::AbstractMatrix, reg::ArrayReg)
+    overlap = eigenvectors'*state(reg)
+    amp_relative = Float64[]
+    bs = Int[]
+
+    for b in basis(overlap)
+        mc = argmax(view(overlap, b+1, :) .|> abs)-1
+        push!(amp_relative, abs2(overlap[b+1, mc+1])/sum(overlap[b+1, :] .|> abs2))
+        push!(bs, mc)
+    end
+    bs, overlap, amp_relative
+end

src/EasyBuild/qft_circuit.jl

@@ -0,0 +1,20 @@
+export qft_circuit
+
+"""
+    cphase(nbits, i, j, θ)
+
+Control-phase gate.
+"""
+cphase(nbits::Int, i::Int, j::Int, θ::T) where T = control(nbits, i, j=>shift(θ))
+
+"""
+    qft_circuit(n)
+
+The quantum Fourer transformation (QFT) circuit.
+
+References
+------------------------
+* [Wiki](https://en.wikipedia.org/wiki/Quantum_Fourier_transform)
+"""
+qft_circuit(n::Int) = chain(n, hcphases(n, i) for i = 1:n)
+hcphases(n, i) = chain(n, i==j ? put(i=>H) : cphase(n, j, i, 2π/(2^(j-i+1))) for j in i:n);