diff --git a/Project.toml b/Project.toml index dca3147..6bbf172 100644 --- a/Project.toml +++ b/Project.toml @@ -4,15 +4,20 @@ authors = ["Roger-luo "] version = "0.1.0" [deps] +BitBasis = "50ba71b6-fa0f-514d-ae9a-0916efc90dcf" DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LuxurySparse = "d05aeea4-b7d4-55ac-b691-9e7fabb07ba2" +MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" Yao = "5872b779-8223-5990-8dd0-5abbb0748c8c" +YaoArrayRegister = "e600142f-9330-5003-8abb-0ebd767abc51" +YaoBlocks = "418bc28f-b43b-5e0b-a6e7-61bbc1a2c1df" [extras] +MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/README.md b/README.md index 37a461a..dbf8a3e 100644 --- a/README.md +++ b/README.md @@ -26,8 +26,9 @@ Please notice, this package is still under development and needs further polish. - [x] QuGAN - [x] QCBM - [x] Hamiltonian Solver -- [ ] QAOA -- [ ] Quantum Chemistry +- [x] QAOA +- [x] Quantum Chemistry +- [x] QuODE ## License diff --git a/examples/Grover.jl b/examples/Grover.jl index 48d2006..891a560 100644 --- a/examples/Grover.jl +++ b/examples/Grover.jl @@ -5,7 +5,7 @@ using QuAlgorithmZoo: groveriter, inference_oracle, prob_match_oracle # ## Target Space and Evidense num_bit = 12 -oracle = matrixgate(Diagonal((v = ones(1<X) for i=1:nbit]) tb = TimeEvolution(HB(3), 0.1) function HC(W::AbstractMatrix) nbit = size(W, 1) - ab = add(nbit) + ab = Any[] for i=1:nbit,j=i+1:nbit if W[i,j] != 0 push!(ab, 0.5*W[i,j]*repeat(nbit, Z, [i,j])) end end - ab + sum(ab) end function qaoa_circuit(W::AbstractMatrix, depth::Int; use_cache::Bool=false) @@ -22,7 +22,7 @@ function qaoa_circuit(W::AbstractMatrix, depth::Int; use_cache::Bool=false) hc = HC(W) use_cache && (hb = hb |> cache; hc = hc |> cache) c = chain(nbit, [repeat(nbit, H, 1:nbit)]) - append!(c, [chain(nbit, [timeevolve(hc, 0.0, tol=1e-5), timeevolve(hb, 0.0, tol=1e-5)]) for i=1:depth]) + append!(c, [chain(nbit, [time_evolve(hc, 0.0, tol=1e-5), time_evolve(hb, 0.0, tol=1e-5)]) for i=1:depth]) end diff --git a/examples/VQE.jl b/examples/VQE.jl index ede7561..d45986e 100644 --- a/examples/VQE.jl +++ b/examples/VQE.jl @@ -1,17 +1,17 @@ -using Yao, Yao.Blocks +using Yao using QuAlgorithmZoo using KrylovKit function ed_groundstate(h::MatrixBlock) E, V = eigsolve(h |> mat, 1, :SR, ishermitian=true) println("Ground State Energy is $(E[1])") - register(V[1]) + ArrayReg(V[1]) end N = 5 c = random_diff_circuit(N, N, [i=>mod(i,N)+1 for i=1:N], mode=:Merged) |> autodiff(:QC) dispatch!(c, :random) hami = heisenberg(N) -ed_groundstate(hami) -vqe_solve(c, hami) +# vqe ground state +vqe_solve!(c, hami) diff --git a/examples/make.jl b/examples/make.jl deleted file mode 100644 index 0daf18e..0000000 --- a/examples/make.jl +++ /dev/null @@ -1,3 +0,0 @@ -using Literate - -Literate.notebook("QCBM.jl", joinpath(@__DIR__, "../notebooks"), execute=false) diff --git a/src/Diff.jl b/src/Diff.jl new file mode 100644 index 0000000..94b63d9 --- /dev/null +++ b/src/Diff.jl @@ -0,0 +1,244 @@ +export Rotor, generator, AbstractDiff, BPDiff, QDiff, backward!, gradient, CPhaseGate, DiffBlock +import Yao: expect, content, chcontent +using StatsBase + +############# General Rotor ############ +const Rotor{N, T} = Union{RotationGate{N, T}, PutBlock{N, <:Any, <:RotationGate, <:Complex{T}}} +const CphaseGate{N, T} = ControlBlock{N,<:ShiftGate{T},<:Any} +const DiffBlock{N, T} = Union{Rotor{N, T}, CphaseGate{N, T}} +""" + generator(rot::Rotor) -> MatrixBlock + +Return the generator of rotation block. +""" +generator(rot::RotationGate) = rot.block +generator(rot::PutBlock{N, C, GT}) where {N, C, GT<:RotationGate} = PutBlock{N}(generator(rot|>content), rot |> occupied_locs) +generator(c::CphaseGate{N}) where N = ControlBlock(N, c.ctrol_locs, ctrl_config, control(2,1,2=>Z), c.locs) + +abstract type AbstractDiff{GT, N, T} <: TagBlock{GT, N, T} end +Base.adjoint(df::AbstractDiff) = Daggered(df) + +istraitkeeper(::AbstractDiff) = Val(true) + +#################### The Basic Diff ################# +""" + QDiff{GT, N, T} <: AbstractDiff{GT, N, Complex{T}} + QDiff(block) -> QDiff + +Mark a block as quantum differentiable. +""" +mutable struct QDiff{GT, N, T} <: AbstractDiff{GT, N, Complex{T}} + block::GT + grad::T + QDiff(block::DiffBlock{N, T}) where {N, T} = new{typeof(block), N, T}(block, T(0)) +end +content(cb::QDiff) = cb.block +chcontent(cb::QDiff, blk::DiffBlock) = QDiff(blk) + +@forward QDiff.block mat, apply! +Base.adjoint(df::QDiff) = QDiff(content(df)') + +function YaoBlocks.print_annotation(io::IO, df::QDiff) + printstyled(io, "[̂∂] "; bold=true, color=:yellow) +end + +#################### The Back Propagation Diff ################# +""" + BPDiff{GT, N, T, PT} <: AbstractDiff{GT, N, Complex{T}} + BPDiff(block, [grad]) -> BPDiff + +Mark a block as differentiable, here `GT`, `PT` is gate type, parameter type. + +Warning: + please don't use the `adjoint` after `BPDiff`! `adjoint` is reserved for special purpose! (back propagation) +""" +mutable struct BPDiff{GT, N, T, PT} <: AbstractDiff{GT, N, T} + block::GT + grad::PT + input::AbstractRegister + BPDiff(block::MatrixBlock{N, T}, grad::PT) where {N, T, PT} = new{typeof(block), N, T, typeof(grad)}(block, grad) +end +BPDiff(block::MatrixBlock) = BPDiff(block, zeros(parameters_eltype(block), nparameters(block))) +BPDiff(block::DiffBlock{N, T}) where {N, T} = BPDiff(block, T(0)) + +content(cb::BPDiff) = cb.block +chcontent(cb::BPDiff, blk::MatrixBlock) = BPDiff(blk) + +@forward BPDiff.block mat +function apply!(reg::AbstractRegister, df::BPDiff) + if isdefined(df, :input) + copyto!(df.input, reg) + else + df.input = copy(reg) + end + apply!(reg, content(df)) + reg +end + +function apply!(δ::AbstractRegister, adf::Daggered{<:BPDiff{<:Rotor}}) + df = adf |> content + apply!(δ, content(df)') + df.grad = -statevec(df.input |> generator(content(df)))' * statevec(δ) |> imag + δ +end + +function YaoBlocks.print_annotation(io::IO, df::BPDiff) + printstyled(io, "[∂] "; bold=true, color=:yellow) +end + + +#### interface ##### +export autodiff, numdiff, opdiff, StatFunctional, statdiff, as_weights + +as_weights(probs::AbstractVector{T}) where T = Weights(probs, T(1)) +""" + autodiff(mode::Symbol, block::AbstractBlock) -> AbstractBlock + autodiff(mode::Symbol) -> Function + +automatically mark differentiable items in a block tree as differentiable. +""" +function autodiff end +autodiff(mode::Symbol) = block->autodiff(mode, block) +autodiff(mode::Symbol, block::AbstractBlock) = autodiff(Val(mode), block) + +# for BP +autodiff(::Val{:BP}, block::DiffBlock) = BPDiff(block) +autodiff(::Val{:BP}, block::AbstractBlock) = block +# Sequential, Roller and ChainBlock can propagate. +function autodiff(mode::Val{:BP}, blk::Union{ChainBlock, Roller, Sequential}) + chsubblocks(blk, autodiff.(mode, subblocks(blk))) +end + +# for QC +autodiff(::Val{:QC}, block::Union{RotationGate, CphaseGate}) = QDiff(block) +# escape control blocks. +autodiff(::Val{:QC}, block::ControlBlock) = block + +function autodiff(mode::Val{:QC}, blk::AbstractBlock) + blks = subblocks(blk) + isempty(blks) ? blk : chsubblocks(blk, autodiff.(mode, blks)) + end + +@inline function _perturb(func, gate::AbstractDiff{<:DiffBlock}, δ::Real) + dispatch!(-, gate, (δ,)) + r1 = func() + dispatch!(+, gate, (2δ,)) + r2 = func() + dispatch!(-, gate, (δ,)) + r1, r2 +end + +@inline function _perturb(func, gate::AbstractDiff{<:Rotor}, δ::Real) # for put + dispatch!(-, gate, (δ,)) + r1 = func() + dispatch!(+, gate, (2δ,)) + r2 = func() + dispatch!(-, gate, (δ,)) + r1, r2 +end + +""" + numdiff(loss, diffblock::AbstractDiff; δ::Real=1e-2) + +Numeric differentiation. +""" +@inline function numdiff(loss, diffblock::AbstractDiff; δ::Real=1e-2) + r1, r2 = _perturb(loss, diffblock, δ) + diffblock.grad = (r2 - r1)/2δ +end + +""" + opdiff(psifunc, diffblock::AbstractDiff, op::MatrixBlock) + +Operator differentiation. +""" +@inline function opdiff(psifunc, diffblock::AbstractDiff, op::MatrixBlock) + r1, r2 = _perturb(()->expect(op, psifunc()) |> real, diffblock, π/2) + diffblock.grad = (r2 - r1)/2 +end + +""" + StatFunctional{N, AT} + StatFunctional(array::AT<:Array) -> StatFunctional{N, <:Array} + StatFunctional{N}(func::AT<:Function) -> StatFunctional{N, <:Function} + +statistic functional, i.e. + * if `AT` is an array, A[i,j,k...], it is defined on finite Hilbert space, which is `∫A[i,j,k...]p[i]p[j]p[k]...` + * if `AT` is a function, F(xᵢ,xⱼ,xₖ...), this functional is `1/C(r,n)... ∑ᵢⱼₖ...F(xᵢ,xⱼ,xₖ...)`, see U-statistics for detail. + +References: + U-statistics, http://personal.psu.edu/drh20/asymp/fall2006/lectures/ANGELchpt10.pdf +""" +struct StatFunctional{N, AT} + data::AT + StatFunctional{N}(data::AT) where {N, AT<:Function} = new{N, AT}(data) + StatFunctional(data::AT) where {N, AT<:AbstractArray{<:Real, N}} = new{N, AT}(data) +end + +@forward StatFunctional.data Base.ndims +Base.parent(stat::StatFunctional) = stat.data + +expect(stat::StatFunctional{2, <:AbstractArray}, px::Weights, py::Weights=px) = px.values' * stat.data * py.values +expect(stat::StatFunctional{1, <:AbstractArray}, px::Weights) = stat.data' * px.values +function expect(stat::StatFunctional{2, <:Function}, xs::AbstractVector{T}) where T + N = length(xs) + res = zero(stat.data(xs[1], xs[1])) + for i = 2:N + for j = 1:i-1 + @inbounds res += stat.data(xs[i], xs[j]) + end + end + res/binomial(N,2) +end +function expect(stat::StatFunctional{2, <:Function}, xs::AbstractVector, ys::AbstractVector) + M = length(xs) + N = length(ys) + ci = CartesianIndices((M, N)) + @inbounds mapreduce(ind->stat.data(xs[ind[1]], ys[ind[2]]), +, ci)/M/N +end +expect(stat::StatFunctional{1, <:Function}, xs::AbstractVector) = mean(stat.data.(xs)) +Base.ndims(stat::StatFunctional{N}) where N = N + +""" + statdiff(probfunc, diffblock::AbstractDiff, stat::StatFunctional{<:Any, <:AbstractArray}; initial::AbstractVector=probfunc()) + statdiff(samplefunc, diffblock::AbstractDiff, stat::StatFunctional{<:Any, <:Function}; initial::AbstractVector=samplefunc()) + +Differentiation for statistic functionals. +""" +@inline function statdiff(probfunc, diffblock::AbstractDiff, stat::StatFunctional{2}; initial::AbstractVector=probfunc()) + r1, r2 = _perturb(()->expect(stat, probfunc(), initial), diffblock, π/2) + diffblock.grad = (r2 - r1)*ndims(stat)/2 +end +@inline function statdiff(probfunc, diffblock::AbstractDiff, stat::StatFunctional{1}) + r1, r2 = _perturb(()->expect(stat, probfunc()), diffblock, π/2) + diffblock.grad = (r2 - r1)*ndims(stat)/2 +end + +""" + backward!(δ::AbstractRegister, circuit::MatrixBlock) -> AbstractRegister + +back propagate and calculate the gradient ∂f/∂θ = 2*Re(∂f/∂ψ*⋅∂ψ*/∂θ), given ∂f/∂ψ*. + +Note: +Here, the input circuit should be a matrix block, otherwise the back propagate may not apply (like Measure operations). +""" +backward!(δ::AbstractRegister, circuit::MatrixBlock) = apply!(δ, circuit') + +""" + gradient(circuit::AbstractBlock, mode::Symbol=:ANY) -> Vector + +collect all gradients in a circuit, mode can be :BP/:QC/:ANY, they will collect `grad` from BPDiff/QDiff/AbstractDiff respectively. +""" +gradient(circuit::AbstractBlock, mode::Symbol=:ANY) = gradient!(circuit, parameters_eltype(circuit)[], mode) + +gradient!(circuit::AbstractBlock, grad, mode::Symbol) = gradient!(circuit, grad, Val(mode)) +function gradient!(circuit::AbstractBlock, grad, mode::Val) + for block in subblocks(circuit) + gradient!(block, grad, mode) + end + grad +end + +gradient!(circuit::BPDiff, grad, mode::Val{:BP}) = append!(grad, circuit.grad) +gradient!(circuit::QDiff, grad, mode::Val{:QC}) = push!(grad, circuit.grad) +gradient!(circuit::AbstractDiff, grad, mode::Val{:ANY}) = append!(grad, circuit.grad) diff --git a/src/Grover.jl b/src/Grover.jl index 3207680..ed369a8 100644 --- a/src/Grover.jl +++ b/src/Grover.jl @@ -14,47 +14,47 @@ inference_oracle(nbit::Int, locs::Vector{Int}) = inference_oracle(locs)(nbit) Return a mask, that disired subspace of an oracle are masked true. """ function target_space(nbit::Int, oracle) - r = register(ones(ComplexF64, 1< oracle real(statevec(r)) .< 0 end -prob_inspace(psi::DefaultRegister, ts) = norm(statevec(psi)[ts])^2 +prob_inspace(psi::ArrayReg, ts) = norm(statevec(psi)[ts])^2 """ prob_match_oracle(psi, oracle) -> Float64 Return the probability that `psi` matches oracle. """ -prob_match_oracle(psi::DefaultRegister, oracle) = prob_inspace(psi, target_space(nqubits(psi), oracle)) +prob_match_oracle(psi::ArrayReg, oracle) = prob_inspace(psi, target_space(nqubits(psi), oracle)) """ - num_grover_step(psi::DefaultRegister, oracle) -> Int + num_grover_step(psi::ArrayReg, oracle) -> Int Return number of grover steps needed to match the oracle. """ -num_grover_step(psi::DefaultRegister, oracle) = _num_grover_step(prob_match_oracle(psi, oracle)) +num_grover_step(psi::ArrayReg, oracle) = _num_grover_step(prob_match_oracle(psi, oracle)) _num_grover_step(prob::Real) = Int(round(pi/4/sqrt(prob)))-1 """ GroverIter{N, T} - GroverIter(oracle, ref::ReflectBlock{N, T}, psi::DefaultRegister, niter::Int) + GroverIter(oracle, ref::ReflectBlock{N, T}, psi::ArrayReg, niter::Int) an iterator that perform Grover operations step by step. An Grover operation consists of applying oracle and Reflection. """ struct GroverIter{N, T} - psi::DefaultRegister + psi::ArrayReg oracle ref::ReflectBlock{N, T} niter::Int end -groveriter(psi::DefaultRegister, oracle, ref::ReflectBlock{N, T}, niter::Int) where {N, T} = GroverIter{N, T}(psi, oracle, ref, niter) -groveriter(psi::DefaultRegister, oracle, niter::Int) = groveriter(psi, oracle, ReflectBlock(psi |> copy), niter) -groveriter(psi::DefaultRegister, oracle) = groveriter(psi, oracle, ReflectBlock(psi |> copy), num_grover_step(psi, oracle)) +groveriter(psi::ArrayReg, oracle, ref::ReflectBlock{N, T}, niter::Int) where {N, T} = GroverIter{N, T}(psi, oracle, ref, niter) +groveriter(psi::ArrayReg, oracle, niter::Int) = groveriter(psi, oracle, ReflectBlock(psi |> copy), niter) +groveriter(psi::ArrayReg, oracle) = groveriter(psi, oracle, ReflectBlock(psi |> copy), num_grover_step(psi, oracle)) function Base.iterate(it::GroverIter, st=1) if it.niter + 1 == st @@ -69,7 +69,7 @@ Base.length(it::GroverIter) = it.niter """ groverblock(oracle, ref::ReflectBlock{N, T}, niter::Int=-1) - groverblock(oracle, psi::DefaultRegister, niter::Int=-1) + groverblock(oracle, psi::ArrayReg, niter::Int=-1) Return a ChainBlock/Sequential as Grover Iteration, the default `niter` will stop at the first optimal step. """ @@ -83,4 +83,4 @@ function groverblock(oracle, ref::ReflectBlock{N, T}, niter::Int=-1) where {N, T sequence(sequence(oracle, ref) for i = 1:niter) end -groverblock(oracle, psi::DefaultRegister, niter::Int=-1) = groverblock(oracle, ReflectBlock(psi |> copy), niter) +groverblock(oracle, psi::ArrayReg, niter::Int=-1) = groverblock(oracle, ReflectBlock(psi |> copy), niter) diff --git a/src/HHL.jl b/src/HHL.jl index d1c8836..dc76331 100644 --- a/src/HHL.jl +++ b/src/HHL.jl @@ -22,14 +22,14 @@ end a, -b, b, a end -function apply!(reg::DefaultRegister, hr::HHLCRot{N, NC, T}) where {N, NC, T} +function apply!(reg::ArrayReg, hr::HHLCRot{N, NC, T}) where {N, NC, T} mask = bmask(hr.ibit) step = 1<<(hr.ibit-1) step_2 = step*2 nbit = nqubits(reg) for j = 0:step_2:size(reg.state, 1)-step for i = j+1:j+step - λ = bfloat(takebit(i-1, hr.cbits...), nbit=nbit-1) + λ = bfloat(readbit(i-1, hr.cbits...), nbits=nbit-1) if λ >= hr.C_value u = hhlrotmat(λ, hr.C_value) u1rows!(state(reg), i, i+step, u...) @@ -40,11 +40,11 @@ function apply!(reg::DefaultRegister, hr::HHLCRot{N, NC, T}) where {N, NC, T} end """ - hhlproject!(all_bit::DefaultRegister, n_reg::Int) -> Vector + hhlproject!(all_bit::ArrayReg, n_reg::Int) -> Vector project to aiming state |1>|00>|u>, and return |u> vector. """ -function hhlproject!(all_bit::DefaultRegister, n_reg::Int) +function hhlproject!(all_bit::ArrayReg, n_reg::Int) all_bit |> focus!(1:(n_reg+1)...) |> select!(1) |> state |> vec end @@ -61,17 +61,17 @@ end """ hhlsolve(A::Matrix, b::Vector) -> Vector - + solving linear system using HHL algorithm. Here, A must be hermitian. """ function hhlsolve(A::Matrix, b::Vector, n_reg::Int, C_value::Real) if !ishermitian(A) throw(ArgumentError("Input matrix not hermitian!")) end - UG = matrixgate(exp(2π*im.*A)) + UG = matblock(exp(2π*im.*A)) # Generating input bits - all_bit = register(b) ⊗ zero_state(n_reg) ⊗ zero_state(1) + all_bit = ArrayReg(b) ⊗ zero_state(n_reg) ⊗ zero_state(1) # Construct HHL circuit. circuit = hhlcircuit(UG, n_reg, C_value) diff --git a/src/HadamardTest.jl b/src/HadamardTest.jl index 0a9e6df..8dd8e8a 100644 --- a/src/HadamardTest.jl +++ b/src/HadamardTest.jl @@ -1,46 +1,33 @@ -export hadamard_test, hadamard_test_circuit, swap_test, swap_test_circuit, singlet_block, state_overlap_circuit +export hadamard_test, hadamard_test_circuit, swap_test_circuit """ see WiKi. """ -function hadamard_test_circuit(U::MatrixBlock{N}) where N +function hadamard_test_circuit(U::MatrixBlock{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::MatrixBlock{N}, reg::AbstractRegister) where N - c = hadamard_test_circuit(U) +function hadamard_test(U::MatrixBlock{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() = hadamard_test_circuit(SWAP) -swap_test(reg::AbstractRegister) = hadamard_test(SWAP, reg) - -function singlet_block(::Type{T}, nbit::Int, i::Int, j::Int) where T - unit = chain(nbit) - push!(unit, put(nbit, i=>chain(XGate{T}(), HGate{T}()))) - push!(unit, control(nbit, -i, j=>XGate{T}())) -end - -singlet_block(nbit::Int, i::Int, j::Int) = singlet_block(ComplexF64, nbit, i, j) -singlet_block() = singlet_block(2,1,2) - """ Estimation of overlap between multiple density matrices. PRL 88.217901 """ -function state_overlap_circuit(nbit::Int, nstate::Int, ϕ::Real) +function swap_test_circuit(nbit::Int, nstate::Int, ϕ::Real) N = nstate*nbit + 1 chain(N, put(N, 1=>H), - put(N, 1=>shift(ϕ)), + 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 - -Yao.mat(ρ::DensityMatrix{1}) = dropdims(state(ρ), dims=3) diff --git a/src/Miscellaneous.jl b/src/Miscellaneous.jl index e077d5c..68659b1 100644 --- a/src/Miscellaneous.jl +++ b/src/Miscellaneous.jl @@ -1,4 +1,4 @@ -export inverselines +export inverselines, singlet_block """ inverselines(nbit::Int; n_reg::Int=nbit) -> ChainBlock @@ -16,4 +16,13 @@ function inverselines(nbit::Int; n_reg::Int=nbit) c end +function singlet_block(::Type{T}, nbit::Int, i::Int, j::Int) where T + unit = chain(nbit) + push!(unit, put(nbit, i=>chain(XGate{T}(), HGate{T}()))) + push!(unit, control(nbit, -i, j=>XGate{T}())) +end + +singlet_block(nbit::Int, i::Int, j::Int) = singlet_block(ComplexF64, nbit, i, j) +singlet_block() = singlet_block(2,1,2) +Yao.mat(ρ::DensityMatrix{1}) = dropdims(state(ρ), dims=3) diff --git a/src/PhaseEstimation.jl b/src/PhaseEstimation.jl index ef8b148..63ba53b 100644 --- a/src/PhaseEstimation.jl +++ b/src/PhaseEstimation.jl @@ -8,7 +8,7 @@ phase estimation circuit. * `n_reg`: the number of bits to store phases, * `n_b`: the number of bits to store vector. """ -function PEBlock(UG::GeneralMatrixGate, n_reg::Int, n_b::Int) +function PEBlock(UG::GeneralMatrixBlock, n_reg::Int, n_b::Int) nbit = n_b + n_reg # Apply Hadamard Gate. hs = repeat(nbit, H, 1:n_reg) @@ -18,22 +18,22 @@ function PEBlock(UG::GeneralMatrixGate, n_reg::Int, n_b::Int) for i = 1:n_reg push!(control_circuit, control(nbit, (i,), (n_reg+1:nbit...,)=>UG)) if i != n_reg - UG = matrixgate(mat(UG) * mat(UG)) + UG = matblock(mat(UG) * mat(UG)) end end # Inverse QFT Block. - iqft = concentrate(nbit, QFTBlock{n_reg}() |> adjoint,[1:n_reg...,]) + iqft = concentrate(nbit, QFTBlock{n_reg}()',[1:n_reg...,]) chain(hs, control_circuit, iqft) end """ - projection_analysis(evec::Matrix, reg::DefaultRegister) -> Tuple + projection_analysis(evec::Matrix, reg::ArrayReg) -> Tuple Analyse using state projection. It returns a tuple of (most probable configuration, the overlap matrix, the relative probability for this configuration) """ -function projection_analysis(evec::Matrix, reg::DefaultRegister) +function projection_analysis(evec::Matrix, reg::ArrayReg) overlap = evec'*state(reg) amp_relative = Float64[] bs = Int[] diff --git a/src/QCBM.jl b/src/QCBM.jl index 8d4296a..9e4b694 100644 --- a/src/QCBM.jl +++ b/src/QCBM.jl @@ -1,4 +1,4 @@ -import Yao.Registers: probs +import Yao: probs export QCBM, QCBMGo!, psi, mmdgrad include("Kernels.jl") @@ -12,7 +12,7 @@ struct QCBM{BT<:AbstractBlock, KT<:AbstractKernel} <: QCOptProblem dbs end function QCBM(circuit::AbstractBlock, kernel::AbstractKernel, ptrain::Vector) - QCBM(circuit, kernel, ptrain, collect(circuit, AbstractDiff)) + QCBM(circuit, kernel, ptrain, collect_blocks(AbstractDiff, circuit)) end # INTERFACES diff --git a/src/QCOptProblem.jl b/src/QCOptProblem.jl index b96b696..9e85749 100644 --- a/src/QCOptProblem.jl +++ b/src/QCOptProblem.jl @@ -1,4 +1,3 @@ -import Yao: gradient export num_gradient, diff_blocks, loss, circuit, QCOptProblem, QCOptGo! """ @@ -36,7 +35,7 @@ function gradient end collection of all differentiable units. """ -diff_blocks(qop::QCOptProblem) = collect(qop |> circuit, AbstractDiff) +diff_blocks(qop::QCOptProblem) = collect_blocks(AbstractDiff, qop |> circuit) """ num_gradient(qop::QCOptProblem) -> Vector diff --git a/src/QFT.jl b/src/QFT.jl index 37123be..16250ea 100644 --- a/src/QFT.jl +++ b/src/QFT.jl @@ -5,7 +5,7 @@ end export QFTCircuit, QFTBlock, invorder_firstdim CRk(i::Int, j::Int, k::Int) = control([i, ], j=>shift(2π/(1<H) : CRk(j, i, j-i+1) for j = i:n) +CRot(n::Int, i::Int) = chain(n, i==j ? kron(i=>H) : CRk(j, i, j-i+1) for j = i:n) QFTCircuit(n::Int) = chain(n, CRot(n, i) for i = 1:n) struct QFTBlock{N} <: PrimitiveBlock{N,ComplexF64} end @@ -20,7 +20,7 @@ isreflexive(q::QFTBlock{N}) where N = N==1 isunitary(q::QFTBlock{N}) where N = true openbox(q::QFTBlock{N}) where N = QFTCircuit(N) -openbox(q::Daggered{<:QFTBlock, N}) where {N} = adjoint(QFTCircuit(N)) +openbox(q::Daggered{<:QFTBlock, N}) where {N} = QFTCircuit(N)' function print_block(io::IO, pb::QFTBlock{N}) where N printstyled(io, "QFT(1-$N)"; bold=true, color=:blue) @@ -36,7 +36,7 @@ function invorder_firstdim(v::Matrix) n_2 = n ÷ 2 mask = [bmask(i, n-i+1) for i in 1:n_2] @simd for b in basis(n) - @inbounds w[breflect(n, b, mask)+1,:] = v[b+1,:] + @inbounds w[breflect(b, mask; nbits=n)+1,:] = v[b+1,:] end w end @@ -48,7 +48,7 @@ function invorder_firstdim(v::Vector) #mask = SVector{n_2, Int}([bmask(i, n-i+1)::Int for i in 1:n_2]) mask = [bmask(i, n-i+1)::Int for i in 1:n_2] @simd for b in basis(n) - @inbounds w[breflect(n, b, mask)+1] = v[b+1] + @inbounds w[breflect(b, mask; nbits=n)+1] = v[b+1] end w end diff --git a/src/QuAlgorithmZoo.jl b/src/QuAlgorithmZoo.jl index 977b995..3b78067 100644 --- a/src/QuAlgorithmZoo.jl +++ b/src/QuAlgorithmZoo.jl @@ -1,13 +1,12 @@ module QuAlgorithmZoo using LuxurySparse, LinearAlgebra -using Yao -using Yao.Intrinsics -using Yao.Registers -using Yao.Blocks -import Yao.Blocks: mat, dispatch!, niparameters, iparameters, setiparameters!, cache_key, print_block, _make_rot_mat, apply!, PrimitiveBlock +using MacroTools: @forward +using Yao, Yao.ConstGate, BitBasis +using YaoArrayRegister: u1rows! +import Yao: mat, dispatch!, niparams, getiparams, setiparams!, cache_key, print_block, apply!, PrimitiveBlock, ishermitian, isunitary, isreflexive +import YaoBlocks: render_params import Base: ==, copy, hash -import Yao.Intrinsics: ishermitian, isreflexive, isunitary export openbox @@ -19,6 +18,8 @@ For a black box, like QFTBlock, you can get its white box (loyal simulation) usi function openbox end include("Miscellaneous.jl") +include("sequence.jl") +include("Diff.jl") include("Adam.jl") include("QFT.jl") include("CircuitBuild.jl") diff --git a/src/QuGAN.jl b/src/QuGAN.jl index 7db37b6..b6bcede 100644 --- a/src/QuGAN.jl +++ b/src/QuGAN.jl @@ -1,6 +1,8 @@ -import Yao.Registers: tracedist +using MacroTools: @forward +import Yao: tracedist export QuGAN, psi, toy_qugan, QuGANGo! + """ Quantum GAN. @@ -8,21 +10,21 @@ Reference: Benedetti, M., Grant, E., Wossnig, L., & Severini, S. (2018). Adversarial quantum circuit learning for pure state approximation, 1–14. """ struct QuGAN{N} <: QCOptProblem - target::DefaultRegister + target::ArrayReg generator::MatrixBlock{N} discriminator::MatrixBlock - reg0::DefaultRegister + reg0::ArrayReg witness_op::MatrixBlock circuit::AbstractBlock gdiffs ddiffs - function QuGAN(target::DefaultRegister, gen::MatrixBlock, dis::MatrixBlock) + function QuGAN(target::ArrayReg, gen::MatrixBlock, dis::MatrixBlock) N = nqubits(target) - c = sequence(gen, addbit(1), dis) + c = Sequence([gen, addbits!(1), dis]) witness_op = put(N+1, (N+1)=>P0) - gdiffs = collect(gen, AbstractDiff) - ddiffs = collect(dis, AbstractDiff) + gdiffs = collect_blocks(AbstractDiff, gen) + ddiffs = collect_blocks(AbstractDiff, dis) new{N}(target, gen, dis, zero_state(N), witness_op, c, gdiffs, ddiffs) end end @@ -52,11 +54,11 @@ psi_disctarget(qg::QuGAN) = copy(qg.target) |> qg.circuit[2:end] tracedist(qg::QuGAN) = tracedist(qg.target, psi(qg))[] """ - toy_qugan(target::DefaultRegister, depth_gen::Int, depth_disc::Int) -> QuGAN + toy_qugan(target::ArrayReg, depth_gen::Int, depth_disc::Int) -> QuGAN Construct a toy qugan. """ -function toy_qugan(target::DefaultRegister, depth_gen::Int, depth_disc::Int) +function toy_qugan(target::ArrayReg, depth_gen::Int, depth_disc::Int) n = nqubits(target) generator = dispatch!(random_diff_circuit(n, depth_gen, pair_ring(n)), :random) |> autodiff(:QC) discriminator = dispatch!(random_diff_circuit(n+1, depth_disc, pair_ring(n+1)), :random) |> autodiff(:QC) diff --git a/src/RotBasis.jl b/src/RotBasis.jl index bdde812..91b1067 100644 --- a/src/RotBasis.jl +++ b/src/RotBasis.jl @@ -10,6 +10,7 @@ mutable struct RotBasis{T} <: PrimitiveBlock{1, Complex{T}} phi::T end +_make_rot_mat(I, block, theta) = I * cos(theta / 2) - im * sin(theta / 2) * block # chain -> * # mat(rb::RotBasis{T}) where T = mat(Ry(-rb.theta))*mat(Rz(-rb.phi)) function mat(x::RotBasis{T}) where T @@ -23,12 +24,13 @@ end copy(block::RotBasis{T}) where T = RotBasis{T}(block.theta, block.phi) dispatch!(block::RotBasis, params::Vector) = ((block.theta, block.phi) = params; block) -iparameters(rb::RotBasis) = (rb.theta, rb.phi) -function setiparameters!(rb::RotBasis, theta::Real, phi::Real) - rb.theta, rb.phi = theta, phi +getiparams(rb::RotBasis) = (rb.theta, rb.phi) +function setiparams!(rb::RotBasis, θ::Real, ϕ::Real) + rb.theta, rb.phi = θ, ϕ rb end -niparameters(::Type{<:RotBasis}) = 2 +niparams(::Type{<:RotBasis}) = 2 +niparams(::RotBasis) = 2 render_params(r::RotBasis, ::Val{:random}) = rand()*π, rand()*2π function print_block(io::IO, R::RotBasis) diff --git a/src/hamiltonian_solvers.jl b/src/hamiltonian_solvers.jl index 7185da3..4fee6f2 100644 --- a/src/hamiltonian_solvers.jl +++ b/src/hamiltonian_solvers.jl @@ -13,12 +13,12 @@ function heisenberg(nbit::Int; periodic::Bool=true) end """ - iter_groundstate!({reg::AbstractRegister}, h::MatrixBlock; niter::Int=100) -> AbstractRegister + iter_groundstate!({reg::AbstractRegister}, h::AbstractBlock; niter::Int=100) -> AbstractRegister project wave function to ground state by iteratively apply -h. """ -iter_groundstate!(h::MatrixBlock; niter::Int=100) = reg -> iter_groundstate!(reg, h, niter=niter) -function iter_groundstate!(reg::AbstractRegister, h::MatrixBlock; niter::Int=100) +iter_groundstate!(h::AbstractBlock; niter::Int=100) = reg -> iter_groundstate!(reg, h, niter=niter) +function iter_groundstate!(reg::AbstractRegister, h::AbstractBlock; niter::Int=100) for i = 1:niter reg |> h i%5 == 0 && reg |> normalize! @@ -27,36 +27,31 @@ function iter_groundstate!(reg::AbstractRegister, h::MatrixBlock; niter::Int=100 end """ - itime_groundstate!({reg::AbstractRegister}, h::MatrixBlock; τ::Int=20, tol=1e-4) -> AbstractRegister + itime_groundstate!({reg::AbstractRegister}, h::AbstractBlock; τ::Int=20, tol=1e-4) -> AbstractRegister Imaginary time evolution method to get ground state, i.e. by projecting wave function to ground state by exp(-hτ). `tol` is for `expmv`. """ -itime_groundstate!(h::MatrixBlock; τ::Real=20, tol=1e-4) = reg -> itime_groundstate!(reg, h; τ=τ, tol=tol) -function itime_groundstate!(reg::AbstractRegister, h::MatrixBlock; τ::Int=20, tol=1e-4) - span = 1 - te = timeevolve(h, -im*span) +itime_groundstate!(h::AbstractBlock; τ::Real=20, tol=1e-4) = reg -> itime_groundstate!(reg, h; τ=τ, tol=tol) +function itime_groundstate!(reg::AbstractRegister, h::AbstractBlock; τ::Int=20, tol=1e-4) + span = 1.0 + te = time_evolve(h, -im*span) for i = 1:τ÷span reg |> te |> normalize! end if τ%span != 0 - reg |> timeevolve(h, τ%span) |> normalize! + reg |> time_evolve(h, τ%span) |> normalize! end reg end -# a patch for Yao.expect, to make it faster -function Yao.expect(op::AddBlock, reg::AbstractRegister{1}) - sum(opi->expect(opi, reg), op) -end - """ - vqe_solve!(circuit::MatrixBlock{N}, hamiltonian::AbstractBlock; niter::Int=100) -> circuit + vqe_solve!(circuit::AbstractBlock{N}, hamiltonian::AbstractBlock; niter::Int=100) -> circuit variational quantum eigensolver, faithful simulation with optimizer Adam(lr=0.01). """ -function vqe_solve!(circuit::MatrixBlock{N}, hamiltonian::AbstractBlock; niter::Int=100) where N +function vqe_solve!(circuit::AbstractBlock{N}, hamiltonian::AbstractBlock; niter::Int=100) where N optimizer = Adam(lr=0.01) - dbs = collect(circuit, AbstractDiff) + dbs = collect_blocks(AbstractDiff, circuit) params = parameters(circuit) for i = 1:niter grad = opdiff.(()->zero_state(N) |> circuit, dbs, Ref(hamiltonian)) diff --git a/src/sequence.jl b/src/sequence.jl new file mode 100644 index 0000000..12924a1 --- /dev/null +++ b/src/sequence.jl @@ -0,0 +1,39 @@ +export Sequence +import YaoBlocks: subblocks, chsubblocks, apply! +using YaoBlocks: _check_size + +struct Sequence <: CompositeBlock{Any, Bool} + blocks::Vector +end + +Sequence(args...) = Sequence(collect(AbstractBlock, args)) + +subblocks(seq::Sequence) = filter(x->x isa AbstractBlock, seq.blocks) +chsubblocks(pb::Sequence, blocks::Vector) = Sequence(blocks) + +function apply!(reg::ArrayReg, seq::Sequence) + for x in seq.blocks + reg |> x + end + reg +end + +for PROP in [:lastindex, :firstindex, :getindex, :length, :eltype, :iterate, :eachindex, :popfirst!, :pop!] + @eval Base.$PROP(c::Sequence, args...; kwargs...) = $PROP(c.blocks, args...; kwargs...) +end + +function Base.:(==)(lhs::Sequence, rhs::Sequence) + (length(lhs.blocks) == length(rhs.blocks)) && all(lhs.blocks .== rhs.blocks) +end + +Base.copy(c::Sequence) = Sequence(copy(c.blocks)) +Base.similar(c::Sequence) = Sequence(empty!(similar(c.blocks))) +Base.getindex(c::Sequence, index::Union{UnitRange, Vector}) = Sequence(getindex(c.blocks, index)) + +Base.setindex!(c::Sequence, val::AbstractBlock{N}, index::Integer) where N = (setindex!(c.blocks, val, index); c) +Base.insert!(c::Sequence, index::Integer, val::AbstractBlock{N}) where N = (insert!(c.blocks, index, val); c) +Base.push!(c::Sequence, m) where N = (push!(c.blocks, m); c) +Base.append!(c::Sequence, list::Vector) where N = (append!(c.blocks, list); c) +Base.append!(c1::Sequence, c2::Sequence) where N = (append!(c1.blocks, c2.blocks); c1) +Base.prepend!(c1::Sequence, list::Vector{<:AbstractBlock{N}}) where N = (prepend!(c1.blocks, list); c1) +Base.prepend!(c1::Sequence, c2::Sequence) where N = (prepend!(c1.blocks, c2.blocks); c1) diff --git a/test/CircuitBuild.jl b/test/CircuitBuild.jl index b87e6c2..d64687c 100644 --- a/test/CircuitBuild.jl +++ b/test/CircuitBuild.jl @@ -1,5 +1,5 @@ using Test -using Yao, Yao.Blocks, QuAlgorithmZoo +using Yao, QuAlgorithmZoo @testset "pairs geometries" begin @test pair_ring(3) == [1=>2,2=>3,3=>1] @@ -21,22 +21,22 @@ end @test length(c) == 45 end -@testset "rotter, collect_rotblocks, num_gradient, opgrad" begin +@testset "rotter, collect_blocks, num_gradient, opgrad" begin @test merged_rotor(true, true) == Rx(0) @test merged_rotor(false, false) == merged_rotor() == chain(Rz(0), Rx(0), Rz(0)) @test merged_rotor(false, true) == chain(Rz(0), Rx(0)) @test merged_rotor(true, false) == chain(Rx(0), Rz(0)) - @test collect(rotorset(:Merged, 5, true, false), RotationGate) |> length == 10 + @test collect_blocks(RotationGate, rotorset(:Merged, 5, true, false)) |> length == 10 @test rotor(5, 2, true, true) isa ChainBlock @test rotor(5, 2, true, true) |> length == 1 @test rotor(5, 2, true, true) |> nqubits == 5 - @test collect(rotorset(:Split, 5, true, false), PutBlock{<:Any, <:Any, <:RotationGate}) |> length == 10 + @test collect_blocks(PutBlock{<:Any, <:Any, <:RotationGate}, rotorset(:Split, 5, true, false)) |> length == 10 end @testset "random diff circuit" begin c = random_diff_circuit(4, 3, [1=>3, 2=>4, 2=>3, 4=>1]) - rots = collect(c, RotationGate) + rots = collect_blocks(RotationGate, c) @test length(rots) == nparameters(c) == 40 @test dispatch!(+, c, ones(40)*0.1) |> parameters == ones(40)*0.1 @test dispatch!(+, c, :random) |> parameters != ones(40)*0.1 @@ -46,7 +46,7 @@ end reg = rand_state(4) dispatch!(c, randn(nparameters(c))) - dbs = collect(c, BPDiff) + dbs = collect_blocks(BPDiff, c) op = kron(4, 1=>Z, 2=>X) loss1z() = expect(op, copy(reg) |> c) # return loss please diff --git a/test/Diff.jl b/test/Diff.jl new file mode 100644 index 0000000..d69d15f --- /dev/null +++ b/test/Diff.jl @@ -0,0 +1,172 @@ +using Yao, QuAlgorithmZoo +using Yao.ConstGate +using LinearAlgebra, Test, Random + +@testset "BP diff" begin + reg = rand_state(4) + block = put(4, 2=>rot(X, 0.3)) + df = BPDiff(block) + @test df.grad == 0 + @test nqubits(df) == 4 + + df2 = BPDiff(rot(CNOT, 0.3)) + @test df2.grad == 0 + @test nqubits(df2) == 2 +end + +@testset "Qi diff" begin + reg = rand_state(4) + df2 = QDiff(rot(CNOT, 0.3)) + @test df2.grad == 0 + @test nqubits(df2) == 2 + + @test df2' isa QDiff + @test mat(df2) == mat(df2')' +end + +""" + loss_expect!(circuit::AbstractBlock, op::AbstractBlock) -> Function + +Return function "loss!(ψ, θ) -> Vector" +""" +function loss_expect!(circuit::AbstractBlock, op::AbstractBlock) + N = nqubits(circuit) + function loss!(ψ::AbstractRegister, θ::Vector) + params = parameters(circuit) + dispatch!(circuit, θ) + ψ |> circuit + popdispatch!(circuit, params) + expect(op, ψ) + end +end + +""" + loss_Z1!(circuit::AbstractBlock; ibit::Int=1) -> Function + +Return the loss function f = (means measuring the ibit-th bit in computation basis). +""" +loss_Z1!(circuit::AbstractBlock; ibit::Int=1) = loss_expect!(circuit, put(nqubits(circuit), ibit=>Z)) + +cnot_entangler(n::Int, pairs) = chain(n, control(n, [ctrl], target=>X) for (ctrl, target) in pairs) + +function rotor(nbit::Int, ibit::Int, noleading::Bool=false, notrailing::Bool=false) + rt = chain(nbit, [put(nbit, ibit=>Rz(0.0)), put(nbit, ibit=>Rx(0.0)), put(nbit, ibit=>Rz(0.0))]) + noleading && popfirst!(rt) + notrailing && pop!(rt) + rt +end + +rset(nbit::Int, noleading::Bool=false, notrailing::Bool=false) = chain(nbit, [rotor(nbit, j, noleading, notrailing) for j=1:nbit]) + +function ibm_diff_circuit(nbit, nlayer, pairs) + circuit = chain(nbit) + + ent = cnot_entangler(nbit, pairs) + for i = 1:(nlayer + 1) + i!=1 && push!(circuit, ent) + push!(circuit, rset(nbit, i==1, i==nlayer+1)) + end + circuit +end + +@testset "BP diff" begin + c = put(4, 3=>Rx(0.5)) |> autodiff(:BP) + cad = c' + @test mat(cad) == mat(c)' + + circuit = chain(4, repeat(4, H, 1:4), put(4, 3=>Rz(0.5)) |> autodiff(:BP), control(4, 2, 1=>shift(0.4)) |> autodiff(:BP), control(2, 1=>X), put(4, 4=>Ry(0.2)) |> autodiff(:BP)) + op = put(4, 3=>Y) + θ = [0.1, 0.2, 0.3] + dispatch!(circuit, θ) + loss! = loss_expect!(circuit, op) + ψ0 = rand_state(4) + ψ = copy(ψ0) |> circuit + + # get gradient + δ = ψ |> op + backward!(δ, circuit) + g1 = gradient(circuit) + + g2 = zero(θ) + η = 0.01 + for i in 1:length(θ) + θ1 = copy(θ) + θ2 = copy(θ) + θ1[i] -= 0.5η + θ2[i] += 0.5η + g2[i] = (loss!(copy(ψ0), θ2) - loss!(copy(ψ0), θ1))/η |> real + end + g3 = opdiff.(() -> copy(ψ0) |> circuit, collect_blocks(BPDiff, circuit), Ref(op)) + @test isapprox.(g1, g2, atol=1e-5) |> all + @test isapprox.(g2, g3, atol=1e-5) |> all +end + +@testset "constructor" begin + @test generator(put(4, 1=>Rx(0.1))) == put(4, 1=>X) + @test generator(Rx(0.1)) == X + circuit = chain(put(4, 1=>Rx(0.1)), control(4, 2, 1=>Ry(0.3))) + c2 = circuit |> autodiff(:BP) + @test c2[1] isa BPDiff + @test !(c2[2] isa BPDiff) +end + +@testset "numdiff & opdiff" begin + @test collect_blocks(XGate, chain([X, Y, Z])) == [X] + + c = chain(put(4, 1=>Rx(0.5))) |> autodiff(:QC) + nd = numdiff(c[1].content) do + expect(put(4, 1=>Z), zero_state(4) |> c) |> real # return loss please + end + + ed = opdiff(c[1].content, put(4, 1=>Z)) do + zero_state(4) |> c # a function get output + end + @test isapprox(nd, ed, atol=1e-4) + + reg = rand_state(4) + c = chain(put(4, 1=>Rx(0.5)), control(4, 1, 2=>Ry(0.5)), control(4, 1, 2=>shift(0.3)), kron(4, 2=>Rz(0.3), 3=>Rx(0.7))) |> autodiff(:QC) + dbs = collect_blocks(QDiff, c) + loss1z() = expect(kron(4, 1=>Z, 2=>X), copy(reg) |> c) |> real # return loss please + nd = numdiff.(loss1z, dbs) + ed = opdiff.(()->copy(reg) |> c, dbs, Ref(kron(4, 1=>Z, 2=>X))) + gd = gradient(c) + @test gradient(c, :QC) == gd + @test gradient(c, :BP) == [] + @test isapprox(nd, ed, atol=1e-4) + @test ed == gd +end + +@testset "stat" begin + nbit = 3 + f(x::Number, y::Number) = Float64(abs(x-y) < 1.5) + x = 0:1<2, 2=>3, 3=>1] + c = ibm_diff_circuit(nbit, 2, prs) |> autodiff(:QC) + dispatch!(c, :random) + dbs = collect_blocks(AbstractDiff,c) + + p0 = zero_state(nbit) |> c |> probs + sample0 = measure(zero_state(nbit) |> c; nshots=5000) + loss0 = expect(V, p0 |> as_weights) + gradsn = numdiff.(()->expect(V, zero_state(nbit) |> c |> probs |> as_weights), dbs) + gradse = statdiff.(()->zero_state(nbit) |> c |> probs |> as_weights, dbs, Ref(V), initial=p0 |> as_weights) + gradsf = statdiff.(()->measure(zero_state(nbit) |> c; nshots=5000), dbs, Ref(VF), initial=sample0) + @test all(isapprox.(gradse, gradsn, atol=1e-4)) + @test norm(gradsf-gradse)/norm(gradsf) <= 0.2 + + # 1D + h = randn(1< autodiff(:QC) + dispatch!(c, :random) + dbs = collect_blocks(AbstractDiff, c) + + p0 = zero_state(nbit) |> c |> probs |> as_weights + loss0 = expect(V, p0 |> as_weights) + gradsn = numdiff.(()->expect(V, zero_state(nbit) |> c |> probs |> as_weights), dbs) + gradse = statdiff.(()->zero_state(nbit) |> c |> probs |> as_weights, dbs, Ref(V)) + @test all(isapprox.(gradse, gradsn, atol=1e-4)) +end diff --git a/test/Grover.jl b/test/Grover.jl index f8e06a1..968808d 100644 --- a/test/Grover.jl +++ b/test/Grover.jl @@ -1,9 +1,9 @@ using Test, Random, LinearAlgebra, SparseArrays +using BitBasis using QuAlgorithmZoo import QuAlgorithmZoo: _num_grover_step -using Yao, Yao.Blocks -using Yao.Intrinsics +using Yao function GroverSearch(oracle, num_bit::Int; psi::DefaultRegister = uniform_state(num_bit)) it = groveriter(psi, oracle) @@ -42,7 +42,7 @@ end @testset "groverblock" begin psi = uniform_state(5) or = inference_oracle(5, [-1,2,5,4,3]) - func_or = FunctionBlock{:Oracle}(reg->apply!(reg, or)) + func_or = or gb = groverblock(or, psi) gb2 = groverblock(func_or, psi) @test apply!(copy(psi), gb) == (for l_psi in groveriter(copy(psi), func_or) psi = l_psi end; psi) @@ -58,7 +58,7 @@ end # the desired subspace basis = collect(UInt, 0:1<0)) + subinds = [itercontrol(num_bit, abs.(evidense), Int.(evidense.>0))...] v_desired = statevec(psi0)[subinds .+ 1] p = norm(v_desired)^2 @@ -67,5 +67,5 @@ end # search the subspace num_iter = _num_grover_step(p) niter, psi = inference(psi0, evidense, num_iter) - @test isapprox((psi.state[subinds .+ 1]'*v_desired) |> abs2, 1, atol=1e-2) + @test isapprox((psi.state[subinds .+ 1]'*v_desired) |> abs2, 1, atol=3e-2) end diff --git a/test/HHL.jl b/test/HHL.jl index 391a4dd..287c9c0 100644 --- a/test/HHL.jl +++ b/test/HHL.jl @@ -1,5 +1,5 @@ using Yao -using Yao.Intrinsics +using BitBasis using QuAlgorithmZoo using Test, LinearAlgebra @@ -11,7 +11,7 @@ function crot(n_reg::Int, C_value::Real) c_bit = Vector(2:n_rot) λ = 0.0 for j = 1:n_reg - if (takebit(i,j) == 0) + if (readbit(i,j) == 0) c_bit[j] = -c_bit[j] end end diff --git a/test/HadamardTest.jl b/test/HadamardTest.jl index 2dba4dd..9cb205b 100644 --- a/test/HadamardTest.jl +++ b/test/HadamardTest.jl @@ -1,7 +1,11 @@ -using Yao, Yao.Blocks +using Yao using Test using LinearAlgebra using QuAlgorithmZoo +using Yao.ConstGate + +single_swap_test_circuit(ϕ::Real) = hadamard_test_circuit(SWAP, ϕ) +single_swap_test(reg::AbstractRegister, ϕ::Real) = hadamard_test(SWAP, reg, ϕ) @testset "state overlap" begin reg1 = rand_state(3) |> focus!(1,2) @@ -11,15 +15,15 @@ using QuAlgorithmZoo reg3 = rand_state(3) |> focus!(1,2) rho3 = reg3 |> ρ desired = tr(mat(rho1)*mat(rho2)) - c = state_overlap_circuit(2, 2, 0) + c = swap_test_circuit(2, 2, 0) res = expect(put(5, 1=>Z), join(join(reg2, reg1), zero_state(1)) |> c) |> tr @test desired ≈ res desired = tr(mat(rho1)*mat(rho2)*mat(rho3)) |> real - c = state_overlap_circuit(2, 3, 0) + c = swap_test_circuit(2, 3, 0) res = expect(put(7, 1=>Z), reduce(⊗, [reg3, reg2, reg1, zero_state(1)]) |> c) |> tr |> real @test desired ≈ res desired = tr(mat(rho1)*mat(rho2)*mat(rho3)) |> imag - c = state_overlap_circuit(2, 3, -π/2) + c = swap_test_circuit(2, 3, -π/2) res = expect(put(7, 1=>Z), reduce(⊗, [reg3, reg2, reg1, zero_state(1)]) |> c) |> tr |> real @test desired ≈ res end @@ -29,13 +33,14 @@ end U = put(nbit, 2=>Rx(0.2)) reg = rand_state(nbit) - @test hadamard_test(U, reg) ≈ real(expect(U, reg)) + @test hadamard_test(U, reg, 0.0) ≈ real(expect(U, reg)) + @test hadamard_test(U, reg, -π/2) ≈ imag(expect(U, reg)) reg = zero_state(2) |> singlet_block() - @test swap_test(reg) ≈ -1 + @test single_swap_test(reg, 0) ≈ -1 reg = zero_state(2) - @test swap_test(reg) ≈ 1 + @test single_swap_test(reg, 0) ≈ 1 reg = product_state(2, 0b11) - @test swap_test(reg) ≈ 1 + @test single_swap_test(reg, 0) ≈ 1 end diff --git a/test/PhaseEstimation.jl b/test/PhaseEstimation.jl index 0db8d35..1aa3f27 100644 --- a/test/PhaseEstimation.jl +++ b/test/PhaseEstimation.jl @@ -1,5 +1,5 @@ using Yao -using Yao.Intrinsics +using BitBasis using Test, LinearAlgebra using QuAlgorithmZoo @@ -27,11 +27,11 @@ end U = V*Diagonal(signs)*V' b = V[:,3] - # Define registers and U operator. + # Define ArrayReg and U operator. M = 6 reg1 = zero_state(M) - reg2 = register(b) - UG = matrixgate(U) + reg2 = ArrayReg(b) + UG = matblock(U) # circuit circuit = PEBlock(UG, M, N) @@ -40,7 +40,7 @@ end reg = apply!(join(reg2, reg1), circuit) # measure - res = breflect(M, measure(focus!(copy(reg), 1:M); nshot=10)[1]) / (1< adjoint) ≈ join(reg2, reg1) @@ -51,11 +51,11 @@ end N = 3 U, b, ϕs, evec = rand_phaseest_setup(N) - # Define registers and U operator. + # Define ArrayReg and U operator. M = 6 reg1 = zero_state(M) - reg2 = register(b) - UG = matrixgate(U); + reg2 = ArrayReg(b) + UG = matblock(U); # run circuit reg= join(reg2, reg1) @@ -64,5 +64,5 @@ end # measure bs, proj, amp_relative = projection_analysis(evec, focus!(reg, M+1:M+N)) - @test isapprox(ϕs, bfloat.(bs, nbit=M), atol=0.05) + @test isapprox(ϕs, bfloat.(bs, nbits=M), atol=0.05) end diff --git a/test/QFT.jl b/test/QFT.jl index e17d0ea..a2961fe 100644 --- a/test/QFT.jl +++ b/test/QFT.jl @@ -1,13 +1,14 @@ using Test, Random, LinearAlgebra, SparseArrays, LuxurySparse using Yao +using YaoArrayRegister: invorder using QuAlgorithmZoo using FFTW @testset "QFT" begin num_bit = 5 fftblock = QFTCircuit(num_bit) - ifftblock = adjoint(fftblock) + ifftblock = fftblock' reg = rand_state(num_bit) rv = copy(statevec(reg)) diff --git a/test/QuGAN.jl b/test/QuGAN.jl index 9fffa5b..17c81e2 100644 --- a/test/QuGAN.jl +++ b/test/QuGAN.jl @@ -1,6 +1,7 @@ using Test using Yao using QuAlgorithmZoo +using Random function run_test(nbit::Int, depth_gen::Int, depth_disc::Int; g_lr=0.1, d_lr=0.2, niter=1000) qg = toy_qugan(rand_state(nbit), depth_gen, depth_disc) @@ -8,7 +9,9 @@ function run_test(nbit::Int, depth_gen::Int, depth_disc::Int; g_lr=0.1, d_lr=0.2 qg end +# to fix @testset "quantum circuit gan - opdiff" begin + Random.seed!(2) N = 3 target = rand_state(N) qcg = toy_qugan(target, 2, 2) @@ -16,7 +19,6 @@ end @test isapprox(grad, num_gradient(qcg), atol=1e-4) qg = run_test(3, 4, 4, g_lr=0.2, d_lr=0.5, niter=300) @test qg |> loss < 0.1 - qg = run_test(3, 4, 4, g_lr=Adam(lr=0.005), d_lr=Adam(lr=0.5), niter=300) + qg = run_test(3, 4, 4, g_lr=Adam(lr=0.005), d_lr=Adam(lr=0.5), niter=1000) @test qg |> loss < 0.1 end - diff --git a/test/RotBasis.jl b/test/RotBasis.jl index 6ad3ed6..3921eeb 100644 --- a/test/RotBasis.jl +++ b/test/RotBasis.jl @@ -18,7 +18,7 @@ using QuAlgorithmZoo rb = roll(1, RotBasis(0.1, 0.3))#rot_basis(1) angles = randpolar(1) # prepair a state in the angles direction. - psi = angles |> polar2u |> register + psi = angles |> polar2u |> ArrayReg # rotate to the same direction for measurements. dispatch!(rb, vec(angles)) diff --git a/test/Sequence.jl b/test/Sequence.jl new file mode 100644 index 0000000..35feb76 --- /dev/null +++ b/test/Sequence.jl @@ -0,0 +1,62 @@ +using Test + +using Yao +using QuAlgorithmZoo + +@testset "constructor" begin + + g = Sequence( + kron(2, X, Y), + kron(2, 1=>phase(0.1)), + ) + + @test g isa Sequence + @test g.blocks == [kron(2, X, Y), kron(2, 1=>phase(0.1))] +end + +@testset "apply" begin + g = Sequence( + kron(2, X, Y), + kron(2, 1=>phase(0.1)), + ) + + reg = rand_state(2) + @test statevec(apply!(copy(reg), g)) ≈ mat(chain(g...)) * reg.state +end + +@testset "iteration" begin + test_list = [X, Y, phase(0.1), rot(X, 0.0)] + g = Sequence(test_list) + + for (src, tg) in zip(g, test_list) + @test src == tg + end + + for (src, tg) in zip(eachindex(g), 1:length(test_list)) + @test src == tg + end +end + +@testset "additional" begin + g = Sequence(X, Y) + push!(g, Z) + @test g[3] == Z + + append!(g, [rot(X, 0.0), rot(Y, 0.0)]) + @test g[4] == rot(X, 0.0) + @test g[5] == rot(Y, 0.0) + + prepend!(g, [phase(0.1)]) + @test g[1] == phase(0.1) + @test g[2] == X + @test g[end] == rot(Y, 0.0) + gg = insert!(g, 4, Z) + @test gg[4] == Z +end + +@testset "traits" begin + # TODO: check traits when primitive blocks' traits are all defined + g = Sequence(X, Y) + @test length(g) == 2 + @test eltype(g) == eltype(g.blocks) +end diff --git a/test/hamiltonian_solvers.jl b/test/hamiltonian_solvers.jl index fa64cab..5b6e5ae 100644 --- a/test/hamiltonian_solvers.jl +++ b/test/hamiltonian_solvers.jl @@ -1,4 +1,4 @@ -using Yao, Yao.Blocks +using Yao using LinearAlgebra using Test using QuAlgorithmZoo diff --git a/test/lin_diffEq_test.jl b/test/lin_diffEq_test.jl index 4cd420b..11a7af8 100644 --- a/test/lin_diffEq_test.jl +++ b/test/lin_diffEq_test.jl @@ -1,5 +1,6 @@ using Yao -using Yao.Intrinsics +using BitBasis +using Random using QuAlgorithmZoo using Test, LinearAlgebra using OrdinaryDiffEq @@ -14,6 +15,7 @@ function diffeq_problem(nbit::Int) end @testset "Linear_differential_equations_HHL" begin + Random.seed!(2) N = 1 h = 0.1 tspan = (0.0,0.6) diff --git a/test/runtests.jl b/test/runtests.jl index e2c4a88..100b854 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -28,6 +28,9 @@ end @testset "HHL" begin include("HHL.jl") end +@testset "diff Eq" begin + include("lin_diffEq_test.jl") +end @testset "QCOptProblem" begin include("QCOptProblem.jl") @@ -40,3 +43,11 @@ end @testset "hadamard test" begin include("HadamardTest.jl") end + +@testset "Sequence" begin + include("Sequence.jl") +end + +@testset "Diff" begin + include("Diff.jl") +end