diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 22db846ac..ab7ed9372 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -39,6 +39,8 @@ jobs: - YaoArrayRegister - YaoBlocks - YaoSym + - YaoPlots + - YaoToEinsum steps: - uses: actions/checkout@v2 - uses: julia-actions/setup-julia@v1 diff --git a/.gitignore b/.gitignore index d79e48dcd..278ef4fb6 100644 --- a/.gitignore +++ b/.gitignore @@ -31,3 +31,5 @@ docs/src/assets/*.pdf *.synctex.gz *.fls *.log +lib/YaoPlots/examples/*.png +lib/YaoPlots/examples/*.svg diff --git a/Makefile b/Makefile new file mode 100644 index 000000000..e412c3289 --- /dev/null +++ b/Makefile @@ -0,0 +1,37 @@ +JL = julia --project + +default: init test + +init: + $(JL) -e 'using Pkg; Pkg.develop([Pkg.PackageSpec(path = joinpath("lib", pkg)) for pkg in ["YaoArrayRegister", "YaoBlocks", "YaoSym", "YaoPlots", "YaoAPI", "YaoToEinsum"]]); Pkg.precompile()' +init-docs: + $(JL) -e 'using Pkg; Pkg.activate("docs"); Pkg.develop([Pkg.PackageSpec(path = joinpath("lib", pkg)) for pkg in ["YaoArrayRegister", "YaoBlocks", "YaoSym", "YaoPlots", "YaoAPI", "YaoToEinsum"]]); Pkg.precompile()' + +update: + $(JL) -e 'using Pkg; Pkg.update(); Pkg.precompile()' +update-docs: + $(JL) -e 'using Pkg; Pkg.activate("docs"); Pkg.update(); Pkg.precompile()' + +test-Yao: + $(JL) -e 'using Pkg; Pkg.test("Yao")' +test-CuYao: + $(JL) -e 'using Pkg; Pkg.activate("ext/CuYao/test"); Pkg.develop(path="."); Pkg.update()' + $(JL) -e 'using Pkg; Pkg.activate("ext/CuYao/test"); include("ext/CuYao/test/runtests.jl")' + +test-%: + $(JL) -e 'using Pkg; Pkg.test("$*")' + +test: + $(JL) -e 'using Pkg; Pkg.test(["YaoAPI", "YaoArrayRegister", "YaoBlocks", "YaoSym", "YaoPlots", "Yao", "YaoToEinsum"])' + +coverage: + $(JL) -e 'using Pkg; Pkg.test(["YaoAPI", "YaoArrayRegister", "YaoBlocks", "YaoSym", "YaoPlots", "Yao", "YaoToEinsum"]; coverage=true)' + +servedocs: + $(JL) -e 'using Pkg; Pkg.activate("docs"); using LiveServer; servedocs(;skip_dirs=["docs/src/assets", "docs/src/generated"])' + +clean: + rm -rf docs/build + find . -name "*.cov" -type f -print0 | xargs -0 /bin/rm -f + +.PHONY: init test diff --git a/Project.toml b/Project.toml index 40f9132dc..f1a7ffd98 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "Yao" uuid = "5872b779-8223-5990-8dd0-5abbb0748c8c" -version = "0.8.8" +version = "0.9.0" [deps] BitBasis = "50ba71b6-fa0f-514d-ae9a-0916efc90dcf" @@ -10,16 +10,28 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69" YaoAPI = "0843a435-28de-4971-9e8b-a9641b2983a8" YaoArrayRegister = "e600142f-9330-5003-8abb-0ebd767abc51" YaoBlocks = "418bc28f-b43b-5e0b-a6e7-61bbc1a2c1df" +YaoPlots = "32cfe2d9-419e-45f2-8191-2267705d8dbc" YaoSym = "3b27209a-d3d6-11e9-3c0f-41eb92b2cb9d" +YaoToEinsum = "9b173c7b-dc24-4dc5-a0e1-adab2f7b6ba9" + +[weakdeps] +CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" + +[extensions] +CuYao = "CUDA" [compat] -BitBasis = "0.8" +BitBasis = "0.8, 0.9" +CUDA = "4, 5" +LinearAlgebra = "1" LuxurySparse = "0.7" Reexport = "1" YaoAPI = "0.4" YaoArrayRegister = "0.9" YaoBlocks = "0.13" +YaoPlots = "0.9" YaoSym = "0.6" +YaoToEinsum = "0.2" julia = "1" [extras] diff --git a/docs/Project.toml b/docs/Project.toml index d6d7fde09..d0ab56514 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -7,6 +7,7 @@ FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" +LiveServer = "16fef848-5104-11e9-1b77-fb7a48bbb589" Optimisers = "3bd65402-5787-11e9-1adc-39752487f4e2" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0" diff --git a/docs/make.jl b/docs/make.jl index 481653a2c..7ca1b4839 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -2,7 +2,7 @@ using Documenter using DocThemeIndigo using Literate using Yao -using Yao: YaoBlocks, YaoArrayRegister, YaoSym, BitBasis, YaoAPI +using Yao: YaoBlocks, YaoArrayRegister, YaoSym, BitBasis, YaoAPI, YaoPlots using YaoBlocks: AD using YaoBlocks: Optimise @@ -15,14 +15,9 @@ attach_notebook_badge(root, name) = str->attach_notebook_badge(root, name, str) function attach_notebook_badge(root, name, str) mybinder_badge_url = "https://mybinder.org/badge_logo.svg" - nbviewer_badge_url = "https://img.shields.io/badge/show-nbviewer-579ACA.svg" - download_badge_url = "https://img.shields.io/badge/download-project-orange" mybinder = "[![]($mybinder_badge_url)](@__BINDER_ROOT_URL__/generated/$root/$name/main.ipynb)" - nbviewer = "[![]($nbviewer_badge_url)](@__NBVIEWER_ROOT_URL__/generated/$root/$name/main.ipynb)" - download = "[![]($download_badge_url)](https://minhaskamal.github.io/DownGit/#/home?url=https://github.com/QuantumBFS/tutorials/tree/gh-pages/dev/generated/$root/$name)" - markdown_only(x) = "#md # " * x - return join(map(markdown_only, (mybinder, nbviewer, download)), "\n") * "\n\n" * str + return join(map(markdown_only, (mybinder, )), "\n") * "\n\n" * str end function build_tutorial(root, name) @@ -30,9 +25,7 @@ function build_tutorial(root, name) generated_abspath = joinpath(@__DIR__, "src", generated_path) source_dir = joinpath(@__DIR__, "src", root, name) source_path = joinpath(source_dir, "main.jl") - Literate.markdown(source_path, generated_abspath; execute=true, name="index", preprocess = attach_notebook_badge(root, name)) - Literate.notebook(source_path, generated_abspath; execute=false, name="main", preprocess = notebook_filter) - Literate.script(source_path, generated_abspath; execute=false, name="main") + Literate.markdown(source_path, generated_abspath; execute=true, name="index", preprocess = attach_notebook_badge(root, name), binder_root_url="") # copy other things for each in readdir(source_dir) if each != "main.jl" @@ -53,6 +46,8 @@ end # download("yaoquantum.org/assets/logo-light.png", "docs/src/assets/logo.png") +example_pages = build("examples") + const PAGES = [ "Home" => "index.md", "Quick Start" => "quick-start.md", @@ -63,18 +58,21 @@ const PAGES = [ "man/registers.md", "man/blocks.md", "man/symbolic.md", + "man/cuda.md", + "man/plot.md", "man/automatic_differentiation.md", + "man/yao2einsum.md", "man/simplification.md", "man/bitbasis.md", ], - "Examples" => build("examples"), + "Examples" => example_pages, "Performance Tips" => "performancetips.md", ] indigo = DocThemeIndigo.install(Yao) makedocs( - modules = [Yao, YaoAPI, YaoArrayRegister, YaoBlocks, BitBasis, YaoSym, AD, Optimise], + modules = [Yao, YaoAPI, YaoArrayRegister, YaoBlocks, BitBasis, YaoSym, YaoPlots, YaoToEinsum, AD, Optimise], format = Documenter.HTML( prettyurls = ("deploy" in ARGS), canonical = ("deploy" in ARGS) ? "https://docs.yaoquantum.org/" : nothing, @@ -88,6 +86,7 @@ makedocs( sitename = "Documentation | Yao", linkcheck = !("skiplinks" in ARGS), pages = PAGES, + warnonly = [:missing_docs, :linkcheck, :cross_references], ) diff --git a/docs/src/assets/images/Riemannian.png b/docs/src/assets/images/Riemannian.png new file mode 100644 index 000000000..7a46dc10d Binary files /dev/null and b/docs/src/assets/images/Riemannian.png differ diff --git a/docs/src/examples/1.prepare-ghz-state/index.md b/docs/src/examples/1.prepare-ghz-state/index.md deleted file mode 100644 index ada7898f4..000000000 --- a/docs/src/examples/1.prepare-ghz-state/index.md +++ /dev/null @@ -1,468 +0,0 @@ -```@meta -EditURL = "/docs/src/quick-start/1.prepare-ghz-state/main.jl" -``` - -[![](https://mybinder.org/badge_logo.svg)](/generated//home/roger/code/julia/Yao/docs/src/quick-start/1.prepare-ghz-state/main.ipynb) -[![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](/generated//home/roger/code/julia/Yao/docs/src/quick-start/1.prepare-ghz-state/main.ipynb) -[![](https://img.shields.io/badge/download-project-orange)](https://minhaskamal.github.io/DownGit/#/home?url=https://github.com/QuantumBFS/tutorials/tree/gh-pages/dev/generated//home/roger/code/julia/Yao/docs/src/quick-start/1.prepare-ghz-state) - -# [Prepare Greenberger–Horne–Zeilinger state with Quantum Circuit](@id example-ghz) - -First, you have to use this package in Julia. - -````julia -using Yao -```` - -Now, we just define the circuit according to the circuit image below: -![ghz](assets/ghz4.png) - -````julia -circuit = chain( - 4, - put(1=>X), - repeat(H, 2:4), - control(2, 1=>X), - control(4, 3=>X), - control(3, 1=>X), - control(4, 3=>X), - repeat(H, 1:4), -) -```` - -```` -nqubits: 4 -chain -├─ put on (1) -│ └─ X -├─ repeat on (2, 3, 4) -│ └─ H -├─ control(2) -│ └─ (1,) X -├─ control(4) -│ └─ (3,) X -├─ control(3) -│ └─ (1,) X -├─ control(4) -│ └─ (3,) X -└─ repeat on (1, 2, 3, 4) - └─ H - -```` - -Let me explain what happens here. - -## Put single qubit gate X to location 1 -we have an `X` gate applied to the first qubit. -We need to tell `Yao` to put this gate on the first qubit by - -````julia -put(4, 1=>X) -```` - -```` -nqubits: 4 -put on (1) -└─ X -```` - -We use Julia's `Pair` to denote the gate and its location in the circuit, -for two-qubit gate, you could also use a tuple of locations: - -````julia -put(4, (1, 2)=>swap(2, 1, 2)) -```` - -```` -nqubits: 4 -put on (1, 2) -└─ put on (1, 2) - └─ SWAP - -```` - -But, wait, why there's no `4` in the definition above? This is because -all the functions in `Yao` that requires to input the number of qubits as its -first argument could be lazy (curried), and let other constructors to infer the total -number of qubits later, e.g - -````julia -put(1=>X) -```` - -```` -(n -> put(n, 1 => X)) -```` - -which will return a lambda that ask for a single argument `n`. - -````julia -put(1=>X)(4) -```` - -```` -nqubits: 4 -put on (1) -└─ X -```` - -## Apply the same gate on different locations - -next we should put Hadmard gates on all locations except the 1st qubits. - -We provide `repeat` to apply the same block repeatly, repeat can take an -iterator of desired locations, and like `put`, we can also leave the total number -of qubits there. - -````julia -repeat(H, 2:4) -```` - -```` -(n -> repeat(n, H, 2:4...)) -```` - -## Define control gates - -In Yao, we could define controlled gates by feeding a gate to `control` - -````julia -control(4, 2, 1=>X) -```` - -```` -nqubits: 4 -control(2) -└─ (1,) X -```` - -Like many others, you could leave the number of total qubits there, and infer it -later. - -````julia -control(2, 1=>X) -```` - -```` -(n -> control(n, 2, 1 => X)) -```` - -## Composite each part together - -This will create a `ControlBlock`, the concept of block in Yao basically -just means quantum operators, since the quantum circuit itself is a quantum operator, -we could create a quantum circuit by composite each part of. - -Here, we use `chain` to chain each part together, a chain of quantum operators -means to apply each operators one by one in the chain. This will create a `ChainBlock`. - -````julia -circuit = chain( - 4, - put(1=>X), - repeat(H, 2:4), - control(2, 1=>X), - control(4, 3=>X), - control(3, 1=>X), - control(4, 3=>X), - repeat(H, 1:4), -) -```` - -```` -nqubits: 4 -chain -├─ put on (1) -│ └─ X -├─ repeat on (2, 3, 4) -│ └─ H -├─ control(2) -│ └─ (1,) X -├─ control(4) -│ └─ (3,) X -├─ control(3) -│ └─ (1,) X -├─ control(4) -│ └─ (3,) X -└─ repeat on (1, 2, 3, 4) - └─ H - -```` - -You can check the type of it with `typeof` - -````julia -typeof(circuit) -```` - -```` -ChainBlock{4} -```` - -## Construct GHZ state from 00...00 - -For simulation, we provide a builtin register type called `ArrayReg`, -we will use the simulated register in this example. - -First, let's create ``|00⋯00⟩``, you can create it with `zero_state` - -````julia -zero_state(4) -```` - -```` -ArrayReg{1, ComplexF64, Array...} - active qubits: 4/4 -```` - -Or we also provide bit string literals to create arbitrary eigen state - -````julia -ArrayReg(bit"0000") -```` - -```` -ArrayReg{1, ComplexF64, Array...} - active qubits: 4/4 -```` - -They will both create a register with Julia's builtin `Array` as storage. - -## Feed Registers to Circuits - -Circuits can be applied to registers with `apply!` - -````julia -apply!(zero_state(4), circuit) -```` - -```` -ArrayReg{1, ComplexF64, Array...} - active qubits: 4/4 -```` - -or you can use pipe operator `|>`, when you want to chain several operations -together, here we measure the state right after the circuit for `1000` times - -````julia -results = zero_state(4) |> circuit |> r->measure(r, nshots=1000) - -using StatsBase, Plots - -hist = fit(Histogram, Int.(results), 0:16) -bar(hist.edges[1] .- 0.5, hist.weights, legend=:none) -```` - -```@raw html - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -``` - -GHZ state will collapse to ``|0000⟩`` or ``|1111⟩``. - ---- - -*This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).* - diff --git a/docs/src/examples/1.prepare-ghz-state/main.ipynb b/docs/src/examples/1.prepare-ghz-state/main.ipynb deleted file mode 100644 index b1056b625..000000000 --- a/docs/src/examples/1.prepare-ghz-state/main.ipynb +++ /dev/null @@ -1,390 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "source": [ - "# Prepare Greenberger–Horne–Zeilinger state with Quantum Circuit" - ], - "metadata": {} - }, - { - "cell_type": "markdown", - "source": [ - "First, you have to use this package in Julia." - ], - "metadata": {} - }, - { - "outputs": [], - "cell_type": "code", - "source": [ - "using Yao" - ], - "metadata": {}, - "execution_count": null - }, - { - "cell_type": "markdown", - "source": [ - "Now, we just define the circuit according to the circuit image below:\n", - "![ghz](assets/ghz4.png)" - ], - "metadata": {} - }, - { - "outputs": [], - "cell_type": "code", - "source": [ - "circuit = chain(\n", - " 4,\n", - " put(1=>X),\n", - " repeat(H, 2:4),\n", - " control(2, 1=>X),\n", - " control(4, 3=>X),\n", - " control(3, 1=>X),\n", - " control(4, 3=>X),\n", - " repeat(H, 1:4),\n", - ")" - ], - "metadata": {}, - "execution_count": null - }, - { - "cell_type": "markdown", - "source": [ - "Let me explain what happens here." - ], - "metadata": {} - }, - { - "cell_type": "markdown", - "source": [ - "## Put single qubit gate X to location 1\n", - "we have an `X` gate applied to the first qubit.\n", - "We need to tell `Yao` to put this gate on the first qubit by" - ], - "metadata": {} - }, - { - "outputs": [], - "cell_type": "code", - "source": [ - "put(4, 1=>X)" - ], - "metadata": {}, - "execution_count": null - }, - { - "cell_type": "markdown", - "source": [ - "We use Julia's `Pair` to denote the gate and its location in the circuit,\n", - "for two-qubit gate, you could also use a tuple of locations:" - ], - "metadata": {} - }, - { - "outputs": [], - "cell_type": "code", - "source": [ - "put(4, (1, 2)=>swap(2, 1, 2))" - ], - "metadata": {}, - "execution_count": null - }, - { - "cell_type": "markdown", - "source": [ - "But, wait, why there's no `4` in the definition above? This is because\n", - "all the functions in `Yao` that requires to input the number of qubits as its\n", - "first arguement could be lazy (curried), and let other constructors to infer the total\n", - "number of qubits later, e.g" - ], - "metadata": {} - }, - { - "outputs": [], - "cell_type": "code", - "source": [ - "put(1=>X)" - ], - "metadata": {}, - "execution_count": null - }, - { - "cell_type": "markdown", - "source": [ - "which will return a lambda that ask for a single arguement `n`." - ], - "metadata": {} - }, - { - "outputs": [], - "cell_type": "code", - "source": [ - "put(1=>X)(4)" - ], - "metadata": {}, - "execution_count": null - }, - { - "cell_type": "markdown", - "source": [ - "## Apply the same gate on different locations" - ], - "metadata": {} - }, - { - "cell_type": "markdown", - "source": [ - "next we should put Hadmard gates on all locations except the 1st qubits." - ], - "metadata": {} - }, - { - "cell_type": "markdown", - "source": [ - "We provide `repeat` to apply the same block repeatly, repeat can take an\n", - "iterator of desired locations, and like `put`, we can also leave the total number\n", - "of qubits there." - ], - "metadata": {} - }, - { - "outputs": [], - "cell_type": "code", - "source": [ - "repeat(H, 2:4)" - ], - "metadata": {}, - "execution_count": null - }, - { - "cell_type": "markdown", - "source": [ - "## Define control gates" - ], - "metadata": {} - }, - { - "cell_type": "markdown", - "source": [ - "In Yao, we could define controlled gates by feeding a gate to `control`" - ], - "metadata": {} - }, - { - "outputs": [], - "cell_type": "code", - "source": [ - "control(4, 2, 1=>X)" - ], - "metadata": {}, - "execution_count": null - }, - { - "cell_type": "markdown", - "source": [ - "Like many others, you could leave the number of total qubits there, and infer it\n", - "later." - ], - "metadata": {} - }, - { - "outputs": [], - "cell_type": "code", - "source": [ - "control(2, 1=>X)" - ], - "metadata": {}, - "execution_count": null - }, - { - "cell_type": "markdown", - "source": [ - "## Composite each part together" - ], - "metadata": {} - }, - { - "cell_type": "markdown", - "source": [ - "This will create a `ControlBlock`, the concept of block in Yao basically\n", - "just means quantum operators, since the quantum circuit itself is a quantum operator,\n", - "we could create a quantum circuit by composite each part of." - ], - "metadata": {} - }, - { - "cell_type": "markdown", - "source": [ - "Here, we use `chain` to chain each part together, a chain of quantum operators\n", - "means to apply each operators one by one in the chain. This will create a `ChainBlock`." - ], - "metadata": {} - }, - { - "outputs": [], - "cell_type": "code", - "source": [ - "circuit = chain(\n", - " 4,\n", - " put(1=>X),\n", - " repeat(H, 2:4),\n", - " control(2, 1=>X),\n", - " control(4, 3=>X),\n", - " control(3, 1=>X),\n", - " control(4, 3=>X),\n", - " repeat(H, 1:4),\n", - ")" - ], - "metadata": {}, - "execution_count": null - }, - { - "cell_type": "markdown", - "source": [ - "You can check the type of it with `typeof`" - ], - "metadata": {} - }, - { - "outputs": [], - "cell_type": "code", - "source": [ - "typeof(circuit)" - ], - "metadata": {}, - "execution_count": null - }, - { - "cell_type": "markdown", - "source": [ - "## Construct GHZ state from 00...00" - ], - "metadata": {} - }, - { - "cell_type": "markdown", - "source": [ - "For simulation, we provide a builtin register type called `ArrayReg`,\n", - "we will use the simulated register in this example." - ], - "metadata": {} - }, - { - "cell_type": "markdown", - "source": [ - "First, let's create ``|00⋯00⟩``, you can create it with `zero_state`" - ], - "metadata": {} - }, - { - "outputs": [], - "cell_type": "code", - "source": [ - "zero_state(4)" - ], - "metadata": {}, - "execution_count": null - }, - { - "cell_type": "markdown", - "source": [ - "Or we also provide bit string literals to create arbitrary eigen state" - ], - "metadata": {} - }, - { - "outputs": [], - "cell_type": "code", - "source": [ - "ArrayReg(bit\"0000\")" - ], - "metadata": {}, - "execution_count": null - }, - { - "cell_type": "markdown", - "source": [ - "They will both create a register with Julia's builtin `Array` as storage." - ], - "metadata": {} - }, - { - "cell_type": "markdown", - "source": [ - "## Feed Registers to Circuits" - ], - "metadata": {} - }, - { - "cell_type": "markdown", - "source": [ - "Circuits can be applied to registers with `apply!`" - ], - "metadata": {} - }, - { - "outputs": [], - "cell_type": "code", - "source": [ - "apply!(zero_state(4), circuit)" - ], - "metadata": {}, - "execution_count": null - }, - { - "cell_type": "markdown", - "source": [ - "or you can use pipe operator `|>`, when you want to chain several operations\n", - "together, here we measure the state right after the circuit for `1000` times" - ], - "metadata": {} - }, - { - "outputs": [], - "cell_type": "code", - "source": [ - "results = zero_state(4) |> circuit |> r->measure(r, nshots=1000)\n", - "\n", - "using StatsBase, Plots\n", - "\n", - "hist = fit(Histogram, Int.(results), 0:16)\n", - "bar(hist.edges[1] .- 0.5, hist.weights, legend=:none)" - ], - "metadata": {}, - "execution_count": null - }, - { - "cell_type": "markdown", - "source": [ - "GHZ state will collapse to ``|0000⟩`` or ``|1111⟩``." - ], - "metadata": {} - }, - { - "cell_type": "markdown", - "source": [ - "---\n", - "\n", - "*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*" - ], - "metadata": {} - } - ], - "nbformat_minor": 3, - "metadata": { - "language_info": { - "file_extension": ".jl", - "mimetype": "application/julia", - "name": "julia", - "version": "1.7.0" - }, - "kernelspec": { - "name": "julia-1.7", - "display_name": "Julia 1.7.0", - "language": "julia" - } - }, - "nbformat": 4 -} diff --git a/docs/src/examples/2.qft-phase-estimation/main.jl b/docs/src/examples/2.qft-phase-estimation/main.jl index 426ddbdf8..50c4e9173 100644 --- a/docs/src/examples/2.qft-phase-estimation/main.jl +++ b/docs/src/examples/2.qft-phase-estimation/main.jl @@ -32,7 +32,7 @@ A(i, j) = control(i, j=>shift(2π/(1<<(i-j+1)))) R4 = A(4, 1) -# If you have read about [preparing GHZ state](@ref example-ghz), +# If you have read about [preparing GHZ state](@ref tutorial-ghz), # you probably know that in Yao, we could just leave the number of qubits, and it # will be evaluated when possible. diff --git a/docs/src/examples/6.quantum-circuit-born-machine/main.jl b/docs/src/examples/6.quantum-circuit-born-machine/main.jl index aa07bd167..3c8b2f6c1 100644 --- a/docs/src/examples/6.quantum-circuit-born-machine/main.jl +++ b/docs/src/examples/6.quantum-circuit-born-machine/main.jl @@ -27,7 +27,7 @@ pg = gaussian_pdf(1:1<<6, 1<<5-0.5, 1<<4); # We can plot the distribution, it looks like -plot(pg) +Plots.plot(pg) # ## Create the Circuit @@ -225,12 +225,12 @@ trained_pg = probs(zero_state(nqubits(qcbm)) |> qcbm) title!("training history") xlabel!("steps"); ylabel!("loss") -plot(history) +Plots.plot(history) # And let's check what we got -fig2 = plot(1:1<<6, trained_pg; label="trained") -plot!(fig2, 1:1<<6, pg; label="target") +fig2 = Plots.plot(1:1<<6, trained_pg; label="trained") +Plots.plot!(fig2, 1:1<<6, pg; label="target") title!("distribution") xlabel!("x"); ylabel!("p") diff --git a/docs/src/examples/8.riemannian-gradient-flow/main.jl b/docs/src/examples/8.riemannian-gradient-flow/main.jl index e1d4ec7a2..93d1bcf6f 100644 --- a/docs/src/examples/8.riemannian-gradient-flow/main.jl +++ b/docs/src/examples/8.riemannian-gradient-flow/main.jl @@ -42,8 +42,8 @@ for i in 1:100 push!(history, real.(expect(h, zero_state(n)=>circuit))) end -plot(history, legend=false) -plot!(1:100, [w[1] for i=1:100]) +Plots.plot(history, legend=false) +Plots.plot!(1:100, [w[1] for i=1:100]) xlabel!("steps") ylabel!("energy") @@ -125,20 +125,25 @@ end; # Finally, let's try it out. # We initialize the state ``|0\rangle`` and apply several optimization steps. -circuit = chain(n) -pauli_strings = generate_2local_pauli_strings(n) -history = Float64[] +# ```julia +# circuit = chain(n) +# pauli_strings = generate_2local_pauli_strings(n) +# history = Float64[] -for i=1:100 - cost = step_and_cost!(n, circuit, h, 0.01, pauli_strings) - push!(history, cost) -end +# for i=1:100 +# cost = step_and_cost!(n, circuit, h, 0.01, pauli_strings) +# push!(history, cost) +# end -plot(history, legend=false) -plot!(1:100, [w[1] for i=1:100]) -xlabel!("steps") -ylabel!("energy") +# Plots.plot(history, legend=false) +# Plots.plot!(1:100, [w[1] for i=1:100]) +# xlabel!("steps") +# ylabel!("energy") +# ``` +# ```@raw html +# Riemannian gradient flow +# ``` # When we compare the final states achieved with the Riemannian gradient flow # optimizer and with the standard VQE we can notice that the former has lower quality. # This is because the Riemannian gradient flow optimizer has only a local view of the cost landscape diff --git a/docs/src/index.md b/docs/src/index.md index 3222982ac..38a98cb29 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -26,7 +26,7 @@ The tutorials are written with [Literate](https://github.com/fredrikekre/Literat ## Pluto Notebooks -There is also a set of Pluto notebooks available in [the notebooks folder](https://github.com/QuantumBFS/Yao.jl/notebooks) +There is also a set of Pluto notebooks available in [the notebooks folder](notebooks) ## Manual diff --git a/docs/src/man/bitbasis.md b/docs/src/man/bitbasis.md index 22afa3956..f36695fcc 100644 --- a/docs/src/man/bitbasis.md +++ b/docs/src/man/bitbasis.md @@ -1,4 +1,5 @@ ```@meta +CurrentModule = BitBasis DocTestSetup = quote using Yao using Yao: YaoBlocks, YaoArrayRegister @@ -16,4 +17,4 @@ For more details please ref to [BitBasis.jl](https://quantumbfs.github.io/BitBas ```@autodocs Modules = [BitBasis] Order = [:macro, :function, :type, :module] -``` +``` \ No newline at end of file diff --git a/docs/src/man/blocks.md b/docs/src/man/blocks.md index 7f5ced0e4..1a38ba412 100644 --- a/docs/src/man/blocks.md +++ b/docs/src/man/blocks.md @@ -1,4 +1,5 @@ ```@meta +CurrentModule = YaoBlocks DocTestSetup = quote using Yao using Yao: YaoBlocks, YaoArrayRegister @@ -10,13 +11,19 @@ end # Blocks **Blocks** are the basic building blocks of a quantum circuit in Yao. -It simply means a quantum operator, thus, all the blocks have matrices in principal and one can get its matrix by [`mat`](@ref). The basic blocks required to build an arbitrary quantum circuit is defined in the component package [`YaoBlocks`](@ref). +It simply means a quantum operator, thus, all the blocks have matrices in principal and one can get its matrix by [`mat`](@ref). The basic blocks required to build an arbitrary quantum circuit is defined in the component package `YaoBlocks`. Block Tree serves as an intermediate representation for Yao to analysis, optimize the circuit, then it will be lowered to instructions like for simulations, blocks will be lowered to [`instruct!`](@ref) calls. The structure of blocks is the same with a small type system, it consists of two basic kinds of blocks: [`CompositeBlock`](@ref) (like composite types), and [`PrimitiveBlock`](@ref) (like primitive types). By combining these two kinds of blocks together, we'll be able to construct a quantum circuit and represent it in a tree data structure. +```@docs +AbstractBlock +PrimitiveBlock +CompositeBlock +``` + ## Primitive Blocks Primitive blocks are subtypes of [`PrimitiveBlock`](@ref), they are the leaf nodes in a block tree, thus primitive types do not have subtypes. @@ -28,6 +35,10 @@ Modules = [YaoBlocks] Filter = t ->(t isa Type && t <: YaoBlocks.PrimitiveBlock) ``` +```@docs +@const_gate +``` + ## Composite Blocks Composite blocks are subtypes of [`CompositeBlock`](@ref), they are the composition of blocks. @@ -39,6 +50,24 @@ Modules = [YaoBlocks] Filter = t -> t isa Type && t <: YaoBlocks.CompositeBlock ``` +## Operations on Blocks +```@docs +unsafe_apply! +apply! +apply + +niparams +getiparams +render_params +nparameters + +occupied_locs +print_block + +apply_back! +mat_back! +``` + ## Error and Exceptions ```@autodocs diff --git a/docs/src/man/cuda.md b/docs/src/man/cuda.md new file mode 100644 index 000000000..98454f1b2 --- /dev/null +++ b/docs/src/man/cuda.md @@ -0,0 +1,71 @@ +```@meta +CurrentModule = Yao +``` + +# CUDA extension - CuYao + +## Tutorial +`CuYao` is a CUDA extension of Yao, which allows you to run Yao circuits on GPU. The usage of `CuYao` is similar to `Yao`, but with some extra APIs to upload and download registers between CPU and GPU: +- `cu(reg)` to upload a registe `reg` to GPU, and +- `cpu(cureg)` to download a register `cureg` from GPU to CPU. + +```julia +julia> using Yao, CUDA + +# create a register on GPU +julia> cureg = rand_state(9; nbatch=1000) |> cu; # or `curand_state(9; nbatch=1000)`. + +# run a circuit on GPU +julia> cureg |> put(9, 2=>Z); + +# measure the register on GPU +julia> measure!(cureg) +1000-element CuArray{DitStr{2, 9, Int64}, 1, CUDA.Mem.DeviceBuffer}: + 110110100 ₍₂₎ + 000100001 ₍₂₎ + 111111001 ₍₂₎ + ⋮ + 010001101 ₍₂₎ + 000100110 ₍₂₎ + +# download the register to CPU +julia> reg = cureg |> cpu; +``` + + +## Features +Supported gates: + +- general U(N) gate +- general U(1) gate +- X, Y, Z gate +- T, S gate +- SWAP gate +- control gates + +Supported register operations: + +- measure!, measure_reset!, measure_remove!, select +- append_qudits!, append_qubits! +- insert_qudit!, insert_qubits! +- focus!, relax! +- join +- density_matrix +- fidelity +- expect + +Autodiff: +- autodiff is supported when the only parameterized gates are rotation gates in a circuit. + +## API +```@docs +cpu +curand_state +cuzero_state +cuproduct_state +cuuniform_state +cughz_state +``` + +!!! note + the `cu` function is not documented in this module, but it is used to upload a register to GPU. \ No newline at end of file diff --git a/docs/src/man/plot.md b/docs/src/man/plot.md new file mode 100644 index 000000000..7e1e91606 --- /dev/null +++ b/docs/src/man/plot.md @@ -0,0 +1,86 @@ +```@meta +CurrentModule = YaoPlots +DocTestSetup = quote + using Yao + using Yao: YaoBlocks, YaoArrayRegister, YaoSym + using YaoBlocks + using YaoArrayRegister + using YaoPlots +end +``` + +# Quantum Circuit Visualization + +`YaoPlots` is the Quantum circuit visualization component for Yao. + + +## Examples +#### Example 1: Visualize a QBIR define in Yao + +```@example plot +using Yao.EasyBuild, YaoPlots + +# show a qft circuit +vizcircuit(qft_circuit(5)) +``` + +If you are using a Pluto/Jupyter notebook, Atom/VSCode editor, you should see the above image in your plotting panel. + +#### Example 2: Visualize a single qubit +```@example plot +using YaoPlots, Yao + +reg = zero_state(1) |> Rx(π/8) |> Rx(π/8) +rho = density_matrix(ghz_state(2), 1) + +bloch_sphere("|ψ⟩"=>reg, "ρ"=>rho; show_projection_lines=true) +``` + +See more [examples](lib/YaoPlots/examples/circuits.jl). + +### Adjusting the plot attributes + +Various attributes of the visualizations can be altered. +The plot can be modified, if we change the following attributes + +- `YaoPlots.CircuitStyles.linecolor[]` for line color, default value being `"#000000"` (black color) +- `YaoPlots.CircuitStyles.gate_bgcolor[]` for background color of square blocks, the default value being `"#FFFFFF"` (white color) +- `YaoPlots.CircuitStyles.textcolor[]` for text color, default value being `"#000000` +- `YaoPlots.CircuitStyles.lw[]` for line width, default value being `1` (pt) +- `YaoPlots.CircuitStyles.textsize[]` for text size, default value being `16` (pt) +- `YaoPlots.CircuitStyles.paramtextsize[]` for parameter text size, for parameterized gates, default value being `10` (pt) + +For example, + +```@example plot +using YaoPlots, Yao +YaoPlots.CircuitStyles.linecolor[] = "pink" +YaoPlots.CircuitStyles.gate_bgcolor[] = "yellow" +YaoPlots.CircuitStyles.textcolor[] = "#000080" # the navy blue color +YaoPlots.CircuitStyles.fontfamily[] = "JuliaMono" +YaoPlots.CircuitStyles.lw[] = 2.5 +YaoPlots.CircuitStyles.textsize[] = 13 +YaoPlots.CircuitStyles.paramtextsize[] = 8 + +vizcircuit(chain(3, put(1=>X), repeat(3, H), put(2=>Y), repeat(3, Rx(π/2)))) +``` + +## Circuit Visualization +```@docs +vizcircuit +plot +``` + +## Bloch Sphere Visualization + +```@docs +CircuitStyles +bloch_sphere +BlochStyles +``` + +## Themes +```@docs +darktheme! +lighttheme! +``` diff --git a/docs/src/man/registers.md b/docs/src/man/registers.md index 29ce6a6d0..51c36f431 100644 --- a/docs/src/man/registers.md +++ b/docs/src/man/registers.md @@ -1,7 +1,9 @@ ```@meta +CurrentModule = YaoArrayRegister DocTestSetup = quote using Yao - using Yao: YaoBlocks, YaoArrayRegister + using BitBasis + using YaoAPI using YaoBlocks using YaoArrayRegister end @@ -9,9 +11,163 @@ end # [Quantum Registers](@id registers) +## Constructing quantum states + A quantum register is a quantum state or a batch of quantum states. -`Yao` provides two types of quantum registers [`ArrayReg`](@ref) and [`BatchedArrayReg`](@ref). +Qubits in a Yao register can be active or inactive. +Only active qubits are visible to quantum operators, which enables applying quantum operators on a subset of qubits. +For example, Suppose we want to run a quantum Fourier transformation circuit of size 4 on qubits `(1, 3, 5, 7)` with the [`focus!`](@ref) function, +we first set these qubits to active qubits the rest to inactive, then we apply the circuit on the active qubits, and finally we switch back to the original configuration with the [`relax!`](@ref) function. + +`Yao` provides two types of quantum registers [`ArrayReg`](@ref) and [`BatchedArrayReg`](@ref). Both use matrices as the storage. +For example, for a quantum register with ``a`` active qubits, ``r`` remaining qubits and batch size ``b``, the storage is as follows. + +![](../assets/images/regstorage.svg) + +The first dimension of size ``2^a`` is for active qubits, only this subset of qubits are allowed to interact with quantum operators. Since we reshaped the state vector into a matrix, applying a quantum operator can be conceptually represented as a matrix-matrix multiplication. + +Various quantum states can be created with the following functions. +```@repl register +using Yao +reg = ArrayReg([0, 1, -1+0.0im, 0]) # a unnormalized Bell state |01⟩ - |10⟩ +statevec(reg) # a quantum state is represented as a vector +print_table(reg) +``` + +```@repl register +reg_zero = zero_state(3) # create a zero state |000⟩ +print_table(reg_zero) +``` + +```@repl register +reg_rand = rand_state(ComplexF32, 3) # a random state +``` + +```@repl register +reg_uniform = uniform_state(ComplexF32, 3) # a uniform state +print_table(reg_uniform) +``` + +```@repl register +reg_prod = product_state(bit"110") # a product state +bit"110"[3] # the bit string is in little-endian format +print_table(reg_prod) +``` + +```@repl register +reg_ghz = ghz_state(3) # a GHZ state +print_table(reg_ghz) +von_neumann_entropy(reg_ghz, (1, 3)) / log(2) # entanglement entropy between qubits (1, 3) and (2,) +``` + +```@repl register +reg_rand3 = rand_state(3, nlevel=3) # a random qutrit state +reg_prod3 = product_state(dit"120;3") # a qudit product state, what follows ";" symbol denotes the number of levels +print_table(reg_prod3) +``` + +```@repl register +reg_batch = rand_state(3; nbatch=2) # a batch of 2 random qubit states +print_table(reg_batch) +reg_view = viewbatch(reg_batch, 1) # view the first state in the batch +print_table(reg_view) +``` + +```@repl register +reg = rand_state(3; nlevel=4, nbatch=2) +nqudits(reg) # the total number of qudits +nactive(reg) # the number of active qubits +nremain(reg) # the number of remaining qubits +nbatch(reg) # the batch size +nlevel(reg) # the number of levels of each qudit +basis(reg) # the basis of the register +focus!(reg, 1:2) # set on the first two qubits as active +nactive(reg) # the number of active qubits +basis(reg) # the basis of the register +relax!(reg) # set all qubits as active +nactive(reg) # the number of active qubits +reorder!(reg, (3,1,2)) # reorder the qubits + +reg1 = product_state(bit"111"); +reg2 = ghz_state(3); +fidelity(reg1, reg2) # the fidelity between two states +tracedist(reg1, reg2) # the trace distance between two states +``` + +## Arithmetic operations + +The list of arithmetic operations for [`ArrayReg`](@ref) include +* `+` +* `-` +* `*` +* `/` (scalar) +* `adjoint` + + +```@repl register +reg1 = rand_state(3) +reg2 = rand_state(3) +reg3 = reg1 + reg2 # addition +normalize!(reg3) # normalize the state +isnormalized(reg3) # check if the state is normalized +reg1 - reg2 # subtraction +reg1 * 2 # scalar multiplication +reg1 / 2 # scalar division +reg1' # adjoint +reg1' * reg1 # inner product +``` + +## Register operations +```@repl register +reg0 = rand_state(3) +append_qudits!(reg0, 2) # append 2 qubits +insert_qudits!(reg0, 2, 2) # insert 2 qubits at the 2nd position +``` + +Comparing with using matrix multiplication for quantum simulation, using specialized instructions are much faster and memory efficient. These instructions are specified with the [`instruct!`](@ref) function. + +```@repl register +reg = zero_state(2) +instruct!(reg, Val(:H), (1,)) # apply a Hadamard gate on the first qubit +print_table(reg) +``` + +## Measurement + +We use the [`measure!`](@ref) function returns the measurement outcome and collapses the state after the measurement. +We also have some "cheating" version [`measure`](@ref) that does not collapse states to facilitate classical simulation. + +```@repl register +measure!(reg0, 1) # measure the qubit, the state collapses +measure!(reg0) # measure all qubits +measure(reg0, 3) # measure the qubit 3 times, the state does not collapse (hacky) +reorder!(reg0, 7:-1:1) # reorder the qubits +measure!(reg0) +invorder!(reg0) # reverse the order of qubits +measure!(reg0) +measure!(RemoveMeasured(), reg0, 2:4) # remove the measured qubits +reg0 + +reg1 = ghz_state(3) +select!(reg1, bit"111") # post-select the |111⟩ state +isnormalized(reg1) # check if the state is normalized +``` + +## Density matrices + +```@repl register +reg = rand_state(3) +rho = density_matrix(reg) # the density matrix of the state +rand_density_matrix(3) # a random density matrix +completely_mixed_state(3) # a completely mixed state +partial_tr(rho, 1) # partial trace on the first qubit +purify(rho) # purify the state +von_neumann_entropy(rho) # von Neumann entropy +mutual_information(rho, 1, 2) # mutual information between qubits 1 and 2 +``` +## API +The constructors and functions for quantum registers are listed below. ```@docs AbstractRegister AbstractArrayReg @@ -19,8 +175,6 @@ ArrayReg BatchedArrayReg ``` -We define some shortcuts to create simulated quantum states easier: - ```@docs arrayreg product_state @@ -32,36 +186,22 @@ ghz_state clone ``` -In a register, qubits are distinguished as active and inactive (or remaining). -The total number of qubits is the number of active qubits plus the number of remaining qubits. -Only active qubits are visible to quantum operators and the number of these qubits are the *size* of a register. -Making this distinction of qubits allows writing reusable quantum circuits. -For example, Suppose we want to run a quantum Fourier transformation circuit of size 4 on qubits `(1, 3, 5, 7)`, -we first set the target qubits to active qubits the reset to inactive, then we apply the circuit on it, finally we unset the inactive qubits. +The following functions are for querying the properties of a quantum register. ```@docs nqudits nqubits nactive nremain -nbatch, -nlevel, +nbatch +nlevel focus! focus relax! -zero_state exchange_sysenv ``` -## Storage - -Both [`ArayReg`](@reef) and [`BatchedArrayReg`](@ref) use matrices as the storage. For example, for a quantum register with ``a`` active qubits, ``r`` remaining qubits and batch size ``b``, the storage is as follows - -![](../assets/images/regstorage.svg) - -The first dimension of size ``2^a`` is for active qubits, only this subset of qubits are allowed to interact with blocks. Since we reshaped the state vector into a matrix, applying a quantum operator can always be represented as a matrix-matrix multiplication . In practice, most gates have in-place implementation that does not require constructing the operator matrix explicitly. - -You can access different views of the storage of an [`ArrayReg`](@ref) with the following functions: +The following functions are for querying the state of a quantum register. ```@docs state @@ -74,23 +214,7 @@ viewbatch transpose_storage ``` -## Operations - -The list of arithmetic operations for [`ArrayReg`](@ref) include -* `+` -* `-` -* `*` -* `/` (scalar) -* `adjoint` - -Then the inner product can be computed as follows. - -```julia -julia> reg = rand_state(3); - -julia> reg' * reg -0.9999999999999998 + 0.0im -``` +The following functions are for arithmetic operations on quantum registers. ```@docs AdjointArrayReg @@ -98,9 +222,9 @@ AdjointArrayReg We also have some faster inplace versions of arithematic operations ```@docs -regadd!, -regsub!, -regscale!, +regadd! +regsub! +regscale! ``` We also define the following functions for state normalization, and distance measurement. @@ -111,33 +235,24 @@ fidelity tracedist ``` -## Resource management and addressing +The following functions are for adding and reordering qubits in a quantum register. ```@docs -add_qudits! -add_qubits! +insert_qudits! +insert_qubits! append_qudits! append_qubits! reorder! invorder! ``` -Only a subset of qubits that does not interact with other qubits can be removed, the best approach is first measuring it in computational basis first. -It can be done with the [`measure!`](@ref) function by setting the first argument to `RemoveMeasured()`. - -## Instruction set - -Although we have matrix representation for Yao blocks, specialized instructions are much faster and memory efficient than using the matrix-matrix product. -These instructions are specified with the `instruct!` function listed bellow. +The `instruct!` function is for applying quantum operators on a quantum register. ```@docs YaoArrayRegister.instruct! ``` -## Measurement - -We have a true measure function `measure!` that collapses the state after the measurement. -We also have some "cheating" functions to facilitate classical simulation. +The following functions are for measurement and post-selection. ```@docs measure! @@ -146,10 +261,10 @@ select! select collapseto! probs -most_probable, +most_probable ``` -## Density matrices +The following functions are for density matrices. ```@docs DensityMatrix @@ -158,6 +273,6 @@ rand_density_matrix completely_mixed_state partial_tr purify -von_neumann_entropy, -mutual_information, +von_neumann_entropy +mutual_information ``` diff --git a/docs/src/man/symbolic.md b/docs/src/man/symbolic.md index 30743c927..b8e67120e 100644 --- a/docs/src/man/symbolic.md +++ b/docs/src/man/symbolic.md @@ -11,11 +11,25 @@ end # Symbolic Computation -Symbolic Computation support for Yao +The symbolic engine of Yao is based on [SymEngine.jl](https://github.com/symengine/SymEngine.jl). It allows one to define quantum circuits with symbolic parameters and perform symbolic computation on them. Two macro/functions play a key role in the symbolic computation: +- `@vars` for defining symbolic variables +- `subs` for substituting symbolic variables with concrete values + +```@repl sym +using Yao +@vars θ +circuit = chain(2, put(1=>H), put(2=>Ry(θ))) +mat(circuit) +new_circuit = subs(circuit, θ=>π/2) +mat(new_circuit) +``` + +## API + +The following functions are for working with symbolic states. ```@docs @ket_str @bra_str -@vars szero_state ``` diff --git a/docs/src/man/yao2einsum.md b/docs/src/man/yao2einsum.md new file mode 100644 index 000000000..f8cce1f6e --- /dev/null +++ b/docs/src/man/yao2einsum.md @@ -0,0 +1,45 @@ +# Tensor network backend + +Simulating quantum circuits using tensor networks has been studied in the literature[^Markov2008][^Pan2022]. The `YaoToEinsum` package provides a convenient way to convert Yao circuits to tensor networks, which can be used for further analysis and optimization. + +## Tutorial +The main function is +```julia +yao2einsum(circuit; initial_state=Dict(), final_state=Dict(), optimizer=TreeSA()) +``` +which transforms a [`Yao`](https://github.com/QuantumBFS/Yao.jl) circuit to a tensor network that generalizes the hyper-graph (einsum notation). The return value is a `TensorNetwork` object. + +* `initial_state` and `final_state` are for specifying the initial state and final state. Left the qubits unspecified if you want to keep them as the open indices. +* `optimizer` is for specifying the contraction order optimizing algorithm of the tensor network. The default value is the `TreeSA()` algorithm that developed in [^Kalachev2021][^Liu2023]. Please check the README of [OMEinsumEinsumContractors.jl](https://github.com/TensorBFS/OMEinsumContractionOrders.jl) for more information. + +In the following example, we show how to convert a quantum Fourier transform circuit to a tensor network and contract it to +- Get the matrix representation of the circuit. +- Get the probability of measuring the zero state after applying the circuit on the zero state. + +```@repl +import Yao +using Yao.EasyBuild: qft_circuit +n = 10; +circuit = qft_circuit(n); # build a quantum Fourier transform circuit +network = Yao.yao2einsum(circuit) # convert this circuit to tensor network +reshape(Yao.contract(network), 1<0 for i=1:n]), final_state=Dict([i=>0 for i=1:n]), + optimizer=Yao.YaoToEinsum.TreeSA(; nslices=3)) # slicing technique +Yao.contract(network)[] ≈ Yao.zero_state(n)' * (Yao.zero_state(n) |> circuit) +``` + +## API +```@docs +yao2einsum +TensorNetwork +optimize_code +contraction_complexity +contract +``` + +## References +[^Pan2022]: Pan, Feng, and Pan Zhang. "Simulation of quantum circuits using the big-batch tensor network method." Physical Review Letters 128.3 (2022): 030501. +[^Kalachev2021]: Kalachev, Gleb, Pavel Panteleev, and Man-Hong Yung. "Recursive multi-tensor contraction for xeb verification of quantum circuits." arXiv preprint arXiv:2108.05665 (2021). +[^Markov2008]: Markov, Igor L., and Yaoyun Shi. "Simulating quantum computation by contracting tensor networks." SIAM Journal on Computing 38.3 (2008): 963-981. +[^Liu2023]: Liu, Jin-Guo, et al. "Computing solution space properties of combinatorial optimization problems via generic tensor networks." SIAM Journal on Scientific Computing 45.3 (2023): A1239-A1270. \ No newline at end of file diff --git a/docs/src/quick-start.md b/docs/src/quick-start.md index ac132eeb4..5807def65 100644 --- a/docs/src/quick-start.md +++ b/docs/src/quick-start.md @@ -11,11 +11,11 @@ is the [`ArrayReg`](@ref), you can create it by feeding a state vector to it, e. ```@repl quick-start using Yao -ArrayReg(rand(ComplexF64, 2^3)) -zero_state(5) -rand_state(5) -product_state(bit"10100") -ghz_state(5) +ArrayReg(randn(ComplexF64, 2^3)) # a random unnormalized 3-qubit state +zero_state(5) # |00000⟩ +rand_state(5) # a random state +product_state(bit"10100") # |10100⟩ +ghz_state(5) # (|00000⟩ + |11111⟩)/√2 ``` the internal quantum state can be accessed via [`statevec`](@ref) method @@ -24,12 +24,11 @@ the internal quantum state can be accessed via [`statevec`](@ref) method statevec(ghz_state(2)) ``` -for more functionalities about registers please refer to the manual of [`registers`](@ref). +for more functionalities about registers please refer to the manual of [Registers](@ref registers). -## Create quantum circuit with Yao blocks +## Create quantum circuit -Yao uses the quantum "block"s to describe quantum circuits, e.g -the following code creates a 2-qubit circuit +Yao introduces an abstract representation for linear maps, called "block"s, which can be used to represent quantum circuits, Hamiltonians, and other quantum operations. The following code creates a 2-qubit circuit ```@repl quick-start chain(2, put(1=>H), put(2=>X)) @@ -39,56 +38,61 @@ where `H` gate is at 1st qubit, `X` gate is at 2nd qubit. A more advanced example is the quantum Fourier transform circuit ```@repl quick-start -A(i, j) = control(i, j=>shift(2π/(1<<(i-j+1)))) +A(i, j) = control(i, j=>shift(2π/(1<<(i-j+1)))) # a cphase gate B(n, k) = chain(n, j==k ? put(k=>H) : A(j, k) for j in k:n) qft(n) = chain(B(n, k) for k in 1:n) -qft(3) +circuit = qft(3) # a 3-qubit QFT circuit +mat(circuit) # the matrix representation of the circuit +apply!(zero_state(3), circuit) # apply the circuit to a zero state ``` -## Create Hamiltonian with Yao blocks +More details about available blocks can be found in the manual of [Blocks](@ref blocks). -the quantum "block"s are expressions on quantum operators, thus, it can -also be used to represent a Hamiltonian, e.g we can create a simple Ising -Hamiltonian on 1D chain as following +## Create Hamiltonian + +We can create a simple Ising Hamiltonian on 1D chain as following ```@repl quick-start -sum(kron(5, i=>Z, mod1(i+1, 5)=>Z) for i in 1:5) +h = sum([kron(5, i=>Z, mod1(i+1, 5)=>Z) for i in 1:5]) # a 5-qubit Ising Hamiltonian +mat(h) # the matrix representation of the Hamiltonian ``` -## Automatic differentiate a Yao block +## Differentiating a quantum circuit Yao has its own automatic differentiation rule implemented, this allows one obtain -gradients of a loss function by simply putting a `'` mark behind [`expect`](@ref) +gradients of a loss function by simply putting a `'` mark following [`expect`](@ref) or [`fidelity`](@ref), e.g +To obtain the gradient of the quantum Fourier transform circuit with respect to its parameters, one can use the following code ```@repl quick-start -expect'(X, zero_state(1)=>Rx(0.2)) +grad_state, grad_circuit_params = expect'(kron(X, X, I2) + kron(I2, X, X), zero_state(3)=>qft(3)) ``` - -or for fiedlity - +where `kron(X, X, I2) + kron(I2, X, X)` is the target Hamiltonian, `zero_state(3)` is the initial state, `qft(3)` is the quantum Fourier transform circuit. +The return value is a vector, each corresponding to the gradient of the loss function with respect to a parameter in the circuit. +The list of parameters can be obtained by [`parameters`](@ref) function. ```@repl quick-start -fidelity'(zero_state(1)=>Rx(0.1), zero_state(1)=>Rx(0.2)) +parameters(qft(3)) ``` -## Combine Yao with ChainRules/Zygote +To obtain the gradient of the fidelity between a state parameterized by a quantum circuit and a target state, one can use the following code +```@repl quick-start +((grad_state1, grad_circuit1), grad_state2) = fidelity'(zero_state(3)=>qft(3), ghz_state(3)) +``` +where `zero_state(3)` is the initial state, `qft(3)` is the quantum Fourier transform circuit, `ghz_state(3)` is the target state. -## Symbolic calculation with Yao block -Yao supports symbolic calculation of quantum circuit via `SymEngine`. We can show - +The automatic differentiation functionality can also be accessed by interfacing with the machine learning libraries [`Zygote`](https://github.com/FluxML/Zygote.jl). ## Plot quantum circuits -The [YaoPlots]() in Yao's ecosystem provides plotting for quantum circuits and ZX diagrams. +The component package `YaoPlots` provides plotting for quantum circuits and ZX diagrams. You can use it to visualize your quantum circuits in [`VSCode`](https://code.visualstudio.com/), [`Jupyter`](https://jupyter.org/) notebook or [`Pluto`](https://github.com/fonsp/Pluto.jl) notebook. ```@example quick-start -using Yao.EasyBuild, YaoPlots +using Yao.EasyBuild, Yao.YaoPlots using Compose # show a qft circuit -Compose.SVG(plot(qft_circuit(5))) +vizcircuit(qft_circuit(5)) ``` -## Convert quantum circuits to tensor network -## Simplify quantum circuit with ZX calculus +More details about the plotting can be found in the manual: [Quantum Circuit Visualization](@ref). \ No newline at end of file diff --git a/ext/CuYao.jl b/ext/CuYao.jl new file mode 100644 index 000000000..1421abcec --- /dev/null +++ b/ext/CuYao.jl @@ -0,0 +1,3 @@ +module CuYao +include("CuYao/src/CuYao.jl") +end diff --git a/ext/CuYao/.gitignore b/ext/CuYao/.gitignore new file mode 100644 index 000000000..56b59ee8d --- /dev/null +++ b/ext/CuYao/.gitignore @@ -0,0 +1,23 @@ +*.jl.cov +*.jl.*.cov +*.jl.mem + +docs/build/ +docs/site/ +docs/src/tutorial/RegisterBasics.md +docs/src/tutorial/BlockBasics.md +docs/src/tutorial/BinaryBasics.md +docs/src/tutorial/QCBM.md + +*.ipynb_checkpoints +**/*.ipynb_checkpoints +**/**/*.ipynb_checkpoints + +_*.dat +*.swp +__pycache__/ +.vscode/ + +Manifest.toml + +_local/ diff --git a/ext/CuYao/README.md b/ext/CuYao/README.md new file mode 100644 index 000000000..becb96326 --- /dev/null +++ b/ext/CuYao/README.md @@ -0,0 +1,8 @@ +## How to run tests + +Open a terminal, switch to Yao folder and run the following commands: +```bash +make init-CuYao +make test-CuYao +``` + diff --git a/ext/CuYao/benchmarks/TestCuda.jl b/ext/CuYao/benchmarks/TestCuda.jl new file mode 100644 index 000000000..8c6a6cd52 --- /dev/null +++ b/ext/CuYao/benchmarks/TestCuda.jl @@ -0,0 +1,54 @@ +using CUDA +using LinearAlgebra +using BenchmarkTools + +function ms!(X::CuArray, s::Number) + function kernel(X, s) + i = (blockIdx().x-1) * blockDim().x + threadIdx().x + @inbounds X[i] *= s + return + end + @cuda blocks=length(X) kernel(X, s) + X +end + +function ms1!(X::CuArray, a, b, c) + k = f(a) + function kernel(X, a, b, c) + i = (blockIdx().x-1) * blockDim().x + threadIdx().x + k(X, i) + k(X, i) + k(X, i) + #@inbounds X[i] *= a + #@inbounds X[i] *= b + #@inbounds X[i] *= c + return + end + @cuda blocks=length(X)÷256 threads=256 kernel(X, a, b, c) + X +end + +@inline function f(s) + @inline function kernel(X, i) + X[i]*=s + end +end +function ms2!(X::CuArray, s::Number) + function kernel(X, s) + i = (blockIdx().x-1) * blockDim().x + threadIdx().x + @inbounds X[i] *= s + return + end + @cuda blocks=length(X)÷256 threads=256 kernel(X, s) + X +end + +a = randn(1<<10) +cua = cu(a) +@benchmark ms1!($cua, 1.001, 1.001, 1.001) +@benchmark ms2!(ms2!(ms2!($cua, 1.001), 1.001), 1.001) +bss = 1:length(cua) +@benchmark f(1.001).(cua, bss) + +@benchmark rmul!($a, 1.01) +@benchmark ms!($cua, 1.0) diff --git a/ext/CuYao/benchmarks/gates.jl b/ext/CuYao/benchmarks/gates.jl new file mode 100644 index 000000000..b63dbd497 --- /dev/null +++ b/ext/CuYao/benchmarks/gates.jl @@ -0,0 +1,19 @@ +using Yao, CuYao, CUDA +using BenchmarkTools + +reg = rand_state(12; nbatch=1000) +creg = reg |> cu +@benchmark CUDA.@sync creg |> put(12, 3=>Z) +@benchmark CUDA.@sync creg |> put(12, 3=>X) +@benchmark reg |> put(12, 3=>Z) +@benchmark CUDA.@sync creg |> control(12, 6, 3=>X) +@benchmark reg |> control(12, 6, 3=>X) +@benchmark CUDA.@sync creg |> put(12, 3=>rot(X, 0.3)) +@benchmark reg |> put(12, 3=>rot(X, 0.3)) + +reg = rand_state(20) +creg = reg |> cu +g = swap(20, 7, 2) +g = put(20, (7, 2)=>matblock(rand_unitary(4))) +@benchmark reg |> g +@benchmark CUDA.@sync creg |> g diff --git a/ext/CuYao/benchmarks/gcompile.jl b/ext/CuYao/benchmarks/gcompile.jl new file mode 100644 index 000000000..62a90ae95 --- /dev/null +++ b/ext/CuYao/benchmarks/gcompile.jl @@ -0,0 +1,92 @@ +using Yao, Yao.Boost, Yao.Intrinsics, StaticArrays, Yao.Blocks +using CuYao, CUDA +using BenchmarkTools, Profile + +nbit = 12 +c = chain(put(nbit, 2=>X), put(nbit, 2=>rot(X, 0.2)), control(nbit, 3, 2=>rot(X,0.3))) +#c = chain(c..., c...,c...) +cc = c |> KernelCompiled +reg = rand_state(nbit) |> cu + +@benchmark $reg |> copy |> $c[1] seconds = 2 +@benchmark $reg |> copy |> $cc seconds = 2 + + +reg = rand_state(9, 1000) +creg = reg |> cu +@benchmark focus!(reg |> copy, 1) do r + measure!(r) + r +end + +@benchmark focus!(creg|>copy, 1) do r + measure!(r) + r +end + +@benchmark focus!(creg, 1) do r + measure_reset!(r) + r +end + +@profile for i=1:10 focus!(creg |> copy, 1) do r + measure!(r) + r +end +end + +a = randn(1<<10) +ca = a |> cu + +gpu_call(ca, (x->x^2, ca)) do state, f, ca + ilin = linear_index(state) + ca[ilin] = f(ca[ilin]) + return +end + +gpu_call(1:1<<10, (x->x^2, ca)) do state, f, ca + ilin = linear_index(state) + ca[ilin] = f(ca[ilin]) + return +end + + +function xx_kernel(bits::Ints) + ctrl = controller(bits[1], 0) + mask = bmask(bits...) + function kernel(state, inds) + i = inds[1] + b = i-1 + ctrl(b) && swaprows!(piecewise(state, inds), i, flip(b, mask) + 1) + end +end + +using CuYao: x_kernel, y_kernel, z_kernel +fx = x_kernel((3,4)) +fy = y_kernel(3) +fz = z_kernel((3,4)) + +a = randn(1<<10) +ca = a |> cu +fs = (fx, fz) + +@benchmark gpu_call(ca, (fs, ca)) do state, f, ca + #ilin = linear_index(state) + ilin = @cartesianidx ca + for fi in fs + fi(ca, ilin) + end + return +end + +@device_code_warntype gpu_call(ca, ((fx, fz), ca)) do state, fs, ca + ilin = @cartesianidx ca + #for fi in fs fi(ca, ilin) end + @Base.Cartesian.@nexprs $(length(fs)) i->fs[i](ca, ilin) + #fs(ca, ilin) + return +end + +a ≈ Vector(ca) + +@edit controller(3,0) diff --git a/ext/CuYao/benchmarks/paralleldot.jl b/ext/CuYao/benchmarks/paralleldot.jl new file mode 100644 index 000000000..38e5dd6e3 --- /dev/null +++ b/ext/CuYao/benchmarks/paralleldot.jl @@ -0,0 +1,11 @@ +function paralleldot(matrices::CuVector, ptrA, ptrB) + @inline function kernel(ctx, matrices) + inds = @cartesianidx state + i = inds[1] + piecewise(state, inds)[i] *= anyone(i-1, mask) ? d : a + return + end + gpu_call(kernel, state, a, d, mask; elements=length(state)) + return state +end + diff --git a/ext/CuYao/benchmarks/statexpect.jl b/ext/CuYao/benchmarks/statexpect.jl new file mode 100644 index 000000000..9e3ca5ddc --- /dev/null +++ b/ext/CuYao/benchmarks/statexpect.jl @@ -0,0 +1,13 @@ +using CuYao +using Yao +using BenchmarkTools + +sf(x, y) = abs(x-y) +a = randn(1024) +ca = a |> cu +b = randn(1024) +cb = b |> cu +@benchmark expect(StatFunctional{2}(sf), a, b) seconds=2 +@benchmark expect(StatFunctional{2}(sf), a) seconds=2 +@benchmark expect(StatFunctional{2}(sf), ca, cb) seconds=2 +@benchmark expect(StatFunctional{2}(sf), ca) seconds=2 diff --git a/ext/CuYao/src/CUDApatch.jl b/ext/CuYao/src/CUDApatch.jl new file mode 100644 index 000000000..aeaf7d367 --- /dev/null +++ b/ext/CuYao/src/CUDApatch.jl @@ -0,0 +1,64 @@ +# TODO +# support norm(view(reshape(A, m, n), :, 1)) +norm2(A::DenseCuArray; dims=1) = mapreduce(abs2, +, A, dims=dims) .|> CUDA.sqrt + +piecewise(state::AbstractVector, inds) = state +piecewise(state::AbstractMatrix, inds) = @inbounds view(state,:,inds[2]) + +function kron(A::DenseCuArray{T1}, B::DenseCuArray{T2}) where {T1, T2} + res = CUDA.zeros(promote_type(T1,T2), (size(A).*size(B))...) + @inline function kernel(ctx, res, A, B) + state = @linearidx res + @inbounds inds = CartesianIndices(res)[state].I + inds_A = (inds.-1) .÷ size(B) .+ 1 + inds_B = (inds.-1) .% size(B) .+ 1 + @inbounds res[state] = A[inds_A...]*B[inds_B...] + return + end + + gpu_call(kernel, res, A, B) + res +end + +""" + kron!(C::CuArray, A, B) + +Computes Kronecker products in-place on the GPU. +The results are stored in 'C', overwriting the existing values of 'C'. +""" +function Yao.YaoArrayRegister.kron!(C::CuArray{T3}, A::DenseCuArray{T1}, B::DenseCuArray{T2}) where {T1, T2, T3} + @boundscheck (size(C) == (size(A,1)*size(B,1), size(A,2)*size(B,2))) || throw(DimensionMismatch()) + CI = Base.CartesianIndices(C) + @inline function kernel(ctx, C, A, B) + state = @linearidx C + @inbounds inds = CI[state].I + inds_A = (inds.-1) .÷ size(B) .+ 1 + inds_B = (inds.-1) .% size(B) .+ 1 + @inbounds C[state] = A[inds_A...]*B[inds_B...] + return + end + + gpu_call(kernel, C, A, B) + C +end + +function getindex(A::DenseCuVector{T}, B::DenseCuArray{<:Integer}) where T + res = CUDA.zeros(T, size(B)...) + @inline function kernel(ctx, res, A, B) + state = @linearidx res + @inbounds res[state] = A[B[state]] + return + end + gpu_call(kernel, res, A, B) + res +end + +function getindex(A::AbstractVector, B::DenseCuArray{<:Integer}) + A[Array(B)] +end + +YaoBlocks.AD.as_scalar(x::DenseCuArray) = Array(x)[] + +# patch for ExponentialUtilities +YaoBlocks.compatible_multiplicative_operand(::CuArray, source::AbstractArray) = CuArray(source) +YaoBlocks.compatible_multiplicative_operand(::CuArray, source::CuArray) = source diff --git a/ext/CuYao/src/CuYao.jl b/ext/CuYao/src/CuYao.jl new file mode 100644 index 000000000..2a442181a --- /dev/null +++ b/ext/CuYao/src/CuYao.jl @@ -0,0 +1,39 @@ +# module CuYao +using Yao.YaoArrayRegister.LuxurySparse, Yao.YaoArrayRegister.StaticArrays, LinearAlgebra, Base.Cartesian +using Yao.YaoArrayRegister.StatsBase, Yao.YaoArrayRegister.BitBasis +using Yao.Reexport +import Yao.YaoArrayRegister.TupleTools +using Yao.YaoArrayRegister.Random + +using Yao.YaoArrayRegister +using Yao +using CUDA +using CUDA.GPUArrays: gpu_call, @linearidx, @cartesianidx, linear_index +using Yao.YaoArrayRegister +using Yao.YaoBlocks +using Yao.ConstGate: SWAPGate +using Yao.ConstGate: S, T, Sdag, Tdag + +import Yao.YaoArrayRegister: insert_qudits!, join +import CUDA: cu +import Yao.YaoArrayRegister: _measure, measure, measure! +import Yao.YaoArrayRegister: batch_normalize! +import Yao.YaoBlocks: BlockedBasis, nblocks, subblock +import Yao: expect +import Yao.YaoArrayRegister: u1rows!, unrows!, autostatic, instruct!, swaprows!, mulrow! +import Yao.YaoArrayRegister.LinearAlgebra: norm +import Base: kron, getindex + +import Yao: cpu, cuzero_state, cuuniform_state, curand_state, cuproduct_state, cughz_state + +const Ints = NTuple{<:Any, Int} + +include("CUDApatch.jl") +include("register.jl") +include("instructs.jl") +include("yao2einsum.jl") + +function __init__() + CUDA.allowscalar(false) +end +# end diff --git a/ext/CuYao/src/instructs.jl b/ext/CuYao/src/instructs.jl new file mode 100644 index 000000000..6ed8e35df --- /dev/null +++ b/ext/CuYao/src/instructs.jl @@ -0,0 +1,312 @@ +# get index +macro idx(shape, grididx=1, ctxsym=:ctx) + quote + x = $(esc(shape)) + i = $linear_index($(esc(ctxsym)), $(esc(grididx))) + i > $prod2(x) && return + @inbounds Base.CartesianIndices(x)[i].I + end +end +replace_first(x::NTuple{2}, v) = (v, x[2]) +replace_first(x::NTuple{1}, v) = (v,) +prod2(x::NTuple{2}) = x[1] * x[2] +prod2(x::NTuple{1}) = x[1] + +gpu_compatible(A::AbstractVecOrMat) = A |> staticize +gpu_compatible(A::StaticArray) = A + +###################### unapply! ############################ +function instruct!(::Val{2}, state::DenseCuVecOrMat, U0::AbstractMatrix, locs::NTuple{M, Int}, clocs::NTuple{C, Int}, cvals::NTuple{C, Int}) where {C, M} + @debug "The generic U(N) matrix of size ($(size(U0))), on: GPU, locations: $(locs), controlled by: $(clocs) = $(cvals)." + nbit = log2dim1(state) + # reorder a unirary matrix. + U = gpu_compatible(all(TupleTools.diff(locs).>0) ? U0 : reorder(U0, collect(locs)|>sortperm)) + locked_bits = [clocs..., locs...] + locked_vals = [cvals..., zeros(Int, M)...] + locs_raw = gpu_compatible([i+1 for i in itercontrol(nbit, setdiff(1:nbit, locs), zeros(Int, nbit-M))]) + configs = itercontrol(nbit, locked_bits, locked_vals) + + len = length(configs) + @inline function kernel(ctx, state, locs_raw, U, configs, len) + CUDA.assume(len > 0) + sz = size(state) + CUDA.assume(length(sz) == 1 || length(sz) == 2) + inds = @idx replace_first(sz, len) + x = @inbounds configs[inds[1]] + @inbounds unrows!(piecewise(state, inds), x .+ locs_raw, U) + return + end + + @inline function kernel_single_entry_diag(ctx, state, loc, val, configs, len) + CUDA.assume(len > 0) + sz = size(state) + CUDA.assume(length(sz) == 1 || length(sz) == 2) + inds = @idx replace_first(sz, len) + x = @inbounds configs[inds[1]] + @inbounds piecewise(state, inds)[x + loc] *= val + return + end + + elements = len*size(state,2) + if U isa Diagonal && count(!isone, U.diag) == 1 + @debug "The single entry diagonal matrix, on: GPU, locations: $(locs), controlled by: $(clocs) = $(cvals)." + k = findfirst(!isone, U.diag) + loc = locs_raw[k] + val = U.diag[k] + gpu_call(kernel_single_entry_diag, state, loc, val, configs, len; elements) + else + gpu_call(kernel, state, locs_raw, U, configs, len; elements) + end + state +end +instruct!(::Val{2}, state::DenseCuVecOrMat, U0::IMatrix, locs::NTuple{M, Int}, clocs::NTuple{C, Int}, cvals::NTuple{C, Int}) where {C, M} = state +instruct!(::Val{2}, state::DenseCuVecOrMat, U0::SDSparseMatrixCSC, locs::NTuple{M, Int}, clocs::NTuple{C, Int}, cvals::NTuple{C, Int}) where {C, M} = instruct!(Val(2), state, U0 |> Matrix, locs, clocs, cvals) + +################## General U1 apply! ################### +function YaoArrayRegister.single_qubit_instruct!(state::DenseCuVecOrMat, U1::SDSparseMatrixCSC, loc::Int) + instruct!(Val(2), state, Matrix(U1), loc, clocs, cval) +end +function YaoArrayRegister.single_qubit_instruct!(state::DenseCuVecOrMat, U1::AbstractMatrix, loc::Int) + @debug "The generic U(2) matrix of size ($(size(U1))), on: GPU, locations: $(loc)." + a, c, b, d = U1 + nbit = log2dim1(state) + step = 1<<(loc-1) + configs = itercontrol(nbit, [loc], [0]) + + len = length(configs) + @inline function kernel(ctx, state, a, b, c, d, len) + inds = @idx replace_first(size(state), len) + i = @inbounds configs[inds[1]]+1 + @inbounds u1rows!(piecewise(state, inds), i, i+step, a, b, c, d) + return + end + gpu_call(kernel, state, a, b, c, d, len; elements=len*size(state,2)) + return state +end + +function YaoArrayRegister.single_qubit_instruct!(state::DenseCuVecOrMat, U1::SDPermMatrix, loc::Int) + @debug "The single qubit permutation matrix of size ($(size(U1))), on: GPU, locations: $(loc)." + nbit = log2dim1(state) + b, c = U1.vals + step = 1<<(loc-1) + configs = itercontrol(nbit, [loc], [0]) + + len = length(configs) + function kernel(ctx, state, b, c, step, len, configs) + inds = @idx replace_first(size(state), len) + x = @inbounds configs[inds[1]] + 1 + @inbounds swaprows!(piecewise(state, inds), x, x+step, c, b) + return + end + gpu_call(kernel, state, b, c, step, len, configs; elements=len*size(state,2)) + return state +end + +function YaoArrayRegister.single_qubit_instruct!(state::DenseCuVecOrMat, U1::SDDiagonal, loc::Int) + @debug "The single qubit diagonal matrix of size ($(size(U1))), on: GPU, locations: $(loc)." + a, d = U1.diag + nbit = log2dim1(state) + mask = bmask(loc) + @inline function kernel(ctx, state, a, d, mask) + inds = @cartesianidx state + i = inds[1] + piecewise(state, inds)[i] *= anyone(i-1, mask) ? d : a + return + end + gpu_call(kernel, state, a, d, mask; elements=length(state)) + return state +end + +YaoArrayRegister.single_qubit_instruct!(state::DenseCuVecOrMat, U::IMatrix, loc::Int) = state + +################## XYZ ############# + +_instruct!(state::DenseCuVecOrMat, ::Val{:X}, locs::NTuple{L,Int}) where {L} = _instruct!(state, Val(:X), locs, (), ()) +function _instruct!(state::DenseCuVecOrMat, ::Val{:X}, locs::NTuple{L,Int}, clocs::NTuple{C, Int}, cvals::NTuple{C, Int}) where {L,C} + @debug "The X gate, on: GPU, locations: $(locs), controlled by: $(clocs) = $(cvals)." + length(locs) == 0 && return state + nbit = log2dim1(state) + configs = itercontrol(nbit, [clocs..., locs[1]], [cvals..., 0]) + mask = bmask(locs...) + len = length(configs) + @inline function kernel(ctx, state, mask, len, configs) + inds = @idx replace_first(size(state), len) + b = @inbounds configs[inds[1]] + @inbounds swaprows!(piecewise(state, inds), b+1, flip(b, mask) + 1) + return + end + gpu_call(kernel, state, mask, len, configs; elements=len*size(state,2)) + return state +end + +function _instruct!(state::DenseCuVecOrMat, ::Val{:Y}, locs::NTuple{C,Int}) where C + @debug "The Y gate, on: GPU, locations: $(locs)." + length(locs) == 0 && return state + nbit = log2dim1(state) + mask = bmask(Int, locs...) + configs = itercontrol(nbit, [locs[1]], [0]) + bit_parity = length(locs)%2 == 0 ? 1 : -1 + factor = (-im)^length(locs) + len = length(configs) + @inline function kernel(ctx, state, mask, bit_parity, configs, len) + inds = @idx replace_first(size(state), len) + b = @inbounds configs[inds[1]] + i_ = flip(b, mask) + 1 + factor1 = count_ones(b&mask)%2 == 1 ? -factor : factor + factor2 = factor1*bit_parity + @inbounds swaprows!(piecewise(state, inds), b+1, i_, factor2, factor1) + return + end + gpu_call(kernel, state, mask, bit_parity, configs, len; elements=len*size(state,2)) + return state +end + +function _instruct!(state::DenseCuVecOrMat, ::Val{:Y}, locs::Tuple{Int}, clocs::NTuple{C, Int}, cvals::NTuple{C, Int}) where C + @debug "The Y gate, on: GPU, locations: $(locs), controlled by: $(clocs) = $(cvals)." + length(locs) == 0 && return state + nbit = log2dim1(state) + configs = itercontrol(nbit, [clocs..., locs...], [cvals..., 0]) + mask = bmask(locs...) + len = length(configs) + @inline function kernel(ctx, state, configs, mask, len) + inds = @idx replace_first(size(state), len) + b = @inbounds configs[inds[1]] + @inbounds swaprows!(piecewise(state, inds), b+1, flip(b, mask) + 1, im, -im) + return + end + gpu_call(kernel, state, configs, mask, len; elements=len*size(state,2)) + return state +end + +function _instruct!(state::DenseCuVecOrMat, ::Val{:Z}, locs::NTuple{C,Int}) where C + @debug "The Z gate, on: GPU, locations: $(locs)." + length(locs) == 0 && return state + nbit = log2dim1(state) + mask = bmask(locs...) + @inline function kernel(ctx, state, mask) + inds = @cartesianidx state + i = inds[1] + piecewise(state, inds)[i] *= count_ones((i-1)&mask)%2==1 ? -1 : 1 + return + end + gpu_call(kernel, state, mask; elements=length(state)) + return state +end + + +for (G, FACTOR) in zip([:Z, :S, :T, :Sdag, :Tdag], [:(-one(Int32)), :(1f0im), :($(exp(im*π/4))), :(-1f0im), :($(exp(-im*π/4)))]) + if G !== :Z + @eval function _instruct!(state::DenseCuVecOrMat, ::Val{$(QuoteNode(G))}, locs::NTuple{C,Int}) where C + @debug "The $($(QuoteNode(G))) gate, on: GPU, locations: $(locs)." + length(locs) == 0 && return state + nbit = log2dim1(state) + mask = bmask(Int32, locs...) + @inline function kernel(ctx, state, mask) + inds = @cartesianidx state + i = inds[1] + piecewise(state, inds)[i] *= $FACTOR ^ count_ones(Int32(i-1)&mask) + return + end + gpu_call(kernel, state, mask; elements=length(state)) + return state + end + end + @eval function _instruct!(state::DenseCuVecOrMat, ::Val{$(QuoteNode(G))}, locs::Tuple{Int}, clocs::NTuple{C, Int}, cvals::NTuple{C, Int}) where C + @debug "The $($(QuoteNode(G))) gate, on: GPU, locations: $(locs), controlled by: $(clocs) = $(cvals)." + ctrl = controller((clocs..., locs...), (cvals..., 1)) + @inline function kernel(ctx, state, ctrl) + inds = @cartesianidx state + i = inds[1] + piecewise(state, inds)[i] *= ctrl(i-1) ? $FACTOR : one($FACTOR) + return + end + gpu_call(kernel, state, ctrl; elements=length(state)) + return state + end +end + +for G in [:X, :Y, :Z, :S, :T, :Sdag, :Tdag] + @eval begin + function YaoArrayRegister.instruct!(::Val{2}, state::DenseCuVecOrMat, g::Val{$(QuoteNode(G))}, locs::NTuple{C,Int}) where C + _instruct!(state, g, locs) + end + + function YaoArrayRegister.instruct!(::Val{2}, state::DenseCuVecOrMat, g::Val{$(QuoteNode(G))}, locs::Tuple{Int}) + _instruct!(state, g, locs) + end + + function YaoArrayRegister.instruct!(::Val{2}, state::DenseCuVecOrMat, g::Val{$(QuoteNode(G))}, locs::Tuple{Int}, clocs::NTuple{C, Int}, cvals::NTuple{C, Int}) where C + _instruct!(state, g, locs, clocs, cvals) + end + + function YaoArrayRegister.instruct!(::Val{2}, state::DenseCuVecOrMat, vg::Val{$(QuoteNode(G))}, locs::Tuple{Int}, cloc::Tuple{Int}, cval::Tuple{Int}) + _instruct!(state, vg, locs, cloc, cval) + end + end + +end + +function instruct!(::Val{2}, state::DenseCuVecOrMat, ::Val{:SWAP}, locs::Tuple{Int,Int}) + @debug "The SWAP gate, on: GPU, locations: $(locs)." + b1, b2 = locs + mask1 = bmask(b1) + mask2 = bmask(b2) + + configs = itercontrol(log2dim1(state), [locs...], [1,0]) + function kf(ctx, state, mask1, mask2) + inds = @idx replace_first(size(state), length(configs)) + + b = configs[inds[1]] + i = b+1 + i_ = b ⊻ (mask1|mask2) + 1 + swaprows!(piecewise(state, inds), i, i_) + nothing + end + gpu_call(kf, state, mask1, mask2; elements=length(configs)*size(state,2)) + state +end + +############## other gates ################ +# parametrized swap gate + +function instruct!(::Val{2}, state::DenseCuVecOrMat, ::Val{:PSWAP}, locs::Tuple{Int, Int}, θ::Real) + @debug "The PSWAP gate, on: GPU, locations: $(locs)." + nbit = log2dim1(state) + mask1 = bmask(locs[1]) + mask2 = bmask(locs[2]) + mask12 = mask1|mask2 + a, c, b_, d = mat(Rx(θ)) + e = exp(-im/2*θ) + configs = itercontrol(nbit, [locs...], [0,0]) + len = length(configs) + @inline function kernel(ctx, state, mask2, mask12, configs, a, b_, c, d) + inds = @idx replace_first(size(state), len) + @inbounds x = configs[inds[1]] + piecewise(state, inds)[x+1] *= e + piecewise(state, inds)[x⊻mask12+1] *= e + y = x ⊻ mask2 + @inbounds u1rows!(piecewise(state, inds), y+1, y⊻mask12+1, a, b_, c, d) + return + end + gpu_call(kernel, state, mask2, mask12, configs, a, b_, c, d; elements=len*size(state,2)) + state +end + +function YaoBlocks._apply_fallback!(r::AbstractCuArrayReg{B,T}, b::AbstractBlock) where {B,T} + YaoBlocks._check_size(r, b) + r.state .= CUDA.adapt(CuArray{T}, mat(T, b)) * r.state + return r +end + +for RG in [:Rx, :Ry, :Rz] + @eval function instruct!( + ::Val{2}, + state::DenseCuVecOrMat{T}, + ::Val{$(QuoteNode(RG))}, + locs::Tuple{Int}, + theta::Number + ) where {T} + YaoArrayRegister.instruct!(Val(2), state, Val($(QuoteNode(RG))), locs, (), (), theta) + return state + end +end diff --git a/ext/CuYao/src/register.jl b/ext/CuYao/src/register.jl new file mode 100644 index 000000000..01c83cbc6 --- /dev/null +++ b/ext/CuYao/src/register.jl @@ -0,0 +1,270 @@ +cu(reg::ArrayReg{D}) where D = ArrayReg{D}(CuArray(reg.state)) +cpu(reg::ArrayReg{D}) where D = ArrayReg{D}(Array(reg.state)) +cu(reg::BatchedArrayReg{D}) where D = BatchedArrayReg{D}(CuArray(reg.state), reg.nbatch) +cpu(reg::BatchedArrayReg{D}) where D = BatchedArrayReg{D}(Array(reg.state), reg.nbatch) +cu(reg::DensityMatrix{D}) where D = DensityMatrix{D}(CuArray(reg.state)) +cpu(reg::DensityMatrix{D}) where D = DensityMatrix{D}(Array(reg.state)) +const AbstractCuArrayReg{D, T, MT} = AbstractArrayReg{D, T, MT} where MT<:DenseCuArray +const CuArrayReg{D, T, MT} = ArrayReg{D, T, MT} where MT<:DenseCuArray +const CuBatchedArrayReg{D, T, MT} = BatchedArrayReg{D, T, MT} where MT<:DenseCuArray +const CuDensityMatrix{D, T, MT} = DensityMatrix{D, T, MT} where MT<:DenseCuMatrix + +function batch_normalize!(s::DenseCuArray, p::Real=2) + p!=2 && throw(ArgumentError("p must be 2!")) + s./=norm2(s, dims=1) + return s +end + +@inline function tri2ij(l::Int) + i = ceil(Int, sqrt(2*l+0.25)-0.5) + j = l-i*(i-1)÷2 + return i+1,j +end + +############### MEASURE ################## +function measure(::ComputationalBasis, reg::ArrayReg{D, T, MT} where MT<:DenseCuArray, ::AllLocs; rng::AbstractRNG=Random.GLOBAL_RNG, nshots::Int=1) where {D,T} + _measure(rng, basis(reg), reg |> probs |> Vector, nshots) +end + +# TODO: optimize the batch dimension using parallel sampling +function measure(::ComputationalBasis, reg::BatchedArrayReg{D, T, MT} where MT<:DenseCuArray, ::AllLocs; rng::AbstractRNG=Random.GLOBAL_RNG, nshots::Int=1) where {D,T} + regm = reg |> rank3 + pl = dropdims(mapreduce(abs2, +, regm, dims=2), dims=2) + return _measure(rng, basis(reg), pl |> Matrix, nshots) +end + +function measure!(::RemoveMeasured, ::ComputationalBasis, reg::AbstractCuArrayReg{D}, ::AllLocs; rng::AbstractRNG=Random.GLOBAL_RNG) where D + regm = reg |> rank3 + B = size(regm, 3) + nregm = similar(regm, D ^ nremain(reg), B) + pl = dropdims(mapreduce(abs2, +, regm, dims=2), dims=2) + pl_cpu = pl |> Matrix + res_cpu = map(ib->_measure(rng, basis(reg), view(pl_cpu, :, ib), 1)[], 1:B) + res = CuArray(res_cpu) + CI = Base.CartesianIndices(nregm) + @inline function kernel(ctx, nregm, regm, res, pl) + state = @linearidx nregm + @inbounds i,j = CI[state].I + @inbounds r = Int(res[j])+1 + @inbounds nregm[i,j] = regm[r,i,j]/CUDA.sqrt(pl[r, j]) + return + end + gpu_call(kernel, nregm, regm, res, pl) + reg.state = reshape(nregm,1,:) + return reg isa ArrayReg ? Array(res)[] : res +end + +function measure!(::NoPostProcess, ::ComputationalBasis, reg::AbstractCuArrayReg{D, T}, ::AllLocs; rng::AbstractRNG=Random.GLOBAL_RNG) where {D, T} + regm = reg |> rank3 + B = size(regm, 3) + pl = dropdims(mapreduce(abs2, +, regm, dims=2), dims=2) + pl_cpu = pl |> Matrix + res_cpu = map(ib->_measure(rng, basis(reg), view(pl_cpu, :, ib), 1)[], 1:B) + res = CuArray(res_cpu) + CI = Base.CartesianIndices(regm) + + @inline function kernel(ctx, regm, res, pl) + state = @linearidx regm + @inbounds k,i,j = CI[state].I + @inbounds rind = Int(res[j]) + 1 + @inbounds regm[k,i,j] = k==rind ? regm[k,i,j]/CUDA.sqrt(pl[k, j]) : T(0) + return + end + gpu_call(kernel, regm, res, pl) + return reg isa ArrayReg ? Array(res)[] : res +end + +function YaoArrayRegister.measure!( + ::NoPostProcess, + bb::BlockedBasis, + reg::AbstractCuArrayReg{D,T}, + ::AllLocs; + rng::AbstractRNG = Random.GLOBAL_RNG, +) where {D,T} + state = @inbounds (reg|>rank3)[bb.perm, :, :] # permute to make eigen values sorted + B = size(state, 3) + pl = dropdims(mapreduce(abs2, +, state, dims=2), dims=2) + pl_cpu = pl |> Matrix + pl_block = zeros(eltype(pl), nblocks(bb), B) + @inbounds for ib = 1:B + for i = 1:nblocks(bb) + for k in subblock(bb, i) + pl_block[i, ib] += pl_cpu[k, ib] + end + end + end + # perform measurements on CPU + res_cpu = Vector{Int}(undef, B) + @inbounds @views for ib = 1:B + ires = sample(rng, 1:nblocks(bb), Weights(pl_block[:, ib])) + # notice ires is `BitStr` type, can be use as indices directly. + range = subblock(bb, ires) + state[range, :, ib] ./= sqrt(pl_block[ires, ib]) + state[1:range.start-1, :, ib] .= zero(T) + state[range.stop+1:size(state, 1), :, ib] .= zero(T) + res_cpu[ib] = ires + end + # undo permute and assign back + _state = reshape(state, 1 << nactive(reg), :) + rstate = reshape(reg.state, 1 << nactive(reg), :) + @inbounds rstate[bb.perm, :] .= _state + return reg isa ArrayReg ? bb.values[res_cpu[]] : CuArray(bb.values[res_cpu]) +end + +function measure!(rst::ResetTo, ::ComputationalBasis, reg::AbstractCuArrayReg{D, T}, ::AllLocs; rng::AbstractRNG=Random.GLOBAL_RNG) where {D, T} + regm = reg |> rank3 + B = size(regm, 3) + pl = dropdims(mapreduce(abs2, +, regm, dims=2), dims=2) + pl_cpu = pl |> Matrix + res_cpu = map(ib->_measure(rng, basis(reg), view(pl_cpu, :, ib), 1)[], 1:B) + res = CuArray(res_cpu) + CI = Base.CartesianIndices(regm) + + @inline function kernel(ctx, regm, res, pl, val) + state = @linearidx regm + @inbounds k,i,j = CI[state].I + @inbounds rind = Int(res[j]) + 1 + @inbounds k==val+1 && (regm[k,i,j] = regm[rind,i,j]/CUDA.sqrt(pl[rind, j])) + CUDA.sync_threads() + @inbounds k!=val+1 && (regm[k,i,j] = 0) + return + end + + gpu_call(kernel, regm, res, pl, rst.x) + return reg isa ArrayReg ? Array(res)[] : res +end + +function YaoArrayRegister.batched_kron(A::DenseCuArray{T1}, B::DenseCuArray{T2}) where {T1 ,T2} + res = CUDA.zeros(promote_type(T1,T2), size(A,1)*size(B, 1), size(A,2)*size(B,2), size(A, 3)) + CI = Base.CartesianIndices(res) + @inline function kernel(ctx, res, A, B) + state = @linearidx res + @inbounds i,j,b = CI[state].I + i_A, i_B = divrem((i-1), size(B,1)) + j_A, j_B = divrem((j-1), size(B,2)) + @inbounds res[state] = A[i_A+1, j_A+1, b]*B[i_B+1, j_B+1, b] + return + end + + gpu_call(kernel, res, A, B) + return res +end + +""" + YaoArrayRegister.batched_kron!(C::CuArray, A, B) + +Performs batched Kronecker products in-place on the GPU. +The results are stored in 'C', overwriting the existing values of 'C'. +""" +function YaoArrayRegister.batched_kron!(C::CuArray{T3, 3}, A::DenseCuArray, B::DenseCuArray) where {T3} + @boundscheck (size(C) == (size(A,1)*size(B,1), size(A,2)*size(B,2), size(A,3))) || throw(DimensionMismatch()) + @boundscheck (size(A,3) == size(B,3) == size(C,3)) || throw(DimensionMismatch()) + CI = Base.CartesianIndices(C) + @inline function kernel(ctx, C, A, B) + state = @linearidx C + @inbounds i,j,b = CI[state].I + i_A, i_B = divrem((i-1), size(B,1)) + j_A, j_B = divrem((j-1), size(B,2)) + @inbounds C[state] = A[i_A+1, j_A+1, b]*B[i_B+1, j_B+1, b] + return + end + + gpu_call(kernel, C, A, B) + return C +end + +function join(reg1::AbstractCuArrayReg{D}, reg2::AbstractCuArrayReg{D}) where {D} + @assert nbatch(reg1) == nbatch(reg2) + s1 = reg1 |> rank3 + s2 = reg2 |> rank3 + state = YaoArrayRegister.batched_kron(s1, s2) + return arrayreg(copy(reshape(state, size(state, 1), :)); nlevel=D, nbatch=nbatch(reg1)) +end + +function Yao.insert_qudits!(reg::AbstractCuArrayReg{D}, loc::Int; nqudits::Int=1) where D + na = nactive(reg) + focus!(reg, 1:loc-1) + reg2 = join(zero_state(nqudits; nbatch=nbatch(reg)) |> cu, reg) |> relax! |> focus!((1:na+nqudits)...) + reg.state = reg2.state + return reg +end + +cuproduct_state(bit_str::BitStr; nbatch::Union{NoBatch,Int} = NoBatch()) = + cuproduct_state(ComplexF64, bit_str; nbatch = nbatch) +cuproduct_state(bit_str::AbstractVector; nbatch::Union{NoBatch,Int} = NoBatch()) = + cuproduct_state(ComplexF64, bit_str; nbatch = nbatch) +cuproduct_state(total::Int, bit_config::Integer; kwargs...) = + cuproduct_state(ComplexF64, total, bit_config; kwargs...) +cuproduct_state(::Type{T}, bit_str::BitStr{N}; kwargs...) where {T,N} = + cuproduct_state(T, N, buffer(bit_str); kwargs...) +cuproduct_state(::Type{T}, bit_configs::AbstractVector; kwargs...) where {T} = + cuproduct_state(T, bit_literal(bit_configs...); kwargs...) +function cuproduct_state( + ::Type{T}, + total::Int, + bit_config::Integer; + nbatch::Union{Int,NoBatch} = NoBatch(), + nlevel::Int=2, +) where {T} + raw = CUDA.zeros(T, nlevel ^ total, YaoArrayRegister._asint(nbatch)) + raw[Int(bit_config)+1,:] .= Ref(one(T)) + return arrayreg(raw; nbatch=nbatch, nlevel=nlevel) +end + +cuzero_state(n::Int; kwargs...) = cuzero_state(ComplexF64, n; kwargs...) +cuzero_state(::Type{T}, n::Int; kwargs...) where {T} = cuproduct_state(T, n, 0; kwargs...) + +curand_state(n::Int; kwargs...) = curand_state(ComplexF64, n; kwargs...) + +function curand_state( + ::Type{T}, + n::Int; + nbatch::Union{Int,NoBatch} = NoBatch(), + nlevel = 2, +) where {T} + raw = CUDA.randn(T, nlevel ^ n, YaoArrayRegister._asint(nbatch)) + return normalize!(arrayreg(raw; nbatch=nbatch, nlevel=nlevel)) +end + +cuuniform_state(n::Int; kwargs...) = cuuniform_state(ComplexF64, n; kwargs...) +function cuuniform_state(::Type{T}, n::Int; + nbatch::Union{Int,NoBatch} = NoBatch(), + nlevel::Int = 2, +) where {T} + raw = CUDA.ones(T, nlevel ^ n, YaoArrayRegister._asint(nbatch)) + return normalize!(arrayreg(raw; nbatch=nbatch, nlevel=nlevel)) +end + +cughz_state(n::Int; kwargs...) = cughz_state(ComplexF64, n; kwargs...) +function cughz_state(::Type{T}, n::Int; kwargs...) where {T} + reg = cuzero_state(T, n; kwargs...) + reg.state[1:1,:] .= Ref(sqrt(T(0.5))) + reg.state[end:end,:] .= Ref(sqrt(T(0.5))) + return reg +end + +#= +for FUNC in [:measure!, :measure!] + @eval function $FUNC(rng::AbstractRNG, op::AbstractBlock, reg::AbstractCuArrayReg, al::AllLocs; kwargs...) where B + E, V = eigen!(mat(op) |> Matrix) + ei = Eigen(E|>cu, V|>cu) + $FUNC(rng::AbstractRNG, ei, reg, al; kwargs...) + end +end +=# + +function YaoBlocks.expect(op::AbstractAdd, dm::CuDensityMatrix) + sum(x->expect(x, dm), subblocks(op)) +end +function YaoBlocks.expect(op::AbstractBlock, dm::CuDensityMatrix{D}) where D + return tr(apply(ArrayReg{D}(dm.state), op).state) +end + +measure( + ::ComputationalBasis, + reg::CuDensityMatrix, + ::AllLocs; + nshots::Int = 1, + rng::AbstractRNG = Random.GLOBAL_RNG, +) = YaoArrayRegister._measure(rng, basis(reg), Array(reg |> probs), nshots) + diff --git a/ext/CuYao/src/yao2einsum.jl b/ext/CuYao/src/yao2einsum.jl new file mode 100644 index 000000000..6d8b0d4e3 --- /dev/null +++ b/ext/CuYao/src/yao2einsum.jl @@ -0,0 +1,3 @@ +function CUDA.cu(tnet::TensorNetwork) + return TensorNetwork(tnet.code, tnet.tensors .|> CuArray) +end diff --git a/ext/CuYao/test/CUDApatch.jl b/ext/CuYao/test/CUDApatch.jl new file mode 100644 index 000000000..136d7aeb7 --- /dev/null +++ b/ext/CuYao/test/CUDApatch.jl @@ -0,0 +1,35 @@ +using Yao +using CUDA +using Test +using Yao.YaoBlocks + +@testset "isapprox-complex" begin + ca = CuArray(randn(ComplexF64,3,3)) + cb = copy(ca) + #@test ca ≈ cb still error! + cb[1:1, 1:1] .+= 1e-7im + @test isapprox(ca, cb, atol=1e-5) + @test !isapprox(ca, cb, atol=1e-9) +end + +@testset "view general" begin + a = randn(5,6,8) |> CuArray + @test view(a, 2:4, 4, [1,4,3]) |> size == (3, 3) +end + +@testset "permutedims vector" begin + ca = randn(ComplexF64,3,4,5,1) + @test permutedims(CuArray(ca), [2,1,4,3]) |> Array ≈ permutedims(ca, [2,1,4,3]) +end + +@testset "Complex pow" begin + for T in [ComplexF64, ComplexF32] + a = CuArray(randn(T, 4, 4)) + @test Array(a .^ Int32(3)) ≈ Array(a).^3 + @test Array(a .^ real(T)(3)) ≈ Array(a).^3 + end +end + +@testset "as_scalar" begin + @test YaoBlocks.AD.as_scalar([2.0] |> CuArray) == 2.0 +end diff --git a/ext/CuYao/test/Project.toml b/ext/CuYao/test/Project.toml new file mode 100644 index 000000000..2e5f158c1 --- /dev/null +++ b/ext/CuYao/test/Project.toml @@ -0,0 +1,23 @@ +[deps] +CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +Reexport = "189a3867-3050-52da-a836-e630ba90ab69" +TupleTools = "9d95972d-f1c8-5527-a6e0-b4b365fa01f6" +Yao = "5872b779-8223-5990-8dd0-5abbb0748c8c" + +[compat] +CUDA = "4, 5" +Reexport = "0.2, 1" +TupleTools = "1" +Yao = "0.8" +YaoBlocks = "0.13.10" +julia = "1" + +[extras] +Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +YaoBlocks = "418bc28f-b43b-5e0b-a6e7-61bbc1a2c1df" + +[targets] +test = ["Test", "Statistics", "YaoBlocks"] diff --git a/ext/CuYao/test/extra.jl b/ext/CuYao/test/extra.jl new file mode 100644 index 000000000..95e759ce9 --- /dev/null +++ b/ext/CuYao/test/extra.jl @@ -0,0 +1,52 @@ +using Yao, Test, CUDA +CUDA.allowscalar(false) + +@testset "gradient" begin + reg1 = rand_state(10) + reg2 = rand_state(10) + c = EasyBuild.variational_circuit(10) + g1, g2 = fidelity'(reg1=>c, reg2) + cg1, cg2 = fidelity'(cu(reg1)=>c, cu(reg2)) + @test g1.first ≈ cpu(cg1.first) + @test g1.second ≈ cg1.second + @test g2 ≈ cpu(cg2) + + h = EasyBuild.heisenberg(10) + g1 = expect'(h, reg1=>c) + cg1 = expect'(h, cu(reg1)=>c) + @test g1.first ≈ cpu(cg1.first) + @test g1.second ≈ cg1.second +end + +@testset "apply density matrix" begin + reg = rand_state(6) + creg = cu(reg) + rho = density_matrix(reg, (3,4)) + crho = density_matrix(creg, (3,4)) + @test rho ≈ cpu(crho) + g = put(2, 1=>Rx(0.3)) + @test cpu(apply(crho, g)) ≈ apply(rho, g) + rho = density_matrix(reg) + crho = density_matrix(creg) + @test rho ≈ cpu(crho) + @test probs(crho) ≈ probs(creg) + g = put(6, 1=>Rx(0.3)) + @test cpu(apply(crho, g)) ≈ apply(rho, g) + + # channel + c = UnitaryChannel([put(6, 1=>Rx(0.3)), put(6, 2=>Z)], [0.4, 0.6]) + @test cpu(apply(crho, c)) ≈ apply(rho, c) +end + +@testset "expect on density matrix" begin + reg = rand_state(6) + rho = density_matrix(reg, (3,4,5)) + crho = cu(rho) + h = EasyBuild.heisenberg(3) + a = expect(h, rho) + b = expect(h, crho) + @test a ≈ b + # fidelity + @test fidelity(crho, crho) ≈ 1 + @test measure(crho; nshots=2) isa Vector +end diff --git a/ext/CuYao/test/instructs.jl b/ext/CuYao/test/instructs.jl new file mode 100644 index 000000000..48ae4c696 --- /dev/null +++ b/ext/CuYao/test/instructs.jl @@ -0,0 +1,135 @@ +using LinearAlgebra, Yao.ConstGate +using Test, Random +using Yao +using Yao.YaoArrayRegister.StaticArrays +using Yao.ConstGate: SWAPGate +using CUDA + +@testset "gpu instruct nbit!" begin + Random.seed!(3) + nbit = 6 + N = 1<Z)), + mat(ComplexF32, put(2,2=>I2)), + mat(ComplexF32, put(2,2=>P0)) + ] + + @test instruct!(Val(2),v1 |> CuArray, UN, (3,1)) |> Vector ≈ instruct!(Val(2),v1 |> copy, UN, (3,1)) + @test instruct!(Val(2),vn |> CuArray, UN, (3,1)) |> Matrix ≈ instruct!(Val(2),vn |> copy, UN, (3,1)) + end + # sparse matrix like P0, P1 et. al. are not implemented. + for g in [mat(ComplexF32, ConstGate.T), mat(ComplexF32, ConstGate.H), mat(ComplexF32, rot(Z, 0.5))] + @test instruct!(Val(2), v1 |> CuArray, g, (3,), (4,), (1,)) |> Vector ≈ instruct!(Val(2), v1 |> copy, g, (3,), (4,), (1,)) + @test instruct!(Val(2), vn |> CuArray, g, (3,), (4,), (1,)) |> Matrix ≈ instruct!(Val(2), vn |> copy, g, (3,), (4,), (1,)) + @test instruct!(Val(2), v1 |> CuArray, g, (3,), (4,), (1,)) |> Vector ≈ instruct!(Val(2), v1 |> copy, g, (3,), (4,), (1,)) + @test instruct!(Val(2), vn |> CuArray, g, (3,), (4,), (1,)) |> Matrix ≈ instruct!(Val(2), vn |> copy, g, (3,), (4,), (1,)) + end +end + +@testset "gpu swapapply!" begin + nbit = 6 + N = 1< CuArray, Val(:SWAP), (3,5)) |> Vector ≈ instruct!(Val(2),v1 |> copy, Val(:SWAP), (3,5)) + @test instruct!(Val(2),vn |> CuArray, Val(:SWAP), (3,5)) |> Matrix ≈ instruct!(Val(2),vn |> copy, Val(:SWAP), (3,5)) +end + +@testset "gpu instruct! 1bit" begin + nbit = 6 + N = 1< CuArray, U1, (3,)) |> Vector ≈ instruct!(Val(2),v1 |> copy, U1, (3,)) + @test instruct!(Val(2),vn |> CuArray, U1, (3,)) |> Matrix ≈ instruct!(Val(2),vn |> copy, U1, (3,)) + end + # sparse matrix like P0, P1 et. al. are not implemented. +end + +@testset "gpu xyz-instruct!" begin + nbit = 6 + N = 1< CuArray, Val(G), (3,)) |> Vector ≈ instruct!(Val(2),v1 |> copy, Val(G), (3,)) + @test instruct!(Val(2),vn |> CuArray, Val(G), (3,)) |> Matrix ≈ instruct!(Val(2),vn |> copy, Val(G), (3,)) + if G != :H + @test instruct!(Val(2),v1 |> CuArray, Val(G), (1,3,4)) |> Vector ≈ instruct!(Val(2),v1 |> copy, Val(G), (1,3,4)) + @test instruct!(Val(2),vn |> CuArray, Val(G),(1,3,4)) |> Matrix ≈ instruct!(Val(2),vn |> copy, Val(G), (1,3,4)) + end + end +end + +@testset "gpu cn-xyz-instruct!" begin + nbit = 6 + N = 1< CuArray, Val(G), (3,), (4,5), (0, 1)) |> Vector ≈ instruct!(Val(2),v1 |> copy, Val(G), (3,), (4,5), (0, 1)) + @test instruct!(Val(2),vn |> CuArray, Val(G), (3,), (4,5), (0, 1)) |> Matrix ≈ instruct!(Val(2),vn |> copy, Val(G), (3,), (4,5), (0, 1)) + @test instruct!(Val(2),v1 |> CuArray, Val(G), (3,), (1,), (1,)) |> Vector ≈ instruct!(Val(2),v1 |> copy, Val(G),(3,), (1,), (1,)) + @test instruct!(Val(2),vn |> CuArray, Val(G), (3,), (1,), (1,)) |> Matrix ≈ instruct!(Val(2),vn |> copy, Val(G),(3,), (1,), (1,)) + end +end + +@testset "pswap" begin + ps = put(6, (2,4)=>rot(SWAP, π/2)) + reg = rand_state(6; nbatch=10) + @test apply!(reg |> cu, ps) |> cpu ≈ apply!(copy(reg), ps) + @test apply!(reg |> cu, ps).state isa CuArray +end + +@testset "regression test: Rx, Ry, Rz, CPHASE" begin + Random.seed!(3) + nbit = 6 + for ps in [put(6, (2,)=>Rx(π/2)), put(6, 2=>Ry(0.5)), put(6, 2=>Rz(0.4))] + reg = rand_state(6; nbatch=10) + @test apply!(reg |> cu, ps) |> cpu ≈ apply!(copy(reg), ps) + @test apply!(reg |> cu, ps).state isa CuArray + end +end + +@testset "density matrix" begin + nbit = 6 + reg = rand_state(nbit) + rho = density_matrix(reg) + c = put(6, (2,)=>Rx(π/3)) + @test apply(rho |> cu, c) |> cpu ≈ apply(rho, c) +end + +@testset "time evolve" begin + g = time_evolve(kron(10, 2=>X, 3=>X), 0.5) + reg = rand_state(10) + @test apply!(copy(reg), g) ≈ apply!(reg |> cu, g) |> cpu +end + +# fix: https://github.com/QuantumBFS/CuYao.jl/issues/81 +@testset "generic sparse (pqc circuit)" begin + # kron of Rx + pqc_circuit = subroutine(10, kron(Rx(0.4), Rx(0.5), Rx(0.6), Rx(0.8)), (1, 2, 6, 5)) + proxy_reg = zero_state(10) + @test apply!(proxy_reg |> cu, pqc_circuit) |> cpu ≈ apply(proxy_reg, pqc_circuit) + + # kron of Rz + pqc_circuit = subroutine(10, kron(Rz(0.4), Rz(0.5), Rz(0.6), Rz(0.8)), (1, 2, 6, 5)) + proxy_reg = zero_state(10) + @test apply!(proxy_reg |> cu, pqc_circuit) |> cpu ≈ apply(proxy_reg, pqc_circuit) +end diff --git a/ext/CuYao/test/register.jl b/ext/CuYao/test/register.jl new file mode 100644 index 000000000..8c6dc1d67 --- /dev/null +++ b/ext/CuYao/test/register.jl @@ -0,0 +1,139 @@ +using Test +using Yao +const CuYao = Base.get_extension(Yao, :CuYao) +using LinearAlgebra +using Yao.YaoArrayRegister.BitBasis +using Statistics: mean +using Yao.YaoArrayRegister.StaticArrays +using CUDA +using Yao.YaoArrayRegister: batched_kron, batched_kron!, batch_normalize! +CUDA.allowscalar(false) + +@testset "basics" begin + a = randn(ComplexF64, 50, 20) + ca = a|>cu + batch_normalize!(ca) + batch_normalize!(a) + @test ca |> Matrix ≈ a + + for l = 1:100 + i, j = CuYao.tri2ij(l) + @test (i-1)*(i-2)÷2+j == l + end +end + +@testset "constructor an measure" begin + reg = rand_state(10) + greg = reg |> cu + @test greg isa CuYao.AbstractCuArrayReg + @test eltype(greg.state) == ComplexF64 + myvec(x) = Vector(x) + myvec(x::Number) = [x] + for reg in [rand_state(10, nbatch=333), rand_state(10)] + greg = reg |> cu + @test size(measure(greg |> copy, nshots=10)) == size(measure(reg, nshots=10)) + @test size(measure!(greg |> copy)) == size(measure!(reg |> copy)) + @test size(measure!(ResetTo(0), greg |> copy)) == size(measure!(ResetTo(0), reg |> copy)) + @test size(measure!(RemoveMeasured(), greg |> copy)) == size(measure!(RemoveMeasured(), reg |> copy)) + @test select(greg |> copy, 12) |> cpu ≈ select(reg, 12) + @test size(measure!(greg |> copy |> focus!(3,4,1))) == size(measure!(reg |> copy |> focus!(3,4,1))) + @test greg |> copy |> focus!(3,4,1) |> relax!(3,4,1) |> cpu ≈ reg + + if nbatch(greg) == 1 + greg1 = greg |> copy |> focus!(1,4,3) + greg0 = copy(greg1) + res = measure!(RemoveMeasured(), greg1) + @test select(greg0, res |> myvec) |> normalize! |> cpu ≈ greg1 |> cpu + end + + greg1 = greg |> copy |> focus!(1,4,3) + greg0 = copy(greg1) + res = measure!(ResetTo(3), greg1) + @test all(measure(greg1, nshots=10) .== 3) + @test greg1 |> isnormalized + @test all(select.(BatchedArrayReg(greg0 |> cpu), res |> myvec) .|> normalize! .≈ select.(BatchedArrayReg(greg1 |> cpu), 3)) + + greg1 = greg |> copy |> focus!(1,4,3) + greg0 = copy(greg1) + res = measure!(greg1) + @test all(select.(BatchedArrayReg(greg0 |> cpu), res |> myvec) .|> normalize! .≈ select.(BatchedArrayReg(greg1 |> cpu), res|>myvec)) + end + + @test join(rand_state(3) |> cu, rand_state(3) |> cu) |> nactive == 6 + @test join(rand_state(3, nbatch=10) |> cu, rand_state(3, nbatch=10) |> cu) |> nactive == 6 +end + +@testset "insert_qubits!" begin + reg = rand_state(5; nbatch=10) + res = insert_qudits!(reg |> cu, 3, 2) |> cpu + @test insert_qudits!(reg, 3, 2) ≈ res + @test append_qudits!(copy(reg) |> cu, 3) |> cpu ≈ append_qudits!(copy(reg), 3) + + reg = rand_state(5, nbatch=10) |>focus!(2,3) + res = insert_qudits!(reg |> cu, 3, 2) |> cpu + @test insert_qudits!(reg, 3, 2) ≈ res + + @test append_qudits!(copy(reg) |> cu, 3) |> cpu ≈ append_qudits!(copy(reg), 3) +end + +@testset "cuda-op-measures" begin + reg = rand_state(8; nbatch=32) |> cu + op = repeat(5, X, 1:5) + + # measure! + reg2 = reg |> copy + res = measure!(op, reg2, 2:6) + res2 = measure!(op, reg2, 2:6) + @test size(res) == (32,) + @test res2 == res + + reg = clone(ArrayReg([1,-1+0im]/sqrt(2.0)), 10) |> cu + @test measure!(X, reg) |> mean ≈ -1 + reg = clone(ArrayReg([1.0,0+0im]), 1000) + @test abs(measure!(X, reg) |> mean) < 0.1 +end + +@testset "cuda kron getindex" begin + a = randn(3,4) + b = randn(4,2) + c = zeros(12,8) + ca, cb, cc = cu(a), cu(b), cu(c) + @test kron(ca, cb) |> Array ≈ kron(a, b) + @test Yao.YaoArrayRegister.kron!(cc, ca, cb) |> Array ≈ kron(a,b) + + Yao.YaoArrayRegister.kron!(c,a,b) + @test cc |> Array ≈ c + + v = randn(100) |> cu + inds = [3,5,2,1,7,1] + @test v[inds] ≈ v[inds |> CuVector] +end + +@testset "cuda batched_kron" begin + a = randn(3,4,5) + b = randn(4,2,5) + c = zeros(12,8,5) + ca, cb, cc = cu(a), cu(b), cu(c) + + @test batched_kron(ca, cb) |> Array ≈ batched_kron(a, b) + @test batched_kron!(cc, ca, cb) |> Array ≈ batched_kron(a, b) + + batched_kron!(c, a, b) + @test cc |> Array ≈ c +end + +@testset "zero_state, et al" begin + for b = [4, NoBatch()] + reg = cuzero_state(3; nbatch=b) + @test cpu(reg) ≈ zero_state(3; nbatch=b) && reg isa CuYao.AbstractCuArrayReg + reg = curand_state(3; nbatch=b) + @test reg isa CuYao.AbstractCuArrayReg + reg = cuuniform_state(3; nbatch=b) + @test cpu(reg) ≈ uniform_state(3; nbatch=b) && reg isa CuYao.AbstractCuArrayReg + reg = cuproduct_state(bit"110"; nbatch=b) + @test cpu(reg) ≈ product_state(bit"110"; nbatch=b) && reg isa CuYao.AbstractCuArrayReg + reg = cughz_state(3; nbatch=b) + @test cpu(reg) ≈ ghz_state(3; nbatch=b) && reg isa CuYao.AbstractCuArrayReg + end +end + diff --git a/ext/CuYao/test/runtests.jl b/ext/CuYao/test/runtests.jl new file mode 100644 index 000000000..8dc8f45bb --- /dev/null +++ b/ext/CuYao/test/runtests.jl @@ -0,0 +1,22 @@ +using CUDA, Yao, Test +CUDA.allowscalar(false) + +@testset "CUDA patch" begin + include("CUDApatch.jl") +end + +@testset "GPU reg" begin + include("register.jl") +end + +@testset "gpu applies" begin + include("instructs.jl") +end + +@testset "extra" begin + include("extra.jl") +end + +@testset "yao2einsum" begin + include("yao2einsum.jl") +end diff --git a/ext/CuYao/test/yao2einsum.jl b/ext/CuYao/test/yao2einsum.jl new file mode 100644 index 000000000..60d599d1a --- /dev/null +++ b/ext/CuYao/test/yao2einsum.jl @@ -0,0 +1,11 @@ +using Yao, CUDA, Test +using Yao.YaoToEinsum: uniformsize + +@testset "Yao Extensions" begin + n = 5 + c = EasyBuild.qft_circuit(n) + net = cu(yao2einsum(c)) + m = reshape(net.code(net.tensors...; size_info=uniformsize(net.code, 2)), 1< expect(op, r) -0.7071067811865474 + 0.0im +0.7071067811865474 ``` """ @interface expect diff --git a/lib/YaoAPI/src/registers.jl b/lib/YaoAPI/src/registers.jl index f66ab4118..d29e3d1dc 100644 --- a/lib/YaoAPI/src/registers.jl +++ b/lib/YaoAPI/src/registers.jl @@ -609,8 +609,8 @@ julia> reg |> probs Return the fidelity between two states. Calcuate the fidelity between `r1` and `r2`, if `r1` or `r2` is not pure state -(`nactive(r) != nqudits(r)`), the fidelity is calcuated by purification. See also -[`pure_state_fidelity`](@ref), [`purification_fidelity`](@ref). +(`nactive(r) != nqudits(r)`), the fidelity is calcuated by purification. See also: +http://iopscience.iop.org/article/10.1088/1367-2630/aa6a4b/meta Obtain the gradient with respect to registers and circuit parameters. For pair input `ψ=>circuit`, the returned gradient is a pair of `gψ=>gparams`, @@ -664,6 +664,8 @@ Trace distance is defined as following: \\frac{1}{2} || A - B ||_{\\rm tr} ``` +It takes values between 0 and 1. + ### Examples ```jldoctest; setup=:(using Yao) @@ -672,7 +674,7 @@ julia> reg1 = uniform_state(3); julia> reg2 = zero_state(3); julia> tracedist(reg1, reg2) -1.8708286933869704 +0.9354143466934852 ``` ### References diff --git a/lib/YaoArrayRegister/Project.toml b/lib/YaoArrayRegister/Project.toml index 26a1abc93..157298c01 100644 --- a/lib/YaoArrayRegister/Project.toml +++ b/lib/YaoArrayRegister/Project.toml @@ -1,6 +1,6 @@ name = "YaoArrayRegister" uuid = "e600142f-9330-5003-8abb-0ebd767abc51" -version = "0.9.5" +version = "0.9.9" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" @@ -18,14 +18,14 @@ TupleTools = "9d95972d-f1c8-5527-a6e0-b4b365fa01f6" YaoAPI = "0843a435-28de-4971-9e8b-a9641b2983a8" [compat] -Adapt = "3" -BitBasis = "0.8" +Adapt = "3, 4" +BitBasis = "0.8, 0.9" DocStringExtensions = "0.8, 0.9" LegibleLambdas = "0.3" LuxurySparse = "0.7" MLStyle = "0.4" StaticArrays = "1" -StatsBase = "0.33" +StatsBase = "0.33 - 0.34" TupleTools = "1.3" YaoAPI = "0.4" julia = "1" diff --git a/lib/YaoArrayRegister/src/density_matrix.jl b/lib/YaoArrayRegister/src/density_matrix.jl index 74c294df2..6e3b02e66 100644 --- a/lib/YaoArrayRegister/src/density_matrix.jl +++ b/lib/YaoArrayRegister/src/density_matrix.jl @@ -77,7 +77,7 @@ end YaoAPI.density_matrix(reg::ArrayReg{D}) where D = DensityMatrix{D}(reg.state * reg.state') YaoAPI.density_matrix(rho::DensityMatrix) = copy(rho) -YaoAPI.tracedist(dm1::DensityMatrix{D}, dm2::DensityMatrix{D}) where {D} = trace_norm(dm1.state .- dm2.state) +YaoAPI.tracedist(dm1::DensityMatrix{D}, dm2::DensityMatrix{D}) where {D} = trace_norm(dm1.state .- dm2.state) / 2.0 """ is_density_matrix(dm::AbstractMatrix; kw...) @@ -106,7 +106,7 @@ end function YaoAPI.purify(r::DensityMatrix{D}; num_env::Int = nactive(r)) where {D} Ne = D ^ num_env Ns = size(r.state, 1) - R, U = eigen!(r.state) + R, U = eigen(r.state) state = view(U, :, Ns-Ne+1:Ns) .* sqrt.(abs.(view(R, Ns-Ne+1:Ns)')) return ArrayReg{D}(state) end diff --git a/lib/YaoArrayRegister/src/instruct.jl b/lib/YaoArrayRegister/src/instruct.jl index a3f85ac3d..39432f131 100644 --- a/lib/YaoArrayRegister/src/instruct.jl +++ b/lib/YaoArrayRegister/src/instruct.jl @@ -221,8 +221,10 @@ function YaoAPI.instruct!(::Val{2}, ::Val{:H}, locs::NTuple{N,Int}, ) where {T, N} - instruct!(Val(2), state, matchtype(T, YaoArrayRegister.Const.H), locs) + instruct!(Val(2), state, _hadamard_matrix(T), locs) end +# this is an interface for future extension (e.g. in YaoSym) +_hadamard_matrix(::Type{T}) where T = matchtype(T, YaoArrayRegister.Const.H) function YaoAPI.instruct!(::Val{2}, state::AbstractVecOrMat{T}, diff --git a/lib/YaoArrayRegister/src/operations.jl b/lib/YaoArrayRegister/src/operations.jl index 0f35c450e..ca6922ccc 100644 --- a/lib/YaoArrayRegister/src/operations.jl +++ b/lib/YaoArrayRegister/src/operations.jl @@ -158,16 +158,28 @@ for op in [:*, :/] end end -function Base.:*(bra::AdjointRegister{D,<:ArrayReg}, ket::ArrayReg{D}) where D - if nremain(bra) == nremain(ket) - return dot(relaxedvec(parent(bra)), relaxedvec(ket)) - elseif nremain(bra) == 0 # |remain> - return ArrayReg{D}(state(bra) * state(ket)) - else - error( - "partially contract ⟨bra|ket⟩ is not supported, expect ⟨bra| to be fully actived. nactive(bra)/nqudits(bra)=$(nactive(bra))/$(nqudits(bra))", +""" +$TYPEDSIGNATURES + +The overlap between `ket` and `bra`, which is only defined for two fully activated equal sized registers. +It is only slightly different from [`inner_product`](@ref) in that it always returns a complex number. + +### Examples +```jldoctest; setup=:(using YaoArrayRegister) +julia> reg1 = ghz_state(3); + +julia> reg2 = uniform_state(3); + +julia> reg1' * reg2 +0.5 + 0.0im +``` +""" +function Base.:*(bra::AdjointRegister{D,<:ArrayReg}, ket::ArrayReg{D})::Number where D + # check the register sizes + nqudits(bra) == nqudits(ket) && nremain(bra) == nremain(ket) || error( + "partially contract ⟨bra|ket⟩ is not supported, expect ⟨bra| and |ket⟩ to have the same size. Got nactive(bra)/nqudits(bra)=$(nactive(bra))/$(nqudits(bra)), nactive(ket)/nqudits(ket)=$(nactive(ket))/$(nqudits(ket))", ) - end + return dot(relaxedvec(parent(bra)), relaxedvec(ket)) end Base.:*(bra::AdjointRegister{D,<:BatchedArrayReg{D}}, ket::BatchedArrayReg{D}) where D = bra .* ket @@ -175,23 +187,15 @@ function Base.:*( bra::AdjointRegister{D,<:BatchedArrayReg{D, T1, <:Transpose}}, ket::BatchedArrayReg{D,T2,<:Transpose}, ) where {D,T1,T2} - if nremain(bra) == nremain(ket) == 0 # all active - A, C = parent(state(parent(bra))), parent(state(ket)) - res = zeros(eltype(promote_type(T1, T2)), nbatch(ket)) - #return mapreduce((x, y) -> conj(x) * y, +, ; dims=2) - for j = 1:size(A, 2) - for i = 1:size(A, 1) - @inbounds res[i] += conj(A[i, j]) * C[i, j] - end - end - res - elseif nremain(bra) == 0 # |remain> - bra .* ket - else - error( - "partially contract ⟨bra|ket⟩ is not supported, expect ⟨bra| to be fully actived. nactive(bra)/nqudits(bra)=$(nactive(bra))/$(nqudits(bra))", + nqudits(bra) == nqudits(ket) && nremain(bra) == nremain(ket) && nbatch(bra) == nbatch(ket) || error( + "partially contract ⟨bra|ket⟩ is not supported, expect ⟨bra| and |ket⟩ to have the same size. Got nactive(bra)/nqudits(bra)/nbatch(bra)=$(nactive(bra))/$(nqudits(bra))/$(nbatch(bra)), nactive(ket)/nqudits(ket)=$(nactive(ket))/$(nqudits(ket))/$(nbatch(ket))", ) + A, C = parent(state(parent(bra))), parent(state(ket)) + res = zeros(eltype(promote_type(T1, T2)), nbatch(ket)) + for j in 1:size(A, 2), i in 1:size(A, 1) + @inbounds res[i] += conj(A[i, j]) * C[i, j] end + return res end # broadcast diff --git a/lib/YaoArrayRegister/src/register.jl b/lib/YaoArrayRegister/src/register.jl index 634e6c20f..292539b9c 100644 --- a/lib/YaoArrayRegister/src/register.jl +++ b/lib/YaoArrayRegister/src/register.jl @@ -1,5 +1,10 @@ import BitBasis: BitStr, BitStr64 +""" + AbstractArrayReg + +Abstract type for quantum registers that are represented by an array. +""" abstract type AbstractArrayReg{D,T,AT} <: AbstractRegister{D} end struct NoBatch end @@ -26,8 +31,7 @@ is the numerical type for each amplitude, it is `ComplexF64` by default. !!! warning `ArrayReg` constructor will not normalize the quantum state. If you need a - normalized quantum state remember to use `normalize!(register)` on the register or - normalize the input raw array with `normalize` or [`batched_normalize!`](@ref). + normalized quantum state remember to use `normalize!(register)` on the register. """ mutable struct ArrayReg{D,T,MT<:AbstractMatrix{T}} <: AbstractArrayReg{D,T,MT} state::MT @@ -63,8 +67,7 @@ is the numerical type for each amplitude, it is `ComplexF64` by default. !!! warning `BatchedArrayReg` constructor will not normalize the quantum state. If you need a - normalized quantum state remember to use `normalize!(register)` on the register or - normalize the input raw array with `normalize` or [`batched_normalize!`](@ref). + normalized quantum state remember to use `normalize!(register)` on the register. """ mutable struct BatchedArrayReg{D,T,MT<:AbstractMatrix{T}} <: AbstractArrayReg{D,T,MT} state::MT @@ -145,6 +148,11 @@ end Adapt.@adapt_structure ArrayReg Adapt.@adapt_structure BatchedArrayReg +""" + AdjointArrayReg{D,T,MT} = AdjointRegister{D,<:AbstractArrayReg{D,T,MT}} + +Adjoint array register type, it is used to represent the bra in the Dirac notation. +""" const AdjointArrayReg{D,T,MT} = AdjointRegister{D,<:AbstractArrayReg{D,T,MT}} const ArrayRegOrAdjointArrayReg{D,T,MT} = Union{AbstractArrayReg{D,T,MT},AdjointArrayReg{D,T,MT}} @@ -618,7 +626,7 @@ Create a uniform state: ```math \\frac{1}{\\sqrt{2^n}} \\sum_{k=0}^{2^{n}-1} |k\\rangle. ``` -This state can also be created by applying [`H`](@ref) (Hadmard gate) on ``|00⋯00⟩`` state. +This state can also be created by applying `H` (Hadmard gate) on ``|00⋯00⟩`` state. ### Example diff --git a/lib/YaoArrayRegister/test/operations.jl b/lib/YaoArrayRegister/test/operations.jl index d236f9d39..57b3e4bca 100644 --- a/lib/YaoArrayRegister/test/operations.jl +++ b/lib/YaoArrayRegister/test/operations.jl @@ -51,9 +51,7 @@ end ket = arrayreg(bit"100") + 2 * arrayreg(bit"110") + 3 * arrayreg(bit"111") focus!(ket, 2:3) - t = bra' * ket - relax!(t, 1) - @test state(t) ≈ [1, 0] + @test_throws ErrorException bra' * ket relax!(ket, 2:3) focus!(ket, 1) @@ -66,7 +64,7 @@ end reg1 = rand_state(2; nbatch = 10) reg2 = rand_state(5; nbatch = 10) focus!(reg2, 2:3) - @test all(reg1' * reg2 .≈ reg1' .* reg2) + @test_throws ErrorException reg1' * reg2 end @testset "inplace funcs" begin diff --git a/lib/YaoBlocks/Project.toml b/lib/YaoBlocks/Project.toml index 649493784..74816f704 100644 --- a/lib/YaoBlocks/Project.toml +++ b/lib/YaoBlocks/Project.toml @@ -1,14 +1,14 @@ name = "YaoBlocks" uuid = "418bc28f-b43b-5e0b-a6e7-61bbc1a2c1df" -version = "0.13.8" +version = "0.13.13" [deps] BitBasis = "50ba71b6-fa0f-514d-ae9a-0916efc90dcf" CacheServers = "a921213e-d44a-5460-ac04-5d720a99ba71" ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" -ExponentialUtilities = "d4d017d3-3776-5f7e-afef-a10c40355c18" InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" +KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" LegibleLambdas = "f1f30506-32fe-5131-bd72-7c197988f9e5" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LuxurySparse = "d05aeea4-b7d4-55ac-b691-9e7fabb07ba2" @@ -22,16 +22,20 @@ YaoAPI = "0843a435-28de-4971-9e8b-a9641b2983a8" YaoArrayRegister = "e600142f-9330-5003-8abb-0ebd767abc51" [compat] -BitBasis = "0.8" +BitBasis = "0.8, 0.9" CacheServers = "0.2" ChainRulesCore = "1.11" DocStringExtensions = "0.8, 0.9" -ExponentialUtilities = "1.5" +InteractiveUtils = "1" +KrylovKit = "0.5, 0.6, 0.7" LegibleLambdas = "0.2, 0.3" +LinearAlgebra = "1" LuxurySparse = "0.7" MLStyle = "0.3, 0.4" +Random = "1" +SparseArrays = "1" StaticArrays = "0.12, 1.0" -StatsBase = "0.32, 0.33" +StatsBase = "0.32 - 0.34" TupleTools = "1.2" YaoAPI = "0.4" YaoArrayRegister = "0.9.5" diff --git a/lib/YaoBlocks/src/YaoBlocks.jl b/lib/YaoBlocks/src/YaoBlocks.jl index ff80f70e9..eace214ee 100644 --- a/lib/YaoBlocks/src/YaoBlocks.jl +++ b/lib/YaoBlocks/src/YaoBlocks.jl @@ -14,7 +14,8 @@ using StatsBase, TupleTools, InteractiveUtils using MLStyle: @match using LinearAlgebra: eigen! using SparseArrays, LuxurySparse -using ExponentialUtilities, Random, CacheServers +using Random, CacheServers +import KrylovKit: exponentiate using DocStringExtensions import StaticArrays: SMatrix @@ -50,7 +51,9 @@ import YaoAPI: subblocks, nparameters, measure!, - measure + measure, + apply_back!, + mat_back! export AbstractBlock, AbstractContainer, @@ -69,6 +72,7 @@ export AbstractBlock, content, dispatch!, dispatch, + sandwich, expect, getiparams, iparams_eltype, diff --git a/lib/YaoBlocks/src/abstract_block.jl b/lib/YaoBlocks/src/abstract_block.jl index c8605aac8..517bb16ca 100644 --- a/lib/YaoBlocks/src/abstract_block.jl +++ b/lib/YaoBlocks/src/abstract_block.jl @@ -339,7 +339,7 @@ render_params(r::AbstractBlock, ::Val{:zero}) = """ cache_type(::Type) -> DataType -Return the element type that a [`CacheFragment`](@ref) +Return the element type that a `CacheFragment` will use. """ cache_type(::Type{<:AbstractBlock}) = Any @@ -420,7 +420,6 @@ end The non-inplace version of applying a block (of quantum circuit) to a quantum register. Check `apply!` for the faster inplace version. """ -# overwrite interface, this one avoids one copy function apply(reg::AbstractRegister{D}, block) where D return apply!(copy(reg), block) end diff --git a/lib/YaoBlocks/src/autodiff/apply_back.jl b/lib/YaoBlocks/src/autodiff/apply_back.jl index c280f829c..32ff2c3aa 100644 --- a/lib/YaoBlocks/src/autodiff/apply_back.jl +++ b/lib/YaoBlocks/src/autodiff/apply_back.jl @@ -164,7 +164,7 @@ end The backward function of `apply!`. Returns a tuple of ((input register, gradient of input register), parameter gradients) """ -function apply_back(st::Tuple{<:AbstractArrayReg,<:AbstractArrayReg}, block::AbstractBlock; kwargs...) +function apply_back(st::Tuple{<:AbstractArrayReg, <:AbstractArrayReg}, block::AbstractBlock; kwargs...) col = [] in, inδ = apply_back!(st, block, col; kwargs...) (in, inδ), col diff --git a/lib/YaoBlocks/src/autodiff/autodiff.jl b/lib/YaoBlocks/src/autodiff/autodiff.jl index 7cbf3b935..a77e0302e 100644 --- a/lib/YaoBlocks/src/autodiff/autodiff.jl +++ b/lib/YaoBlocks/src/autodiff/autodiff.jl @@ -15,6 +15,7 @@ using BitBasis, YaoArrayRegister, YaoAPI using ..YaoBlocks import ChainRulesCore: rrule, @non_differentiable, NoTangent, Tangent, backing, AbstractTangent, ZeroTangent +import YaoAPI: mat_back!, apply_back! using SparseArrays, LuxurySparse, LinearAlgebra include("NoParams.jl") diff --git a/lib/YaoBlocks/src/autodiff/chainrules_patch.jl b/lib/YaoBlocks/src/autodiff/chainrules_patch.jl index 5cab81016..6895feccd 100644 --- a/lib/YaoBlocks/src/autodiff/chainrules_patch.jl +++ b/lib/YaoBlocks/src/autodiff/chainrules_patch.jl @@ -79,17 +79,18 @@ end function rrule(::typeof(apply), reg::AbstractArrayReg, block::AbstractBlock) out = apply(reg, block) out, function (outδ) - (in, inδ), paramsδ = apply_back((copy(out), outδ), block) + (in, inδ), paramsδ = apply_back((copy(out), tangent_to_reg(typeof(out), outδ)), block) return (NoTangent(), inδ, create_circuit_tangent(block, paramsδ)) end end function rrule(::typeof(apply), reg::AbstractArrayReg, block::AbstractAdd) out = apply(reg, block) out, function (outδ) - (in, inδ), paramsδ = apply_back((copy(out), outδ), block; in = reg) + (in, inδ), paramsδ = apply_back((copy(out), tangent_to_reg(typeof(out), outδ)), block; in = reg) return (NoTangent(), inδ, create_circuit_tangent(block, paramsδ)) end end +tangent_to_reg(::Type{T}, reg) where T<:AbstractArrayReg = reg isa Tangent ? T(reg.state) : reg function rrule(::typeof(dispatch), block::AbstractBlock, params) diff --git a/lib/YaoBlocks/src/blocktools.jl b/lib/YaoBlocks/src/blocktools.jl index 261fa7ceb..df3bdcfe8 100644 --- a/lib/YaoBlocks/src/blocktools.jl +++ b/lib/YaoBlocks/src/blocktools.jl @@ -74,9 +74,9 @@ collect_blocks(::Type{T}, x::AbstractBlock) where {T<:AbstractBlock} = #expect(op::AbstractBlock, dm::DensityMatrix) = mapslices(x->sum(mat(op).*x)[], dm.state, dims=[1,2]) |> vec """ - expect(op::AbstractBlock, reg) -> Vector - expect(op::AbstractBlock, reg => circuit) -> Vector - expect(op::AbstractBlock, density_matrix) -> Vector + expect(op::AbstractBlock, reg) -> Real + expect(op::AbstractBlock, reg => circuit) -> Real + expect(op::AbstractBlock, density_matrix) -> Real Get the expectation value of an operator, the second parameter can be a register `reg` or a pair of input register and circuit `reg => circuit`. @@ -92,13 +92,21 @@ For register input, the return value is a register. For batched register, `expect(op, reg=>circuit)` returns a vector of size number of batch as output. However, one can not differentiate over a vector loss, so `expect'(op, reg=>circuit)` accumulates the gradient over batch, rather than returning a batched gradient of parameters. """ +function expect(op::AbstractBlock, reg::AbstractRegister) + # NOTE: broadcast because the input register can be a batched one + return safe_real.(sandwich(reg, op, reg)) +end +function expect(op, plan::Pair{<:AbstractRegister,<:AbstractBlock}) + expect(op, copy(plan.first) |> plan.second) +end + function expect(op::AbstractBlock, dm::DensityMatrix) # NOTE: we use matrix form here because the matrix size is known to be small, # while applying a circuit on a reduced density matrix might take much more than constructing the matrix. mop = mat(op) # TODO: switch to `IterNz` # sum(x->dm.state[x[2],x[1]]*x[3], IterNz(mop)) - return sum(transpose(dm.state) .* mop) + return safe_real(sum(transpose(dm.state) .* mop)) end function expect(op::AbstractAdd, reg::DensityMatrix) # NOTE: this is faster in e.g. when the op is Heisenberg @@ -108,23 +116,28 @@ function expect(op::Scale, reg::DensityMatrix) factor(op) * expect(content(op), reg) end -# NOTE: assume an register has a bra. Can we define it for density matrix? -expect(op::AbstractBlock, reg::AbstractRegister) = reg' * apply!(copy(reg), op) +""" + sandwich(bra::AbstractRegister, op::AbstractBlock, ket::AbstracRegister) -> Complex + +Compute the sandwich function ⟨bra|op|ket⟩. +""" +sandwich(bra::AbstractArrayReg, op::AbstractBlock, reg::AbstractArrayReg) = bra' * apply!(copy(reg), op) -function expect(op::AbstractBlock, reg::BatchedArrayReg) +function sandwich(bra::BatchedArrayReg, op::AbstractBlock, reg::BatchedArrayReg) + @assert nbatch(bra) == nbatch(reg) B = YaoArrayRegister._asint(nbatch(reg)) ket = apply!(copy(reg), op) - if !(reg.state isa Transpose) # not-transposed storage + if !(bra.state isa Transpose) # not-transposed storage C = reshape(ket.state, :, B) - A = reshape(reg.state, :, B) + A = reshape(bra.state, :, B) # reduce over the 1st dimension conjsumprod1(A, C) - elseif size(reg.state, 2) == B # transposed storage, no environment qubits + elseif size(bra.state, 2) == B # transposed storage, no environment qubits # reduce over the second dimension - conjsumprod2(reg.state.parent, ket.state.parent) + conjsumprod2(bra.state.parent, ket.state.parent) else - C = reshape(ket.state.parent, :, B, size(reg.state, 1)) - A = reshape(reg.state.parent, :, B, size(reg.state, 1)) + C = reshape(ket.state.parent, :, B, size(bra.state, 1)) + A = reshape(bra.state.parent, :, B, size(bra.state, 1)) # reduce over the 1st and 3rd dimension conjsumprod13(A, C) end @@ -161,19 +174,15 @@ function conjsumprod13(A::AbstractArray, C::AbstractArray) res end -for REG in [:AbstractRegister, :BatchedArrayReg] - @eval function expect(op::AbstractAdd, reg::$REG) - sum(opi -> expect(opi, reg), op) +for REG in [:AbstractArrayReg, :BatchedArrayReg] + @eval function sandwich(bra::$REG, op::AbstractAdd, reg::$REG) + sum(opi -> sandwich(bra, opi, reg), op) end - @eval function expect(op::Scale, reg::$REG) - factor(op) * expect(content(op), reg) + @eval function sandwich(bra::$REG, op::Scale, reg::$REG) + factor(op) * sandwich(bra, content(op), reg) end end -function expect(op, plan::Pair{<:AbstractRegister,<:AbstractBlock}) - expect(op, copy(plan.first) |> plan.second) -end - # obtaining Dense Matrix of a block LinearAlgebra.Matrix(blk::AbstractBlock) = Matrix(mat(blk)) @@ -183,7 +192,7 @@ LinearAlgebra.Matrix(blk::AbstractBlock) = Matrix(mat(blk)) Operator fidelity defined as ```math -F^2 = \\frac{1}{d^2}\\left[{\\rm Tr}(b1^\\dagger b2)\\right] +F^2 = \\frac{1}{d}\\left[{\\rm Tr}(b1^\\dagger b2)\\right] ``` Here, `d` is the size of the Hilbert space. Note this quantity is independant to global phase. diff --git a/lib/YaoBlocks/src/channel/unitary_channel.jl b/lib/YaoBlocks/src/channel/unitary_channel.jl index ab89d1657..dac24bcc1 100644 --- a/lib/YaoBlocks/src/channel/unitary_channel.jl +++ b/lib/YaoBlocks/src/channel/unitary_channel.jl @@ -50,6 +50,7 @@ function YaoAPI.unsafe_apply!(r::DensityMatrix{D,T}, for (locs, block) in zip(k.locs, k.blocks) YaoAPI.unsafe_apply!(r, put(k.n, locs => block)) end + return r end function mat(::Type{T}, x::UnitaryChannel) where {T} diff --git a/lib/YaoBlocks/src/composite/control.jl b/lib/YaoBlocks/src/composite/control.jl index 8e87fea72..12f000da5 100644 --- a/lib/YaoBlocks/src/composite/control.jl +++ b/lib/YaoBlocks/src/composite/control.jl @@ -1,5 +1,16 @@ export ControlBlock, control, cnot, cz +""" + $(TYPEDSIGNATURES) + +A control block is a composite block that applies a block when the control qubits are all ones. + +!!! note + If control qubit index is negative, it means the inverse control, i.e., the block is applied when the control qubit is zero. + +### Fields +$(TYPEDFIELDS) +""" struct ControlBlock{BT<:AbstractBlock,C,M} <: AbstractContainer{BT,2} n::Int ctrl_locs::NTuple{C,Int} diff --git a/lib/YaoBlocks/src/composite/kron.jl b/lib/YaoBlocks/src/composite/kron.jl index 34f653b4f..094081d44 100644 --- a/lib/YaoBlocks/src/composite/kron.jl +++ b/lib/YaoBlocks/src/composite/kron.jl @@ -218,7 +218,7 @@ Base.length(k::KronBlock) = length(k.blocks) Base.eachindex(k::KronBlock) = k.locs function Base.:(==)(lhs::KronBlock, rhs::KronBlock) - return nqudits(lhs) == nqudits(rhs) && all(lhs.locs .== rhs.locs) && all(lhs.blocks .== rhs.blocks) + return nqudits(lhs) == nqudits(rhs) && length(lhs.locs) == length(rhs.locs) && all(lhs.locs .== rhs.locs) && all(lhs.blocks .== rhs.blocks) end Base.adjoint(blk::KronBlock) = KronBlock(blk.n, blk.locs, map(adjoint, blk.blocks)) diff --git a/lib/YaoBlocks/src/composite/put_block.jl b/lib/YaoBlocks/src/composite/put_block.jl index cc01e8aa3..493dba57f 100644 --- a/lib/YaoBlocks/src/composite/put_block.jl +++ b/lib/YaoBlocks/src/composite/put_block.jl @@ -119,7 +119,24 @@ function YaoAPI.iscommute(x::PutBlock{D}, y::PutBlock{D}) where {D} end end +""" + Swap = PutBlock{2,2,G} where {G<:ConstGate.SWAPGate} + Swap(n::Int, locs::Tuple{Int,Int}) + +Swap gate, which swaps two qubits. +""" const Swap = PutBlock{2,2,G} where {G<:ConstGate.SWAPGate} + +""" + PSwap = PutBlock{2,2,RotationGate{2,T,G}} where {G<:ConstGate.SWAPGate} + PSwap(n::Int, locs::Tuple{Int,Int}, θ::Real) + +Parametrized swap gate that swaps two qubits with a phase, defined as + +```math +{\\rm SWAP}(θ) = e^{-iθ{\\rm SWAP}/2} +``` +""" const PSwap{T} = PutBlock{2,2,RotationGate{2,T,G}} where {G<:ConstGate.SWAPGate} Swap(n::Int, locs::Tuple{Int,Int}) = PutBlock(n, ConstGate.SWAPGate(), locs) PSwap(n::Int, locs::Tuple{Int,Int}, θ::Real) = diff --git a/lib/YaoBlocks/src/composite/repeated.jl b/lib/YaoBlocks/src/composite/repeated.jl index a49006ff5..eed9e79c2 100644 --- a/lib/YaoBlocks/src/composite/repeated.jl +++ b/lib/YaoBlocks/src/composite/repeated.jl @@ -12,13 +12,13 @@ struct RepeatedBlock{D,C,GT<:AbstractBlock} <: AbstractContainer{GT,D} end function RepeatedBlock(n::Int, block::AbstractBlock{D}, locs::NTuple{C,Int}) where {D,C} - @assert_locs_safe n Tuple(i:i+nqudits(block)-1 for i in locs) + @assert_locs_safe n [i:i+nqudits(block)-1 for i in locs] nqudits(block) > 1 && throw( ArgumentError("RepeatedBlock does not support multi-qubit content for the moment."), ) # sort the locations - locs = TupleTools.sort(locs) - return RepeatedBlock{D,C,typeof(block)}(n, block, locs) + slocs = sort([locs...]) + return RepeatedBlock{D,C,typeof(block)}(n, block, (slocs...,)) end function RepeatedBlock(n::Int, block::AbstractBlock{D}, locs::UnitRange{Int}) where {D} @@ -105,7 +105,7 @@ julia> @time apply!(reg, chain([put(20, i=>X) for i=1:20])); Base.repeat(n::Int, x::AbstractBlock, locs::Int...) = repeat(n, x, locs) Base.repeat(n::Int, x::AbstractBlock, locs::NTuple{C,Int}) where {C} = RepeatedBlock(n, x, locs) -Base.repeat(n::Int, x::AbstractBlock, locs) = repeat(n, x, locs...) +Base.repeat(n::Int, x::AbstractBlock, locs) = repeat(n, x, (locs...,)) Base.repeat(n::Int, x::AbstractBlock, locs::UnitRange) = RepeatedBlock(n, x, locs) Base.repeat(n::Int, x::AbstractBlock) = RepeatedBlock(n, x) Base.repeat(x::AbstractBlock) = @λ(n -> repeat(n, x)) @@ -128,9 +128,8 @@ mat(::Type{T}, rb::RepeatedBlock{D}) where {T,D} = mat(::Type{T}, rb::RepeatedBlock{D,0,GT}) where {T,D,GT} = IMatrix{T}(D^nqudits(rb)) function YaoAPI.unsafe_apply!(r::AbstractRegister, rp::RepeatedBlock) - m = mat_matchreg(r, rp.content) for addr in rp.locs - instruct!(r, m, Tuple(addr:addr+nqudits(rp.content)-1)) + YaoAPI.unsafe_apply!(r, put(nqudits(rp), addr:addr+nqudits(rp.content)-1=>rp.content)) end return r end diff --git a/lib/YaoBlocks/src/layout.jl b/lib/YaoBlocks/src/layout.jl index 35eb5f930..9cfa764df 100644 --- a/lib/YaoBlocks/src/layout.jl +++ b/lib/YaoBlocks/src/layout.jl @@ -77,7 +77,7 @@ Print the block tree. # Keywords - `maxdepth`: max tree depth to print -- `charset`: default is ('├','└','│','─'). See also [`BlockTreeCharSet`](@ref). +- `charset`: default is ('├','└','│','─'). See also `YaoBlocks.BlockTreeCharSet`. - `title`: control whether to print the title, `true` or `false`, default is `true` """ function print_tree( diff --git a/lib/YaoBlocks/src/primitive/reflect.jl b/lib/YaoBlocks/src/primitive/reflect.jl index 7b1268aa4..f305fd2a2 100644 --- a/lib/YaoBlocks/src/primitive/reflect.jl +++ b/lib/YaoBlocks/src/primitive/reflect.jl @@ -1,18 +1,22 @@ export ReflectGate, reflect -ReflectGate{D, T, Tt, AT<:AbstractArrayReg{D, T}} = TimeEvolution{D,Tt,Projector{D,T,AT}} - """ -$(TYPEDSIGNATURES) + ReflectGate{D, T, Tt, AT<:AbstractArrayReg{D, T}} = TimeEvolution{D,Tt,Projector{D,T,AT}} -Create a [`ReflectGate`](@ref) with respect to an quantum state vector `v`. -It defines the following gate operation. +Let `|v⟩` be a quantum state vector, a reflection gate is a unitary operator that defined as the following operation. ```math |v⟩ → 1 - (1-exp(-iθ)) |v⟩⟨v| ``` When ``θ = π``, it defines a standard reflection gate ``1-2|v⟩⟨v|``. +""" +ReflectGate{D, T, Tt, AT<:AbstractArrayReg{D, T}} = TimeEvolution{D,Tt,Projector{D,T,AT}} + +""" +$(TYPEDSIGNATURES) + +Create a [`ReflectGate`](@ref) with respect to an quantum state vector `v`. ### Example diff --git a/lib/YaoBlocks/src/primitive/time_evolution.jl b/lib/YaoBlocks/src/primitive/time_evolution.jl index 63c0ff3f5..7f11df33c 100644 --- a/lib/YaoBlocks/src/primitive/time_evolution.jl +++ b/lib/YaoBlocks/src/primitive/time_evolution.jl @@ -87,10 +87,12 @@ function LinearAlgebra.mul!(y::AbstractVector, A::BlockMap{T,GT}, x::AbstractVec return y end +compatible_multiplicative_operand(::AbstractArray, source::AbstractArray) = source + function YaoAPI.unsafe_apply!(reg::AbstractArrayReg{D,T}, te::TimeEvolution) where {D,T} if isdiagonal(te.H) # `compatible_multiplicative_operand` is used to ensure GPU compatibility. - reg.state .*= exp.((-im * te.dt) .* ExponentialUtilities.compatible_multiplicative_operand(reg.state, diag(mat(T, te.H)))) + reg.state .*= exp.((-im * te.dt) .* compatible_multiplicative_operand(reg.state, diag(mat(T, te.H)))) return reg end st = state(reg) @@ -98,10 +100,14 @@ function YaoAPI.unsafe_apply!(reg::AbstractArrayReg{D,T}, te::TimeEvolution) whe A = BlockMap(T, te.H) @inbounds for j = 1:size(st, 2) v = view(st, :, j) - # TODO: opnorm could be estimated by inspecting the block operation. - # We should fix it in the future in improve the accuracy. - Ks = arnoldi(A, v; tol = te.tol, ishermitian = true, opnorm = 1.0) - expv!(v, dt, Ks) + v_out, info = exponentiate(A, dt, v, + tol=te.tol, + krylovdim=min(1000, size(A,1)), + ishermitian=true, + eager=true,) + # info.converged is 1 if it converged and 0 otherwise + Bool(info.converged) || @warn "Application of $(typeof(te)) did not converge. error = $(info.normres) but tol=$(te.tol)" + v .= v_out end return reg end @@ -124,11 +130,8 @@ function Base.adjoint(te::TimeEvolution) end Base.copy(te::TimeEvolution) = TimeEvolution(te.H, te.dt, tol = te.tol) -YaoAPI.isdiagonal(r::TimeEvolution) = isdiagonal(r.H) -function YaoAPI.isunitary(te::TimeEvolution) - iszero(imag(te.dt)) || return false - return true -end +YaoAPI.isdiagonal(te::TimeEvolution) = isdiagonal(te.H) +YaoAPI.isunitary(te::TimeEvolution) = iszero(imag(te.dt)) iparams_range(::TimeEvolution{D,T}) where {D,T} = ((typemin(T), typemax(T)),) diff --git a/lib/YaoBlocks/src/treeutils/optimise.jl b/lib/YaoBlocks/src/treeutils/optimise.jl index 6bb934438..5ecbba8f1 100644 --- a/lib/YaoBlocks/src/treeutils/optimise.jl +++ b/lib/YaoBlocks/src/treeutils/optimise.jl @@ -214,7 +214,7 @@ const __default_simplification_rules__ = simplify(block[; rules=__default_simplification_rules__]) Simplify a block tree accroding to given rules, default to use -[`__default_simplification_rules__`](@ref). +`YaoBlocks.Optimise.__default_simplification_rules__`. """ function simplify(ex::AbstractBlock; rules = __default_simplification_rules__) out1 = simplify_pass(rules, ex) diff --git a/lib/YaoBlocks/src/utils.jl b/lib/YaoBlocks/src/utils.jl index e7aac8342..357f61178 100644 --- a/lib/YaoBlocks/src/utils.jl +++ b/lib/YaoBlocks/src/utils.jl @@ -311,3 +311,13 @@ end SparseArrays.sparse(et::EntryTable) = SparseVector(et) Base.vec(et::EntryTable) = Vector(et) + +# convert a (maybe complex) number x to real number. +function safe_real(x::Complex{T}) where T + img = imag(x) + if !iszero(img) && !(hasmethod(eps, Tuple{Type{T}}) && isapprox(x - im*img, x; rtol=10*eps(T), atol=10*eps(T))) + @warn "Convert number $x to real is unsafe due to the nonzero imaginary part: $img. Maybe you want to use the `sandwich` function?" + end + return real(x) +end +safe_real(x) = x \ No newline at end of file diff --git a/lib/YaoBlocks/test/autodiff/chainrules_patch.jl b/lib/YaoBlocks/test/autodiff/chainrules_patch.jl index d58acfe34..baba15b9d 100644 --- a/lib/YaoBlocks/test/autodiff/chainrules_patch.jl +++ b/lib/YaoBlocks/test/autodiff/chainrules_patch.jl @@ -272,4 +272,20 @@ end return real((reg' * apply(reg, cost))^2) end @test Zygote.gradient(_theta -> loss(_theta, alpha, var, 6), params)[1] ≈ gt1 +end + +@testset "zygote apply back" begin + function loss3(c, θs) + # the assign is nessesary! + c = dispatch(c, θs) + reg = zero_state(nqubits(c)) + reg = apply(reg, c) + v = vec(reg.state) + real(sum(abs2.(v .* v))) + end + N = 5 + c = chain(rand([[control(N, i, mod1(i+1,N)=>Rx(0.0)) for i=1:N]..., [put(N, i=>Ry(0.0)) for i=1:N]..., [put(N, i=>Rx(0.0)) for i=1:N]...], 6)) + loss3(c, rand(nparameters(c))) + zygote_grad = Zygote.gradient(θs->loss3(c, θs), rand(nparameters(c)))[1] + @test zygote_grad isa Vector end \ No newline at end of file diff --git a/lib/YaoBlocks/test/blocktools.jl b/lib/YaoBlocks/test/blocktools.jl index b818e7ac5..1328a70b2 100644 --- a/lib/YaoBlocks/test/blocktools.jl +++ b/lib/YaoBlocks/test/blocktools.jl @@ -96,4 +96,4 @@ end @test qc |> gatecount |> values |> sum == 6 res = gatecount(repeat(5, X, (2, 3))) @test res |> values |> sum == 2 -end +end \ No newline at end of file diff --git a/lib/YaoBlocks/test/channel/unitary_channel.jl b/lib/YaoBlocks/test/channel/unitary_channel.jl index acfc7dad9..b2c4a46df 100644 --- a/lib/YaoBlocks/test/channel/unitary_channel.jl +++ b/lib/YaoBlocks/test/channel/unitary_channel.jl @@ -49,4 +49,14 @@ end ops = unitary_channel([put(2, 1=>chain(Rx(0.4), Ry(0.5), Rz(0.4))), put(2,2=>Y), put(2,1=>Z)], [0.3, 0.1, 0.6]) channel2 = put(6, (2,3)=>ops) @test apply(r, channel1) ≈ apply(r, channel2) -end \ No newline at end of file +end + +@testset "kron and repeat" begin + dm = density_matrix(ghz_state(5), 1:3) + dpolarizing = single_qubit_depolarizing_channel(0.1) + dmkron = apply(dm, kron(3, 1=>dpolarizing, 2=>dpolarizing)) + dmrepeat = apply(dm, repeat(3, dpolarizing, (1, 2))) + dmput = apply(apply(dm, put(3, 1=>dpolarizing)), put(3, 2=>dpolarizing)) + @test dmkron ≈ dmrepeat + @test dmkron ≈ dmput +end diff --git a/lib/YaoBlocks/test/measure_ops.jl b/lib/YaoBlocks/test/measure_ops.jl index 6b26682d6..9baa6c98a 100644 --- a/lib/YaoBlocks/test/measure_ops.jl +++ b/lib/YaoBlocks/test/measure_ops.jl @@ -51,7 +51,7 @@ end @show op @test isapprox( sum(measure(op, reg2; nshots = 100000)) / 100000, - expect(op, reg), + sandwich(reg, op, reg), rtol = 0.1, ) @test reg ≈ reg2 @@ -60,7 +60,7 @@ end reg2 = copy(reg) @test isapprox( dropdims(sum(measure(op, reg2; nshots = 100000), dims = 1), dims = 1) / 100000, - expect(op, reg), + sandwich(reg, op, reg), rtol = 0.1, ) @test reg ≈ reg2 @@ -91,7 +91,7 @@ end @show op @test isapprox( sum(measure(op, reg2, locs; nshots = 100000)) / 100000, - expect(put(Nbit, locs => op), reg), + sandwich(reg, put(Nbit, locs => op), reg), rtol = 0.2, ) @test reg ≈ reg2 diff --git a/lib/YaoBlocks/test/utils.jl b/lib/YaoBlocks/test/utils.jl index 6a703381b..ae3ba0cee 100644 --- a/lib/YaoBlocks/test/utils.jl +++ b/lib/YaoBlocks/test/utils.jl @@ -45,4 +45,13 @@ end @test isclean(cet) @test cet[bit"000"] == 1.0+0im @test cet[bit"111"] == 0.0im +end + +@testset "#issue 508" begin + @test YaoBlocks.safe_real(0 + 1e-323im) ≈ 0 + @test YaoBlocks.safe_real(0+0im) ≈ 0 + @test_warn () YaoBlocks.safe_real(0+1e-15im) + @test YaoBlocks.safe_real(0+1e-17im) == 0 + @test YaoBlocks.safe_real(1e10+1e-10im) == 1e10 + @test_warn () YaoBlocks.safe_real(0+1im) end \ No newline at end of file diff --git a/lib/YaoPlots/LICENSE b/lib/YaoPlots/LICENSE new file mode 100644 index 000000000..4eedc0116 --- /dev/null +++ b/lib/YaoPlots/LICENSE @@ -0,0 +1,201 @@ +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 [yyyy] [name of copyright owner] + + 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/YaoPlots/Project.toml b/lib/YaoPlots/Project.toml new file mode 100644 index 000000000..7e53dc1e5 --- /dev/null +++ b/lib/YaoPlots/Project.toml @@ -0,0 +1,24 @@ +name = "YaoPlots" +uuid = "32cfe2d9-419e-45f2-8191-2267705d8dbc" +version = "0.9.2" + +[deps] +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Luxor = "ae8d54c2-7ccd-5906-9d76-62fc9837b5bc" +Thebes = "8b424ff8-82f5-59a4-86a6-de3761897198" +YaoArrayRegister = "e600142f-9330-5003-8abb-0ebd767abc51" +YaoBlocks = "418bc28f-b43b-5e0b-a6e7-61bbc1a2c1df" + +[compat] +LinearAlgebra = "1" +Luxor = "3" +Thebes = "1" +YaoArrayRegister = "0.9" +YaoBlocks = "0.13" +julia = "1" + +[extras] +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[targets] +test = ["Test"] diff --git a/lib/YaoPlots/README.md b/lib/YaoPlots/README.md new file mode 100644 index 000000000..ec5d9cd7a --- /dev/null +++ b/lib/YaoPlots/README.md @@ -0,0 +1,61 @@ +# YaoPlots + +[![CI](https://github.com/QuantumBFS/YaoPlots.jl/actions/workflows/CI.yml/badge.svg)](https://github.com/QuantumBFS/YaoPlots.jl/actions/workflows/CI.yml) +[![Coverage](https://codecov.io/gh/QuantumBFS/YaoPlots.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/QuantumBFS/YaoPlots.jl) + +## Example 1: Visualize a QBIR define in Yao + +```julia +using Yao.EasyBuild, YaoPlots + +# show a qft circuit +vizcircuit(qft_circuit(5)) +``` + +If you are using a Pluto/Jupyter notebook, Atom/VSCode editor, you should see the following image in your plotting panel. + +![qft](examples/qft.png) + +## Example 2: Visualize a single qubit +```julia +using YaoPlots, Yao + +reg = zero_state(1) |> Rx(π/8) |> Rx(π/8) +rho = density_matrix(ghz_state(2), 1) + +bloch_sphere("|ψ⟩"=>reg, "ρ"=>rho; show_projection_lines=true) +``` + +Similarly, you will see +![bloch](examples/bloch.png) + +See more [examples](examples/circuits.jl). + +### Adjusting the plot attributes + +Various attributes of the visualizations can be altered. +The plot can be modified, if we change the following attributes + +- `YaoPlots.CircuitStyles.linecolor[]` for line color, default value being `"#000000"` (black color) +- `YaoPlots.CircuitStyles.gate_bgcolor[]` for background color of square blocks, the default value being `"#FFFFFF"` (white color) +- `YaoPlots.CircuitStyles.textcolor[]` for text color, default value being `"#000000` +- `YaoPlots.CircuitStyles.lw[]` for line width, default value being `1` (pt) +- `YaoPlots.CircuitStyles.textsize[]` for text size, default value being `16` (pt) +- `YaoPlots.CircuitStyles.paramtextsize[]` for parameter text size, for parameterized gates, default value being `10` (pt) + +For example, + +```julia +using YaoPlots, Yao +YaoPlots.CircuitStyles.linecolor[] = "pink" +YaoPlots.CircuitStyles.gate_bgcolor[] = "yellow" +YaoPlots.CircuitStyles.textcolor[] = "#000080" # the navy blue color +YaoPlots.CircuitStyles.fontfamily[] = "JuliaMono" +YaoPlots.CircuitStyles.lw[] = 2.5 +YaoPlots.CircuitStyles.textsize[] = 13 +YaoPlots.CircuitStyles.paramtextsize[] = 8 + +vizcircuit(chain(3, put(1=>X), repeat(3, H), put(2=>Y), repeat(3, Rx(π/2)))) +``` + +![attribute_example_2](examples/attr_circuit_2.svg) \ No newline at end of file diff --git a/lib/YaoPlots/examples/bloch.jl b/lib/YaoPlots/examples/bloch.jl new file mode 100644 index 000000000..3ad1e129b --- /dev/null +++ b/lib/YaoPlots/examples/bloch.jl @@ -0,0 +1,54 @@ +using YaoPlots, Yao +using LinearAlgebra: axpy! + +function commutator(x, y; ishermitian=false) + res = x * y + if ishermitian + return res - res' + else + return res - y * x + end +end +function anticommutator(x, y; ishermitian=false) + res = x * y + if ishermitian + return res + res' + else + return res + y * x + end +end + +# single step master equation +function mestep(rho::DensityMatrix{D}, h, Ls, dt) where D + res = copy(rho.state) + # The im*[ρ, H] term. + # NOTE: transposed storage is faster + reg = arrayreg(copy(rho.state)'; nlevel=size(rho.state, 2)) + apply!(reg, h) + axpy!(dt * im, reg.state' - reg.state, res) + for L in Ls + # the LρL' term + axpy!(dt, apply(rho, L).state, res) + # the -(ρL'L + L'Lρ)/2 term + reg = arrayreg(copy(rho.state)'; nlevel=size(rho.state, 2)) + apply!(reg, L' * L) + axpy!(-0.5dt, reg.state, res) + axpy!(-0.5dt, reg.state', res) + end + return DensityMatrix(res) +end + +function simulate(t; dt=0.02, dissiplation=0.2, filename=nothing) + reg0 = zero_state(1) |> H + rho = density_matrix(reg0) + h = Z # rotate around Z + Ls = [sqrt(dissiplation) * ConstGate.Pd] # dissipation + states = ["|+⟩"=>rho] + for t = 0:dt:t + rho = mestep(rho, h, Ls, dt) + push!(states, ""=>rho) + end + return bloch_sphere(states...; filename) +end + +simulate(π/2; filename=joinpath(@__DIR__, "mesolve.png")) \ No newline at end of file diff --git a/lib/YaoPlots/examples/circuits.jl b/lib/YaoPlots/examples/circuits.jl new file mode 100644 index 000000000..eb64245e6 --- /dev/null +++ b/lib/YaoPlots/examples/circuits.jl @@ -0,0 +1,50 @@ +using Yao.EasyBuild, YaoPlots, Yao + +lighttheme!() +YaoPlots.CircuitStyles.gate_bgcolor[] = "white" # default is transparent +# qft circuit +vizcircuit(qft_circuit(5); filename=joinpath(@__DIR__, "qft.png")) + +# labeled and time evolution +vizcircuit(chain(control(5, 3, (2,4)=>matblock(rand_unitary(4); tag="label")), + put(5, (2,4)=>matblock(rand_unitary(4); tag="labellabel")), time_evolve(+(put(5, 2=>X), kron(5, 2=>Z, 3=>Y)), 0.2)), filename=joinpath(@__DIR__, "labelled.png")) + +# variational circuit +vizcircuit(variational_circuit(5), filename=joinpath(@__DIR__, "variational.png")) +# vizcircuit(variational_circuit(5; mode=:Merged)) + +# general U4 gate +vizcircuit(general_U4(), filename=joinpath(@__DIR__, "u4.png")) + +# quantum supremacy circuit +vizcircuit(rand_supremacy2d(2, 2, 8), filename=joinpath(@__DIR__, "supremacy2d.png")) + +# google 52 qubit +vizcircuit(rand_google53(5); filename=joinpath(@__DIR__, "google53.png")) + +# control blocks +vizcircuit(chain(control(5, (2,-3), 4=>X), control(5, (-4, -2), 1=>Z)), filename=joinpath(@__DIR__, "controls.png")) + +# controlled kron +vizcircuit(control(4, 2, (1, 3)=>kron(X, X)), filename=joinpath(@__DIR__, "cxx.png")) + +vizcircuit(control(4, -collect(1:4-1), 4=>-Z), filename=joinpath(@__DIR__, "reflect.png")) + +vizcircuit(chain(5, [put(5, 2=>X), Yao.Measure(5; locs=(2,3)), Yao.Measure(5;locs=(2,)), Yao.Measure(5; resetto=bit"00110")]), filename=joinpath(@__DIR__, "measure.png")) + +vizcircuit(chain(5, [put(5, 2=>ConstGate.Sdag), put(5, 3=>ConstGate.Tdag), + put(5, (2,3)=>ConstGate.CNOT), put(5, (1,4)=>ConstGate.CZ), put(5, (1,2,5)=>ConstGate.Toffoli), + put(5, (2,3)=>ConstGate.SWAP), put(5, (1,)=>ConstGate.P0), put(5, (1,)=>ConstGate.I2), + put(5, (2,)=>ConstGate.P1), put(5, (1,)=>ConstGate.Pu), put(5, (1,)=>ConstGate.Pd), + put(5, (2,)=>ConstGate.T), + put(5, (2,)=>phase(0.4π)), + put(5, (2,)=>shift(0.4π)), + ]), filename=joinpath(@__DIR__, "constgates.png")) + +vizcircuit(chain(5, [put(5, (2,3)=>matblock(Matrix(SWAP), tag="SWAP")'), put(5, 2=>matblock(mat(I2), tag="id")), put(5, 2=>addlabel(X, "Xr")), control(5, (5,3), (2,4,1)=>put(3, (1,3)=>addlabel(SWAP, "SWAP")))]), filename=joinpath(@__DIR__, "multiqubit.png")) + +YaoPlots.darktheme!() +YaoPlots.CircuitStyles.gate_bgcolor[] = "transparent" # default is transparent +# qft circuit +vizcircuit(qft_circuit(5), filename=joinpath(@__DIR__, "qft-white.png")) + diff --git a/lib/YaoPlots/examples/rotation_gates.jl b/lib/YaoPlots/examples/rotation_gates.jl new file mode 100644 index 000000000..acb6eb198 --- /dev/null +++ b/lib/YaoPlots/examples/rotation_gates.jl @@ -0,0 +1,7 @@ +using Yao, YaoPlots, Yao.EasyBuild +N = 4 +d = 1 +ising(nbit, i, j) = rot(kron(nbit, i=>X, j=>X), 0.0) +circuit = dispatch!(variational_circuit(N, d, pair_ring(N), entangler=ising), :zero) + +circuit |> vizcircuit(scale=0.7) diff --git a/lib/YaoPlots/src/YaoPlots.jl b/lib/YaoPlots/src/YaoPlots.jl new file mode 100644 index 000000000..baa890904 --- /dev/null +++ b/lib/YaoPlots/src/YaoPlots.jl @@ -0,0 +1,24 @@ +module YaoPlots + +using YaoBlocks +using YaoBlocks.DocStringExtensions +using YaoArrayRegister +import Luxor +import Thebes +using Luxor: @layer, Point +using Thebes: Point3D, project +using LinearAlgebra: tr + +export CircuitStyles, vizcircuit, darktheme!, lighttheme! +export bloch_sphere, BlochStyles +export plot + +"""An alias of `vizcircuit`""" +plot(;kwargs...) = x->plot(x;kwargs...) +plot(blk::AbstractBlock; kwargs...) = vizcircuit(blk; kwargs...) + +include("helperblock.jl") +include("vizcircuit.jl") +include("bloch.jl") + +end diff --git a/lib/YaoPlots/src/bloch.jl b/lib/YaoPlots/src/bloch.jl new file mode 100644 index 000000000..a681bf7ec --- /dev/null +++ b/lib/YaoPlots/src/bloch.jl @@ -0,0 +1,288 @@ +""" + BlochStyles + +The module to define the default styles for bloch sphere drawing. +To change the default styles, you can modify the values in this module, e.g. +```julia +using YaoPlots +YaoPlots.BlochStyles.lw[] = 2.0 +``` + +### Style variables +#### Generic +- `lw`: the line width of the drawing +- `textsize`: the size of the text +- `fontfamily`: the font family of the text +- `background_color`: the background color of the drawing +- `color`: the color of the drawing + +#### Sphere +- `ball_size`: the size of the ball +- `dot_size`: the size of the dot +- `eye_point`: the eye point of the drawing + +#### Axis +- `axes_lw`: the line width of the axes +- `axes_colors`: the colors of the axes +- `axes_texts`: the texts of the axes, default to `["x", "y", "z"]` + +#### State display +- `show_projection_lines`: whether to show the projection lines +- `show_angle_texts`: whether to show the angle texts +- `show_line`: whether to show the line +- `show01`: whether to show the 0 and 1 states +""" +module BlochStyles + using Luxor + # generic config + const lw = Ref(1.0) + const textsize = Ref(16.0) + const fontfamily = Ref("monospace") + const background_color = Ref("transparent") + const color = Ref("#000000") + + # bloch sphere config + const ball_size = Ref(100) + const dot_size = Ref(3) + const eye_point = Ref((500, 200, 200)) + + # axes config + const axes_lw = Ref(1.0) + const axes_colors = ["#000000", "#000000", "#000000"] + const axes_texts = ["x", "y", "z"] + + # state display config + const show_projection_lines = Ref(false) + const show_angle_texts = Ref(false) + const show_line = Ref(true) + const show01 = Ref(false) +end + +""" +$TYPEDSIGNATURES + +Draw a bloch sphere, with the inputs being a list of `string => state` pairs, +where the string is a label for the state and a state can be a complex vector of size 2, a Yao register or `DensityMatrix`. +If you want to get a raw drawing, use `draw_bloch_sphere` instead. + +### Keyword Arguments +Note: The default values can be specified in submodule `BlochStyles`. + +- `textsize`: the size of the text +- `color`: the color of the drawing +- `drawing_size`: the size of the drawing +- `offset_x`: the offset of the drawing in x direction +- `offset_y`: the offset of the drawing in y direction +- `filename`: the filename of the output file, if not specified, a temporary file will be used +- `format`: the format of the output file, if not specified, the format will be inferred from the filename +- `fontfamily`: the font family of the text +- `background_color`: the background color of the drawing +- `lw`: the line width of the drawing +- `eye_point`: the eye point of the drawing +- `extra_kwargs`: extra keyword arguments passed to `draw_bloch_sphere` + - `dot_size`: the size of the dot + - `ball_size`: the size of the ball + - `show_projection_lines`: whether to show the projection lines + - `show_angle_texts`: whether to show the angle texts + - `show_line`: whether to show the line + - `show01`: whether to show the 0 and 1 states + - `colors`: the colors of the states + - `axes_lw`: the line width of the axes + - `axes_textsize`: the size of the axes texts + - `axes_colors`: the colors of the axes + - `axes_texts`: the texts of the axes + +### Examples + +```jldoctest +julia> using YaoPlots, YaoArrayRegister + +julia> bloch_sphere("|ψ⟩"=>rand_state(1), "ρ"=>density_matrix(rand_state(2), 1)); +``` +""" +function bloch_sphere(states...; + textsize=BlochStyles.textsize[], + color = BlochStyles.color[], + drawing_size = 300, + offset_x = 0, + offset_y = 0, + filename = nothing, + format = :svg, + fontfamily = BlochStyles.fontfamily[], + background_color = BlochStyles.background_color[], + lw = BlochStyles.lw[], + eye_point = BlochStyles.eye_point[], + extra_kwargs...) + + # file format + if filename === nothing + if format == :pdf + _format = tempname()*".pdf" + else + _format = format + end + else + _format = filename + end + # Set up the drawing canvas + Luxor.Drawing(drawing_size, drawing_size, _format) + Luxor.origin(drawing_size/2 + offset_x, drawing_size/2 + offset_y) + Luxor.background(background_color) + Luxor.sethue(color) + Luxor.fontsize(textsize) + fontfamily !== nothing && Luxor.fontface(fontfamily) + Luxor.setline(lw) + Thebes.eyepoint(eye_point...) + + draw_bloch_sphere(states...; eye_point, extra_kwargs...) + + # Save the drawing to a file + Luxor.finish() + Luxor.preview() +end + +# draw bloch sphere at the origin +function draw_bloch_sphere(states::Pair{<:AbstractString}...; + dot_size=BlochStyles.dot_size[], + ball_size=BlochStyles.ball_size[], + eye_point=BlochStyles.eye_point[], + show_projection_lines = BlochStyles.show_projection_lines[], + show_angle_texts = BlochStyles.show_angle_texts[], + show_line = BlochStyles.show_line[], + show01 = BlochStyles.show01[], + colors = fill(BlochStyles.color[], length(states)), + extra_kwargs... + ) + # get coordinate of a state + getcoo(x) = Point3D(ball_size .* state_to_cartesian(x)) + + # ball + Luxor.circle(Point(0, 0), ball_size, :stroke) + + # equator + nstep = 100 + equator_points = map(LinRange(0, 2π*(1-1/nstep), nstep)) do ϕ + project(Point3D(ball_size .* polar_to_cartesian(1.0, π/2, ϕ))) + end + Luxor.line.(equator_points[1:2:end], equator_points[2:2:end], :stroke) + + # show axes + axes3D(ball_size*3 ÷ 2; extra_kwargs...) + + # show 01 states + if show01 + for (txt, point) in [("|0⟩", [1, 0.0im]), ("|1⟩", [0.0im, 1])] + p = getcoo(point) + @layer begin + Luxor.sethue(BlochStyles.color[]) + if Thebes.distance(Point3D(0, 0, 0), Point3D(eye_point...)) < Thebes.distance(p, Point3D(eye_point...)) + Luxor.setopacity(0.3) + end + show_point(txt, project(p); dot_size, text_offset=Point(10, 0), show_line=false) + end + end + end + + # show points + for ((txt, point), color) in zip(states, colors) + p = getcoo(point) + @layer begin + Luxor.sethue(color) + if Thebes.distance(Point3D(0, 0, 0), Point3D(eye_point...)) < Thebes.distance(p, Point3D(eye_point...)) + Luxor.setopacity(0.3) + end + show_point(txt, project(p); dot_size, text_offset=Point(10, 0), show_line=show_line) + end + if show_projection_lines + # show θ + ratio = 0.2 + sz = project(Point3D(0, 0, ball_size*ratio)) + if show_angle_texts + Luxor.move(sz) + Luxor.arc2r(Point(0, 0), sz, project(p) * ratio, :stroke) + Luxor.text("θ", sz - Point(0, ball_size*0.07)) + end + # show equator projection and ϕ + equatorp = Point3D(p[1], p[2], 0) + sx = project(Point3D(ball_size*ratio, 0, 0)) + + if show_angle_texts + Luxor.move(sx) + Luxor.carc2r(Point(0, 0), sx, project(equatorp) * ratio, :stroke) + Luxor.text("ϕ", sx - Point(ball_size*0.12, 0)) + end + + @layer begin + Luxor.setdash("dot") + Luxor.setline(1) + Luxor.line(project(p), project(equatorp), :stroke) + Luxor.line(Point(0, 0), project(equatorp), :stroke) + end + end + end +end + +function show_point(txt, p; dot_size, text_offset, show_line) + Luxor.circle(p, dot_size, :fill) + Luxor.text(txt, p + text_offset) + show_line && Luxor.line(Point(0, 0), p, :stroke) +end + +function polar_to_cartesian(r, θ, ϕ) + x = r * sin(θ) * cos(ϕ) + y = r * sin(θ) * sin(ϕ) + z = r * cos(θ) + return x, y, z +end + +function cartesian_to_polar(x, y, z) + r = sqrt(x^2 + y^2 + z^2) + θ = acos(z/r) + ϕ = atan(y, x) + return r, θ, ϕ +end + +function state_to_polar(state::AbstractVector{Complex{T}}) where T + @assert length(state) == 2 + r = norm(state) + ϕ = iszero(state[1]) ? zero(T) : angle(state[2]/state[1]) + θ = 2 * atan(abs(state[2]), abs(state[1])) + return r, θ, ϕ +end +state_to_cartesian(state) = polar_to_cartesian(state_to_polar(state)...) + +# Draw labelled 3D axes with length `n`. +function axes3D(n::Int; + axes_lw = BlochStyles.axes_lw[], + axes_textsize = BlochStyles.textsize[], + axes_colors = BlochStyles.axes_colors, + axes_texts = BlochStyles.axes_texts, + ) + @layer begin + Luxor.fontsize(axes_textsize) + Luxor.setline(axes_lw) + for i = 1:3 + axis1 = project(Point3D(0.1, 0.1, 0.1)) + axis2 = [0.1, 0.1, 0.1] + axis2[i] = n + axis2 = project(Point3D(axis2...)) + Luxor.sethue(axes_colors[i]) + if (axis1 !== nothing) && (axis2 !== nothing) && !isapprox(axis1, axis2) + Luxor.arrow(axis1, axis2) + Luxor.label(axes_texts[i], :N, axis2, offset=10) + end + end + end +end + +# Interace to Yao +function state_to_polar(reg::AbstractRegister{2}) + @assert nqudits(reg) == 1 "Only single qubit register is allowed to plot on bloch sphere. If you want to plot a subsystem as a mixed state, please construct a density matrix first with `density_matrix`." + @assert nbatch(reg) == 1 || nbatch(reg) == YaoBlocks.NoBatch() "Only single batch register is allowed to plot on bloch sphere." + return state_to_polar(statevec(reg)) +end + +function state_to_cartesian(reg::DensityMatrix{2}) + @assert nqudits(reg) == 1 "Only single qubit density matrix is allowed to plot on bloch sphere" + return real(tr(reg.state * mat(X))), real(tr(reg.state * mat(Y))), real(tr(reg.state * mat(Z))) +end diff --git a/lib/YaoPlots/src/helperblock.jl b/lib/YaoPlots/src/helperblock.jl new file mode 100644 index 000000000..61aaa9418 --- /dev/null +++ b/lib/YaoPlots/src/helperblock.jl @@ -0,0 +1,50 @@ +using YaoBlocks +export LabelBlock + +""" + LabelBlock{BT,D} <: TagBlock{BT,D} + +A marker to mark a circuit applying on a continous block for better plotting. +""" +struct LabelBlock{BT<:AbstractBlock,D} <: TagBlock{BT,D} + content::BT + name::String +end + +YaoBlocks.content(cb::LabelBlock) = cb.content +function LabelBlock(x::BT, name::String) where {D,BT<:AbstractBlock{D}} + LabelBlock{BT,D}(x, name) +end + +function is_continuous_chunk(x) + length(x) == 0 && return true + return length(x) == maximum(x)-minimum(x)+1 +end + +YaoBlocks.PropertyTrait(::LabelBlock) = YaoBlocks.PreserveAll() +YaoBlocks.mat(::Type{T}, blk::LabelBlock) where {T} = mat(T, content(blk)) +YaoBlocks.unsafe_apply!(reg::YaoBlocks.AbstractRegister, blk::LabelBlock) = YaoBlocks.unsafe_apply!(reg, content(blk)) +YaoBlocks.chsubblocks(blk::LabelBlock, target::AbstractBlock) = LabelBlock(target, blk.name) + +Base.adjoint(x::LabelBlock) = LabelBlock(adjoint(content(x)), endswith(x.name, "†") ? x.name[1:end-1] : x.name*"†") +Base.copy(x::LabelBlock) = LabelBlock(copy(content(x)), x.name) +YaoBlocks.Optimise.to_basictypes(block::LabelBlock) = block + +export addlabel +addlabel(b::AbstractBlock, str::String) = LabelBlock(b, str) + +# to fix issue +function YaoBlocks.print_tree( + io::IO, + root::AbstractBlock, + node::LabelBlock, + depth::Int = 1, + islast::Bool = false, + active_levels = (); + maxdepth = 5, + charset = YaoBlocks.BlockTreeCharSet(), + title = true, + compact = false, +) + print(io, node.name) +end diff --git a/lib/YaoPlots/src/vizcircuit.jl b/lib/YaoPlots/src/vizcircuit.jl new file mode 100644 index 000000000..d4d28cb90 --- /dev/null +++ b/lib/YaoPlots/src/vizcircuit.jl @@ -0,0 +1,545 @@ +""" + CircuitStyles + +A module to define the styles of the circuit visualization. +To change the styles, please modify the variables in this module, e.g. +```julia +julia> using YaoPlots + +julia> YaoPlots.CircuitStyles.unit[] = 40 +40 +``` + +### Style variables +#### Sizes +* `unit` is the number of pixels in a unit. +* `r` is the size of nodes. +* `lw` is the line width. + +#### Texts +* `textsize` is the text size. +* `paramtextsize` is the text size for longer texts. +* `fontfamily` is the font family. + +#### Colors +* `linecolor` is the line color. +* `gate_bgcolor` is the gate background color. +* `textcolor` is the text color. +""" +module CircuitStyles + using Luxor + const unit = Ref(60) # number of pixels in a unit + const barrier_for_chain = Ref(false) + const r = Ref(0.2) + const lw = Ref(1.0) + const textsize = Ref(16.0) + const paramtextsize = Ref(10.0) + const fontfamily = Ref("monospace") + #const fontfamily = Ref("Dejavu Sans") + const linecolor = Ref("#000000") + const gate_bgcolor = Ref("transparent") + const textcolor = Ref("#000000") + + abstract type Gadget end + struct Box{FT} <: Gadget + height::FT + width::FT + end + struct Cross <: Gadget end + struct Dot <: Gadget end + struct NDot <: Gadget end + struct OPlus <: Gadget end + struct MeasureBox <: Gadget end + struct Phase <: Gadget + text::String + end + struct Text{FT} <: Gadget + fontsize::FT + end + struct Line <: Gadget end + get_width(::Cross) = 0.0 + get_width(::Phase) = r[]/2.5 + get_width(::Dot) = r[]/2.5 + get_width(::NDot) = r[]/2.5 + get_width(::OPlus) = r[]*1.4 + boxsize(b::Gadget) = (w = get_width(b); (w, w)) + function boxsize(b::Box) + return b.width, b.height + end + function boxsize(::MeasureBox) + return 2 * r[], 2 * r[] + end + + function render(b::Box, loc) + setcolor(gate_bgcolor[]) + Luxor.box(Point(loc)*unit[], b.width*unit[], b.height*unit[], :fill) + setcolor(linecolor[]) + setline(lw[]) + Luxor.box(Point(loc)*unit[], b.width*unit[], b.height*unit[], :stroke) + end + + function render(d::Dot, loc) + setcolor(linecolor[]) + circle(Point(loc)*unit[], get_width(d)*unit[]/2, :fill) + end + function render(d::Phase, loc) + x0 = Point(loc)*unit[] + setcolor(linecolor[]) + circle(x0, get_width(d)*unit[]/2, :fill) + setcolor(textcolor[]) + fontsize(textsize[]) + fontface(fontfamily[]) + text(d.text, x0+Point(8,8); valign=:middle, halign=:center) + end + function render(d::NDot, loc) + setcolor(gate_bgcolor[]) + circle(Point(loc)*unit[], get_width(d)*unit[]/2, :fill) + setcolor(linecolor[]) + setline(lw[]) + circle(Point(loc)*unit[], get_width(d)*unit[]/2, :stroke) + end + function render(::Cross, loc) + setline(lw[]) + setcolor(linecolor[]) + line(Point(loc[1]-r[]/sqrt(2), loc[2]-r[]/sqrt(2))*unit[], Point(loc[1]+r[]/sqrt(2), loc[2]+r[]/sqrt(2))*unit[], :stroke) + line(Point(loc[1]-r[]/sqrt(2), loc[2]+r[]/sqrt(2))*unit[], Point(loc[1]+r[]/sqrt(2), loc[2]-r[]/sqrt(2))*unit[], :stroke) + end + function render(d::OPlus, loc) + w = get_width(d) / 2 + setcolor(gate_bgcolor[]) + circle(Point(loc)*unit[], w*unit[], :fill) + setcolor(linecolor[]) + setline(lw[]) + x0 = Point(loc)*unit[] + circle(x0, w*unit[], :stroke) + line(x0+Point(-w*unit[], 0.0), x0+Point(w*unit[], 0.0*unit[]), :stroke) + line(x0+Point(0.0, -w*unit[]), x0+Point(0.0*unit[], w*unit[]), :stroke) + end + function render(::MeasureBox, loc) + x0 = Point(loc)*unit[] + setcolor(gate_bgcolor[]) + box(x0, 2*r[]*unit[], 2*r[]*unit[], :fill) + setcolor(linecolor[]) + setline(lw[]) + box(x0, 2*r[]*unit[], 2*r[]*unit[], :stroke) + move(x0+Point(-0.8*r[], 0.5*r[])*unit[]) + curve( + x0+Point(-0.8*r[], -0.6*r[])*unit[], + x0+Point(0.8*r[], -0.6*r[])*unit[], + x0+Point(0.8*r[], 0.5*r[])*unit[]) + strokepath() + line(x0+Point(0.0, 0.5*r[])*unit[], x0+Point(0.7*r[], -0.4*r[])*unit[], :stroke) + end + + function render(::Line, locs) + setcolor(linecolor[]) + setline(lw[]) + line(Point(locs[1])*unit[], Point(locs[2])*unit[], :stroke) + end + function render(t::Text, loctxt) + loc, txt, width, height = loctxt + fontsize(t.fontsize) + fontface(fontfamily[]) + #fontface("Dejavu Sans") + setcolor(textcolor[]) + if contains(txt, '\n') + for (i, txt) in enumerate(split(loctxt[2], "\n")) + text(txt, Point(loc)*unit[]+i*Point(0, t.fontsize)-Point((width-0.1)*unit[]/2, height*unit[]/2); halign=:left, valign=:middle) + end + else + text(txt, Point(loc)*unit[]; halign=:center, valign=:middle) + end + end + + Base.@kwdef struct GateStyles + g = Box(2*r[], 2*r[]) + c = Dot() + x = Cross() + nc = NDot() + not = OPlus() + measure = MeasureBox() + + # other styles + line = Line() + text = Text(textsize[]) + smalltext = Text(paramtextsize[]) + end +end + +struct CircuitGrid + frontier::Vector{Int} + w_depth::Float64 + w_line::Float64 + gatestyles::CircuitStyles.GateStyles +end + +nline(c::CircuitGrid) = length(c.frontier) +depth(c::CircuitGrid) = frontier(c, 1, nline(c)) +Base.getindex(c::CircuitGrid, i, j) = (c.w_depth*i, c.w_line*j) +Base.typed_vcat(c::CircuitGrid, ij1, ij2) = (c[ij1...], c[ij2...]) + +function CircuitGrid(nline::Int; w_depth=1.0, w_line=1.0, gatestyles=CircuitStyles.GateStyles()) + CircuitGrid(zeros(Int, nline), w_depth, w_line, gatestyles) +end + +function frontier(c::CircuitGrid, args...) + maximum(i->c.frontier[i], min(args..., nline(c)):max(args..., 1)) +end + +function _draw!(c::CircuitGrid, loc_brush_texts) + isempty(loc_brush_texts) && return + # a loc can be a integer, or a range + loc_brush_texts = sort(loc_brush_texts, by=x->maximum(x[1])) + locs = Iterators.flatten(getindex.(loc_brush_texts, 1)) |> collect + + # get the gate width and the circuit depth to draw + boxwidths, boxheights = Float64[], Float64[] + loc_brush_texts = map(loc_brush_texts) do (j, b, txt) + if length(j) != 0 + wspace, _ = text_width_and_size(txt) + hspace = (maximum(j)-minimum(j)) * c.w_line + # make box larger + b = b isa CircuitStyles.Box ? CircuitStyles.Box(b.height + hspace, b.width + wspace) : b + boxwidth, boxheight = CircuitStyles.boxsize(b) + push!(boxwidths, boxwidth) + push!(boxheights, boxheight) + end + (j, b, txt) + end + max_width = maximum(boxwidths) + ncolumn = max(1, ceil(Int, max_width/c.w_depth + 0.2)) # 0.1 is the minimum gap between two columns + + ipre = frontier(c, minimum(locs):maximum(locs)...) + i = ipre + ncolumn/2 + local jpre + for (k, ((j, b, txt), boxheight)) in enumerate(zip(loc_brush_texts, boxheights)) + length(j) == 0 && continue + _, fontsize = text_width_and_size(txt) + jmid = (minimum(j)+maximum(j))/2 + CircuitStyles.render(b, c[i, jmid]) + CircuitStyles.render(CircuitStyles.Text(fontsize), (c[i,jmid], txt, CircuitStyles.boxsize(b)...)) + # use line to connect blocks in the same gate + if k!=1 + CircuitStyles.render(c.gatestyles.line, c[(i, jmid-boxheight/2/c.w_line); (i, jpre)]) + end + jpre = jmid + boxheight/2/c.w_line + end + + #jmin, jmax = min(locs..., nline(c)), max(locs..., 1) + # connect horizontal lines + for (width, (j, b, txt)) in zip(boxwidths, loc_brush_texts) + for jj in j + CircuitStyles.render(c.gatestyles.line, c[(c.frontier[jj], jj); (i-width/2/c.w_depth, jj)]) + CircuitStyles.render(c.gatestyles.line, c[(i+width/2/c.w_depth, jj); (ipre+ncolumn, jj)]) + c.frontier[jj] = ipre + ncolumn + end + end + for j in setdiff(minimum(locs):maximum(locs), locs) + CircuitStyles.render(c.gatestyles.line, c[(c.frontier[j], j); (ipre+ncolumn, j)]) + c.frontier[j] = ipre + ncolumn + end +end + +function text_width_and_size(text) + lines = split(text, "\n") + widths = map(x->textwidth(x), lines) + (mw, loc) = findmax(widths) + fontsize = mw > 3 ? CircuitStyles.paramtextsize[] : CircuitStyles.textsize[] + # -2 because the gate has a default size + #width = max(W - 4, 0) * fontsize * 0.016 # mm to cm + Luxor.fontsize(fontsize) + Luxor.fontface(CircuitStyles.fontfamily[]) + width, height = Luxor.textextents(lines[loc])[3:4] + w = max(width / CircuitStyles.unit[] - CircuitStyles.r[]*2 + 0.2, 0.0) + return w, fontsize +end + +function initialize!(c::CircuitGrid; starting_texts, starting_offset) + starting_texts !== nothing && for j=1:nline(c) + CircuitStyles.render(c.gatestyles.text, (c[starting_offset, j], string(starting_texts[j]), c.w_depth, c.w_line)) + end +end + +function finalize!(c::CircuitGrid; show_ending_bar, ending_offset, ending_texts) + i = frontier(c, 1, nline(c)) + for j=1:nline(c) + show_ending_bar && CircuitStyles.render(c.gatestyles.line, c[(i, j-0.2); (i, j+0.2)]) + CircuitStyles.render(c.gatestyles.line, c[(i, j); (c.frontier[j], j)]) + ending_texts !== nothing && CircuitStyles.render(c.gatestyles.text, (c[i+ending_offset, j], string(ending_texts[j]), c.w_depth, c.w_line)) + end +end + +# elementary +function draw!(c::CircuitGrid, p::PrimitiveBlock, address, controls) + bts = length(controls)>=1 ? get_cbrush_texts(c, p) : get_brush_texts(c, p) + _draw!(c, [controls..., (getindex.(Ref(address), occupied_locs(p)), bts[1], bts[2])]) +end +function draw!(c::CircuitGrid, p::Daggered{<:PrimitiveBlock}, address, controls) + bts = length(controls)>=1 ? get_cbrush_texts(c, content(p)) : get_brush_texts(c, content(p)) + _draw!(c, [controls..., (getindex.(Ref(address), occupied_locs(p)), bts[1], bts[2]*"'")]) +end + +function draw!(c::CircuitGrid, p::Scale, address, controls) + fp = YaoBlocks.factor(p) + if !(abs(fp) ≈ 1) + error("can not visualize non-phase factor.") + end + draw!(c, YaoBlocks.phase(angle(fp)), [first(address)], controls) + draw!(c, p.content, address, controls) +end +# Special primitive gates +function draw!(c::CircuitGrid, ::I2Gate, address, controls) + return +end +function draw!(c::CircuitGrid, ::IdentityGate, address, controls) + return +end +function draw!(c::CircuitGrid, p::ConstGate.SWAPGate, address, controls) + bts = [(c.gatestyles.x, ""), (c.gatestyles.x, "")] + _draw!(c, [controls..., [(address[l], bt...) for (l, bt) in zip(occupied_locs(p), bts)]...]) +end +function draw!(c::CircuitGrid, p::ConstGate.CNOTGate, address, controls) + bts = [(c.gatestyles.c, ""), (c.gatestyles.x, "")] + _draw!(c, [controls..., [(address[l], bt...) for (l, bt) in zip(occupied_locs(p), bts)]...]) +end +function draw!(c::CircuitGrid, p::ConstGate.CZGate, address, controls) + bts = [(c.gatestyles.c, ""), (c.gatestyles.c, "")] + _draw!(c, [controls..., [(address[l], bt...) for (l, bt) in zip(occupied_locs(p), bts)]...]) +end +function draw!(c::CircuitGrid, p::ConstGate.ToffoliGate, address, controls) + bts = [(c.gatestyles.c, ""), (c.gatestyles.c, ""), (c.gatestyles.x, "")] + _draw!(c, [controls..., [(address[l], bt...) for (l, bt) in zip(occupied_locs(p), bts)]...]) +end + +# composite +function draw!(c::CircuitGrid, p::ChainBlock, address, controls) + for block in subblocks(p) + draw!(c, block, address, controls) + CircuitStyles.barrier_for_chain[] && set_barrier!(c, Int[address..., controls...]) + end +end + +function set_barrier!(c::CircuitGrid, locs::AbstractVector{Int}) + front = maximum(c.frontier[locs]) + for loc in locs + if c.frontier[loc] < front + CircuitStyles.render(c.gatestyles.line, c[(c.frontier[loc], loc); (front, loc)]) + c.frontier[loc] = front + end + end +end + +function draw!(c::CircuitGrid, p::PutBlock, address, controls) + locs = [address[i] for i in p.locs] + draw!(c, p.content, locs, controls) +end + +function draw!(c::CircuitGrid, m::YaoBlocks.Measure, address, controls) + @assert length(controls) == 0 + if m.postprocess isa RemoveMeasured + error("can not visualize post-processing: `RemoveMeasured`.") + end + if !(m.operator isa ComputationalBasis) + error("can not visualize measure blocks for operators") + end + locs = m.locations isa AllLocs ? collect(1:nqudits(m)) : [address[i] for i in m.locations] + for (i, loc) in enumerate(locs) + _draw!(c, [(loc, c.gatestyles.measure, "")]) + if m.postprocess isa ResetTo + # read i-th bit value + val = Int(m.postprocess.x)>>(i-1) & 1 + _draw!(c, [(loc, c.gatestyles.g, val == 1 ? "P₁" : "P₀")]) + end + end +end + +function draw!(c::CircuitGrid, cb::ControlBlock{GT,C}, address, controls) where {GT,C} + ctrl_locs = [address[i] for i in cb.ctrl_locs] + locs = [address[i] for i in cb.locs] + mycontrols = [(loc, (bit == 1 ? c.gatestyles.c : c.gatestyles.nc), "") for (loc, bit)=zip(ctrl_locs, cb.ctrl_config)] + draw!(c, cb.content, locs, [controls..., mycontrols...]) +end + +for B in [:LabelBlock, :GeneralMatrixBlock, :Add] + @eval function draw!(c::CircuitGrid, cb::$B, address, controls) + length(address) == 0 && return + #is_continuous_chunk(address) || error("address not continuous in a block marked as continous.") + _draw!(c, [controls..., (address, c.gatestyles.g, string(cb))]) + end +end + +# [:KronBlock, :RepeatedBlock, :CachedBlock, :Subroutine, :(YaoBlocks.AD.NoParams)] +function draw!(c::CircuitGrid, p::CompositeBlock, address, controls) + barrier_style = CircuitStyles.barrier_for_chain[] + CircuitStyles.barrier_for_chain[] = false + draw!(c, YaoBlocks.Optimise.to_basictypes(p), address, controls) + CircuitStyles.barrier_for_chain[] = barrier_style +end +for (GATE, SYM) in [(:XGate, :Rx), (:YGate, :Ry), (:ZGate, :Rz)] + @eval get_brush_texts(c, b::RotationGate{D,T,<:$GATE}) where {D,T} = (c.gatestyles.g, "$($(SYM))($(pretty_angle(b.theta)))") +end + +pretty_angle(theta) = string(theta) +function pretty_angle(theta::AbstractFloat) + c = continued_fraction(theta/π, 10) + if c.den < 100 + res = if c.num == 1 + "π" + elseif c.num==0 + "0" + elseif c.num==-1 + "-π" + else + "$(c.num)π" + end + if c.den != 1 + res *= "/$(c.den)" + end + res + else + "$(round(theta; digits=2))" + end +end +""" + continued_fraction(ϕ, n::Int) -> Rational + +Obtain `s` and `r` from `ϕ` that satisfies `|s//r - ϕ| ≦ 1/2r²` +""" +function continued_fraction(fl, n::Int) + if n == 1 || abs(mod(fl, 1)) < 1e-10 + Rational(floor(Int, fl), 1) + else + floor(Int, fl) + 1//continued_fraction(1/mod(fl, 1), n-1) + end +end + +get_brush_texts(c, ::ConstGate.SdagGate) = (c.gatestyles.g, "S'") +get_brush_texts(c, ::ConstGate.TdagGate) = (c.gatestyles.g, "T'") +get_brush_texts(c, ::ConstGate.PuGate) = (c.gatestyles.g, "P+") +get_brush_texts(c, ::ConstGate.PdGate) = (c.gatestyles.g, "P-") +get_brush_texts(c, ::ConstGate.P0Gate) = (c.gatestyles.g, "P₀") +get_brush_texts(c, ::ConstGate.P1Gate) = (c.gatestyles.g, "P₁") +get_brush_texts(c, b::PrimitiveBlock) = (c.gatestyles.g, string(b)) +get_brush_texts(c, b::TimeEvolution) = (c.gatestyles.g, string(b)) +get_brush_texts(c, b::ShiftGate) = (c.gatestyles.g, "ϕ($(pretty_angle(b.theta)))") +get_brush_texts(c, b::PhaseGate) = (CircuitStyles.Phase("$(pretty_angle(b.theta))"), "") +function get_brush_texts(c, b::T) where T<:ConstantGate + namestr = string(T.name.name) + if endswith(namestr, "Gate") + namestr = namestr[1:end-4] + end + # Fix! + (c.gatestyles.g, namestr) +end + +get_cbrush_texts(c, b::PrimitiveBlock) = get_brush_texts(c, b) +get_cbrush_texts(c, ::XGate) = (c.gatestyles.not, "") +get_cbrush_texts(c, ::ZGate) = (c.gatestyles.c, "") + +# front end +""" + vizcircuit(circuit; w_depth=0.85, w_line=0.75, format=:svg, filename=nothing, + show_ending_bar=false, starting_texts=nothing, starting_offset=-0.3, + ending_texts=nothing, ending_offset=0.3) + +Visualize a `Yao` quantum circuit. + +### Keyword Arguments +* `w_depth` is the circuit column width. +* `w_line` is the circuit row width. +* `format` can be `:svg`, `:png` or `:pdf`. +* `filename` can be `"*.svg"`, `"*.png"`, `"*.pdf"` or nothing (not saving to a file). +* `starting_texts` and `ending_texts` are texts shown before and after the circuit. +* `starting_offset` and `end_offset` are offsets (real values) for starting and ending texts. +* `show_ending_bar` is a boolean switch to show ending bar. + +### Styles +To change the gates styles like colors and lines, please modify the constants in submodule `CircuitStyles`. +They are defined as: + +* CircuitStyles.unit = Ref(60) # number of points in a unit +* CircuitStyles.r = Ref(0.2) # size of nodes +* CircuitStyles.lw = Ref(1.0) # line width +* CircuitStyles.textsize = Ref(16.0) # text size +* CircuitStyles.paramtextsize = Ref(10.0) # text size (longer texts) +* CircuitStyles.fontfamily = Ref("monospace") # font family +* CircuitStyles.linecolor = Ref("#000000") # line color +* CircuitStyles.gate_bgcolor = Ref("transparent") # gate background color +* CircuitStyles.textcolor = Ref("#000000") # text color +""" +function vizcircuit(blk::AbstractBlock; w_depth=0.85, w_line=0.75, format=:svg, filename=nothing, + show_ending_bar=false, starting_texts=nothing, starting_offset=-0.3, + ending_texts=nothing, ending_offset=0.3, gatestyles=CircuitStyles.GateStyles()) + img = circuit_canvas(nqudits(blk); w_depth, w_line, show_ending_bar, starting_texts, starting_offset, + ending_texts, ending_offset, gatestyles, format, filename) do c + addblock!(c, blk) + end + return img +end + +addblock!(c::CircuitGrid, blk::AbstractBlock) = draw!(c, blk, collect(1:nqudits(blk)), []) +addblock!(c::CircuitGrid, blk::Function) = addblock!(c, blk(nline(c))) + +function circuit_canvas(f, nline::Int; format=:svg, filename=nothing, w_depth=0.85, w_line=0.75, + show_ending_bar=false, starting_texts=nothing, starting_offset=-0.3, ending_texts=nothing, + ending_offset=0.3, gatestyles=CircuitStyles.GateStyles()) + # the first time to estimate the canvas size + Luxor.Drawing(50, 50, :png) + c = CircuitGrid(nline; w_depth, w_line, gatestyles) + initialize!(c; starting_texts, starting_offset) + f(c) + finalize!(c; show_ending_bar, ending_texts, ending_offset) + Luxor.finish() + # the second time draw + u = CircuitStyles.unit[] + a, b = ceil(Int, (depth(c)+1)*w_depth*u), ceil(Int, nline*w_line*u) + _luxor(a, b, w_depth/2*u, -w_line/2*u; format, filename) do + c = CircuitGrid(nline; w_depth, w_line, gatestyles) + initialize!(c; starting_texts, starting_offset) + f(c) + finalize!(c; show_ending_bar, ending_texts, ending_offset) + end +end + +function _luxor(f, Dx, Dy, offsetx, offsety; format, filename) + if filename === nothing + if format == :pdf + _format = tempname()*".pdf" + else + _format = format + end + else + _format = filename + end + Luxor.Drawing(round(Int,Dx), round(Int,Dy), _format) + Luxor.origin(offsetx, offsety) + f() + Luxor.finish() + Luxor.preview() +end + +vizcircuit(; kwargs...) = c->vizcircuit(c; kwargs...) + +""" + darktheme!() + +Change the default theme to dark. +""" +function darktheme!() + const CircuitStyles.linecolor[] = "#FFFFFF" + const CircuitStyles.textcolor[] = "#FFFFFF" + const BlochStyles.color[] = "#FFFFFF" + BlochStyles.axes_colors .= ["#FFFFFF", "#FFFFFF", "#FFFFFF"] +end + +""" + lighttheme!() + +Change the default theme to light. +""" +function lighttheme!() + const CircuitStyles.linecolor[] = "#000000" + const CircuitStyles.textcolor[] = "#000000" + const BlochStyles.color[] = "#000000" + BlochStyles.axes_colors .= ["#000000", "#000000", "#000000"] +end diff --git a/lib/YaoPlots/test/bloch.jl b/lib/YaoPlots/test/bloch.jl new file mode 100644 index 000000000..e26c641cc --- /dev/null +++ b/lib/YaoPlots/test/bloch.jl @@ -0,0 +1,43 @@ +using YaoPlots, Test +using LinearAlgebra: normalize!, eigen +using Luxor: Drawing +using YaoBlocks +using YaoArrayRegister + +@testset "spherical coo" begin + for i=1:10 + x = randn(3) + @test collect(YaoPlots.polar_to_cartesian(YaoPlots.cartesian_to_polar(x...)...)) ≈ x + end +end + +@testset "state to polar" begin + nx, ny, nz = normalize!(randn(3)) + sigma = nx * X + ny * Y + nz * Z + m = Matrix(sigma) + E, U = eigen(m) + @test collect(YaoPlots.state_to_cartesian(U[:, 2])) ≈ [nx, ny, nz] + @test collect(YaoPlots.state_to_cartesian(U[:, 1])) ≈ -[nx, ny, nz] + + reg = rand_state(1) + @test collect(YaoPlots.state_to_cartesian(reg)) ≈ collect(YaoPlots.state_to_cartesian(density_matrix(reg))) +end + +@testset "bloch" begin + @test bloch_sphere("|ψ⟩"=>[2.2, 0.3im+0.3]; show_projection_lines=true, + show_angle_texts=true, show_line=true, show01=true) isa Drawing + # dark theme + darktheme!() + @test bloch_sphere("|ψ⟩"=>[2.2, 0.3im+0.3]; show_projection_lines=true, + show_angle_texts=true, show_line=true, show01=true) isa Drawing +end + +@testset "draw reg" begin + @test bloch_sphere("|ψ⟩"=>rand_state(1)) isa Drawing +end + +@testset "draw density matrix" begin + rho = density_matrix(rand_state(2), 1) + @test bloch_sphere("ρ"=>rho) isa Drawing + bloch_sphere("|ψ⟩"=>rand_state(1), "ρ"=>density_matrix(rand_state(2), 1)) +end \ No newline at end of file diff --git a/lib/YaoPlots/test/helperblock.jl b/lib/YaoPlots/test/helperblock.jl new file mode 100644 index 000000000..85baaf701 --- /dev/null +++ b/lib/YaoPlots/test/helperblock.jl @@ -0,0 +1,29 @@ +using YaoPlots, YaoBlocks, Luxor +using YaoArrayRegister +using Test + +@testset "LabelBlock" begin + x = put(5, (2,3)=>matblock(rand_unitary(4))) + cb = LabelBlock(x, "x") + @test mat(copy(cb)) == mat(cb) + @test isunitary(cb) + @test ishermitian(cb) == ishermitian(x) + @test isreflexive(cb) == isreflexive(x) + @test mat(cb) == mat(x) + reg = rand_state(5) + @test apply!(copy(reg), cb) ≈ apply!(copy(reg), x) + @test cb' isa LabelBlock && mat(cb') ≈ mat(cb)' + @test (cb').name == "x†" && (cb'').name == "x" + + y = put(5, (3,4)=>matblock(rand_unitary(4))) + cc = chsubblocks(cb, y) + @test YaoPlots.is_continuous_chunk([1,2,3]) == true + @test YaoPlots.is_continuous_chunk([1,2,4]) == false + @test YaoPlots.is_continuous_chunk([3,2,4]) == true + + c1 = chain(5, [put(5, (2,3)=>addlabel(SWAP, "SWAP")), put(5, 2=>addlabel(I2, "id")), put(5, 2=>addlabel(X, "X")), control(5, (5,3), (2,4,1)=>put(3, (1,3)=>addlabel(SWAP, "SWAP")))]) + c2 = chain(5, [put(5, (2,3)=>addlabel(SWAP, "SWAP")), put(5, 2=>addlabel(I2, "id")), put(5, 2=>addlabel(X, "X")), control(5, (5,3), (2,4,1)=>put(3, (1,2)=>addlabel(SWAP, "SWAP")))]) + + @test vizcircuit(c1) isa Drawing + @test vizcircuit(c2) isa Drawing +end diff --git a/lib/YaoPlots/test/runtests.jl b/lib/YaoPlots/test/runtests.jl new file mode 100644 index 000000000..27a917f48 --- /dev/null +++ b/lib/YaoPlots/test/runtests.jl @@ -0,0 +1,14 @@ +using YaoPlots +using Test + +@testset "helperblock" begin + include("helperblock.jl") +end + +@testset "vizcircuit" begin + include("vizcircuit.jl") +end + +@testset "bloch" begin + include("bloch.jl") +end \ No newline at end of file diff --git a/lib/YaoPlots/test/vizcircuit.jl b/lib/YaoPlots/test/vizcircuit.jl new file mode 100644 index 000000000..bf063aee6 --- /dev/null +++ b/lib/YaoPlots/test/vizcircuit.jl @@ -0,0 +1,89 @@ +using YaoPlots +using Test +using YaoBlocks +using YaoPlots: addblock!, CircuitGrid, circuit_canvas +using YaoArrayRegister +using Luxor + +@testset "gate styles" begin + c = YaoPlots.CircuitGrid(1) + @test YaoPlots.get_brush_texts(c, X)[2] == "X" + @test YaoPlots.get_brush_texts(c, Rx(0.5))[2] == "Rx(0.5)" + @test YaoPlots.get_brush_texts(c, shift(0.5))[2] == "ϕ(0.5)" + @test YaoPlots.get_brush_texts(c, YaoBlocks.phase(0.5))[2] == "" +end + +@testset "circuit canvas" begin + c = CircuitGrid(5) + @test YaoPlots.nline(c) == 5 + @test YaoPlots.frontier(c, 2, 3) == 0 + @test YaoPlots.depth(c) == 0 + circuit_canvas(5) do c + YaoPlots.draw!(c, put(5, 3=>X), 1:5, []) + @test YaoPlots.frontier(c, 1, 2) == 0 + @test YaoPlots.frontier(c, 3, 5) == 1 + @test YaoPlots.depth(c) == 1 + end + + gg = circuit_canvas(5) do c + addblock!(c, put(3=>X)) + addblock!(c, put(3=>matblock(rand_unitary(2); tag="xxx"))) + addblock!(c, put(3=>igate(1))) + addblock!(c, time_evolve(put(3,3=>igate(1)), 0.1)) + addblock!(c, control(2, 3=>X)) + addblock!(c, control(2, 3=>shift(0.5))) + addblock!(c, chain(5, control(2, 3=>X), put(1=>X))) + @test YaoPlots.depth(c) >= 6 + end + @test gg isa Drawing + + g = put(5, (3, 4)=>SWAP) |> vizcircuit(; w_line=0.8, w_depth=0.9) + @test g isa Drawing + @test vizcircuit(put(5, (3,4)=>kron(X, Y)); w_line=0.8, w_depth=0.9) isa Drawing + + @test vizcircuit(control(10, (2, -3), 6=>X)) isa Drawing + @test vizcircuit(control(10, (2, -3), 6=>im*X)) isa Drawing + @test plot(put(7, (2,3)=>matblock(randn(4,4)))) isa Drawing + # qudit + @test vizcircuit(put(10, 6=>matblock(rand_unitary(3); nlevel=3, tag="3 level"))) isa Drawing +end + +@testset "fix #3" begin + @test (control(4, 2, (1,3)=>kron(X, X)) |> vizcircuit) isa Drawing +end + +@testset "pretty_angle" begin + @test YaoPlots.pretty_angle(π) == "π" + @test YaoPlots.pretty_angle(π*1.0) == "π" + @test YaoPlots.pretty_angle(-π*1.0) == "-π" + @test YaoPlots.pretty_angle(-π*0.0) == "0" + @test YaoPlots.pretty_angle(-π*0.5) == "-π/2" + @test YaoPlots.pretty_angle(1.411110) == "1.41" +end + +@testset "readme" begin + YaoPlots.CircuitStyles.linecolor[] = "pink" + YaoPlots.CircuitStyles.gate_bgcolor[] = "yellow" + YaoPlots.CircuitStyles.textcolor[] = "#000080" # the navy blue color + YaoPlots.CircuitStyles.fontfamily[] = "JuliaMono" + YaoPlots.CircuitStyles.lw[] = 2.5 + YaoPlots.CircuitStyles.textsize[] = 13 + YaoPlots.CircuitStyles.paramtextsize[] = 8 + + @test plot(chain(3, put(1=>X), repeat(3, H), put(2=>Y), repeat(3, Rx(π/2)))) isa Drawing + darktheme!() + @test plot(chain(3, put(1=>X), repeat(3, H), put(2=>Y), repeat(3, Rx(π/2)))) isa Drawing + lighttheme!() + @test plot(chain(3, put(1=>X), repeat(3, H), put(2=>Y), repeat(3, Rx(π/2)))) isa Drawing +end + +@testset "readme" begin + circuit = chain( + 4, + kron(X, H, H, H), + kron(1=>Y, 4=>H), + put(2=>Y), + ) + YaoPlots.CircuitStyles.barrier_for_chain[] = true + @test vizcircuit(circuit) isa Drawing +end \ No newline at end of file diff --git a/lib/YaoSym/Project.toml b/lib/YaoSym/Project.toml index 2d58194d9..96a87a533 100644 --- a/lib/YaoSym/Project.toml +++ b/lib/YaoSym/Project.toml @@ -1,28 +1,26 @@ name = "YaoSym" uuid = "3b27209a-d3d6-11e9-3c0f-41eb92b2cb9d" -version = "0.6.4" +version = "0.6.8" [deps] BitBasis = "50ba71b6-fa0f-514d-ae9a-0916efc90dcf" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LuxurySparse = "d05aeea4-b7d4-55ac-b691-9e7fabb07ba2" -Requires = "ae029012-a4dd-5104-9daa-d747884805df" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +SymEngine = "123dc426-2d89-5057-bbad-38513e3affd8" YaoArrayRegister = "e600142f-9330-5003-8abb-0ebd767abc51" YaoBlocks = "418bc28f-b43b-5e0b-a6e7-61bbc1a2c1df" [compat] -BitBasis = "0.8" +BitBasis = "0.8, 0.9" LuxurySparse = "0.7" -Requires = "1" -SymEngine = "0.8" +SymEngine = "0.11" YaoArrayRegister = "0.9" YaoBlocks = "0.13" julia = "1" [extras] -SymEngine = "123dc426-2d89-5057-bbad-38513e3affd8" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Test", "SymEngine"] +test = ["Test"] diff --git a/lib/YaoSym/src/YaoSym.jl b/lib/YaoSym/src/YaoSym.jl index 6d3af978b..53397e32f 100644 --- a/lib/YaoSym/src/YaoSym.jl +++ b/lib/YaoSym/src/YaoSym.jl @@ -1,11 +1,20 @@ module YaoSym -using Requires +using SparseArrays, LuxurySparse, LinearAlgebra +using BitBasis, YaoArrayRegister, YaoBlocks +import YaoArrayRegister: parametric_mat +using SymEngine +using SymEngine: @vars, Basic, BasicType, BasicOp, BasicTrigFunction, BasicComplexNumber -include("register.jl") +# SymEngine APIs +export @vars, Basic, subs, expand, simplify_expi + +# Symbolic registers +export @ket_str, @bra_str +export SymReg, AdjointSymReg, SymRegOrAdjointSymReg +export szero_state -function __init__() - @require SymEngine = "123dc426-2d89-5057-bbad-38513e3affd8" include("symengine/backend.jl") -end +include("register.jl") +include("symengine/symengine.jl") end # module diff --git a/lib/YaoSym/src/register.jl b/lib/YaoSym/src/register.jl index 1e80db8cb..87cab1d62 100644 --- a/lib/YaoSym/src/register.jl +++ b/lib/YaoSym/src/register.jl @@ -1,6 +1,3 @@ -using SparseArrays, BitBasis, YaoArrayRegister -export @ket_str, @bra_str - function parse_str(s::String) v = 0 k = 1 @@ -38,14 +35,14 @@ Create a ket register. See also [`@bra_str`](@ref). a symbolic quantum state can be created simply by -```jldoctest; setup=:(using Yao, SymEngine) +```jldoctest; setup=:(using Yao) julia> ket"110" + 2ket"111" |110⟩ + 2.0|111⟩ ``` qubits can be partially actived by [`focus!`](@ref) -```jldoctest; setup=:(using Yao, SymEngine) +```jldoctest; setup=:(using Yao) julia> ket"100" + ket"111" |> focus!(1:2) |100⟩ + |111⟩ ``` @@ -65,7 +62,7 @@ Create a bra register. See also [`@ket_str`](@ref). Similar to `@ket_str` literal, a symbolic quantum state can be created by -```jldoctest; setup=:(using Yao, SymEngine) +```jldoctest; setup=:(using Yao) julia> bra"111" + 2bra"101" 2.0⟨101| + ⟨111| diff --git a/lib/YaoSym/src/symengine/blocks.jl b/lib/YaoSym/src/symengine/blocks.jl index e2052f9df..93024ac3b 100644 --- a/lib/YaoSym/src/symengine/blocks.jl +++ b/lib/YaoSym/src/symengine/blocks.jl @@ -1,63 +1,37 @@ -using YaoBlocks -using LuxurySparse -using LinearAlgebra -using ..SymEngine -using ..SymEngine: BasicType, BasicOp, BasicTrigFunction - op_types = [:Mul, :Add, :Pow] const BiVarOp = Union{[SymEngine.BasicType{Val{i}} for i in op_types]...} -export @vars - simag = SymFunction("Im") sreal = SymFunction("Re") -sabs = SymFunction("abs") +sangle = SymFunction("angle") Base.promote_rule(::Type{Bool}, ::Type{Basic}) = Basic # NOTE: need to annotate the output because otherwise it is type unstable! -Base.conj(x::Basic)::Basic = Basic(conj(SymEngine.BasicType(x))) -Base.conj(x::BasicType) = real(x) - im * imag(x) +Base.abs(x::Basic)::Basic = Basic(abs(SymEngine.BasicType(x))) +function Base.abs(x::BasicType{Val{:Pow}}) + a, b = get_args(x.x) + abs(a)^real(b) * exp(-angle(a) * imag(b)) +end + +Base.angle(x::Basic)::Basic = Basic(angle(SymEngine.BasicType(x))) +Base.angle(x::BasicType)::Basic = sangle(x.x) +Base.angle(::BasicType{Val{:Symbol}})::Basic = Basic(0) +Base.angle(::BasicType{Val{:Constant}})::Basic = Basic(0) +Base.angle(::SymEngine.BasicRealNumber)::Basic = Basic(0) + Base.conj(x::BiVarOp) = juliafunc(x)(conj.(get_args(x.x))...) Base.conj(x::BasicTrigFunction) = juliafunc(x)(conj.(get_args(x.x)...)...) -# WARNING: symbols and constants are assumed real! -Base.imag(x::BasicType{Val{:Constant}}) = Basic(0) -Base.imag(x::BasicType{Val{:Symbol}}) = Basic(0) -Base.abs(x::Basic) = Basic(abs(SymEngine.BasicType(x))) -Base.abs(x::BasicType{Val{:Constant}}) = x -Base.abs(x::BasicType{Val{:Symbol}}) = x + +Base.imag(::BasicType{Val{:Constant}}) = Basic(0) +Base.imag(::BasicType{Val{:Symbol}}) = Basic(0) function Base.imag(x::BasicType{Val{:Add}}) args = get_args(x.x) mapreduce(imag, +, args) end - -function Base.real(x::BasicType{Val{:Add}}) - args = get_args(x.x) - mapreduce(real, +, args) -end - -function Base.abs(x::BasicType{Val{:Add}}) - args = get_args(x.x) - mapreduce(abs, +, args) -end - function Base.imag(x::BasicType{Val{:Mul}}) args = (get_args(x.x)...,) get_mul_imag(args) end - -function Base.real(x::BasicType{Val{:Pow}}) - a, b = get_args(x.x) - if imag(a) == 0 && imag(b) == 0 - return x.x - else - if imag(a) == 0 - return a^real(b) * cos(log(a) * imag(b)) - else - return sreal(x.x) - end - end -end - function Base.imag(x::BasicType{Val{:Pow}}) a, b = get_args(x.x) if imag(a) == 0 && imag(b) == 0 @@ -70,53 +44,51 @@ function Base.imag(x::BasicType{Val{:Pow}}) end end end - -function Base.abs(x::BasicType{Val{:Pow}}) - a, b = get_args(x.x) - abs(a)^real(b) -end - -function Base.real(x::BasicType{Val{:Mul}}) - args = (get_args(x.x)...,) - get_mul_real(args) +function Base.imag(x::BasicTrigFunction) + a, = get_args(x.x) + if imag(a) == 0 + return Basic(0) + else + return simag(x.x) + end end - function get_mul_imag(args::NTuple{N,Any}) where {N} imag(args[1]) * get_mul_real(args[2:end]) + real(args[1]) * get_mul_imag(args[2:end]) end get_mul_imag(args::Tuple{Basic}) = imag(args[1]) -function get_mul_real(args::NTuple{N,Any}) where {N} - real(args[1]) * get_mul_real(args[2:end]) - imag(args[1]) * get_mul_imag(args[2:end]) +function Base.real(x::BasicType{Val{:Add}}) + args = get_args(x.x) + mapreduce(real, +, args) end -get_mul_real(args::Tuple{Basic}) = real(args[1]) - -function Base.real(x::BasicTrigFunction) - a, = get_args(x.x) - if imag(a) == 0 +function Base.real(x::BasicType{Val{:Pow}}) + a, b = get_args(x.x) + if imag(a) == 0 && imag(b) == 0 return x.x else - return sreal(x.x) + if imag(a) == 0 + return a^real(b) * cos(log(a) * imag(b)) + else + return sreal(x.x) + end end end - -function Base.abs(x::BasicTrigFunction) +function Base.real(x::BasicType{Val{:Mul}}) + args = (get_args(x.x)...,) + get_mul_real(args) +end +function Base.real(x::BasicTrigFunction) a, = get_args(x.x) if imag(a) == 0 return x.x else - return sabs(x.x) + return sreal(x.x) end end - -function Base.imag(x::BasicTrigFunction) - a, = get_args(x.x) - if imag(a) == 0 - return Basic(0) - else - return simag(x.x) - end +function get_mul_real(args::NTuple{N,Any}) where {N} + real(args[1]) * get_mul_real(args[2:end]) - imag(args[1]) * get_mul_imag(args[2:end]) end +get_mul_real(args::Tuple{Basic}) = real(args[1]) @generated function juliafunc(x::BasicType{Val{T}}) where {T} SymEngine.map_fn(T, SymEngine.fn_map) @@ -132,7 +104,8 @@ YaoBlocks.shift(θ::SymReal) = ShiftGate(θ) YaoBlocks.mat(::Type{Basic}, gate::GT) where GT<:ConstantGate = _pretty_basic.(mat(gate)) YaoBlocks.mat(::Type{Basic}, gate::ConstGate.TGate) = Diagonal(Basic[1, exp(Basic(im)*Basic(π)/4)]) YaoBlocks.mat(::Type{Basic}, gate::ConstGate.TdagGate) = Diagonal(Basic[1, exp(-Basic(im)*Basic(π)/4)]) -YaoBlocks.mat(::Type{Basic}, ::HGate) = 1 / sqrt(Basic(2)) * Basic[1 1; 1 -1] +YaoArrayRegister._hadamard_matrix(::Type{Basic}) = 1 / sqrt(Basic(2)) * Basic[1 1; 1 -1] +YaoBlocks.mat(::Type{Basic}, ::HGate) = YaoArrayRegister._hadamard_matrix(Basic) YaoBlocks.mat(::Type{Basic}, gate::ShiftGate) = Diagonal([1, exp(im * gate.theta)]) YaoBlocks.mat(::Type{Basic}, gate::PhaseGate) = exp(im * gate.theta) * IMatrix(2) function YaoBlocks.mat(::Type{Basic}, R::RotationGate{D}) where {D} @@ -154,7 +127,6 @@ YaoBlocks.PSwap(n::Int, locs::Tuple{Int,Int}, θ::SymReal) = YaoBlocks.pswap(n::Int, i::Int, j::Int, α::SymReal) = PSwap(n, (i, j), α) YaoBlocks.pswap(i::Int, j::Int, α::SymReal) = n -> pswap(n, i, j, α) -export subs SymEngine.subs(c::AbstractBlock, args...; kwargs...) = subs(Basic, c, args...; kwargs...) function SymEngine.subs(::Type{T}, c::AbstractBlock, args...; kwargs...) where {T} c = setiparams(c, map(x -> T(subs(x, args...; kwargs...)), getiparams(c))...) diff --git a/lib/YaoSym/src/symengine/instruct.jl b/lib/YaoSym/src/symengine/instruct.jl index bae33ba4e..8fa89a345 100644 --- a/lib/YaoSym/src/symengine/instruct.jl +++ b/lib/YaoSym/src/symengine/instruct.jl @@ -1,6 +1,3 @@ -using ..SymEngine -import YaoArrayRegister: parametric_mat - parametric_mat(::Type{T}, ::Val{:Rx}, theta::Basic) where {T} = Basic[cos(theta / 2) -im*sin(theta / 2); -im*sin(theta / 2) cos(theta / 2)] parametric_mat(::Type{T}, ::Val{:Ry}, theta::Basic) where {T} = @@ -26,29 +23,29 @@ for G in [:Rx, :Ry, :Rz, :CPHASE] # forward single gates @eval function YaoArrayRegister.instruct!(::Val{2}, - state::AbstractVecOrMat{T}, + state::AbstractVecOrMat, g::Val{$(QuoteNode(G))}, locs::NTuple{N1,Int}, theta::Basic, - ) where {T,N1} + ) where {N1} instruct!(Val(2), state, g, locs, (), (), theta) end @eval function YaoArrayRegister.instruct!(::Val{2}, - state::AbstractVecOrMat{T}, + state::AbstractVecOrMat, g::Val{$(QuoteNode(G))}, locs::Tuple{Int}, theta::Basic, - ) where {T,N1} + ) instruct!(Val(2), state, g, locs, (), (), theta) end end @eval function YaoArrayRegister.instruct!(::Val{2}, - state::AbstractVecOrMat{T}, + state::AbstractVecOrMat, g::Val{:PSWAP}, locs::Tuple{Int,Int}, theta::Basic, -) where {T,N1} +) instruct!(Val(2), state, g, locs, (), (), theta) end diff --git a/lib/YaoSym/src/symengine/patch.jl b/lib/YaoSym/src/symengine/patch.jl index b22fa5aa9..470241f57 100644 --- a/lib/YaoSym/src/symengine/patch.jl +++ b/lib/YaoSym/src/symengine/patch.jl @@ -1,6 +1,3 @@ -using ..SymEngine: Basic -export simplify_expi - function Base.iszero(x::Basic) isempty(free_symbols(x)) && iszero(N(x)) end diff --git a/lib/YaoSym/src/symengine/register.jl b/lib/YaoSym/src/symengine/register.jl index 4599a9f5c..50ffcd79d 100644 --- a/lib/YaoSym/src/symengine/register.jl +++ b/lib/YaoSym/src/symengine/register.jl @@ -1,9 +1,3 @@ -using SparseArrays, BitBasis, YaoArrayRegister -using ..SymEngine -export @ket_str, @bra_str -export SymReg, AdjointSymReg, SymRegOrAdjointSymReg, expand -export szero_state - YaoArrayRegister._warn_type(raw::AbstractMatrix{Basic}) = nothing const SymReg{D,MT} = AbstractArrayReg{D,Basic,MT} where {MT<:AbstractMatrix{Basic}} @@ -34,7 +28,7 @@ function _pretty_basic(x::Complex) end end -function ket_m(s) +function ket_m(s::AbstractString) v, N = parse_str(s) st = spzeros(Basic, 1 << N, 1) st[v+1] = 1 diff --git a/lib/YaoSym/src/symengine/backend.jl b/lib/YaoSym/src/symengine/symengine.jl similarity index 53% rename from lib/YaoSym/src/symengine/backend.jl rename to lib/YaoSym/src/symengine/symengine.jl index 14c46ae01..407334c5a 100644 --- a/lib/YaoSym/src/symengine/backend.jl +++ b/lib/YaoSym/src/symengine/symengine.jl @@ -1,7 +1,3 @@ -using ..SymEngine -using ..SymEngine: @vars, Basic -export @vars, Basic, subs - include("register.jl") include("instruct.jl") include("blocks.jl") diff --git a/lib/YaoSym/test/symengine/blocks.jl b/lib/YaoSym/test/symengine/blocks.jl index 78bc4e022..37cc8bdb4 100644 --- a/lib/YaoSym/test/symengine/blocks.jl +++ b/lib/YaoSym/test/symengine/blocks.jl @@ -1,10 +1,10 @@ using YaoSym, YaoBlocks, YaoArrayRegister, SymEngine -using YaoSym: simplify_expi +using YaoSym: simplify_expi, BasicType using SymEngine using Test @testset "imag and conj" begin - @vars a b c x + @vars a b c d x @test imag(Basic(2 + im * a)) == a @test imag(Basic((2 + im) * (im * a + 1))) == 2 * a + 1 @test real(sin(a)) == sin(a) @@ -14,8 +14,9 @@ using Test @test real(exp(im * a)) == cos(a) @test imag(exp(im * a)) == sin(a) @test abs(exp(im * a)) == 1 - @test abs(sin(a)) == sin(a) - @test abs(sin(a) + 2) == sin(a) + 2 + @test BasicType(imag(cos(im * a + b))) isa BasicType{Val{:FunctionSymbol}} + @test BasicType(real((im*c + d)^(im * a + b))) isa BasicType{Val{:FunctionSymbol}} + @test BasicType(real(sin(im * a + b))) isa BasicType{Val{:FunctionSymbol}} end @testset "mat" begin @@ -124,3 +125,12 @@ end @test test_check_dumpload(time_evolve(X, θ)) @test test_check_dumpload(rot(X, θ)) end + +@testset "angle" begin + @vars x + @test_broken Float64(angle(Basic(1 + im))) ≈ π / 4 + @test angle(x) == Basic(0) + @test BasicType(angle((1+x)^x)) isa BasicType{Val{:FunctionSymbol}} + @test angle(Basic(1)) == Basic(0) + @test angle(exp(Basic(1))) == Basic(0) +end \ No newline at end of file diff --git a/lib/YaoToEinsum/CITATION.bib b/lib/YaoToEinsum/CITATION.bib new file mode 100644 index 000000000..42462489b --- /dev/null +++ b/lib/YaoToEinsum/CITATION.bib @@ -0,0 +1,61 @@ +@article{Luo2020, + title = {Yao.jl: {Extensible}, {Efficient} {Framework} for {Quantum} {Algorithm} {Design}}, + volume = {4}, + shorttitle = {Yao.jl}, + url = {https://quantum-journal.org/papers/q-2020-10-11-341/}, + doi = {10.22331/q-2020-10-11-341}, + abstract = {Xiu-Zhe Luo, Jin-Guo Liu, Pan Zhang, and Lei Wang, +Quantum 4, 341 (2020). +We introduce \${\textbackslash}texttt\{Yao\}\$, an extensible, efficient open-source framework for quantum algorithm design. \${\textbackslash}texttt\{Yao\}\$ features generic and differentiable programming of quantum circuits. It a…}, + language = {en-GB}, + urldate = {2023-03-23}, + journal = {Quantum}, + author = {Luo, Xiu-Zhe and Liu, Jin-Guo and Zhang, Pan and Wang, Lei}, + month = oct, + year = {2020}, + note = {Publisher: Verein zur Förderung des Open Access Publizierens in den Quantenwissenschaften}, + pages = {341}, +} + +@article{Pan2022, + title = {Simulation of {Quantum} {Circuits} {Using} the {Big}-{Batch} {Tensor} {Network} {Method}}, + volume = {128}, + url = {https://link.aps.org/doi/10.1103/PhysRevLett.128.030501}, + doi = {10.1103/PhysRevLett.128.030501}, + abstract = {We propose a tensor network approach to compute amplitudes and probabilities for a large number of correlated bitstrings in the final state of a quantum circuit. As an application, we study Google’s Sycamore circuits, which are believed to be beyond the reach of classical supercomputers and have been used to demonstrate quantum supremacy. By employing a small computational cluster containing 60 graphical processing units (GPUs), we compute exact amplitudes and probabilities of 2×106 correlated bitstrings with some entries fixed (which span a subspace of the output probability distribution) for the Sycamore circuit with 53 qubits and 20 cycles. The obtained results verify the Porter-Thomas distribution of the large and deep quantum circuits of Google, provide datasets and benchmarks for developing approximate simulation methods, and can be used for spoofing the linear cross entropy benchmark of quantum supremacy. Then we extend the proposed big-batch method to a full-amplitude simulation approach that is more efficient than the existing Schrödinger method on shallow circuits and the Schrödinger-Feynman method in general, enabling us to obtain the state vector of Google’s simplifiable circuit with n=43 qubits and m=14 cycles using only one GPU. We also manage to obtain the state vector for Google’s simplifiable circuits with n=50 qubits and m=14 cycles using a small GPU cluster, breaking the previous record on the number of qubits in full-amplitude simulations. Our method is general in computing bitstring probabilities for a broad class of quantum circuits and can find applications in the verification of quantum computers. We anticipate that our method will pave the way for combining tensor network–based classical computations and near-term quantum computations for solving challenging problems in the real world.}, + number = {3}, + urldate = {2023-02-09}, + journal = {Physical Review Letters}, + author = {Pan, Feng and Zhang, Pan}, + month = jan, + year = {2022}, + note = {Publisher: American Physical Society}, + pages = {030501}, +} + +@article{Kalachev2021, + title = {Recursive {Multi}-{Tensor} {Contraction} for {XEB} {Verification} of {Quantum} {Circuits}}, + url = {http://arxiv.org/abs/2108.05665}, + abstract = {The computational advantage of noisy quantum computers have been demonstrated by sampling the bitstrings of quantum random circuits. An important issue is how the performance of quantum devices could be quantified in the so-called "supremacy regime". The standard approach is through the linear cross entropy (XEB), where the theoretical value of the probability is required for each bitstring. However, the computational cost of XEB grows exponentially. So far, random circuits of the 53-qubit Sycamore chip was verified up to 10 cycles of gates only; the XEB fidelities of deeper circuits were approximated with simplified circuits instead. Here we present a multi-tensor contraction algorithm for speeding up the calculations of XEB of quantum circuits, where the computational cost can be significantly reduced through a recursive manner with some form of memoization. As a demonstration, we analyzed the experimental data of the 53-qubit Sycamore chip and obtained the exact values of the corresponding XEB fidelities up to 16 cycles using only moderate computing resources (few GPUs). If the algorithm was implemented on the Summit supercomputer, we estimate that for the 20-cycles supremacy circuits, it would only cost 7.5 days, which is several orders of magnitudes lower than previously estimated in the literature.}, + author = {Kalachev, Gleb and Panteleev, Pavel and Yung, Man-Hong}, + year = {2021}, + note = {arXiv: 2108.05665}, + pages = {1--9}, +} + +@article{Markov2008, + title = {Simulating {Quantum} {Computation} by {Contracting} {Tensor} {Networks}}, + volume = {38}, + issn = {0097-5397}, + url = {https://epubs.siam.org/doi/abs/10.1137/050644756}, + doi = {10.1137/050644756}, + abstract = {Adiabatic quantum computation has recently attracted attention in the physics and computer science communities, but its computational power was unknown. We describe an efficient adiabatic simulation of any given quantum algorithm, which implies that the adiabatic computation model and the conventional quantum computation model are polynomially equivalent. Our result can be extended to the physically realistic setting of particles arranged on a two‐dimensional grid with nearest neighbor interactions. The equivalence between the models allows stating the main open problems in quantum computation using well‐studied mathematical objects such as eigenvectors and spectral gaps of sparse matrices.}, + number = {3}, + urldate = {2023-10-05}, + journal = {SIAM Journal on Computing}, + author = {Markov, Igor L. and Shi, Yaoyun}, + month = jan, + year = {2008}, + note = {Publisher: Society for Industrial and Applied Mathematics}, + pages = {963--981}, +} diff --git a/lib/YaoToEinsum/LICENSE b/lib/YaoToEinsum/LICENSE new file mode 100644 index 000000000..75029780e --- /dev/null +++ b/lib/YaoToEinsum/LICENSE @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2021 GiggleLiu and contributors + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/lib/YaoToEinsum/Project.toml b/lib/YaoToEinsum/Project.toml new file mode 100644 index 000000000..200d6e320 --- /dev/null +++ b/lib/YaoToEinsum/Project.toml @@ -0,0 +1,21 @@ +name = "YaoToEinsum" +uuid = "9b173c7b-dc24-4dc5-a0e1-adab2f7b6ba9" +authors = ["GiggleLiu and contributors"] +version = "0.2.3" + +[deps] +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +OMEinsum = "ebe7aa44-baf0-506c-a96f-8464559b3922" +YaoBlocks = "418bc28f-b43b-5e0b-a6e7-61bbc1a2c1df" + +[compat] +OMEinsum = "0.8" +YaoBlocks = "0.13" +julia = "1.9" + +[extras] +SymEngine = "123dc426-2d89-5057-bbad-38513e3affd8" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[targets] +test = ["Test", "SymEngine"] diff --git a/lib/YaoToEinsum/README.md b/lib/YaoToEinsum/README.md new file mode 100644 index 000000000..f689653c5 --- /dev/null +++ b/lib/YaoToEinsum/README.md @@ -0,0 +1,14 @@ +# YaoToEinsum + +Convert [Yao](https://github.com/QuantumBFS/Yao.jl) circuit to tensor networks (einsum). + +## Installation + +`YaoToEinsum` is a [Julia language](https://julialang.org/) package. To install `YaoToEinsum`, please [open Julia's interactive session (known as REPL)](https://docs.julialang.org/en/v1/manual/getting-started/) and press ] key in the REPL to use the package mode, then type the following command + +```julia +pkg> add YaoToEinsum +``` + +## Using +Please refer to the [documentation](). \ No newline at end of file diff --git a/lib/YaoToEinsum/src/Core.jl b/lib/YaoToEinsum/src/Core.jl new file mode 100644 index 000000000..3d36cc22e --- /dev/null +++ b/lib/YaoToEinsum/src/Core.jl @@ -0,0 +1,62 @@ +""" + TensorNetwork + +A (generalized) tensor network representation of a quantum circuit. + +### Fields +* `code::AbstractEinsum`: The einsum code. +* `tensors::Vector`: The tensors in the network. +""" +struct TensorNetwork + code::AbstractEinsum + tensors::Vector +end +function Base.show(io::IO, c::TensorNetwork) + print(io, "TensorNetwork") + print(io, "\n") + print(io, contraction_complexity(c)) +end +function Base.show(io::IO, ::MIME"text/plain", c::TensorNetwork) + Base.show(io, c) +end +function Base.iterate(c::TensorNetwork, state=1) + if state > 2 + return nothing + elseif state == 1 + return (c.code, 2) + else + return (c.tensors, 3) + end +end + +""" + contract(c::TensorNetwork) + +Contract the tensor network, and return the result tensor. +""" +function contract(c::TensorNetwork) + return c.code(c.tensors...; size_info=uniformsize(c.code, 2)) +end + +""" + optimize_code(c::TensorNetwork, optimizer=TreeSA()) + +Optimize the code of the tensor network. + +### Arguments +* `c::TensorNetwork`: The tensor network. +* `optimizer::Optimizer`: The optimizer to use, default is `TreeSA()`. Please check [OMEinsumContractors.jl](https://github.com/TensorBFS/OMEinsumContractionOrders.jl) for more information. +""" +function OMEinsum.optimize_code(c::TensorNetwork, args...) + optcode = optimize_code(c.code, uniformsize(c.code, 2), args...) + return TensorNetwork(optcode, c.tensors) +end + +""" + contraction_complexity(c::TensorNetwork) + +Return the contraction complexity of the tensor network. +""" +function OMEinsum.contraction_complexity(c::TensorNetwork) + return contraction_complexity(c.code, uniformsize(c.code, 2)) +end \ No newline at end of file diff --git a/lib/YaoToEinsum/src/YaoToEinsum.jl b/lib/YaoToEinsum/src/YaoToEinsum.jl new file mode 100644 index 000000000..1b10d22b3 --- /dev/null +++ b/lib/YaoToEinsum/src/YaoToEinsum.jl @@ -0,0 +1,13 @@ +module YaoToEinsum + +using YaoBlocks, YaoBlocks.YaoArrayRegister, OMEinsum +using LinearAlgebra + +export yao2einsum +export TensorNetwork, optimize_code, contraction_complexity, contract +export TreeSA + +include("Core.jl") +include("circuitmap.jl") + +end diff --git a/lib/YaoToEinsum/src/circuitmap.jl b/lib/YaoToEinsum/src/circuitmap.jl new file mode 100644 index 000000000..c61901077 --- /dev/null +++ b/lib/YaoToEinsum/src/circuitmap.jl @@ -0,0 +1,215 @@ +struct EinBuilder{T} + slots::Vector{Int} + labels::Vector{Vector{Int}} + tensors::Vector{AbstractArray{T}} + maxlabel::Base.RefValue{Int} +end + +YaoBlocks.nqubits(eb::EinBuilder) = length(eb.slots) +function add_tensor!(eb::EinBuilder{T}, tensor::AbstractArray{T,N}, labels::Vector{Int}) where {N,T} + @assert N == length(labels) + push!(eb.tensors, tensor) + push!(eb.labels, labels) +end + +function EinBuilder(::Type{T}, n::Int) where T + EinBuilder(collect(1:n), Vector{Int}[], AbstractArray{T}[], Ref(n)) +end +newlabel!(eb::EinBuilder) = (eb.maxlabel[] += 1; eb.maxlabel[]) + +function add_gate!(eb::EinBuilder{T}, b::PutBlock{D,C}) where {T,D,C} + return add_matrix!(eb, C, mat(T, b.content), collect(b.locs)) +end +# general and diagonal gates +function add_matrix!(eb::EinBuilder{T}, k::Int, m::AbstractMatrix, locs::Vector) where T + if isdiag(m) + add_tensor!(eb, reshape(Vector{T}(diag(m)), fill(2, k)...), eb.slots[locs]) + elseif m isa YaoBlocks.OuterProduct # low rank + nlabels = [newlabel!(eb) for _=1:k] + K = rank(m) + if K == 1 # projector + add_tensor!(eb, reshape(Vector{T}(m.right), fill(2, k)...), [eb.slots[locs]...]) + add_tensor!(eb, reshape(Vector{T}(m.left), fill(2, k)...), [nlabels...]) + eb.slots[locs] .= nlabels + else + midlabel = newlabel!(eb) + add_tensor!(eb, reshape(Matrix{T}(m.right), fill(2, k)..., K), [eb.slots[locs]..., midlabel]) + add_tensor!(eb, reshape(Matrix{T}(m.left), fill(2, k)..., K), [nlabels..., midlabel]) + eb.slots[locs] .= nlabels + end + else + nlabels = [newlabel!(eb) for _=1:k] + add_tensor!(eb, reshape(Matrix{T}(m), fill(2, 2k)...), [nlabels..., eb.slots[locs]...]) + eb.slots[locs] .= nlabels + end + return eb +end +# swap gate +function add_gate!(eb::EinBuilder{T}, b::PutBlock{2,2,ConstGate.SWAPGate}) where {T} + lj = eb.slots[b.locs[2]] + eb.slots[b.locs[2]] = eb.slots[b.locs[1]] + eb.slots[b.locs[1]] = lj + return eb +end + +# projection gate, todo: generalize to arbitrary low rank gate +function add_gate!(eb::EinBuilder{T}, b::PutBlock{2,1,ConstGate.P0Gate}) where {T} + add_matrix!(eb, 1, YaoBlocks.OuterProduct(T[1, 0], T[1, 0]), collect(b.locs)) + return eb +end + +# projection gate, todo: generalize to arbitrary low rank gate +function add_gate!(eb::EinBuilder{T}, b::PutBlock{2,1,ConstGate.P1Gate}) where {T} + add_matrix!(eb, 1, YaoBlocks.OuterProduct(T[0, 1], T[0, 1]), collect(b.locs)) + return eb +end + + +# control gates +function add_gate!(eb::EinBuilder{T}, b::ControlBlock{BT,C,M}) where {T, BT,C,M} + return add_controlled_matrix!(eb, M, mat(T, b.content), collect(b.locs), collect(b.ctrl_locs), collect(b.ctrl_config)) +end +function add_controlled_matrix!(eb::EinBuilder{T}, k::Int, m::AbstractMatrix, locs::Vector, control_locs, control_vals) where T + if length(control_locs) == 0 + return add_matrix!(eb, k, m, locs) + end + sig = eb.slots[control_locs[1]] + val = control_vals[1] + for i=1:length(control_locs)-1 + newsig = newlabel!(eb) + add_tensor!(eb, and_gate(T, control_vals[i+1], val), [newsig,eb.slots[control_locs[i+1]],sig]) + sig = newsig + val = 1 + end + if !isdiag(m) + t1 = reshape(Matrix{T}(m), fill(2, 2k)...) + t2 = reshape(Matrix{T}(I, 1<1, 2=>1, 3=>0, 4=>1)` specifies a product state `|1⟩⊗|1⟩⊗|0⟩⊗|1⟩`. + - In the second interface, a state is specified as an `ArrayReg`, e.g. `Dict(1=>rand_state(1), 2=>rand_state(1))`. +If any qubit in initial state or final state is not specified, it will be treated as a free leg in the tensor network. +* `optimizer` is the optimizer used to optimize the tensor network. The default is `TreeSA()`. +Please check [OMEinsumContractors.jl](https://github.com/TensorBFS/OMEinsumContractionOrders.jl) for more information. + + +```jldoctest +julia> using Yao + +julia> c = chain(3, put(3, 2=>X), put(3, 1=>Y), control(3, 1, 3=>Y)) +nqubits: 3 +chain +├─ put on (2) +│ └─ X +├─ put on (1) +│ └─ Y +└─ control(1) + └─ (3,) Y + + +julia> yao2einsum(c; initial_state=Dict(1=>0, 2=>1), final_state=Dict(1=>ArrayReg([0.6, 0.8im]), 2=>1)) +TensorNetwork +Time complexity: 2^4.700439718141093 +Space complexity: 2^2.0 +Read-write complexity: 2^6.0 +``` +""" +function yao2einsum(circuit::AbstractBlock{D}; initial_state::Dict=Dict{Int,Int}(), final_state::Dict=Dict{Int,Int}(), optimizer=TreeSA()) where {D} + T = promote_type(ComplexF64, dict_regtype(initial_state), dict_regtype(final_state), YaoBlocks.parameters_eltype(circuit)) + vec_initial_state = Dict{keytype(initial_state),ArrayReg{D,T}}([k=>render_single_qubit_state(T, v) for (k, v) in initial_state]) + vec_final_state = Dict{keytype(final_state),ArrayReg{D,T}}([k=>render_single_qubit_state(T, v) for (k, v) in final_state]) + yao2einsum(circuit, vec_initial_state, vec_final_state, optimizer) +end +dict_regtype(d::Dict) = promote_type(_regtype.(values(d))...) +_regtype(::ArrayReg{D,VT}) where {D,VT} = VT +_regtype(::Int) = ComplexF64 +render_single_qubit_state(::Type{T}, x::Int) where T = x == 0 ? zero_state(T, 1) : product_state(T, bit"1") +render_single_qubit_state(::Type{T}, x::ArrayReg) where T = ArrayReg(collect(T, statevec(x))) + +function yao2einsum(circuit::AbstractBlock{D}, initial_state::Dict{Int,<:ArrayReg{D,T}}, final_state::Dict{Int,<:ArrayReg{D,T}}, optimizer) where {D,T} + v_initial_state = Dict{Vector{Int}, ArrayReg{D,T}}([[k]=>v for (k, v) in initial_state]) + v_final_state = Dict{Vector{Int}, ArrayReg{D, T}}([[k]=>v for (k, v) in final_state]) + yao2einsum(circuit, v_initial_state, v_final_state, optimizer) +end +function yao2einsum(circuit::AbstractBlock{D}, initial_state::Dict{Vector{Int},<:ArrayReg{D,T}}, final_state::Dict{Vector{Int},<:ArrayReg{D,T}}, optimizer) where {D,T} + n = nqubits(circuit) + eb = EinBuilder(T, n) + openindices = add_states!(eb, initial_state) + add_gate!(eb, circuit) + openindices2 = add_states!(eb, final_state; conjugate=true) + network = build_einsum(eb, vcat(openindices2, openindices)) + return optimizer === nothing ? network : optimize_code(network, optimizer, MergeVectors()) +end +function check_state_spec(state::Dict{Vector{Int},<:ArrayReg{D,T}}, n::Int) where {D,T} + iks = collect(Int, vcat(keys(state)...)) + @assert length(unique(iks)) == length(iks) "state qubit indices must be unique" + @assert all(1 .<= iks .<= n) "state qubit indices must be in the range 1 to $n" + return iks +end +function add_states!(eb::EinBuilder{T}, states::Dict; conjugate=false) where {T} + n = nqubits(eb) + unique_indices = check_state_spec(states, n) + openindices = eb.slots[setdiff(1:n, unique_indices)] + for (k, state) in states + add_tensor!(eb, (conjugate ? conj : identity).(dropdims(hypercubic(state); dims=length(k)+1)), eb.slots[k]) + end + return openindices +end + +function build_einsum(eb::EinBuilder, openindices) + return TensorNetwork(EinCode(eb.labels, openindices), eb.tensors) +end + diff --git a/lib/YaoToEinsum/test/circuitmap.jl b/lib/YaoToEinsum/test/circuitmap.jl new file mode 100644 index 000000000..56b86cfd3 --- /dev/null +++ b/lib/YaoToEinsum/test/circuitmap.jl @@ -0,0 +1,350 @@ +module CircuitMapTest +using YaoToEinsum +using Test, OMEinsum +using YaoBlocks, YaoBlocks.YaoArrayRegister +using SymEngine + +const CPhaseGate{T} = ControlBlock{<:ShiftGate{T},<:Any} + +@const_gate ISWAP = PermMatrix([1,3,2,4], [1,1.0im,1.0im,1]) +@const_gate SqrtX = [0.5+0.5im 0.5-0.5im; 0.5-0.5im 0.5+0.5im] +@const_gate SqrtY = [0.5+0.5im -0.5-0.5im; 0.5+0.5im 0.5+0.5im] +# √W is a non-Clifford gate +@const_gate SqrtW = mat(rot((X+Y)/sqrt(2), π/2)) + +""" + singlet_block(θ::Real, ϕ::Real) + +The circuit block for initialzing a singlet state. +""" +singlet_block() = chain(put(2, 1=>chain(X, H)), control(2, -1, 2=>X)) + + +mutable struct FSimGate{T<:Number} <: PrimitiveBlock{2} + theta::T + phi::T +end +YaoBlocks.nqudits(fs::FSimGate) = 2 +YaoBlocks.print_block(io::IO, block::FSimGate) = print(io, "FSim(θ=$(block.theta), ϕ=$(block.phi))") + +function Base.:(==)(fs1::FSimGate, fs2::FSimGate) + return fs1.theta == fs2.theta && fs1.phi == fs2.phi +end + +function YaoBlocks.mat(::Type{T}, fs::FSimGate) where T + θ, ϕ = fs.theta, fs.phi + T[1 0 0 0; + 0 cos(θ) -im*sin(θ) 0; + 0 -im*sin(θ) cos(θ) 0; + 0 0 0 exp(-im*ϕ)] +end + +YaoBlocks.iparams_eltype(::FSimGate{T}) where T = T +YaoBlocks.getiparams(fs::FSimGate{T}) where T = (fs.theta, fs.phi) +function YaoBlocks.setiparams!(fs::FSimGate{T}, θ, ϕ) where T + fs.theta = θ + fs.phi = ϕ + return fs +end + +YaoBlocks.@dumpload_fallback FSimGate FSimGate +YaoBlocks.Optimise.to_basictypes(fs::FSimGate) = fsim_block(fs.theta, fs.phi) + +""" + fsim_block(θ::Real, ϕ::Real) + +The circuit representation of FSim gate. +""" +function fsim_block(θ::Real, ϕ::Real) + if θ ≈ π/2 + return cphase(2,2,1,-ϕ)*SWAP*rot(kron(Z,Z), -π/2)*put(2,1=>phase(-π/4)) + else + return cphase(2,2,1,-ϕ)*rot(SWAP,2*θ)*rot(kron(Z,Z), -θ)*put(2,1=>phase(θ/2)) + end +end + +cphase(nbits::Int, i::Int, j::Int, θ::T) where T = control(nbits, i, j=>shift(θ)) +qft_circuit(::Type{T}, n::Int) where T = chain(n, hcphases(T, n, i) for i = 1:n) +hcphases(::Type{T}, n, i) where T = chain(n, i==j ? put(i=>H) : cphase(n, j, i, 2T(π)/T(2^(j-i+1))) for j in i:n); +qft_circuit(n::Int) = qft_circuit(Float64, n) + +entangler_google53(::Type{T}, nbits::Int, i::Int, j::Int) where T = put(nbits, (i,j)=>FSimGate(T(π)/2, T(π)/6)) + +struct Lattice53 + labels::Matrix{Int} +end + +function Lattice53(;nbits::Int=53) + config = ones(Bool, 5, 12) + config[end,2:2:end] .= false + config[1, 7] = false + labels = zeros(Int, 5, 12) + k = 0 + for (i,c) in enumerate(config) + if c + k += 1 + labels[i] = k + k>=nbits && break + end + end + return Lattice53(labels) +end + +nbits(lattice::Lattice53) = maximum(lattice.labels) + +function Base.getindex(lattice::Lattice53, i, j) + 1<=i<=size(lattice.labels, 1) && 1<=j<=size(lattice.labels, 2) ? lattice.labels[i,j] : 0 +end +upperleft(lattice::Lattice53,i,j) = lattice[i-j%2,j-1] +lowerleft(lattice::Lattice53,i,j) = lattice[i+(j-1)%2,j-1] +upperright(lattice::Lattice53,i,j) = lattice[i-j%2,j+1] +lowerright(lattice::Lattice53,i,j) = lattice[i+(j-1)%2,j+1] + +function pattern53(lattice::Lattice53, chr::Char) + res = Tuple{Int,Int}[] + # i0, di, j0, dj and direction + di = 1 + (chr>'D') + dj = 2 - (chr>'D') + j0 = 1 + min(dj-1, mod(chr-'A',2)) + direction = 'C'<=chr<='F' ? lowerright : upperright + for j=j0:dj:12 + i0 = chr>'D' ? mod((chr-'D') + (j-(chr>='G'))÷2, 2) : 1 + for i = i0:di:5 + src = lattice[i, j] + dest = direction(lattice, i, j) + src!=0 && dest !=0 && push!(res, (src, dest)) + end + end + return res +end + +function print_lattice53(lattice, pattern) + for i_=1:10 + i = (i_+1)÷2 + for j=1:12 + if i_%2 == j%2 && lattice[i,j]!=0 + print(" ∘ ") + else + print(" ") + end + end + println() + for j=1:12 + if i_%2 == j%2 && lattice[i,j]!=0 + hasll = (lowerleft(lattice, i, j), lattice[i,j]) in pattern + haslr = (lattice[i,j], lowerright(lattice, i, j)) in pattern + print(hasll ? "/ " : " ") + print(haslr ? " \\" : " ") + else + print(" ") + end + end + println() + end +end + +""" + rand_google53([T=Float64], depth::Int; nbits=53) -> AbstactBlock + +Google supremacy circuit with 53 qubits, also know as the Sycamore quantum supremacy circuits. `T` is the parameter type. + +References +------------------------- +* Arute, Frank, et al. "Quantum supremacy using a programmable superconducting processor." Nature 574.7779 (2019): 505-510. +""" +rand_google53(depth::Int; nbits::Int=53) = rand_google53(Float64, depth; nbits) +function rand_google53(::Type{T}, depth::Int; nbits::Int=53) where T + c = chain(nbits) + lattice = Lattice53(nbits=nbits) + k = 0 + for pattern in Iterators.cycle(['A', 'B', 'C', 'D', 'C', 'D', 'A', 'B']) + push!(c, rand_google53_layer(T, lattice, pattern)) + k += 1 + k>=depth && break + end + return c +end + +function rand_google53_layer(::Type{T}, lattice, pattern) where T + nbit = nbits(lattice) + chain(nbit, chain(nbit, [put(nbit, i=>rand([SqrtW, SqrtX, SqrtY])) for i=1:nbit]), + chain(nbit, [entangler_google53(T, nbit,i,j) for (i,j) in pattern53(lattice, pattern)]) + ) +end + +""" + pair_ring(n::Int) -> Vector + +Pair ring entanglement layout. +""" +pair_ring(n::Int) = [i=>mod(i, n)+1 for i=1:n] + +""" + pair_square(m::Int, n::Int) -> Vector + +Pair square entanglement layout. +""" +function pair_square(m::Int, n::Int; periodic=false) + res = Vector{Pair{Int, Int}}(undef, (m-!periodic)*n+m*(n-!periodic)) + li = LinearIndices((m, n)) + k = 1 + for i = 1:2:m, j=1:n + if periodic || i li[i%m+1, j] + k+=1 + end + end + for i = 2:2:m, j=1:n + if periodic || i li[i%m+1, j] + k+=1 + end + end + for i = 1:m, j=1:2:n + if periodic || j li[i, j%n+1] + k+=1 + end + end + for i = 1:m, j=2:2:n + if periodic || j li[i, j%n+1] + k+=1 + end + end + res +end + +###################### rotor and rotorset ##################### +""" + merged_rotor(noleading::Bool=false, notrailing::Bool=false) -> ChainBlock{1, ComplexF64} + +Single qubit arbitrary rotation unit, set parameters notrailing, noleading true to remove trailing and leading Z gates. + +!!! note + + Here, `merged` means `Rz(η)⋅Rx(θ)⋅Rz(ξ)` are multiplied first, this kind of operation if now allowed in differentiable + circuit with back-propagation (`:BP`) mode (just because we are lazy to implement it!). + But is a welcoming component in quantum differentiation. +""" +merged_rotor(::Type{T}, noleading::Bool=false, notrailing::Bool=false) where T = noleading ? (notrailing ? Rx(zero(T)) : chain(Rx(zero(T)), Rz(zero(T)))) : (notrailing ? chain(Rz(zero(T)), Rx(zero(T))) : chain(Rz(zero(T)), Rx(zero(T)), Rz(zero(T)))) + +""" + rotor(nbit::Int, ibit::Int, noleading::Bool=false, notrailing::Bool=false) -> ChainBlock{nbit, ComplexF64} + +Arbitrary rotation unit (put in `nbit` space), set parameters notrailing, noleading true to remove trailing and leading Z gates. +""" +function rotor(::Type{T}, nbit::Int, ibit::Int, noleading::Bool=false, notrailing::Bool=false) where T + rt = chain(nbit, [put(nbit, ibit=>Rz(zero(T))), put(nbit, ibit=>Rx(zero(T))), put(nbit, ibit=>Rz(zero(T)))]) + noleading && popfirst!(rt) + notrailing && pop!(rt) + rt +end + +rotorset(::Type{T}, ::Val{:Merged}, nbit::Int, noleading::Bool=false, notrailing::Bool=false) where T = chain(nbit, [put(nbit, j=>merged_rotor(T, noleading, notrailing)) for j=1:nbit]) +rotorset(::Type{T}, ::Val{:Split}, nbit::Int, noleading::Bool=false, notrailing::Bool=false) where T = chain(nbit, [rotor(T, nbit, j, noleading, notrailing) for j=1:nbit]) +rotorset(::Type{T}, mode::Symbol, nbit::Int, noleading::Bool=false, notrailing::Bool=false) where T = rotorset(T, Val(mode), nbit, noleading, notrailing) + +""" + variational_circuit([T=Float64], nbit[, nlayer][, pairs]; mode=:Split, do_cache=false, entangler=cnot) + +A kind of widely used differentiable quantum circuit, angles in the circuit is randomely initialized. + +### Arguments + +* `T` is the parameter type. +* `pairs` is list of `Pair`s for entanglers in a layer, default to `pair_ring` structure, +* `mode` can be :Split or :Merged, +* `do_cache` decides whether cache the entangler matrix or not, +* `entangler` is a constructor returns a two qubit gate, `f(n,i,j) -> gate`. + The default value is `cnot(n,i,j)`. + +### References + +1. Kandala, A., Mezzacapo, A., Temme, K., Takita, M., Chow, J. M., & Gambetta, J. M. (2017). Hardware-efficient Quantum Optimizer for Small Molecules and Quantum Magnets. Nature Publishing Group, 549(7671), 242–246. https://doi.org/10.1038/nature23879. +""" +function variational_circuit(::Type{T}, nbit, nlayer, pairs; mode=:Split, do_cache=false, entangler=cnot) where T + circuit = chain(nbit) + + ent = chain(nbit, entangler(nbit, i, j) for (i, j) in pairs) + if do_cache + ent = ent |> cache + end + has_param = nparameters(ent) != 0 + for i = 1:(nlayer + 1) + i!=1 && push!(circuit, has_param ? deepcopy(ent) : ent) + push!(circuit, rotorset(T, mode, nbit, i==1, i==nlayer+1)) + end + circuit +end + +variational_circuit(::Type{T}, n::Int; kwargs...) where T = variational_circuit(T, n, 3, pair_ring(n); kwargs...) + +variational_circuit(::Type{T}, nbit::Int, nlayer::Int; kwargs...) where T = variational_circuit(T, nbit, nlayer, pair_ring(nbit); kwargs...) +variational_circuit(nbit::Int; kwargs...) = variational_circuit(Float64, nbit; kwargs...) +variational_circuit(nbit::Int, nlayer::Int; kwargs...) = variational_circuit(Float64, nbit, nlayer; kwargs...) + +@testset "YaoToEinsum.jl" begin + n = 5 + a = rand_unitary(4)[:, 1:2] + a1 = rand_unitary(4)[:, 2] + mb = matblock(OuterProduct(conj.(a), a)) + mb1 = matblock(OuterProduct(conj.(a1), a1)) + for c in [put(n, 2=>Y), put(n, 2=>ConstGate.P0), put(n, (3,2)=>mb1), put(n, (2, 1)=>mb), put(n, 2=>ConstGate.P1), put(n, (5,3)=>SWAP), put(n, (4,2)=>ConstGate.CNOT), put(n, (2,3,1)=>kron(ConstGate.CNOT, X)), + put(n, 2=>Z), control(n, -3, 2=>X), control(n, 3, 2=>X), control(n, (2, -1), 3=>Y), control(n, (4,1,-2), 5=>Z)] + @show c + C = chain([put(n, i=>Rx(rand()*2π)) for i=1:n]..., c) + code, xs = yao2einsum(C; optimizer=nothing) + optcode = optimize_code(code, uniformsize(code, 2), GreedyMethod()) + @test reshape(optcode(xs...; size_info=uniformsize(code, 2)), 1<rand_state(1) for i=1:n]) + reg = join([initial_state[i] for i=n:-1:1]...) + reg |> c + inner = (2,3) + focus!(reg, inner) + for final_state in [Dict([i=>rand_state(1) for i in inner]), Dict([i=>1 for i in inner])] + freg = join(YaoToEinsum.render_single_qubit_state(ComplexF64, final_state[3]), YaoToEinsum.render_single_qubit_state(ComplexF64, final_state[2])) + net = yao2einsum(c; initial_state=initial_state, final_state=final_state, optimizer=TreeSA(nslices=3)) + println(net) + @test vec(contract(net)) ≈ vec(statevec(freg)' * state(reg)) + end +end + +@testset "symbolic" begin + n = 5 + c = qft_circuit(n) + initial_state = Dict([i=>zero_state(Basic, 1) for i=1:n]) + code, xs = yao2einsum(c; initial_state=initial_state) + @test eltype(xs) == AbstractArray{Basic} +end + +@testset "fix to basic type" begin + c = chain(kron(X,X)) + @test (yao2einsum(c) |> first) isa OMEinsum.SlicedEinsum +end + +@testset "multiple qubit states" begin + n = 4 + reg1 = rand_state(2) + reg2 = rand_state(2) + reg3 = rand_state(1) + reg4 = rand_state(3) + c = chain(4) + code, xs = yao2einsum(c; initial_state=Dict([1, 2]=>reg1, [3, 4]=>reg2), final_state=Dict([1]=>reg3, [2,3,4]=>reg4)) + @test code(xs...; size_info=uniformsize(code, 2))[] ≈ join(reg4, reg3)' * join(reg2, reg1) +end +end diff --git a/lib/YaoToEinsum/test/runtests.jl b/lib/YaoToEinsum/test/runtests.jl new file mode 100644 index 000000000..fd313091d --- /dev/null +++ b/lib/YaoToEinsum/test/runtests.jl @@ -0,0 +1,5 @@ +using YaoToEinsum, Test + +@testset "circuitmap" begin + include("circuitmap.jl") +end \ No newline at end of file diff --git a/notebooks/getting-started/assignment1.jl b/notebooks/getting-started/assignment1.jl index 6733e9078..d819b7156 100755 --- a/notebooks/getting-started/assignment1.jl +++ b/notebooks/getting-started/assignment1.jl @@ -88,7 +88,7 @@ Matrix(bellcircuit) == Matrix(chain(2, put(1=>H), control(1,2=>X))) ? md"✅" : # ╔═╡ a59c3f42-0667-11eb-34e2-c9ab7e9080f2 begin reversebellcircuit = chain(2) #complete the reverse bell circuit - plot(bellcircuit) + plot(reversebellcircuit) end # ╔═╡ c2a9bb30-067d-11eb-0594-5d416dc1e763 @@ -101,7 +101,7 @@ begin end # ╔═╡ 5d52aca4-067f-11eb-3174-a966da065d76 -Matrix(bell_and_reverse_bell_circuit) == Matrix(chain(2, put(1=>H), control(1, 2=>H), control(1, 2=>H), put(1=>H))) ? md"✅" : md"❌" +Matrix(bell_and_reverse_bell_circuit) == Matrix(chain(2, put(1=>H), control(1, 2=>X), control(1, 2=>X), put(1=>H))) ? md"✅" : md"❌" # ╔═╡ d87106c6-067d-11eb-3beb-8bfea8d168a3 begin @@ -135,7 +135,7 @@ md"_**Note:**_ If you're using ` ArrayReg(bit\" \") ` to create the qubits, the md"_**Assignment 3:**_" # ╔═╡ 131c831a-0681-11eb-2bab-3d4cf5b23006 -md"Suppose that Alice and Bob have 2 pairs of entangled qubits. Both the pairs are in the state $ \frac{|00> \;+\; |11>}{\sqrt2} $. Suppose Alice has an extra qubit, with completely random state." +md"Suppose that Alice and Bob have 2 pairs of entangled qubits. Both the pairs are in the state `` \frac{|00> \;+\; |11>}{\sqrt2} ``. Suppose Alice has an extra qubit, with completely random state." # ╔═╡ c89d84b4-0681-11eb-3c22-4def72fc3477 md"1. _Make the circuit for quantum teleportation_ @@ -170,7 +170,7 @@ end # ╔═╡ 88f1d456-06df-11eb-1827-912c481bb7a4 input = 0 -#Remove the "0". Pass the above created qubit through the teleportation circuit, measuring the first two qubits using measure_remove!() function. Hint: qubit |> circuit |> r->measure_remove!(r, m:n) +#Remove the "0". Pass the above created qubit through the teleportation circuit, measuring the first two qubits using RemoveMeasured() function. Hint: qubit |> circuit |> r->measure!(RemoveMeasured(), r, m:n) # ╔═╡ 79280796-06e3-11eb-1e5c-b5cea995a88a (input == bit"00" || input == bit"11" || input==bit"10" || input==bit"01") && (typeof(input) != Int64) ? md"✅" : md"❌" @@ -180,13 +180,13 @@ md"The above measurement will act as the information for superdense coding circu # ╔═╡ 26cc4008-06e0-11eb-0cdf-51954af7254d begin - if(input==bit"00") + if input==bit"00" superdense_coding_circuit = chain(2) - elseif(input==bit"01") + elseif input==bit"01" superdense_coding_circuit = chain(2) - elseif(input==bit"10") + elseif input==bit"10" superdense_coding_circuit = chain(2) - elseif(input==bit"11") + elseif input==bit"11" superdense_coding_circuit = chain(2) end #Remember the bits are read from right to left @@ -195,17 +195,15 @@ begin end # ╔═╡ dc0b86be-06e4-11eb-24db-4f56cb821e0f -begin - if(input==bit"00") - Matrix(superdense_coding_circuit) == Matrix(chain(2)) - elseif(input==bit"01") - Matrix(superdense_coding_circuit) == Matrix(chain(2, put(1=>Z))) - elseif(input==bit"10") - Matrix(superdense_coding_circuit) == Matrix(chain(2, put(1=>X))) - elseif(input==bit"11") - Matrix(superdense_coding_circuit) == Matrix(chain(2, put(1=>Y))) - end && typeof(input) != Int64 ? md"✅" : md"❌" -end +if input==bit"00" + Matrix(superdense_coding_circuit) == Matrix(chain(2)) +elseif input==bit"01" + Matrix(superdense_coding_circuit) == Matrix(chain(2, put(1=>Z))) +elseif input==bit"10" + Matrix(superdense_coding_circuit) == Matrix(chain(2, put(1=>X))) +elseif input==bit"11" + Matrix(superdense_coding_circuit) == Matrix(chain(2, put(1=>Y))) +end && typeof(input) != Int64 ? md"✅" : md"❌" # ╔═╡ 28ccce76-06e1-11eb-1c70-b3d7779a062e Bobs_part = ((Alice_and_Bobs_second_entangled_qubit |> superdense_coding_circuit) |> reversebellcircuit) |> r->measure(r, nshots=1000) @@ -213,29 +211,28 @@ Bobs_part = ((Alice_and_Bobs_second_entangled_qubit |> superdense_coding_circuit # ╔═╡ 2aaa0910-06e2-11eb-1d6c-0fae35763d48 #Bobs_part now contains the information Alice wanted to convey to Bob regarding his qubits. Use this information to make Bob's qubit's state jump to Alice's random qubit's state begin - if(Bobs_part[1]==bit"00") + if Bobs_part[1]==bit"00" Bobs_qubit = input_to_teleportation_circuit |> chain(1) - elseif(Bobs_part[1]==bit"01") + elseif Bobs_part[1]==bit"01" Bobs_qubit = input_to_teleportation_circuit |> chain(1) #Complete the circuit - elseif(Bobs_part[1]==bit"10") + elseif Bobs_part[1]==bit"10" Bobs_qubit = input_to_teleportation_circuit |> chain(1) #Complete the circuit - elseif(Bobs_part[1]==bit"11") + elseif Bobs_part[1]==bit"11" Bobs_qubit = input_to_teleportation_circuit |> chain(1) #Complete the circuit end state(Bobs_qubit) end # ╔═╡ b1392da2-06e5-11eb-3632-7d96a1b81a24 -begin - log(2,length(state(input_to_teleportation_circuit))) == 2 && ((Bobs_part[1]==bit"00" && Bobs_qubit == (input_to_teleportation_circuit |> chain(1))) || (Bobs_part[1]==bit"01" && Bobs_qubit == (input_to_teleportation_circuit |> chain(1, put(1=>Z)))) || ( Bobs_part[1]==bit"10" && Bobs_qubit == (input_to_teleportation_circuit |> chain(1, put(1=>X))) ) || ( Bobs_part[1]==bit"11" && Bobs_qubit == (input_to_teleportation_circuit |> chain(1, put(1=>Y))))) && (sum(Bobs_part .== Bobs_part[1]) == 1000) ? md"✅" : md"❌" - #=elseif() - Bobs_qubit == - elseif() - Bobs_qubit == - elseif() - Bobs_qubit == =# - -end +if Bobs_part[1]==bit"00" + Bobs_qubit == input_to_teleportation_circuit +elseif Bobs_part[1]==bit"01" + Bobs_qubit == (input_to_teleportation_circuit |> ZGate()) +elseif Bobs_part[1]==bit"10" + Bobs_qubit == (input_to_teleportation_circuit |> XGate()) +elseif Bobs_part[1]==bit"11" + Bobs_qubit == (input_to_teleportation_circuit |> YGate()) +end && (length(state(input_to_teleportation_circuit)) == 2) && (sum(Bobs_part .== Bobs_part[1]) == 1000) ? md"✅" : md"❌" # ╔═╡ 83df1840-06e7-11eb-2502-173dac5d4963 md"Seems to work... Although, wouldn't using the information from superdense coding make teleportation pointless." diff --git a/src/EasyBuild/block_extension/FSimGate.jl b/src/EasyBuild/block_extension/FSimGate.jl index 738aaa4ec..dae21cd27 100644 --- a/src/EasyBuild/block_extension/FSimGate.jl +++ b/src/EasyBuild/block_extension/FSimGate.jl @@ -51,3 +51,5 @@ function fsim_block(θ::Real, ϕ::Real) return cphase(2,2,1,-ϕ)*rot(SWAP,2*θ)*rot(kron(Z,Z), -θ)*put(2,1=>phase(θ/2)) end end + +YaoPlots.get_brush_texts(c, b::FSimGate) = (c.gatestyles.g, "FSim($(pretty_angle(b.theta)), $(pretty_angle(b.phi)))") diff --git a/src/EasyBuild/block_extension/shortcuts.jl b/src/EasyBuild/block_extension/shortcuts.jl index c7c31c41d..465e7f861 100644 --- a/src/EasyBuild/block_extension/shortcuts.jl +++ b/src/EasyBuild/block_extension/shortcuts.jl @@ -14,4 +14,9 @@ const CPhaseGate{T} = ControlBlock{<:ShiftGate{T},<:Any} The circuit block for initialzing a singlet state. """ -singlet_block() = chain(put(2, 1=>chain(X, H)), control(2, -1, 2=>X)) \ No newline at end of file +singlet_block() = chain(put(2, 1=>chain(X, H)), control(2, -1, 2=>X)) + +YaoPlots.get_brush_texts(c, ::SqrtWGate) = (c.gatestyles.g, "√W") +YaoPlots.get_brush_texts(c, ::SqrtXGate) = (c.gatestyles.g, "√X") +YaoPlots.get_brush_texts(c, ::SqrtYGate) = (c.gatestyles.g, "√Y") + diff --git a/src/EasyBuild/easybuild.jl b/src/EasyBuild/easybuild.jl index 9e8a673c1..9780918ee 100644 --- a/src/EasyBuild/easybuild.jl +++ b/src/EasyBuild/easybuild.jl @@ -1,6 +1,7 @@ module EasyBuild using YaoBlocks, YaoBlocks.LuxurySparse, YaoBlocks.YaoAPI, YaoBlocks.YaoArrayRegister using YaoBlocks.LinearAlgebra +import YaoPlots using YaoArrayRegister: sparse include("block_extension/blocks.jl") @@ -12,4 +13,4 @@ include("variational_circuit.jl") include("supremacy_circuit.jl") include("google53.jl") include("hadamardtest.jl") -end \ No newline at end of file +end diff --git a/src/Yao.jl b/src/Yao.jl index 2a402bce8..5b29208a5 100644 --- a/src/Yao.jl +++ b/src/Yao.jl @@ -15,8 +15,7 @@ Extensible Framework for Quantum Algorithm Design for Humans. const 幺 = Yao using Reexport -@reexport using YaoArrayRegister, YaoBlocks, YaoSym -export EasyBuild +@reexport using YaoArrayRegister, YaoBlocks, YaoSym, YaoPlots, YaoToEinsum using YaoArrayRegister.BitBasis, YaoAPI using YaoBlocks: @@ -28,8 +27,12 @@ using YaoBlocks: print_title, print_block +export EasyBuild +# CUDA APIs +export cpu, cuzero_state, cuuniform_state, curand_state, cuproduct_state, cughz_state -include("deprecations.jl") +include("cudainterfaces.jl") include("EasyBuild/easybuild.jl") +include("deprecations.jl") end # module diff --git a/src/cudainterfaces.jl b/src/cudainterfaces.jl new file mode 100644 index 000000000..335ba000b --- /dev/null +++ b/src/cudainterfaces.jl @@ -0,0 +1,41 @@ +""" + cuproduct_state([T=ComplexF64], total::Int, bit_config::Integer; nbatch=NoBatch()) + +The GPU version of [`product_state`](@ref). +""" +function cuproduct_state end + +""" + curand_state([T=ComplexF64], n::Int; nbatch=1) + +The GPU version of [`rand_state`](@ref). +""" +function curand_state end + +""" + cuzero_state([T=ComplexF64], n::Int; nbatch=1) + +The GPU version of [`zero_state`](@ref). +""" +function cuzero_state end + +""" + cuuniform_state([T=ComplexF64], n::Int; nbatch=1) + +The GPU version of [`uniform_state`](@ref). +""" +function cuuniform_state end + +""" + cughz_state([T=ComplexF64], n::Int; nbatch=1) + +The GPU version of [`ghz_state`](@ref). +""" +function cughz_state end + +""" + cpu(cureg) + +Download the register state from GPU to CPU. +""" +function cpu end \ No newline at end of file diff --git a/test/easybuild/block_extension/shortcuts.jl b/test/easybuild/block_extension/shortcuts.jl index c0f652363..9ff12c6f5 100644 --- a/test/easybuild/block_extension/shortcuts.jl +++ b/test/easybuild/block_extension/shortcuts.jl @@ -1,6 +1,7 @@ using Test, Yao.EasyBuild, YaoBlocks using YaoBlocks.Optimise: replace_block, to_basictypes, simplify using YaoBlocks: parse_ex +using YaoPlots: vizcircuit, Luxor @testset "gates" begin @test isunitary(FSimGate(0.5, 0.6)) @@ -21,3 +22,7 @@ using YaoBlocks: parse_ex @test Matrix(c) ≈ Matrix(simplify(c_, rules=[to_basictypes])) end + +@testset "regression" begin + @test vizcircuit(put(10, (8,2,3)=>heisenberg(3)), starting_texts=string.(1:10), ending_texts=string.(1:10), show_ending_bar=true) isa Luxor.Drawing +end \ No newline at end of file diff --git a/test/easybuild/hadamardtest.jl b/test/easybuild/hadamardtest.jl index d7053c020..2c533e520 100644 --- a/test/easybuild/hadamardtest.jl +++ b/test/easybuild/hadamardtest.jl @@ -32,8 +32,8 @@ end 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)) + @test hadamard_test(U, reg, 0.0) ≈ real(sandwich(reg, U, reg)) + @test hadamard_test(U, reg, -π/2) ≈ imag(sandwich(reg, U, reg)) reg = zero_state(2) |> EasyBuild.singlet_block() @test single_swap_test(reg, 0) ≈ -1 diff --git a/test/runtests.jl b/test/runtests.jl index 0142031eb..5ec31bf96 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,6 +3,7 @@ using Yao using YaoAPI using YaoArrayRegister using YaoBlocks +using YaoToEinsum using YaoSym using BitBasis using Documenter @@ -12,10 +13,12 @@ using Random include("easybuild/easybuild.jl") end -DocMeta.setdocmeta!(Yao, :DocTestSetup, :(using Yao, YaoAPI, YaoArrayRegister, YaoBlocks, YaoSym, BitBasis); recursive=true) +DocMeta.setdocmeta!(Yao, :DocTestSetup, :(using Yao, YaoAPI, YaoArrayRegister, YaoBlocks, YaoPlots, YaoSym, YaoToEinsum, BitBasis); recursive=true) Documenter.doctest(YaoAPI; manual=false) Documenter.doctest(YaoArrayRegister; manual=false) Documenter.doctest(YaoBlocks; manual=false) Documenter.doctest(YaoSym; manual=false) +Documenter.doctest(YaoPlots; manual=false) +Documenter.doctest(YaoToEinsum; manual=false) Documenter.doctest(Yao; manual=false)