From d86cccca7f856fe8256305734adc4037cd587692 Mon Sep 17 00:00:00 2001 From: Jinguo Liu Date: Tue, 6 Aug 2024 06:04:23 +0800 Subject: [PATCH] Switch plotting backend from PythonCall to CairoMakie (#673) * Switch plotting backend from PythonCall to CairoMakie * fix test and broken documents * fix a marker issue * fix test --- .ci/Project.toml | 1 - .ci/src/example.jl | 2 - CondaPkg.toml | 2 - Project.toml | 3 +- docs/Project.toml | 1 - docs/src/install.md | 4 - docs/src/waveform.md | 27 +++-- examples/1.blockade/Project.toml | 1 - examples/1.blockade/main.jl | 43 ++++---- examples/2.adiabatic/Project.toml | 1 - examples/2.adiabatic/main.jl | 81 ++++++-------- examples/3.quantum-scar/Project.toml | 3 +- examples/3.quantum-scar/main.jl | 78 ++++++-------- examples/4.LGT/Project.toml | 1 - examples/4.LGT/main.jl | 81 ++++++-------- examples/5.MIS/Project.toml | 1 - examples/5.MIS/main.jl | 115 +++++++++----------- examples/6.MWIS/Project.toml | 1 - examples/6.MWIS/main.jl | 33 +++--- examples/7.QMC/Project.toml | 8 +- examples/7.QMC/main.jl | 35 +++--- examples/8.simulating noise/Project.toml | 8 +- examples/8.simulating noise/main.jl | 130 ++++++++++++++--------- lib/BloqadeNoisy/src/problem.jl | 10 +- src/Bloqade.jl | 10 +- src/plots/bitstring_hist.jl | 45 ++++---- src/plots/waveform.jl | 15 ++- 27 files changed, 345 insertions(+), 395 deletions(-) delete mode 100644 CondaPkg.toml diff --git a/.ci/Project.toml b/.ci/Project.toml index 5d94f0d8c..e1748090d 100644 --- a/.ci/Project.toml +++ b/.ci/Project.toml @@ -5,7 +5,6 @@ version = "0.1.0" [deps] Comonicon = "863f3e99-da2a-4334-8734-de3dacbe5542" -CondaPkg = "992eb4ea-22a4-4c89-a5bb-47a3300528ab" CoverageTools = "c36e975a-824b-4404-a568-ef97ca766997" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" LiveServer = "16fef848-5104-11e9-1b77-fb7a48bbb589" diff --git a/.ci/src/example.jl b/.ci/src/example.jl index bbd60f5e8..006ab6898 100644 --- a/.ci/src/example.jl +++ b/.ci/src/example.jl @@ -141,7 +141,6 @@ in parallel. example_dir = root_dir("examples") script = """ using Pkg - using CondaPkg using Literate for name in readdir(\"$example_dir\") project_dir = joinpath(\"$example_dir\", name) @@ -149,7 +148,6 @@ in parallel. Pkg.activate(project_dir) Pkg.instantiate() - CondaPkg.resolve() @info "building" project_dir Literate.$target( diff --git a/CondaPkg.toml b/CondaPkg.toml deleted file mode 100644 index 255866961..000000000 --- a/CondaPkg.toml +++ /dev/null @@ -1,2 +0,0 @@ -[deps] -matplotlib = "3.5.1" diff --git a/Project.toml b/Project.toml index 4beb2b54d..428673598 100644 --- a/Project.toml +++ b/Project.toml @@ -9,11 +9,11 @@ BloqadeLattices = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe4" BloqadeMIS = "bd27d05e-4ce1-5e74-84dd-c5d7d508bbe2" BloqadeODE = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe5" BloqadeWaveforms = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe7" +CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4" Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7" -PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" Yao = "5872b779-8223-5990-8dd0-5abbb0748c8c" YaoSubspaceArrayReg = "bd27d05e-4ce1-5e79-84dd-c5d7d508ade2" @@ -29,7 +29,6 @@ ColorSchemes = "3" Colors = "0.12" ForwardDiff = "0.10" Measurements = "2" -PythonCall = "0.8,0.9" Reexport = "1" Yao = "0.9" YaoSubspaceArrayReg = "0.2" diff --git a/docs/Project.toml b/docs/Project.toml index a879c23f1..2e879db8c 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -16,7 +16,6 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" -PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" Yao = "5872b779-8223-5990-8dd0-5abbb0748c8c" YaoAPI = "0843a435-28de-4971-9e8b-a9641b2983a8" diff --git a/docs/src/install.md b/docs/src/install.md index 10fa5fb18..2231bd2c0 100644 --- a/docs/src/install.md +++ b/docs/src/install.md @@ -57,8 +57,6 @@ Built on top of Ubuntu Server 20.04 LTS, this image includes - [Yao.jl](https://yaoquantum.org/) - [Revise.jl](https://github.com/timholy/Revise.jl) - [BenchmarkTools.jl](https://juliaci.github.io/BenchmarkTools.jl/stable/) -- [PythonCall.jl](https://cjdoris.github.io/PythonCall.jl/stable/) -- Conda package manager, provided by [Miniconda](https://docs.conda.io/en/latest/miniconda.html) #### Bloqade CUDA AMI @@ -227,8 +225,6 @@ Built on top of Ubuntu Server 20.04 LTS, this image includes - [Yao.jl](https://yaoquantum.org/) - [Revise.jl](https://github.com/timholy/Revise.jl) - [BenchmarkTools.jl](https://juliaci.github.io/BenchmarkTools.jl/stable/) -- [PythonCall.jl](https://cjdoris.github.io/PythonCall.jl/stable/) -- Conda package manager, provided by [Mamba](https://mamba.readthedocs.io/en/latest/index.html) - Jupyter Lab interface with dedicated Julia and Python kernels - Integrated Terminal for interactive command-line sessions diff --git a/docs/src/waveform.md b/docs/src/waveform.md index 7d3134fbd..c02c8f8db 100644 --- a/docs/src/waveform.md +++ b/docs/src/waveform.md @@ -17,14 +17,11 @@ which is a created by providing a callable object and a real number `duration`: Bloqade gives users the flexibility to specify general waveforms by inputting functions. The following code constructs a sinusoidal waveform with a time duration of 2 μs: ```@example waveform -using Bloqade -using PythonCall -plt = pyimport("matplotlib.pyplot") +using Bloqade, Bloqade.CairoMakie waveform = Waveform(t->2.2*2π*sin(2π*t), duration = 2); Bloqade.plot(waveform) ``` -In our documentation, we use the -python package [`matplotlib`](https://matplotlib.org) for plotting. +In our documentation, we use the [`CairoMakie`](https://docs.makie.org/) for plotting. Bloqade supports built-in waveforms for convenience (see References below). For example, the codes below create different waveform shapes with a single line: @@ -63,7 +60,9 @@ wf1 = piecewise_linear(;clocks, values=values1); values2 = 2π*rand(length(clocks)-1) wf2 = piecewise_constant(;clocks, values=values2); -fig, (ax1, ax2) = plt.subplots(figsize=(12, 4), ncols=2) +fig = Figure(size=(960, 400)) +ax1 = Axis(fig[1, 1]) +ax2 = Axis(fig[1, 2]) Bloqade.plot!(ax1, wf1) Bloqade.plot!(ax2, wf2) fig @@ -104,7 +103,9 @@ waveform: wf = piecewise_linear(clocks=[0.0, 2.0, 3.0, 4.0], values=2π*[0.0, 3.0, 1.1, 2.2]); swf = smooth(wf;kernel_radius=0.1); -fig, (ax1, ax2) = plt.subplots(figsize=(12, 4), ncols=2) +fig = Figure(size=(960, 400)) +ax1 = Axis(fig[1, 1]) +ax2 = Axis(fig[1, 2]) Bloqade.plot!(ax1, wf) Bloqade.plot!(ax2, swf) fig @@ -120,7 +121,9 @@ wf2 = sinusoidal(duration = 2.2, amplitude = 2.2*2π); wf3 = wf1 + wf2; wf4 = wf1 - wf2; -fig, (ax1, ax2) = plt.subplots(figsize=(12, 4), ncols=2) +fig = Figure(size=(960, 400)) +ax1 = Axis(fig[1, 1]) +ax2 = Axis(fig[1, 2]) Bloqade.plot!(ax1, wf3) Bloqade.plot!(ax2, wf4) fig @@ -133,7 +136,9 @@ To increase the strength of a waveform by some factors, we can directly use `*`: wf = linear_ramp(;duration=2.2, start_value=0.0, stop_value=1.0*2π); wf_t = 3 * wf; -fig, (ax1, ax2) = plt.subplots(figsize=(12, 4), ncols=2) +fig = Figure(size=(960, 400)) +ax1 = Axis(fig[1, 1]) +ax2 = Axis(fig[1, 2]) Bloqade.plot!(ax1, wf) Bloqade.plot!(ax2, wf_t) fig @@ -144,7 +149,9 @@ Such operations can also be broadcasted by using `.*`: ```@example waveform wf2, wf3 = [2.0, 3.0] .* wf1; -fig, (ax1, ax2) = plt.subplots(figsize=(12, 4), ncols=2) +fig = Figure(size=(960, 400)) +ax1 = Axis(fig[1, 1]) +ax2 = Axis(fig[1, 2]) Bloqade.plot!(ax1, wf2) Bloqade.plot!(ax2, wf3) fig diff --git a/examples/1.blockade/Project.toml b/examples/1.blockade/Project.toml index 6de2daf32..d0d02f70a 100644 --- a/examples/1.blockade/Project.toml +++ b/examples/1.blockade/Project.toml @@ -7,7 +7,6 @@ BloqadeMIS = "bd27d05e-4ce1-5e74-84dd-c5d7d508bbe2" BloqadeODE = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe5" BloqadeWaveforms = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe7" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" -PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" Yao = "5872b779-8223-5990-8dd0-5abbb0748c8c" YaoArrayRegister = "e600142f-9330-5003-8abb-0ebd767abc51" diff --git a/examples/1.blockade/main.jl b/examples/1.blockade/main.jl index 251b84bc5..17079e1ec 100644 --- a/examples/1.blockade/main.jl +++ b/examples/1.blockade/main.jl @@ -345,26 +345,29 @@ for _ in TimeChoiceIterator(integrator2, 0.0:dt:Tmax) end # Plot the data: -using PythonCall # Use matplotlib to generate plots -matplotlib = pyimport("matplotlib") -plt = pyimport("matplotlib.pyplot") - -ax = plt.subplot(1, 1, 1) -plt.plot(times, real(densities), "k", label = "Full space") -plt.plot(times, real(densities2), "r--", label = "Subspace") -ax.axis([0, Tmax, 0, 0.45]) -plt.xlabel("Time (us)") -plt.ylabel("Rydberg density") -plt.tight_layout() -plt.legend() - -inset_axes = pyimport("mpl_toolkits.axes_grid1.inset_locator") -ax2 = inset_axes.inset_axes(ax, width = "20%", height = "30%", loc = "lower right", borderpad = 1) -plt.plot(times, real(densities - densities2)) -plt.axis([0, 0.5, -0.001, 0.003]) -plt.ylabel("Difference", fontsize = 12) -plt.yticks(LinRange(-0.001, 0.003, 5), fontsize = 12); -plt.xticks([0, 0.2, 0.4, 0.6], fontsize = 12); +using Bloqade.CairoMakie # Use CairoMakie to generate plots + +fig = Figure() +ax = Axis(fig[1, 1], xlabel = "Time (us)", ylabel = "Rydberg density", limits=(0.0, Tmax, 0.0, 0.45)) +lines!(ax, times, real(densities), label = "Full space", color=:black) +lines!(ax, times, real(densities2), label = "Subspace", color=:red, linestyle=:dash) +axislegend(ax, position = :rt) +fig + +# The inset axis +ax2 = Axis(fig[1, 1], + width=Relative(0.2), + height=Relative(0.3), + halign=0.9, + valign=0.1, + limits=(0.0, 0.5, -0.001, 0.003), + ylabel="Difference", + yticks=LinRange(-0.001, 0.003, 5), + xticks=[0, 0.2, 0.4, 0.6], + backgroundcolor=:lightgray) + +lines!(ax2, times, real(densities - densities2), color = :black) +fig; # ![RydbergBlockadeSubspace](../../../assets/RydbergBlockadeSubspace.png) diff --git a/examples/2.adiabatic/Project.toml b/examples/2.adiabatic/Project.toml index 604de30b8..6a814e65b 100644 --- a/examples/2.adiabatic/Project.toml +++ b/examples/2.adiabatic/Project.toml @@ -7,7 +7,6 @@ BloqadeMIS = "bd27d05e-4ce1-5e74-84dd-c5d7d508bbe2" BloqadeODE = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe5" BloqadeWaveforms = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe7" KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" -PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d" YaoSubspaceArrayReg = "bd27d05e-4ce1-5e79-84dd-c5d7d508ade2" [extras] diff --git a/examples/2.adiabatic/main.jl b/examples/2.adiabatic/main.jl index 109502f71..bfaab0c4a 100644 --- a/examples/2.adiabatic/main.jl +++ b/examples/2.adiabatic/main.jl @@ -19,12 +19,10 @@ # Let's start by importing the required libraries: using Bloqade -using PythonCall +using Bloqade.CairoMakie using KrylovKit using SparseArrays -plt = pyimport("matplotlib.pyplot"); - # # Ground State Properties # We start by probing the ground state properties of the Rydberg Hamiltonian in a 1D system. @@ -59,21 +57,15 @@ end # To compare, we first plot the density profile when ``\Delta= -2π \times 10`` MHz: -fig, ax = plt.subplots(figsize = (10, 4)) -ax.bar(1:nsites, density_g[1, :]) -ax.set_xticks(1:nsites) -ax.set_xlabel("Sites") -ax.set_ylabel("Rydberg density") -ax.set_title("Density Profile: 1D Chain, Δ = -2π * 10 MHz") +fig = Figure(size=(800, 320)) +ax = Axis(fig[1, 1], xlabel = "Sites", ylabel = "Rydberg density", title = "Density Profile: 1D Chain, Δ = -2π * 10 MHz", xticks = 1:nsites) +barplot!(ax, 1:nsites, density_g[1, :]) fig # We can see that the Rydberg densities in this case are close to 0 for all sites. In contrast, for ``\Delta= 2π \times 10`` MHz, the density shows a clear ``Z_2`` ordered profile: -fig, ax = plt.subplots(figsize = (10, 4)) -ax.bar(1:nsites, density_g[30, :]) -ax.set_xticks(1:nsites) -ax.set_xlabel("Sites") -ax.set_ylabel("Rydberg density") -ax.set_title("Density Profile: 1D Chain, Δ = 2π * 10 MHz") +fig = Figure(size=(800, 320)) +ax = Axis(fig[1, 1], xlabel = "Sites", ylabel = "Rydberg density", title = "Density Profile: 1D Chain, Δ = 2π * 10 MHz", xticks = 1:nsites) +barplot!(ax, 1:nsites, density_g[30, :]) fig # More generally, we can plot an order parameter as a function of ``\Delta`` to clearly see the onset of phase transition. @@ -83,13 +75,10 @@ order_para = map(1:Δ_step) do ii return sum(density_g[ii, 1:2:nsites]) - sum(density_g[ii, 2:2:nsites]) # density on odd sites - density on even sites end -fig, ax = plt.subplots(figsize = (10, 4)) -ax.plot(Δ / 2π, order_para) -ax.set_xlabel("Δ/2π (MHz) ") -ax.set_ylabel("Order parameter") +fig = Figure(size=(800, 320)) +ax = Axis(fig[1, 1], xlabel = "Δ/2π (MHz)", ylabel = "Order parameter") +lines!(ax, Δ / (2π), order_para) fig - -# From the density profile of ground states and the change in the order parameter, we can observe a phase transition with changing ``\Delta``. # Below, we show that by slowly changing the parameters of the Hamiltonian, we can follow the trajectory of the ground states and adiabatically evolve the atoms from the ground state to the ``Z_2`` # ordered state. @@ -106,13 +95,12 @@ U2 = 2π * 10; Δ = piecewise_linear(clocks = [0.0, 0.6, 2.1, total_time], values = [U1, U1, U2, U2]); # We plot the two waveforms: -fig, (ax1, ax2) = plt.subplots(ncols = 2, figsize = (12, 4)) -Bloqade.plot!(ax1, Ω) -ax1.set_ylabel("Ω/2π (MHz)") -Bloqade.plot!(ax2, Δ) -ax2.set_ylabel("Δ/2π (MHz)") +fig = Figure(size=(960, 320)) +ax1 = Axis(fig[1, 1], ylabel = "Ω/2π (MHz)") +ax2 = Axis(fig[1, 2], ylabel = "Δ/2π (MHz)") +Bloqade.plot!(ax1, Ω / (2π)) +Bloqade.plot!(ax2, Δ / (2π)) fig - # We generate the positions for a 1D atomic chain again: nsites = 9 @@ -136,19 +124,13 @@ densities = [] for _ in TimeChoiceIterator(integrator, 0.0:1e-3:total_time) push!(densities, rydberg_density(reg)) end -D = hcat(densities...); - -# and finally plot the time-dependent dynamics of Rydberg density for each site: -fig, ax = plt.subplots(figsize = (10, 4)) -shw = ax.imshow(real(D), interpolation = "nearest", aspect = "auto", extent = [0, total_time, 0.5, nsites + 0.5]) -ax.set_xlabel("time (μs)") -ax.set_ylabel("site") -ax.set_xticks(0:0.2:total_time) -ax.set_yticks(1:nsites) -bar = fig.colorbar(shw) -fig +D = hcat(densities...)'; -# We can clearly see that a ``Z_2`` ordered state has been generated by the specified adiabatic pulse sequence. +fig = Figure(size=(800, 320)) +ax = Axis(fig[1, 1], xlabel = "time (μs)", ylabel = "site", xticks=0:0.2:total_time, yticks=(0.5:1:nsites, string.(1:nsites))) +shw = heatmap!(ax, LinRange(0.0, total_time, size(D, 1)), LinRange(0.5, nsites - 0.5, size(D, 2)), real(D)) +Colorbar(fig[1, 2], shw, label = "Rydberg density") +fig # We can also confirm it by plotting the bitstring distribution at the final time step: bitstring_hist(reg; nlargest = 20) @@ -194,11 +176,11 @@ total_time = 2.9 U = 2π * 15.0 Δ = piecewise_linear(clocks = [0.0, 0.3, 2.6, total_time], values = [-U, -U, U, U]); -fig, (ax1, ax2) = plt.subplots(ncols = 2, figsize = (10, 4)) +fig = Figure(size=(800, 320)) +ax1 = Axis(fig[1, 1], ylabel = "Ω/2π (MHz)") +ax2 = Axis(fig[1, 2], ylabel = "Δ/2π (MHz)") Bloqade.plot!(ax1, Ω) -ax1.set_ylabel("Ω/2π (MHz)") Bloqade.plot!(ax2, Δ) -ax2.set_ylabel("Δ/2π (MHz)") fig # Then, we use the waveforms and atom positions to create a Hamiltonian and define a time evolution problem: @@ -213,13 +195,10 @@ densities = []; for _ in TimeChoiceIterator(integrator, 0.0:1e-3:total_time) push!(densities, rydberg_density(reg)) end -D = hcat(densities...) - -fig, ax = plt.subplots(figsize = (10, 4)) -shw = ax.imshow(real(D), interpolation = "nearest", aspect = "auto", extent = [0, total_time, 0.5, nsites + 0.5]) -ax.set_xlabel("time (μs)") -ax.set_ylabel("site") -ax.set_xticks(0:0.2:total_time) -ax.set_yticks(1:nsites) -bar = fig.colorbar(shw) +D = hcat(densities...)' + +fig = Figure(size=(800, 320)) +ax = Axis(fig[1, 1], xlabel = "time (μs)", ylabel = "site", xticks=0:0.2:total_time, yticks=(0.5:1:nsites, string.(1:nsites))) +shw = heatmap!(ax, LinRange(0.0, total_time, size(D, 1)), LinRange(0.5, nsites - 0.5, size(D, 2)), real(D)) +Colorbar(fig[1, 2], shw, label = "Rydberg density") fig diff --git a/examples/3.quantum-scar/Project.toml b/examples/3.quantum-scar/Project.toml index 9c8f6f348..6b49bc0fa 100644 --- a/examples/3.quantum-scar/Project.toml +++ b/examples/3.quantum-scar/Project.toml @@ -1,5 +1,4 @@ [deps] -CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" Bloqade = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe1" BloqadeExpr = "bd27d05e-4ce1-5e79-84dd-c5d7d508abe2" BloqadeKrylov = "bd27d05e-4cd1-5e79-84dd-c5d7d508ade2" @@ -7,8 +6,8 @@ BloqadeLattices = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe4" BloqadeMIS = "bd27d05e-4ce1-5e74-84dd-c5d7d508bbe2" BloqadeODE = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe5" BloqadeWaveforms = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe7" +CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" -PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" Yao = "5872b779-8223-5990-8dd0-5abbb0748c8c" YaoArrayRegister = "e600142f-9330-5003-8abb-0ebd767abc51" diff --git a/examples/3.quantum-scar/main.jl b/examples/3.quantum-scar/main.jl index f172a7871..bbb0922df 100644 --- a/examples/3.quantum-scar/main.jl +++ b/examples/3.quantum-scar/main.jl @@ -20,10 +20,9 @@ # To start, we first import the required libraries: using Bloqade -using PythonCall +using Bloqade.CairoMakie using Random -plt = pyimport("matplotlib.pyplot"); # # Many-Body Rabi Oscillations with Rydberg Blockade # We first demonstrate that the strong Rydberg interactions have important effects on the Rabi oscillations of Rydberg atoms. @@ -57,11 +56,9 @@ end # Here, we use a [`KrylovEvolution`](@ref) object to simulate the dynamics for a time-independent Hamiltonian. # One can also use ODE to simulate the dynamics. For an example, see the [Adiabatic Evolution](@ref) tutorial. # The Rydberg density of this atom exhibits Rabi oscillations as a function of time, shown by the plot below: -fig, ax = plt.subplots() -ax.plot(clocks, density1[1, :]) -ax.set_xlabel("Time (μs)") -ax.set_ylabel("Single Rydberg Probability") -ax.set_title("Rabi Oscillation: Single Atom Case") +fig = Figure() +ax = Axis(fig[1,1], xlabel = "Time (μs)", ylabel = "Single Rydberg Probability", title = "Rabi Oscillation: Single Atom Case") +lines!(ax, clocks, density1[1, :], label = "1 atom") fig # For the case of 2 and 3 atoms, if they are separated far enough with negligible interactions, the total Rydberg excitation densities are simply the sum of the Rydberg density for each atom. @@ -92,14 +89,12 @@ density3 = sum(density3, dims = 1); # where ``N`` is the number of atoms. # For more information, please refer to [H. Bernien, et al. (10.1038/nature24622)](https://www.nature.com/articles/nature24622). # The total Rydberg density for the 1-, 2-, and 3-atom system is plotted below: -fig, ax = plt.subplots() -ax.plot(clocks, density1[1, :]) -ax.plot(clocks, density2[1, :]) -ax.plot(clocks, density3[1, :]) -ax.set_xlabel("Time (μs)") -ax.set_ylabel("Rydberg Probability") -ax.set_title("Many-body Rabi Oscillation for 1-, 2-, and 3-atom system") -ax.legend(["1 atom", "2 atoms", "3 atoms"], loc = "lower right") +fig = Figure() +ax = Axis(fig[1,1], xlabel = "Time (μs)", ylabel = "Rydberg Probability", title = "Many-body Rabi Oscillation for 1-, 2-, and 3-atom system") +lines!(ax, clocks, density1[1, :], label = "1 atom") +lines!(ax, clocks, density2[1, :], label = "2 atoms") +lines!(ax, clocks, density3[1, :], label = "3 atoms") +axislegend(position = :rb) fig # From this plot, we can see that the total Rydberg density for 2 (3) atom case does not exceed 1. This is because @@ -134,11 +129,13 @@ atoms = generate_sites(ChainLattice(), nsites, scale = 5.72) Ω_tot = append(Ω1, Ω2); Δ_tot = append(Δ1, Δ2); -fig, (ax1, ax2) = plt.subplots(ncols = 2, figsize = (12, 4)) +fig = Figure(size=(960, 320)) +ax1 = Axis(fig[1, 1], ylabel = "Ω/2π (MHz)") +ax2 = Axis(fig[1, 2], ylabel = "Δ/2π (MHz)") Bloqade.plot!(ax1, Ω_tot) Bloqade.plot!(ax2, Δ_tot) -ax1.set_ylabel("Ω/2π (MHz)") -ax2.set_ylabel("Δ/2π (MHz)") +fig[1, 1] = ax1 +fig[1, 2] = ax2 fig # Note that the total evolution time is 4.2 μs. @@ -171,28 +168,23 @@ end # We first plot the Rydberg density for each site as a function of time: clocks = 0:1e-3:total_time -D = hcat(densities...) - -fig, ax = plt.subplots(figsize = (10, 4)) -shw = ax.imshow(real(D), interpolation = "nearest", aspect = "auto", extent = [0, total_time, 0.5, nsites + 0.5]) -ax.set_xlabel("time (μs)") -ax.set_ylabel("site") -ax.set_xticks(0:0.4:total_time) -ax.set_yticks(1:nsites) -bar = fig.colorbar(shw) -fig +D = hcat(densities...)' +fig = Figure(size=(800, 320)) +ax = Axis(fig[1, 1]; xlabel = "Time (μs)", ylabel = "Site", xticks = 0:0.4:total_time, yticks = (0.5:1:nsites, string.(1:nsites))) +shw = heatmap!(ax, LinRange(0, total_time, size(D, 1)), LinRange(0.5, nsites - 0.5, size(D, 2)), real(D), colormap = :viridis) +bar = Colorbar(fig[1, 2], shw) +fig # We can see that the state evolves to a Neel state after the first part of the pulse (time around 2.2 μs). # After that, there are clear oscillations between the two patterns of the Rydberg density, which is a signature of the quantum scar. # We can also plot the entanglement as a function of time: -fig, ax = plt.subplots(figsize = (10, 4)) -ax.plot(clocks, entropy) -ax.set_xlabel("time (μs)") -ax.set_ylabel("entanglement entropy") +fig = Figure(size=(800, 320)) +ax = Axis(fig[1, 1]; xlabel = "time (μs)", ylabel = "entanglement entropy") +line = lines!(ax, clocks, entropy) +fig[1, 1] = ax fig - # # A Different Initial State # In order to see that the revivals depends strongly on the initial state, @@ -204,26 +196,18 @@ clocks = 0.0:1e-2:total_time; init_d = product_state(bit"100000101") prob_d = KrylovEvolution(init_d, clocks, hd) -density_mat_d = zeros(nsites, length(clocks)) +density_mat_d = zeros(length(clocks), nsites) for info in prob_d for i in 1:nsites - density_mat_d[i, info.step] = rydberg_density(info.reg, i) + density_mat_d[info.step, i] = rydberg_density(info.reg, i) end end -fig, ax = plt.subplots(figsize = (10, 4)) -shw = ax.imshow( - real(density_mat_d), - interpolation = "nearest", - aspect = "auto", - extent = [0, total_time, 0.5, nsites + 0.5], -) -ax.set_xlabel("time (μs)") -ax.set_ylabel("site") -ax.set_xticks(0:0.2:total_time) -ax.set_yticks(1:nsites) -bar = fig.colorbar(shw) +fig = Figure(size=(800, 320)) +ax = Axis(fig[1, 1]; xlabel = "time (μs)", ylabel = "site", xticks = 0:0.2:total_time, yticks = (0.5:1.0:nsites, string.(1:nsites))) +shw = heatmap!(ax, LinRange(0.0, total_time, size(density_mat_d, 1)), LinRange(0.5, nsites - 0.5, size(density_mat_d, 2)), density_mat_d, colormap = :viridis) +bar = Colorbar(fig[1, 2], shw) fig # From this figure, we see that the density does not show long-lived oscillations. diff --git a/examples/4.LGT/Project.toml b/examples/4.LGT/Project.toml index 604de30b8..6a814e65b 100644 --- a/examples/4.LGT/Project.toml +++ b/examples/4.LGT/Project.toml @@ -7,7 +7,6 @@ BloqadeMIS = "bd27d05e-4ce1-5e74-84dd-c5d7d508bbe2" BloqadeODE = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe5" BloqadeWaveforms = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe7" KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" -PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d" YaoSubspaceArrayReg = "bd27d05e-4ce1-5e79-84dd-c5d7d508ade2" [extras] diff --git a/examples/4.LGT/main.jl b/examples/4.LGT/main.jl index 227fbdac5..e9fca0173 100644 --- a/examples/4.LGT/main.jl +++ b/examples/4.LGT/main.jl @@ -13,9 +13,7 @@ # We first import the required packages using Bloqade -using PythonCall -plt = pyimport("matplotlib.pyplot"); - +using Bloqade.CairoMakie # ## Mapping between the Rydberg system and the LGT @@ -59,16 +57,13 @@ total_time = 3.5; Ω1 = piecewise_linear(clocks = [0.0, 0.2, total_time], values = [0.0, Ωmax, Ωmax]); # The waveforms for the detuning and the Rabi frequency are shown below -fig1, (ax1, ax2) = plt.subplots(ncols = 2, figsize = (12, 4)) +fig1 = Figure(size=(960, 320)) +ax1 = Axis(fig1[1, 1], ylabel="Ω1/2π (MHz)", xgridvisible=true, ygridvisible=true) +ax2 = Axis(fig1[1, 2], ylabel="Δ1/2π (MHz)", xgridvisible=true, ygridvisible=true) Bloqade.plot!(ax1, Ω1) Bloqade.plot!(ax2, Δ1) -ax1.set_ylabel("Ω1/2π (MHz)") -ax2.set_ylabel("Δ1/2π (MHz)") -ax1.grid() -ax2.grid() fig1 - # In order to simulate the gauge theory dynamics, we define the function # `get_average_rydberg_densities` which takes in a detuning # and Rabi frequency $(Δ, Ω)$, and returns the final state and the Rydberg density: @@ -93,11 +88,11 @@ end; # by simulating the dynamics governed by the waveforms, followed by plotting the density profile, as shown below: dens1 = get_average_rydberg_densities(Δ1, Ω1) ; -fig2, ax = plt.subplots(figsize = (10, 4)) ; -ax.bar(1:N, dens1[end]) ; -ax.set_xlabel("site index") -ax.set_ylabel("Average Rydberg densities") -fig2 + +fig2 = Figure(size = (800, 320)) +ax = Axis(fig2[1, 1], xlabel="site index", ylabel="Average Rydberg densities") +barplot!(ax, 1:N, dens1[end]) +fig2 # Recalling the mapping between the Rydberg chain and the LGT illustrated above, # we see that the final state is an approximation of the anti-string state of the LGT, or a $Z_2$ ordered state of the Rydberg chain. @@ -148,7 +143,10 @@ end ; # As an example, for the case with a single defect, we show the detuning for the central site, which is the defect, and those for other sites separately below: -fig3, (ax1, ax2) = plt.subplots(ncols = 2, figsize = (12, 4)) +fig3 = Figure(size=(960, 320)) +ax1 = Axis(fig3[1, 1], title="Detuning for the central site", xgridvisible=true, ygridvisible=true) +ax2 = Axis(fig3[1, 2], title="Detunings for other sites", xgridvisible=true, ygridvisible=true) + for idx in 1 : length(atoms) if idx == floor(Int, N/2)+1 Bloqade.plot!(ax1, Δ2_single_defect[idx]) @@ -156,10 +154,6 @@ for idx in 1 : length(atoms) Bloqade.plot!(ax2, Δ2_single_defect[idx]) end end -ax1.grid() -ax2.grid() -ax1.set_title("Detuning for the central site") -ax2.set_title("Detunings for other sites") fig3 @@ -169,14 +163,13 @@ fig3 dens2 = get_average_rydberg_densities(Δ2_single_defect, Ω2) dens3 = get_average_rydberg_densities(Δ2_two_defects, Ω2) -fig4, (ax1, ax2) = plt.subplots(nrows = 2, figsize = (10, 4), frameon=false) -ax1.bar(1:N, dens2[end]) -ax2.bar(1:N, dens3[end]) -fig4.supxlabel("site index") -fig4.supylabel("Average Rydberg densities", x=0.06) +fig4 = Figure(size=(800, 320)) +ax1 = Axis(fig4[1, 1], xlabel="site index", ylabel="Average Rydberg densities") +ax2 = Axis(fig4[2, 1], xlabel="site index", ylabel="Average Rydberg densities") +barplot!(ax1, 1:N, dens2[end]) +barplot!(ax2, 1:N, dens3[end]) fig4 - # Again, we see that the Rydberg density at the defects are not exactly zero, but the prepared states, # as we shall see below, serve as good initial states to study the propagation of particle-antiparticle pairs in LGT. @@ -197,20 +190,15 @@ end # Again, as an example, for the case with a single defect, we show the detuning for the central site, # which is the defect, and those for other sites separately below: -fig5, (ax1, ax2) = plt.subplots(ncols = 2, figsize = (12, 4)) -for idx in 1 : length(atoms) - if idx == floor(Int, N/2)+1 - Bloqade.plot!(ax1, Δ3_single_defect[idx]) - else - Bloqade.plot!(ax2, Δ3_single_defect[idx]) - end -end -Bloqade.plot!(ax1, Ω3) -ax1.grid() -ax2.grid() -ax1.legend(["Detuning for the central site", "Rabi frequency for all sites"]) -ax2.legend(["Detunings for other sites"], loc="center left") +fig5 = Figure(size=(960, 320)) +ax1 = Axis(fig5[1, 1], xlabel="Time", ylabel="Detuning", xgridvisible=true, ygridvisible=true) +ax2 = Axis(fig5[1, 2], xlabel="Time", ylabel="Detuning", xgridvisible=true, ygridvisible=true) +Bloqade.plot!(ax1, Δ3_single_defect[floor(Int, N/2)+1], label="Detuning for the central site") +Bloqade.plot!(ax2, Δ3_single_defect[1], label="Detunings for other sites") +Bloqade.plot!(ax1, Ω3; label="Rabi frequency for all sites") +axislegend(ax1, position=:rb) +axislegend(ax2, position=:rb) fig5 # ### Simulating Particle-Antiparticle Pairs in LGT Dynamics @@ -235,16 +223,15 @@ clocks = clocks[ind0: end]; # Then we plot the Rydberg density as a function of time, where the two panels correspond to the cases with single and two defects respectively: yticks = range(clocks[1], stop=clocks[end], length=10); yticks = [string(ytick)[1:4] for ytick in yticks][end:-1:1]; -fig6, (ax1, ax2) = plt.subplots(ncols = 2, figsize = (12, 8), sharey=true) -ax1.imshow(transpose(D_single_defect)[end:-1:1,:], aspect="auto", interpolation="none") -ax1.set_xlabel("sites") -ax1.set_ylabel("time (μs)") -ax1.set_yticks(range(1, stop = length(clocks), length=10), yticks) -im = ax2.imshow(transpose(D_two_defects)[end:-1:1,:], aspect="auto", interpolation="none") -ax2.set_xlabel("sites") -fig6.colorbar(im, ax=[ax1, ax2]) -fig6 +fig6 = Figure(size=(960, 640)) +ax1 = Axis(fig6[1, 1], xlabel = "sites", ylabel = "time (μs)", yticks = (range(1, stop = length(clocks), length = 10), string.(round.(range(clocks[1], stop=clocks[end], length=10); digits=2)))) +ax2 = Axis(fig6[1, 2], xlabel = "sites", yticklabelsvisible=false) +clims = (0.0, 1.0) +img1 = image!(ax1, D_single_defect[end:-1:1, :], colormap = :inferno, colorrange=clims, interpolate=false) +img2 = image!(ax2, D_two_defects[end:-1:1, :], colormap = :inferno, colorrange=clims, interpolate=false) +Colorbar(fig6[1, 3]; limits=clims) +fig6 # From the left panel, we can observe a light-cone-shaped region originating from the particle-antiparticle pair in the vacuum. # At the right panel, we show the interference of two light cones, which produces an additional change of periodicity corresponding to the elastic scattering. # When the particle or antiparticle reaches the boundary of the chain, it will be scattered back as observed. diff --git a/examples/5.MIS/Project.toml b/examples/5.MIS/Project.toml index c345c5027..16ee9abaa 100644 --- a/examples/5.MIS/Project.toml +++ b/examples/5.MIS/Project.toml @@ -10,7 +10,6 @@ BloqadeWaveforms = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe7" GenericTensorNetworks = "3521c873-ad32-4bb4-b63d-f4f178f42b49" Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" Optim = "429524aa-4258-5aef-a3af-852621145aeb" -PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" YaoSubspaceArrayReg = "bd27d05e-4ce1-5e79-84dd-c5d7d508ade2" diff --git a/examples/5.MIS/main.jl b/examples/5.MIS/main.jl index 418a74042..8266bb856 100644 --- a/examples/5.MIS/main.jl +++ b/examples/5.MIS/main.jl @@ -34,8 +34,7 @@ using Bloqade using Random using GenericTensorNetworks using Optim -using PythonCall -plt = pyimport("matplotlib.pyplot"); +using Bloqade.CairoMakie # ## Setting Up the Problem @@ -59,7 +58,8 @@ Bloqade.plot(atoms, blockade_radius = 7.5) # the [generic tensor network algorithm (J. Liu et al. (10.48550/arXiv.2205.03718))](https://arxiv.org/pdf/2205.03718.pdf) # in the package [`GenericTensorNetworks`](https://github.com/QuEraComputing/GenericTensorNetworks.jl): graph = BloqadeMIS.unit_disk_graph(atoms, 7.5) -mis_size_and_counting = GenericTensorNetworks.solve(IndependentSet(graph), CountingMax())[] +independent_set_problem = IndependentSet(graph) +mis_size_and_counting = GenericTensorNetworks.solve(GenericTensorNetwork(independent_set_problem), CountingMax())[] # The `solve` function takes a graph instance and a solution space property as inputs, # where the graph instance is generated by the [`unit_disk_graph`](@ref) function in the Bloqade submodule `BloqadeMIS`, @@ -80,11 +80,11 @@ T_max = 0.6 Δ_end = 2π * 11 Δ = piecewise_linear(clocks = [0.0, 0.1, 0.5, T_max], values = [Δ_start, Δ_start, Δ_end, Δ_end]) -fig, (ax1, ax2) = plt.subplots(ncols = 2, figsize = (12, 4)) +fig = Figure(size=(960, 320)) +ax1 = Axis(fig[1, 1]; ylabel="Ω/2π (MHz)") +ax2 = Axis(fig[1, 2]; ylabel="Δ/2π (MHz)") Bloqade.plot!(ax1, Ω) -ax1.set_ylabel("Ω/2π (MHz)") Bloqade.plot!(ax2, Δ) -ax2.set_ylabel("Δ/2π (MHz)") fig # Here, the total time is fixed to `T_max`, the adiabatic evolution path is specified by the [`piecewise_linear`](@ref) function. @@ -111,7 +111,7 @@ bitstring_hist(prob.reg; nlargest = 20) # The correctness of the output can be verified by comparing it to the classical solution: best_bit_strings = most_probable(prob.reg, 2) -all_optimal_configs = GenericTensorNetworks.solve(IndependentSet(graph), ConfigsMax())[] +all_optimal_configs = GenericTensorNetworks.solve(GenericTensorNetwork(independent_set_problem), ConfigsMax())[] @assert all(bs -> GenericTensorNetworks.StaticBitVector([bs...]) ∈ all_optimal_configs.c, best_bit_strings) # We can also visualize these atoms and check them visually: @@ -175,11 +175,12 @@ clocks = [0, cumsum(durations)...] Ω2 = piecewise_constant(; clocks = clocks, values = repeat([Ω_max, 0.0], 3)) Δ2 = piecewise_constant(; clocks = clocks, values = repeat([0.0, Δ_end], 3)) -fig, (ax1, ax2) = plt.subplots(ncols = 2, figsize = (12, 4)) + +fig = Figure(size=(960, 320)) +ax1 = Axis(fig[1, 1]; ylabel="Ω/2π (MHz)") +ax2 = Axis(fig[1, 2]; ylabel="Δ/2π (MHz)") Bloqade.plot!(ax1, Ω2) -ax1.set_ylabel("Ω/2π (MHz)") Bloqade.plot!(ax2, Δ2) -ax2.set_ylabel("Δ/2π (MHz)") fig # The `piecewise_constant` pulses can be more accurately simulated with the [`KrylovEvolution`](@ref) solver. @@ -231,8 +232,8 @@ end # Let us check the loss function again using the initial point above: x0 = durations -rydberg_density, reg1 = loss_piecewise_constant(atoms, x0) -rydberg_density +density, reg1 = loss_piecewise_constant(atoms, x0) +density # The most probable bitstrings are: @@ -267,11 +268,11 @@ bitstring_hist(reg1_final; nlargest = 20) pulse_piecewise_linear = piecewise_linear(clocks = [0.0, 0.05, 0.1, 0.5, 0.55, T_max], values = [0, 0, 0.4, 0.4, 0, 0]); pulse_smooth = smooth(pulse_piecewise_linear; kernel_radius = 0.02); -fig, ax = plt.subplots() -Bloqade.plot!(ax, pulse_piecewise_linear) -Bloqade.plot!(ax, pulse_smooth) -ax.set_ylabel("strength") -ax.legend(["piecewise linear", "smoothened piecewise linear"], loc = "lower right") +fig = Figure() +ax = Axis(fig[1, 1]; ylabel="strength") +Bloqade.plot!(ax, pulse_piecewise_linear; label = "piecewise linear") +Bloqade.plot!(ax, pulse_smooth; label = "smoothened piecewise linear") +axislegend(ax; position=:rb) fig # Here, the function [`smooth`](@ref) takes a `kernel_radius` keyword parameter as the Gaussian kernel parameter. @@ -317,15 +318,15 @@ T_max = 0.6 values = [Δ_start, Δ_start, Δ0 * x0[1], Δ0 * x0[2], Δ0 * x0[3], Δ_end, Δ_end], ) -rydberg_density, reg2, Δ_initial_smooth = loss_piecewise_linear(atoms, x0) -rydberg_density +density, reg2, Δ_initial_smooth = loss_piecewise_linear(atoms, x0) +density # and plot the smoothened waveform: -fig, ax = plt.subplots() -Bloqade.plot!(ax, Δ_initial) -Bloqade.plot!(ax, Δ_initial_smooth) -ax.set_ylabel("Δ/2π (MHz)") -ax.legend(["piecewise linear", "smoothened piecewise linear"], loc = "lower right") +fig = Figure() +ax = Axis(fig[1, 1]; ylabel="Δ/2π (MHz)") +Bloqade.plot!(ax, Δ_initial; label = "piecewise linear") +Bloqade.plot!(ax, Δ_initial_smooth; label = "smoothened piecewise linear") +axislegend(ax; position=:rb) fig # Let's plot the distribution: @@ -344,11 +345,11 @@ bitstring_hist(reg_final; nlargest = 20) # We can also plot out the final optimized waveform for Δ # and compare it with the initial waveform: -fig, ax = plt.subplots() -Bloqade.plot!(ax, Δ_initial_smooth) -Bloqade.plot!(ax, Δ_final) -ax.set_ylabel("Δ/2π (MHz)") -ax.legend(["initial", "optimized"], loc = "lower right") +fig = Figure() +ax = Axis(fig[1, 1]; ylabel="Δ/2π (MHz)") +Bloqade.plot!(ax, Δ_initial_smooth; label = "initial") +Bloqade.plot!(ax, Δ_final; label = "optimized") +axislegend(ax; position=:rb) fig # ## Applications in VLSI Chip Manufacturing @@ -380,35 +381,28 @@ fig using Bloqade using GenericTensorNetworks using GenericTensorNetworks: unit_disk_graph -using PythonCall +using Bloqade.CairoMakie using Random using Optim -plt = pyimport("matplotlib.pyplot"); -patches = pyimport("matplotlib.patches"); -np = pyimport("numpy") # ### Setting up the Problem # We can visualize this problem as a grid of 64K RAM chips, where some of the grid cells are removed. -fig, ax = plt.subplots(figsize=(5, 5)) - G = 5 # grid size (G x G) defects = [] for i in range(1, G*G//15) push!(defects, ((rand(Int)%G+10) % G, (rand(Int)%G + 10) % G)) end -ax.grid(visible=true) -plt.xticks(1:G) -plt.yticks(1:G) - +fig = Figure(size=(400, 400)) +ax = Axis(fig[1, 1]; ygridvisible=true, xgridvisible=true, xticks=1:G, yticks=1:G, aspect=DataAspect(), limits=(0, G, 0, G)) grid = zeros((G, G)) function plot_defects(defects) for (x, y) in defects - ax.add_patch(patches.Rectangle((x, y), 1, 1, color="red")) + CairoMakie.poly!(Rect(x, y, 1.0, 1.0), color = :red) grid[x+1, y+1] = 1 end end @@ -446,9 +440,9 @@ for x in 1:G-1 end R_opt = sqrt(2sqrt(2))/2 +circles = [] for circ in circs - patch = patches.Circle(circ, R_opt, fill=false) - ax.add_patch(patch) + push!(circles, CairoMakie.poly!(Circle(CairoMakie.Point(Float64.(circ)), R_opt), color = :transparent, strokewidth=1)) end fig @@ -476,7 +470,8 @@ Bloqade.plot(atoms, blockade_radius = R) # We will first use classical algorithms to find the optimal solution to the problem: unit_graph = unit_disk_graph(atoms, R) -configs = GenericTensorNetworks.solve(IndependentSet(unit_graph), ConfigsMax())[] +independent_set_problem = IndependentSet(unit_graph) +configs = GenericTensorNetworks.solve(GenericTensorNetwork(independent_set_problem), ConfigsMax())[] MIS_config = configs.c[1] # We can visualize the solution on the atom lattice: @@ -485,24 +480,10 @@ Bloqade.plot(atoms, blockade_radius = R; colors = [iszero(b) ? "white" : "red" f # We can also visualize the solution on the chip grid: -function clear_circ() - for i in 1:5 - for patch in ax.patches - try - patch.radius - patch.remove() - catch - print(patch) - end - end - end -end - -clear_circ() +delete!.(Ref(ax.scene), circles) for (circ, used) in zip(circs, MIS_config) if used == 1 - patch = patches.Circle(circ, R_opt, fill=false) - ax.add_patch(patch) + CairoMakie.poly!(Circle(CairoMakie.Point(Float64.(circ)), R_opt), color = :transparent, strokewidth=1) end end @@ -521,12 +502,12 @@ T_max = 0.6 Δ_end = 2π * 11 Δ = piecewise_linear(clocks = [0.0, 0.1, 0.5, T_max], values = [Δ_start, Δ_start, Δ_end, Δ_end]) -plot, (ax1, ax2) = plt.subplots(ncols = 2, figsize = (12, 4)) +fig = Figure(size=(960, 320)) +ax1 = Axis(fig[1, 1]; ylabel="Ω/2π (MHz)") +ax2 = Axis(fig[1, 2]; ylabel="Δ/2π (MHz)") Bloqade.plot!(ax1, Ω) -ax1.set_ylabel("Ω/2π (MHz)") Bloqade.plot!(ax2, Δ) -ax2.set_ylabel("Δ/2π (MHz)") -plot +fig # And construct the Hamiltonian: @@ -586,12 +567,12 @@ end function plot_waves(params::Vector{Float64}) Ω, Δ = get_waves(params) - graph, (ax1, ax2) = plt.subplots(ncols = 2, figsize = (12, 4)) + fig = Figure(size=(960, 320)) + ax1 = Axis(fig[1, 1]; ylabel="Ω/2π (MHz)") + ax2 = Axis(fig[1, 2]; ylabel="Δ/2π (MHz)") Bloqade.plot!(ax1, Ω) - ax1.set_ylabel("Ω/2π (MHz)") Bloqade.plot!(ax2, Δ) - ax2.set_ylabel("Δ/2π (MHz)") - return graph + return fig end x0 = [0.1, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8] diff --git a/examples/6.MWIS/Project.toml b/examples/6.MWIS/Project.toml index a20f21a7a..34714e734 100644 --- a/examples/6.MWIS/Project.toml +++ b/examples/6.MWIS/Project.toml @@ -11,6 +11,5 @@ GenericTensorNetworks = "3521c873-ad32-4bb4-b63d-f4f178f42b49" Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" -PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" YaoSubspaceArrayReg = "bd27d05e-4ce1-5e79-84dd-c5d7d508ade2" diff --git a/examples/6.MWIS/main.jl b/examples/6.MWIS/main.jl index 5f498fe01..d3199258f 100644 --- a/examples/6.MWIS/main.jl +++ b/examples/6.MWIS/main.jl @@ -16,9 +16,7 @@ Random.seed!(42) using Graphs using GenericTensorNetworks using GenericTensorNetworks: unit_disk_graph -using Bloqade -using PythonCall -plt = pyimport("matplotlib.pyplot"); +using Bloqade, Bloqade.CairoMakie # We now specify the atom locations and construct an example # unit disk graph on a square lattice with nearest-neighbor connections. @@ -26,7 +24,7 @@ plt = pyimport("matplotlib.pyplot"); # and all vertices closer than a distance 1.5 are connected by an edge: locs = [(1, -1), (4, 0), (1, 1), (2, 0), (0, 0), (2, 2), (2, -2), (3, 1), (3, -1)]; g = unit_disk_graph(locs, 1.5) -show_graph(g; locs = locs, vertex_colors = ["white" for i in 1:nv(g)]) +show_graph(g, map(x->60 .* x, locs); vertex_colors = ["white" for i in 1:nv(g)]) # We then assign random weights to each vertex for this example problem: weights = [rand() for i in 1:nv(g)]; @@ -34,9 +32,10 @@ weights = [rand() for i in 1:nv(g)]; # We can solve the MWIS problem classically for this graph # using the [GenericTensorNetworks](https://github.com/QuEraComputing/GenericTensorNetworks.jl) package. # The MWIS is shown in red. -configs_mapped = GenericTensorNetworks.solve(IndependentSet(g; weights = weights), ConfigsMax())[] +wmis_problem = IndependentSet(g, weights) +configs_mapped = GenericTensorNetworks.solve(GenericTensorNetwork(wmis_problem), ConfigsMax())[] MIS_config = configs_mapped.c[1] -show_graph(g; locs = locs, vertex_colors = [iszero(MIS_config[i]) ? "white" : "red" for i in 1:nv(g)]) +show_graph(g, map(x->60 .* x, locs), vertex_colors = [iszero(MIS_config[i]) ? "white" : "red" for i in 1:nv(g)]) # # Quantum Adiabatic Algorithm to Solve the MWIS Problem @@ -85,13 +84,13 @@ end t_max = 1.5 H, Ω, Δ = build_adiabatic_sweep(g, Ω_max, Δ_max, t_max, weights); -fig, (ax1, ax2) = plt.subplots(ncols = 2, figsize = (14, 4)) +fig = Figure(; size=(700, 200)) +ax1 = Axis(fig[1, 1]; ylabel = "Ω/2π (MHz)") +ax2 = Axis(fig[1, 2]; ylabel = "Δ/2π (MHz)") Bloqade.plot!(ax1, Ω) -ax1.set_ylabel("Ω/2π (MHz)") for i in 1:nv(g) Bloqade.plot!(ax2, Δ[i]) end -ax2.set_ylabel("Δ/2π (MHz)") fig # ## Compute the MIS Probability and the Adiabatic Timescale @@ -130,14 +129,12 @@ a, b = linear_fit(t_list[P_MWIS.>0.9], y) T_LZ = -1 / b; # Finally, we plot the results: -fig, (ax1, ax2) = plt.subplots(ncols = 2, figsize = (16, 6)) -ax1.scatter(t_list, P_MWIS) -ax1.set_ylabel("MWIS Probability") -ax1.set_xlabel("Time (μs)") - -ax2.scatter(t_list, broadcast(log, 1 .- P_MWIS)) -ax2.plot(t_list, a .+ b .* t_list) -ax2.set_xlabel("Time (μs)") -ax2.set_ylabel("log(1 - MWIS Probability)") +fig = Figure(; size=(800, 300)) +ax1 = Axis(fig[1, 1]; ylabel = "MWIS Probability", xlabel = "Time (μs)") +ax2 = Axis(fig[1, 2]; ylabel = "log(1 - MWIS Probability)", xlabel = "Time (μs)") + +scatter!(ax1, t_list, P_MWIS) +scatter!(ax2, t_list, broadcast(log, 1 .- P_MWIS)) +lines!(ax2, t_list, a .+ b .* t_list) fig diff --git a/examples/7.QMC/Project.toml b/examples/7.QMC/Project.toml index 199fb4b6d..717d54772 100644 --- a/examples/7.QMC/Project.toml +++ b/examples/7.QMC/Project.toml @@ -1,10 +1,16 @@ [deps] BinningAnalysis = "b7192094-8e58-5052-a244-180a858778ee" Bloqade = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe1" +BloqadeExpr = "bd27d05e-4ce1-5e79-84dd-c5d7d508abe2" +BloqadeKrylov = "bd27d05e-4cd1-5e79-84dd-c5d7d508ade2" +BloqadeLattices = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe4" +BloqadeMIS = "bd27d05e-4ce1-5e74-84dd-c5d7d508bbe2" +BloqadeODE = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe5" BloqadeQMC = "2c0eea08-e05e-11ec-33df-b7e33af02477" +BloqadeWaveforms = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe7" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7" -Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Yao = "5872b779-8223-5990-8dd0-5abbb0748c8c" +YaoSubspaceArrayReg = "bd27d05e-4ce1-5e79-84dd-c5d7d508ade2" diff --git a/examples/7.QMC/main.jl b/examples/7.QMC/main.jl index 44dbec2ef..a1431c98e 100644 --- a/examples/7.QMC/main.jl +++ b/examples/7.QMC/main.jl @@ -27,11 +27,10 @@ using BloqadeQMC using Random -using Plots # In addition, we will import some libraries that can be used later to check our QMC against ED results for small system sizes: -using Bloqade +using Bloqade, Bloqade.CairoMakie using Yao: mat, ArrayReg using LinearAlgebra using Measurements @@ -184,12 +183,9 @@ end # Let us plot the results. -using Plots: bar - - -results_plot = bar(densities_QMC, label="") -xlabel!("Site number") -ylabel!("Occupation density") +fig = Figure() +ax = Axis(fig[1, 1]; xlabel="Site number", ylabel="Occupation density") +results_plot = barplot!(ax, densities_QMC) results_plot # As expected, we see a $\mathbb{Z}_2$ pattern has emerged, just as we saw using the exact diagonalization method. @@ -302,11 +298,13 @@ for ii in 1:Δ_step energy_ED[ii] = sum(w .* energies) / sum(w) end -fig_energy = scatter(Δ/2π, energy_ED, label="ED", marker=:x); -scatter!(Δ/2π, value.(energy_QMC); yerror=uncertainty.(energy_QMC), label="QMC", marker=:x) -xlabel!("Δ/2π") -ylabel!("Energy") -fig_energy +fig = Figure() +ax = Axis(fig[1, 1]; xlabel="Δ/2π", ylabel="Energy") +scatter!(ax, Δ/2π, energy_ED, label="ED", marker=:x); +errorbars!(ax, Δ/2π, value.(energy_QMC), uncertainty.(energy_QMC)) +scatter!(ax, Δ/2π, value.(energy_QMC); label="QMC", marker=:x) +axislegend(ax, position=:rt) +fig # We see that using the QMC, we have achieved the same results as for the ED with high accuracy. @@ -345,8 +343,9 @@ for ii in 1:Δ_step append!(order_param_QMC, measurement(mean(BD), std_error(BD)*sqrt(ratio)) ) end - -fig_order = scatter(Δ/2π, value.(order_param_QMC); yerror=uncertainty.(order_param_QMC), label="", marker=:x); -xlabel!("Δ/2π (MHz)") -ylabel!("Stag mag") -fig_order +fig = Figure() +ax = Axis(fig[1, 1]; xlabel="Δ/2π", ylabel="Stag mag") +errorbars!(ax, Δ/2π, value.(order_param_QMC), uncertainty.(order_param_QMC)) +fig_order = scatter!(ax, Δ/2π, value.(order_param_QMC); label="", marker=:x) +axislegend(ax, position=:rt) +fig diff --git a/examples/8.simulating noise/Project.toml b/examples/8.simulating noise/Project.toml index 9a0ca79f7..cec82edd4 100644 --- a/examples/8.simulating noise/Project.toml +++ b/examples/8.simulating noise/Project.toml @@ -1,13 +1,19 @@ [deps] Bloqade = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe1" +BloqadeExpr = "bd27d05e-4ce1-5e79-84dd-c5d7d508abe2" +BloqadeKrylov = "bd27d05e-4cd1-5e79-84dd-c5d7d508ade2" +BloqadeLattices = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe4" +BloqadeMIS = "bd27d05e-4ce1-5e74-84dd-c5d7d508bbe2" BloqadeNoisy = "7534646f-0cda-4b9e-a311-3a9166d831c9" +BloqadeODE = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe5" +BloqadeWaveforms = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe7" CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" ProgressBars = "49802e3a-d2f1-5c88-81d8-b72133a6f568" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" Yao = "5872b779-8223-5990-8dd0-5abbb0748c8c" +YaoSubspaceArrayReg = "bd27d05e-4ce1-5e79-84dd-c5d7d508ade2" diff --git a/examples/8.simulating noise/main.jl b/examples/8.simulating noise/main.jl index e3ac75106..3c0f61917 100644 --- a/examples/8.simulating noise/main.jl +++ b/examples/8.simulating noise/main.jl @@ -40,13 +40,11 @@ using Yao using CSV using DataFrames using JSON -using LaTeXStrings using LinearAlgebra -using Plots +using Bloqade.CairoMakie using StatsBase using Printf using ProgressBars -pythonplot() # ## Noisy Single-Qubit Rabi Oscillations @@ -117,21 +115,24 @@ sol = emulate_noisy( report_error = true ) -plot( +fig = Figure() +ax = Axis(fig[1, 1]; xlabel = "Rabi periods", ylabel = L"\langle \hat{n} \rangle", title="Noisy Rabi Oscillation") +ylims!(ax, 0, 1) +x_plot = lines!(ax, save_times, sol.expectations[1], - ribbon = sol.twosigma[1], - title = "Noisy Rabi Oscillation", - xlabel = "Rabi periods", - ylabel = L"\langle \hat n \rangle", - label = "simulated" + label = "simulated", + color=:black ) -plot!( +band!(ax, save_times, sol.expectations[1] .- sol.twosigma[1], sol.expectations[1] .+ sol.twosigma[1], alpha=0.2, color=:cyan) +plot!(ax, save_times, 1/2 .- 1/2*exp.(-4γ * save_times) .* cos.(Ω * save_times), - label = "analytic" + label = "analytic", + color=:blue ) -ylims!(0,1) +axislegend(ax, position = :rt) +fig # ![BloqadeNoisy]("../../../assets/BloqadeNoisy_tutorial/noisy_depolarizingrabi.png") @@ -172,16 +173,21 @@ em = ErrorModel( prob = NoisySchrodingerProblem(zero_state(1), times, h, em) sim = emulate_noisy(prob, 2000, [mat(Z)]; report_error = true); -plot( - times, sim.expectations[1], ribbon = sim.twosigma[1], - label = "simulated", ylabel = L"\langle Z \rangle", - xlabel = "Rabi periods", title = "Coherent noise" +fig = Figure() +ax = Axis(fig[1, 1]; xlabel = "Rabi periods", ylabel = L"\langle Z \rangle", title="Coherent noise") +lines!(ax, + times, sim.expectations[1], + label = "simulated", + color=:black ) -plot!( +band!(ax, times, sim.expectations[1] .- sim.twosigma[1], sim.expectations[1] .+ sim.twosigma[1], alpha=0.2, color=:cyan) +scatter!(ax, times, (t->cos(Ω*t)*exp.(-σ^2*t^2/2 - 4γ*t)).(times), label = "analytic" ) +axislegend(ax, position = :rt) +fig # ![BloqadeNoisy]("../../../assets/BloqadeNoisy_tutorial/noisy_coherentnoise_rabi.png") @@ -214,7 +220,7 @@ plot!( # ### Experimental validation # The experimental data used to calibrate the noise model is taken from the Aquila whitepaper [3]. Readout error can be added to expectation values in the computational basis by passing `readout = true` to `emulate`. The operators must be of type `Diagonal`. Below, we show the estimation of a noisy expectation value incorporating the effect of readout error and compare to the experimental whitepaper data. -whitepaper_data = CSV.read("/Bloqade.jl/lib/BloqadeNoisy/examples/whitepaper_comparison/data/15MHz_long.csv", DataFrame, delim = ",", header = false) +whitepaper_data = CSV.read(pkgdir(Bloqade, "lib", "BloqadeNoisy", "examples", "whitepaper_comparison", "data", "15MHz_long.csv"), DataFrame; delim = ",", header = false) times = collect(whitepaper_data[1,:]) data = collect(whitepaper_data[2,:]) save_times = LinRange(0, last(times), 400) @@ -226,12 +232,15 @@ H = rydberg_h([(0.0,)]; Ω = Ω, Δ = Δ) prob = NoisySchrodingerProblem(zero_state(1), save_times, H, Aquila()) sim = emulate_noisy(prob, 1000, [mat(Op.n)]; readout_error = true); -plot(times, data, marker = :diamond, +fig = Figure() +ax = Axis(fig[1, 1]; xlabel = "Rabi periods", ylabel = L"\langle \hat{n} \rangle", title="Simulation vs Data") +scatterlines!(ax, times, data, marker = :diamond, markersize = 4, color = :green, linestyle = :dash, - xlabel = L"t \ (\mu s)", ylabel = L"\langle \hat n \rangle", - label = "Data", title = "Simulation vs Data" + label = "Data" ) -plot!(save_times, sim[1], color = :blue, label = "Sim") +plot!(ax, save_times, sim[1], color = :blue, label = "Sim") +axislegend(ax, position = :rt) +fig # ![BloqadeNoisy]("../../../assets/BloqadeNoisy_tutorial/noisy_whitepapercomparison.png") @@ -257,7 +266,7 @@ new_error_model = load_error_model(default_error_model); Ω = 2π C = 862690 * 2π Rb = (C/Ω)^(1/6) -N = 14 +N = 12 atoms = generate_sites(ChainLattice(), N, scale = .65*Rb) Δ_NNN = 1/(2N) * sum( @@ -289,21 +298,25 @@ H = rydberg_h(atoms; Ω = Ω, Δ = Δ) expt_times = LinRange(t_prep, t_prep+t_quench, 70) prob = NoisySchrodingerProblem(zero_state(N), expt_times, H, Aquila()) -@time sim = emulate_noisy(prob, 500, [mat(chain(igate(N)-put(N, i=>Op.n) for i in 1:N))]; +@time sim = emulate_noisy(prob, 10, [mat(chain(igate(N)-put(N, i=>Op.n) for i in 1:N))]; readout_error = true, report_error = true, shots = 1000, ensemble_algo = EnsembleThreads() ) -plot(times, [abs(p[1])^2 for p in ψ.u], +fig = Figure() +ax = Axis(fig[1, 1]; xlabel = L"t (\mu s)", ylabel = L"\langle (1-\hat{n})^{\otimes N}\rangle", title="Many-body scar simulation") +lines!(ax, times, [abs(p[1])^2 for p in ψ.u], label = "ideal", - color = :blue, title = "Many-body scar simulation" + color = :blue +) +errorbars!(ax, expt_times, sim.expectations[1], + sqrt.(sim.shot_error[1].^2 .+ sim.propagated_err[1].^2), + color = :red ) -scatter!(expt_times, sim.expectations, - yerr = sqrt.(sim.shot_error[1].^2 .+ sim.propagated_err[1].^2), - xlabel = L"t \ (\mu s)", - ylabel = L"\langle (1-\hat n)^{\otimes N}\rangle", +scatter!(ax, expt_times, sim.expectations[1], label = "noisy", color = :red, ) +fig # ![BloqadeNoisy]("../../../assets/BloqadeNoisy_tutorial/noisy_manybodyscar.png") @@ -330,7 +343,7 @@ tend = 10 Ω = Waveform(t-> 2.5*2π*(sin(sqrt(10)*2π*t)+cos(sqrt(2)*2π*t)), tend) R = (862690/2.5)^(1/6) -N = 12 +N = 10 D = 2^N atoms = [(R*i+.05*randn(),) for i in 1:N] @@ -348,7 +361,7 @@ function sim_depol(γ) h, depol(γ, N) ) - sim = emulate_noisy(prob, 1000, + sim = emulate_noisy(prob, 20, sol -> [ abs.(sol[end]).^2, abs(sol[end]' * ψ)^2 @@ -365,19 +378,23 @@ p = [sim_depol(γ) for γ in depol_strengths]; colors = [:blue, :green, :purple] -plot(xlabel = L"p", ylabel = L"P(p)", title = "Probability distribution vs noise strength") -stephist!(abs.(ψ).^2, yaxis =:log, normalize =true, color = :black, label = "γ = 0") -plot!(x->D*exp(-D*x), 0, .002, color = :black, linestyle = :dash, label = "") +fig = Figure() +ps = 0:1e-4:0.006 +ax = Axis(fig[1, 1]; xlabel = L"p", ylabel = L"P(p)", title = "Probability distribution vs noise strength", yscale=log10) +xlims!(ax, 0, .006) +ylims!(ax, 3, 1f4) +stephist!(normalize!(abs2.(ψ), 1), color = :black, label = "γ = 0") +lines!(ax, ps, D .* exp.(-D .* ps), color = :black, linestyle = :dash) for (p, c, r) in zip(p, colors, depol_strengths) F = exp(-r * 3N * tend) - stephist!(p[1][1], yaxis =:log, normalize =true, color = c, label = "γ = $(@sprintf("%.3f", r))") - plot!(x->D/F*exp((1-F)/F)*exp(-D* x/F), (1-F)/D, .002, color = c, linestyle = :dash, label = "") - plot!([(1-F)/D, (1-F)/D], [0, D/F], color = c, linestyle = :dash, label = "") + ps = (1-F)/D:1e-4:0.006 + stephist!(normalize(p[1][1], 1), color = c, label = "γ = $(@sprintf("%.3f", r))") + lines!(ax, ps, D ./ F .* exp((1-F)/F) .* exp.(-D .* ps ./ F), color = c, linestyle = :dash) + lines!(ax, [(1-F)/D, (1-F)/D], [0, D/F], color = c, linestyle = :dash) end - -ylims!(10, 2f4) -xlims!(0, .0015) +axislegend(ax, position = :rt) +fig # ![BloqadeNoisy]("../../../assets/BloqadeNoisy_tutorial/noisy_depol_dists.png") @@ -399,9 +416,14 @@ xlims!(0, .0015) # ``` # Substituting these in gives ``\text{XEB} = D^2(2/D^2)F + D^2(1-F)/D^2 - 1 = F``. Thus all three quantities estimate ``F``. -plot(depol_strengths, [p[i][1][2] for i in 1:3], title = "Depolarizing fidelity and cross-entropy", xlabel = L"\gamma", ylabel = "fidelity", label = L"F", yerr = [p[i][2][2] for i in 1:3]) -plot!(depol_strengths, exp.(-3N * tend .* depol_strengths), label = L"e^{-3\gamma N T}") -plot!(depol_strengths, [(2^N*sum(p[i][1][1] .* abs.(ψ).^2)-1) for i in 1:3], label = "XEB") +fig = Figure() +ax = Axis(fig[1, 1]; title = "Depolarizing fidelity and cross-entropy", xlabel = L"\gamma", ylabel = "fidelity") +errorbars!(ax, depol_strengths, [p[i][1][2] for i in 1:3], [p[i][2][2] for i in 1:3]) +scatter!(ax, depol_strengths, [p[i][1][2] for i in 1:3], label = L"F") +lines!(ax, depol_strengths, exp.(-3N * tend .* depol_strengths), label = L"e^{-3\gamma N T}") +lines!(ax, depol_strengths, [(2^N*sum(p[i][1][1] .* abs.(ψ).^2)-1) for i in 1:3], label = "XEB") +axislegend(ax, position = :rt) +fig # ![BloqadeNoisy]("../../../assets/BloqadeNoisy_tutorial/noisy_fidelitymeasures.png") @@ -409,7 +431,7 @@ plot!(depol_strengths, [(2^N*sum(p[i][1][1] .* abs.(ψ).^2)-1) for i in 1:3], la # When simulating large systems, memory can be an issue. The `NoisySchrodingerProblem` can also be used directly with the `DifferentialEquations` interface to simulate each trajectory manually if more control is required. The `randomize` function reinitializes the trajectory with a new sample from the specified distirbution of Hamiltonian parameters and chooses a random initial condition. R = (862690/2.5)^(1/6) -atoms = generate_sites(SquareLattice(), 4, 4, scale = .85*R) +atoms = generate_sites(SquareLattice(), 4, 3, scale = .85*R) N = length(atoms) tend = 3.0 Ω = 2.5*2π @@ -419,24 +441,26 @@ h = rydberg_h(atoms; Ω = Ω, Δ = 0) ψ = solve(SchrodingerProblem(zero_state(N), tend, h), DP8()).u; p = zeros(2^N) -ntraj = 250 +ntraj = 100 prob = NoisySchrodingerProblem(zero_state(N), [tend], h, Aquila()) for i in ProgressBar(1:ntraj) sample = randomize(prob) int = init(sample, DP8()) solve!(int) - p .+= abs.(int.u).^2/ntraj + p .+= abs2.(int.u)/ntraj end -scatter( +fig = Figure() +ax = Axis(fig[1, 1]; xlabel = "state", ylabel = "probability", title = "Memory-constrained simulation") +scatterlines!(ax, [real(ψ[end]' * mat(put(N, i=>Op.n)) * ψ[end]) for i in 1:N], - marker = :circle, markersize = 6, title = "Noisy vs noiseless excitation number", - xlabel = "site #", - ylabel = L"\langle \hat n_i \rangle", + marker = :circle, markersize = 6, label = "ideal" ) -scatter!([expectation_value_noisy(Aquila(), p, mat(put(N, i=>Op.n))) for i in 1:N], marker = :square, markersize = 6, label = "noisy") -ylims!(0,1) +scatterlines!(ax, [expectation_value_noisy(Aquila(), p, mat(put(N, i=>Op.n))) for i in 1:N], marker = :rect, markersize = 6, label = "noisy") +ylims!(ax, 0,1) +axislegend(ax, position = :rt) +fig # ![BloqadeNoisy]("../../../assets/BloqadeNoisy_tutorial/noisy_numberdensity.png") diff --git a/lib/BloqadeNoisy/src/problem.jl b/lib/BloqadeNoisy/src/problem.jl index 875aa4c22..e315038f7 100644 --- a/lib/BloqadeNoisy/src/problem.jl +++ b/lib/BloqadeNoisy/src/problem.jl @@ -435,9 +435,9 @@ function simulation_series_mean(sim; index=false) ntraj = length(sim) times = length(sim[:, 1]) if index == false - [mean([sim[t, i] for i in 1:ntraj]) for t in 1:times] + [mean([sim.u[i][t] for i in 1:ntraj]) for t in 1:times] else - [mean([sim[t, i][index] for i in 1:ntraj]) for t in 1:times] + [mean([sim.u[i][t][index] for i in 1:ntraj]) for t in 1:times] end end @@ -453,11 +453,11 @@ Convenience method to estimate the sampling error in the ensemble solution """ function simulation_series_err(sim; index=false, factor=2) ntraj = length(sim) - times = length(sim[:, 1]) + times = length(sim.u[1]) if index == false - [factor * std([sim[t, i] for i in 1:ntraj]) / sqrt(ntraj) for t in 1:times] + [factor * std([sim.u[i][t] for i in 1:ntraj]) / sqrt(ntraj) for t in 1:times] else - [factor * std([sim[t, i][index] for i in 1:ntraj]) / sqrt(ntraj) for t in 1:times] + [factor * std([sim.u[i][t][index] for i in 1:ntraj]) / sqrt(ntraj) for t in 1:times] end end diff --git a/src/Bloqade.jl b/src/Bloqade.jl index 8509fadd1..da92fb7f9 100644 --- a/src/Bloqade.jl +++ b/src/Bloqade.jl @@ -42,15 +42,7 @@ using Reexport export rydberg_density, rydberg_corr, bitstring_hist, bitstring_hist!, get_average_rydberg_densities -using PythonCall -const plt = PythonCall.pynew() - -function __init__() - # copied from PyPlotCall.jl - PythonCall.pycopy!(plt, pyimport("matplotlib.pyplot")) - return -end - +using CairoMakie using Colors, ColorSchemes include("plots/plots.jl") diff --git a/src/plots/bitstring_hist.jl b/src/plots/bitstring_hist.jl index 30f387756..3878c5993 100644 --- a/src/plots/bitstring_hist.jl +++ b/src/plots/bitstring_hist.jl @@ -11,27 +11,24 @@ end """ bitstring_hist!(ax, register; nlargest::Int, title="", kw...) -Plot the bitstring histgram. +Plot the bitstring histogram. # Arguments -- `ax`: the axis object from `matplotlib.pyplot`. +- `ax`: the axis object from `CairoMakie`. - `register`: the register to plot. # Keyword Arguments - `nlargest`: plot the first `nlargest` bitstrings. - `title`: title of the plot. -- `kw`: other keyword supported by `matplotlib.bar` function. +- `kw`: other keyword arguments supported by `CairoMakie.bars!` function. """ function bitstring_hist!(ax, r::ArrayReg; nlargest::Int, title = "", kw...) ps = probs(r) indices = find_largest(ps, nlargest) - obj = ax.bar(1:length(indices), ps[indices], kw...) - ax.set_xlabel("bitstring") - ax.set_ylabel("probability") - ax.set_xticks(1:length(indices), string.(indices .- 1; base = 2, pad = nqubits(r)), rotation = 60, ha = "right") + obj = barplot!(ax, 1:length(indices), ps[indices]; kw...) return obj end @@ -40,22 +37,14 @@ function bitstring_hist!(ax, r::SubspaceArrayReg; nlargest::Int, title = "", kw. probs = abs2.(state) indices = find_largest(probs, nlargest) - obj = ax.bar(1:length(indices), probs[indices]; kw...) - ax.set_xlabel("bitstring") - ax.set_ylabel("probability") - ax.set_xticks( - 1:length(indices), - string.(space(r).subspace_v[indices]; base = 2, pad = nqubits(r)), - rotation = 60, - ha = "right", - ) + obj = barplot!(ax, 1:length(indices), probs[indices]; kw...) return obj end """ bitstring_hist(r; kw...) -Plot the bitstring histgram. +Plot the bitstring histogram. # Arguments @@ -65,10 +54,24 @@ Plot the bitstring histgram. - `nlargest`: plot the first `nlargest` bitstrings. - `title`: title of the plot. -- `kw...`: other keyword supported by `matplotlib.bar` function. +- `kw...`: other keyword arguments supported by `CairoMakie.bars!` function. """ -function bitstring_hist(r; kw...) - fig, ax = plt.subplots() - bitstring_hist!(ax, r; kw...) +function bitstring_hist(r; nlargest::Int, kw...) + fig = CairoMakie.Figure() + xticks_string = _xticks_string(r, nlargest) + ax = Axis(fig[1, 1], xlabel="bitstring", ylabel="probability", xticks = (1:length(xticks_string), xticks_string), xticklabelrotation=π/3) + + bitstring_hist!(ax, r; nlargest, kw...) return fig end +function _xticks_string(r::ArrayReg, nlargest::Int) + ps = probs(r) + indices = find_largest(ps, nlargest) + return string.(indices .- 1; base = 2, pad = nqubits(r)) +end +function _xticks_string(r::SubspaceArrayReg, nlargest::Int) + state = statevec(r) + probs = abs2.(state) + indices = find_largest(probs, nlargest) + return string.(space(r).subspace_v[indices]; base = 2, pad = nqubits(r)) +end \ No newline at end of file diff --git a/src/plots/waveform.jl b/src/plots/waveform.jl index 4ec37a309..d305dea9b 100644 --- a/src/plots/waveform.jl +++ b/src/plots/waveform.jl @@ -1,13 +1,12 @@ -function plot!(ax, wf::Waveform) +function plot!(ax, wf::Waveform; kwargs...) clocks = sample_clock(wf) - fig = ax.plot(clocks, BloqadeWaveforms._rm_err.(sample_values(wf, clocks) ./ (2π));) - ax.set_xlabel("time (μs)") - ax.set_ylabel("value / 2π (MHz)") + lines!(ax, clocks, BloqadeWaveforms._rm_err.(sample_values(wf, clocks) ./ (2π)); kwargs...) return ax end -function plot(wf::Waveform) - fig, ax = plt.subplots() - plot!(ax, wf) +function plot(wf::Waveform; kwargs...) + fig = CairoMakie.Figure() + ax = Axis(fig[1, 1]; xlabel="time (μs)", ylabel = "value / 2π (MHz)") + plot!(ax, wf; kwargs...) return fig -end +end \ No newline at end of file