diff --git a/tutorials/images/steane_triangle.jpg b/tutorials/images/steane_triangle.jpg new file mode 100644 index 0000000..0a17c34 Binary files /dev/null and b/tutorials/images/steane_triangle.jpg differ diff --git a/tutorials/steane_code.ipynb b/tutorials/steane_code.ipynb new file mode 100644 index 0000000..b7e0c47 --- /dev/null +++ b/tutorials/steane_code.ipynb @@ -0,0 +1,899 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# **Exploring the Steane Code in pyQuEST**\n", + "\n", + "---\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The Steane code is a famous quantum error-correcting code that encodes a single logical qubit in seven physical qubits. Also known as a **[[7, 1, 3]]** CSS code, the Steane code is the smallest qubit CSS code capable of correcting a single error, whether a phase-flip or bit-flip.\n", + "\n", + "Quantum error correction is essential for protecting quantum information from errors due to decoherence, imperfect gate operations, and other noise sources. The Steane code is a simple example of how such errors might be handled in quantum systems, and using pyQuEST, we can probe this quantitatively. **In particular, we can use pyQuEST to determine the rate of logical errors, or how often the logical encoding will change states.**

\n", + "\n", + "#### **Objectives**\n", + "\n", + "To delve deeper into the Steane code's capabilities and practical implementations, we will achieve the following objectives:\n", + "\n", + "1. **Review the composition of the Steane code and understand its error-correcting capabilities.**\n", + "\n", + "2. **Construct the Steane code using pyQuEST for a Monte Carlo simulation.**\n", + "\n", + "3. **Construct the Steane code using pyQuEST for a density matrix simulation.**\n", + "\n", + "In our pyQuEST simulations, we will investigate the performance of the code for storing a logical state in the presence of noise.

\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## **Composition of the Steane Code**" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "At the core of the Steane code (and indeed quantum error-correction codes more generally) are a set of **stabilisers**. In the Steane code, the sets of $X$ and $Z$ stabilisers is given by:

\n", + "\n", + "$S_X = \\{X_0 X_1 X_2 X_3, X_1 X_2 X_4 X_5, X_2 X_3 X_5 X_6 \\}$\n", + "\n", + "$S_Z = \\{Z_0 Z_1 Z_2 Z_3, Z_1 Z_2 Z_4 Z_5, Z_2 Z_3 Z_5 Z_6 \\}$

