Merge branch 'master' of github.com:QuantumBFS/Yao.jl

docs/README.md

@@ -3,5 +3,5 @@
 use the following command to build locally.
 
 ```sh
-> julia make.jl loacl
+> julia make.jl local
 ```

docs/make.jl

@@ -18,6 +18,7 @@ makedocs(
         "Tutorial" => Any[
             "tutorial/GHZ.md",
             "tutorial/QFT.md",
+            "tutorial/Grover.md",
             "tutorial/QCBM.md",
         ],
         "Manual" => Any[

docs/src/tutorial/Grover.md

@@ -1,162 +1,172 @@
-# Grover Search and Inference
+# Grover Search and Quantum Inference
 
 ## Grover Search
 ![grover](../assets/figures/grover.png)
 
+First, we construct the reflection block ``R(|\psi\rangle) = 2|\psi\rangle\langle\psi|-1``, given we know how to construct ``|\psi\rangle=A|0\rangle``.
+Then it equivalent to construct $R(|\psi\rangle) = A(2|0\rangle\langle 0|-1)A^\dagger$
 ```@example Grover
 using Yao
+using Yao.Blocks
+using Compat
+using Compat.Test
+using StatsBase
+
+"""
+A way to construct oracle, e.g. inference_oracle([1,2,-3,5]) will
+invert the sign when a qubit configuration matches: 1=>1, 2=>1, 3=>0, 5=>1.
+"""
+function inference_oracle(locs::Vector{Int})
+    control(locs[1:end-1], abs(locs[end]) => (locs[end]>0 ? Z : chain(phase(π), Z)))
+end
 
-# Control-R(k) gate in block-A
-A(i::Int, j::Int, k::Int) = control([i, ], j=>shift(2π/(1<<k)))
-# block-B
-B(n::Int, i::Int) = chain(i==j ? kron(i=>H) : A(j, i, j-i+1) for j = i:n)
-QFT(n::Int) = chain(n, B(n, i) for i = 1:n)
-
-# define QFT and IQFT block.
-num_bit = 5
-qft = QFT(num_bit)
-iqft = adjoint(qft)
-```
-
-The basic building block - controled phase shift gate is defined as
-
-```math
-R(k)=\begin{bmatrix}
-1 & 0\\
-0 & \exp\left(\frac{2\pi i}{2^k}\right)
-\end{bmatrix}
-```
-In Yao, factory methods for blocks will be loaded lazily. For example, if you missed the total
-number of qubits of `chain`, then it will return a function that requires an input of an integer.
-So the following two statements are equivalent
-```@example QFT
-control([4, ], 1=>shift(-2π/(1<<4)))(5) == control(5, [4, ], 1=>shift(-2π/(1<<4)))
-```
-Both of then will return a `ControlBlock` instance. If you missed the total number of qubits. It is OK. Just go on, it will be filled when its possible.
-
-Once you have construct a block, you can inspect its matrix using `mat` function.
-Let's construct the circuit in dashed box A, and see the matrix of ``R_4`` gate
-```julia
-julia> a = A(4, 1, 4)(5)
-Total: 5, DataType: Complex{Float64}
-control(4)
-└─ 1=>Phase Shift Gate:-0.39269908169872414
+function reflectblock(A::MatrixBlock{N}) where N
+    chain(N, A |> adjoint, inference_oracle(-collect(1:N)), A)
+end
 
+nbit = 12
+A = repeat(nbit, H)
+ref = reflectblock(A)
 
-julia> mat(a.block)
-2×2 Diagonal{Complex{Float64}}:
- 1.0+0.0im          ⋅         
-     ⋅      0.92388-0.382683im
+@testset "test reflect" begin
+    reg = rand_state(nbit)
+    ref_vec = apply!(zero_state(nbit), A) |> statevec
+    v0 = reg |> statevec
+    @test -2*(ref_vec'*v0)*ref_vec + v0 ≈ apply!(copy(reg), ref) |> statevec
+end
 ```
+Then we define the oracle and target state
 
-Similarly, you can use `put` and `chain` to construct `PutBlock` (basic placement of a single gate) and `ChainBlock` (sequential application of `MatrixBlock`s) instances. `Yao.jl` view every component in a circuit as an `AbstractBlock`, these blocks can be integrated to perform higher level functionality.
-
-You can check the result using classical `fft`
-```@example QFT
-# if you're using lastest julia, you need to add the fft package.
-@static if VERSION >= v"0.7-"
-    using FFTW
+```@example Grover
+# first, construct the oracle with desired state in the range 100-105.
+oracle!(reg::DefaultRegister) = (reg.state[100:105,:]*=-1; reg)
+
+# transform it into a function block, so it can be put inside a `Sequential`.
+fb_oracle = FunctionBlock{:Oracle}(reg->oracle!(reg))
+
+"""
+ratio of components in a wavefunction that flip sign under oracle.
+"""
+function prob_match_oracle(psi::DefaultRegister, oracle)
+    fliped_reg = apply!(register(ones(Complex128, 1<<nqubits(psi))), oracle)
+    match_mask = fliped_reg |> statevec |> real .< 0
+    norm(statevec(psi)[match_mask])^2
 end
-using Compat.Test
 
-@test chain(num_bit, qft, iqft) |> mat ≈ eye(2^num_bit)
+# uniform state as initial state
+psi0 = apply!(zero_state(nbit), A)
 
-# define a register and get its vector representation
-reg = rand_state(num_bit)
-rv = reg |> statevec |> copy
+# the number of grover steps that can make it reach first maximum overlap.
+num_grover_step(prob::Real) = Int(round(pi/4/sqrt(prob)))-1
+niter = num_grover_step(prob_match_oracle(psi0, fb_oracle))
 
-# test fft
-reg_qft = apply!(copy(reg) |>invorder!, qft)
-kv = ifft(rv)*sqrt(length(rv))
-@test reg_qft |> statevec ≈ kv
-
-# test ifft
-reg_iqft = apply!(copy(reg), iqft)
-kv = fft(rv)/sqrt(length(rv))
-@test reg_iqft |> statevec ≈ kv |> invorder
+# construct the whole circuit
+gb = sequence(sequence(fb_oracle, ref) for i = 1:niter);
 ```
 
-QFT and IQFT are different from FFT and IFFT in three ways,
+Now, let's start training
+```@example Grover
+for (i, blk) in enumerate(gb)
+    apply!(psi0, blk)
+    overlap = prob_match_oracle(psi0, fb_oracle)
+    println("step $i, overlap = $overlap")
+end
+```
 
-1. they are different by a factor of ``\sqrt{2^n}`` with ``n`` the number of qubits.
-2. the little end and big end will exchange after applying QFT or IQFT.
-3. due to the convention, QFT is more related to IFFT rather than FFT.
+The above is the standard Grover Search algorithm, it can find target state in $O(\sqrt N)$ time, with $N$ the size of an unordered database.
+Similar algorithm can be used in more useful applications, like inference, i.e. get conditional probability distribution $p(x|y)$ given $p(x, y)$.
 
+```@example Grover
+function rand_circuit(nbit::Int, ngate::Int)
+    circuit = chain(nbit)
+    gate_list = [X, H, Ry(0.3), CNOT]
+    for i = 1:ngate
+        gate = rand(gate_list)
+        push!(circuit, put(nbit, (sample(1:nbit, nqubits(gate),replace=false)...,)=>gate))
+    end
+    circuit
+end
+A = rand_circuit(nbit, 200)
+psi0 = apply!(zero_state(nbit), A)
+
+# now we want to search the subspace with [1,3,5,8,9,11,12]
+# fixed to 1 and [4,6] fixed to 0.
+evidense = [1, 3, -4, 5, -6, 8, 9, 11, 12]
+
+"""
+Doing Inference, psi is the initial state,
+the target is to search target space with specific evidense.
+e.g. evidense [1, -3, 6] means the [1, 3, 6]-th bits take value [1, 0, 1].
+"""
+oracle_infer = inference_oracle(evidense)(nqubits(psi0))
+
+niter = num_grover_step(prob_match_oracle(psi0, oracle_infer))
+gb_infer = chain(nbit, chain(oracle_infer, reflectblock(A)) for i = 1:niter);
+```
 
-## Phase Estimation
-Since we have QFT and IQFT blocks we can then use them to realize phase estimation circuit, what we want to realize is the following circuit
-![phase estimation](../assets/figures/phaseest.png)
+Now, let's start training
+```@example Grover
+for (i, blk) in enumerate(gb_infer)
+    apply!(psi0, blk)
+    p_target = prob_match_oracle(psi0, oracle_infer)
+    println("step $i, overlap^2 = $p_target")
+end
+```
 
-In the following simulation, we use equivalent `QFTBlock` in the Yao.`Zoo` module rather than the above chain block,
-it is faster than the above construction because it hides all the simulation details (yes, we are cheating :D) and get the equivalent output.
+Here is an application, suppose we have constructed some digits and stored it in a wave vector.
 
-```@example PhaseEstimation
-using Yao
-using Yao.Zoo
-using Yao.Blocks
+```@example Grover
 using Yao.Intrinsics
 
-function phase_estimation(reg1::DefaultRegister, reg2::DefaultRegister, U::GeneralMatrixGate{N}, nshot::Int=1) where {N}
-    M = nqubits(reg1)
-    iqft = QFTBlock{M}() |> adjoint
-    HGates = rollrepeat(M, H)
-
-    control_circuit = chain(M+N)
-    for i = 1:M
-        push!(control_circuit, control(M+N, (i,), (M+1:M+N...,)=>U))
-        if i != M
-            U = matrixgate(mat(U) * mat(U))
-        end
-    end
+x1 = [0 1 0; 0 1 0; 0 1 0; 0 1 0; 0 1 0]
+x2 = [1 1 1; 0 0 1; 1 1 1; 1 0 0; 1 1 1]
+x0 = [1 1 1; 1 0 1; 1 0 1; 1 0 1; 1 1 1]
 
-    # calculation
-    # step1 apply hadamard gates.
-    apply!(reg1, HGates)
-    # join two registers
-    reg = join(reg1, reg2)
-    # using iqft to read out the phase
-    apply!(reg, sequence(control_circuit, focus(1:M...), iqft))
-    # measure the register (on focused bits), if the phase can be exactly represented by M qubits, only a single shot is needed.
-    res = measure(reg, nshot)
-    # inverse the bits in result due to the exchange of big and little ends, so that we can get the correct phase.
-    breflect.(M, res)./(1<<M), reg
+nbit = 15
+v = zeros(1<<nbit)
+
+# they occur with different probabilities.
+for (x, p) in [(x0, 0.7), (x1, 0.29), (x2,0.01)]
+    v[(x |> vec |> BitArray |> packbits)+1] = sqrt(p)
 end
 ```
-Here, `reg1` (``Q_{1-5}``) is used as the output space to store phase ϕ, and `reg2` (``Q_{6-8}``) is the input state which corresponds to an eigenvector of oracle matrix `U`.
-The algorithm detials can be found [here](https://en.wikipedia.org/wiki/Quantum_phase_estimation_algorithm).
+Plot them, you will see these digits
 
-In this function, `HGates` corresponds to circuit block in dashed box `A`, `control_circuit` corresponds to block in dashed box `B`.
-`matrixgate` is a factory function for `GeneralMatrixGate`.
+![digits](../assets/figures/digits012.png)
 
-Here, the only difficult concept is `focus`, `focus` returns a `FunctionBlock`, that will make focused bits the active bits.
-An operator sees only active bits, and operating active space is more efficient, most importantly, it becomes much easier to integrate blocks.
-However, it has the potential ability to change line orders, for safety consideration, you may also need safer [`Concentrator`](@ref).
+Then we construct the inference circuit.
+Here, we choose to use `reflect` to construct a [`ReflectBlock`](@ref),
+instead of constructing it explicitly.
+```@example Grover
+rb = reflect(copy(v))
+psi0 = register(v)
 
-```@example PhaseEstimation
-r = rand_state(6)
-apply!(r, focus(4,1,2))  # or equivalently using focus!(r, [4,1,2])
-nactive(r)
-```
+# we want to find the digits with the first 5 qubits [1, 0, 1, 1, 1].
+evidense = [1, -2, 3, 4, 5]
+oracle_infer = inference_oracle(evidense)(nbit)
 
-Then we will have a check to above function
+niter = num_grover_step(prob_match_oracle(psi0, oracle_infer))
+gb_infer = chain(nbit, chain(oracle_infer, rb) for i = 1:niter)
+```
 
-```@example PhaseEstimation
-rand_unitary(N::Int) = qr(randn(N, N))[1]
+Now, let's start training
+```@example Grover
+for (i, blk) in enumerate(gb_infer)
+    apply!(psi0, blk)
+    p_target = prob_match_oracle(psi0, oracle_infer)
+    println("step $i, overlap^2 = $p_target")
+end
+```
 
-M = 5
-N = 3
+The result is
+```@example Grover
+pl = psi0 |> probs
+config = findn(pl.>0.5)[] - 1 |> bitarray(nbit)
+res = reshape(config, 5,3)
+```
 
-# prepair oracle matrix U
-V = rand_unitary(1<<N)
-phases = rand(1<<N)
-ϕ = Int(0b11101)/(1<<M)
-phases[3] = ϕ  # set the phase of the 3rd eigenstate manually.
-signs = exp.(2pi*im.*phases)
-U = V*Diagonal(signs)*V'  # notice U is unitary
+It is 2 ~
 
-# the state with phase ϕ
-psi = U[:,3]
+![infer](../assets/figures/digit2.png)
 
-res, reg = phase_estimation(zero_state(M), register(psi), GeneralMatrixGate(U))
-println("Phase is 2π * $(res[]), the exact value is 2π * $ϕ")
-```
+Congratuations! You get state of art quantum inference circuit!

examples/Grover.jl

@@ -1,13 +1,13 @@
 using Yao
-using Yao.Zoo: GroverIter, inference_oracle, prob_match_oracle
+using Yao.Zoo: groveriter!, inference_oracle, prob_match_oracle
 
 num_bit = 12
 oracle(reg::DefaultRegister) = (reg.state[100:101,:]*=-1; reg)
 target_state = zeros(1<<num_bit); target_state[100:101] = sqrt(0.5)
 
 # then solve the above problem
-it = GroverIter(oracle, uniform_state(num_bit))
-for (i, psi) in it
+it = groveriter!(uniform_state(num_bit), oracle)
+for (i, psi) in enumerate(it)
     overlap = abs(statevec(psi)'*target_state)
     println("step $(i-1), overlap = $overlap")
 end
@@ -23,8 +23,8 @@ Doing Inference, psi is the initial state, we want to search target space with s
 e.g. evidense [1, -3, 6] means the [1, 3, 6]-th bits take value [1, 0, 1].
 """
 oracle_infer = inference_oracle(evidense)(nqubits(psi))
-it = GroverIter(oracle_infer, psi)
-for (i, psi) in it
+it = groveriter!(psi, oracle_infer)
+for (i, psi) in enumerate(it)
     p_target = prob_match_oracle(psi, oracle_infer)
     println("step $(i-1), overlap^2 = $p_target")
 end

src/Blocks/ReflectBlock.jl

@@ -6,20 +6,21 @@ export ReflectBlock
 Householder reflection with respect to some target state, ``|\\psi\\rangle = 2|s\\rangle\\langle s|-1``.
 """
 struct ReflectBlock{N, T} <: PrimitiveBlock{N, T}
-    state :: Vector{T}
+    psi :: DefaultRegister{1, T}
 end
-ReflectBlock(state::Vector{T}) where T = ReflectBlock{log2i(length(state)), T}(state)
-ReflectBlock(psi::DefaultRegister) = ReflectBlock(statevec(psi))
+ReflectBlock(psi::DefaultRegister{1, T}) where T = ReflectBlock{nqubits(psi), T}(psi)
+ReflectBlock(state::Vector{T}) where T = ReflectBlock(register(state))
 
 function apply!(r::DefaultRegister, g::ReflectBlock)
-    r.state[:,:] .= 2.* (g.state'*r.state) .* reshape(g.state, :, 1) - r.state
+    v = g.psi |> state
+    r.state[:,:] .= 2 .* (v'*r.state) .* v - r.state
     r
 end
 
-==(A::ReflectBlock, B::ReflectBlock) = A.state == B.state
-copy(r::ReflectBlock) = ReflectBlock(r.state)
+==(A::ReflectBlock, B::ReflectBlock) = A.psi == B.psi
+copy(r::ReflectBlock) = ReflectBlock(r.psi)
 
-mat(r::ReflectBlock) = 2*r.state*r.state' - IMatrix(length(r.state))
+mat(r::ReflectBlock) = (v = r.psi |> statevec; 2*v*v' - IMatrix(length(v)))
 isreflexive(::ReflectBlock) = true
 ishermitian(::ReflectBlock) = true
 isunitary(::ReflectBlock) = true

src/Interfaces/Sequential.jl

@@ -12,6 +12,3 @@ function sequence end
 sequence() = Sequential([])
 sequence(blocks::AbstractBlock...) = Sequential(blocks...)
 sequence(blocks) = sequence(blocks...)
-
-# lazy constructors
-sequence(blocks...) = n->sequence([parse_block(n, each) for each in blocks])

src/Zoo/Differential.jl

@@ -1,4 +1,4 @@
-export diff_circuit, num_gradient, rotter, cnot_entangler, opgrad, collect_rotblocks
+export diff_circuit, num_gradient, rotter, cnot_entangler, opgrad, collect_rotblocks, perturb
 
 """
     rotter(noleading::Bool=false, notrailing::Bool=false) -> ChainBlock{1, ComplexF64}