diff --git a/lib/QuAlgorithmZoo/.github/workflows/TagBot.yml b/lib/QuAlgorithmZoo/.github/workflows/TagBot.yml new file mode 100644 index 000000000..d77d3a0c3 --- /dev/null +++ b/lib/QuAlgorithmZoo/.github/workflows/TagBot.yml @@ -0,0 +1,11 @@ +name: TagBot +on: + schedule: + - cron: 0 * * * * +jobs: + TagBot: + runs-on: ubuntu-latest + steps: + - uses: JuliaRegistries/TagBot@v1 + with: + token: ${{ secrets.GITHUB_TOKEN }} diff --git a/lib/QuAlgorithmZoo/.gitignore b/lib/QuAlgorithmZoo/.gitignore new file mode 100644 index 000000000..334a91acc --- /dev/null +++ b/lib/QuAlgorithmZoo/.gitignore @@ -0,0 +1,18 @@ +*.jl.cov +*.jl.*.cov +*.jl.mem + +build/ +docs/build/ +docs/site/ +docs/src/tutorial/ + +*.ipynb_checkpoints +**/*.ipynb_checkpoints +**/**/*.ipynb_checkpoints + +_*.dat +*.swp +__pycache__/ + +Manifest.toml diff --git a/lib/QuAlgorithmZoo/LICENSE.md b/lib/QuAlgorithmZoo/LICENSE.md new file mode 100644 index 000000000..a714ea550 --- /dev/null +++ b/lib/QuAlgorithmZoo/LICENSE.md @@ -0,0 +1,206 @@ +The QuAlgorithmZoo.jl package is licensed under the Apache License, Version 2.0: + +> Copyright (c) 2018: QuantumBFS. +> +> Apache License +> Version 2.0, January 2004 +> http://www.apache.org/licenses/ +> +> TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION +> +> 1. Definitions. +> +> "License" shall mean the terms and conditions for use, reproduction, +> and distribution as defined by Sections 1 through 9 of this document. +> +> "Licensor" shall mean the copyright owner or entity authorized by +> the copyright owner that is granting the License. +> +> "Legal Entity" shall mean the union of the acting entity and all +> other entities that control, are controlled by, or are under common +> control with that entity. For the purposes of this definition, +> "control" means (i) the power, direct or indirect, to cause the +> direction or management of such entity, whether by contract or +> otherwise, or (ii) ownership of fifty percent (50%) or more of the +> outstanding shares, or (iii) beneficial ownership of such entity. +> +> "You" (or "Your") shall mean an individual or Legal Entity +> exercising permissions granted by this License. +> +> "Source" form shall mean the preferred form for making modifications, +> including but not limited to software source code, documentation +> source, and configuration files. +> +> "Object" form shall mean any form resulting from mechanical +> transformation or translation of a Source form, including but +> not limited to compiled object code, generated documentation, +> and conversions to other media types. +> +> "Work" shall mean the work of authorship, whether in Source or +> Object form, made available under the License, as indicated by a +> copyright notice that is included in or attached to the work +> (an example is provided in the Appendix below). +> +> "Derivative Works" shall mean any work, whether in Source or Object +> form, that is based on (or derived from) the Work and for which the +> editorial revisions, annotations, elaborations, or other modifications +> represent, as a whole, an original work of authorship. For the purposes +> of this License, Derivative Works shall not include works that remain +> separable from, or merely link (or bind by name) to the interfaces of, +> the Work and Derivative Works thereof. +> +> "Contribution" shall mean any work of authorship, including +> the original version of the Work and any modifications or additions +> to that Work or Derivative Works thereof, that is intentionally +> submitted to Licensor for inclusion in the Work by the copyright owner +> or by an individual or Legal Entity authorized to submit on behalf of +> the copyright owner. For the purposes of this definition, "submitted" +> means any form of electronic, verbal, or written communication sent +> to the Licensor or its representatives, including but not limited to +> communication on electronic mailing lists, source code control systems, +> and issue tracking systems that are managed by, or on behalf of, the +> Licensor for the purpose of discussing and improving the Work, but +> excluding communication that is conspicuously marked or otherwise +> designated in writing by the copyright owner as "Not a Contribution." +> +> "Contributor" shall mean Licensor and any individual or Legal Entity +> on behalf of whom a Contribution has been received by Licensor and +> subsequently incorporated within the Work. +> +> 2. Grant of Copyright License. Subject to the terms and conditions of +> this License, each Contributor hereby grants to You a perpetual, +> worldwide, non-exclusive, no-charge, royalty-free, irrevocable +> copyright license to reproduce, prepare Derivative Works of, +> publicly display, publicly perform, sublicense, and distribute the +> Work and such Derivative Works in Source or Object form. +> +> 3. Grant of Patent License. Subject to the terms and conditions of +> this License, each Contributor hereby grants to You a perpetual, +> worldwide, non-exclusive, no-charge, royalty-free, irrevocable +> (except as stated in this section) patent license to make, have made, +> use, offer to sell, sell, import, and otherwise transfer the Work, +> where such license applies only to those patent claims licensable +> by such Contributor that are necessarily infringed by their +> Contribution(s) alone or by combination of their Contribution(s) +> with the Work to which such Contribution(s) was submitted. If You +> institute patent litigation against any entity (including a +> cross-claim or counterclaim in a lawsuit) alleging that the Work +> or a Contribution incorporated within the Work constitutes direct +> or contributory patent infringement, then any patent licenses +> granted to You under this License for that Work shall terminate +> as of the date such litigation is filed. +> +> 4. Redistribution. You may reproduce and distribute copies of the +> Work or Derivative Works thereof in any medium, with or without +> modifications, and in Source or Object form, provided that You +> meet the following conditions: +> +> (a) You must give any other recipients of the Work or +> Derivative Works a copy of this License; and +> +> (b) You must cause any modified files to carry prominent notices +> stating that You changed the files; and +> +> (c) You must retain, in the Source form of any Derivative Works +> that You distribute, all copyright, patent, trademark, and +> attribution notices from the Source form of the Work, +> excluding those notices that do not pertain to any part of +> the Derivative Works; and +> +> (d) If the Work includes a "NOTICE" text file as part of its +> distribution, then any Derivative Works that You distribute must +> include a readable copy of the attribution notices contained +> within such NOTICE file, excluding those notices that do not +> pertain to any part of the Derivative Works, in at least one +> of the following places: within a NOTICE text file distributed +> as part of the Derivative Works; within the Source form or +> documentation, if provided along with the Derivative Works; or, +> within a display generated by the Derivative Works, if and +> wherever such third-party notices normally appear. The contents +> of the NOTICE file are for informational purposes only and +> do not modify the License. You may add Your own attribution +> notices within Derivative Works that You distribute, alongside +> or as an addendum to the NOTICE text from the Work, provided +> that such additional attribution notices cannot be construed +> as modifying the License. +> +> You may add Your own copyright statement to Your modifications and +> may provide additional or different license terms and conditions +> for use, reproduction, or distribution of Your modifications, or +> for any such Derivative Works as a whole, provided Your use, +> reproduction, and distribution of the Work otherwise complies with +> the conditions stated in this License. +> +> 5. Submission of Contributions. Unless You explicitly state otherwise, +> any Contribution intentionally submitted for inclusion in the Work +> by You to the Licensor shall be under the terms and conditions of +> this License, without any additional terms or conditions. +> Notwithstanding the above, nothing herein shall supersede or modify +> the terms of any separate license agreement you may have executed +> with Licensor regarding such Contributions. +> +> 6. Trademarks. This License does not grant permission to use the trade +> names, trademarks, service marks, or product names of the Licensor, +> except as required for reasonable and customary use in describing the +> origin of the Work and reproducing the content of the NOTICE file. +> +> 7. Disclaimer of Warranty. Unless required by applicable law or +> agreed to in writing, Licensor provides the Work (and each +> Contributor provides its Contributions) on an "AS IS" BASIS, +> WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +> implied, including, without limitation, any warranties or conditions +> of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A +> PARTICULAR PURPOSE. You are solely responsible for determining the +> appropriateness of using or redistributing the Work and assume any +> risks associated with Your exercise of permissions under this License. +> +> 8. Limitation of Liability. In no event and under no legal theory, +> whether in tort (including negligence), contract, or otherwise, +> unless required by applicable law (such as deliberate and grossly +> negligent acts) or agreed to in writing, shall any Contributor be +> liable to You for damages, including any direct, indirect, special, +> incidental, or consequential damages of any character arising as a +> result of this License or out of the use or inability to use the +> Work (including but not limited to damages for loss of goodwill, +> work stoppage, computer failure or malfunction, or any and all +> other commercial damages or losses), even if such Contributor +> has been advised of the possibility of such damages. +> +> 9. Accepting Warranty or Additional Liability. While redistributing +> the Work or Derivative Works thereof, You may choose to offer, +> and charge a fee for, acceptance of support, warranty, indemnity, +> or other liability obligations and/or rights consistent with this +> License. However, in accepting such obligations, You may act only +> on Your own behalf and on Your sole responsibility, not on behalf +> of any other Contributor, and only if You agree to indemnify, +> defend, and hold each Contributor harmless for any liability +> incurred by, or claims asserted against, such Contributor by reason +> of your accepting any such warranty or additional liability. +> +> END OF TERMS AND CONDITIONS +> +> APPENDIX: How to apply the Apache License to your work. +> +> To apply the Apache License to your work, attach the following +> boilerplate notice, with the fields enclosed by brackets "{}" +> replaced with your own identifying information. (Don't include +> the brackets!) The text should be enclosed in the appropriate +> comment syntax for the file format. We also recommend that a +> file or class name and description of purpose be included on the +> same "printed page" as the copyright notice for easier +> identification within third-party archives. +> +> Copyright [year] [fullname] +> +> Licensed under the Apache License, Version 2.0 (the "License"); +> you may not use this file except in compliance with the License. +> You may obtain a copy of the License at +> +> http://www.apache.org/licenses/LICENSE-2.0 +> +> Unless required by applicable law or agreed to in writing, software +> distributed under the License is distributed on an "AS IS" BASIS, +> WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +> See the License for the specific language governing permissions and +> limitations under the License. +> diff --git a/lib/QuAlgorithmZoo/Project.toml b/lib/QuAlgorithmZoo/Project.toml new file mode 100644 index 000000000..db3d4ec44 --- /dev/null +++ b/lib/QuAlgorithmZoo/Project.toml @@ -0,0 +1,26 @@ +name = "QuAlgorithmZoo" +uuid = "65c24e16-9b0a-11e8-1353-efc5bc5f6586" +version = "0.1.0" + +[deps] +BitBasis = "50ba71b6-fa0f-514d-ae9a-0916efc90dcf" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +Yao = "5872b779-8223-5990-8dd0-5abbb0748c8c" +YaoArrayRegister = "e600142f-9330-5003-8abb-0ebd767abc51" +YaoBlocks = "418bc28f-b43b-5e0b-a6e7-61bbc1a2c1df" +YaoExtensions = "7a06699c-c960-11e9-3c98-9f78548b5f0f" + +[compat] +BitBasis = "0.6, 0.7" +Yao = "0.6" +YaoArrayRegister = "0.6, 0.7" +YaoBlocks = "0.8, 0.9, 0.10, 0.11" +YaoExtensions = "0.2" +julia = "1" + +[extras] +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[targets] +test = ["Test"] diff --git a/lib/QuAlgorithmZoo/README.md b/lib/QuAlgorithmZoo/README.md new file mode 100644 index 000000000..8c9ef9e4a --- /dev/null +++ b/lib/QuAlgorithmZoo/README.md @@ -0,0 +1,39 @@ +# QuAlgorithmZoo + +A curated implementation of quantum algorithms with [Yao.jl](https://github.com/QuantumBFS/Yao.jl) + +## Installation + +QuAlgorithmZoo.jl is not registered, please use the following command: + +```julia +pkg> dev https://github.com/QuantumBFS/QuAlgorithmZoo.jl.git +``` + +Then open directory `.julia/dev/QuAlgorithmZoo/examples` to find algorithms. + +## Contents + +- [QFT](https://github.com/QuantumBFS/YaoExtensions.jl) +- Phase Estimation +- Imaginary Time Evolution Quantum Eigensolver +- Variational Quantum Eigensolver +- Hadamard Test +- State Overlap Algorithms +- Quantum SVD + +In examples folder, you will find + +- HHL +- QAOA +- Quantum Circuit Born Machine +- QuGAN +- Shor +- Grover search + +- [QuODE](https://github.com/QuantumBFS/QuDiffEq.jl) +- [TensorNetwork Inspired Circuits](https://github.com/GiggleLiu/QuantumPEPS.jl) + +## License + +QuAlgorithmZoo.jl is released under Apache License 2.0. diff --git a/lib/QuAlgorithmZoo/docs/Project.toml b/lib/QuAlgorithmZoo/docs/Project.toml new file mode 100644 index 000000000..df7e87967 --- /dev/null +++ b/lib/QuAlgorithmZoo/docs/Project.toml @@ -0,0 +1,12 @@ +[deps] +BitBasis = "50ba71b6-fa0f-514d-ae9a-0916efc90dcf" +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" +Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c" +KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" +Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316" +Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" +Weave = "44d3d7a6-8a23-5bf8-98c5-b353f8df5ec9" +YaoArrayRegister = "e600142f-9330-5003-8abb-0ebd767abc51" +YaoBase = "a8f54c17-34bc-5a9d-b050-f522fe3f755f" +YaoBlocks = "418bc28f-b43b-5e0b-a6e7-61bbc1a2c1df" diff --git a/lib/QuAlgorithmZoo/docs/README.md b/lib/QuAlgorithmZoo/docs/README.md new file mode 100644 index 000000000..259000068 --- /dev/null +++ b/lib/QuAlgorithmZoo/docs/README.md @@ -0,0 +1,7 @@ +# QuAlgorithmZoo Documentation + +To build the documentation locally: + +```julia +julia make.jl local +``` diff --git a/lib/QuAlgorithmZoo/docs/make.jl b/lib/QuAlgorithmZoo/docs/make.jl new file mode 100644 index 000000000..862e9b4b1 --- /dev/null +++ b/lib/QuAlgorithmZoo/docs/make.jl @@ -0,0 +1,61 @@ +using Pkg + +pkg""" +add Documenter Literate Plots Yao +""" + +using Documenter, Literate, QuAlgorithmZoo + +const Examples = ["Grover", "VQE", "Shor"] + +const PATH = ( + tutorial = joinpath(@__DIR__, "src/tutorial"), + examples = joinpath(@__DIR__, "..", "examples") +) + +function process_literate_scripts() + TUTORIALS = [] + for token in Examples + file = "$token.jl" + filepath = joinpath(PATH.examples, token, file) + Literate.markdown(filepath, PATH.tutorial) + + filename, _ = splitext(file) + mdfile = join([filename, ".md"]) + # TODO: use PATH.tutorial rather then manual path + push!(TUTORIALS, relpath(joinpath("tutorial", mdfile))) + end + TUTORIALS +end + +#----------------------------------------------- + +function generate(islocal::Bool="local" in ARGS) + makedocs( + modules = [QuAlgorithmZoo], + clean = false, + format = :html, + sitename = "Quantum Algorithm Zoo", + linkcheck = !("skiplinks" in ARGS), + analytics = "UA-89508993-1", + pages = [ + "Home" => "index.md", + "Algorithms" => process_literate_scripts(), + "Manual" => Any[ + "man/zoo.md", + ], + ], + html_prettyurls = !islocal, + html_canonical = "https://quantumbfs.github.io/QuAlgorithmZoo.jl/latest/", + ) + + deploydocs( + repo = "github.com/QuantumBFS/QuAlgorithmZoo.jl.git", + target = "build", + julia = "1.0", + deps = nothing, + make = nothing, + ) +end + +generate(false) diff --git a/lib/QuAlgorithmZoo/docs/src/index.md b/lib/QuAlgorithmZoo/docs/src/index.md new file mode 100644 index 000000000..bacf06ca4 --- /dev/null +++ b/lib/QuAlgorithmZoo/docs/src/index.md @@ -0,0 +1,30 @@ +```@meta +CurrentModule = QuAlgorithmZoo +``` + +# Quantum Algorithm Zoo + +A curated implementation of quantum algorithms with [Yao.jl](https://github.com/QuantumBFS/Yao.jl) + +## Tutorial +```@contents +Pages = [ + "tutorial/Grover.md", + "tutorial/VQE.md", + "tutorial/Shor.md", +] +Depth = 1 +``` + +Some examples are not yet documented, please refer the [README](https://github.com/QuantumBFS/QuAlgorithmZoo.jl) for a full list of algorithms. + +## Manual + +The manual is a references of utilities defined in QuAlgorithmZoo, +like `Adam` optimizer and `NumberTheory` submodule. +```@contents +Pages = [ + "man/zoo.md", +] +Depth = 1 +``` diff --git a/lib/QuAlgorithmZoo/docs/src/man/zoo.md b/lib/QuAlgorithmZoo/docs/src/man/zoo.md new file mode 100644 index 000000000..672f4fe89 --- /dev/null +++ b/lib/QuAlgorithmZoo/docs/src/man/zoo.md @@ -0,0 +1,11 @@ +# QuAlgorithmZoo + +```@autodocs +Modules = [QuAlgorithmZoo] +Order = [:module, :constant, :type, :macro, :function] +``` + +```@autodocs +Modules = [QuAlgorithmZoo.NumberTheory] +Order = [:module, :constant, :type, :macro, :function] +``` diff --git a/lib/QuAlgorithmZoo/examples/GateLearning/gate_learning.jl b/lib/QuAlgorithmZoo/examples/GateLearning/gate_learning.jl new file mode 100644 index 000000000..bbc76fcb8 --- /dev/null +++ b/lib/QuAlgorithmZoo/examples/GateLearning/gate_learning.jl @@ -0,0 +1,27 @@ +using YaoExtensions, Yao +using Test, Random +using Optim: LBFGS, optimize +using Optim + +""" + learn_u4(u::AbstractMatrix; niter=100) + +Learn a general U4 gate. The optimizer is LBFGS. +""" +function learn_u4(u::AbstractBlock; niter=100) + ansatz = general_U4() + params = parameters(ansatz) + println("initial loss = $(operator_fidelity(u,ansatz))") + optimize(x->-operator_fidelity(u, dispatch!(ansatz, x)), + (G, x) -> (G .= -operator_fidelity'(u, dispatch!(ansatz, x))[2]), + parameters(ansatz), + LBFGS(), + Optim.Options(iterations=niter)) + println("final fidelity = $(operator_fidelity(u,ansatz))") + return ansatz +end + +using Random +Random.seed!(2) +u = matblock(rand_unitary(4)) +c = learn_u4(u; niter=150) diff --git a/lib/QuAlgorithmZoo/examples/Grover/Grover.jl b/lib/QuAlgorithmZoo/examples/Grover/Grover.jl new file mode 100644 index 000000000..852597eca --- /dev/null +++ b/lib/QuAlgorithmZoo/examples/Grover/Grover.jl @@ -0,0 +1,119 @@ +# # [Grover Search](@id Grover) +using Yao +using YaoExtensions: variational_circuit +using LinearAlgebra + +# ## Grover Step +# A single grover step is consist of applying oracle circuit and reflection circuit. +# The `reflection_circuit` function takes the wave function generator `U` as the input and returns `U|0><0|U'`. +function grover_step!(reg::AbstractRegister, oracle, U::AbstractBlock) + apply!(reg |> oracle, reflect_circuit(U)) +end + +function reflect_circuit(gen::AbstractBlock{N}) where N + reflect0 = control(N, -collect(1:N-1), N=>-Z) + chain(gen', reflect0, gen) +end + +# Compute the propotion of target states to estimate the number of iterations, +# which requires computing the output state. +function solution_state(oracle, gen::AbstractBlock{N}) where N + reg= zero_state(N) |> gen + reg.state[real.(statevec(ArrayReg(ones(ComplexF64, 1< oracle)) .> 0] .= 0 + normalize!(reg) +end + +function num_grover_step(oracle, gen::AbstractBlock{N}) where N + reg = zero_state(N) |> gen + ratio = abs2(solution_state(oracle, gen)'*reg) + Int(round(pi/4/sqrt(ratio)))-1 +end + +# #### Run +# First, we define the problem by an oracle, it finds bit string `bit"000001100100"`. +num_bit = 12 +oracle = matblock(Diagonal((v = ones(ComplexF64, 1< gen + +target_state = solution_state(oracle, gen) + +for i = 1:num_grover_step(oracle, gen) + grover_step!(reg, oracle, gen) + overlap = abs(reg'*target_state) + println("step $(i-1), overlap = $overlap") +end + +# ## Rejection Sampling +# In practise, it is often not possible to determine the number of iterations before actual running. +# we can use rejection sampling technique to avoid estimating the number of grover steps. + +using Random; Random.seed!(2) #src + +# In a single try, we `apply` the grover algorithm for `nstep` times. +function single_try(oracle, gen::AbstractBlock{N}, nstep::Int; nbatch::Int) where N + reg = zero_state(N+1; nbatch=nshot) + focus!(reg, 1:N) do r + r |> gen + for i = 1:nstep + grover_step!(r, oracle, gen) + end + return r + end + reg |> checker + res = measure!(RemoveMeasured(), reg, (N+1)) + return res, reg +end + +# After running the grover search, we have a checker program that flips the ancilla qubit +# if the output is the desired value, we assume the checker program can be implemented in polynomial time. +# to gaurante the output is correct. +# We contruct a checker "program", if the result is correct, flip the ancilla qubit +ctrl = -collect(1:num_bit); ctrl[[3,6,7]] *= -1 +checker = control(num_bit+1,ctrl, num_bit+1=>X) + +# The register is batched, with batch dimension `nshot`. +# [`focus!`](@ref Yao.focus!) views the first 1-N qubts as system. +# For a batched register, [`measure!`](@ref Yao.measure!) +# returns a vector of bitstring as output. + +# #### Run +maxtry = 100 +nshot = 3 + +for nstep = 0:maxtry + println("number of iter = $nstep") + res, reg = single_try(oracle, gen, nstep; nbatch=3) + + ## success! + if any(==(1), res) + overlap_final = viewbatch(reg, findfirst(==(1), res))'*target_state + println("success, overlap = $(overlap_final)") + break + end +end + +# The final state has an overlap of `1` with the target state. + +# ## Amplitude Amplification +# Given a circuit to generate a state, +# now we want to project out the subspace with [1,3,5,8,9,11,12] fixed to 1 and [4,6] fixed to 0. +# We can construct an oracle +evidense = [1, 3, -4, 5, -6, 8, 9, 11, 12] +function inference_oracle(nbit::Int, locs::Vector{Int}) + control(nbit, locs[1:end-1], abs(locs[end]) => (locs[end]>0 ? Z : -Z)) +end +oracle = inference_oracle(nqubits(reg), evidense) + +# We use a variational circuit generator defined in `YaoExtensions` +gen = dispatch!(variational_circuit(num_bit), :random) +reg = zero_state(num_bit) |> gen + +# #### Run +solution = solution_state(oracle, gen) +for i = 1:num_grover_step(oracle, gen) + grover_step!(reg, oracle, gen) + println("step $(i-1), overlap = $(abs(reg'*solution))") +end diff --git a/lib/QuAlgorithmZoo/examples/Grover/tests.jl b/lib/QuAlgorithmZoo/examples/Grover/tests.jl new file mode 100644 index 000000000..a9ba1e048 --- /dev/null +++ b/lib/QuAlgorithmZoo/examples/Grover/tests.jl @@ -0,0 +1,67 @@ +include("Grover.jl") +using Test, BitBasis + +"""traditional grover search algorithm.""" +function grover_search(oracle::AbstractBlock{N}, gen::AbstractBlock{N}=repeat(N,H,1:N)) where N + reg = zero_state(N) |> gen + for i = 1:num_grover_step(oracle, gen) + grover_step!(reg, oracle, gen) + end + return reg +end + +function grover_circuit(oracle::AbstractBlock{N}, gen::AbstractBlock{N}, niter::Int=num_grover_step(oracle, gen)) where {N} + chain(N, chain(oracle, reflect_circuit(gen)) for i = 1:niter) +end + +#################### Tests ################## + +@testset "oracle" begin + oracle = inference_oracle(3, [2,-1,3]) + # ≈ method for Identity/PermuteMultiply/Sparse + # add and mimus between sparse matrices. + # alway use sorted CSC format. + v = ones(1<<3) + v[Int(0b110)+1] *= -1 + @test mat(oracle) ≈ Diagonal(v) + @test mat(chain(phase(π), Z)) ≈ -mat(Z) +end + +@testset "Grover Search" begin + ####### Construct Grover Search Using Reflection Block + num_bit = 12 + oracle = inference_oracle(num_bit, push!(collect(Int, 1:num_bit-1), num_bit)) + + psi = grover_search(oracle) + target_state = zeros(1< gen, gb) == grover_search(or) +end + +@testset "test inference" begin + Random.seed!(2) + num_bit = 12 + gen = dispatch!(variational_circuit(num_bit), :random) + psi0 = zero_state(num_bit) |> gen + #psi0 = uniform_state(num_bit) + evidense = [1, 2, 3, 4, 5, 6, 7, 8, 9] + #evidense = collect(1:num_bit) + + # the desired subspace + basis = collect(UInt, 0:1<0))...] + + v_desired = statevec(psi0)[subinds .+ 1] + p = norm(v_desired)^2 + v_desired[:] ./= sqrt(p) + + # search the subspace + psi = grover_search(inference_oracle(num_bit, evidense), gen) + @test isapprox((psi.state[subinds .+ 1]'*v_desired) |> abs2, 1, atol=3e-2) +end diff --git a/lib/QuAlgorithmZoo/examples/HHL/HHL.jl b/lib/QuAlgorithmZoo/examples/HHL/HHL.jl new file mode 100644 index 000000000..2c4901389 --- /dev/null +++ b/lib/QuAlgorithmZoo/examples/HHL/HHL.jl @@ -0,0 +1,90 @@ +using Yao, YaoBlocks +using BitBasis +using YaoArrayRegister +using Test, LinearAlgebra +using QuAlgorithmZoo: PEBlock + +include("HHLlib.jl") + +function crot(n_reg::Int, C_value::Real) + n_rot = n_reg + 1 + rot = chain(n_rot) + θ = zeros(1< 1 + return println("C_value = $C_value, λ = $λ, sinθ = $sin_value > 1, please lower C_value.\n") + end + θ[(i)] = 2.0*asin(C_value / λ) + push!(rot, control(c_bit, 1=>Ry(θ[(i)]))) + end + rot +end + +@testset "HHLCRot" begin + hr = HHLCRot{4}([4,3,2], 1, 0.01) + reg = rand_state(4) + @test reg |> copy |> hr |> isnormalized + + hr2 = crot(3, 0.01) + reg1 = reg |> copy |> hr + reg2 = reg |> copy |> hr2 + @test fidelity(reg1, reg2)[] ≈ 1 +end + +""" + hhl_problem(nbit::Int) -> Tuple + +Returns (A, b), where + * `A` is positive definite, hermitian, with its maximum eigenvalue λ_max < 1. + * `b` is normalized. +""" +function hhl_problem(nbit::Int) + siz = 1< = |b>. + ## signs: Diagonal Matrix of eigen values of A. + ## base_space: the eigen space of A. + ## x: |x>. + using Random + Random.seed!(2) + N = 3 + A, b = hhl_problem(N) + x = A^(-1)*b # base_i = base_space[:,i] ϕ1 = (A*base_i./base_i)[1] + + ## n_b : number of bits for |b>. + ## n_reg: number of PE register. + ## n_all: number of all bits. + n_reg = 12 + + ## C_value: value of constant C in control rotation. + ## It should be samller than the minimum eigen value of A. + C_value = minimum(eigvals(A) .|> abs)*0.25 + #C_value = 1.0/(1<|00>|u>. + @test isapprox.(x, res, atol=0.5) |> all +end diff --git a/lib/QuAlgorithmZoo/examples/HHL/HHLlib.jl b/lib/QuAlgorithmZoo/examples/HHL/HHLlib.jl new file mode 100644 index 000000000..1ff8146a0 --- /dev/null +++ b/lib/QuAlgorithmZoo/examples/HHL/HHLlib.jl @@ -0,0 +1,84 @@ +export hhlcircuit, hhlproject!, hhlsolve, HHLCRot + +""" + HHLCRot{N, NC} <: PrimitiveBlock{N} + +Controlled rotation gate used in HHL algorithm, applied on N qubits. + + * cbits: control bits. + * ibit:: the ancilla bit. + * C_value:: the value of constant "C", should be smaller than the spectrum "gap". +""" +struct HHLCRot{N, NC, T} <: PrimitiveBlock{N} + cbits::Vector{Int} + ibit::Int + C_value::T + HHLCRot{N}(cbits::Vector{Int}, ibit::Int, C_value::T) where {N, T} = new{N, length(cbits), T}(cbits, ibit, C_value) +end + +@inline function hhlrotmat(λ::Real, C_value::Real) + b = C_value/λ + a = sqrt(1-b^2) + a, -b, b, a +end + +function YaoBlocks._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(readbit(i-1, hr.cbits...), nbits=nbit-1) + if λ >= hr.C_value + u = hhlrotmat(λ, hr.C_value) + YaoArrayRegister.u1rows!(state(reg), i, i+step, u...) + end + end + end + reg +end + +""" + hhlproject!(all_bit::ArrayReg, n_reg::Int) -> Vector + +project to aiming state |1>|00>|u>, and return |u> vector. +""" +function hhlproject!(all_bit::ArrayReg, n_reg::Int) + all_bit |> focus!(1:(n_reg+1)...) |> select!(1) |> state |> vec +end + +""" +Function to build up a HHL circuit. +""" +function hhlcircuit(UG, n_reg::Int, C_value::Real) + n_b = nqubits(UG) + n_all = 1 + n_reg + n_b + pe = PEBlock(UG, n_reg, n_b) + cr = HHLCRot{n_reg+1}([2:n_reg+1...], 1, C_value) + chain(n_all, subroutine(n_all, pe, [2:n_all...,]), subroutine(n_all, cr, [1:(n_reg+1)...,]), subroutine(n_all, pe', [2:n_all...,])) +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 = matblock(exp(2π*im.*A)) + + # Generating input bits + all_bit = join(ArrayReg(b), zero_state(n_reg), zero_state(1)) + + # Construct HHL circuit. + circuit = hhlcircuit(UG, n_reg, C_value) + + # Apply bits to the circuit. + apply!(all_bit, circuit) + + # Get state of aiming state |1>|00>|u>. + hhlproject!(all_bit, n_reg) ./ C_value +end diff --git a/lib/QuAlgorithmZoo/examples/HHL/lin_diffEq_HHL.jl b/lib/QuAlgorithmZoo/examples/HHL/lin_diffEq_HHL.jl new file mode 100644 index 000000000..11a7af89a --- /dev/null +++ b/lib/QuAlgorithmZoo/examples/HHL/lin_diffEq_HHL.jl @@ -0,0 +1,52 @@ +using Yao +using BitBasis +using Random +using QuAlgorithmZoo +using Test, LinearAlgebra +using OrdinaryDiffEq + +function diffeq_problem(nbit::Int) + siz = 1< all + + res = solve(qprob, QuLeapfrog(), h, n_reg) + r = res[(N_t + 1)*2 + 2^N - 1: (N_t + 1)*2 + 2^N + N_t - 3] # range of relevant values in the obtained state. + @test isapprox.(s, r, atol = 0.3) |> all + + res = solve(qprob, QuAB2(), h,n_reg) + r = res[(N_t + 1)*2 + 2^N - 1: (N_t + 1)*2 + 2^N + N_t - 3] # range of relevant values in the obtained state. + @test isapprox.(s, r, atol = 0.3) |> all + + res = solve(qprob, QuAB3(), h,n_reg) + r = res[(N_t + 1)*2 + 2^N - 1: (N_t + 1)*2 + 2^N + N_t - 3] # range of relevant values in the obtained state. + @test isapprox.(s, r, atol = 0.3) |> all + + res = solve(qprob, QuAB4(), h ,n_reg) + r = res[(N_t + 1)*2 + 2^N - 1: (N_t + 1)*2 + 2^N + N_t - 3] # range of relevant values in the obtained state. + @test isapprox.(s, r, atol = 0.3) |> all +end; diff --git a/lib/QuAlgorithmZoo/examples/HHL/lin_diffEq_HHLlib.jl b/lib/QuAlgorithmZoo/examples/HHL/lin_diffEq_HHLlib.jl new file mode 100644 index 000000000..3e57b6a81 --- /dev/null +++ b/lib/QuAlgorithmZoo/examples/HHL/lin_diffEq_HHLlib.jl @@ -0,0 +1,157 @@ +export array_qudiff, prepare_init_state, LDEMSAlgHHL, bval, aval +export QuEuler, QuLeapfrog, QuAB2, QuAB3, QuAB4 +export QuLDEMSProblem + +using DiffEqBase +""" + Based on : arxiv.org/abs/1010.2745v2 + + * array_qudiff(N_t,N,h,A) - generates matrix for k-step solver + * prepare_init_state(b,x,h,N_t) - generates inital states + * solve_qudiff - solver + + x' = Ax + b + + * A - input matrix. + * b - input vector. + * x - inital vector + * N - dimension of b (as a power of 2). + * h - step size. + * tspan - time span. +""" + +""" + LDEMSAlgHHL + * step - step for multistep method + * α - coefficients for xₙ + * β - coefficent for xₙ' +""" +abstract type QuODEAlgorithm <: DiffEqBase.AbstractODEAlgorithm end +#abstract type QuODEProblem{uType,tType,isinplace} <: DiffEqBase.AbstractODEProblem{uType,tType,isinplace} end +abstract type LDEMSAlgHHL <: QuODEAlgorithm end + +struct QuLDEMSProblem{F,C,U,T} #<: QuODEProblem{uType,tType,isinplace} + A::F + b::C + u0::U + tspan::NTuple{2,T} + + #function QuLDEMSProblem(A,b,u0,tspan) + # new{typeof(u0),typeof(tspan),false,typeof(A),typeof(b)}(A,b,u0,tspan) + #end +end + +""" + Explicit Linear Multistep Methods +""" +struct QuEuler{T}<:LDEMSAlgHHL + step::Int + α::Vector{T} + β::Vector{T} + + QuEuler(::Type{T} = Float64) where {T} = new{T}(1,[1.0,],[1.0,]) +end + +struct QuLeapfrog{T}<:LDEMSAlgHHL + step::Int + α::Vector{T} + β::Vector{T} + + QuLeapfrog(::Type{T} = Float64) where {T} = new{T}(2,[0, 1.0],[2.0, 0]) +end +struct QuAB2{T}<:LDEMSAlgHHL + step::Int + α::Vector{T} + β::Vector{T} + + QuAB2(::Type{T} = Float64) where {T} = new{T}(2,[1.0, 0], [1.5, -0.5]) +end +struct QuAB3{T}<:LDEMSAlgHHL + step::Int + α::Vector{T} + β::Vector{T} + QuAB3(::Type{T} = Float64) where {T} = new{T}(3,[1.0, 0, 0], [23/12, -16/12, 5/12]) +end +struct QuAB4{T}<:LDEMSAlgHHL + step::Int + α::Vector{T} + β::Vector{T} + + QuAB4(::Type{T} = Float64) where {T} = new{T}(4,[1.0, 0, 0, 0], [55/24, -59/24, 37/24, -9/24]) +end + +function bval(alg::LDEMSAlgHHL,t,h,g::Function) + b = zero(g(1)) + for i in 1:(alg.step) + b += alg.β[i]*g(t-(i-1)*h) + end + return b +end + +function aval(alg::LDEMSAlgHHL,t,h,g::Function) + sz, = size(g(1)) + A = Array{ComplexF64}(undef,sz,(alg.step + 1)*sz) + i_mat = Matrix{Float64}(I, size(g(1))) + A[1:sz,sz*(alg.step) + 1:sz*(alg.step + 1)] = i_mat + for i in 1:alg.step + A[1:sz,sz*(i - 1) + 1: sz*i] = -1*(alg.α[alg.step - i + 1]*i_mat + h*alg.β[alg.step - i + 1]*g(t - (alg.step - i)*h)) + end + return A +end + +function prepare_init_state(tspan::NTuple{2, Float64},x::Vector,h::Float64,g::Function,alg::LDEMSAlgHHL) + N_t = round(Int, (tspan[2] - tspan[1])/h + 1) #number of time steps + N = nextpow(2,2*N_t + 1) # To ensure we have a power of 2 dimension for matrix + sz, = size(g(1)) + init_state = zeros(ComplexF64,2*(N)*sz) + #inital value + init_state[1:sz] = x + for i in 2:N_t + b = bval(alg,h*(i - 1) + tspan[1],h,g) + init_state[Int(sz*(i - 1) + 1):Int(sz*(i))] = h*b + end + return init_state +end + +function array_qudiff(tspan::NTuple{2, Float64},h::Float64,g::Function,alg::LDEMSAlgHHL) + sz, = size(g(1)) + i_mat = Matrix{Float64}(I, size(g(1))) + N_t = round(Int, (tspan[2] - tspan[1])/h + 1) #number of time steps + N = nextpow(2,2*N_t + 1) # To ensure we have a power of 2 dimension for matrix + A_ = zeros(ComplexF64, N*sz, N*sz) + # Generates First two rows + @inbounds A_[1:sz, 1:sz] = i_mat + @inbounds A_[sz + 1:2*sz, 1:sz] = -1*(i_mat + h*g(tspan[1])) + @inbounds A_[sz + 1:2*sz,sz+1:sz*2] = i_mat + #Generates additional rows based on k - step + for i in 3:alg.step + @inbounds A_[sz*(i - 1) + 1:sz*i, sz*(i - 3) + 1:sz*i] = aval(QuAB2(),(i-2)*h + tspan[1],h,g) + end + for i in alg.step + 1:N_t + @inbounds A_[sz*(i - 1) + 1:sz*(i), sz*(i - alg.step - 1) + 1:sz*i] = aval(alg,(i - 2)*h + tspan[1],h,g) + end + #Generates half mirroring matrix + for i in N_t + 1:N + @inbounds A_[sz*(i - 1) + 1:sz*(i), sz*(i - 2) + 1:sz*(i - 1)] = -1*i_mat + @inbounds A_[sz*(i - 1) + 1:sz*(i), sz*(i - 1) + 1:sz*i] = i_mat + end + A_ = [zero(A_) A_;A_' zero(A_)] + return A_ +end + +function DiffEqBase.solve(prob::QuLDEMSProblem{F,C,U,T}, alg::LDEMSAlgHHL, dt = (prob.tspan[2]-prob.tspan[1])/100, n_reg::Int = 12) where {F,C,U,T} + A = prob.A + b = prob.b + tspan = prob.tspan + x = prob.u0 + + mat = array_qudiff(tspan, dt, A, alg) + state = prepare_init_state(tspan, x, dt, b, alg) + λ = maximum(eigvals(mat)) + C_value = minimum(eigvals(mat) .|> abs)*0.01; + mat = 1/(λ*2)*mat + state = state*1/(2*λ) |> normalize! + res = hhlsolve(mat,state, n_reg, C_value) + res = res/λ + return res +end; diff --git a/lib/QuAlgorithmZoo/examples/PortQuantumInfromation/yao_port_qi.ipynb b/lib/QuAlgorithmZoo/examples/PortQuantumInfromation/yao_port_qi.ipynb new file mode 100644 index 000000000..3c3a438d4 --- /dev/null +++ b/lib/QuAlgorithmZoo/examples/PortQuantumInfromation/yao_port_qi.ipynb @@ -0,0 +1,374 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Porting Yao.jl with QuantumInformation.jl\n", + "### GiggleLiu" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# overview\n", + "\n", + " [`Yao`](https://github.com/QuantumBFS/Yao.jl) is a powerful tool for quantum circuit based simulation, but it does not support many density matrix related operations. This is why we need to port `Yao.jl` with [`QuantumInformation (QI)`](https://github.com/QuantumBFS/QuantumInformation.jl) sometimes (e.g. for computing entanglement entropy).\n", + " \n", + "* `Yao.jl` Documentation: https://quantumbfs.github.io/Yao.jl/latest/ (paper is comming out)\n", + "* `QuantumInformation.jl` paper: https://arxiv.org/abs/1806.11464\n", + " \n", + "### `Yao` provides\n", + "* high performance quantum circuit based simulation\n", + " * parameter management\n", + " * gradients\n", + " * batched regiser\n", + "* operator matrix representation and arithmatics\n", + "* [quantum algorithms](https://github.com/QuantumBFS/QuAlgorithmZoo.jl)\n", + "* [GPU support](https://github.com/QuantumBFS/CuYao.jl)\n", + "\n", + "### `QI` provides\n", + "\n", + "* Compute entropy from density matrices\n", + "* Quantum channels, four types of channel representations\n", + " * Kraus Operator\n", + " * Super operator\n", + " * Dynamic matrices\n", + " * Stinespring representation\n", + "* Compute norm, distance and distingushability between \"states\" (density matrices)\n", + " * Hilbert–Schmidt norm and distance\n", + " * trace norm and *distance*\n", + " * diamond norm\n", + " * Bures distane and Bures angles\n", + " * *fidelity* and superfidelity\n", + " * KL-divergence\n", + " * JS-distance\n", + "* Compute the amount of entanglement\n", + " * negativity\n", + " * positive partial trace\n", + " * concurrence\n", + "* POVM measurements" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [], + "source": [ + "import Yao\n", + "using Yao: ArrayReg, ρ, mat, ConstGate, purify, exchange_sysenv, @bit_str, statevec\n", + "import QuantumInformation; const QI = QuantumInformation\n", + "using QuantumInformation: ket\n", + "using LinearAlgebra\n", + "using Test" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Obtain reduced density matrices in Yao\n", + "-------------------------\n", + "The memory layout of `Yao` register and `QI` ket are similar, their basis are both [little endian](https://en.wikipedia.org/wiki/Endianness), despite they have different representation powers\n", + "\n", + "* `Yao` support batch,\n", + "* `QI` is not limited to qubits.\n", + "\n", + "\n", + "`Yao` does not have much operations defined on density matrices, but purified states with environment,\n", + "On the other side, most operations in `QI` are defined on **(density) matrices**, they can be easily obtained in `Yao`." + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "\u001b[32m\u001b[1mTest Passed\u001b[22m\u001b[39m" + ] + }, + "execution_count": 17, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# construct a product state, notice the indexing in `QI` starts from `1`\n", + "@test QI.ket(3, 1<<4) ≈ statevec(ArrayReg(bit\"0010\"))\n", + "\n", + "# join two registers, notice little endian convension is used here.\n", + "reg = join(ArrayReg(bit\"10\"), ArrayReg(bit\"11\"))\n", + "v = QI.:⊗(QI.ket(0b10+1,1<<2), QI.ket(0b11+1,1<<2))\n", + "@test statevec(reg) ≈ v" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "data": { + "text/plain": [ + "16×16 Array{Complex{Float64},2}:\n", + " 0.0668399+0.0im … 0.0048149+0.00800254im \n", + " -0.00683079+0.00430075im 0.00271044+0.013467im \n", + " -0.00405524-0.00233655im 0.00489161-0.00506099im \n", + " 0.0041184-0.00690317im -0.00724508+0.00433365im \n", + " 0.000248112-0.00614303im -0.00169715-0.0060107im \n", + " -0.00638715+0.00343611im … -0.00346919+0.0104737im \n", + " -0.0032589-0.00594789im -0.00502371+0.00889227im \n", + " 0.0053714-0.00448422im 0.000149836+0.00490488im \n", + " -0.00485418+0.00190183im -0.00707738-0.0117206im \n", + " -0.00185245-0.0113168im -0.00100021+0.00456715im \n", + " 0.000202351+0.00648573im … 9.29962e-5-0.00362312im \n", + " 0.0038004-0.00408768im 0.00290617+0.0109155im \n", + " -0.00488166-0.00699333im -0.00471523+0.000137239im\n", + " 0.00485705+0.00532262im 0.00956895-0.00457732im \n", + " 0.00756613-0.00569826im 0.0032851+0.0014402im \n", + " 0.0048149-0.00800254im … 0.0786938+0.0im " + ] + }, + "execution_count": 18, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# convert a Yao register to density matrix in QI\n", + "reg2dm(reg::ArrayReg{1}) = reg |> ρ |> Matrix\n", + "\n", + "# e.g. obtain a reduced denstiy matrix for subsystem 1,2,3,4\n", + "reg = Yao.rand_state(10)\n", + "freg = Yao.focus!(reg, 1:4) # make qubits 1-4 active\n", + "reg2dm(freg)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "One can also convert a density matrix to a a quantum state through **purification**" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "\u001b[32m\u001b[1mTest Passed\u001b[22m\u001b[39m" + ] + }, + "execution_count": 19, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# e.g. purify a state and recover it\n", + "reg = Yao.rand_state(6) |> Yao.focus!(1:4)\n", + "reg_p = purify(reg |> ρ; nbit_env=2)\n", + "@test Yao.fidelity(reg, reg_p)[] ≈ 1" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "entanglement & state distance\n", + "----------------\n" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "\u001b[32m\u001b[1mTest Passed\u001b[22m\u001b[39m" + ] + }, + "execution_count": 20, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "reg1 = Yao.rand_state(10)\n", + "freg1 = Yao.focus!(reg1, 1:4)\n", + "reg2 = Yao.rand_state(6)\n", + "freg2 = Yao.focus!(reg2, 1:4)\n", + "dm1, dm2 = freg1 |> reg2dm, freg2 |> reg2dm\n", + "\n", + "# trace distance between two registers (different by a factor 2)\n", + "@test Yao.tracedist(freg1, freg2)[]/2 ≈ QI.trace_distance(dm1, dm2)" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "QI.vonneumann_entropy(dm1) = 2.6568839293081608\n", + "QI.vonneumann_entropy(dm2) = 1.3245543916726097\n" + ] + }, + { + "data": { + "text/plain": [ + "1.3245543916726097" + ] + }, + "execution_count": 21, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# get the entanglement entropy between system and env\n", + "@show QI.vonneumann_entropy(dm1)\n", + "@show QI.vonneumann_entropy(dm2)" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "1.5694621854109723" + ] + }, + "execution_count": 22, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# KL-divergence (or relative entropy)\n", + "QI.kl_divergence(dm2, dm1)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Note: you can defined many distances and entropies in a similar way, we don't enumerate it." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Quantum Operations/Quantum Gates\n", + "------------------------\n", + "\n", + "A quantum gate in `Yao` is equivalent to a unitary channel in `QI`, matrix representations of blocks in `Yao` can be used to construct channels." + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "QuantumInformation.KrausOperators{Array{Complex{Float64},2}}\n", + " dimensions: (2, 2)\n", + " Complex{Float64}[1.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im]\n", + " Complex{Float64}[0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 1.0 + 0.0im]\n", + " Complex{Float64}[0.0 + 0.0im 1.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im]" + ] + }, + "execution_count": 23, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# construct a Kraus Operator\n", + "QI.KrausOperators([Matrix(ConstGate.P0), Matrix(ConstGate.P1), Matrix(ConstGate.Pu)])" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + ":(#= In[24]:9 =# @test (copy(reg) |> Yao.chain(b1, b2)) |> reg2dm ≈ (c2 ∘ c1)(reg |> reg2dm))" + ] + }, + "execution_count": 24, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# applying a rotation gate\n", + "b1 = Yao.put(2,2=>Yao.Rx(0.3π))\n", + "c1 = QI.UnitaryChannel(mat(b1))\n", + "b2 = Yao.put(2,2=>Yao.Ry(0.3π))\n", + "c2 = QI.UnitaryChannel(mat(b2))\n", + "\n", + "reg = Yao.rand_state(2)\n", + "@test copy(reg) |> b1 |> reg2dm ≈ c1(reg |> reg2dm)\n", + ":@test copy(reg) |> Yao.chain(b1,b2) |> reg2dm ≈ (c2∘c1)(reg |> reg2dm)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "@webio": { + "lastCommId": null, + "lastKernelId": null + }, + "kernelspec": { + "display_name": "Julia 1.3.0", + "language": "julia", + "name": "julia-1.3" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.3.0" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/lib/QuAlgorithmZoo/examples/PortZygote/chainrules_patch.jl b/lib/QuAlgorithmZoo/examples/PortZygote/chainrules_patch.jl new file mode 100644 index 000000000..76b706fc3 --- /dev/null +++ b/lib/QuAlgorithmZoo/examples/PortZygote/chainrules_patch.jl @@ -0,0 +1,66 @@ +import ChainRulesCore: rrule, @non_differentiable, NoTangent +using Yao, Yao.AD + +function rrule(::typeof(apply), reg::ArrayReg, block::AbstractBlock) + out = apply(reg, block) + out, function (outδ) + (in, inδ), paramsδ = apply_back((copy(out), outδ), block) + return (NoTangent(), inδ, paramsδ) + end +end + +function rrule(::typeof(dispatch), block::AbstractBlock, params) + out = dispatch(block, params) + out, function (outδ) + (NoTangent(), NoTangent(), outδ) + end +end + +function rrule(::typeof(expect), op::AbstractBlock, reg::AbstractRegister{B}) where {B} + out = expect(op, reg) + out, function (outδ) + greg = Yao.AD.expect_g(op, reg) + for b=1:B + viewbatch(greg, b).state .*= 2*outδ[b] + end + return (NoTangent(), NoTangent(), greg) + end +end + +function rrule(::Type{Matrix}, block::AbstractBlock) + out = Matrix(block) + out, function (outδ) + paramsδ = mat_back(block, outδ) + return (NoTangent(), paramsδ) + end +end + +function rrule(::Type{ArrayReg{B}}, raw::AbstractArray) where B + ArrayReg{B}(raw), adjy->(NoTangent(), reshape(adjy.state, size(raw))) +end + +function rrule(::Type{ArrayReg{B}}, raw::ArrayReg) where B + ArrayReg{B}(raw), adjy->(NoTangent(), adjy) +end + +function rrule(::Type{ArrayReg}, raw::AbstractArray) + ArrayReg(raw), adjy->(NoTangent(), reshape(adjy.state, size(raw))) +end + +function rrule(::typeof(copy), reg::ArrayReg) where B + copy(reg), adjy->(NoTangent(), adjy) +end + +_totype(::Type{T}, x::AbstractArray{T}) where T = x +_totype(::Type{T}, x::AbstractArray{T2}) where {T,T2} = convert.(T, x) +rrule(::typeof(state), reg::ArrayReg{B,T}) where {B,T} = state(reg), adjy->(NoTangent(), ArrayReg(_totype(T, adjy))) +rrule(::typeof(statevec), reg::ArrayReg{B,T}) where {B,T} = statevec(reg), adjy->(NoTangent(), ArrayReg(_totype(T, adjy))) +rrule(::typeof(state), reg::AdjointArrayReg{B,T}) where {B,T} = state(reg), adjy->(NoTangent(), ArrayReg(_totype(T, adjy)')') +rrule(::typeof(statevec), reg::AdjointArrayReg{B,T}) where {B,T} = statevec(reg), adjy->(NoTangent(), ArrayReg(_totype(adjy)')') +rrule(::typeof(parent), reg::AdjointArrayReg) = parent(reg), adjy->(NoTangent(), adjy') +rrule(::typeof(Base.adjoint), reg::ArrayReg) = Base.adjoint(reg), adjy->(NoTangent(), parent(adjy)) +@non_differentiable Yao.nparameters(::Any) +@non_differentiable Yao.zero_state(args...) +@non_differentiable Yao.rand_state(args...) +@non_differentiable Yao.uniform_state(args...) +@non_differentiable Yao.product_state(args...) diff --git a/lib/QuAlgorithmZoo/examples/PortZygote/gate_learning.jl b/lib/QuAlgorithmZoo/examples/PortZygote/gate_learning.jl new file mode 100644 index 000000000..2deb7fbf1 --- /dev/null +++ b/lib/QuAlgorithmZoo/examples/PortZygote/gate_learning.jl @@ -0,0 +1,33 @@ +using YaoExtensions, Yao +using Test, Random +using Optim: LBFGS, optimize +using Optim + +# port the `Matrix` function to Yao's AD. +using Zygote +include("chainrules_patch.jl") + +function loss(u, ansatz) + m = Matrix(ansatz) + sum(abs.(u .- m)) +end + +""" + learn_u4(u::AbstractMatrix; niter=100) + +Learn a general U4 gate. The optimizer is LBFGS. +""" +function learn_u4(u::AbstractMatrix; niter=100) + ansatz = general_U4() * put(2, 1=>phase(0.0)) # initial values are 0, here, we attach a global phase. + params = parameters(ansatz) + g!(G, x) = (ansatz=dispatch(ansatz, x); G .= Zygote.gradient(ansatz->loss(u, ansatz), ansatz)[1]) + optimize(x->(ansatz=dispatch(ansatz, x); loss(u, ansatz)), g!, parameters(ansatz), + LBFGS(), Optim.Options(iterations=niter)) + println("final loss = $(loss(u,ansatz))") + return ansatz +end + +using Random +Random.seed!(3) +u = rand_unitary(4) +c = learn_u4(u; niter=150) diff --git a/lib/QuAlgorithmZoo/examples/PortZygote/shared_parameters.jl b/lib/QuAlgorithmZoo/examples/PortZygote/shared_parameters.jl new file mode 100644 index 000000000..02602d824 --- /dev/null +++ b/lib/QuAlgorithmZoo/examples/PortZygote/shared_parameters.jl @@ -0,0 +1,40 @@ +using Zygote +include("chainrules_patch.jl") + +import YaoExtensions, Random + +c = YaoExtensions.variational_circuit(5) +h = YaoExtensions.heisenberg(5) + +function loss(h, c, θ) where N + # the assign is nessesary! + c = dispatch(c, fill(θ, nparameters(c))) + reg = apply(zero_state(nqubits(c)), c) + real(expect(h, reg)) +end + +reg0 = zero_state(5) +zygote_grad = Zygote.gradient(θ->loss(h, c, θ), 0.5)[1] + + +# check gradients +using Test +dispatch!(c, fill(0.5, nparameters(c))) +greg, gparams = expect'(h, zero_state(5)=>c) +true_grad = sum(gparams) + +@test zygote_grad ≈ true_grad + +# the batched version +function loss2(h, c, θ) where N + # the assign is nessesary! + c = dispatch(c, fill(θ, nparameters(c))) + reg = zero_state(nqubits(c),nbatch=2) + reg = apply(reg, c) + sum(real(expect(h, reg))) +end + +reg0 = zero_state(5) +zygote_grad2 = Zygote.gradient(θ->loss2(h, c, θ), 0.5)[1] +true_grad2 = (loss2(h, c, 0.5+1e-5) - loss2(h, c, 0.5-1e-5)) / 2e-5 +@test true_grad2 ≈ zygote_grad2 \ No newline at end of file diff --git a/lib/QuAlgorithmZoo/examples/PortZygote/simple_example.jl b/lib/QuAlgorithmZoo/examples/PortZygote/simple_example.jl new file mode 100644 index 000000000..f3993c95f --- /dev/null +++ b/lib/QuAlgorithmZoo/examples/PortZygote/simple_example.jl @@ -0,0 +1,18 @@ +using Zygote +include("chainrules_patch.jl") + +import YaoExtensions, Random + +c = YaoExtensions.variational_circuit(5) +dispatch!(c, :random) + +function loss(reg::AbstractRegister, circuit::AbstractBlock{N}) where N + #copy(reg) |> circuit + reg = apply(copy(reg), circuit) + st = state(reg) + sum(real(st.*st)) +end + +reg0 = zero_state(5) +paramsδ = gradient(c->loss(reg0, c), c)[1] +regδ = gradient(reg->loss(reg, c), reg0)[1] diff --git a/lib/QuAlgorithmZoo/examples/QAOA/QAOA.jl b/lib/QuAlgorithmZoo/examples/QAOA/QAOA.jl new file mode 100644 index 000000000..c0d39427d --- /dev/null +++ b/lib/QuAlgorithmZoo/examples/QAOA/QAOA.jl @@ -0,0 +1,77 @@ +# required packages: NLopt, SCS and Convex +using Yao, BitBasis +using NLopt + +include("maxcut_gw.jl") + +HB(nbit::Int) = sum([put(nbit, i=>X) for i=1:nbit]) +tb = TimeEvolution(HB(3) |> cache, 0.1; check_hermicity=false) +function HC(W::AbstractMatrix) + nbit = size(W, 1) + 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 + sum(ab) +end + +function qaoa_circuit(W::AbstractMatrix, depth::Int; use_cache::Bool=false) + nbit = size(W, 1) + hb = HB(nbit) + hc = HC(W) + use_cache && (hb = hb |> cache; hc = hc |> cache) + c = chain(nbit, [repeat(nbit, H, 1:nbit)]) + append!(c, [chain(nbit, [time_evolve(hc, 0.0, tol=1e-5, check_hermicity=false), time_evolve(hb, 0.0, tol=1e-5, check_hermicity=false)]) for i=1:depth]) +end + + +function cobyla_optimize(circuit::AbstractBlock{N}, hc::AbstractBlock; niter::Int) where N + function f(params, grad) + reg = zero_state(N) |> dispatch!(circuit, params) + loss = expect(hc, reg) |> real + #println(loss) + loss + end + opt = Opt(:LN_COBYLA, nparameters(circuit)) + min_objective!(opt, f) + maxeval!(opt, niter) + cost, params, info = optimize(opt, parameters(circuit)) + pl = zero_state(N) |> circuit |> probs + cost, params, pl +end + +using Random, Test +@testset "qaoa circuit" begin + Random.seed!(2) + nbit = 5 + W = [0 5 2 1 0; + 5 0 3 2 0; + 2 3 0 0 0; + 1 2 0 0 4; + 0 0 0 4 0] + + # the exact solution + exact_cost, sets = goemansWilliamson(W) + + hc = HC(W) + @test ishermitian(hc) + # the actual loss is - + sum(wij)/2 + @test -expect(hc, product_state(5, bmask(sets[1]...)))+sum(W)/4 ≈ exact_cost + + # build a QAOA circuit and start training. + qc = qaoa_circuit(W, 20; use_cache=true) + opt_pl = nothing + opt_cost = Inf + for i = 1:10 + dispatch!(qc, :random) + cost, params, pl = cobyla_optimize(qc, hc; niter=2000) + @show cost + cost < opt_cost && (opt_pl = pl) + end + + # check the correctness of result + config = argmax(opt_pl)-1 + @test config in [bmask(set...) for set in sets] +end diff --git a/lib/QuAlgorithmZoo/examples/QAOA/imag_time_evolve.jl b/lib/QuAlgorithmZoo/examples/QAOA/imag_time_evolve.jl new file mode 100644 index 000000000..299a2b309 --- /dev/null +++ b/lib/QuAlgorithmZoo/examples/QAOA/imag_time_evolve.jl @@ -0,0 +1,51 @@ +"""imaginary time evolution solution to maximum cut problem.""" +using Yao + +@const_gate ZZ = mat(kron(Z,Z)) + +function maxcut_circuit(W::AbstractMatrix, τ) + nbit = size(W, 1) + ab = chain(nbit) + for i=1:nbit,j=i+1:nbit + if W[i,j] != 0 + #push!(ab, put(nbit, (i,j)=>0.5*W[i,j]*ZZ)) + push!(ab, put(nbit, (i,j)=>rot(ZZ, -im*W[i,j]*τ))) + end + end + return ab +end + +function maxcut_h(W::AbstractMatrix) + nbit = size(W, 1) + ab = Add{nbit}() + for i=1:nbit,j=i+1:nbit + if W[i,j] != 0 + #push!(ab, put(nbit, (i,j)=>0.5*W[i,j]*ZZ)) + push!(ab, put(nbit, (i,j)=>0.5*W[i,j]*ZZ)) + end + end + return ab +end + +using Random, Test +@testset "iqaoa circuit" begin + Random.seed!(2) + nbit = 5 + # exact solution is [(1,3,4), (2,5)] + W = [0 5 2 1 0; + 5 0 3 2 0; + 2 3 0 0 0; + 1 2 0 0 4; + 0 0 0 4 0] + + c = maxcut_circuit(W, 10) + h = maxcut_h(W) + reg = uniform_state(nbit) |> c + + # check the correctness of result + opt_pl = reg |> probs + config = argmax(opt_pl)-1 + @show BitStr64{nbit}(config), opt_pl[config+1] + @test config == bit"01101" || config == bit"10010" + @test isapprox(expect(h, reg |> normalize!), -5.5) +end diff --git a/lib/QuAlgorithmZoo/examples/QAOA/maxcut_gw.jl b/lib/QuAlgorithmZoo/examples/QAOA/maxcut_gw.jl new file mode 100644 index 000000000..dff4f928b --- /dev/null +++ b/lib/QuAlgorithmZoo/examples/QAOA/maxcut_gw.jl @@ -0,0 +1,83 @@ +using Convex +using SCS +using LinearAlgebra + +" +The Goemans-Williamson algorithm for the MAXCUT problem. + +From: https://github.com/ericproffitt/MaxCut +" + +function goemansWilliamson(W::Matrix{T}; tol::Real=1e-1, iter::Int=100) where T<:Real + "Partition a graph into two disjoint sets such that the sum of the edge weights + which cross the partition is as large as possible (known to be NP-hard)." + + "A cut of a graph can be produced by assigning either 1 or -1 to each vertex. The Goemans-Williamson + algorithm relaxes this binary condition to allow for vector assignments drawn from the (n-1)-sphere + (choosing an n-1 dimensional space will ensure seperability). This relaxation can then be written as + an SDP. Once the optimal vector assignments are found, origin centered hyperplanes are generated and + their corresponding cuts evaluated. After 'iter' trials, or when the desired tolerance is reached, + the hyperplane with the highest corresponding binary cut is used to partition the vertices." + + "W: Adjacency matrix." + "tol: Maximum acceptable distance between a cut and the MAXCUT upper bound." + "iter: Maximum number of hyperplane iterations before a cut is chosen." + + LinearAlgebra.checksquare(W) + @assert LinearAlgebra.issymmetric(W) "Adjacency matrix must be symmetric." + @assert all(W .>= 0) "Entries of the adjacency matrix must be nonnegative." + @assert all(diag(W) .== 0) "Diagonal entries of adjacency matrix must be zero." + @assert tol > 0 "The tolerance 'tol' must be positive." + @assert iter > 0 "The number of iterations 'iter' must be a positive integer." + + "This is the standard SDP Relaxation of the MAXCUT problem, a reference can be found at + http://www.sfu.ca/~mdevos/notes/semidef/GW.pdf." + k = size(W, 1) + S = Semidefinite(k) + + expr = dot(W, S) + constr = [S[i,i] == 1.0 for i in 1:k] + problem = minimize(expr, constr...) + solve!(problem, SCS.Optimizer(verbose=0)) + + ### Ensure symmetric positive-definite. + A = 0.5 * (S.value + S.value') + A += (max(0, -eigmin(A)) + eps(1e3))*I + + X = Matrix(cholesky(A)) + + ### A non-trivial upper bound on MAXCUT. + upperbound = (sum(W) - dot(W, S.value)) / 4 + + "Random origin-centered hyperplanes, generated to produce partitions of the graph." + maxcut = 0 + maxpartition = nothing + + for i in 1:iter + gweval = X' * randn(k) + partition = (findall(x->x>0, gweval), findall(x->x<0, gweval)) + cut = sum(W[partition...]) + + if cut > maxcut + maxpartition = partition + maxcut = cut + end + + upperbound - maxcut < tol && break + i == iter && println("Max iterations reached.") + end + return round(maxcut; digits=3), maxpartition +end + +using Test +@testset "max cut" begin + W = [0 5 2 1 0; + 5 0 3 2 0; + 2 3 0 0 0; + 1 2 0 0 4; + 0 0 0 4 0] + + maxcut, maxpartition = goemansWilliamson(W) + @test maxcut == 14 + @test [2,5] in maxpartition +end diff --git a/lib/QuAlgorithmZoo/examples/QCBM/QCBM.jl b/lib/QuAlgorithmZoo/examples/QCBM/QCBM.jl new file mode 100644 index 000000000..1dbce8ac8 --- /dev/null +++ b/lib/QuAlgorithmZoo/examples/QCBM/QCBM.jl @@ -0,0 +1,54 @@ +# # Quantum Circuit Born Machine +using Yao, YaoExtensions +import Yao: probs +using QuAlgorithmZoo: Adam, update! + +struct QCBM{BT<:AbstractBlock, MT<:MMD} + circuit::BT + mmd::MT +end + +"""generated output probability distribution""" +function probs(qcbm::QCBM) + zero_state(qcbm.circuit |> nqubits) |> qcbm.circuit |> probs +end + +function loss(qcbm::QCBM) + expect(qcbm.mmd, probs(qcbm) |> as_weights) +end + +function getgrad(qcbm::QCBM) + expect'(qcbm.mmd, zero_state(nqubits(qcbm.circuit))=>qcbm.circuit).second +end + +# ## DATA: Target Probability Distribution +# The gaussian probability disctribution in phase space of 2^6 +nbit = 6 +N = 1<put(nbit, i=>Z), +, (1:Na..., Na+Nc+1:Na+Nb...)) + params = parameters(c) + for i = 1:maxiter + grad = expect'(obs, reg => c).second + QuAlgorithmZoo.update!(params, grad, optimizer) + println("Iter $i, Loss = $(Na+expect(obs, copy(reg) |> c))") + dispatch!(c, params) + end +end + +"""build the circuit for quantum SVD training.""" +function circuit_qsvd(circuit_a::AbstractBlock{Na}, circuit_b::AbstractBlock{Nb}, Nc::Int) where {Na, Nb} + nbit = Na+Nb + cnots = chain(control(nbit, i+Na, i=>X) for i=1:Nc) + c = chain(subroutine(nbit, circuit_a, 1:Na), subroutine(nbit, circuit_b, Na+1:nbit), cnots) +end + +"""read QSVD results""" +function readout_qsvd(reg::AbstractRegister, circuit_a::AbstractBlock{Na}, circuit_b::AbstractBlock{Nb}, Nc::Int) where {Na, Nb} + reg = copy(reg) |> subroutine(Na+Nb, circuit_a, 1:Na) |> subroutine(Na+Nb, circuit_b, Na+1:Na+Nb) + _S = [select(reg, b|b< abs2, I, atol=1e-2) + @test isapprox(V'*V_exact .|> abs2, I, atol=1e-2) +end diff --git a/lib/QuAlgorithmZoo/examples/QSVD/README.md b/lib/QuAlgorithmZoo/examples/QSVD/README.md new file mode 100644 index 000000000..8d259f71f --- /dev/null +++ b/lib/QuAlgorithmZoo/examples/QSVD/README.md @@ -0,0 +1,2 @@ +# Quantum SVD +## Reference: https://arxiv.org/abs/1905.01353 diff --git a/lib/QuAlgorithmZoo/examples/QuGAN/QuGAN.jl b/lib/QuAlgorithmZoo/examples/QuGAN/QuGAN.jl new file mode 100644 index 000000000..7bb4b968c --- /dev/null +++ b/lib/QuAlgorithmZoo/examples/QuGAN/QuGAN.jl @@ -0,0 +1,78 @@ +using Yao +using YaoExtensions: variational_circuit, Sequence, faithful_grad, numdiff +using QuAlgorithmZoo: Adam, update! +import Yao: tracedist + +""" +Quantum GAN. + +Reference: + Benedetti, M., Grant, E., Wossnig, L., & Severini, S. (2018). Adversarial quantum circuit learning for pure state approximation, 1–14. +""" +struct QuGAN{N} + target::ArrayReg + generator::AbstractBlock{N} + discriminator::AbstractBlock + reg0::ArrayReg + witness_op::AbstractBlock + circuit::AbstractBlock + + function QuGAN(target::ArrayReg, gen::AbstractBlock, dis::AbstractBlock) + N = nqubits(target) + c = Sequence([gen, addbits!(1), dis]) + witness_op = put(N+1, (N+1)=>ConstGate.P0) + new{N}(target, gen, dis, zero_state(N), witness_op, c) + end +end + +# INTERFACES +circuit(qg::QuGAN) = qg.circuit +loss(qg::QuGAN) = p0t(qg) - p0g(qg) + +function gradient(qg::QuGAN) + grad_gen = faithful_grad(qg.witness_op, qg.reg0 => qg.circuit) + grad_tar = faithful_grad(qg.witness_op, qg.target => qg.circuit[2:end]) + ngen = nparameters(qg.generator) + [-grad_gen[1:ngen]; grad_tar - grad_gen[ngen+1:end]] +end + +"""probability to get evidense qubit 0 on generation set.""" +p0g(qg::QuGAN) = expect(qg.witness_op, qg.reg0 => qg.circuit) |> real +"""probability to get evidense qubit 0 on target set.""" +p0t(qg::QuGAN) = expect(qg.witness_op, qg.target => qg.circuit[2:end]) |> real +"""generated wave function""" +outputψ(qg::QuGAN) = copy(qg.reg0) |> qg.generator + +"""tracedistance between target and generated wave function""" +tracedist(qg::QuGAN) = tracedist(qg.target, outputψ(qg))[] + +using Test, Random +Random.seed!(2) + +nbit = 3 +depth_gen = 4 +depth_dis = 4 + +# define a QuGAN +target = rand_state(nbit) +generator = dispatch!(variational_circuit(nbit, depth_gen), :random) +discriminator = dispatch!(variational_circuit(nbit+1, depth_dis), :random) +qg = QuGAN(target, generator, discriminator) + +# check the gradient +grad = gradient(qg) +numgrad = numdiff(c->loss(qg), qg.circuit) +@test isapprox(grad, numgrad, atol=1e-4) + +# learning rates for the generator and discriminator +g_lr = 0.2 +d_lr = 0.5 +for i=1:300 + ng = nparameters(qg.generator) + grad = gradient(qg) + dispatch!(-, qg.generator, grad[1:ng]*g_lr) + dispatch!(-, qg.discriminator, -grad[ng+1:end]*d_lr) + println("Step $i, trace distance = $(tracedist(qg))") +end + +@test qg |> loss < 0.1 diff --git a/lib/QuAlgorithmZoo/examples/Shor/Shor.jl b/lib/QuAlgorithmZoo/examples/Shor/Shor.jl new file mode 100644 index 000000000..87bca87d2 --- /dev/null +++ b/lib/QuAlgorithmZoo/examples/Shor/Shor.jl @@ -0,0 +1,125 @@ +# # [Shor's Algorithm](@id Shor) + +# ## References +# * [Neilsen](https://aapt.scitation.org/doi/abs/10.1119/1.1463744?journalCode=ajp) +# * [An Insightful Blog](https://algassert.com/post/1718) + +# ## Main Program +# The main program of a Shor's algorithm can be summrized in several lines of code. +# For the theory part, please refer the reference materials above. +# It factorize an integer `L`, and returns one of the factors. +# Here, the input `ver` can be either `Val(:quantum)` or `Val(:classical)`, +# where the classical version is for comparison. + +using Yao, BitBasis +using YaoExtensions: KMod, QFTCircuit +using QuAlgorithmZoo: NumberTheory + +function shor(L::Int, ver=Val(:quantum); maxtry=100) + L%2 == 0 && return 2 + + ## find short cut solutions like `a^b` + res = NumberTheory.factor_a_power_b(L) + res !== nothing && return res[1] + + for i in 1:maxtry + ## step 1 + x = NumberTheory.rand_primeto(L) + + ## step 2 + r = get_order(ver, x, L; ) + if r%2 == 0 && powermod(x, r÷2, L) != L-1 + ## step 3 + f1, f2 = gcd(powermod(x, r÷2, L)-1, L), gcd(powermod(x, r÷2, L)+1, L) + if f1!=1 + return f1 + elseif f2!=1 + return f2 + else + error("Algorithm Fail!") + end + end + end +end + +# Except some shortcuts, in each try, the main program can be summarized in several steps +# 1. randomly pick a number that prime to the input numebr `L`, i.e. `gcd(x, L) = 1`. +# The complexity of this algorithm is polynoial. +# 2. get the order `x`, i.e. finding a number `r` that satisfies `mod(x^r, L) = 1`. +# If `r` is even and `x^(r÷2)` is non-trivil, go on, otherwise start another try. +# Here, trivil means equal to `L-1 (mod L)`. +# 3. According to Theorem 5.2 in Neilsen book, +# one of `gcd(x^(r÷2)-1, L)` and `gcd(x^(r÷2)+1, L)` must be a non-trivil (`!=1`) factor of `L`. +# Notice `powermod(x, r÷2, L)` must be `-1` rather than `1`, +# otherwise the order should be `r/2` according to definition. + +# The only difference between classical and quantum version is the order finding algorithm. + +# ## Order Finding +# We provided a classical order finding algorithm in `NumberTheory`, +# here we focus on the quantum version. +# The algorithm is consisted +# 1. run the circuit to get a bitstring, +# 2. interpret this bitstring in output register as a rational number `s/r`. +# To achieve this, we first interpret it as a floating point number, +# then the continued fraction algorithm can find the best match for us. +# +# When using the quantum version, we have the flexibility to set key word arguments `nshot`, +# `nbit` (size of input data register) and `ncbit` (size of control register, or output register). +# `nbit` can be simply chosen as the minimum register size to store input, +# while `ncbit` can be estimated with the following function +"""estimate the required size of the output register.""" +estimate_ncbit(nbit::Int, ϵ::Real) = 2*nbit + 1 + ceil(Int,log2(2+1/2ϵ)) + +get_order(::Val{:classical}, x::Int, L::Int; kwargs...) = NumberTheory.find_order(x, L) +function get_order(::Val{:quantum}, x::Int, L::Int; nshots::Int=10, + nbit::Int=bit_length(L-1), ncbit::Int=estimate_ncbit(nbit, 0.25)) + c = order_finding_circuit(x, L; nbit=nbit, ncbit=ncbit) + reg = join(product_state(nbit, 1), zero_state(ncbit)) + + res = measure(copy(reg) |> c; nshots=nshots) + for r in res + ## split bit string b into lower bits `k` and higher bits `r`. + mask = bmask(1:ncbit) + k,i = r&mask, r>>ncbit + ## get s/r + ϕ = bfloat(k) # + ϕ == 0 && continue + + ## order_from_float: given a floating point number, + ## return the closest rational number with bounded number of continued fraction steps. + order = NumberTheory.order_from_float(ϕ, x, L) + if order === nothing + continue + else + return order + end + end + return nothing +end + +# #### The circuit used for finding order +""" + order_finding_circuit(x::Int, L::Int; nbit::Int=bit_length(L-1), ncbit::Int=estimate_ncbit(nbit, 0.25)) -> AbstractBlock + +Returns the circuit for finding the order of `x` to `L`, +feeding input `|1>⊗|0>` will get the resulting quantum register with the desired "phase" information. +""" +function order_finding_circuit(x::Int, L::Int; nbit::Int, ncbit::Int) + N = nbit+ncbit + chain(N, repeat(N, H, 1:ncbit), KMod{N, ncbit}(x, L), + subroutine(N, QFTCircuit(ncbit)', 1:ncbit)) +end + +# The circuit for order finding is consisted of three parts +# 1. Hadamard gates, +# 2. `KMod` that computes a classical function `mod(a^k*x, L)`. +# `k` is the integer stored in first `K` (or `ncbit`) qubits and the rest `N-K` qubits stores `a`. +# Notice it is not a basic gate, it should have been compiled to multiple gates, which is not implemented in `Yao` for the moment. +# To learn more about implementing arithmatics on a quantum circuit, please read [this paper](https://arxiv.org/abs/1805.12445). +# 3. Inverse quantum fourier transformation. + +# ## Run +# Factorizing `15`, you should see `3` or `5`, please report a bug if it is not... +using Random; Random.seed!(129) #src +shor(15, Val(:quantum)) diff --git a/lib/QuAlgorithmZoo/examples/Shor/tests.jl b/lib/QuAlgorithmZoo/examples/Shor/tests.jl new file mode 100644 index 000000000..8c0c104ab --- /dev/null +++ b/lib/QuAlgorithmZoo/examples/Shor/tests.jl @@ -0,0 +1,44 @@ +include("Shor.jl") +using Test + +"""Euler theorem states that the order is a devisor of Eulerφ (or the size of Z* group of `N`)""" +function check_Euler_theorem(N::Int) + Z = NumberTheory.Z_star(N) + Nz = length(Z) # Eulerφ + for x in Z + @test powermod(x,Nz,N) == 1 # the order is a devisor of Eulerφ + end +end + +@testset "Euler" begin + check_Euler_theorem(150) +end + +@testset "shor_classical" begin + Random.seed!(129) + L = 35 + f = shor(L, Val(:classical)) + @test f == 5 || f == 7 + + L = 25 + f = shor(L, Val(:classical)) + @test f == 5 + + L = 7*11 + f = shor(L, Val(:classical)) + @test f == 7 || f == 11 + + L = 14 + f = shor(L, Val(:classical)) + @test f == 2 || f == 7 + + @test NumberTheory.factor_a_power_b(25) == (5, 2) + @test NumberTheory.factor_a_power_b(15) == nothing +end + +@testset "shor quantum" begin + Random.seed!(129) + L = 15 + f = shor(L, Val(:quantum)) + @test f == 5 || f == 3 +end diff --git a/lib/QuAlgorithmZoo/examples/StructureOpt/StructureOpt.jl b/lib/QuAlgorithmZoo/examples/StructureOpt/StructureOpt.jl new file mode 100644 index 000000000..7fb5da1fb --- /dev/null +++ b/lib/QuAlgorithmZoo/examples/StructureOpt/StructureOpt.jl @@ -0,0 +1,46 @@ +using Yao, YaoExtensions +using YaoBlocks.Optimise: replace_block +using YaoBlocks.ConstGate +using YaoBlocks +using Test +using LinearAlgebra + +""" + TODO: + * optimize to optimal value, + * if this optimal value's loss is close to `θ = 0.0`, kill it. +""" +performance(params, gradient) = abs.(params .* gradient) + +nbit = 5 +c = variational_circuit(nbit, 50) + +dispatch!(c, :random) +h = heisenberg(nbit) +cc = replace_block(x->x isa RotationGate ? Bag(x) : x, c) + +function train(h, c; kill_rate=0.1, kill_step=10, niter=100) + bags = blockfilter(b->b isa Bag && isenabled(b), c) + for i=1:niter + regδ, paramsδ = expect'(h, zero_state(nbit) => c) + dispatch!(-, c, 0.01 .* paramsδ) + loss = expect(h, zero_state(nbit)=>c) |> real + + Nk = round(Int, kill_rate*length(paramsδ)) + + if i%kill_step == 0 && kill_rate != 0 + ps = performance(parameters(c), paramsδ) + cut = sort(ps)[Nk] + disable_block!.(bags[ps .<= cut]) + bags = bags[ps .> cut] + end + @show i,nparameters(c), loss + end +end + + +train(h, cc; kill_rate=0.01, niter=200) + +cr = variational_circuit(nbit, 5) +dispatch!(cr, :random) +train(h, cr; kill_rate=0.0, niter=200) diff --git a/lib/QuAlgorithmZoo/examples/TraceOperator/README.md b/lib/QuAlgorithmZoo/examples/TraceOperator/README.md new file mode 100644 index 000000000..523b5f3f2 --- /dev/null +++ b/lib/QuAlgorithmZoo/examples/TraceOperator/README.md @@ -0,0 +1,8 @@ +# Tracing an operator + +To compute $\Re[{\rm Tr}(Ue^{-i\theta})]$, we utilize the Hadamard test. + +The circuit is +![trace circuit](traceop.png) + +This result can be proved from its tensor network representation. \ No newline at end of file diff --git a/lib/QuAlgorithmZoo/examples/TraceOperator/TraceOperator.jl b/lib/QuAlgorithmZoo/examples/TraceOperator/TraceOperator.jl new file mode 100644 index 000000000..992ff2f85 --- /dev/null +++ b/lib/QuAlgorithmZoo/examples/TraceOperator/TraceOperator.jl @@ -0,0 +1,17 @@ +using Yao +using QuAlgorithmZoo +using YaoBlocks.ConstGate +using Test +using LinearAlgebra + +@testset "hadamard test" begin + n = 2 + U = chain(put(n, 2=>Rx(0.2)), put(n, 1=>Rz(1.2)), put(n, 1=>phase(0.4))) + US = chain(2n, put(2n, (3,4)=>U), + chain(2n, [swap(2n,i,i+n) for i=1:n])) + reg = zero_state(2n) + reg |> repeat(2n, H, 1:n) |> chain(2n, [cnot(2n,i,i+n) for i=1:n]) + @show real(tr(mat(U)))/(1<Z) + Z2 = put(2,2=>Z) + X1 = put(2,1=>X) + X2 = put(2,2=>X) + 0.011280*Z1*Z2 + 0.397936*Z1 + 0.397936*Z2 + 0.180931*X1*X2 +end + +using Flux: Optimise +function train!(circ, hamiltonian; optimizer, niter::Int=100) + params = parameters(circ) + dispatch!(circ, :random) + for i=1:niter + _, grad = expect'(hamiltonian, zero_state(nqubits(circ)) => circ) + Optimise.update!(optimizer, params, grad) + dispatch!(circ, params) + println("Energy = $(expect(hamiltonian, zero_state(nqubits(hamiltonian)) |> circ) |> real)") + end + return expect(hamiltonian, zero_state(nqubits(hamiltonian)) |> circ) +end + +h = hydrogen_hamiltonian() +c = variational_circuit(2) +emin_vqe = train!(c, h; optimizer=Optimise.ADAM(0.1)) + +using LinearAlgebra +emin = eigvals(Matrix(mat(h)))[1] +@assert isapprox(emin, emin, atol=1e-1) diff --git a/lib/QuAlgorithmZoo/examples/VQE/H2_OpenFermion.jl b/lib/QuAlgorithmZoo/examples/VQE/H2_OpenFermion.jl new file mode 100644 index 000000000..b924fd9b1 --- /dev/null +++ b/lib/QuAlgorithmZoo/examples/VQE/H2_OpenFermion.jl @@ -0,0 +1,110 @@ +# A more precise example of $\text H_2$ with 4 qubits hamiltonian. + +using Yao +using YaoExtensions + +# Here we we invoke [OpenFermion](https://github.com/quantumlib/OpenFermion) in Python +# to get the hamiltonian under [Jordan-Weigner transformation](https://en.wikipedia.org/wiki/Jordan%E2%80%93Wigner_transformation) +# (with 2 electrons and 4 sto-3g basis set orbitals). + +using PyCall + +py""" +import numpy as np +from openfermion.config import * +from openfermion.hamiltonians import MolecularData +from openfermion.transforms import jordan_wigner + +def make_H2_mol(bond_len): + + ## Returns: + ## molecule(openfermion.hamiltonians.MolecularData): + ## MolecularData for H2 with certain bond length + + ## load molecule + atom_1 = 'H' + atom_2 = 'H' + basis = 'sto-3g' + multiplicity = 1 + charge = 0 + coordinate_1 = (0.0, 0.0, 0.0) + coordinate_2 = (0.0, 0.0, bond_len) + geometry = [(atom_1, coordinate_1), (atom_2, coordinate_2)] + molecule = MolecularData(geometry, basis, multiplicity, + charge, description=str(bond_len)) + ## molecule = run_pyscf(molecule,run_scf=1,run_ccsd=1,run_fci=1) + ## can calculate on arbitrary bond length if package openfermionpyscf is available, + ## only existing molecule data can be used + molecule.load() + return molecule + +def get_hamiltonian(bond_len): + ## get the coefficients and pauli terms of Hamiltonian of H_2 with given bond length + ## Returns: + ## (PyArray,PyArray) + molecule = make_H2_mol(bond_len) + qubit_hamiltonian = jordan_wigner(molecule.get_molecular_hamiltonian()) + qubit_hamiltonian.compress() + return list(qubit_hamiltonian.terms.keys()),list(qubit_hamiltonian.terms.values()) + +""" + +# Define the function to transfer OpenFermion hamiltonian to Yao circuit + +function get_fermion_hamiltonian(n::Int,terms::Array,coefs::Vector) + gates=Dict("Z"=>Z,"X"=>X,"Y"=>Y) + to_pauli(t::Tuple{Int,String})=put(n,t[1]+1=>get(gates,t[2],ErrorException("Invalid"))) + return sum(coefs[2:end].*map(x->reduce(*,map(to_pauli,x)),terms[2:end])) +end + +# Use VQE code given in the previous Yao example + +using Flux: Optimise + +function train!(circ, hamiltonian; optimizer, niter::Int=200,verbose::Bool=true) + params = parameters(circ) + dispatch!(circ, :random) + for i=1:niter + _, grad = expect'(hamiltonian, zero_state(nqubits(circ)) => circ) + Optimise.update!(optimizer, params, grad) + dispatch!(circ, params) + if verbose + println("Energy = + $(expect(hamiltonian, zero_state(nqubits(hamiltonian)) |> circ) |> real)") + end + end + return expect(hamiltonian, zero_state(nqubits(hamiltonian)) |> circ) +end + +# Train VQE to search ground energy on various bond lengths. + +bond_lens = Array(0.2:0.1:1.5)#\AA +pes = Vector{Real}() + +for l in bond_lens + terms, coefs = py"get_hamiltonian"(l) + e = coefs[1] + + + h = get_fermion_hamiltonian(4,terms,coefs) + c = variational_circuit(4) + emin_vqe = train!(c, h; optimizer=Optimise.ADAM(0.1),verbose=false) + println("Estimated minimum energy at bond length $(l): $(e+emin_vqe|>real)") + push!(pes,e+emin_vqe|>real) + + using LinearAlgebra + emin = eigvals(Matrix(mat(h)))[1] + @assert isapprox(emin, emin_vqe, atol=1e-2) +end + +# Plot the potential energy surface (PES) + +using Plots + +plot(bond_lens,pes, markershape = :circle, linestyle = :dash, label = "VQE",ylabel = "Ground State Energy (Hartree)", xlabel = "Bond Length (Å) ") + +# Since OpenFermion only generates molecule data at discrete bond length with minimum step to be 0.1 Å, +# +# the VQE gives lowest energy -1.1361862609 Hartree at 0.7 Å. +# +# As a comparison, the ground energy at stable bond length given by FCI (Full Configuration Interaction) is -1.13730688 Hartree at 0.734832 Å. diff --git a/lib/QuAlgorithmZoo/examples/VQE/README.md b/lib/QuAlgorithmZoo/examples/VQE/README.md new file mode 100644 index 000000000..a0fb31738 --- /dev/null +++ b/lib/QuAlgorithmZoo/examples/VQE/README.md @@ -0,0 +1,32 @@ +# Variational Quantum Eigensolver + +## Start + +1. install [Julia 1.1](https://julialang.org/downloads/) +2. type `]` in julia REPL to enter `pkg` mode, and install packages with +```julia +pkg> add Yao KrylovKit +pkg> dev git@github.com:QuantumBFS/YaoExtensions.jl.git +``` +3. run H2 example with +```bash +julia examples/VQE/H2.jl +``` +## More +4. set your Python environment in Julia +``` +julia> ENV["PYTHON"] = "... path of the python executable ..." +``` +5. install [PyCall](https://github.com/JuliaPy/PyCall.jl) +``` +pkg> add PyCall +pkg> build PyCall +``` +6. install [OpenFermion](https://github.com/quantumlib/OpenFermion) +```bash +pip install openfermion +``` +7. run H2_OpenFermion example +```bash +julia examples/VQE/H2_OpenFermion.jl +``` diff --git a/lib/QuAlgorithmZoo/examples/VQE/VQE.jl b/lib/QuAlgorithmZoo/examples/VQE/VQE.jl new file mode 100644 index 000000000..0d27a585a --- /dev/null +++ b/lib/QuAlgorithmZoo/examples/VQE/VQE.jl @@ -0,0 +1,64 @@ +# # [Variational Quantum Eigensolver](@id VQE) + +# ## References +# * [A variational eigenvalue solver on a quantum processor](https://arxiv.org/abs/1304.3061) +# * [Variational Quantum Eigensolver with Fewer Qubits](https://arxiv.org/abs/1902.02663) + +# ## Define a hamiltonian + +# construct a 5-site heisenberg hamiltonian + +using Yao, YaoExtensions + +N = 5 +hami = heisenberg(N) + +# The ground state can be obtained by a sparse matrix ground state solver. +# The high performance `mat` function in `Yao.jl` makes computation time lower than `10s` +# to construct a `20` site Heisenberg hamltonian +using KrylovKit: eigsolve + +function ed_groundstate(h::AbstractBlock) + E, V = eigsolve(h |> mat, 1, :SR, ishermitian=true) + E[1], V[1] +end + +ed_groundstate(hami) + +# Here we use the `heisenberg` hamiltonian that defined in [`YaoExtensions.jl`](https://github.com/QuantumBFS/YaoExtensions.jl), +# for tutorial purpose, we pasted the code for construction here. +# ```julia +# function heisenberg(nbit::Int; periodic::Bool=true) +# sx = i->put(nbit, i=>X) +# sy = i->put(nbit, i=>Y) +# sz = i->put(nbit, i=>Z) +# map(1:(periodic ? nbit : nbit-1)) do i +# j=i%nbit+1 +# sx(i)*sx(j)+sy(i)*sy(j)+sz(i)*sz(j) +# end |> sum +# end +# ``` + +# ## Define an ansatz +# As an ansatz, we use the canonical circuit for demonstration [`variational_circuit`](@ref YaoExtensions.variational_circuit) +# defined in [`YaoExtensions.jl`](https://github.com/QuantumBFS/YaoExtensions.jl). +c = variational_circuit(N) +dispatch!(c, :random) + +# ## Run +# Use the [`Adam`](@ref) optimizer for parameter optimization, +# we provide a poorman's implementation in `QuAlgorithmZoo` +using QuAlgorithmZoo: Adam, update! + +optimizer = Adam(lr=0.01) +params = parameters(c) +niter = 100 + +for i = 1:niter + ## `expect'` gives the gradient of an observable. + grad_input, grad_params = expect'(hami, zero_state(N) => c) + + ## feed the gradients into the circuit. + dispatch!(c, update!(params, grad_params, optimizer)) + println("Step $i, Energy = $(expect(hami, zero_state(N) |> c))") +end diff --git a/lib/QuAlgorithmZoo/examples/run_examples.jl b/lib/QuAlgorithmZoo/examples/run_examples.jl new file mode 100644 index 000000000..222510dc2 --- /dev/null +++ b/lib/QuAlgorithmZoo/examples/run_examples.jl @@ -0,0 +1,9 @@ +# run this file to make sure all examples work +# Note: running all these examples are too heavy for CI. +include("../examples/Grover/Grover.jl") +include("../examples/QAOA/QAOA.jl") +include("../examples/HHL/HHL.jl") +include("../examples/QCBM/QCBM.jl") +include("../examples/QuGAN/QuGAN.jl") +include("../examples/Shor/Shor.jl") +include("../examples/VQE/VQE.jl") diff --git a/lib/QuAlgorithmZoo/src/Adam.jl b/lib/QuAlgorithmZoo/src/Adam.jl new file mode 100644 index 000000000..3183e70a6 --- /dev/null +++ b/lib/QuAlgorithmZoo/src/Adam.jl @@ -0,0 +1,40 @@ +export Adam + +mutable struct Adam + lr::AbstractFloat + gclip::AbstractFloat + beta1::AbstractFloat + beta2::AbstractFloat + eps::AbstractFloat + t::Int + fstm + scndm +end + +Adam(; lr=0.001, gclip=0, beta1=0.9, beta2=0.999, eps=1e-8)=Adam(lr, gclip, beta1, beta2, eps, 0, nothing, nothing) + +function update!(w, g, p::Adam) + gclip!(g, p.gclip) + if p.fstm===nothing; p.fstm=zero(w); p.scndm=zero(w); end + p.t += 1 + lmul!(p.beta1, p.fstm) + BLAS.axpy!(1-p.beta1, g, p.fstm) + lmul!(p.beta2, p.scndm) + BLAS.axpy!(1-p.beta2, g .* g, p.scndm) + fstm_corrected = p.fstm / (1 - p.beta1 ^ p.t) + scndm_corrected = p.scndm / (1 - p.beta2 ^ p.t) + BLAS.axpy!(-p.lr, @.(fstm_corrected / (sqrt(scndm_corrected) + p.eps)), w) +end + +function gclip!(g, gclip) + if gclip == 0 + g + else + gnorm = vecnorm(g) + if gnorm <= gclip + g + else + BLAS.scale!(gclip/gnorm, g) + end + end +end diff --git a/lib/QuAlgorithmZoo/src/HadamardTest.jl b/lib/QuAlgorithmZoo/src/HadamardTest.jl new file mode 100644 index 000000000..1d823930e --- /dev/null +++ b/lib/QuAlgorithmZoo/src/HadamardTest.jl @@ -0,0 +1,33 @@ +export hadamard_test, hadamard_test_circuit, swap_test_circuit + +""" +see WiKi. +""" +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 + +""" +Estimation of overlap between multiple density matrices. + +PRL 88.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 diff --git a/lib/QuAlgorithmZoo/src/PhaseEstimation.jl b/lib/QuAlgorithmZoo/src/PhaseEstimation.jl new file mode 100644 index 000000000..5aef22073 --- /dev/null +++ b/lib/QuAlgorithmZoo/src/PhaseEstimation.jl @@ -0,0 +1,47 @@ +export PEBlock, projection_analysis +""" + PEBlock(UG, n_reg, n_b) -> ChainBlock + +phase estimation circuit. + + * `UG`: the input unitary matrix. + * `n_reg`: the number of bits to store phases, + * `n_b`: the number of bits to store vector. +""" +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) + + # Construct a control circuit. + control_circuit = chain(nbit) + for i = 1:n_reg + push!(control_circuit, control(nbit, (i,), (n_reg+1:nbit...,)=>UG)) + if i != n_reg + UG = matblock(mat(UG) * mat(UG)) + end + end + + # Inverse QFT Block. + iqft = subroutine(nbit, QFTBlock{n_reg}()',[1:n_reg...,]) + chain(hs, control_circuit, iqft) +end + +""" + 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::ArrayReg) + overlap = evec'*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 diff --git a/lib/QuAlgorithmZoo/src/QuAlgorithmZoo.jl b/lib/QuAlgorithmZoo/src/QuAlgorithmZoo.jl new file mode 100644 index 000000000..168ee3fd7 --- /dev/null +++ b/lib/QuAlgorithmZoo/src/QuAlgorithmZoo.jl @@ -0,0 +1,15 @@ +module QuAlgorithmZoo + +using LinearAlgebra +using Yao, BitBasis +using YaoExtensions + +include("Adam.jl") +include("PhaseEstimation.jl") +include("hamiltonian_solvers.jl") +include("HadamardTest.jl") +include("number_theory.jl") + +@deprecate random_diff_circuit variational_circuit + +end # module diff --git a/lib/QuAlgorithmZoo/src/hamiltonian_solvers.jl b/lib/QuAlgorithmZoo/src/hamiltonian_solvers.jl new file mode 100644 index 000000000..bb33b65ac --- /dev/null +++ b/lib/QuAlgorithmZoo/src/hamiltonian_solvers.jl @@ -0,0 +1,49 @@ +export iter_groundstate!, itime_groundstate!, vqe_solve! + +""" + iter_groundstate!({reg::AbstractRegister}, h::AbstractBlock; niter::Int=100) -> AbstractRegister + +project wave function to ground state by iteratively apply -h. +""" +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! + end + reg |> normalize! +end + +""" + 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::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 |> time_evolve(h, τ%span) |> normalize! + end + reg +end + +""" + 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::AbstractBlock{N}, hamiltonian::AbstractBlock; niter::Int=100) where N + optimizer = Adam(lr=0.01) + params = parameters(circuit) + for i = 1:niter + grad = expect'(hamiltonian, zero_state(N) => circuit).second + dispatch!(circuit, update!(params, grad, optimizer)) + println("Step $i, Energy = $(expect(hamiltonian, zero_state(N) |> circuit))") + end + circuit +end diff --git a/lib/QuAlgorithmZoo/src/number_theory.jl b/lib/QuAlgorithmZoo/src/number_theory.jl new file mode 100644 index 000000000..355e281ae --- /dev/null +++ b/lib/QuAlgorithmZoo/src/number_theory.jl @@ -0,0 +1,101 @@ +export NumberTheory + +module NumberTheory + +export Z_star, Eulerφ, continued_fraction, mod_inverse, rand_primeto, factor_a_power_b +export is_order, order_from_float, find_order + +""" + Z_star(N::Int) -> Vector + +returns the Z* group elements of `N`, i.e. {x | gcd(x, N) == 1} +""" +Z_star(N::Int) = filter(i->gcd(i, N)==1, 0:N-1) +Eulerφ(N) = length(Z_star(N)) + +""" + continued_fraction(ϕ, niter::Int) -> Rational + +obtain `s` and `r` from `ϕ` that satisfies `|s/r - ϕ| ≦ 1/2r²` +""" +continued_fraction(ϕ, niter::Int) = niter==0 ? floor(Int, ϕ) : floor(Int, ϕ) + 1//continued_fraction(1/mod(ϕ, 1), niter-1) +continued_fraction(ϕ::Rational, niter::Int) = niter==0 || ϕ.den==1 ? floor(Int, ϕ) : floor(Int, ϕ) + 1//continued_fraction(1/mod(ϕ, 1), niter-1) + +""" + mod_inverse(x::Int, N::Int) -> Int + +Return `y` that `(x*y)%N == 1`, notice the `(x*y)%N` operations in Z* forms a group and this is the definition of inverse. +""" +function mod_inverse(x::Int, N::Int) + for i=1:N + (x*i)%N == 1 && return i + end + throw(ArgumentError("Can not find the inverse, $x is probably not in Z*($N)!")) +end + +""" + is_order(r, x, N) -> Bool + +Returns true if `r` is the order of `x`, i.e. `r` satisfies `x^r % N == 1`. +""" +is_order(r, x, N) = powermod(x, r, N) == 1 + +""" + find_order(x::Int, N::Int) -> Int + +Find the order of `x` by brute force search. +""" +function find_order(x::Int, N::Int) + findfirst(r->is_order(r, x, N), 1:N) +end + +""" + rand_primeto(N::Int) -> Int + +Returns a random number `2 ≦ x < N` that is prime to `N`. +""" +function rand_primeto(N::Int) + while true + x = rand(2:N-1) + d = gcd(x, N) + if d == 1 + return x + end + end +end + +""" + order_from_float(ϕ, x, L) -> Int + +Estimate the order of `x` to `L`, `r`, from a floating point number `ϕ ∼ s/r` using the continued fraction method. +""" +function order_from_float(ϕ, x, L) + k = 1 + rnum = continued_fraction(ϕ, k) + while rnum.den < L + r = rnum.den + if is_order(r, x, L) + return r + end + k += 1 + rnum = continued_fraction(ϕ, k) + end + return nothing +end + +""" + factor_a_power_b(N::Int) -> (Int, Int) or nothing + +Factorize `N` into the power form `a^b`. +""" +function factor_a_power_b(N::Int) + y = log2(N) + for b = 2:ceil(Int, y) + x = 2^(y/b) + u1 = floor(Int, x) + u1^b == N && return (u1, b) + (u1+1)^b == N && return (u1+1, b) + end +end + +end diff --git a/lib/QuAlgorithmZoo/test/HadamardTest.jl b/lib/QuAlgorithmZoo/test/HadamardTest.jl new file mode 100644 index 000000000..d04ab038e --- /dev/null +++ b/lib/QuAlgorithmZoo/test/HadamardTest.jl @@ -0,0 +1,46 @@ +using Yao +using Test +using LinearAlgebra +using QuAlgorithmZoo +using YaoBlocks.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) + rho1 = reg1 |> ρ + reg2 = rand_state(3) |> focus!(1,2) + rho2 = reg2 |> ρ + reg3 = rand_state(3) |> focus!(1,2) + rho3 = reg3 |> ρ + desired = tr(Matrix(rho1)*Matrix(rho2)) + 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(Matrix(rho1)*Matrix(rho2)*Matrix(rho3)) |> real + c = swap_test_circuit(2, 3, 0) + res = expect(put(7, 1=>Z), reduce(join, [reg3, reg2, reg1, zero_state(1)]) |> c) |> tr |> real + @test desired ≈ res + desired = tr(Matrix(rho1)*Matrix(rho2)*Matrix(rho3)) |> imag + c = swap_test_circuit(2, 3, -π/2) + res = expect(put(7, 1=>Z), reduce(join, [reg3, reg2, reg1, zero_state(1)]) |> c) |> tr |> real + @test desired ≈ res +end + +@testset "hadamard test" begin + nbit = 4 + U = put(nbit, 2=>Rx(0.2)) + reg = rand_state(nbit) + + @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) |> YaoExtensions.singlet_block() + @test single_swap_test(reg, 0) ≈ -1 + + reg = zero_state(2) + @test single_swap_test(reg, 0) ≈ 1 + reg = product_state(2, 0b11) + @test single_swap_test(reg, 0) ≈ 1 +end diff --git a/lib/QuAlgorithmZoo/test/PhaseEstimation.jl b/lib/QuAlgorithmZoo/test/PhaseEstimation.jl new file mode 100644 index 000000000..1aa3f2713 --- /dev/null +++ b/lib/QuAlgorithmZoo/test/PhaseEstimation.jl @@ -0,0 +1,68 @@ +using Yao +using BitBasis +using Test, LinearAlgebra +using QuAlgorithmZoo + +""" +random phase estimation problem setup. +""" +function rand_phaseest_setup(N::Int) + U = rand_unitary(1< adjoint) ≈ join(reg2, reg1) +end + +@testset "phaseest, non-eigen" begin + # Generate a random matrix. + N = 3 + U, b, ϕs, evec = rand_phaseest_setup(N) + + # Define ArrayReg and U operator. + M = 6 + reg1 = zero_state(M) + reg2 = ArrayReg(b) + UG = matblock(U); + + # run circuit + reg= join(reg2, reg1) + pe = PEBlock(UG, M, N) + apply!(reg, pe) + + # measure + bs, proj, amp_relative = projection_analysis(evec, focus!(reg, M+1:M+N)) + @test isapprox(ϕs, bfloat.(bs, nbits=M), atol=0.05) +end diff --git a/lib/QuAlgorithmZoo/test/hamiltonian_solvers.jl b/lib/QuAlgorithmZoo/test/hamiltonian_solvers.jl new file mode 100644 index 000000000..740618340 --- /dev/null +++ b/lib/QuAlgorithmZoo/test/hamiltonian_solvers.jl @@ -0,0 +1,32 @@ +using Yao +using LinearAlgebra +using Test +using QuAlgorithmZoo, YaoExtensions +using YaoBlocks: ConstGate + +@testset "solving hamiltonian" begin + nbit = 8 + h = heisenberg(nbit) |> cache + @test ishermitian(h) + reg = rand_state(nqubits(h)) + E0 = expect(h, reg)/nbit + reg |> iter_groundstate!(h, niter=1000) + EG = expect(h, reg)/nbit/4 + @test isapprox(EG, -0.4564, atol=1e-4) + + # using Time Evolution + reg = rand_state(nqubits(h)) + reg |> itime_groundstate!(h, τ=20) + EG = expect(h, reg)/nbit/4 + @test isapprox(EG, -0.4564, atol=1e-4) + + # using VQE + N = 4 + h = heisenberg(N) + E = eigen(h |> mat |> Matrix).values[1] + c = YaoExtensions.variational_circuit(N, 5) + dispatch!(c, :random) + vqe_solve!(c, h) + E2 = expect(h, zero_state(N) |> c) + @test isapprox(E, E2, rtol=1e-1) +end diff --git a/lib/QuAlgorithmZoo/test/runtests.jl b/lib/QuAlgorithmZoo/test/runtests.jl new file mode 100644 index 000000000..9634f07f8 --- /dev/null +++ b/lib/QuAlgorithmZoo/test/runtests.jl @@ -0,0 +1,16 @@ +using Test, Random, LinearAlgebra, SparseArrays +using Yao +using YaoExtensions +using QuAlgorithmZoo + +@testset "PhaseEstimation" begin + include("PhaseEstimation.jl") +end + +@testset "hamiltonian solvers" begin + include("hamiltonian_solvers.jl") +end + +@testset "hadamard test" begin + include("HadamardTest.jl") +end