Scalable aces

README.md

@@ -1,6 +1,6 @@
 # QuantumACES.jl
 
-<!--[![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://evanhockings.github.io/QuantumACES.jl/stable/)-->
+[![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://evanhockings.github.io/QuantumACES.jl/stable/)
 [![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://evanhockings.github.io/QuantumACES.jl/dev/)
 [![Build Status](https://github.com/evanhockings/QuantumACES.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/evanhockings/QuantumACES.jl/actions/workflows/CI.yml?query=branch%3Amain)
 [![Coverage](https://codecov.io/gh/evanhockings/QuantumACES.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/evanhockings/QuantumACES.jl)
@@ -9,7 +9,7 @@
 `QuantumACES.jl` is a package for designing and simulating scalable and performant Pauli noise characterisation experiments for stabiliser circuits with averaged circuit eigenvalue sampling (ACES).
 It is particularly interested in characterising the noise associated with fault-tolerant gadgets in the context of topological quantum error correcting codes, such as surface code syndrome extraction circuits.
 
-The methods used in this package are detailed in [arXiv:2404.06545](https://arxiv.org/abs/2404.06545).<!--, and the code generating the data for this paper can be found in the `scalable_aces` folder on the [scalable_aces](https://github.com/evanhockings/QuantumACES.jl/tree/scalable_aces) branch.{:/comment}-->
+The methods used in this package are detailed in [arXiv:2404.06545](https://arxiv.org/abs/2404.06545), and the code generating the data for this paper can be found in the `scalable_aces` folder on the [scalable_aces](https://github.com/evanhockings/QuantumACES.jl/tree/scalable_aces) branch.
 These methods build on the original ACES protocol presented in [arXiv:2108.05803](https://arxiv.org/abs/2108.05803).
 This package relies on [Stim](https://github.com/quantumlib/Stim) for stabiliser circuit simulations.
 
@@ -35,7 +35,13 @@ rotated_param = get_rotated_param(dist)
 rotated_planar = get_circuit(rotated_param, dep_param)
 ```
 
-Optimise an experimental design for these parameters.
+Next, generate an experimental design for this circuit.
+
+```julia
+d = generate_design(rotated_planar)
+```
+
+Alternatively, optimise an experimental design to improve its sample efficiency.
 
 ```julia
 d = optimise_design(rotated_planar, options = OptimOptions(; ls_type = :wls, seed = seed))
@@ -94,30 +100,31 @@ fgls_z_scores_big =
 
 ## Installation and setup
 
-This is not currently a registered package, so to add it you can run
+To install this package, run the following command in the Julia REPL.
 
 ```
-julia> # press ] to enter the Pkg REPL
-
-pkg> add https://github.com/evanhockings/QuantumACES.jl
+] add QuantumACES
 ```
 
-This package relies on the Python package [Stim](https://github.com/quantumlib/Stim) to perform stabiliser simulations.
-It calls stim with [PythonCall](https://github.com/JuliaPy/PythonCall.jl), which can be a little tricky to set up.
-One helpful method for managing Python versions is [pyenv](https://github.com/pyenv/pyenv), or for Windows, [pyenv-win](https://github.com/pyenv-win/pyenv-win), which is analogous to [Juliaup](https://github.com/JuliaLang/juliaup) for Julia.
+This package relies on the Python package [Stim](https://github.com/quantumlib/Stim) to perform stabiliser circuit simulations.
+It calls Stim with [PythonCall](https://github.com/JuliaPy/PythonCall.jl).
+By default, PythonCall creates its own Python environment, but you may wish to use an existing Python installation.
+
+One helpful method for managing Python versions is [pyenv](https://github.com/pyenv/pyenv), or for Windows, [pyenv-win](https://github.com/pyenv-win/pyenv-win); these are analogous to [Juliaup](https://github.com/JuliaLang/juliaup) for Julia.
+The following assumes you are using pyenv or pyenv-win.
 
 On Windows, to instruct PythonCall to use the Python version set by pyenv, configure PythonCall's environment variables by adding the following to your `~/.julia/config/startup.jl` file
 
-```
+```julia
 ENV["JULIA_CONDAPKG_BACKEND"] = "Null"
 python_exe = readchomp(`cmd /C pyenv which python`)
 ENV["JULIA_PYTHONCALL_EXE"] = python_exe
 ```
 
-On Unix systems, shell commands are parsed directly by Julia and appear to be unaware of your PATH variable.
-I am not sure how to fix this, so you may need to manually supply `python_exe` for the Python version `<version>` as
+On Unix systems, shell commands are parsed directly by Julia and appear to be unaware of your PATH variable, and I am not sure how to work around this.
+Therefore, you may need to manually supply `python_exe` for the Python version `<version>` as
 
-```
+```julia
 python_exe = homedir() * "/.pyenv/versions/<version>/bin/python"
 ```
 

docs/example_creation.jl

@@ -8,7 +8,7 @@ struct ExampleParameters <: AbstractCircuitParameters
 end
 
 # Construct the example parameters
-function get_example_parameters(;
+function get_example_param(;
     pad_identity = true,
     single_qubit_time::Float64 = 29.0,
     two_qubit_time::Float64 = 29.0,
@@ -27,7 +27,7 @@ function get_example_parameters(;
     if pad_identity != true
         circuit_name *= "_pad_identity_$(pad_identity)"
     end
-    example_param = ExampleParameters(pad_identity, layer_time_dict, "example_circuit")
+    example_param = ExampleParameters(pad_identity, layer_time_dict, circuit_name)
     return example_param::ExampleParameters
 end
 
@@ -42,7 +42,7 @@ function example_circuit(example_param::ExampleParameters)
     qubit_num = 3
     circuit = [
         Layer([Gate("CZ", 0, [2; 3])], qubit_num),
-        Layer([Gate("CX", 0, [1; 2]), Gate("H", 0, [3])], qubit_num),
+        Layer([Gate("CZ", 0, [1; 2]), Gate("H", 0, [3])], qubit_num),
         Layer([Gate("H", 0, [1]), Gate("S", 0, [2]), Gate("H", 0, [3])], qubit_num),
     ]
     layer_types = [two_qubit_type, two_qubit_type, single_qubit_type]
@@ -171,19 +171,18 @@ r_m = m
 phen_param = get_phen_param(p, m)
 dep_param = get_dep_param(r_1, r_2, r_m)
 # Generate the circuit
-example_param = get_example_parameters()
+example_param = get_example_param()
 circuit_example = get_circuit(example_param, dep_param)
 # Optimise the experimental design
 seed = UInt(0)
 d = optimise_design(circuit_example; options = OptimOptions(; ls_type = :gls, seed = seed))
+pretty_print(d)
 # Updat the noise to the phenomenological noise model
 d_phen = update_noise(d, phen_param)
-merit_dep = calc_gls_merit(d)
-merit_phen = calc_gls_merit(d_phen)
+merit_set_dep = calc_merit_set(d)
+merit_set_phen = calc_merit_set(d_phen)
 # Simulate ACES experiments
 budget_set = [10^6; 10^7; 10^8]
 repetitions = 20
 aces_data = simulate_aces(d_phen, budget_set; repetitions = repetitions, seed = seed)
-fgls_z_scores_phen =
-    (aces_data.fgls_gate_norm_coll[:, 3] .- merit_phen.expectation) /
-    sqrt(merit_phen.variance)
+pretty_print(aces_data, merit_set_phen)

docs/src/manual.md

@@ -18,30 +18,31 @@ Typical usage of this package involves a few steps:
 
 ## Installation and setup
 
-This is not currently a registered package, so to add it you can run
+To install this package, run the following command in the Julia REPL.
 
 ```
-julia> # press ] to enter the Pkg REPL
-
-pkg> add https://github.com/evanhockings/QuantumACES.jl
+] add QuantumACES
 ```
 
-This package relies on the Python package [Stim](https://github.com/quantumlib/Stim) to perform stabiliser simulations.
-It calls stim with [PythonCall](https://github.com/JuliaPy/PythonCall.jl), which can be a little tricky to set up.
-One helpful method for managing Python versions is [pyenv](https://github.com/pyenv/pyenv), or for Windows, [pyenv-win](https://github.com/pyenv-win/pyenv-win), which is analogous to [Juliaup](https://github.com/JuliaLang/juliaup) for Julia.
+This package relies on the Python package [Stim](https://github.com/quantumlib/Stim) to perform stabiliser circuit simulations.
+It calls Stim with [PythonCall](https://github.com/JuliaPy/PythonCall.jl).
+By default, PythonCall creates its own Python environment, but you may wish to use an existing Python installation.
+
+One helpful method for managing Python versions is [pyenv](https://github.com/pyenv/pyenv), or for Windows, [pyenv-win](https://github.com/pyenv-win/pyenv-win); these are analogous to [Juliaup](https://github.com/JuliaLang/juliaup) for Julia.
+The following assumes you are using pyenv or pyenv-win.
 
 On Windows, to instruct PythonCall to use the Python version set by pyenv, configure PythonCall's environment variables by adding the following to your `~/.julia/config/startup.jl` file
 
-```
+```julia
 ENV["JULIA_CONDAPKG_BACKEND"] = "Null"
 python_exe = readchomp(`cmd /C pyenv which python`)
 ENV["JULIA_PYTHONCALL_EXE"] = python_exe
 ```
 
-On Unix systems, shell commands are parsed directly by Julia and appear to be unaware of your PATH variable.
-I am not sure how to fix this, so you may need to manually supply `python_exe` for the Python version `<version>` as
+On Unix systems, shell commands are parsed directly by Julia and appear to be unaware of your PATH variable, and I am not sure how to work around this.
+Therefore, you may need to manually supply `python_exe` for the Python version `<version>` as
 
-```
+```julia
 python_exe = homedir() * "/.pyenv/versions/<version>/bin/python"
 ```
 
@@ -74,7 +75,13 @@ rotated_param = get_rotated_param(dist)
 rotated_planar = get_circuit(rotated_param, dep_param)
 ```
 
-Optimise an experimental design for these parameters, configuring the optimisation with the parameters associated with [`OptimOptions`](@ref).
+Next, generate an experimental design for this circuit.
+
+```julia
+d = generate_design(rotated_planar)
+```
+
+Alternatively, optimise an experimental design to improve its sample efficiency, configuring the optimisation with the parameters associated with [`OptimOptions`](@ref).
 
 ```julia
 d = optimise_design(rotated_planar, options = OptimOptions(; ls_type = :wls, seed = seed))
@@ -149,7 +156,7 @@ end
 We need a function to construct the parameter struct.
 
 ```julia
-function get_example_parameters(;
+function get_example_param(;
     pad_identity = true,
     single_qubit_time::Float64 = 29.0,
     two_qubit_time::Float64 = 29.0,
@@ -168,7 +175,7 @@ function get_example_parameters(;
     if pad_identity != true
         circuit_name *= "_pad_identity_$(pad_identity)"
     end
-    example_param = ExampleParameters(pad_identity, layer_time_dict, "example_circuit")
+    example_param = ExampleParameters(pad_identity, layer_time_dict, circuit_name)
     return example_param::ExampleParameters
 end
 ```
@@ -186,7 +193,7 @@ function example_circuit(example_param::ExampleParameters)
     qubit_num = 3
     circuit = [
         Layer([Gate("CZ", 0, [2; 3])], qubit_num),
-        Layer([Gate("CX", 0, [1; 2]), Gate("H", 0, [3])], qubit_num),
+        Layer([Gate("CZ", 0, [1; 2]), Gate("H", 0, [3])], qubit_num),
         Layer([Gate("H", 0, [1]), Gate("S", 0, [2]), Gate("H", 0, [3])], qubit_num),
     ]
     layer_types = [two_qubit_type, two_qubit_type, single_qubit_type]
@@ -343,8 +350,8 @@ dep_param = get_dep_param(r_1, r_2, r_m)
 
 Then construct the circuit
 
-```
-example_param = get_example_parameters()
+```julia
+example_param = get_example_param()
 circuit_example = get_circuit(example_param, dep_param)
 ```
 
@@ -354,30 +361,26 @@ This is because the circuit acts on only three qubits and, unlike the surface co
 ```julia
 seed = UInt(0)
 d = optimise_design(circuit_example; options = OptimOptions(; ls_type = :gls, seed = seed))
+pretty_print(d)
 ```
 
 Create a copy of the optimised design that associates phenomenological noise with the circuit, and compare the predicted performance of the experimental design with depolarising and phenomenological noise.
+In particular, we can predict the expectation and mean of the normalised root-mean-square (RMS) error between the estimated and true gate eigenvalues.
 
 ```julia
 d_phen = update_noise(d, phen_param)
-merit_dep = calc_gls_merit(d)
-merit_phen = calc_gls_merit(d_phen)
+merit_set_dep = calc_merit_set(d)
+merit_set_phen = calc_merit_set(d_phen)
 ```
 
-We can also simulate the performance of the experimental design with phenomenological noise.
+We can also simulate noise characterisation experiments with this experimental design and phenomenological noise, and compare the performance to predictions by computing z-scores for the normalised RMS error with respect to the predicted expectation and variance.
+Note that the generalised least squares (GLS) estimator is the most performant and interesting here, and is implemented as an iterative feasible generalised least squares (FGLS) method.
 
 ```julia
 budget_set = [10^6; 10^7; 10^8]
 repetitions = 20
 aces_data = simulate_aces(d_phen, budget_set; repetitions = repetitions, seed = seed)
-```
-
-Finally, compare the performance to predictions at the largest measurement budget.
-
-```julia
-fgls_z_scores_phen =
-    (aces_data.fgls_gate_norm_coll[:, 3] .- merit_phen.expectation) /
-    sqrt(merit_phen.variance)
+pretty_print(aces_data, merit_set_phen)
 ```
 
 As before, note that the distribution of the normalised RMS error between the estimated and true gate eigenvalues is not quite normally distributed.

scalable_aces/Project.toml

@@ -0,0 +1,32 @@
+[deps]
+ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+GLM = "38e38edf-8417-5370-95a0-9cbb8c7f171a"
+JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
+LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Optim = "429524aa-4258-5aef-a3af-852621145aeb"
+PGFPlotsX = "8314cec4-20b6-5062-9cdb-752b83310925"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
+QuantumACES = "57b1db7a-c9fc-4178-b66e-6b8364be52bb"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd"
+
+[compat]
+ColorSchemes = "3.24"
+DataFrames = "1.6"
+GLM = "1.9"
+JLD2 = "0.4"
+LaTeXStrings = "1.3"
+LinearAlgebra = "1.10"
+Optim = "1.9"
+PGFPlotsX = "1.6"
+Plots = "1.40"
+PrettyTables = "2.3"
+QuantumACES = "0.1"
+Random = "1.10"
+StatsBase = "0.34"
+StatsPlots = "0.15"
+julia = "1.10"

scalable_aces/README.md

@@ -0,0 +1,15 @@
+# Scalable noise characterisation of syndrome extraction circuits with averaged circuit eigenvalue sampling
+
+This contains the code used to generate the results for [arXiv:2404.06545](https://arxiv.org/abs/2404.06545).
+Beware: Stim does not guarantee consistency for random seeds between versions, so do not expect to obtain exactly the same results.
+
+For these files, the prefix `rot` refers to the rotated surface code, whereas the prefix `unrot` refers to the unrotated surface code, and the Jupyter notebooks corresponding to each file plots and displays the results.
+The file endings correspond to the following functionalities:
+
+  - `optimise` files optimise designs at a range of depolarising noise strengths.
+  - `scaling` files calculate the scaling behaviour of the optimised design for depolarising and log-normal Pauli noise.
+  - `simulate` files perform many small-scale simulations for optimised and basic designs, both for depolarising noise and the seed-0 instance of log-normal Pauli noise.
+  - `simulate_big` files perform large-scale simulation for optimised and basic designs, both for depolarising noise and the seed-0 instance of log-normal Pauli noise.
+  - `runfiles` files run all of the above files.
+
+Moreover, the `toy_design.ipynb` notebook examines relative precision estimation in the context of a toy experimental design, and the `google_data.ipynb` notebook examines the error probabilities of the Google quantum device in `Suppressing quantum errors by scaling a surface code logical qubit` by Google Quantum AI (2023).