\n", + "\n", + "The Steane code can be nicely represented with a triangle, with each coloured face (or plaquette) supporting an $X$ and a $Z$ stabiliser:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + " \"\n", + "
" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can see how the stabilisers above can be matched to the triangle — each of the pairs of $X$ and $Z$ stabilisers can be placed on a distinct plaquette. Using pyQuEST, we can build these stabilisers.\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### **Stabilisers in pyQuEST**" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Before we begin, we need to define a qubit register — in the case of the Steane code, we have seven data qubits. We also need an ancilla that we will recycle throughout our code." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Import the libraries for constructing a register and a circuit with relevant gates\n", + "from pyquest import Circuit, Register\n", + "from pyquest.unitaries import H, X\n", + "\n", + "\n", + "# Define an eight-qubit register - seven data qubits and one ancilla\n", + "reg = Register(8)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's now implement the $X$ stabilisers — these will help us to correct phase-flip errors. We can define a general circuit that contains the operators we need:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Set the ancilla qubit\n", + "ancilla = 7\n", + "\n", + "# Set the stabiliser qubits - e.g. for the red plaquette\n", + "qubits = [0, 1, 2, 3]\n", + "\n", + "# Define a circuit with the operators for an X stabiliser\n", + "x_stabiliser_circ = Circuit([\n", + " H(ancilla), X(qubits[0], controls=ancilla), X(qubits[1], controls=ancilla), \n", + " X(qubits[2], controls=ancilla), X(qubits[3], controls=ancilla), H(ancilla)\n", + "])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can similarly implement the $Z$ stabilisers — these will help us to correct bit-flip errors:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Define a circuit with the operators for a Z stabiliser\n", + "z_stabiliser_circ = Circuit([\n", + " X(ancilla, controls=qubits[0]), X(ancilla, controls=qubits[1]), \n", + " X(ancilla, controls=qubits[2]), X(ancilla, controls=qubits[3]),\n", + "])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "After each stabiliser circuit is run, we must include measurement operators that allow us to identify errors. \n", + "\n", + ">Note that in pyQuEST, the results of a particular measurement operator is returned as an array — this makes it easy to access the results of an operator that is applied over several qubits. If the measurement is applied as part of a circuit, results are returned as a 2D array, with each 1D array containing the results from a measurement operator in the circuit, in order. Since we plan to include the measurement operator as part of the circuit and want to measure only the ancilla qubit, the result of our measurement must be accessed using 2D indexing.\n", + "\n", + "Given that we are only using one ancilla, we must also include a 'correction' that resets the ancilla if the measurement result is -1 (indicating a 1 state). Putting everything together in functions:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Import the measurement operator\n", + "from pyquest.gates import M\n", + "\n", + "\n", + "# Define a function for implementing the X stabiliser\n", + "def x_stabiliser(ancilla, qubits):\n", + " circ = Circuit([\n", + " H(ancilla), X(qubits[0], controls=ancilla), X(qubits[1], controls=ancilla), \n", + " X(qubits[2], controls=ancilla), X(qubits[3], controls=ancilla), H(ancilla),\n", + " M([ancilla])\n", + " ])\n", + "\n", + " # Apply the circuit - the output will contain the state of the ancilla\n", + " measurement = reg.apply_circuit(circ)\n", + "\n", + " # Access the measurement result using 2D indexing\n", + " # If the ancilla is in a 1 state, reset the ancilla and return -1\n", + " if measurement[0][0] == 1:\n", + " reg.apply_operator(X(ancilla))\n", + " return -1\n", + " # Else if the ancilla is in a 0 state, just return +1\n", + " else:\n", + " return 1\n", + " \n", + "\n", + "def z_stabiliser(ancilla, qubits):\n", + " circ = Circuit([\n", + " X(ancilla, controls=qubits[0]), X(ancilla, controls=qubits[1]), \n", + " X(ancilla, controls=qubits[2]), X(ancilla, controls=qubits[3]),\n", + " M([ancilla])\n", + " ])\n", + "\n", + " measurement = reg.apply_circuit(circ)\n", + "\n", + " if measurement[0][0] == 1:\n", + " reg.apply_operator(X(ancilla))\n", + " return -1\n", + " else:\n", + " return 1" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now that we have the stabiliser framework, we can get to work on using this in a Monte Carlo simulation.

" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## **Monte Carlo Implementation**" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In a Monte Carlo approach, we aim to run a large number of trials in order to obtain useful statistical estimates of our code performance, particularly in the presence of noise. Before generating noise models, let's first explore how we can put together a simulation that can repeatedly runs our stabilisers. \n", + "\n", + "Above, we defined a single set of qubits corresponding to one of the plaquettes. We can deifne a 2D array containing all three qubit combinations:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Define the sets of qubits needed for the Steane code stabiliser measurements\n", + "chosen_qubits = [[0, 1, 2, 3], [1, 2, 4, 5], [2, 3, 5, 6]]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For our simulations, we will consider an encoded 0 logical state. We will need to prepare a perfect such state before any live runs — to do so, we will need a **decoder**. Not only will this be useful for preparing a perfect logical state, but the decoder will also be central to correcting real errors when we start to consider noise models.\n", + "\n", + "When first running our $Z$ stabilisers, we should obtain a perfect syndrome — i.e. all +1 stabiliser results. When first running our $X$ stabilisers however, the syndrome may not be perfect — each stabiliser has a 50% chance of returning +1, and a 50% chance of returning -1. We will want to perform 'corrections' that will not affect our logical state but will allow us to start with a system that always returns +1 stabiliser results in the absence of noise.\n", + "\n", + "We can build our own decoder using a lookup table, stored as a Python dictionary. If we consider all possible error syndromes for one stabiliser type (there are three results with two possibilities, so $2^3 = 8$ combinations) and the corresponding qubit error locations, we can devise the lookup table below:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Match the possible syndromes to the corresponding error locations\n", + "syndrome_lookup = {\n", + " \"[-1, 1, 1]\" : 0,\n", + " \"[-1, -1, 1]\" : 1,\n", + " \"[-1, -1, -1]\" : 2,\n", + " \"[-1, 1, -1]\" : 3,\n", + " \"[1, -1, 1]\" : 4,\n", + " \"[1, -1, -1]\" : 5,\n", + " \"[1, 1, -1]\" : 6\n", + "}" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "vscode": { + "languageId": "plaintext" + } + }, + "source": [ + "Note that the order of the stabiliser results as seen in the syndrome is given by the order in **chosen_qubits** — in the Steane triangle, this is the red plaquette, then the blue plaquette, and then finally the green plaquette. **[1, 1, 1]** is not stored as this syndrome requires no corrections. \n", + "\n", + "To see how we obtained the above, choose a data qubit and consider the effect of an error on the stabiliser results: if there was an error on data qubit 4, only the blue plaquette (or stabiliser number two) would indicate an error, leading to a -1 result for that stabiliser. If there was an error on data qubit 2 instead, all three plaquettes would indicate an error, leading to -1 results for all three stablisers.\n", + "\n", + "With this lookup table, we can define a function that takes a syndrome (as a string) and a correction operator, and performs the desired correction:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def perform_correction(results, gate):\n", + " # Check if the syndrome is in the lookup table\n", + " if results in syndrome_lookup:\n", + " # If it is in the lookup table, find the target qubit and apply the given correction operator\n", + " target_qubit = syndrome_lookup[results]\n", + " reg.apply_operator(gate(target_qubit))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "With these pieces in place, not only can we prepare a perfect logical 0 state, but we have now also laid the groundwork for the stabiliser cycles that will be the core of our simulations.

" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### **Putting it Together**" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Our Monte Carlo simulation will involve preparing a perfect logical 0 state, and running stabiliser cycles with corrections a large number of times. We can run the code below before our live runs to prepare the state and within our live loop to perform our stabiliser cycles:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Import the Z operator\n", + "from pyquest.unitaries import Z\n", + "\n", + "\n", + "# Create arrays to store X and Z stabiliser results\n", + "x_stabiliser_results = []\n", + "z_stabiliser_results = []\n", + "\n", + "# Perform all three X stabilisers, and all three Z stabilisers\n", + "# Store the measurement results in the arrays above\n", + "for i in range(3):\n", + " x_stabiliser_results.append(x_stabiliser(7, chosen_qubits[i]))\n", + "for j in range(3):\n", + " z_stabiliser_results.append(z_stabiliser(7, chosen_qubits[j]))\n", + " \n", + "# Apply phase-flips given the X stabiliser results, and bit-flips given the Z stabiliser results \n", + "perform_correction(str(x_stabiliser_results), Z)\n", + "perform_correction(str(z_stabiliser_results), X)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We will then need to construct a loop for our live runs. If, for example, we want to run 1,000,000 stabiliser cycles in our simulation, we would run:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Set the number of runs to 1,000,000\n", + "num_runs = 1000000\n", + "\n", + "for run in range(num_runs):\n", + " \n", + " x_stabiliser_results = []\n", + " z_stabiliser_results = []\n", + "\n", + " for i in range(3):\n", + " x_stabiliser_results.append(x_stabiliser(7, chosen_qubits[i]))\n", + " for j in range(3):\n", + " z_stabiliser_results.append(z_stabiliser(7, chosen_qubits[j]))\n", + "\n", + " perform_correction(str(x_stabiliser_results), Z)\n", + " perform_correction(str(z_stabiliser_results), X)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "So far, our code has run only clean stabiliser cycles: every set of stabilisers returns perfect syndromes. This does not offer any real insights into the Steane code's performance.\n", + "\n", + "By introducing noise, however, we will be able to observe how well the code does its job of preserving the correct logical state.

" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### **Exploring Noise**" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In real quantum computing systems, noise can be quite complex and involve several poorly understood degrees of freedom. For our purposes, we will only consider noise that is 'environmental,' acting on each data qubit independently and only once per cycle before all stabilisers are performed. \n", + "\n", + ">Note that we could also explore a wide range of other noise models, with pyQuEST offering a high degree of customisability.\n", + "\n", + "Such a noise model can be implemented using a random number generator, where an error on a particular data qubit is only applied if the random number is less than a given error probability, for instance 0.1%:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Import the Y operator and Python's 'random' library\n", + "from pyquest.unitaries import Y\n", + "import random\n", + "\n", + "\n", + "# Store the possible error operators\n", + "error_gates = [X, Y, Z]\n", + "\n", + "# Set the desired error probability\n", + "error_prob = 0.001\n", + "\n", + "# Loop through every data qubit\n", + "for qubit in range(7):\n", + "\n", + " # Create an array for storing the chosen errors\n", + " chosen_errors = []\n", + "\n", + " # If the random number is less than the error probability, choose an error and store it\n", + " if random.random() < error_prob:\n", + " chosen_errors.append((random.choice(error_gates))(qubit))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In order for a logical error to occur — where our encoded logical state flips from 0 to 1 — our system will need to face at least two errors at a time. This is because the Steane code can handle one error perfectly — only when there is a combination of errors disguised as a single error in the syndrome can a logical error occur. As a result, we can optimise our code so that the stabiliser cycle only runs when two or more errors occur. \n", + "\n", + "We also need to keep track of any logical errors that occur during our simulation — since performing a logical measurement will destroy quantum information, we do not want to perform any such measurements to the qubits in our simulation. Instead, in pyQuEST, we can make an exact copy of our register and perform destructive measurements on that copy, discarding it afterwards. In the Steane code, measuring along any edge of the triangle will give us the logical state of the qubit, so we can choose to measure qubits 0, 1, and 4, and check the parity of that measurement. **A logical error will be counted if the current logical state does not match the logical state of the previous stabiliser cycle.**\n", + "\n", + "Putting our overall simulation loop together:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Create a variable to store the number of logical errors\n", + "num_log_errors = 0\n", + "\n", + "# Create a variable to store the logical state of the previous cycle\n", + "# We will use this to check for logical flips in the current cycle\n", + "prev_logical = 0\n", + "\n", + "for run in range(num_runs):\n", + "\n", + " chosen_errors = []\n", + "\n", + " for qubit in range(7):\n", + "\n", + " if random.random() < error_prob:\n", + " chosen_errors.append((random.choice(error_gates))(qubit))\n", + "\n", + " # If the number of errors is less than 2, skip the run\n", + " if len(chosen_errors) < 2:\n", + " continue\n", + " # Else store the errors and continue with the stabiliser cycle\n", + " else:\n", + " for error in chosen_errors:\n", + " reg.apply_operator(error)\n", + " \n", + "\n", + " x_stabiliser_results = []\n", + " z_stabiliser_results = []\n", + "\n", + " for i in range(3):\n", + " x_stabiliser_results.append(x_stabiliser(7, chosen_qubits[i]))\n", + " for j in range(3):\n", + " z_stabiliser_results.append(z_stabiliser(7, chosen_qubits[j]))\n", + "\n", + " perform_correction(str(x_stabiliser_results), Z)\n", + " perform_correction(str(z_stabiliser_results), X)\n", + "\n", + "\n", + " # Creates a copy of the register\n", + " measure_reg = Register(copy_reg = reg)\n", + "\n", + " # Performs a destructive logical measurement\n", + " # The parity of the states measured gives the logical state\n", + " logical_result = sum(measure_reg.apply_operator(M([0, 1, 4]))) % 2\n", + "\n", + " # Checks if there has been a logical flip\n", + " if prev_logical != logical_result:\n", + " num_log_errors += 1\n", + " # Stores the new logical state for comparison in the next cycle\n", + " prev_logical = logical_result" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Finally, if we wanted to find the logical error rate, we can just divide the number of logical errors by the number of live runs:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "logical_rate = num_log_errors/num_runs\n", + "\n", + "print(\"Logical Error Rate:\", logical_rate)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We have now constructed a Monte Carlo approach to simulating the performance of the Steane code — running this over a range of physical error rates (the error probability given to the simulation) will allow you to analyse how this performance varies as the error rate increases.

" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## **Density Matrix Implementation**" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In the density matrix approach, we need to:\n", + "\n", + "- consider all possible outcomes in a stabiliser cycle;\n", + "- keep track of the probabilities of each outcome;\n", + "- and sum together the weighted density matrices for all of these outcomes.\n", + "\n", + "This means a few important differences compared to the Monte Carlo approach.

" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### **Setting Up**" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We will again need to define a register in pyQuEST, this time with the **density_matrix** parameter set to True. As before, we will also need to define functions for our stabilisers — in this approach, however, we will not include measurement operators in our circuit constructions. Using pyQuEST's built-in features, we will later make non-destructive measurements of the probabilities defining a qubit's state.\n", + "\n", + "Note that our stabiliser functions will also take **register** as an argument — this will be important for when we perform operations on copy density matrices as we consider all possible stabiliser outcomes." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from pyquest import Circuit, Register\n", + "from pyquest.unitaries import H, X\n", + "\n", + "# Define an eight-qubit register - make it a density_matrix\n", + "reg = Register(8, density_matrix = True)\n", + "\n", + "\n", + "# Define a function for implementing an X stabiliser (without measurement), taking register as an argument\n", + "def x_stabiliser(ancilla, qubits, register):\n", + " circ = Circuit([\n", + " H(ancilla), X(qubits[0], controls=ancilla), X(qubits[1], controls=ancilla), \n", + " X(qubits[2], controls=ancilla), X(qubits[3], controls=ancilla), H(ancilla)\n", + " ])\n", + "\n", + " register.apply_circuit(circ)\n", + " \n", + "\n", + "def z_stabiliser(ancilla, qubits, register):\n", + " circ = Circuit([\n", + " X(ancilla, controls=qubits[0]), X(ancilla, controls=qubits[1]), \n", + " X(ancilla, controls=qubits[2]), X(ancilla, controls=qubits[3])\n", + " ])\n", + " \n", + " register.apply_circuit(circ)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "While the lookup table we will use is similar to the Monte Carlo approach, we will need to replace the measurement outcomes, +1 and -1, with the actual measured states, 0 and 1 respectively. This will allow us to work with stabiliser outcomes using binary representations.\n", + "\n", + "Our correction function will remain the same, just taking the additional register argument as in the stabiliser functions above:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Match the possible stabiliser outcomes to the corresponding error locations\n", + "syndrome_lookup = {\n", + " \"[1, 0, 0]\" : 0,\n", + " \"[1, 1, 0]\" : 1,\n", + " \"[1, 1, 1]\" : 2,\n", + " \"[1, 0, 1]\" : 3,\n", + " \"[0, 1, 0]\" : 4,\n", + " \"[0, 1, 1]\" : 5,\n", + " \"[0, 0, 1]\" : 6\n", + "}\n", + "\n", + "\n", + "def perform_correction(results, register, gate):\n", + " # Check if the outcome set is in the lookup table\n", + " if results in syndrome_lookup:\n", + " # If it is in the lookup table, find the target qubit and apply the given correction operator\n", + " target_qubit = syndrome_lookup[results]\n", + " register.apply_operator(gate(target_qubit))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can initialise a perfect logical 0 state in a similar way to our Monte Carlo approach:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from pyquest.unitaries import Z\n", + "\n", + "\n", + "# Define the sets of qubits needed for the Steane code stabiliser measurements\n", + "chosen_qubits = [[0, 1, 2, 3], [1, 2, 4, 5], [2, 3, 5, 6]]\n", + "\n", + "# Create an array to store the stabiliser measurements\n", + "measurement = []\n", + "\n", + "\n", + "# Perform all three X stabilisers, and then all three Z stabilisers\n", + "# Store the measurement results for each in the array above\n", + "for x_round in range(0,3):\n", + " x_stabiliser(7, chosen_qubits[x_round], reg)\n", + " measurement.append(reg.apply_operator(M([7]))[0])\n", + " # Reset the ancilla if the most recent measurement was 1\n", + " if measurement[x_round] == 1:\n", + " reg.apply_operator(X(7))\n", + "\n", + "\n", + "for z_round in range(3,6):\n", + " # Adjust the indexing for chosen_qubits to cycle through all sets again\n", + " z_stabiliser(7, chosen_qubits[z_round - 3], reg)\n", + " measurement.append(reg.apply_operator(M([7]))[0])\n", + " if measurement[z_round] == 1:\n", + " reg.apply_operator(X(7))\n", + "\n", + "\n", + "# Extract the outcomes from the X stabilisers, and then from the Z stabilisers\n", + "# Pass the outcome string into the correction function\n", + "perform_correction(str(measurement[0:3]), reg, Z)\n", + "perform_correction(str(measurement[3:6]), reg, X)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "With our logical state prepared, we can now focus on setting up the core of our simulation.

" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### **Constructing the Stabiliser Cycle Loop**" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Unlike in the Monte Carlo approach, in the density matrix approach we only need to run our stabiliser cycle a small number of times in order to understand the performance of the code for a particular physical error rate. The slope when plotting the probability of measuring a logical 1 state against runs gives the logical error rate, so only a few runs — for instance, less than ten — provide sufficient information to understand the code's performance.\n", + "\n", + "For every stabiliser cycle, we need to accumulate the density matrices of all possible stabiliser cycle outcomes. We therefore need to create a blank density matrix (where all the probability amplitudes are 0) that will become the weighted sum of all possible density matrix outcomes.\n", + "\n", + "Just as in the Monte Carlo approach, we want to investigate the code under noise — we can use a similar environmental noise model, applying **depolarising** noise to each data qubit at the start of every stabiliser cycle." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Import the depolarising noise operator\n", + "from pyquest.decoherence import Depolarising\n", + "\n", + "\n", + "# Set the number of full stabiliser cycles - e.g. 5\n", + "num_runs = 5\n", + "\n", + "\n", + "for run in range(num_runs):\n", + "\n", + " # Apply depolarising noise to each data qubit\n", + " for i in range(7):\n", + " reg.apply_operator(Depolarising(i, prob=error_prob))\n", + " \n", + " # Create an accumulator to store the weighted sum of density matrices\n", + " accumulator_reg = Register(8, density_matrix = True)\n", + " accumulator_reg.init_blank_state()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Since there are six stabiliser measurements, each with two possible state measurements, there are $2^6 = 64$ different measurement combinations. We can consider all 64 using binary, taking the following steps for each combination:

\n", + "\n", + "1. **Create a copy of the current or 'master' density matrix and perform subsequent operations on this.**\n", + "2. **For each bit in the binary representation in turn, apply the corresponding stabiliser ($X$ or $Z$), update a running total with the probability of measuring the state given by the value of the bit, and force the stabiliser into 'measuring' that outcome (thereby collapsing the state).**\n", + "3. **Use the lookup table to perform corrections — just as in a normal stabiliser cycle — based on the overall measurement outcome (equivalent to the current binary representation).**\n", + "4. **Add the copy register multiplied by the running total — representing the overall probability of seeing the particular measurement outcome — and add this to the accumulator.**

\n", + "\n", + ">Note that pyQuEST allows an indirect (non-destructive) measurement of the probabilities of measuring 0 or 1 on a particular qubit with **register.prob_of_all_outcomes([qubits])**. The probability of 0 can be accessed at the 0th index, and the probability of 1 can be accessed at the 1st index.\n", + "\n", + "

Once all 64 measurement combinations have been worked through, update the master density matrix with the accumulator. The master density matrix now contains the weighted sum of all possible measurement outcomes, and will be used as the basis for the next stabiliser cycle.\n", + "\n", + "Putting these ideas together, the code within each **run** loop will look like:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Loop through every measurement combination\n", + "for k in range (64):\n", + "\n", + "\n", + " # Initialise a running total that will store the overall measurement probability\n", + " meas_prob = 1\n", + "\n", + " # Create a binary representation of the loop number, fixed at six bits (for the six stabiliser measurements)\n", + " binary_rep = [int(bit) for bit in f'{k:06b}'] \n", + "\n", + " # Make a copy of the master density matrix\n", + " # Use this for operations applied for the current measurement combination\n", + " copy_reg = Register(copy_reg = reg)\n", + "\n", + "\n", + " # Loop through the three X stabilisers\n", + " for x_qubit in range(0,3):\n", + "\n", + " # Perform the X stabiliser\n", + " x_stabiliser(7, chosen_qubits[x_qubit], copy_reg)\n", + " \n", + " # Find the probability of measuring the outcome specified by the corresponding bit in the binary number\n", + " meas_prob *= (copy_reg.prob_of_all_outcomes([7])[binary_rep[x_qubit]])\n", + "\n", + " # Force a measurement of that specified outcome\n", + " copy_reg.apply_operator(M([7], force=binary_rep[x_qubit]))\n", + "\n", + " # Reset the ancilla if the measurement was of state 1\n", + " if binary_rep[x_qubit] == 1:\n", + " copy_reg.apply_operator(X(7))\n", + "\n", + "\n", + " for z_qubit in range(3,6):\n", + " \n", + " z_stabiliser(7, chosen_qubits[z_qubit - 3], copy_reg)\n", + "\n", + " meas_prob *= (copy_reg.prob_of_all_outcomes([7])[binary_rep[z_qubit]])\n", + "\n", + " copy_reg.apply_operator(M([7], force=binary_rep[z_qubit]))\n", + "\n", + " if binary_rep[z_qubit] == 1:\n", + " copy_reg.apply_operator(X(7))\n", + "\n", + "\n", + " # Perform corrections based on the stabiliser results\n", + " perform_correction(str(binary_rep[0:3]), copy_reg, Z)\n", + " perform_correction(str(binary_rep[3:6]), copy_reg, X)\n", + "\n", + " # Add the weighted density matrix to the accumulator\n", + " accumulator_reg += meas_prob * copy_reg\n", + " \n", + "\n", + "# After all combinations, set the master density matrix to the accumulator\n", + "reg = accumulator_reg" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "After updating the register with the accumulator, we will also need to measure the probability of measuring a logical 1 state. We do not want to destructively measure this out — we can instead apply CNOTs that outputs the logical state onto the ancilla, make an indirect measurement of the state probabilities, and then reverse the circuit such that the system is left unchanged.\n", + "\n", + "The measurement function will look like:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def measure_logical(ancilla, qubits, register):\n", + "\n", + " circ = Circuit([\n", + " X(ancilla, controls=qubits[0]),\n", + " X(ancilla, controls=qubits[1]),\n", + " X(ancilla, controls=qubits[2])\n", + " ])\n", + "\n", + " # Apply the measurement circuit\n", + " register.apply_circuit(circ)\n", + " # Store the probability of measuring a 1\n", + " prob = register.prob_of_all_outcomes([ancilla])\n", + " # Apply the inverse of the measurement circuit to leave the system unchanged\n", + " register.apply_circuit(circ.inverse) \n", + "\n", + " return prob" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "All that is left to do is to add a logical check at the end of each stabiliser cycle, adding the probability of measuring a logical 1 to a storage array:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# BEFORE THE OVERALL RUNS LOOP\n", + "\n", + "# Create an array to store the probabilities of a logical 1 for every cycle\n", + "logical_prob = []\n", + "\n", + "\n", + "# ......\n", + "\n", + "\n", + "# INSIDE THE OVERALL RUNS LOOP, AT THE END\n", + "\n", + "# Use a set of edge qubits, as in the Monte Carlo carlo approach\n", + "logical_prob.append(measure_logical(7, [0, 1, 4], reg)[1])\n", + "print(\"Probability of Logical 1:\", logical_prob[run])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## **Conclusion**" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "There we have it! In pyQuEST, we have seen how to implement and analyse a fundamental quantum error-correcting code through both Monte Carlo and density matrix simulations. These methods give us valuable insights into the code's effectiveness at handling errors and preserving the integrity of logical qubits. \n", + "\n", + "See if you can implement other codes and explore their performance under various noise models!" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.6" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}