Welcome to the Docs for the Forest SDK!

The Rigetti Forest Software Development Kit includes pyQuil, the Rigetti Quil Compiler (quilc), and the Quantum Virtual Machine (qvm).

Longtime users of Rigetti Forest will notice a few changes. First, the SDK now contains a downloadable compiler and a QVM. Second, the SDK contains pyQuil 2.0, with significant updates to previous versions. As a result, programs written using previous versions of the Forest toolkit will need to be updated to pyQuil 2.0 to be compatible with the QVM or compiler.

After installing the SDK and updating pyQuil in Installation and Getting Started, see Forest 2.0: Migration Guide to get caught up on what’s new!

Quantum Cloud Services will provide users with a dedicated Quantum Machine Image, which will come prepackaged with the Forest SDK. We’re releasing a Preview to the Forest SDK now, so current users can begin migrating code (and share feedback with us early and often!). Longtime Forest users should start with the Migration Guide which outlines key changes in this SDK Preview release.

If you’re new to Forest, we hope this documentation will provide everything you need to get up and running with the toolkit. Once you’ve oriented yourself here, proceed to the section Installation and Getting Started to get started. If you’re new to quantum computing, you also go to our section on Introduction to Quantum Computing. There, you’ll learn the basic concepts needed to write quantum software. You can also work through an introduction to quantum computing in a jupyter notebook; launch the notebook from the source folder in pyquil’s docs:

cd pyquil/docs/source
jupyter notebook intro_to_qc.ipynb

A few terms to orient you as you get started with Forest:

  1. pyQuil: An open source Python library to help you write and run quantum programs. The source is hosted on github.
  2. Quil: The Quantum Instruction Language standard. Instructions written in Quil can be executed on any implementation of a quantum abstract machine, such as the quantum virtual machine (QVM), or on a real quantum processing unit (QPU). More details regarding Quil can be found in the whitepaper, A Practical Quantum Instruction Set Architecture.
  3. QVM: The Quantum Virtual Machine is an implementation of a quantum abstract machine on classical hardware. The QVM lets you use a regular computer to simulate a small quantum computer and execute Quil programs.
  4. QPU: Quantum processing unit. This refers to the physical hardware chip which we run quantum programs on.
  5. Quil Compiler: The compiler, quilc, compiles Quil written for one quantum abstract machine (QAM) to another. Our compiler will take arbitrary Quil and compile it for the given QAM, according to its supported instruction set architecture.
  6. Forest SDK: Our software development kit, optimized for near-term quantum computers that operate as coprocessors, working in concert with traditional processors to run hybrid quantum-classical algorithms. For references on problems addressable with near-term quantum computers, see Quantum Computing in the NISQ era and beyond.

Our flagship product Quantum Cloud Services offers users an on-premise, dedicated access point to our quantum computers. This access point is a fully-configured VM, which we call a Quantum Machine Image. A QMI is bundled with the same downloadable SDK mentioned above, and a command line interface (CLI), which is used for scheduling compute time on our quantum computers. To sign up for our waitlist, please click the link above. If you’d like to access to our quantum computers for research, please email support@rigetti.com.

Note

To join our user community, connect to the Rigetti Slack workspace at https://rigetti-forest.slack.com.

Contents

Installation and Getting Started

To make full use of the Rigetti Forest SDK, you will need pyQuil, the QVM, and the Quil Compiler. On this page, we will take you through the process of installing all three of these. We also step you through running a basic pyQuil program.

Note

If you’re running from a Quantum Machine Image, installation has been completed for you. Continue to Getting Started.

Upgrading or Installing pyQuil

PyQuil 2.0 is our library for generating and executing Quil programs on the Rigetti Forest platform.

Before you install, we recommend that you activate a Python 3.6+ virtual environment. Then, install pyQuil using pip:

pip install pyquil

For those of you that already have pyQuil, you can upgrade with:

pip install --upgrade pyquil

If you would like to stay up to date with the latest changes and bug fixes, you can also opt to install pyQuil from the source here.

Note

PyQuil requires Python 3.6 or later.

Downloading the QVM and Compiler

The Forest 2.0 Downloadable SDK Preview currently contains:

  • The Rigetti Quantum Virtual Machine (qvm) which allows high-performance simulation of Quil programs
  • The Rigetti Quil Compiler (quilc) which allows compilation and optimization of Quil programs to native gate sets

The QVM and the compiler are packed as program binaries that are accessed through the command line. Both of them provide support for direct command-line interaction, as well as a server mode. The server mode is required for use with pyQuil.

Request the Forest SDK here. You’ll receive an email right away with the download links for macOS, Linux (.deb), Linux (.rpm), and Linux (bare-bones).

All installation mechanisms, except the bare-bones package, require administrative privileges to install. To use the QVM and Quil Compiler from the bare-bones package, you will have to install the prerequisite dependencies on your own.

Installing on macOS

Mount the file forest-sdk.dmg by double clicking on it in your email. From there, open forest-sdk.pkg by double-clicking on it. Follow the installation instructions.

Upon successful installation, one should be able to open a new terminal window and run the following two commands:

qvm --version
quilc --version

To uninstall, delete the following files:

/usr/local/bin/qvm
/usr/local/bin/quilc
/usr/local/share/man/man1/qvm.1
/usr/local/share/man/man1/quilc.1
Installing the QVM and Compiler on Linux (deb)

Download the Debian distribution by clicking on the link in your email. Unpack the tarball and change to that directory by doing:

tar -xf forest-sdk-linux-deb.tar.bz2
cd forest-sdk-2.0rc2-linux-deb

From here, run the following command:

sudo ./forest-sdk-2.0rc2-linux-deb.run

Upon successful installation, one should be able to run the following two commands:

qvm --version
quilc --version

To uninstall, type:

sudo apt remove forest-sdk
Installing the QVM and Compiler on Linux (rpm)

Download the RPM-based distribution by clicking on the link in your email. Unpack the tarball and change to that directory by doing:

tar -xf forest-sdk-linux-rpm.tar.bz2
cd forest-sdk-2.0rc2-linux-rpm

From here, run the following command:

sudo ./forest-sdk-2.0rc2-linux-rpm.run

Upon successful installation, one should be able to run the following two commands:

qvm --version
quilc --version

To uninstall, type:

sudo rpm -e forest-sdk
# or
sudo yum uninstall forest-sdk
Installing the QVM and Compiler on Linux (bare-bones)

The bare-bones installation only contains the executable binaries and manual pages, and doesn’t contain any of the requisite dynamic libraries. As such, installation doesn’t require administrative or sudo privileges.

First, unpack the tarball and change to that directory by doing:

tar -xf forest-sdk-linux-barebones.tar.bz2
cd forest-sdk-2.0rc2-linux-barebones

From here, run the following command:

./forest-sdk-2.0rc2-linux-barebones.run

Upon successful installation, this will have created a new directory rigetti in your home directory that contains all of the binary and documentation artifacts.

This method of installation requires one, through whatever means, to install shared libraries for BLAS, LAPACK, and libffi. On a Debian-derivative system, this could be accomplished with

sudo apt-get install liblapack-dev libblas-dev libffi-dev

To uninstall, remove the directory ~/rigetti.

Getting Started

To get started using the SDK, you can either interact with the QVM and the compiler directly from the command line, or you can run them in server mode and use them with pyQuil. In this section, we’re going to explain how to do the latter.

For more information about directly interacting with the QVM and the compiler, refer to their respective manual pages. After installation, you can read the manual pages by opening a new terminal window and typing man qvm (for the QVM) or man quilc (for the compiler). Quit out of the manual page by typing q.

Setting Up Server Mode for PyQuil

Note

This set up is only necessary to run pyQuil locally. If you’re running in a QMI, this has already been done for you.

It’s easy to start up local servers for the QVM and quilc on your laptop. You should have two terminal windows open to run in the background. We recommend using a resource such as tmux for running and managing multiple programs in one terminal.

### CONSOLE 1
$ qvm -S

Welcome to the Rigetti QVM
(Configured with 10240 MiB of workspace and 8 workers.)
[2018-09-20 15:39:50] Starting server on port 5000.


### CONSOLE 2
$ quilc -S

Welcome to the Rigetti Quil Compiler
[2018-09-19 11:22:37] Starting server: 0.0.0.0 : 6000.

That’s it! You’re all set up to run pyQuil locally. Your programs will make requests to these server endpoints to compile your Quil programs to native Quil, and to simulate those programs on the QVM.

Run Your First Program

Now that our local endpoints are up and running, we can start running pyQuil programs! We will run a simple program on the Quantum Virtual Machine (QVM).

The program we will create prepares a fully entangled state between two qubits, called a Bell State. This state is in an equal superposition between |00⟩ and |11⟩, meaning that it is equally likely that a measurement will result in measuring both qubits in the ground state or both qubits in the excited state. For more details about the physics behind these concepts, see Introduction to Quantum Computing.

To begin, start up python however you like. You can open a jupyter notebook (type jupyter notebook in your terminal), open an interactive python notebook in your terminal (with ipython3), or simply launch python in your terminal (type python3). Recall that you need Python 3.6+ to use pyQuil.

Import a few things from pyQuil:

from pyquil import Program, get_qc
from pyquil.gates import *

The Program object allows us to build up a Quil program. get_qc() connects us to a QuantumComputer object, which specifies what our program should run on (see: The Quantum Computer). We’ve also imported all (*) gates from the pyquil.gates module, which allows us to add operations to our program (Programs and Gates).

Next, let’s construct our Bell State.

# construct a Bell State program
p = Program(H(0), CNOT(0, 1))

We’ve accomplished this by driving qubit 0 into a superposition state (that’s what the “H” gate does), and then creating an entangled state between qubits 0 and 1 (that’s what the “CNOT” gate does). Finally, we’ll want to run our program:

# run the program on a QVM
qc = get_qc('9q-square-qvm')
result = qc.run_and_measure(p, trials=10)
print(result[0])
print(result[1])

Compare the two arrays of measurement results. The results will be correlated between the qubits and random from shot to shot.

The qc is a simulated quantum computer. By specifying we want to .run_and_measure, we’ve told our QVM to run the program specified above, collapse the state with a measurement, and return the results to us. trials refers to the number of times we run the whole program.

The call to run_and_measure will make a request to the two servers we started up in the previous section: first, to the quilc server instance to compile the Quil program into native Quil, and then to the qvm server instance to simulate and return measurement results of the program 10 times. If you open up the terminal windows where your servers are running, you should see output printed to the console regarding the requests you just made.

In the following sections, we’ll cover gates, program construction & execution, and go into detail about our Quantum Virtual Machine, our QPUs, noise models and more. If you’ve used pyQuil before, continue on to our Forest 2.0: Migration Guide. Once you’re set with that, jump to Programs and Gates to continue.

Forest 2.0: Migration Guide

The goals of this guide are to cover changes to the Forest SDK (containing pyquil 2.0, new Quil, Quil Compiler, and QVM), and to go through an example of migrating a VQE program from Forest 1.3 (pyQuil 1.9, Quil 1.0) to be able to run on the new Forest SDK.

Note

For installation & setup, follow the download instructions in the section Installation and Getting Started at the top of the page.

What’s changed

With the new Forest SDK, users will be able to run pyQuil programs on a downloadable QVM and Quil Compiler!

In the following section, we’ll cover the main changes to pyQuil, Quil, the Quil Compiler, and the QVM.

Overview of Updates to Quil and pyQuil

The primary differences in the programming language Quil 1.0 (as appearing in pyQuil 1.3) and Quil 2 (as appearing in 2.0) amount to an enhanced memory model. Whereas the classical memory model in Quil 1.0 amounted to an flat bit array of indefinite size, the memory model in Quil 2 is segmented into typed, sized, named regions.

In terms of compatibility with Quil 1.0, this primarily changes how MEASURE instructions are formulated, since their classical address targets must be modified to fit the new framework. In terms of new functionality, this allows angle values to be read in from classical memory.

Quil 2 also introduces easier ways to manipulate gates by using gate modifiers. Two gate modifiers are supported currently, DAGGER and CONTROLLED.

DAGGER can be written before a gate to refer to its inverse. For instance

DAGGER RX(pi/3) 0

would have the same effect as

RX(-pi/3) 0

DAGGER can be applied to any gate, but also circuits defined with DEFCIRCUIT. This allows for easy reversal of unitary circuits:

DEFCIRCUIT BELL:
    H 0
    CNOT 0 1

# construct a Bell state
BELL
# disentangle, bringing us back to identity
DAGGER BELL

Parametric programs

The main benefit for users of declared memory regions in Quil is that angle values for parametric gates can be loaded at execution time on the QPU. Consider the following simple QAOA instance:

DECLARE ro BIT[2]
DECLARE beta REAL
DECLARE gamma REAL

H 0
RZ(beta) 0
H 0
H 1
RZ(beta) 1
H 1

CNOT 0 1
RZ(gamma) 1
CNOT 0 1

MEASURE 0 ro[0]
MEASURE 1 ro[1]

To generate a “landscape” plot as beta and gamma range, it was previously required to generate a different program for each possible pair of values, substitute that pair in, send it to the compiler, and send the resulting compiled program to the QPU for execution (and hence generate the expectation values). With Quil 2, this exact program can be sent to the compiler, which returns a nativized Quil program that still has parametric gates with parameters referencing the classical memory regions beta and gamma. This program can then be loaded onto the QPU for repeated execution with different values of beta and gamma, without recompilation in between.

Details of updates to Quil

Classical memory regions must be explicitly requested and named by a Quil program using DECLARE directive. A generic DECLARE directive has the following syntax:

DECLARE region-name type([count])? (SHARING parent-region-name (OFFSET (offset-count offset-type)+))?

The non-keyword items have the following allowable values:

  • region-name: any non-keyword formal name.
  • type: one of REAL, BIT, OCTET, or INTEGER
  • parent-region-name: any non-keyword formal name previously used as region-name in a different DECLARE statement.
  • offset-count: a nonnegative integer.
  • offset-type: the same allowable values as type.

Here are some examples:

DECLARE beta REAL[32]
DECLARE ro BIT[128]
DECLARE beta-bits BIT[1436] SHARING beta
DECLARE fourth-bit-in-beta1 BIT SHARING beta OFFSET 1 REAL 4 BIT

In order, the intention of these DECLARE statements is:

  • Allocate an array called beta of length 32, each entry of which is a REAL number.
  • Allocate an array called ro of length 128, each entry of which is a BIT.
  • Name an array called beta-bits, which is an overlay onto the existing array beta, so that the bit representations of elements of beta can be directly examined and manipulated.
  • Name a single BIT called fourth-bit-in-beta1 which overlays the fourth bit of the bit representation of the REAL value beta[1].

Backwards compatibility

Quil 1.0 is not compatible with Quil 2 in the following ways:

  • The unnamed memory references [n] and [n-m] have no direct equivalent in Quil 2 and must be replaced by named memory references. (This primarily affects MEASURE instructions.)
  • The classical memory manipulation instructions have been modified: the operands of AND have been reversed (so that in Quil 2, the left operand is the target address) and OR has been replaced by IOR and its operands reversed (so that, again, in Quil 2 the left operand is the target address).

In all other instances, Quil 1.0 will operate identically with Quil 2.

When confronted with program text conforming to Quil 1.0, pyQuil 2.0 will automatically rewrite MEASURE q [n] to MEASURE q ro[n] and insert a DECLARE statement which allocates a BIT-array of the appropriate size named ro.

Details of pyQuil and Forest updates

Updates to Forest

  • In Forest 1.3, job submission to the QPU was done from your workstation and the ability was gated by on user ID. In Forest 2.0, job submission to the QPU must be done from your remote virtual machine, called a QMI (Quantum Machine Image).
  • In Forest 1.3, user data persisted indefinitely in cloud storage and could be accessed using the assigned job ID. In Forest 2.0, user data is stored only transiently, and it is the user’s responsibility to handle long-term data storage on their QMI.

Updates to pyQuil

  • In pyQuil 1.9, API calls were organized by endpoint (e.g., all simulation calls were passed to a QVMConnection object). In pyQuil 2.0, API calls are organized by type (e.g., run calls are sent to a QuantumComputer but wavefunction calls are sent to a WavefunctionSimulator).
  • In pyQuil 1.9, quantum program evaluation was asynchronous on the QPU and a mix of synchronuous or asynchronous on the QVM. In pyQuil 2.0, all quantum program evaluation is synchronous.
  • In pyQuil 1.9, each quantum program execution call started from scratch. In pyQuil 2.0, compiled program objects can be reused.

Backwards compatibility and migration

PyQuil 2.0 is not backwards compatible with pyQuil 1.9. However, the new API objects available in pyQuil 2.0 have compatibility methods that make migration to pyQuil 2.0 easier.

Note

Users writing new programs from scratch are encouraged to use the bare pyQuil 2.0 programming model over the compatibility methods. It is not possible to use the fanciest new features of Forest 2.0 (e.g., parametric execution of parametric programs) from within the compatibility model.

Whereas pyQuil 1.9 organized API calls around “connection objects” (viz., CompilerConnection, QPUConnection, and QVMConnection), pyQuil 2.0 organizes API calls around function, so that QVM- and QPU-based objects can be more easily swapped. These API objects fall into two groups:

  • QuantumComputer: This wrapper object houses the typical ingredients for execution of a hybrid classical-quantum algorithm: an interface to a compiler, an interface to a quantum computational device, and some optional wrapper routines. QuantumComputer objects themselves can be manually initialized with these ingredients, or they can be requested by name from the Forest 2.0 service, which will populate these subfields with the appropriate objects for execution on a particular quantum device, real or simulated.
  • AbstractCompiler: An interface to a compiler service. Compilers are responsible for two tasks: converting arbitrary Quil to “native” (or “device-specific”) Quil, and converting native Quil to control system binaries.
  • QAM: An interface to a quantum computational device. This can be populated by a connection to an actual QPU, or it can be populated by a connection to a QVM (Quantum Virtual Machine).
  • Wrapper routines: Execution of programs in pyQuil 1.9 was typically done with a single API call (e.g., .run()). QuantumComputer exposes a near-identical interface for single runs of quantum programs, which wraps and hides the more low-level pyQuil 2.0 infrastructure.
  • WavefunctionSimulator: This wrapper object houses the typical ingredients used for the debug process of wavefunction inspection. This is inherently not a procedure natively available on a quantum computational device, and so this wrapper either calls out to a QVM or functions as a repeated sampling wrapper from a physical quantum computational device.
Example: Computing the bond energy of molecular hydrogen, pyQuil 1.9 vs 2.0

By way of example, let’s consider the following pyQuil 1.9 program, which computes the natural bond distance in molecular hydrogen using a VQE-type algorithm:

from pyquil.api import QVMConnection
from pyquil.quil import Program


def setup_forest_objects():
    qvm = QVMConnection()
    return qvm


def build_wf_ansatz_prep(theta):
    program = Program(f"""
# set up initial state
X 0
X 1

# build the exponentiated operator
RX(pi/2) 0
H 1
H 2
H 3

CNOT 0 1
CNOT 1 2
CNOT 2 3
RZ({theta}) 3
CNOT 2 3
CNOT 1 2
CNOT 0 1

RX(-pi/2) 0
H 1
H 2
H 3

# measure out the results
MEASURE 0 [0]
MEASURE 1 [1]
MEASURE 2 [2]
MEASURE 3 [3]""")
    return program


# some constants
bond_step, bond_min, bond_max = 0.05, 0, 200
angle_step, angle_min, angle_max = 0.1, 0, 63
convolution_coefficients = [0.1698845197777728, 0.16988451977777283, -0.2188630663199042,
                            -0.2188630663199042]
shots = 1000

# set up the Forest object
qvm = setup_forest_objects()

# get all the unweighted expectations for all the sample wavefunctions
occupations = list(range(angle_min, angle_max))
indices = list(range(4))
for offset in occupations:
    # set up the Program object, each time we have a new parameter
    program = build_wf_ansatz_prep(angle_min + offset * angle_step)
    bitstrings = qvm.run(program, indices, trials=shots)

    totals = [0, 0, 0, 0]
    for bitstring in bitstrings:
        for index in indices:
            totals[index] += bitstring[index]
    occupations[offset] = [t / shots for t in totals]

# compute minimum energy as a function of bond length
min_energies = list(range(bond_min, bond_max))
for bond_length in min_energies:
    energies = []
    for offset in range(angle_min, angle_max):
        energy = 0
        for j in range(4):
            energy += occupations[offset][j] * convolution_coefficients[j]
        energies.append(energy)

    min_energies[bond_length] = min(energies)

min_index = min_energies.index(min(min_energies))
min_energy, relaxed_length = min_energies[min_index], min_index * bond_step

In order to port this code to pyQuil 2.0, we need change only one thing: the part referencing QVMConnection should be replaced by an equivalent part referencing a QuantumComputer connected to a QVM. Specifically, the following snippet

from pyquil.api import QVMConnection

def setup_forest_objects():
    qvm = QVMConnection()
    return qvm

can be changed to

from pyquil.api import get_qc

def setup_forest_objects():
    qc = get_qc("9q-square-qvm")
    return qc

and the references to qvm in the main body are changed to qc instead. Since the QuantumComputer object also exposes a run routine and pyQuil itself automatically rewrites 1.9-style MEASURE instructions into 2.0-style instructions, this is all we need to do.

If we are willing to be more intrusive, we can also take advantage of pyQuil 2.0’s classical memory and parametric programs. The first piece to change is the Quil program itself: we remove the argument theta from the Python function build_wf_ansatz_prep, with the intention of letting the QPU fill it in later. In turn, we modify the Quil program itself to have a REAL memory parameter named theta. We also declare a few BITs for our MEASURE instructions to target.

def build_wf_ansatz_prep():
    program = Program("""
# set up memory
DECLARE ro BIT[4]
DECLARE theta REAL

# set up initial state
X 0
X 1

# build the exponentiated operator
RX(pi/2) 0
H 1
H 2
H 3

CNOT 0 1
CNOT 1 2
CNOT 2 3
RZ(theta) 3
CNOT 2 3
CNOT 1 2
CNOT 0 1

RX(-pi/2) 0
H 1
H 2
H 3

# measure out the results
MEASURE 0 ro[0]
MEASURE 1 ro[1]
MEASURE 2 ro[2]
MEASURE 3 ro[3]""")
    return program

Next, we modify the execution loop. Rather than reformulating the Program object each time, we build and compile it once, then use the .load() method to transfer the parametric program to the (simulated) quantum device. We then set only the angle value within the inner loop, and we change to using .run() and .wait() methods to manage control between us and the quantum device.

More specifically, the old execution loop

# get all the unweighted expectations for all the sample wavefunctions
occupations = list(range(angle_min, angle_max))
indices = list(range(4))
for offset in occupations:
    # set up the Program object, each time we have a new parameter
    program = build_wf_ansatz_prep(angle_min + offset * angle_step)
    bitstrings = qvm.run(program, indices, trials=shots)

    totals = [0, 0, 0, 0]
    for bitstring in bitstrings:
        for index in indices:
            totals[index] += bitstring[index]
    occupations[offset] = [t / shots for t in totals]

becomes

# set up the Program object, ONLY ONCE
program = build_wf_ansatz_prep()
program.wrap_in_numshots_loop(shots=shots)
nq_program = qc.compiler.quil_to_native_quil(program)
binary = qc.compiler.native_quil_to_executable(nq_program)
qc.qam.load(binary)

# get all the unweighted expectations for all the sample wavefunctions
occupations = list(range(angle_min, angle_max))
indices = list(range(4))
for offset in occupations:
    qc.qam.write_memory(region_name='theta', value=angle_min + offset * angle_step)
    qc.qam.run()
    qc.qam.wait()
    bitstrings = qc.qam.read_memory(region_name="ro")

    totals = [0, 0, 0, 0]
    for bitstring in bitstrings:
        for index in indices:
            totals[index] += bitstring[index]
    occupations[offset] = [t / shots for t in totals]

Overall, the resulting program looks like this:

from pyquil.api import get_qc
from pyquil.quil import Program


def setup_forest_objects():
    qc = get_qc("9q-square-qvm")
    return qc


def build_wf_ansatz_prep():
    program = Program("""
# set up memory
DECLARE ro BIT[4]
DECLARE theta REAL

# set up initial state
X 0
X 1

# build the exponentiated operator
RX(pi/2) 0
H 1
H 2
H 3

CNOT 0 1
CNOT 1 2
CNOT 2 3
RZ(theta) 3
CNOT 2 3
CNOT 1 2
CNOT 0 1

RX(-pi/2) 0
H 1
H 2
H 3

# measure out the results
MEASURE 0 ro[0]
MEASURE 1 ro[1]
MEASURE 2 ro[2]
MEASURE 3 ro[3]""")
    return program


# some constants
bond_step, bond_min, bond_max = 0.05, 0, 200
angle_step, angle_min, angle_max = 0.1, 0, 63
convolution_coefficients = [0.1698845197777728, 0.16988451977777283, -0.2188630663199042,
                            -0.2188630663199042]
shots = 1000

# set up the Forest object
qc = setup_forest_objects()

# set up the Program object, ONLY ONCE
program = build_wf_ansatz_prep()
program.wrap_in_numshots_loop(shots=shots)
nq_program = qc.compiler.quil_to_native_quil(program)
binary = qc.compiler.native_quil_to_executable(nq_program)
qc.qam.load(binary)

# get all the unweighted expectations for all the sample wavefunctions
occupations = list(range(angle_min, angle_max))
indices = list(range(4))
for offset in occupations:
    qc.qam.write_memory(region_name='theta', value=angle_min + offset * angle_step)
    qc.qam.run()
    qc.qam.wait()
    bitstrings = qc.qam.read_memory(region_name="ro")

    totals = [0, 0, 0, 0]
    for bitstring in bitstrings:
        for index in indices:
            totals[index] += bitstring[index]
    occupations[offset] = [t / shots for t in totals]

# compute minimum energy as a function of bond length
min_energies = list(range(bond_min, bond_max))
for bond_length in min_energies:
    energies = []
    for offset in range(angle_min, angle_max):
        energy = 0
        for j in range(4):
            energy += occupations[offset][j] * convolution_coefficients[j]
        energies.append(energy)

    min_energies[bond_length] = min(energies)

min_index = min_energies.index(min(min_energies))
min_energy, relaxed_length = min_energies[min_index], min_index * bond_step

Miscellanea

Quil promises that a BIT is 1 bit and that an OCTET is 8 bits. Quil does not promise, however, the size or layout of INTEGER or REAL. These are implementation-dependent.

On the QPU, INTEGER refers to an unsigned integer stored in a 48-bit wide little-endian word, and REAL refers to a 48-bit wide little-endian fixed point number of type <0.48>. In general, these datatypes are implementation-dependent. OCTET always refers to an 8-bit wide unsigned integer and is independent of implementation.

Memory regions are all “global”: DECLARE directives cannot appear in the body of a DEFCIRCUIT.

On the QVM, INTEGER is a two’s complement signed 64-bit integer. REAL is an IEEE-754 double-precision floating-point number.

Error reporting

Because the Forest 2.0 execution model is no longer asynchronous, our error reporting model has also changed. Rather than writing to technical support with a job ID, users will need to provide all pertinent details to how they produced an error.

PyQuil 2.0 makes this task easy with the function decorator @pyquil_protect, found in the module pyquil.api. By decorating a failing function (or a function that has the potential to fail), any unhandled exceptions will cause an error log to be written to disk (at a user-specifiable location). For example, the nonsense code block

from pyquil.api import pyquil_protect

...

@pyquil_protect
def my_function():
    ...
    qc.qam.load(qc)
    ...

my_function()

causes the following error to be printed:

>>> PYQUIL_PROTECT <<<
An uncaught exception was raised in a function wrapped in pyquil_protect.  We are writing out a
log file to "/Users/your_name/Documents/pyquil/pyquil_error.log".

Along with a description of what you were doing when the error occurred, send this file to Rigetti Computing
support by email at support@rigetti.com for assistance.
>>> PYQUIL_PROTECT <<<

as well as the following log file to be written to disk at the indicated location:

{
  "stack_trace": [
    {
      "name": "pyquil_protect_wrapper",
      "filename": "/Users/your_name/Documents/pyquil/pyquil/error_reporting.py",
      "line_number": 197,
      "locals": {
        "e": "TypeError('quil_binary argument must be a QVMExecutableResponse. This error is typically triggered by
            forgetting to pass (nativized) Quil to native_quil_to_executable or by using a compiler meant to be used
            for jobs bound for a QPU.',)",
        "old_filename": "'pyquil_error.log'",
        "kwargs": "{}",
        "args": "()",
        "log_filename": "'pyquil_error.log'",
        "func": "<function my_function at 0x106dc4510>"
      }
    },
    {
      "name": "my_function",
      "filename": "<stdin>",
      "line_number": 10,
      "locals": {
        "offset": "0",
        "occupations": "[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24,
            25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50,
            51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62]"
      }
    },
    {
      "name": "wrapper",
      "filename": "/Users/your_name/Documents/pyquil/pyquil/error_reporting.py",
      "line_number": 228,
      "locals": {
        "pre_entry": "CallLogValue(timestamp_in=datetime.datetime(2018, 9, 11, 18, 40, 19, 65538),
            timestamp_out=None, return_value=None)",
        "key": "run('<pyquil.api._qvm.QVM object at 0x1027e3940>', )",
        "kwargs": "{}",
        "args": "(<pyquil.api._qvm.QVM object at 0x1027e3940>,)",
        "func": "<function QVM.run at 0x106db4e18>"
      }
    },
    {
      "name": "run",
      "filename": "/Users/your_name/Documents/pyquil/pyquil/api/_qvm.py",
      "line_number": 376,
      "locals": {
        "self": "<pyquil.api._qvm.QVM object at 0x1027e3940>",
        "__class__": "<class 'pyquil.api._qvm.QVM'>"
      }
    }
  ],
  "timestamp": "2018-09-11T18:40:19.253286",
  "call_log": {
    "__init__('<pyquil.api._qvm.QVM object at 0x1027e3940>', '<pyquil.api._base_connection.ForestConnection object at
        0x1027e3588>', )": {
      "timestamp_in": "2018-09-11T18:40:18.967750",
      "timestamp_out": "2018-09-11T18:40:18.968170",
      "return_value": "None"
    },
    "run('<pyquil.api._qvm.QVM object at 0x1027e3940>', )": {
      "timestamp_in": "2018-09-11T18:40:19.065538",
      "timestamp_out": null,
      "return_value": null
    }
  },
  "exception": "TypeError('quil_binary argument must be a QVMExecutableResponse. This error is typically triggered
    by forgetting to pass (nativized) Quil to native_quil_to_executable or by using a compiler meant to be used for
    jobs bound for a QPU.',)",
  "system_info": {
    "python_version": "3.6.3 (default, Jan 25 2018, 13:55:02) \n[GCC 4.2.1 Compatible Apple LLVM 9.0.0
        (clang-900.0.39.2)]",
    "pyquil_version": "2.0.0-internal.1"
  }
}

Please attach such a logfile to any request for support.

Parametric Programs

In PyQuil 1.x, there was an object named ParametricProgram:

# This function returns a quantum circuit with different rotation angles on a gate on qubit 0
def rotator(angle):
    return Program(RX(angle, 0))

from pyquil.parametric import ParametricProgram
par_p = ParametricProgram(rotator) # This produces a new type of parameterized program object

This object has been removed from PyQuil 2. Please consider simply using a Python function for the above functionality:

par_p = rotator

Or using declared classical memory:

p = Program()
angle = p.declare('angle', 'REAL')
p += RX(angle, 0)

Programs and Gates

Note

If you’re running locally, remember set up the QVM and quilc in server mode before trying to use them: Setting Up Server Mode for PyQuil.

Introduction

Quantum programs are written in Forest using the Program object. This Program abstraction will help us compose Quil programs.

from pyquil import Program

Programs are constructed by adding quantum gates to it, which are defined in the gates module. We can import all standard gates with the following:

from pyquil.gates import *

Let’s instantiate a Program and add an operation to it. We will act an X gate on qubit 0.

p = Program()
p += X(0)

All qubits begin in the ground state. This means that if we measure a qubit without applying operations on it, we expect to receive a 0 result. The X gate will rotate qubit 0 from the ground state to the excited state, so a measurement immediately after should return a 1 result. More details about gate operations are explained in Introduction to Quantum Computing.

We can print our pyQuil program (print(p)) to see the equivalent Quil representation:

X 0

This isn’t going to be very useful to us without measurements. In pyQuil 2.0, we have to DECLARE a memory space to read measurement results, which we call “readout results” and abbreviate as ro. With measurement, our whole program looks like this:

from pyquil import Program
from pyquil.gates import *

p = Program()
ro = p.declare('ro', 'BIT', 1)
p += X(0)
p += MEASURE(0, ro[0])

print(p)

DECLARE ro BIT[1]
X 0
MEASURE 0 ro[0]

We’ve instantiated a program, declared a memory space named ro with one single bit of memory, applied an X gate on qubit 0, and finally measured qubit 0 into the zeroth index of the memory space named ro.

Awesome! That’s all we need to get results back. Now we can actually see what happens if we run this program on the Quantum Virtual Machine (QVM). We just have to add a few lines to do this.

from pyquil import get_qc

...

qc = get_qc('1q-qvm')  # You can make any 'nq-qvm' this way for any reasonable 'n'
executable = qc.compile(p)
result = qc.run(executable)
print(result)

Congratulations! You just ran your program on the QVM. The returned value should be:

[[1]]

For more information on what the above result means, and on executing quantum programs on the QVM in general, see The Quantum Computer. The remainder of this section of the docs will be dedicated to constructing programs in detail, an essential part of becoming fluent in quantum programming.

The Standard Gate Set

The following gates methods come standard with Quil and gates.py:

  • Pauli gates I, X, Y, Z
  • Hadamard gate: H
  • Phase gates: PHASE(theta), S, T
  • Controlled phase gates: CZ, CPHASE00(alpha), CPHASE01(alpha), CPHASE10(alpha), CPHASE(alpha)
  • Cartesian rotation gates: RX(theta), RY(theta), RZ(theta)
  • Controlled \(X\) gates: CNOT, CCNOT
  • Swap gates: SWAP, CSWAP, ISWAP, PSWAP(alpha)

The parameterized gates take a real or complex floating point number as an argument.

Declaring Memory

Classical memory regions must be explicitly requested and named by a Quil program using the DECLARE directive. Details about the Quil directive can be found in Details of updates to Quil.

In pyQuil, we declare memory with the .declare method on a Program. Let’s inspect the function signature

# pyquil.quil.Program

def declare(self, name, memory_type='BIT', memory_size=1, shared_region=None, offsets=None):

and break down each argument:

  • name is any name you want to give this memory region.
  • memory_type is one of 'REAL', 'BIT', 'OCTET', or 'INTEGER' (given as a string). Only BIT and OCTET always have a determined size, which is 1 bit and 8 bits respectively.
  • memory_size is the number of elements of that type to reserve.
  • shared_region and offsets allow you to alias memory regions. For example, you might want to name the third bit in your readout array as q3_ro. SHARING is currently disallowed for our QPUs, so we won’t focus on this here.

Now we can get into an example.

from pyquil import Program

p = Program()
ro = p.declare('ro', 'BIT', 16)
theta = p.declare('theta', 'REAL')

Warning

.declare cannot be chained, since it doesn’t return a modified Program object.

Notice that the .declare method returns a reference to the memory we’ve just declared. We will need this reference to make use of these memory spaces again. Let’s see how the Quil is looking so far:

DECLARE ro BIT[16]
DECLARE theta REAL[1]

That’s all we have to do to declare the memory. Continue to the next section on Measurement to learn more about using ro to store measured readout results. Check out Parametric Compilation to see how you might use theta to compile gate parameters dynamically.

Measurement

There are several ways you can handle measurements in your program. We will start with the simplest method – letting the QuantumComputer abstraction do it for us.

from pyquil import Program, get_qc
from pyquil.gates import H, CNOT

# Get our QuantumComputer instance, with a Quantum Virutal Machine (QVM) backend
qc = get_qc("8q-qvm")

# Construct a simple Bell State
p = Program(H(0), CNOT(0, 1))

results = qc.run_and_measure(p, trials=10)
print(results)

{0: array([1, 1, 0, 1, 0, 0, 1, 1, 0, 1]),
 1: array([1, 1, 0, 1, 0, 0, 1, 1, 0, 1]),
 2: array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0]),
 3: array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0]),
 4: array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0]),
 5: array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0]),
 6: array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0]),
 7: array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0])}

The method .run_and_measure will handle declaring memory for readout results, adding MEASURE instructions for each qubit in the QVM, telling the QVM how many trials to run, running and returning the measurement results.

You might sometimes want finer grained control. In this case, we’re probably only interested in the results on qubits 0 and 1, but .run_and_measure returns the results for all eight qubits in the QVM. We can change our program to be more particular about what we want.

from pyquil import Program
from pyquil.gates import *

p = Program()
ro = p.declare('ro', 'BIT', 2)
p += H(0)
p += CNOT(0, 1)
p += MEASURE(0, ro[0])
p += MEASURE(1, ro[1])

In the last two lines, we’ve added our MEASURE instructions, saying that we want to store the result of qubit 0 into the 0th bit of ro, and the result of qubit 1 into the 1st bit of ro. The following snippet could be a useful way to measure many qubits, in particular, on a lattice that doesn’t start at qubit 0 (although you can use the compiler to re-index your qubits):

qubits = [5, 6, 7]
# ...
for i, q in enumerate(qubits):
    p += MEASURE(q, ro[i])

Note

The QPU can only handle MEASURE final programs. You can’t operate gates after measurements.

Specifying the number of trials

Quantum computing is inherently probabilistic. We often have to repeat the same experiment many times to get the results we need. Sometimes we expect the results to all be the same, such as when we apply no gates, or only an X gate. When we prepare a superposition state, we expect probabilistic outcomes, such as a 50% probability measuring 0 or 1.

The number of shots (also called trials) is the number of times to execute a program at once. This determines the length of the results that are returned.

As we saw above, the .run_and_measure method of the QuantumComputer object can handle multiple executions of a program. If you would like more explicit control for representing multi-shot execution, another way to do this is with .wrap_in_numshots_loop. This puts the number of shots to be run in the representation of the program itself, as opposed to in the arguments list of the execution method itself. Below, we specify that our program should be executed 1000 times.

p = Program()
...   # build up your program here...
p.wrap_in_numshots_loop(1000)

Note

Did You Know?

The word “shot” comes from experimental physics where an experiment is performed many times, and each result is called a shot.

Parametric Compilation

Modern quantum algorithms are often parametric, following a hybrid model. In this hybrid model, the program ansatz (template of gates) is fixed, and iteratively updated with new parameters. These new parameters are often determined by an update given by a classical optimizer. Depending on the complexity of the algorithm, problem of interest, and capabilities of the classical optimizer, this loop may need to run many times. In order to efficiently operate within this hybrid model, parametric compilation can be used.

Parametric compilation allows one to compile the program ansatz just once. Making use of declared memory regions, we can load values to the parametric gates at execution time, after compilation. Taking the compiler out of the execution loop for programs like this offers a huge performance improvement compared to compiling the program each time a parameter update is required. (Some more details about this and an example are found :ref:here <parametric>.)

The first step is to build our parametric program, which functions like a template for all the precise programs we will run. Below we create a simple example program to illustrate, which puts the qubit onto the equator of the Bloch Sphere and then rotates it around the Z axis for some variable angle theta before applying another X pulse and measuring.

import numpy as np

from pyquil import Program
from pyquil.gates import RX, RZ, MEASURE

qubit = 0

p = Program()
ro = p.declare("ro", "BIT", 1)
theta_ref = p.declare("theta", "REAL")

p += RX(np.pi / 2, qubit)
p += RZ(theta, qubit)
p += RX(-np.pi / 2, qubit)

p += MEASURE(qubit, ro[0])

Note

The example program, although simple, is actually more than just a toy example. It is similar to an experiment which measures the qubit frequency.

Notice how theta hasn’t been specified yet. The next steps will have to involve a QuantumComputer or a compiler implementation. For simplicity, we will demonstrate with a QuantumComputer instance.

from pyquil import get_qc

# Get a Quantum Virtual Machine to simulate execution
qc = get_qc("1q-qvm")
executable = qc.compile(p)

We are able to compile our program, even with theta still not specified. Now we want to run our program with theta filled in for, say, 200 values between \(0\) and \(2\pi\). We demonstrate this below.

# Load the executable onto the QVM.
# QAM is a quantum abstract machine, which can be a QVM or a QPU
qc.qam.load(executable)

# Somewhere to store each list of results
parametric_measurements = []

for theta in np.linspace(0, 2 * np.pi, 200):
    # Write the value we want to execute with
    qc.qam.write_memory(region_name='theta', value=theta)
    # Get the results of the run
    bitstrings = qc.qam.run().wait().read_memory(region_name='ro')
    # Store our results
    parameteric_measurements.append(bitstrings)

Gate Modifiers

Gate applications in Quil can be preceded by a gate modifier. There are two supported modifiers: DAGGER and CONTROLLED. The DAGGER modifier represents the dagger of the gate. For instance,

DAGGER RX(pi/3) 0

would have an equivalent effect to RX(-pi/3) 0.

The CONTROLLED modifier takes a gate and makes it a controlled gate. For instance, one could write the Toffoli gate in any of the three following ways:

CCNOT 0 1 2
CONTROLLED CNOT 0 1 2
CONTROLLED CONTROLLED X 0 1 2

Note

The letter C in the gate name has no semantic significance in Quil. To make a controlled Y gate, one cannot write CY, but rather one has to write CONTROLLED Y.

Defining New Gates

New gates can be easily added inline to Quil programs. All you need is a matrix representation of the gate. For example, below we define a \(\sqrt{X}\) gate.

import numpy as np

from pyquil import Program
from pyquil.quil import DefGate

# First we define the new gate from a matrix
sqrt_x = np.array([[ 0.5+0.5j,  0.5-0.5j],
                   [ 0.5-0.5j,  0.5+0.5j]])

# Get the Quil definition for the new gate
sqrt_x_definition = DefGate("SQRT-X", sqrt_x)
# Get the gate constructor
SQRT_X = sqrt_x_definition.get_constructor()

# Then we can use the new gate
p = Program()
p += sqrt_x_definition
p += SQRT_X(0)
print(p)

DEFGATE SQRT-X:
    0.5+0.5i, 0.5-0.5i
    0.5-0.5i, 0.5+0.5i

SQRT-X 0

Below we show how we can define \(X_0\otimes \sqrt{X_1}\) as a single gate.

# A multi-qubit defgate example
x_gate_matrix = np.array(([0.0, 1.0], [1.0, 0.0]))
sqrt_x = np.array([[ 0.5+0.5j,  0.5-0.5j],
                [ 0.5-0.5j,  0.5+0.5j]])
x_sqrt_x = np.kron(x_gate_matrix, sqrt_x)

Now we can use this gate in the same way that we used SQRT_X, but we will pass it two arguments rather than one, since it operates on two qubits.

x_sqrt_x_definition = DefGate("X-SQRT-X", x_sqrt_x)
X_SQRT_X = x_sqrt_x_definition.get_constructor()

# Then we can use the new gate
p = Program(x_sqrt_x_definition, X_SQRT_X(0, 1))

Tip

To inspect the wavefunction that will result from applying your new gate, you can use the Wavefunction Simulator (e.g. print(WavefunctionSimulator().wavefunction(p))).

Defining Parametric Gates

Let’s say we want to have a controlled RX gate. Since RX is a parametric gate, we need a slightly different way of defining it than in the previous section.

from pyquil import Program, WavefunctionSimulator
from pyquil.parameters import Parameter, quil_sin, quil_cos
from pyquil.quilbase import DefGate
import numpy as np

# Define the new gate from a matrix
theta = Parameter('theta')
crx = np.array([
    [1, 0, 0, 0],
    [0, 1, 0, 0],
    [0, 0, quil_cos(theta / 2), -1j * quil_sin(theta / 2)],
    [0, 0, -1j * quil_sin(theta / 2), quil_cos(theta / 2)]
])

gate_definition = DefGate('CRX', crx, [theta])
CRX = gate_definition.get_constructor()

# Create our program and use the new parametric gate
p = Program()
p += gate_definition
p += H(0)
p += CRX(np.pi/2)(0, 1)

quil_sin and quil_cos work as the regular sines and cosines, but they support the parametrization. Parametrized functions you can use with pyQuil are: quil_sin, quil_cos, quil_sqrt, quil_exp, and quil_cis.

Tip

To inspect the wavefunction that will result from applying your new gate, you can use the Wavefunction Simulator (e.g. print(WavefunctionSimulator().wavefunction(p))).

Pragmas

PRAGMA directives give users more control over how Quil programs are processed or simulated but generally do not change the semantics of the Quil program itself. As a general rule of thumb, deleting all PRAGMA directives in a Quil program should leave a valid and semantically equivalent program.

In pyQuil, PRAGMA directives play many roles, such as controlling the behavior of gates in noisy simulations, or commanding the Quil compiler to perform actions in a certain way. Here, we will cover the basics of two very common use cases for including a PRAGMA in your program: qubit rewiring and delays. For a more comprehensive review of what pragmas are and what the compiler supports, check out The Quil Compiler. For more information about PRAGMA in Quil, see A Practical Quantum ISA, and Simulating Quantum Processor Errors.

Specifying A Qubit Rewiring Scheme

Qubit rewiring is one of the most powerful features of the Quil compiler. We are able to write Quil programs which are agnostic to the topology of the chip, and the compiler will intelligently relabel our qubits to give better performance.

When we intend to run a program on the QPU, sometimes we write programs which use specific qubits targeting a specific device topology, perhaps to achieve a high-performance program. Other times, we write programs that are agnostic to the underlying topology, thereby making the programs more portable. Qubit rewiring accommodates both use cases in an automatic way.

Consider the following program.

from pyquil import Program
from pyquil.gates import *

p = Program(X(3))

We’ve tested this on the QVM, and we’ve reserved a lattice on the QPU which has qubits 4, 5, and 6, but not qubit 3. Rather than rewrite our program for each reservation, we modify our program to tell the compiler to do this for us.

from pyquil.quil import Pragma

p = Program(Pragma('INITIAL_REWIRING', ['"GREEDY"']))
p += X(3)

Now, when we pass our program through the compiler (such as with QuantumComputer.compile()) we will get native Quil with the qubit reindexed to one of 4, 5, or 6. If qubit 3 is available, and we don’t want that pulse to be applied to any other qubit, we would instead use Pragma('INITIAL_REWIRING', ['"NAIVE"']]. Detailed information about the available options is here.

Note

In general, we assume that the qubits you’re supplying as input are also the ones which you prefer to operate on, and so NAIVE rewiring is the default.

Asking for a Delay

At times, we may want to add a delay in our program. Usually this is associated with qubit characterization. Delays are not regular gate operations, and they do not affect the abstract semantics of the Quil program, so they’re implemented with a PRAGMA directive.

#  ...
# qubit index and time in seconds must be defined and provided
p += Pragma('DELAY', [qubit], str(time))

Warning

Keep in mind, the program duration is currently capped at 15 seconds, and the length of the program is multiplied by the number of shots. If you have a 1000 shot program, where each shot contains a 100ms delay, you won’t be able to execute it.

Ways to Construct Programs

PyQuil supports a variety of methods for constructing programs however you prefer. Multiple instructions can be applied at once, and programs can be added together. PyQuil can also produce a Program by interpreting raw Quil text. You can still use the more pyQuil 1.X style of using the .inst method to add instruction gates. Thus, the following are all valid programs:

# Preferred method
p = Program()
p += X(0)
p += Y(1)
print(p)

# Multiple instructions in declaration
print(Program(X(0), Y(1)))

# A composition of two programs
print(Program(X(0)) + Program(Y(1)))

# Raw Quil with newlines
print(Program("X 0\nY 1"))

# Raw Quil comma separated
print(Program("X 0", "Y 1"))

# Chained inst; less preferred
print(Program().inst(X(0)).inst(Y(1)))

All of the above methods will produce the same output:

X 0
Y 1

The pyquil.parser submodule provides a front-end to other similar parser functionality.

Fixing a Mistaken Instruction

If an instruction was appended to a program incorrectly, you can pop it off.

p = Program(X(0), Y(1))
print(p)

print("We can fix by popping:")
p.pop()
print(p)

X 0
Y 1

We can fix by popping:
X 0

QPU-allowable Quil

Apart from DECLARE and PRAGMA directives, a program must break into the following three regions, each optional:

  1. A RESET command.
  2. A sequence of quantum gate applications.
  3. A sequence of MEASURE commands.

The only memory that is writeable is the region named ro, and only through MEASURE instructions. All other memory is read-only.

The keyword SHARING is disallowed.

Compilation is unavailable for invocations of DEFGATEs with parameters read from classical memory.

The Quantum Computer

The Quantum Virtual Machine (QVM)

The Rigetti Quantum Virtual Machine is an implementation of the Quantum Abstract Machine from A Practical Quantum Instruction Set Architecture. [1] It is implemented in ANSI Common LISP and executes programs specified in the Quantum Instruction Language (Quil).

Quil is an opinionated quantum instruction language: its basic belief is that in the near term quantum computers will operate as coprocessors, working in concert with traditional CPUs. This means that Quil is designed to execute on a Quantum Abstract Machine that has a shared classical/quantum architecture at its core.

The QVM is a wavefunction simulation of unitary evolution with classical control flow and shared quantum classical memory.

Using the QVM

After downloading the SDK, the QVM is available on your local machine. You can initialize a local QVM server by typing qvm -S into your terminal. You should see the following message.

$ qvm -S
******************************
* Welcome to the Rigetti QVM *
******************************
Copyright (c) 2018 Rigetti Computing.

(Configured with 2048 MiB of workspace and 8 workers.)

[2018-11-06 18:18:18] Starting server on port 5000.

For a detailed description of how to use the qvm from the command line, see The QVM manual page.

Once the QVM is serving requests, we can run the following pyQuil program to get a QuantumComputer object which will use the QVM.

from pyquil import get_qc, Program
from pyquil.gates import *
qc = get_qc('9q-square-qvm')

One executes quantum programs on the QVM using the .run(...) method, intended to closely mirror how one will execute programs on a real QPU. We also offer a Wavefunction Simulator (formerly a part of the QVM object), which allows users to contruct and inspect wavefunctions of quantum programs. Learn more about the Wavefunction Simulator here. For information on constructing quantum programs, please refer back to Programs and Gates.

The .run(...) method

The .run(...) method takes in a compiled program. You are responsible for compiling your program before running it. Remember to also start up a quilc compiler server, too, with quilc -S.

p = Program()
ro = p.declare('ro', 'BIT', 1)
p += X(0)
p += MEASURE(0, ro[0])
p += MEASURE(1, ro[1])
p.wrap_in_numshots_loop(5)
executable = qc.compile(p)
results = qvm.run(executable)
print(results)

The results returned are a list of lists of integers. In the above case, that’s

[[1, 0], [1, 0], [1, 0], [1, 0], [1, 0]]

Let’s unpack this. The outer list is an enumeration over the trials; the argument given to wrap_in_numshots_loop will match the length of results.

The inner list, on the other hand, is an enumeration over the results stored in the memory region named ro, which we use as our readout register. We see that the result of this program is that the memory region ro[0] now stores the state of qubit 0, which should be 1 after an \(X\)-gate. See Declaring Memory and Measurement for more details about declaring and accessing classical memory regions.

[1]https://arxiv.org/abs/1608.03355
Simulating the QPU using the QVM

The QVM is a powerful tool for testing quantum programs before executing them on the QPU.

qc = get_qc("QuantumComputerName")
qc = get_qc("QuantumComputerName-qvm")

By simply providing -qvm in the device name, all programs executed on this QVM will, have the same topology as the named QPU. To learn how to add noise models to your virtual QuantumComputer instance, check out Noise and Quantum Computation.

The Quantum Processing Unit

Coming soon: Detailed information about how to use get_qc() to target a QPU.

The Wavefunction Simulator

Formerly a part of the QVM object in pyQuil, the Wavefunction Simulator allows you to directly inspect the wavefunction of a quantum state prepared by your program. Because of the probabilistic nature of quantum information, the programs you’ll be running on the QPU can give a distribution of outputs. When running on the QPU or QVM, you would aggregate results (anywhere from tens of trials to 100k+!) that you can sample to get back a distribution.

With the Wavefunction Simulator, you can look at the distribution without having to collect samples from your program. This can save a lot of time for small programs. Let’s walk through a basic example of using WavefunctionSimulator:

from pyquil import Program
from pyquil.gates import *
from pyquil.api import WavefunctionSimulator
wf_sim = WavefunctionSimulator()
coin_flip = Program((H(0))
wf_sim.wavefunction(coin_flip)

<pyquil.wavefunction.Wavefunction at 0x1088a2c10>

The return value is a Wavefunction object that stores the amplitudes of the quantum state. We can print this object

coin_flip = Program(H(0))
wavefunction = wf_sim.wavefunction(coin_flip)
print(wavefunction)

(0.7071067812+0j)|0> + (0.7071067812+0j)|1>

to see the amplitudes listed as a sum of computational basis states. We can index into those amplitudes directly or look at a dictionary of associated outcome probabilities.

assert wavefunction[0] == 1 / np.sqrt(2)
# The amplitudes are stored as a numpy array on the Wavefunction object
print(wavefunction.amplitudes)
prob_dict = wavefunction.get_outcome_probs() # extracts the probabilities of outcomes as a dict
print(prob_dict)
prob_dict.keys() # these store the bitstring outcomes
assert len(wavefunction) == 1 # gives the number of qubits

[ 0.70710678+0.j  0.70710678+0.j]
{'1': 0.49999999999999989, '0': 0.49999999999999989}

It is important to remember that this wavefunction method is a useful debugging tool for small quantum systems, and cannot be feasibly obtained on a quantum processor.

Multi-Qubit Basis Enumeration

The WavefunctionSimulator enumerates bitstrings such that qubit 0 is the least significant bit (LSB) and therefore on the right end of a bitstring as shown in the table below which contains some examples.

bitstring qubit_(n-1) … qubit_2 qubit_1 qubit_0
1…101 1 … 1 0 1
0…110 0 … 1 1 0

This convention is counter to that often found in the quantum computing literature where bitstrings are often ordered such that the lowest-index qubit is on the left. The vector representation of a wavefunction assumes the “canonical” ordering of basis elements. I.e., for two qubits this order is 00, 01, 10, 11. In the typical Dirac notation for quantum states, the tensor product of two different degrees of freedom is not always explicitly understood as having a fixed order of those degrees of freedom. This is in contrast to the kronecker product between matrices which uses the same mathematical symbol and is clearly not commutative. This, however, becomes important when writing things down as coefficient vectors or matrices:

\[\begin{split}\ket{0}_0 \otimes \ket{1}_1 = \ket{1}_1 \otimes \ket{0}_0 = \ket{10}_{1,0} \equiv \begin{pmatrix} 0 \\ 0 \\ 1 \\ 0 \end{pmatrix}\end{split}\]

As a consequence there arise some subtle but important differences in the ordering of wavefunction and multi-qubit gate matrix coefficients. According to our conventions the matrix

\[\begin{split}U_{\rm CNOT(1,0)} \equiv \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 \end{pmatrix}\end{split}\]

corresponds to the Quil instruction CNOT(1, 0) which is counter to how most other people in the field order their tensor product factors (or more specifically their kronecker products). In this convention CNOT(0, 1) is given by

\[\begin{split}U_{\rm CNOT(0,1)} \equiv \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 \\ 0 & 1 & 0 & 0 \end{pmatrix}\end{split}\]

For additional information why we decided on this basis ordering check out our note Someone shouts, “|01000>!” Who is Excited?.

The Quil Compiler

Expectations for Program Contents

The QPUs have much more limited natural gate sets than the standard gate set offered by pyQuil: on Rigetti QPUs, the gate operators are constrained to lie in RZ(θ), RX(k*π/2), and CZ; and the gates are required to act on physically available hardware (for single-qubit gates, this means acting only on live qubits, and for qubit-pair gates, this means acting on neighboring qubits). However, as a programmer, it is often (though not always) desirable to to be able to write programs which don’t take these details into account. These generally leads to more portable code if one isn’t tied to a specific set of gates or QPU architecture. To ameliorate these limitations, the Rigetti software toolkit contains an optimizing compiler that translates arbitrary Quil to native Quil and native ProtoQuil to executables suitable for Rigetti hardware.

Interacting with the Compiler

After downloading the SDK, the Quil Compiler, quilc is available on your local machine. You can initialize a local quilc server by typing quilc -S into your terminal. You should see the following message.

$ quilc -S
+-----------------+
|  W E L C O M E  |
|   T O   T H E   |
|  R I G E T T I  |
|     Q U I L     |
| C O M P I L E R |
+-----------------+
Copyright (c) 2018 Rigetti Computing.

This is a part of the Forest SDK. By using this program
you agree to the End User License Agreement (EULA) supplied
with this program. If you did not receive the EULA, please
contact <support@rigetti.com>.

[2018-11-06 10:59:22] Starting server: 0.0.0.0 : 6000.

To get a description of quilc, and options and examples of its command line use, see QUILC Man Page.

A QuantumComputer object supplied by the function pyquil.api.get_qc() comes equipped with a connection to your local Rigetti Quil compiler. This can be accessed using the instance method .compile(), as in the following:

from pyquil.quil import Pragma, Program
from pyquil.api import get_qc
from pyquil.gates import CNOT, H

qc = get_qc("9q-square-qvm")

ep = qc.compile(Program(H(0), CNOT(0,1), CNOT(1,2)))

print(ep.program) # here ep is of type PyquilExecutableResponse, which is not always inspectable

with output

PRAGMA EXPECTED_REWIRING "#(0 1 2 3 4 5 6 7)"
RZ(pi/2) 0
RX(pi/2) 0
RZ(-pi/2) 1
RX(pi/2) 1
CZ 1 0
RX(-pi/2) 1
RZ(-pi/2) 2
RX(pi/2) 2
CZ 2 1
RZ(-pi/2) 0
RZ(-pi/2) 1
RX(-pi/2) 2
RZ(pi/2) 2
PRAGMA CURRENT_REWIRING "#(0 1 2 3 4 5 6 7)"

The compiler connection is also available directly via the property qc.compiler. The precise class of this object changes based on context (e.g., QPUCompiler, QVMCompiler, or LocalQVMCompiler), but it always conforms to the interface laid out by pyquil.api._qac:

  • compiler.quil_to_native_quil(program): This method converts a Quil program into native Quil, according to the ISA that the compiler is initialized with. The input parameter is specified as a Program object, and the output is given as a new Program object, equipped with a .metadata property that gives extraneous information about the compilation output (e.g., gate depth, as well as many others). This call blocks until Quil compilation finishes.
  • compiler.native_quil_to_executable(nq_program): This method converts a ProtoQuil program, which is promised to consist only of native gates for a given ISA, into an executable suitable for submission to one of a QVM or a QPU. This call blocks until the executable is generated.

The instance method qc.compile described above is a combination of these two methods: first the incoming Quil is nativized, and then that is immediately turned into an executable. Accordingly, the previous example snippet is identical to the following:

from pyquil.quil import Pragma, Program
from pyquil.api import get_qc
from pyquil.gates import CNOT, H

qc = get_qc("9q-square-qvm")

p = Program(H(0), CNOT(0,1), CNOT(1,2))
np = qc.compiler.quil_to_native_quil(p)
ep = qc.compiler.native_quil_to_executable(np)

print(ep.program) # here ep is of type PyquilExecutableResponse, which is not always inspectable

Region-specific compiler features through PRAGMA

The Quil compiler can also be communicated with through PRAGMA commands embedded in the Quil program.

Note

The interface to the Quil compiler from pyQuil is under construction, and some of the PRAGMA directives will soon be replaced by finer-grained method calls.

Preserved regions

The compiler can be circumvented in user-specified regions. The start of such a region is denoted by PRAGMA PRESERVE_BLOCK, and the end is denoted by PRAGMA END_PRESERVE_BLOCK. The Quil compiler promises not to modify any instructions contained in such a region.

Warning

If a preserved block is not legal QPU input, then it is not guaranteed to execute or it may produced unexpected results.

The following is an example of a program that prepares a Bell state on qubits 0 and 1, then performs a time delay to invite noisy system interaction before measuring the qubits. The time delay region is marked by PRAGMA PRESERVE_BLOCK and PRAGMA END_PRESERVE_BLOCK; without these delimiters, the compiler will remove the identity gates that serve to provide the time delay. However, the regions outside of the PRAGMA region will still be compiled, converting the Bell state preparation to the native gate set.

DECLARE ro BIT[2]

#   prepare a Bell state
H 0
CNOT 0 1

#   wait a while
PRAGMA PRESERVE_BLOCK
I 0
I 1
I 0
I 1
# ...
I 0
I 1
PRAGMA END_PRESERVE_BLOCK

#   and read out the results
MEASURE 0 ro[0]
MEASURE 1 ro[1]
Parallelizable regions

The compiler can sometimes arrange gate sequences more cleverly if the user gives it hints about sequences of gates that commute. A region containing commuting sequences is bookended by PRAGMA COMMUTING_BLOCKS and PRAGMA END_COMMUTING_BLOCKS; within such a region, a given commuting sequence is bookended by PRAGMA BLOCK and PRAGMA END_BLOCK.

Warning

Lying to the compiler about what blocks can commute can cause incorrect results.

The following snippet demonstrates this hinting syntax in a context typical of VQE-type algorithms: after a first stage of performing some state preparation on individual qubits, there is a second stage of “mixing operations” that both re-use qubit resources and mutually commute, followed by a final rotation and measurement. The following program is naturally laid out on a ring with vertices (read either clockwise or counterclockwise) as 0, 1, 2, 3. After scheduling the first round of preparation gates, the compiler will use the hinting to schedule the first and third blocks (which utilize qubit pairs 0-1 and 2-3) before the second and fourth blocks (which utilize qubit pairs 1-2 and 0-3), resulting in a reduction in circuit depth by one half. Without hinting, the compiler will instead execute the blocks in their written order.

DECLARE ro BIT[4]

# Stage one
H 0
H 1
H 2
H 3

# Stage two
PRAGMA COMMUTING_BLOCKS
PRAGMA BLOCK
CNOT 0 1
RZ(0.4) 1
CNOT 0 1
PRAGMA END_BLOCK
PRAGMA BLOCK
CNOT 1 2
RZ(0.6) 2
CNOT 1 2
PRAGMA END_BLOCK
PRAGMA BLOCK
CNOT 2 3
RZ(0.8) 3
CNOT 2 3
PRAGMA END_BLOCK
PRAGMA BLOCK
CNOT 0 3
RZ(0.9) 3
CNOT 0 3
PRAGMA END_BLOCK
PRAGMA END_COMMUTING_BLOCKS

# Stage three
H 0
H 1
H 2
H 3

MEASURE 0 ro[0]
MEASURE 1 ro[1]
MEASURE 2 ro[2]
MEASURE 3 ro[3]
Rewirings

When a Quil program contains multi-qubit instructions that do not name qubit-qubit links present on a target device, the compiler will rearrange the qubits so that execution becomes possible. In order to help the user understand what rearrangement may have been done, the compiler emits two forms of PRAGMA: PRAGMA EXPECTED_REWIRING and PRAGMA CURRENT_REWIRING. From the perspective of the user, both PRAGMA instructions serve the same purpose: PRAGMA ..._REWIRING "#(n0 n1 ... nk)" indicates that the logical qubit labeled j in the program has been assigned to lie on the physical qubit labeled nj on the device. This is strictly for human-readability: user-supplied instructions of the form PRAGMA [EXPECTED|CURRENT]_REWIRING are discarded and have no effect.

In addition, you have some control over how the compiler constructs its rewiring, which is controlled by PRAGMA INITIAL_REWIRING. The syntax is as follows.

# <type> can be NAIVE, RANDOM, PARTIAL, or GREEDY
#
# The double quotes are required.
PRAGMA INITIAL_REWIRING "<type>"

Including this before any non-pragmas will allow the compiler to alter its rewiring behavior. The possible options are:

  • NAIVE (default): The compiler will start with an identity mapping as the initial rewiring. In particular, qubits will not be rewired unless the program requests a qubit-qubit interaction not natively available on the QPU.
  • PARTIAL: The compiler will start with nothing assigned to each physical qubit. Then, it will fill in the logical-to-physical mapping as it encounters new qubits in the program, making its best guess for where they should be placed.
  • RANDOM: the compiler will start with a random permutation.
  • GREEDY: the compiler will make a guess for the initial rewiring based on a quick initial scan of the entire program.

Note

NAIVE rewiring is the default, and for the most part, it follows the “Do What I Mean” (DWIM) principle. It is the least sophisticated, but attempts to follow what the user has constructed with their program. Choosing another rewiring, such as PARTIAL, may lead to higher-performing programs because the compiler has more freedom to optimize the layout of the gates on the qubits.

Common Error Messages

The compiler itself is subject to some limitations, and some of the more commonly observed errors follow:

  • ! ! ! Error: Matrices do not lie in the same projective class. The compiler attempted to decompose an operator as native Quil instructions, and the resulting instructions do not match the original operator. This can happen when the original operator is not a unitary matrix, and could indicate an invalid DEFGATE block. In some rare circumstances, it can also happen due to floating point precision issues.

Noise and Quantum Computation

Modeling Noisy Quantum Gates

Pure States vs. Mixed States

Errors in quantum computing can introduce classical uncertainty in what the underlying state is. When this happens we sometimes need to consider not only wavefunctions but also probabilistic sums of wavefunctions when we are uncertain as to which one we have. For example, if we think that an X gate was accidentally applied to a qubit with a 50-50 chance then we would say that there is a 50% chance we have the \(\ket{0}\) state and a 50% chance that we have a \(\ket{1}\) state. This is called an “impure” or “mixed”state in that it isn’t just a wavefunction (which is pure) but instead a distribution over wavefunctions. We describe this with something called a density matrix, which is generally an operator. Pure states have very simple density matrices that we can write as an outer product of a ket vector \(\ket{\psi}\) with its own bra version \(\bra{\psi}=\ket{\psi}^\dagger\). For a pure state the density matrix is simply

\[\rho_\psi = \ket{\psi}\bra{\psi}.\]

The expectation value of an operator for a mixed state is given by

\[\langle X \rangle_\rho = \tr{X \rho}\]

where \(\tr{\cdot}\) is the trace of an operator, which is the sum of its diagonal elements, which is independent of choice of basis. Pure state density matrices satisfy

\[\rho \text{ is pure } \Leftrightarrow \rho^2 = \rho\]

which you can easily verify for \(\rho_\psi\) assuming that the state is normalized. If we want to describe a situation with classical uncertainty between states \(\rho_1\) and \(\rho_2\), then we can take their weighted sum

\[\rho = p \rho_1 + (1-p) \rho_2\]

where \(p\in [0,1]\) gives the classical probability that the state is \(\rho_1\).

Note that classical uncertainty in the wavefunction is markedly different from superpositions. We can represent superpositions using wavefunctions, but use density matrices to describe distributions over wavefunctions. You can read more about density matrices here [DensityMatrix].

[DensityMatrix]https://en.wikipedia.org/wiki/Density_matrix
Quantum Gate Errors

For a quantum gate given by its unitary operator \(U\), a “quantum gate error” describes the scenario in which the actually induces transformation deviates from \(\ket{\psi} \mapsto U\ket{\psi}\). There are two basic types of quantum gate errors:

  1. coherent errors are those that preserve the purity of the input state, i.e., instead of the above mapping we carry out a perturbed, but unitary operation \(\ket{\psi} \mapsto \tilde{U}\ket{\psi}\), where \(\tilde{U} \neq U\).

  2. incoherent errors are those that do not preserve the purity of the input state, in this case we must actually represent the evolution in terms of density matrices. The state \(\rho := \ket{\psi}\bra{\psi}\) is then mapped as

    \[\rho \mapsto \sum_{j=1}^n K_j\rho K_j^\dagger,\]

    where the operators \(\{K_1, K_2, \dots, K_m\}\) are called Kraus operators and must obey \(\sum_{j=1}^m K_j^\dagger K_j = I\) to conserve the trace of \(\rho\). Maps expressed in the above form are called Kraus maps. It can be shown that every physical map on a finite dimensional quantum system can be represented as a Kraus map, though this representation is not generally unique. You can find more information about quantum operations here

In a way, coherent errors are in principle amendable by more precisely calibrated control. Incoherent errors are more tricky.

Why Do Incoherent Errors Happen?

When a quantum system (e.g., the qubits on a quantum processor) is not perfectly isolated from its environment it generally co-evolves with the degrees of freedom it couples to. The implication is that while the total time evolution of system and environment can be assumed to be unitary, restriction to the system state generally is not.

Let’s throw some math at this for clarity: Let our total Hilbert space be given by the tensor product of system and environment Hilbert spaces: \(\mathcal{H} = \mathcal{H}_S \otimes \mathcal{H}_E\). Our system “not being perfectly isolated” must be translated to the statement that the global Hamiltonian contains a contribution that couples the system and environment:

\[H = H_S \otimes I + I \otimes H_E + V\]

where \(V\) non-trivally acts on both the system and the environment. Consequently, even if we started in an initial state that factorized over system and environment \(\ket{\psi}_{S,0}\otimes \ket{\psi}_{E,0}\) if everything evolves by the Schrödinger equation

\[\ket{\psi_t} = e^{-i \frac{Ht}{\hbar}} \left(\ket{\psi}_{S,0}\otimes \ket{\psi}_{E,0}\right)\]

the final state will generally not admit such a factorization.

A Toy Model

In this (somewhat technical) section we show how environment interaction can corrupt an identity gate and derive its Kraus map. For simplicity, let us assume that we are in a reference frame in which both the system and environment Hamiltonian’s vanish \(H_S = 0, H_E = 0\) and where the cross-coupling is small even when multiplied by the duration of the time evolution \(\|\frac{tV}{\hbar}\|^2 \sim \epsilon \ll 1\) (any operator norm \(\|\cdot\|\) will do here). Let us further assume that \(V = \sqrt{\epsilon} V_S \otimes V_E\) (the more general case is given by a sum of such terms) and that the initial environment state satisfies \(\bra{\psi}_{E,0} V_E\ket{\psi}_{E,0} = 0\). This turns out to be a very reasonable assumption in practice but a more thorough discussion exceeds our scope.

Then the joint system + environment state \(\rho = \rho_{S,0} \otimes \rho_{E,0}\) (now written as a density matrix) evolves as

\[\rho \mapsto \rho' := e^{-i \frac{Vt}{\hbar}} \rho e^{+i \frac{Vt}{\hbar}}\]

Using the Baker-Campbell-Hausdorff theorem we can expand this to second order in \(\epsilon\)

\[\rho' = \rho - \frac{it}{\hbar} [V, \rho] - \frac{t^2}{2\hbar^2} [V, [V, \rho]] + O(\epsilon^{3/2})\]

We can insert the initially factorizable state \(\rho = \rho_{S,0} \otimes \rho_{E,0}\) and trace over the environmental degrees of freedom to obtain

\[\begin{split}\begin{align} \rho_S' := \tr{\rho'}_E & = \rho_{S,0} \underbrace{\tr{\rho_{E,0}}_{1}} - \frac{i\sqrt{\epsilon} t}{\hbar} \underbrace{\left[ V_S \rho_{S,0} \underbrace{\tr{V_E\rho_{E,0}}}_{\bra{\psi}_{E,0} V_E\ket{\psi}_{E,0} = 0} - \rho_{S,0}V_S \underbrace{\tr{\rho_{E,0}V_E}}_{\bra{\psi}_{E,0} V_E\ket{\psi}_{E,0} = 0} \right]}_0 \\ & - \frac{\epsilon t^2}{2\hbar^2} \left[ V_S^2\rho_{S,0}\tr{V_E^2 \rho_{E,0}} + \rho_{S,0} V_S^2 \tr{\rho_{E,0}V_E^2} - 2 V_S\rho_{S,0}V_S\tr{V_E \rho_{E,0}V_E}\right] \\ & = \rho_{S,0} - \frac{\gamma}{2} \left[ V_S^2\rho_{S,0} + \rho_{S,0} V_S^2 - 2 V_S\rho_{S,0}V_S\right] \end{align}\end{split}\]

where the coefficient in front of the second part is by our initial assumption very small \(\gamma := \frac{\epsilon t^2}{2\hbar^2}\tr{V_E^2 \rho_{E,0}} \ll 1\). This evolution happens to be approximately equal to a Kraus map with operators \(K_1 := I - \frac{\gamma}{2} V_S^2, K_2:= \sqrt{\gamma} V_S\):

\[\begin{align} \rho_S \to \rho_S' &= K_1\rho K_1^\dagger + K_2\rho K_2^\dagger = \rho - \frac{\gamma}{2}\left[ V_S^2 \rho + \rho V_S^2\right] + \gamma V_S\rho_S V_S + O(\gamma^2) \end{align}\]

This agrees to \(O(\epsilon^{3/2})\) with the result of our derivation above. This type of derivation can be extended to many other cases with little complication and a very similar argument is used to derive the Lindblad master equation.

Noisy Gates on the Rigetti QVM

As of today, users of our Forest SDK can annotate their QUIL programs by certain pragma statements that inform the QVM that a particular gate on specific target qubits should be replaced by an imperfect realization given by a Kraus map.

The QVM propagates pure states — so how does it simulate noisy gates? It does so by yielding the correct outcomes in the average over many executions of the QUIL program: When the noisy version of a gate should be applied the QVM makes a random choice which Kraus operator is applied to the current state with a probability that ensures that the average over many executions is equivalent to the Kraus map. In particular, a particular Kraus operator \(K_j\) is applied to \(\ket{\psi}_S\)

\[\ket{\psi'}_S = \frac{1}{\sqrt{p_j}} K_j \ket{\psi}_S\]

with probability \(p_j:= \bra{\psi}_S K_j^\dagger K_j \ket{\psi}_S\). In the average over many execution \(N \gg 1\) we therefore find that

\[\begin{split}\begin{align} \overline{\rho_S'} & = \frac{1}{N} \sum_{n=1}^N \ket{\psi'_n}_S\bra{\psi'_n}_S \\ & = \frac{1}{N} \sum_{n=1}^N p_{j_n}^{-1}K_{j_n}\ket{\psi'}_S \bra{\psi'}_SK_{j_n}^\dagger \end{align}\end{split}\]

where \(j_n\) is the chosen Kraus operator label in the \(n\)-th trial. This is clearly a Kraus map itself! And we can group identical terms and rewrite it as

\[\begin{align} \overline{\rho_S'} & = \sum_{\ell=1}^n \frac{N_\ell}{N} p_{\ell}^{-1}K_{\ell}\ket{\psi'}_S \bra{\psi'}_SK_{\ell}^\dagger \end{align}\]

where \(N_{\ell}\) is the number of times that Kraus operator label \(\ell\) was selected. For large enough \(N\) we know that \(N_{\ell} \approx N p_\ell\) and therefore

\[\begin{align} \overline{\rho_S'} \approx \sum_{\ell=1}^n K_{\ell}\ket{\psi'}_S \bra{\psi'}_SK_{\ell}^\dagger \end{align}\]

which proves our claim. The consequence is that noisy gate simulations must generally be repeated many times to obtain representative results.

Getting Started

  1. Come up with a good model for your noise. We will provide some examples below and may add more such examples to our public repositories over time. Alternatively, you can characterize the gate under consideration using Quantum Process Tomography or Gate Set Tomography and use the resulting process matrices to obtain a very accurate noise model for a particular QPU.

  2. Define your Kraus operators as a list of numpy arrays kraus_ops = [K1, K2, ..., Km].

  3. For your QUIL program p, call:

    p.define_noisy_gate("MY_NOISY_GATE", [q1, q2], kraus_ops)

    where you should replace MY_NOISY_GATE with the gate of interest and q1, q2 the indices of the qubits.

Scroll down for some examples!

from __future__ import print_function
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import binom
import matplotlib.colors as colors
%matplotlib inline

from pyquil.quil import Program, MEASURE
from pyquil.api import QVMConnection
from pyquil.gates import CZ, H, I, X
from scipy.linalg import expm

cxn = QVMConnection()
Example 1: Amplitude Damping

Amplitude damping channels are imperfect identity maps with Kraus operators

\[\begin{split}K_1 = \begin{pmatrix} 1 & 0 \\ 0 & \sqrt{1-p} \end{pmatrix} \\ K_2 = \begin{pmatrix} 0 & \sqrt{p} \\ 0 & 0 \end{pmatrix}\end{split}\]

where \(p\) is the probability that a qubit in the \(\ket{1}\) state decays to the \(\ket{0}\) state.

def damping_channel(damp_prob=.1):
    """
    Generate the Kraus operators corresponding to an amplitude damping
    noise channel.

    :params float damp_prob: The one-step damping probability.
    :return: A list [k1, k2] of the Kraus operators that parametrize the map.
    :rtype: list
    """
    damping_op = np.sqrt(damp_prob) * np.array([[0, 1],
                                                [0, 0]])

    residual_kraus = np.diag([1, np.sqrt(1-damp_prob)])
    return [residual_kraus, damping_op]

def append_kraus_to_gate(kraus_ops, g):
    """
    Follow a gate `g` by a Kraus map described by `kraus_ops`.

    :param list kraus_ops: The Kraus operators.
    :param numpy.ndarray g: The unitary gate.
    :return: A list of transformed Kraus operators.
    """
    return [kj.dot(g) for kj in kraus_ops]


def append_damping_to_gate(gate, damp_prob=.1):
    """
    Generate the Kraus operators corresponding to a given unitary
    single qubit gate followed by an amplitude damping noise channel.

    :params np.ndarray|list gate: The 2x2 unitary gate matrix.
    :params float damp_prob: The one-step damping probability.
    :return: A list [k1, k2] of the Kraus operators that parametrize the map.
    :rtype: list
    """
    return append_kraus_to_gate(damping_channel(damp_prob), gate)

%%time

# single step damping probability
damping_per_I = 0.02

# number of program executions
trials = 200

results = []
outcomes = []
lengths = np.arange(0, 201, 10, dtype=int)
for jj, num_I in enumerate(lengths):

    print("{}/{}, ".format(jj, len(lengths)), end="")


    p = Program(X(0))
    # want increasing number of I-gates
    p.inst([I(0) for _ in range(num_I)])
    p.inst(MEASURE(0, [0]))

    # overload identity I on qc 0
    p.define_noisy_gate("I", [0], append_damping_to_gate(np.eye(2), damping_per_I))
    cxn.random_seed = int(num_I)
    res = cxn.run(p, [0], trials=trials)
    results.append([np.mean(res), np.std(res) / np.sqrt(trials)])

results = np.array(results)

0/21, 1/21, 2/21, 3/21, 4/21, 5/21, 6/21, 7/21, 8/21, 9/21, 10/21, 11/21, 12/21, 13/21, 14/21, 15/21, 16/21, 17/21, 18/21, 19/21, 20/21, CPU times: user 138 ms, sys: 19.2 ms, total: 157 ms
Wall time: 6.4 s

dense_lengths = np.arange(0, lengths.max()+1, .2)
survival_probs = (1-damping_per_I)**dense_lengths
logpmf = binom.logpmf(np.arange(trials+1)[np.newaxis, :], trials, survival_probs[:, np.newaxis])/np.log(10)

DARK_TEAL = '#48737F'
FUSCHIA = "#D6619E"
BEIGE = '#EAE8C6'
cm = colors.LinearSegmentedColormap.from_list('anglemap', ["white", FUSCHIA, BEIGE], N=256, gamma=1.5)

plt.figure(figsize=(14, 6))
plt.pcolor(dense_lengths, np.arange(trials+1)/trials, logpmf.T, cmap=cm, vmin=-4, vmax=logpmf.max())
plt.plot(dense_lengths, survival_probs, c=BEIGE, label="Expected mean")
plt.errorbar(lengths, results[:,0], yerr=2*results[:,1], c=DARK_TEAL,
             label=r"noisy qvm, errorbars $ = \pm 2\hat{\sigma}$", marker="o")
cb = plt.colorbar()
cb.set_label(r"$\log_{10} \mathrm{Pr}(n_1; n_{\rm trials}, p_{\rm survival}(t))$", size=20)

plt.title("Amplitude damping model of a single qubit", size=20)
plt.xlabel(r"Time $t$ [arb. units]", size=14)
plt.ylabel(r"$n_1/n_{\rm trials}$", size=14)
plt.legend(loc="best", fontsize=18)
plt.xlim(*lengths[[0, -1]])
plt.ylim(0, 1)
_images/GateNoiseModels_14_1.png
Example 2: Dephased CZ-gate

Dephasing is usually characterized through a qubit’s \(T_2\) time. For a single qubit the dephasing Kraus operators are

\[\begin{split}K_1(p) = \sqrt{1-p} I_2 \\ K_2(p) = \sqrt{p} \sigma_Z\end{split}\]

where \(p = 1 - \exp(-T_2/T_{\rm gate})\) is the probability that the qubit is dephased over the time interval of interest, \(I_2\) is the \(2\times 2\)-identity matrix and \(\sigma_Z\) is the Pauli-Z operator.

For two qubits, we must construct a Kraus map that has four different outcomes:

  1. No dephasing
  2. Qubit 1 dephases
  3. Qubit 2 dephases
  4. Both dephase

The Kraus operators for this are given by

where we assumed a dephasing probability \(p\) for the first qubit and \(q\) for the second.

Dephasing is a diagonal error channel and the CZ gate is also diagonal, therefore we can get the combined map of dephasing and the CZ gate simply by composing \(U_{\rm CZ}\) the unitary representation of CZ with each Kraus operator

\[\begin{split}\begin{align} K^{\rm CZ}_1(p,q) = K_1(p)\otimes K_1(q)U_{\rm CZ} \\ K^{\rm CZ}_2(p,q) = K_2(p)\otimes K_1(q)U_{\rm CZ} \\ K^{\rm CZ}_3(p,q) = K_1(p)\otimes K_2(q)U_{\rm CZ} \\ K^{\rm CZ}_4(p,q) = K_2(p)\otimes K_2(q)U_{\rm CZ} \end{align}\end{split}\]

Note that this is not always accurate, because a CZ gate is often achieved through non-diagonal interaction Hamiltonians! However, for sufficiently small dephasing probabilities it should always provide a good starting point.

def dephasing_kraus_map(p=.1):
    """
    Generate the Kraus operators corresponding to a dephasing channel.

    :params float p: The one-step dephasing probability.
    :return: A list [k1, k2] of the Kraus operators that parametrize the map.
    :rtype: list
    """
    return [np.sqrt(1-p)*np.eye(2), np.sqrt(p)*np.diag([1, -1])]

def tensor_kraus_maps(k1, k2):
    """
    Generate the Kraus map corresponding to the composition
    of two maps on different qubits.

    :param list k1: The Kraus operators for the first qubit.
    :param list k2: The Kraus operators for the second qubit.
    :return: A list of tensored Kraus operators.
    """
    return [np.kron(k1j, k2l) for k1j in k1 for k2l in k2]


def append_kraus_to_gate(kraus_ops, g):
    """
    Follow a gate `g` by a Kraus map described by `kraus_ops`.

    :param list kraus_ops: The Kraus operators.
    :param numpy.ndarray g: The unitary gate.
    :return: A list of transformed Kraus operators.
    """
    return [kj.dot(g) for kj in kraus_ops]

%%time
# single step damping probabilities
ps = np.linspace(.001, .5, 200)

# number of program executions
trials = 500

results = []

for jj, p in enumerate(ps):

    corrupted_CZ = append_kraus_to_gate(
    tensor_kraus_maps(
        dephasing_kraus_map(p),
        dephasing_kraus_map(p)
    ),
    np.diag([1, 1, 1, -1]))


    print("{}/{}, ".format(jj, len(ps)), end="")

    # make Bell-state
    p = Program(H(0), H(1), CZ(0,1), H(1))

    p.inst(MEASURE(0, [0]))
    p.inst(MEASURE(1, [1]))

    # overload identity I on qc 0
    p.define_noisy_gate("CZ", [0, 1], corrupted_CZ)
    cxn.random_seed = jj
    res = cxn.run(p, [0, 1], trials=trials)
    results.append(res)

results = np.array(results)

0/200, 1/200, 2/200, 3/200, 4/200, 5/200, 6/200, 7/200, 8/200, 9/200, 10/200, 11/200, 12/200, 13/200, 14/200, 15/200, 16/200, 17/200, 18/200, 19/200, 20/200, 21/200, 22/200, 23/200, 24/200, 25/200, 26/200, 27/200, 28/200, 29/200, 30/200, 31/200, 32/200, 33/200, 34/200, 35/200, 36/200, 37/200, 38/200, 39/200, 40/200, 41/200, 42/200, 43/200, 44/200, 45/200, 46/200, 47/200, 48/200, 49/200, 50/200, 51/200, 52/200, 53/200, 54/200, 55/200, 56/200, 57/200, 58/200, 59/200, 60/200, 61/200, 62/200, 63/200, 64/200, 65/200, 66/200, 67/200, 68/200, 69/200, 70/200, 71/200, 72/200, 73/200, 74/200, 75/200, 76/200, 77/200, 78/200, 79/200, 80/200, 81/200, 82/200, 83/200, 84/200, 85/200, 86/200, 87/200, 88/200, 89/200, 90/200, 91/200, 92/200, 93/200, 94/200, 95/200, 96/200, 97/200, 98/200, 99/200, 100/200, 101/200, 102/200, 103/200, 104/200, 105/200, 106/200, 107/200, 108/200, 109/200, 110/200, 111/200, 112/200, 113/200, 114/200, 115/200, 116/200, 117/200, 118/200, 119/200, 120/200, 121/200, 122/200, 123/200, 124/200, 125/200, 126/200, 127/200, 128/200, 129/200, 130/200, 131/200, 132/200, 133/200, 134/200, 135/200, 136/200, 137/200, 138/200, 139/200, 140/200, 141/200, 142/200, 143/200, 144/200, 145/200, 146/200, 147/200, 148/200, 149/200, 150/200, 151/200, 152/200, 153/200, 154/200, 155/200, 156/200, 157/200, 158/200, 159/200, 160/200, 161/200, 162/200, 163/200, 164/200, 165/200, 166/200, 167/200, 168/200, 169/200, 170/200, 171/200, 172/200, 173/200, 174/200, 175/200, 176/200, 177/200, 178/200, 179/200, 180/200, 181/200, 182/200, 183/200, 184/200, 185/200, 186/200, 187/200, 188/200, 189/200, 190/200, 191/200, 192/200, 193/200, 194/200, 195/200, 196/200, 197/200, 198/200, 199/200, CPU times: user 1.17 s, sys: 166 ms, total: 1.34 s
Wall time: 1min 49s

Z1s = (2*results[:,:,0]-1.)
Z2s = (2*results[:,:,1]-1.)
Z1Z2s = Z1s * Z2s

Z1m = np.mean(Z1s, axis=1)
Z2m = np.mean(Z2s, axis=1)
Z1Z2m = np.mean(Z1Z2s, axis=1)

plt.figure(figsize=(14, 6))
plt.axhline(y=1.0, color=FUSCHIA, alpha=.5, label="Bell state")

plt.plot(ps, Z1Z2m, "x", c=FUSCHIA, label=r"$\overline{Z_1 Z_2}$")
plt.plot(ps, 1-2*ps, "--", c=FUSCHIA, label=r"$\langle Z_1 Z_2\rangle_{\rm theory}$")

plt.plot(ps, Z1m, "o", c=DARK_TEAL, label=r"$\overline{Z}_1$")
plt.plot(ps, 0*ps, "--", c=DARK_TEAL, label=r"$\langle Z_1\rangle_{\rm theory}$")

plt.plot(ps, Z2m, "d", c="k", label=r"$\overline{Z}_2$")
plt.plot(ps, 0*ps, "--", c="k", label=r"$\langle Z_2\rangle_{\rm theory}$")

plt.xlabel(r"Dephasing probability $p$", size=18)
plt.ylabel(r"$Z$-moment", size=18)
plt.title(r"$Z$-moments for a Bell-state prepared with dephased CZ", size=18)
plt.xlim(0, .5)
plt.legend(fontsize=18)
_images/GateNoiseModels_20_1.png

Adding Decoherence Noise

In this example, we investigate how a program might behave on a near-term device that is subject to T1- and T2-type noise using the convenience function pyquil.noise.add_decoherence_noise(). The same module also contains some other useful functions to define your own types of noise models, e.g., pyquil.noise.tensor_kraus_maps() for generating multi-qubit noise processes, pyquil.noise.combine_kraus_maps() for describing the succession of two noise processes and pyquil.noise.append_kraus_to_gate() which allows appending a noise process to a unitary gate.

from pyquil.quil import Program
from pyquil.paulis import PauliSum, PauliTerm, exponentiate, exponential_map, trotterize
from pyquil.gates import MEASURE, H, Z, RX, RZ, CZ
import numpy as np
The Task

We want to prepare \(e^{i \theta XY}\) and measure it in the \(Z\) basis.

from numpy import pi
theta = pi/3
xy = PauliTerm('X', 0) * PauliTerm('Y', 1)
The Idiomatic PyQuil Program
prog = exponential_map(xy)(theta)
print(prog)

H 0
RX(pi/2) 1
CNOT 0 1
RZ(2*pi/3) 1
CNOT 0 1
H 0
RX(-pi/2) 1
The Compiled Program

To run on a real device, we must compile each program to the native gate set for the device. The high-level noise model is similarly constrained to use a small, native gate set. In particular, we can use

  • \(I\)
  • \(RZ(\theta)\)
  • \(RX(\pm \pi/2)\)
  • \(CZ\)

For simplicity, the compiled program is given below but generally you will want to use a compiler to do this step for you.

def get_compiled_prog(theta):
    return Program([
        RZ(-pi/2, 0),
        RX(-pi/2, 0),
        RZ(-pi/2, 1),
        RX( pi/2, 1),
        CZ(1, 0),
        RZ(-pi/2, 1),
        RX(-pi/2, 1),
        RZ(theta, 1),
        RX( pi/2, 1),
        CZ(1, 0),
        RX( pi/2, 0),
        RZ( pi/2, 0),
        RZ(-pi/2, 1),
        RX( pi/2, 1),
        RZ(-pi/2, 1),
    ])
Scan Over Noise Parameters

We perform a scan over three levels of noise each at 20 theta points.

Specifically, we investigate T1 values of 1, 3, and 10 us. By default, T2 = T1 / 2, 1 qubit gates take 50 ns, and 2 qubit gates take 150 ns.

In alignment with the device, \(I\) and parametric \(RZ\) are noiseless while \(RX\) and \(CZ\) gates experience 1q and 2q gate noise, respectively.

from pyquil.api import QVMConnection
cxn = QVMConnection()

t1s = np.logspace(-6, -5, num=3)
thetas = np.linspace(-pi, pi, num=20)
t1s * 1e6 # us

array([  1.        ,   3.16227766,  10.        ])

from pyquil.noise import add_decoherence_noise
records = []
for theta in thetas:
    for t1 in t1s:
        prog = get_compiled_prog(theta)
        noisy = add_decoherence_noise(prog, T1=t1).inst([
            MEASURE(0, 0),
            MEASURE(1, 1),
        ])
        bitstrings = np.array(cxn.run(noisy, [0,1], 1000))

        # Expectation of Z0 and Z1
        z0, z1 = 1 - 2*np.mean(bitstrings, axis=0)

        # Expectation of ZZ by computing the parity of each pair
        zz = 1 - (np.sum(bitstrings, axis=1) % 2).mean() * 2

        record = {
            'z0': z0,
            'z1': z1,
            'zz': zz,
            'theta': theta,
            't1': t1,
        }
        records += [record]
Plot the Results

Note that to run the code below you will need to install the pandas and seaborn packages.

%matplotlib inline
from matplotlib import pyplot as plt
import seaborn as sns
sns.set(style='ticks', palette='colorblind')

import pandas as pd
df_all = pd.DataFrame(records)
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(12,4))

for t1 in t1s:
    df = df_all.query('t1 == @t1')

    ax1.plot(df['theta'], df['z0'], 'o-')
    ax2.plot(df['theta'], df['z1'], 'o-')
    ax3.plot(df['theta'], df['zz'], 'o-', label='T1 = {:.0f} us'.format(t1*1e6))

ax3.legend(loc='best')

ax1.set_ylabel('Z0')
ax2.set_ylabel('Z1')
ax3.set_ylabel('ZZ')
ax2.set_xlabel(r'$\theta$')
fig.tight_layout()
_images/NoiseInvestigation_14_0.png

Modeling Readout Noise

Qubit-Readout can be corrupted in a variety of ways. The two most relevant error mechanisms on the Rigetti QPU right now are:

  1. Transmission line noise that makes a 0-state look like a 1-state or vice versa. We call this classical readout bit-flip error. This type of readout noise can be reduced by tailoring optimal readout pulses and using superconducting, quantum limited amplifiers to amplify the readout signal before it is corrupted by classical noise at the higher temperature stages of our cryostats.
  2. T1 qubit decay during readout (our readout operations can take more than a µsecond unless they have been specially optimized), which leads to readout signals that initially behave like 1-states but then collapse to something resembling a 0-state. We will call this T1-readout error. This type of readout error can be reduced by achieving shorter readout pulses relative to the T1 time, i.e., one can try to reduce the readout pulse length, or increase the T1 time or both.
Qubit Measurements

This section provides the necessary theoretical foundation for accurately modeling noisy quantum measurements on superconducting quantum processors. It relies on some of the abstractions (density matrices, Kraus maps) introduced in our notebook on gate noise models.

The most general type of measurement performed on a single qubit at a single time can be characterized by some set \(\mathcal{O}\) of measurement outcomes, e.g., in the simplest case \(\mathcal{O} = \{0, 1\}\), and some unnormalized quantum channels (see notebook on gate noise models) that encapsulate 1. the probability of that outcome 2. how the qubit state is affected conditional on the measurement outcome.

Here the outcome is understood as classical information that has been extracted from the quantum system.

Projective, Ideal Measurement

The simplest case that is usually taught in introductory quantum mechanics and quantum information courses are Born’s rule and the projection postulate which state that there exist a complete set of orthogonal projection operators

\[P_{\mathcal{O}} := \{\Pi_x \text{ Projector }\mid x \in \mathcal{O}\},\]

i.e., one for each measurement outcome. Any projection operator must satisfy \(\Pi_x^\dagger = \Pi_x = \Pi_x^2\) and for an orthogonal set of projectors any two members satisfy

\[\begin{split}\Pi_x\Pi_y = \delta_{xy} \Pi_x = \begin{cases} 0 & \text{ if } x \ne y \\ \Pi_x & \text{ if } x = y \end{cases}\end{split}\]

and for a complete set we additionally demand that \(\sum_{x\in\mathcal{O}} \Pi_x = 1\). Following our introduction to gate noise, we write quantum states as density matrices as this is more general and in closer correspondence with classical probability theory.

With these the probability of outcome \(x\) is given by \(p(x) = \tr{\Pi_x \rho \Pi_x} = \tr{\Pi_x^2 \rho} = \tr{\Pi_x \rho}\) and the post measurement state is

\[\rho_x = \frac{1}{p(x)} \Pi_x \rho \Pi_x,\]

which is the projection postulate applied to mixed states.

If we were a sloppy quantum programmer and accidentally erased the measurement outcome then our best guess for the post measurement state would be given by something that looks an awful lot like a Kraus map:

\[\rho_{\text{post measurement}} = \sum_{x\in\mathcal{O}} p(x) \rho_x = \sum_{x\in\mathcal{O}} \Pi_x \rho \Pi_x.\]

The completeness of the projector set ensures that the trace of the post measurement is still 1 and the Kraus map form of this expression ensures that \(\rho_{\text{post measurement}}\) is a positive (semi-)definite operator.

Classical Readout Bit-Flip Error

Consider now the ideal measurement as above, but where the outcome \(x\) is transmitted across a noisy classical channel that produces a final outcome \(x'\in \mathcal{O}' = \{0', 1'\}\) according to some conditional probabilities \(p(x'|x)\) that can be recorded in the assignment probability matrix

\[\begin{split}P_{x'|x} = \begin{pmatrix} p(0 | 0) & p(0 | 1) \\ p(1 | 0) & p(1 | 1) \end{pmatrix}\end{split}\]

Note that this matrix has only two independent parameters as each column must be a valid probability distribution, i.e. all elements are non-negative and each column sums to 1.

This matrix allows us to obtain the probabilities \(\mathbf{p}' := (p(x'=0), p(x'=1))^T\) from the original outcome probabilities \(\mathbf{p} := (p(x=0), p(x=1))^T\) via \(\mathbf{p}' = P_{x'|x}\mathbf{p}\). The difference relative to the ideal case above is that now an outcome \(x' = 0\) does not necessarily imply that the post measurement state is truly \(\Pi_{0} \rho \Pi_{0} / p(x=0)\). Instead, the post measurement state given a noisy outcome \(x'\) must be

\[\begin{split}\begin{aligned} \rho_{x'} & = \sum_{x\in \mathcal{O}} p(x|x') \rho_x \\ & = \sum_{x\in \mathcal{O}} p(x'|x)\frac{p(x)}{p(x')} \rho_x \\ & = \frac{1}{p(x')}\sum_{x\in \mathcal{O}} p(x'|x) \Pi_x \rho \Pi_x \end{aligned}\end{split}\]

where

\[\begin{split}\begin{aligned} p(x') & = \sum_{x\in\mathcal{O}} p(x'|x) p(x) \\ & = \tr{\sum_{x\in \mathcal{O}} p(x'|x) \Pi_x \rho \Pi_x} \\ & = \tr{\rho \sum_{x\in \mathcal{O}} p(x'|x)\Pi_x} \\ & = \tr{\rho E_x'}. \end{aligned}\end{split}\]

where we have exploited the cyclical property of the trace \(\tr{ABC}=\tr{BCA}\) and the projection property \(\Pi_x^2 = \Pi_x\). This has allowed us to derive the noisy outcome probabilities from a set of positive operators

\[E_{x'} := \sum_{x\in \mathcal{O}} p(x'|x)\Pi_x \ge 0\]

that must sum to 1:

\[\sum_{x'\in\mathcal{O}'} E_x' = \sum_{x\in\mathcal{O}}\underbrace{\left[\sum_{x'\in\mathcal{O}'} p(x'|x)\right]}_{1}\Pi_x = \sum_{x\in\mathcal{O}}\Pi_x = 1.\]

The above result is a type of generalized Bayes’ theorem that is extremely useful for this type of (slightly) generalized measurement and the family of operators \(\{E_{x'}| x' \in \mathcal{O}'\}\) whose expectations give the probabilities is called a positive operator valued measure (POVM). These operators are not generally orthogonal nor valid projection operators but they naturally arise in this scenario. This is not yet the most general type of measurement, but it will get us pretty far.

How to Model \(T_1\) Error

T1 type errors fall outside our framework so far as they involve a scenario in which the quantum state itself is corrupted during the measurement process in a way that potentially erases the pre-measurement information as opposed to a loss of purely classical information. The most appropriate framework for describing this is given by that of measurement instruments, but for the practical purpose of arriving at a relatively simple description, we propose describing this by a T1 damping Kraus map followed by the noisy readout process as described above.

Further Reading

Chapter 3 of John Preskill’s lecture notes http://www.theory.caltech.edu/people/preskill/ph229/notes/chap3.pdf

Working with Readout Noise

  1. Come up with a good guess for your readout noise parameters \(p(0|0)\) and \(p(1|1)\), the off-diagonals then follow from the normalization of \(P_{x'|x}\). If your assignment fidelity \(F\) is given, and you assume that the classical bit flip noise is roughly symmetric, then a good approximation is to set \(p(0|0)=p(1|1)=F\).

  2. For your QUIL program p, and a qubit index q call:

    p.define_noisy_readout(q, p00, p11)

    where you should replace p00 and p11 with the assumed probabilities.

Scroll down for some examples!

from __future__ import print_function, division
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

from pyquil.quil import Program, MEASURE, Pragma
from pyquil.api.qvm import QVMConnection
from pyquil.gates import I, X, RX, H, CNOT
from pyquil.noise import (estimate_bitstring_probs, correct_bitstring_probs,
                          bitstring_probs_to_z_moments, estimate_assignment_probs)

DARK_TEAL = '#48737F'
FUSCHIA = '#D6619E'
BEIGE = '#EAE8C6'

cxn = QVMConnection()
Example 1: Rabi Sequence with Noisy Readout
%%time

# number of angles
num_theta = 101

# number of program executions
trials = 200

thetas = np.linspace(0, 2*np.pi, num_theta)

p00s = [1., 0.95, 0.9, 0.8]

results_rabi = np.zeros((num_theta, len(p00s)))

for jj, theta in enumerate(thetas):
    for kk, p00 in enumerate(p00s):
        cxn.random_seed = hash((jj, kk))
        p = Program(RX(theta, 0))
        # assume symmetric noise p11 = p00
        p.define_noisy_readout(0, p00=p00, p11=p00)
        p.measure(0, 0)
        res = cxn.run(p, [0], trials=trials)
        results_rabi[jj, kk] = np.sum(res)

CPU times: user 1.2 s, sys: 73.6 ms, total: 1.27 s
Wall time: 3.97 s

plt.figure(figsize=(14, 6))
for jj, (p00, c) in enumerate(zip(p00s, [DARK_TEAL, FUSCHIA, "k", "gray"])):
    plt.plot(thetas, results_rabi[:, jj]/trials, c=c, label=r"$p(0|0)=p(1|1)={:g}$".format(p00))
plt.legend(loc="best")
plt.xlim(*thetas[[0,-1]])
plt.ylim(-.1, 1.1)
plt.grid(alpha=.5)
plt.xlabel(r"RX angle $\theta$ [radian]", size=16)
plt.ylabel(r"Excited state fraction $n_1/n_{\rm trials}$", size=16)
plt.title("Effect of classical readout noise on Rabi contrast.", size=18)

<matplotlib.text.Text at 0x104314250>
_images/ReadoutNoise_10_1.png
Example 2: Estimate the Assignment Probabilities

Here we will estimate \(P_{x'|x}\) ourselves! You can run some simple experiments to estimate the assignment probability matrix directly from a QPU.

On a perfect quantum computer

estimate_assignment_probs(0, 1000, cxn, Program())

array([[ 1.,  0.],
       [ 0.,  1.]])

On an imperfect quantum computer

cxn.seed = None
header0 = Program().define_noisy_readout(0, .85, .95)
header1 = Program().define_noisy_readout(1, .8, .9)
header2 = Program().define_noisy_readout(2, .9, .85)

ap0 = estimate_assignment_probs(0, 100000, cxn, header0)
ap1 = estimate_assignment_probs(1, 100000, cxn, header1)
ap2 = estimate_assignment_probs(2, 100000, cxn, header2)

print(ap0, ap1, ap2, sep="\n")

[[ 0.84967  0.04941]
 [ 0.15033  0.95059]]
[[ 0.80058  0.09993]
 [ 0.19942  0.90007]]
[[ 0.90048  0.14988]
 [ 0.09952  0.85012]]
Example 3: Correct for Noisy Readout
3a) Correcting the Rabi Signal from Above
ap_last = np.array([[p00s[-1], 1 - p00s[-1]],
                    [1 - p00s[-1], p00s[-1]]])
corrected_last_result = [correct_bitstring_probs([1-p, p], [ap_last])[1] for p in results_rabi[:, -1] / trials]

plt.figure(figsize=(14, 6))
for jj, (p00, c) in enumerate(zip(p00s, [DARK_TEAL, FUSCHIA, "k", "gray"])):
    if jj not in [0, 3]:
        continue
    plt.plot(thetas, results_rabi[:, jj]/trials, c=c, label=r"$p(0|0)=p(1|1)={:g}$".format(p00), alpha=.3)
plt.plot(thetas, corrected_last_result, c="red", label=r"Corrected $p(0|0)=p(1|1)={:g}$".format(p00s[-1]))
plt.legend(loc="best")
plt.xlim(*thetas[[0,-1]])
plt.ylim(-.1, 1.1)
plt.grid(alpha=.5)
plt.xlabel(r"RX angle $\theta$ [radian]", size=16)
plt.ylabel(r"Excited state fraction $n_1/n_{\rm trials}$", size=16)
plt.title("Corrected contrast", size=18)

<matplotlib.text.Text at 0x1055e7310>
_images/ReadoutNoise_19_1.png

We find that the corrected signal is fairly noisy (and sometimes exceeds the allowed interval \([0,1]\)) due to the overall very small number of samples \(n=200\).

3b) Corrupting and Correcting GHZ State Correlations

In this example we will create a GHZ state \(\frac{1}{\sqrt{2}}\left[\left|000\right\rangle + \left|111\right\rangle \right]\) and measure its outcome probabilities with and without the above noise model. We will then see how the Pauli-Z moments that indicate the qubit correlations are corrupted (and corrected) using our API.

ghz_prog = Program(H(0), CNOT(0, 1), CNOT(1, 2),
                   MEASURE(0, 0), MEASURE(1, 1), MEASURE(2, 2))
print(ghz_prog)
results = cxn.run(ghz_prog, [0, 1, 2], trials=10000)

H 0
CNOT 0 1
CNOT 1 2
MEASURE 0 [0]
MEASURE 1 [1]
MEASURE 2 [2]

header = header0 + header1 + header2
noisy_ghz = header + ghz_prog
print(noisy_ghz)
noisy_results = cxn.run(noisy_ghz, [0, 1, 2], trials=10000)

PRAGMA READOUT-POVM 0 "(0.85 0.050000000000000044 0.15000000000000002 0.95)"
PRAGMA READOUT-POVM 1 "(0.8 0.09999999999999998 0.19999999999999996 0.9)"
PRAGMA READOUT-POVM 2 "(0.9 0.15000000000000002 0.09999999999999998 0.85)"
H 0
CNOT 0 1
CNOT 1 2
MEASURE 0 [0]
MEASURE 1 [1]
MEASURE 2 [2]
Uncorrupted probability for \(\left|000\right\rangle\) and \(\left|111\right\rangle\) 
probs = estimate_bitstring_probs(results)
probs[0, 0, 0], probs[1, 1, 1]

(0.50419999999999998, 0.49580000000000002)

As expected the outcomes 000 and 111 each have roughly probability \(1/2\).

Corrupted probability for \(\left|011\right\rangle\) and \(\left|100\right\rangle\) 
noisy_probs = estimate_bitstring_probs(noisy_results)
noisy_probs[0, 0, 0], noisy_probs[1, 1, 1]

(0.30869999999999997, 0.3644)

The noise-corrupted outcome probabilities deviate significantly from their ideal values!

Corrected probability for \(\left|011\right\rangle\) and \(\left|100\right\rangle\) 
corrected_probs = correct_bitstring_probs(noisy_probs, [ap0, ap1, ap2])
corrected_probs[0, 0, 0], corrected_probs[1, 1, 1]

(0.50397601453064977, 0.49866843912900716)

The corrected outcome probabilities are much closer to the ideal value.

Estimate \(\langle Z_0^{j} Z_1^{k} Z_2^{\ell}\rangle\) for \(jkl=100, 010, 001\) from non-noisy data

We expect these to all be very small

zmoments = bitstring_probs_to_z_moments(probs)
zmoments[1, 0, 0], zmoments[0, 1, 0], zmoments[0, 0, 1]

(0.0083999999999999631, 0.0083999999999999631, 0.0083999999999999631)
Estimate \(\langle Z_0^{j} Z_1^{k} Z_2^{\ell}\rangle\) for \(jkl=110, 011, 101\) from non-noisy data

We expect these to all be close to 1.

zmoments[1, 1, 0], zmoments[0, 1, 1], zmoments[1, 0, 1]

(1.0, 1.0, 1.0)
Estimate \(\langle Z_0^{j} Z_1^{k} Z_2^{\ell}\rangle\) for \(jkl=100, 010, 001\) from noise-corrected data 
zmoments_corr = bitstring_probs_to_z_moments(corrected_probs)
zmoments_corr[1, 0, 0], zmoments_corr[0, 1, 0], zmoments_corr[0, 0, 1]

(0.0071476770049732075, -0.0078641261685578612, 0.0088462563282706852)
Estimate \(\langle Z_0^{j} Z_1^{k} Z_2^{\ell}\rangle\) for \(jkl=110, 011, 101\) from noise-corrected data 
zmoments_corr[1, 1, 0], zmoments_corr[0, 1, 1], zmoments_corr[1, 0, 1]

(0.99477496902638118, 1.0008376440216553, 1.0149652015905912)

Overall the correction can restore the contrast in our multi-qubit observables, though we also see that the correction can lead to slightly non-physical expectations. This effect is reduced the more samples we take.

Advanced Usage

Note

If you’re running locally, remember set up the QVM and quilc in server mode before trying to use them: Setting Up Server Mode for PyQuil.

PyQuil Configuration Files

Network endpoints for the Rigetti Forest infrastructure and information pertaining to QPU access are stored in a pair of configuration files. These files are located by default at ~/.qcs_config and ~/.forest_config. The location can be changed by setting the environment variables QCS_CONFIG or FOREST_CONFIG to point to the new location.

When running on a QMI, the values in these configuration files are automatically managed so as to point to the correct endpoints. When running locally, configuration files are not necessary. Thus, the average user will not have to do any work to get their configuration files set up.

If for some reason you want to use an atypical configuration, you may need to modify these files.

The default QCS config file on any QMI looks similar to the following:

# .qcs_config
[Rigetti Forest]
url = https://forest-server.qcs.rigetti.com
key = 4fd12391-11eb-52ec-35c2-262765ae4c4f
user_id = 4fd12391-11eb-52ec-35c2-262765ae4c4f

[QPU]
exec_on_engage = bash exec_on_engage.sh

where

  • url is the endpoint that pyQuil hits for device information and for the 2.0 endpoints,
  • key stores the Forest 1.X API key,
  • user_id stores a Forest 2.0 user ID, and
  • exec_on_engage specifies the shell command that the QMI will launch when the QMI becomes QPU-engaged. It would have no effect if you are running locally, but is important if you are running on the QMI. By default, it runs the exec_on_engage.sh shell script. It’s best to leave the configuration as is, and edit that script. More documentation about exec_on_engage.sh can be found in the QCS docs here.

The Forest config file on any QMI has these contents, with specific IP addresses filled in:

# .forest_config
[Rigetti Forest]
qpu_endpoint_address = None
qvm_address = http://10.1.165.XX:5000
compiler_server_address = tcp://10.1.165.XX:5555

where

  • qpu_endpoint_address is the endpoint where pyQuil will try to communicate with the QPU orchestrating service during QPU-engagement. It may not appear until your QMI engages, and furthermore will have no effect if you are running locally. It’s best to leave this alone. If you obtain access to one of our QPUs, we will fill it in for you.
  • qvm_address is the endpoint where pyQuil will try to communicate with the Rigetti Quantum Virtual Machine. On a QMI, this points to the provided QVM instance. On a local installation, this should be set to the server endpoint for a locally running QVM instance. However, pyQuil will use the default value http://localhost:5000 if this file isn’t found, which is the correct endpoint when you run the QVM locally with qvm -S.
  • compiler_server_address: This is the endpoint where pyQuil will try to communicate with the compiler server. On a QMI, this points to a provided compiler server instance. On a local installation, this should be set to the server endpoint for a locally running quilc instance. However, pyQuil will use the default value http://localhost:6000 if this isn’t set, which is the correct endpoint when you run quilc locally with quilc -S.

Note

PyQuil itself reads these values out using the helper class pyquil._config.PyquilConfig. PyQuil users should not ever need to touch this class directly.

Using Qubit Placeholders

Note

The functionality provided inline by QubitPlaceholders is similar to writing a function which returns a Program, with qubit indices taken as arguments to the function.

In pyQuil, we typically use integers to identify qubits

from pyquil import Program
from pyquil.gates import CNOT, H
print(Program(H(0), CNOT(0, 1)))

H 0
CNOT 0 1

However, when running on real, near-term QPUs we care about what particular physical qubits our program will run on. In fact, we may want to run the same program on an assortment of different qubits. This is where using QubitPlaceholders comes in.

from pyquil.quilatom import QubitPlaceholder
q0 = QubitPlaceholder()
q1 = QubitPlaceholder()
p = Program(H(q0), CNOT(q0, q1))
print(p)

H {q4402789176}
CNOT {q4402789176} {q4402789120}

If you try to use this program directly, it will not work

print(p.out())

RuntimeError: Qubit q4402789176 has not been assigned an index

Instead, you must explicitly map the placeholders to physical qubits. By default, the function address_qubits will address qubits from 0 to N.

from pyquil.quil import address_qubits
print(address_qubits(p))

H 0
CNOT 0 1

The real power comes into play when you provide an explicit mapping:

print(address_qubits(prog, qubit_mapping={
    q0: 14,
    q1: 19,
}))

H 14
CNOT 14 19
Register

Usually, your algorithm will use an assortment of qubits. You can use the convenience function QubitPlaceholder.register() to request a list of qubits to build your program.

qbyte = QubitPlaceholder.register(8)
p_evens = Program(H(q) for q in qbyte)
print(address_qubits(p_evens, {q: i*2 for i, q in enumerate(qbyte)}))

H 0
H 2
H 4
H 6
H 8
H 10
H 12
H 14

Classical Control Flow

Note

Classical control flow is not yet supported on the QPU.

Here are a couple quick examples that show how much richer a Quil program can be with classical control flow. In this first example, we create a while loop by following these steps:

  1. Declare a register called flag_register to use as a boolean test for looping.
  2. Initialize this register to 1 program so our while loop will execute. This is often called the loop preamble or loop initialization.
  3. Write the body of the loop in its own Program. This will be a program that applies an \(X\) gate followed by a \(H\) gate on our qubit.
  4. Using the while_do() method to add control flow.
from pyquil import Program
from pyquil.gates import *

# Initialize the Program and declare a 1 bit memory space for our boolean flag
outer_loop = Program()
flag_register = outer_loop.declare('flag_register', 'BIT')

# Set the initial flag value to 1
outer_loop += MOVE(flag_register, 1)

# Define the body of the loop with a new Program
inner_loop = Program()
inner_loop += Program(X(0), H(0))
inner_loop += MEASURE(0, flag_register)

# Run inner_loop in a loop until flag_register is 0
outer_loop.while_do(flag_register, inner_loop)

print(outer_loop)

DECLARE flag_register BIT[1]
MOVE flag_register 1
LABEL @START1
JUMP-UNLESS @END2 flag_register
X 0
H 0
MEASURE 0 flag_register
JUMP @START1
LABEL @END2

Notice that the outer_loop program applied a Quil instruction directly to a classical register. There are several classical commands that can be used in this fashion:

  • NOT which flips a classical bit
  • AND which operates on two classical bits
  • IOR which operates on two classical bits
  • MOVE which moves the value of a classical bit at one classical address into another
  • EXCHANGE which swaps the value of two classical bits

In this next example, we show how to do conditional branching in the form of the traditional if construct as in many programming languages. Much like the last example, we construct programs for each branch of the if, and put it all together by using the if_then() method.

# Declare our memory spaces
branching_prog = Program()
test_register = branching_prog.declare('test_register', 'BIT')
ro = branching_prog.declare('ro', 'BIT')

# Construct each branch of our if-statement. We can have empty branches
# simply by having empty programs.
then_branch = Program(X(0))
else_branch = Program()

# Construct our program so that the result in test_register is equally likely to be a 0 or 1
branching_prog += H(1)
branching_prog += MEASURE(1, test_register)

# Add the conditional branching
branching_prog.if_then(test_register, then_branch, else_branch)

# Measure qubit 0 into our readout register
branching_prog += MEASURE(0, ro)

print(branching_prog)

DECLARE test_register BIT[1]
DECLARE ro BIT[1]
H 1
MEASURE 1 test_register
JUMP-WHEN @THEN1 test_register
JUMP @END2
LABEL @THEN1
X 0
LABEL @END2
MEASURE 0 ro

We can run this program a few times to see what we get in the readout register ro.

from pyquil import get_qc

qc = get_qc("2q-qvm")
branching_prog.wrap_in_numshots_loop(10)
qc.run(branching_prog)

[[1], [1], [1], [0], [1], [0], [0], [1], [1], [0]]

Parametric Depolarizing Noise

The Rigetti QVM has support for emulating certain types of noise models. One such model is parametric Pauli noise, which is defined by a set of 6 probabilities:

  • The probabilities \(P_X\), \(P_Y\), and \(P_Z\) which define respectively the probability of a Pauli \(X\), \(Y\), or \(Z\) gate getting applied to each qubit after every gate application. These probabilities are called the gate noise probabilities.
  • The probabilities \(P_X'\), \(P_Y'\), and \(P_Z'\) which define respectively the probability of a Pauli \(X\), \(Y\), or \(Z\) gate getting applied to the qubit being measured before it is measured. These probabilities are called the measurement noise probabilities.

We can instantiate a noisy QVM by creating a new connection with these probabilities specified.

# 20% chance of a X gate being applied after gate applications and before measurements.
gate_noise_probs = [0.2, 0.0, 0.0]
meas_noise_probs = [0.2, 0.0, 0.0]
noisy_qvm = qvm(gate_noise=gate_noise_probs, measurement_noise=meas_noise_probs)

We can test this by applying an \(X\)-gate and measuring. Nominally, we should always measure 1.

p = Program().inst(X(0)).measure(0, 0)
print("Without Noise: {}".format(qvm.run(p, [0], 10)))
print("With Noise   : {}".format(noisy_qvm.run(p, [0], 10)))

Without Noise: [[1], [1], [1], [1], [1], [1], [1], [1], [1], [1]]
With Noise   : [[0], [0], [0], [0], [0], [1], [1], [1], [1], [0]]

Pauli Operator Algebra

Many algorithms require manipulating sums of Pauli combinations, such as \(\sigma = \frac{1}{2}I - \frac{3}{4}X_0Y_1Z_3 + (5-2i)Z_1X_2,\) where \(G_n\) indicates the gate \(G\) acting on qubit \(n\). We can represent such sums by constructing PauliTerm and PauliSum. The above sum can be constructed as follows:

from pyquil.paulis import ID, sX, sY, sZ

# Pauli term takes an operator "X", "Y", "Z", or "I"; a qubit to act on, and
# an optional coefficient.
a = 0.5 * ID()
b = -0.75 * sX(0) * sY(1) * sZ(3)
c = (5-2j) * sZ(1) * sX(2)

# Construct a sum of Pauli terms.
sigma = a + b + c
print(f"sigma = {sigma}")

sigma = (0.5+0j)*I + (-0.75+0j)*X0*Y1*Z3 + (5-2j)*Z1*X2

Right now, the primary thing one can do with Pauli terms and sums is to construct the exponential of the Pauli term, i.e., \(\exp[-i\beta\sigma]\). This is accomplished by constructing a parameterized Quil program that is evaluated when passed values for the coefficients of the angle \(\beta\).

Related to exponentiating Pauli sums we provide utility functions for finding the commuting subgroups of a Pauli sum and approximating the exponential with the Suzuki-Trotter approximation through fourth order.

When arithmetic is done with Pauli sums, simplification is automatically done.

The following shows an instructive example of all three.

from pyquil.paulis import exponential_map

sigma_cubed = sigma * sigma * sigma
print(f"Simplified: {sigma_cubed}\n")

# Produce Quil code to compute exp[iX]
H = -1.0 * sX(0)
print(f"Quil to compute exp[iX] on qubit 0:\n"
       f"{exponential_map(H)(1.0)}")

Simplified: (32.46875-30j)*I + (-16.734375+15j)*X0*Y1*Z3 + (71.5625-144.625j)*Z1*X2

Quil to compute exp[iX] on qubit 0:
H 0
RZ(-2.0) 0
H 0

exponential_map returns a function allowing you to fill in a multiplicative constant later. This commonly occurs in variational algorithms. The function exponential_map is used to compute \(\exp[-i \alpha H]\) without explicitly filling in a value for \(\alpha\).

expH = exponential_map(H)
print(f"0:\n{expH(0.0)}\n")
print(f"1:\n{expH(1.0)}\n")
print(f"2:\n{expH(2.0)}")

0:
H 0
RZ(0) 0
H 0

1:
H 0
RZ(-2.0) 0
H 0

2:
H 0
RZ(-4.0) 0
H 0

To take it one step further, you can use Parametric Compilation with exponential_map. For instance:

ham = sZ(0) * sZ(1)
prog = Program()
theta = prog.declare('theta', 'REAL')
prog += exponential_map(ham)(theta)

Exercises

Exercise 1: Quantum Dice

Write a quantum program to simulate throwing an 8-sided die. The Python function you should produce is:

def throw_octahedral_die():
    # return the result of throwing an 8 sided die, an int between 1 and 8, by running a quantum program

Next, extend the program to work for any kind of fair die:

def throw_polyhedral_die(num_sides):
    # return the result of throwing a num_sides sided die by running a quantum program

Exercise 2: Controlled Gates

We can use the full generality of NumPy to construct new gate matrices.

  1. Write a function controlled which takes a \(2\times 2\) matrix \(U\) representing a single qubit operator, and makes a \(4\times 4\) matrix which is a controlled variant of \(U\), with the first argument being the control qubit.
  2. Write a Quil program to define a controlled-\(Y\) gate in this manner. Find the wavefunction when applying this gate to qubit 1 controlled by qubit 0.

Exercise 3: Grover’s Algorithm

Write a quantum program for the single-shot Grover’s algorithm. The Python function you should produce is:

# data is an array of 0's and 1's such that there are exactly three times as many
# 0's as 1's
def single_shot_grovers(data):
    # return an index that contains the value 1

As an example: single_shot_grovers([0,0,1,0]) should return 2.

HINT - Remember that the Grover’s diffusion operator is:

\[\begin{split}\begin{pmatrix} 2/N - 1 & 2/N & \cdots & 2/N \\ 2/N & & &\\ \vdots & & \ddots & \\ 2/N & & & 2/N-1 \end{pmatrix}\end{split}\]

Exercise 4: Prisoner’s Dilemma

A classic strategy game is the prisoner’s dilemma where two prisoners get the minimal penalty if they collaborate and stay silent, get zero penalty if one of them defects and the other collaborates (incurring maximum penalty) and get intermediate penalty if they both defect. This game has an equilibrium where both defect and incur intermediate penalty.

However, things change dramatically when we allow for quantum strategies leading to the Quantum Prisoner’s Dilemma.

Can you design a program that simulates this game?

Exercise 5: Quantum Fourier Transform

The quantum Fourier transform (QFT) is a quantum implementation of the discrete Fourier transform. The Fourier transform can be used to transform a function from the time domain into the frequency domain.

Compute the discrete Fourier transform of [0, 1, 0, 0, 0, 0, 0, 0], using pyQuil:
  1. Write a state preparation quantum program.
  2. Write a function to make a 3-qubit QFT program, taking qubit indices as arguments.
  3. Combine your solutions to part a and b into one program and use the WavefunctionSimulator to get the solution.

Note

For a more challenging initial state, try 01100100.

Solution
Part a: Prepare the initial state

We are going to apply the QFT on the amplitudes of the states.

We want to prepare a state that corresponds to the sequence for which we want to compute the discrete Fourier transform. As the exercise hinted in part b, we need 3 qubits to transform an 8 bit sequence. It is simplest to understand if we think of the qubits as three digits in a binary string (aka bitstring). There are 8 possible values the bitstring can have, and in our quantum state, each of these possibilities has an amplitude. Our 8 indices in the QFT sequence label each of these states. For clarity:

\(|000\rangle\) => 10000000

\(|001\rangle\) => 01000000

…

\(|111\rangle\) -> 00000001

The sequence we want to compute is 01000000, so our initial state is simply \(|001\rangle\). For a bitstring with more than one 1, we would want an equal superposition over all the selected states. (E.g. 01100000 would be an equal superposition of \(|001\rangle\) and \(|010\rangle\)).

To set up the \(|001\rangle\) state, we only have to apply one \(X\)-gate to the zeroth qubit.

from pyquil import Program
from pyquil.gates import *

state_prep = Program(X(0))

We can verify that this works by computing its wavefunction with the Wavefunction Simulator. However, we need to add some “dummy” qubits, because otherwise wavefunction would return a two-element vector for only qubit 0.

from pyquil.api import WavefunctionSimulator

add_dummy_qubits = Program(I(1), I(2))  # The identity gate I has no affect

wf_sim = WavefunctionSimulator()
wavefunction = wf_sim.wavefunction(state_prep + add_dummy_qubits)
print(wavefunction)

(1+0j)|001>

We’ll need wf_sim for part c, too.

Part b: Three qubit QFT program

In this part, we define a function, qft3, to make a 3-qubit QFT quantum program. The algorithm is nicely described on this page. It is a mix of Hadamard and CPHASE gates, with a SWAP gate for bit reversal correction.

from math import pi

def qft3(q0, q1, q2):
    p = Program()
    p += [SWAP(q0, q2),
          H(q0),
          CPHASE(-pi / 2.0, q0, q1),
          H(q1),
          CPHASE(-pi / 4.0, q0, q2),
          CPHASE(-pi / 2.0, q1, q2),
          H(q2)]
    return p

There is a very important detail to recognize here: The function qft3 doesn’t compute the QFT, but rather it makes a quantum program to compute the QFT on qubits q0, q1, and q2.

We can see what this program looks like in Quil notation with print(qft(0, 1, 2)).

SWAP 0 2
H 0
CPHASE(-pi/2) 0 1
H 1
CPHASE(-pi/4) 0 2
CPHASE(-pi/2) 1 2
H 2
Part c: Execute the QFT

Combining parts a and b:

compute_qft_prog = state_prep + qft3(0, 1, 2)
wavefunction = wf_sim.wavefunction(compute_qft_prog)
print(wavefunction.amplitudes)

array([ 3.53553391e-01+0.j        ,  2.50000000e-01-0.25j      ,
        2.16489014e-17-0.35355339j, -2.50000000e-01-0.25j      ,
       -3.53553391e-01+0.j        , -2.50000000e-01+0.25j      ,
       -2.16489014e-17+0.35355339j,  2.50000000e-01+0.25j      ])

We can verify this works by computing the inverse FFT on the output with NumPy and seeing that we get back our input (with some floating point error).

from numpy.fft import ifft
ifft(wavefunction.amplitudes, norm="ortho")

array([0.+0.00000000e+00j, 1.+9.38127079e-17j, 0.+0.00000000e+00j,
       0.-1.53080850e-17j, 0.+0.00000000e+00j, 0.-6.31965379e-17j,
       0.+0.00000000e+00j, 0.-1.53080850e-17j])

After ignoring the terms that are on the order of 1e-17, we get [0, 1, 0, 0, 0, 0, 0, 0], which was our input!

Example: The Meyer-Penny Game

To create intuition for quantum algorithms, it is useful (and fun) to play with the abstraction that the software provides.

The Meyer-Penny Game [1] is a simple example we’ll use from quantum game theory. The interested reader may want to read more about quantum game theory in the article Toward a general theory of quantum games [2]. The Meyer-Penny Game goes as follows:

The Starship Enterprise, during one of its deep-space missions, is facing an immediate calamity at the edge of a wormhole, when a powerful alien suddenly appears. The alien, named Q, offers to help Picard, the captain of the Enterprise, under the condition that Picard beats Q in a simple game of heads or tails.

The rules

Picard is to place a penny heads up into an opaque box. Then Picard and Q take turns to flip or not flip the penny without being able to see it; first Q then P then Q again. After this the penny is revealed; Q wins if it shows heads (H), while tails (T) makes Picard the winner.

Picard vs. Q

Picard quickly estimates that his chance of winning is 50% and agrees to play the game. He loses the first round and insists on playing again. To his surprise Q agrees, and they continue playing several rounds more, each of which Picard loses. How is that possible?

What Picard did not anticipate is that Q has access to quantum tools. Instead of flipping the penny, Q puts the penny into a superposition of heads and tails proportional to the quantum state \(|H\rangle+|T\rangle\). Then no matter whether Picard flips the penny or not, it will stay in a superposition (though the relative sign might change). In the third step Q undoes the superposition and always finds the penny to show heads.

Let’s see how this works!

To simulate the game we first construct the corresponding quantum circuit, which takes two qubits: one to simulate Picard’s choice whether or not to flip the penny, and the other to represent the penny. The initial state for all qubits is \(|0\rangle\) (which is mapped to \(|T\rangle\), tails). To simulate Picard’s decision, we assume that he chooses randomly whether or not to flip the coin, in agreement with the optimal strategy for the classic penny-flip game. This random choice can be created by putting one qubit into an equal superposition, e.g. with the Hadamard gate \(H\), and then measure its state. The measurement will show heads or tails with equal probability \(p_h=p_t=0.5\).

To simulate the penny flip game we take the second qubit and put it into its excited state \(|1\rangle\) (which is mapped to \(|H\rangle\), heads) by applying the X (or NOT) gate. Q’s first move is to apply the Hadamard gate H. Picard’s decision about the flip is simulated as a CNOT operation where the control bit is the outcome of the random number generator described above. Finally Q applies a Hadamard gate again, before we measure the outcome. The full circuit is shown in the figure below.

_images/MeyerPennyGame.png
In pyQuil

We first import and initialize the necessary tools [3]

from pyquil import Program
from pyquil.api import WavefunctionSimulator
from pyquil.gates import *

wf_sim = WavefunctionSimulator()
p = Program()

and then wire it all up into the overall measurement circuit; remember that qubit 0 is the penny, and qubit 1 represents Picard’s choice.

p += X(0)
p += H(0)
p += H(1)
p += CNOT(1, 0)
p += H(0)

We use the quantum mechanics principle of deferred measurement to keep all the measurement logic separate from the gates. Our method call to the WavefunctionSimulator will handle measuring for us [4].

Finally, we play the game several times. (Remember to run your qvm server.)

wf_sim.run_and_measure(p, trials=10)

array([[1, 1],
       [1, 1],
       [1, 1],
       [1, 1],
       [1, 1],
       [1, 0],
       [1, 1],
       [1, 1],
       [1, 1],
       [1, 0]])

In each trial, the first number is the outcome of the game, whereas the second number represents Picard’s choice to flip or not flip the penny.

Inspecting the results, we see that no matter what Picard does, Q will always win!

[1]https://link.aps.org/doi/10.1103/PhysRevLett.82.1052
[2]https://arxiv.org/abs/quant-ph/0611234
[3]See more: Programs and Gates
[4]More about measurements and run_and_measure: Measurement

Source Code Documentation

pyquil.device

class pyquil.device.AbstractDevice[source]

Bases: abc.ABC

get_isa(oneq_type='Xhalves', twoq_type='CZ')[source]

Construct an ISA suitable for targeting by compilation.

This will raise an exception if the requested ISA is not supported by the device.

Parameters:
  • oneq_type – The family of one-qubit gates to target
  • twoq_type – The family of two-qubit gates to target
Return type:

ISA
get_specs()[source]

Construct a Specs object required by compilation

Return type:Specs
qubit_topology()[source]

The connectivity of qubits in this device given as a NetworkX graph.

Return type:Graph
qubits()[source]

A sorted list of qubits in the device topology.

class pyquil.device.Device(name, raw)[source]

Bases: pyquil.device.AbstractDevice

A device (quantum chip) that can accept programs.

Only devices that are online will actively be accepting new programs. In addition to the self._raw attribute, two other attributes are optionally constructed from the entries in self._raw – isa and noise_model – which should conform to the dictionary format required by the .from_dict() methods for ISA and NoiseModel, respectively.

Variables:
  • _raw (dict) – Raw JSON response from the server with additional information about the device.
  • isa (ISA) – The instruction set architecture (ISA) for the device.
  • noise_model (NoiseModel) – The noise model for the device.
Parameters:
  • name – name of the device
  • raw – raw JSON response from the server with additional information about this device.
get_isa(oneq_type='Xhalves', twoq_type='CZ')[source]

Construct an ISA suitable for targeting by compilation.

This will raise an exception if the requested ISA is not supported by the device.

Parameters:
  • oneq_type – The family of one-qubit gates to target
  • twoq_type – The family of two-qubit gates to target
Return type:

ISA
get_specs()[source]

Construct a Specs object required by compilation

isa
qubit_topology()[source]

The connectivity of qubits in this device given as a NetworkX graph.

Return type:Graph
qubits()[source]

A sorted list of qubits in the device topology.

class pyquil.device.Edge(targets, type, dead)

Bases: tuple

Create new instance of Edge(targets, type, dead)

dead

Alias for field number 2

targets

Alias for field number 0

type

Alias for field number 1

pyquil.device.EdgeSpecs

alias of pyquil.device._QubitQubitSpecs

class pyquil.device.ISA[source]

Bases: pyquil.device._ISA

Basic Instruction Set Architecture specification.

Variables:
  • qubits (Sequence[Qubit]) – The qubits associated with the ISA.
  • edges (Sequence[Edge]) – The multi-qubit gates.

Create new instance of _ISA(qubits, edges)

static from_dict(d)[source]

Re-create the ISA from a dictionary representation.

Parameters:d (Dict[str,Any]) – The dictionary representation.
Returns:The restored ISA.
Return type:ISA
to_dict()[source]

Create a JSON-serializable representation of the ISA.

The dictionary representation is of the form:

{
    "1Q": {
        "0": {
            "type": "Xhalves"
        },
        "1": {
            "type": "Xhalves",
            "dead": True
        },
        ...
    },
    "2Q": {
        "1-4": {
            "type": "CZ"
        },
        "1-5": {
            "type": "CZ"
        },
        ...
    },
    ...
}
Returns:A dictionary representation of self.
Return type:Dict[str, Any]
class pyquil.device.NxDevice(topology)[source]

Bases: pyquil.device.AbstractDevice

A shim over the AbstractDevice API backed by a NetworkX graph.

A Device holds information about the physical device. Specifically, you might want to know about connectivity, available gates, performance specs, and more. This class implements the AbstractDevice API for devices not available via get_devices(). Instead, the user is responsible for constructing a NetworkX graph which represents a chip topology.

edges()[source]
Return type:List[Tuple[int, int]]
get_isa(oneq_type='Xhalves', twoq_type='CZ')[source]

Construct an ISA suitable for targeting by compilation.

This will raise an exception if the requested ISA is not supported by the device.

Parameters:
  • oneq_type – The family of one-qubit gates to target
  • twoq_type – The family of two-qubit gates to target
get_specs()[source]

Construct a Specs object required by compilation

qubit_topology()[source]

The connectivity of qubits in this device given as a NetworkX graph.

qubits()[source]

A sorted list of qubits in the device topology.

Return type:List[int]
class pyquil.device.Qubit(id, type, dead)

Bases: tuple

Create new instance of Qubit(id, type, dead)

dead

Alias for field number 2

id

Alias for field number 0

type

Alias for field number 1

pyquil.device.QubitSpecs

alias of pyquil.device._QubitSpecs

class pyquil.device.Specs[source]

Bases: pyquil.device._Specs

Basic specifications for the device, such as gate fidelities and coherence times.

Variables:
  • qubits_specs (List[QubitSpecs]) – The specs associated with individual qubits.
  • edges_specs (List[EdgesSpecs]) – The specs associated with edges, or qubit-qubit pairs.

Create new instance of _Specs(qubits_specs, edges_specs)

T1s()[source]

Get a dictionary of T1s (in seconds) from the specs, keyed by qubit index.

Returns:A dictionary of T1s, in seconds.
Return type:Dict[int, float]
T2s()[source]

Get a dictionary of T2s (in seconds) from the specs, keyed by qubit index.

Returns:A dictionary of T2s, in seconds.
Return type:Dict[int, float]
f1QRBs()[source]

Get a dictionary of single-qubit randomized benchmarking fidelities (normalized to unity) from the specs, keyed by qubit index.

Returns:A dictionary of 1QRBs, normalized to unity.
Return type:Dict[int, float]
fActiveResets()[source]

Get a dictionary of single-qubit active reset fidelities (normalized to unity) from the specs, keyed by qubit index.

Returns:A dictionary of reset fidelities, normalized to unity.
fBellStates()[source]

Get a dictionary of two-qubit Bell state fidelities (normalized to unity) from the specs, keyed by targets (qubit-qubit pairs).

Returns:A dictionary of Bell state fidelities, normalized to unity.
Return type:Dict[tuple(int, int), float]
fCPHASEs()[source]

Get a dictionary of CPHASE fidelities (normalized to unity) from the specs, keyed by targets (qubit-qubit pairs).

Returns:A dictionary of CPHASE fidelities, normalized to unity.
Return type:Dict[tuple(int, int), float]
fCZs()[source]

Get a dictionary of CZ fidelities (normalized to unity) from the specs, keyed by targets (qubit-qubit pairs).

Returns:A dictionary of CZ fidelities, normalized to unity.
Return type:Dict[tuple(int, int), float]
fROs()[source]

Get a dictionary of single-qubit readout fidelities (normalized to unity) from the specs, keyed by qubit index.

Returns:A dictionary of RO fidelities, normalized to unity.
Return type:Dict[int, float]
static from_dict(d)[source]

Re-create the Specs from a dictionary representation.

Parameters:Any] d (Dict[str,) – The dictionary representation.
Returns:The restored Specs.
Return type:Specs
to_dict()[source]

Create a JSON-serializable representation of the device Specs.

The dictionary representation is of the form:

{
    '1Q': {
        "0": {
            "f1QRB": 0.99,
            "T1": 20e-6,
            ...
        },
        "1": {
            "f1QRB": 0.989,
            "T1": 19e-6,
            ...
        },
        ...
    },
    '2Q': {
        "1-4": {
            "fBellState": 0.93,
            "fCZ": 0.92,
            "fCPHASE": 0.91
        },
        "1-5": {
            "fBellState": 0.9,
            "fCZ": 0.89,
            "fCPHASE": 0.88
        },
        ...
    },
    ...
}
Returns:A dctionary representation of self.
Return type:Dict[str, Any]
pyquil.device.THETA = Parameter('theta')

Used as the symbolic parameter in RZ, CPHASE gates.

pyquil.device.gates_in_isa(isa)[source]

Generate the full gateset associated with an ISA.

Parameters:isa (ISA) – The instruction set architecture for a QPU.
Returns:A sequence of Gate objects encapsulating all gates compatible with the ISA.
Return type:Sequence[Gate]
pyquil.device.isa_from_graph(graph, oneq_type='Xhalves', twoq_type='CZ')[source]

Generate an ISA object from a NetworkX graph.

Parameters:
  • graph (Graph) – The graph
  • oneq_type – The type of 1-qubit gate. Currently ‘Xhalves’
  • twoq_type – The type of 2-qubit gate. One of ‘CZ’ or ‘CPHASE’.
Return type:

ISA
pyquil.device.isa_to_graph(isa)[source]

Construct a NetworkX qubit topology from an ISA object.

This discards information about supported gates.

Parameters:isa (ISA) – The ISA.
Return type:Graph
pyquil.device.specs_from_graph(graph)[source]

Generate a Specs object from a NetworkX graph with placeholder values for the actual specs.

Parameters:graph (Graph) – The graph

pyquil.gates

pyquil.gates.I(qubit)[source]

Produces the I instruction.

I = [1, 0]
[0, 1]

This gate is a single qubit identity gate. Note that this gate is different that the NOP instruction as noise channels are typically still applied during the duration of identity gates. Identities will also block parallelization like any other gate.

Parameters:qubit – The qubit apply the gate to.
Returns:A Gate object.
pyquil.gates.X(qubit)[source]

Produces the X instruction.

X = [[0, 1],
[1, 0]]

This gate is a single qubit X-gate.

Parameters:qubit – The qubit apply the gate to.
Returns:A Gate object.
pyquil.gates.Y(qubit)[source]

Produces the Y instruction.

Y = [[0, 0 - 1j],
[0 + 1j, 0]]

This gate is a single qubit Y-gate.

Parameters:qubit – The qubit apply the gate to.
Returns:A Gate object.
pyquil.gates.Z(qubit)[source]

Produces the Z instruction.

Z = [[1, 0],
[0, -1]]

This gate is a single qubit Z-gate.

Parameters:qubit – The qubit apply the gate to.
Returns:A Gate object.
pyquil.gates.H(qubit)[source]
H = (1 / sqrt(2)) * [[1, 1],
[1, -1]]

Produces the H instruction. This gate is a single qubit Hadamard gate.

Parameters:qubit – The qubit apply the gate to.
Returns:A Gate object.
pyquil.gates.S(qubit)[source]

Produces the S instruction.

S = [[1, 0],
[0, 1j]]

This gate is a single qubit S-gate.

Parameters:qubit – The qubit apply the gate to.
Returns:A Gate object.
pyquil.gates.T(qubit)[source]

Produces the T instruction.

T = [[1, 0],
[0, exp(1j * pi / 4)]]

This gate is a single qubit T-gate. It is the same as RZ(pi/4).

Parameters:qubit – The qubit apply the gate to.
Returns:A Gate object.
pyquil.gates.PHASE(angle, qubit)[source]

Produces the PHASE instruction.

PHASE(phi) = [[1, 0],
[0, exp(1j * phi)]]

This is the same as the RZ gate.

Parameters:
  • angle – The angle to rotate around the z-axis on the bloch sphere.
  • qubit – The qubit apply the gate to.
Returns:

A Gate object.
pyquil.gates.RX(angle, qubit)[source]

Produces the RX instruction.

RX(phi) = [[cos(phi / 2), -1j * sin(phi / 2)],
[-1j * sin(phi / 2), cos(phi / 2)]]

This gate is a single qubit X-rotation.

Parameters:
  • angle – The angle to rotate around the x-axis on the bloch sphere.
  • qubit – The qubit apply the gate to.
Returns:

A Gate object.
pyquil.gates.RY(angle, qubit)[source]

Produces the RY instruction.

RY(phi) = [[cos(phi / 2), -sin(phi / 2)],
[sin(phi / 2), cos(phi / 2)]]

This gate is a single qubit Y-rotation.

Parameters:
  • angle – The angle to rotate around the y-axis on the bloch sphere.
  • qubit – The qubit apply the gate to.
Returns:

A Gate object.
pyquil.gates.RZ(angle, qubit)[source]

Produces the RZ instruction.

RZ(phi) = [[cos(phi / 2) - 1j * sin(phi / 2), 0]
[0, cos(phi / 2) + 1j * sin(phi / 2)]]

This gate is a single qubit Z-rotation.

Parameters:
  • angle – The angle to rotate around the z-axis on the bloch sphere.
  • qubit – The qubit apply the gate to.
Returns:

A Gate object.
pyquil.gates.CZ(control, target)[source]

Produces a CZ instruction.

CZ = [[1, 0, 0, 0],
[0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, -1]]

This gate applies to two qubit arguments to produce the controlled-Z gate instruction.

Parameters:
  • control – The control qubit.
  • target – The target qubit. The target qubit has an Z-gate applied to it if the control qubit is in the excited state.
Returns:

A Gate object.
pyquil.gates.CNOT(control, target)[source]

Produces a CNOT instruction.

CNOT = [[1, 0, 0, 0],
[0, 1, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]]

This gate applies to two qubit arguments to produce the controlled-not gate instruction.

Parameters:
  • control – The control qubit.
  • target – The target qubit. The target qubit has an X-gate applied to it if the control qubit is in the excited state.
Returns:

A Gate object.
pyquil.gates.CCNOT(control1, control2, target)[source]

Produces a CCNOT instruction.

CCNOT = [[1, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0, 0, 0], [0, 0, 1, 0, 0, 0, 0, 0], [0, 0, 0, 1, 0, 0, 0, 0], [0, 0, 0, 0, 1, 0, 0, 0], [0, 0, 0, 0, 0, 1, 0, 0], [0, 0, 0, 0, 0, 0, 0, 1], [0, 0, 0, 0, 0, 0, 1, 0]]

This gate applies to three qubit arguments to produce the controlled-controlled-not gate instruction.

Parameters:
  • control1 – The first control qubit.
  • control2 – The second control qubit.
  • target – The target qubit. The target qubit has an X-gate applied to it if both control qubits are in the excited state.
Returns:

A Gate object.
pyquil.gates.CPHASE00(angle, control, target)[source]

Produces a CPHASE00 instruction.

CPHASE00(phi) = diag([exp(1j * phi), 1, 1, 1])

This gate applies to two qubit arguments to produce the variant of the controlled phase instruction that affects the state 00.

Parameters:
  • angle – The input phase angle to apply when both qubits are in the ground state.
  • control – Qubit 1.
  • target – Qubit 2.
Returns:

A Gate object.
pyquil.gates.CPHASE01(angle, control, target)[source]

Produces a CPHASE01 instruction.

CPHASE01(phi) = diag([1.0, exp(1j * phi), 1.0, 1.0])

This gate applies to two qubit arguments to produce the variant of the controlled phase instruction that affects the state 01.

Parameters:
  • angle – The input phase angle to apply when q1 is in the excited state and q2 is in the ground state.
  • control – Qubit 1.
  • target – Qubit 2.
Returns:

A Gate object.
pyquil.gates.CPHASE10(angle, control, target)[source]

Produces a CPHASE10 instruction.

CPHASE10(phi) = diag([1, 1, exp(1j * phi), 1])

This gate applies to two qubit arguments to produce the variant of the controlled phase instruction that affects the state 10.

Parameters:
  • angle – The input phase angle to apply when q2 is in the excited state and q1 is in the ground state.
  • control – Qubit 1.
  • target – Qubit 2.
Returns:

A Gate object.
pyquil.gates.CPHASE(angle, control, target)[source]

Produces a CPHASE instruction, which is a synonym for CPHASE11.

CPHASE(phi) = diag([1, 1, 1, exp(1j * phi)])

This gate applies to two qubit arguments to produce the variant of the controlled phase instruction that affects the state 11.

Parameters:
  • angle – The input phase angle to apply when both qubits are in the excited state.
  • control – Qubit 1.
  • target – Qubit 2.
Returns:

A Gate object.
pyquil.gates.SWAP(q1, q2)[source]

Produces a SWAP instruction.

SWAP = [[1, 0, 0, 0],
[0, 0, 1, 0], [0, 1, 0, 0], [0, 0, 0, 1]]

This gate swaps the state of two qubits.

Parameters:
  • q1 – Qubit 1.
  • q2 – Qubit 2.
Returns:

A Gate object.
pyquil.gates.CSWAP(control, target_1, target_2)[source]
CSWAP = [[1, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0, 0, 0], [0, 0, 1, 0, 0, 0, 0, 0], [0, 0, 0, 1, 0, 0, 0, 0], [0, 0, 0, 0, 1, 0, 0, 0], [0, 0, 0, 0, 0, 0, 1, 0], [0, 0, 0, 0, 0, 1, 0, 0], [0, 0, 0, 0, 0, 0, 0, 1]]

Produces a CSWAP instruction. This gate swaps the state of two qubits.

Parameters:
  • control – The control qubit.
  • target-1 – The first target qubit.
  • target-2 – The second target qubit. The two target states are swapped if the control is in the excited state.
pyquil.gates.ISWAP(q1, q2)[source]

Produces an ISWAP instruction.

ISWAP = [[1, 0, 0, 0],
[0, 0, 1j, 0], [0, 1j, 0, 0], [0, 0, 0, 1]]

This gate swaps the state of two qubits, applying a -i phase to q1 when it is in the excited state and a -i phase to q2 when it is in the ground state.

Parameters:
  • q1 – Qubit 1.
  • q2 – Qubit 2.
Returns:

A Gate object.
pyquil.gates.PSWAP(angle, q1, q2)[source]

Produces a PSWAP instruction.

PSWAP(phi) = [[1, 0, 0, 0],
[0, 0, exp(1j * phi), 0], [0, exp(1j * phi), 0, 0], [0, 0, 0, 1]]

This is a parameterized swap gate.

Parameters:
  • angle – The angle of the phase to apply to the swapped states. This phase is applied to q1 when it is in the excited state and to q2 when it is in the ground state.
  • q1 – Qubit 1.
  • q2 – Qubit 2.
Returns:

A Gate object.
pyquil.gates.WAIT = <pyquil.quilbase.Wait object>

This instruction tells the quantum computation to halt. Typically these is used while classical memory is being manipulated by a CPU in a hybrid classical/quantum algorithm.

Returns:A Wait object.
pyquil.gates.RESET(qubit_index=None)[source]

Reset all qubits or just a specific qubit at qubit_index.

Parameters:qubit_index (Optional[int]) – The address of the qubit to reset. If None, reset all qubits.
Returns:A Reset or ResetQubit Quil AST expression corresponding to a global or targeted reset, respectively.
Return type:Union[Reset, ResetQubit]
pyquil.gates.NOP = <pyquil.quilbase.Nop object>

This instruction applies no operation at that timestep. Typically these are ignored in error-models.

Returns:A Nop object.
pyquil.gates.HALT = <pyquil.quilbase.Halt object>

This instruction ends the program.

Returns:A Halt object.
pyquil.gates.MEASURE(qubit, classical_reg=None)[source]

Produce a MEASURE instruction.

Parameters:
  • qubit – The qubit to measure.
  • classical_reg – The classical register to measure into, or None.
Returns:

A Measurement instance.
pyquil.gates.TRUE(classical_reg)[source]

Produce a TRUE instruction.

Parameters:classical_reg – A classical register to modify.
Returns:An instruction object representing the equivalent MOVE.
pyquil.gates.FALSE(classical_reg)[source]

Produce a FALSE instruction.

Parameters:classical_reg – A classical register to modify.
Returns:An instruction object representing the equivalent MOVE.
pyquil.gates.NOT(classical_reg)[source]

Produce a NOT instruction.

Parameters:classical_reg – A classical register to modify.
Returns:A ClassicalNot instance.
pyquil.gates.AND(classical_reg1, classical_reg2)[source]

Produce an AND instruction.

NOTE: The order of operands was reversed in pyQuil <=1.9 .

Parameters:
  • classical_reg1 – The first classical register, which gets modified.
  • classical_reg2 – The second classical register or immediate value.
Returns:

A ClassicalAnd instance.
pyquil.gates.OR(classical_reg1, classical_reg2)[source]

Produce an OR instruction.

NOTE: Deprecated. Use IOR instead.

Parameters:
  • classical_reg1 – The first classical register.
  • classical_reg2 – The second classical register, which gets modified.
Returns:

A ClassicalOr instance.
pyquil.gates.MOVE(classical_reg1, classical_reg2)[source]

Produce a MOVE instruction.

Parameters:
  • classical_reg1 – The first classical register, which gets modified.
  • classical_reg2 – The second classical register or immediate value.
Returns:

A ClassicalMove instance.
pyquil.gates.EXCHANGE(classical_reg1, classical_reg2)[source]

Produce an EXCHANGE instruction.

Parameters:
  • classical_reg1 – The first classical register, which gets modified.
  • classical_reg2 – The second classical register, which gets modified.
Returns:

A ClassicalExchange instance.
pyquil.gates.IOR(classical_reg1, classical_reg2)[source]

Produce an inclusive OR instruction.

Parameters:
  • classical_reg1 – The first classical register, which gets modified.
  • classical_reg2 – The second classical register or immediate value.
Returns:

A ClassicalOr instance.
pyquil.gates.XOR(classical_reg1, classical_reg2)[source]

Produce an exclusive OR instruction.

Parameters:
  • classical_reg1 – The first classical register, which gets modified.
  • classical_reg2 – The second classical register or immediate value.
Returns:

A ClassicalOr instance.
pyquil.gates.NEG(classical_reg)[source]

Produce a NEG instruction.

Parameters:classical_reg – A classical memory address to modify.
Returns:A ClassicalNeg instance.
pyquil.gates.ADD(classical_reg, right)[source]

Produce an ADD instruction.

Parameters:
  • classical_reg – Left operand for the arithmetic operation. Also serves as the store target.
  • right – Right operand for the arithmetic operation.
Returns:

A ClassicalAdd instance.
pyquil.gates.SUB(classical_reg, right)[source]

Produce a SUB instruction.

Parameters:
  • classical_reg – Left operand for the arithmetic operation. Also serves as the store target.
  • right – Right operand for the arithmetic operation.
Returns:

A ClassicalSub instance.
pyquil.gates.MUL(classical_reg, right)[source]

Produce a MUL instruction.

Parameters:
  • classical_reg – Left operand for the arithmetic operation. Also serves as the store target.
  • right – Right operand for the arithmetic operation.
Returns:

A ClassicalMul instance.
pyquil.gates.DIV(classical_reg, right)[source]

Produce an DIV instruction.

Parameters:
  • classical_reg – Left operand for the arithmetic operation. Also serves as the store target.
  • right – Right operand for the arithmetic operation.
Returns:

A ClassicalDiv instance.
pyquil.gates.EQ(classical_reg1, classical_reg2, classical_reg3)[source]

Produce an EQ instruction.

Parameters:
  • classical_reg1 – Memory address to which to store the comparison result.
  • classical_reg2 – Left comparison operand.
  • classical_reg3 – Right comparison operand.
Returns:

A ClassicalEqual instance.
pyquil.gates.GT(classical_reg1, classical_reg2, classical_reg3)[source]

Produce an GT instruction.

Parameters:
  • classical_reg1 – Memory address to which to store the comparison result.
  • classical_reg2 – Left comparison operand.
  • classical_reg3 – Right comparison operand.
Returns:

A ClassicalGreaterThan instance.
pyquil.gates.GE(classical_reg1, classical_reg2, classical_reg3)[source]

Produce an GE instruction.

Parameters:
  • classical_reg1 – Memory address to which to store the comparison result.
  • classical_reg2 – Left comparison operand.
  • classical_reg3 – Right comparison operand.
Returns:

A ClassicalGreaterEqual instance.
pyquil.gates.LE(classical_reg1, classical_reg2, classical_reg3)[source]

Produce an LE instruction.

Parameters:
  • classical_reg1 – Memory address to which to store the comparison result.
  • classical_reg2 – Left comparison operand.
  • classical_reg3 – Right comparison operand.
Returns:

A ClassicalLessEqual instance.
pyquil.gates.LT(classical_reg1, classical_reg2, classical_reg3)[source]

Produce an LT instruction.

Parameters:
  • classical_reg1 – Memory address to which to store the comparison result.
  • classical_reg2 – Left comparison operand.
  • classical_reg3 – Right comparison operand.
Returns:

A ClassicalLessThan instance.
pyquil.gates.LOAD(target_reg, region_name, offset_reg)[source]

Produce a LOAD instruction.

Parameters:
  • target_reg – LOAD storage target.
  • region_name – Named region of memory to load from.
  • offset_reg – Offset into region of memory to load from. Must be a MemoryReference.
Returns:

A ClassicalLoad instance.
pyquil.gates.STORE(region_name, offset_reg, source)[source]

Produce a STORE instruction.

Parameters:
  • region_name – Named region of memory to store to.
  • offset_reg – Offset into memory region. Must be a MemoryReference.
  • source – Source data. Can be either a MemoryReference or a constant.
Returns:

A ClassicalStore instance.
pyquil.gates.CONVERT(classical_reg1, classical_reg2)[source]

Produce a CONVERT instruction.

Parameters:
  • classical_reg1 – MemoryReference to store to.
  • classical_reg2 – MemoryReference to read from.
Returns:

A ClassicalCONVERT instance.
class pyquil.gates.Gate(name, params, qubits)[source]

Bases: pyquil.quilbase.AbstractInstruction

This is the pyQuil object for a quantum gate instruction.

get_qubits(indices=True)[source]
out()[source]

pyquil.noise

Module for creating and verifying noisy gate and readout definitions.

pyquil.noise.INFINITY = inf

Used for infinite coherence times.

class pyquil.noise.KrausModel[source]

Bases: pyquil.noise._KrausModel

Encapsulate a single gate’s noise model.

Variables:
  • gate (str) – The name of the gate.
  • params (Sequence[float]) – Optional parameters for the gate.
  • targets (Sequence[int]) – The target qubit ids.
  • kraus_ops (Sequence[np.array]) – The Kraus operators (must be square complex numpy arrays).
  • fidelity (float) – The average gate fidelity associated with the Kraus map relative to the ideal operation.

Create new instance of _KrausModel(gate, params, targets, kraus_ops, fidelity)

static from_dict(d)[source]

Recreate a KrausModel from the dictionary representation.

Parameters:d (dict) – The dictionary representing the KrausModel. See to_dict for an example.
Returns:The deserialized KrausModel.
Return type:KrausModel
to_dict()[source]

Create a dictionary representation of a KrausModel.

For example:

{
    "gate": "RX",
    "params": np.pi,
    "targets": [0],
    "kraus_ops": [            # In this example single Kraus op = ideal RX(pi) gate
        [[[0,   0],           # element-wise real part of matrix
          [0,   0]],
          [[0, -1],           # element-wise imaginary part of matrix
          [-1, 0]]]
    ],
    "fidelity": 1.0
}
Returns:A JSON compatible dictionary representation.
Return type:Dict[str,Any]
static unpack_kraus_matrix(m)[source]

Helper to optionally unpack a JSON compatible representation of a complex Kraus matrix.

Parameters:m (Union[list,np.array]) – The representation of a Kraus operator. Either a complex square matrix (as numpy array or nested lists) or a JSON-able pair of real matrices (as nested lists) representing the element-wise real and imaginary part of m.
Returns:A complex square numpy array representing the Kraus operator.
Return type:np.array
class pyquil.noise.NoiseModel[source]

Bases: pyquil.noise._NoiseModel

Encapsulate the QPU noise model containing information about the noisy gates.

Variables:
  • gates (Sequence[KrausModel]) – The tomographic estimates of all gates.
  • assignment_probs (Dict[int,np.array]) – The single qubit readout assignment probability matrices keyed by qubit id.

Create new instance of _NoiseModel(gates, assignment_probs)

static from_dict(d)[source]

Re-create the noise model from a dictionary representation.

Parameters:d (Dict[str,Any]) – The dictionary representation.
Returns:The restored noise model.
Return type:NoiseModel
gates_by_name(name)[source]

Return all defined noisy gates of a particular gate name.

Parameters:name (str) – The gate name.
Returns:A list of noise models representing that gate.
Return type:Sequence[KrausModel]
to_dict()[source]

Create a JSON serializable representation of the noise model.

For example:

{
    "gates": [
        # list of embedded dictionary representations of KrausModels here [...]
    ]
    "assignment_probs": {
        "0": [[.8, .1],
              [.2, .9]],
        "1": [[.9, .4],
              [.1, .6]],
    }
}
Returns:A dictionary representation of self.
Return type:Dict[str,Any]
exception pyquil.noise.NoisyGateUndefined[source]

Bases: Exception

Raise when user attempts to use noisy gate outside of currently supported set.

pyquil.noise.add_decoherence_noise(prog, T1=3e-05, T2=3e-05, gate_time_1q=5e-08, gate_time_2q=1.5e-07, ro_fidelity=0.95)[source]

Add generic damping and dephasing noise to a program.

This high-level function is provided as a convenience to investigate the effects of a generic noise model on a program. For more fine-grained control, please investigate the other methods available in the pyquil.noise module.

In an attempt to closely model the QPU, noisy versions of RX(+-pi/2) and CZ are provided; I and parametric RZ are noiseless, and other gates are not allowed. To use this function, you need to compile your program to this native gate set.

The default noise parameters

  • T1 = 30 us
  • T2 = 30 us
  • 1q gate time = 50 ns
  • 2q gate time = 150 ns

are currently typical for near-term devices.

This function will define new gates and add Kraus noise to these gates. It will translate the input program to use the noisy version of the gates.

Parameters:
  • prog – A pyquil program consisting of I, RZ, CZ, and RX(+-pi/2) instructions
  • T1 (Union[Dict[int,float],float]) – The T1 amplitude damping time either globally or in a dictionary indexed by qubit id. By default, this is 30 us.
  • T2 (Union[Dict[int,float],float]) – The T2 dephasing time either globally or in a dictionary indexed by qubit id. By default, this is also 30 us.
  • gate_time_1q (float) – The duration of the one-qubit gates, namely RX(+pi/2) and RX(-pi/2). By default, this is 50 ns.
  • gate_time_2q (float) – The duration of the two-qubit gates, namely CZ. By default, this is 150 ns.
  • ro_fidelity (Union[Dict[int,float],float]) – The readout assignment fidelity \(F = (p(0|0) + p(1|1))/2\) either globally or in a dictionary indexed by qubit id.
Returns:

A new program with noisy operators.
pyquil.noise.append_kraus_to_gate(kraus_ops, gate_matrix)[source]

Follow a gate gate_matrix by a Kraus map described by kraus_ops.

Parameters:
  • kraus_ops (list) – The Kraus operators.
  • gate_matrix (numpy.ndarray) – The unitary gate.
Returns:

A list of transformed Kraus operators.
pyquil.noise.apply_noise_model(prog, noise_model)[source]

Apply a noise model to a program and generated a ‘noisy-fied’ version of the program.

Parameters:
  • prog (Program) – A Quil Program object.
  • noise_model (NoiseModel) – A NoiseModel, either generated from an ISA or from a simple decoherence model.
Returns:

A new program translated to a noisy gateset and with noisy readout as described by the noisemodel.
Return type:

Program
pyquil.noise.bitstring_probs_to_z_moments(p)[source]

Convert between bitstring probabilities and joint Z moment expectations.

Parameters:p (np.array) – An array that enumerates bitstring probabilities. When flattened out p = [p_00...0, p_00...1, ...,p_11...1]. The total number of elements must therefore be a power of 2. The canonical shape has a separate axis for each qubit, such that p[i,j,...,k] gives the estimated probability of bitstring ij...k.
Returns:z_moments, an np.array with one length-2 axis per qubit which contains the expectations of all monomials in {I, Z_0, Z_1, ..., Z_{n-1}}. The expectations of each monomial can be accessed via:
<Z_0^j_0 Z_1^j_1 ... Z_m^j_m> = z_moments[j_0,j_1,...,j_m]
Return type:np.array
pyquil.noise.combine_kraus_maps(k1, k2)[source]

Generate the Kraus map corresponding to the composition of two maps on the same qubits with k1 being applied to the state after k2.

Parameters:
  • k1 (list) – The list of Kraus operators that are applied second.
  • k2 (list) – The list of Kraus operators that are applied first.
Returns:

A combinatorially generated list of composed Kraus operators.
pyquil.noise.correct_bitstring_probs(p, assignment_probabilities)[source]

Given a 2d array of corrupted bitstring probabilities (outer axis iterates over shots, inner axis over bits) and a list of assignment probability matrices (one for each bit in the readout) compute the corrected probabilities.

Parameters:
  • p (np.array) – An array that enumerates bitstring probabilities. When flattened out p = [p_00...0, p_00...1, ...,p_11...1]. The total number of elements must therefore be a power of 2. The canonical shape has a separate axis for each qubit, such that p[i,j,...,k] gives the estimated probability of bitstring ij...k.
  • assignment_probabilities (List[np.array]) –

    A list of assignment probability matrices per qubit. Each assignment probability matrix is expected to be of the form:

    [[p00 p01]
 [p10 p11]]
Returns:

p_corrected an array with as many dimensions as there are qubits that contains the noisy-readout-corrected estimated probabilities for each measured bitstring, i.e., p[i,j,...,k] gives the estimated probability of bitstring ij...k.
Return type:

np.array
pyquil.noise.corrupt_bitstring_probs(p, assignment_probabilities)[source]

Given a 2d array of true bitstring probabilities (outer axis iterates over shots, inner axis over bits) and a list of assignment probability matrices (one for each bit in the readout, ordered like the inner axis of results) compute the corrupted probabilities.

Parameters:
  • p (np.array) – An array that enumerates bitstring probabilities. When flattened out p = [p_00...0, p_00...1, ...,p_11...1]. The total number of elements must therefore be a power of 2. The canonical shape has a separate axis for each qubit, such that p[i,j,...,k] gives the estimated probability of bitstring ij...k.
  • assignment_probabilities (List[np.array]) –

    A list of assignment probability matrices per qubit. Each assignment probability matrix is expected to be of the form:

    [[p00 p01]
 [p10 p11]]
Returns:

p_corrected an array with as many dimensions as there are qubits that contains the noisy-readout-corrected estimated probabilities for each measured bitstring, i.e., p[i,j,...,k] gives the estimated probability of bitstring ij...k.
Return type:

np.array
pyquil.noise.damping_after_dephasing(T1, T2, gate_time)[source]

Generate the Kraus map corresponding to the composition of a dephasing channel followed by an amplitude damping channel.

Parameters:
  • T1 (float) – The amplitude damping time
  • T2 (float) – The dephasing time
  • gate_time (float) – The gate duration.
Returns:

A list of Kraus operators.
pyquil.noise.damping_kraus_map(p=0.1)[source]

Generate the Kraus operators corresponding to an amplitude damping noise channel.

Parameters:p (float) – The one-step damping probability.
Returns:A list [k1, k2] of the Kraus operators that parametrize the map.
Return type:list
pyquil.noise.decoherence_noise_with_asymmetric_ro(gates, p00=0.975, p11=0.911)[source]

Similar to :py:func`_decoherence_noise_model`, but with asymmetric readout.

For simplicity, we use the default values for T1, T2, gate times, et al. and only allow the specification of readout fidelities.

pyquil.noise.dephasing_kraus_map(p=0.1)[source]

Generate the Kraus operators corresponding to a dephasing channel.

Params float p:The one-step dephasing probability.
Returns:A list [k1, k2] of the Kraus operators that parametrize the map.
Return type:list
pyquil.noise.estimate_assignment_probs(q, trials, cxn, p0=None)[source]

Estimate the readout assignment probabilities for a given qubit q. The returned matrix is of the form:

[[p00 p01]
 [p10 p11]]
Parameters:
  • q (int) – The index of the qubit.
  • trials (int) – The number of samples for each state preparation.
  • cxn (Union[QVMConnection,QPUConnection]) – The quantum abstract machine to sample from.
  • p0 (Program) – A header program to prepend to the state preparation programs.
Returns:

The assignment probability matrix
Return type:

np.array
pyquil.noise.estimate_bitstring_probs(results)[source]

Given an array of single shot results estimate the probability distribution over all bitstrings.

Parameters:results (np.array) – A 2d array where the outer axis iterates over shots and the inner axis over bits.
Returns:An array with as many axes as there are qubit and normalized such that it sums to one. p[i,j,...,k] gives the estimated probability of bitstring ij...k.
Return type:np.array
pyquil.noise.get_noisy_gate(gate_name, params)[source]

Look up the numerical gate representation and a proposed ‘noisy’ name.

Parameters:
  • gate_name (str) – The Quil gate name
  • params (Tuple[float]) – The gate parameters.
Returns:

A tuple (matrix, noisy_name) with the representation of the ideal gate matrix and a proposed name for the noisy version.
Return type:

Tuple[np.array, str]
pyquil.noise.pauli_kraus_map(probabilities)[source]

Generate the Kraus operators corresponding to a pauli channel.

Params list|floats probabilities:
 The 4^num_qubits list of probabilities specifying the desired pauli channel.

There should be either 4 or 16 probabilities specified in the order I, X, Y, Z for 1 qubit or II, IX, IY, IZ, XI, XX, XY, etc for 2 qubits.

For example:

The d-dimensional depolarizing channel \Delta parameterized as
\Delta(\rho) = p \rho + [(1-p)/d] I
is specified by the list of probabilities
[p + (1-p)/d, (1-p)/d,  (1-p)/d), ... , (1-p)/d)]
Returns:A list of the 4^num_qubits Kraus operators that parametrize the map.
Return type:list
pyquil.noise.tensor_kraus_maps(k1, k2)[source]

Generate the Kraus map corresponding to the composition of two maps on different qubits.

Parameters:
  • k1 (list) – The Kraus operators for the first qubit.
  • k2 (list) – The Kraus operators for the second qubit.
Returns:

A list of tensored Kraus operators.

pyquil.parser

Module for parsing Quil programs from text into PyQuil objects

pyquil.parser.parse(quil)[source]

Parse a raw Quil program and return a corresponding list of PyQuil objects.

Parameters:quil (str) – a single or multiline Quil program
Returns:list of instructions
pyquil.parser.parse_program(quil)[source]

Parse a raw Quil program and return a PyQuil program.

Parameters:quil (str) – a single or multiline Quil program
Returns:PyQuil Program object

pyquil.paulis

Module for working with Pauli algebras.

pyquil.paulis.HASH_PRECISION = 1000000.0

The precision used when hashing terms to check equality. The simplify() method uses np.isclose() for coefficient comparisons to 0 which has its own default precision. We can’t use np.isclose() for hashing terms though.

pyquil.paulis.ID()[source]

The identity Pauli Term.

class pyquil.paulis.PauliSum(terms)[source]

Bases: object

A sum of one or more PauliTerms.

Parameters:terms (Sequence) – A Sequence of PauliTerms.
get_programs()[source]

Get a Pyquil Program corresponding to each term in the PauliSum and a coefficient for each program

Returns:(programs, coefficients)
get_qubits()[source]

The support of all the operators in the PauliSum object.

Returns:A list of all the qubits in the sum of terms.
Return type:list
simplify()[source]

Simplifies the sum of Pauli operators according to Pauli algebra rules.

class pyquil.paulis.PauliTerm(op, index, coefficient=1.0)[source]

Bases: object

A term is a product of Pauli operators operating on different qubits.

Create a new Pauli Term with a Pauli operator at a particular index and a leading coefficient.

Parameters:
  • op (string) – The Pauli operator as a string “X”, “Y”, “Z”, or “I”
  • index (int) – The qubit index that that operator is applied to.
  • coefficient (float) – The coefficient multiplying the operator, e.g. 1.5 * Z_1
copy()[source]

Properly creates a new PauliTerm, with a completely new dictionary of operators

classmethod from_list(terms_list, coefficient=1.0)[source]

Allocates a Pauli Term from a list of operators and indices. This is more efficient than multiplying together individual terms.

Parameters:terms_list (list) – A list of tuples, e.g. [(“X”, 0), (“Y”, 1)]
Returns:PauliTerm
get_qubits()[source]

Gets all the qubits that this PauliTerm operates on.

id(sort_ops=True)[source]

Returns an identifier string for the PauliTerm (ignoring the coefficient).

Don’t use this to compare terms. This function will not work with qubits that aren’t sortable.

Parameters:sort_ops – Whether to sort operations by qubit. This is True by default for backwards compatibility but will change in pyQuil 2.0. Callers should never rely on comparing id’s for testing equality. See operations_as_set instead.
Returns:A string representation of this term’s operations.
Return type:string
operations_as_set()[source]

Return a frozenset of operations in this term.

Use this in place of id() if the order of operations in the term does not matter.

Returns:frozenset of strings representing Pauli operations
pauli_string(qubits=None)[source]

Return a string representation of this PauliTerm mod its phase, as a concatenation of the string representation of the >>> p = PauliTerm(“X”, 0) * PauliTerm(“Y”, 1, 1.j) >>> p.pauli_string() “XY” >>> p.pauli_string([0]) “X” >>> p.pauli_string([0, 2]) “XI”

Parameters:qubits (list) – The list of qubits to represent, given as ints. If None, defaults to all qubits in this PauliTerm.
Returns:The string representation of this PauliTerm, modulo its phase.
Return type:String
program
exception pyquil.paulis.UnequalLengthWarning(*args, **kwargs)[source]

Bases: Warning

pyquil.paulis.ZERO()[source]

The zero Pauli Term.

pyquil.paulis.check_commutation(pauli_list, pauli_two)[source]

Check if commuting a PauliTerm commutes with a list of other terms by natural calculation. Derivation similar to arXiv:1405.5749v2 fo the check_commutation step in the Raesi, Wiebe, Sanders algorithm (arXiv:1108.4318, 2011).

Parameters:
  • pauli_list (list) – A list of PauliTerm objects
  • pauli_two_term (PauliTerm) – A PauliTerm object
Returns:

True if pauli_two object commutes with pauli_list, False otherwise
Return type:

bool
pyquil.paulis.commuting_sets(pauli_terms)[source]

Gather the Pauli terms of pauli_terms variable into commuting sets

Uses algorithm defined in (Raeisi, Wiebe, Sanders, arXiv:1108.4318, 2011) to find commuting sets. Except uses commutation check from arXiv:1405.5749v2

Parameters:pauli_terms (PauliSum) – A PauliSum object
Returns:List of lists where each list contains a commuting set
Return type:list
pyquil.paulis.exponential_map(term)[source]

Creates map alpha -> exp(-1j*alpha*term) represented as a Program.

Parameters:term (PauliTerm) – Tests is a PauliTerm is the identity operator
Returns:Program
Return type:Function
pyquil.paulis.exponentiate(term)[source]

Creates a pyQuil program that simulates the unitary evolution exp(-1j * term)

Parameters:term (PauliTerm) – Tests is a PauliTerm is the identity operator
Returns:A Program object
Return type:Program
pyquil.paulis.exponentiate_commuting_pauli_sum(pauli_sum)[source]

Returns a function that maps all substituent PauliTerms and sums them into a program. NOTE: Use this function with care. Substituent PauliTerms should commute.

Parameters:pauli_sum (PauliSum) – PauliSum to exponentiate.
Returns:A function that parametrizes the exponential.
Return type:function
pyquil.paulis.integer_types = (<class 'int'>, <class 'numpy.int64'>, <class 'numpy.int32'>, <class 'numpy.int16'>, <class 'numpy.int8'>)

Explicitly include numpy integer dtypes (for python 3).

pyquil.paulis.is_identity(term)[source]

Check if Pauli Term is a scalar multiple of identity

Parameters:term (PauliTerm) – A PauliTerm object
Returns:True if the PauliTerm is a scalar multiple of identity, false otherwise
Return type:bool
pyquil.paulis.is_zero(pauli_object)[source]

Tests to see if a PauliTerm or PauliSum is zero.

Parameters:pauli_object – Either a PauliTerm or PauliSum
Returns:True if PauliTerm is zero, False otherwise
Return type:bool
pyquil.paulis.sI(q)[source]

A function that returns the identity operator on a particular qubit.

Parameters:qubit_index (int) – The index of the qubit
Returns:A PauliTerm object
Return type:PauliTerm
pyquil.paulis.sX(q)[source]

A function that returns the sigma_X operator on a particular qubit.

Parameters:qubit_index (int) – The index of the qubit
Returns:A PauliTerm object
Return type:PauliTerm
pyquil.paulis.sY(q)[source]

A function that returns the sigma_Y operator on a particular qubit.

Parameters:qubit_index (int) – The index of the qubit
Returns:A PauliTerm object
Return type:PauliTerm
pyquil.paulis.sZ(q)[source]

A function that returns the sigma_Z operator on a particular qubit.

Parameters:qubit_index (int) – The index of the qubit
Returns:A PauliTerm object
Return type:PauliTerm
pyquil.paulis.simplify_pauli_sum(pauli_sum)[source]
pyquil.paulis.suzuki_trotter(trotter_order, trotter_steps)[source]

Generate trotterization coefficients for a given number of Trotter steps.

U = exp(A + B) is approximated as exp(w1*o1)exp(w2*o2)… This method returns a list [(w1, o1), (w2, o2), … , (wm, om)] of tuples where o=0 corresponds to the A operator, o=1 corresponds to the B operator, and w is the coefficient in the exponential. For example, a second order Suzuki-Trotter approximation to exp(A + B) results in the following [(0.5/trotter_steps, 0), (1/trotteri_steps, 1), (0.5/trotter_steps, 0)] * trotter_steps.

Parameters:
  • trotter_order (int) – order of Suzuki-Trotter approximation
  • trotter_steps (int) – number of steps in the approximation
Returns:

List of tuples corresponding to the coefficient and operator type: o=0 is A and o=1 is B.
Return type:

list
pyquil.paulis.term_with_coeff(term, coeff)[source]

Change the coefficient of a PauliTerm.

Parameters:
  • term (PauliTerm) – A PauliTerm object
  • coeff (Number) – The coefficient to set on the PauliTerm
Returns:

A new PauliTerm that duplicates term but sets coeff
Return type:

PauliTerm
pyquil.paulis.trotterize(first_pauli_term, second_pauli_term, trotter_order=1, trotter_steps=1)[source]

Create a Quil program that approximates exp( (A + B)t) where A and B are PauliTerm operators.

Parameters:
  • first_pauli_term (PauliTerm) – PauliTerm denoted A
  • second_pauli_term (PauliTerm) – PauliTerm denoted B
  • trotter_order (int) – Optional argument indicating the Suzuki-Trotter approximation order–only accepts orders 1, 2, 3, 4.
  • trotter_steps (int) – Optional argument indicating the number of products to decompose the exponential into.
Returns:

Quil program
Return type:

Program

pyquil.quilbase

Contains the core pyQuil objects that correspond to Quil instructions.

class pyquil.quilbase.AbstractInstruction[source]

Bases: object

Abstract class for representing single instructions.

out()[source]
class pyquil.quilbase.ArithmeticBinaryOp(left, right)[source]

Bases: pyquil.quilbase.AbstractInstruction

The abstract class for binary arithmetic classical instructions.

out()[source]
class pyquil.quilbase.ClassicalAdd(left, right)[source]

Bases: pyquil.quilbase.ArithmeticBinaryOp

The ADD instruction.

op = 'ADD'
class pyquil.quilbase.ClassicalAnd(left, right)[source]

Bases: pyquil.quilbase.LogicalBinaryOp

WARNING: The operand order for ClassicalAnd has changed. In pyQuil versions <= 1.9, AND had signature

AND %source %target

Now, AND has signature

AND %target %source
op = 'AND'
class pyquil.quilbase.ClassicalComparison(target, left, right)[source]

Bases: pyquil.quilbase.AbstractInstruction

Abstract class for ternary comparison instructions.

out()[source]
class pyquil.quilbase.ClassicalConvert(left, right)[source]

Bases: pyquil.quilbase.AbstractInstruction

The CONVERT instruction.

op = 'CONVERT'
out()[source]
class pyquil.quilbase.ClassicalDiv(left, right)[source]

Bases: pyquil.quilbase.ArithmeticBinaryOp

The DIV instruction.

op = 'DIV'
class pyquil.quilbase.ClassicalEqual(target, left, right)[source]

Bases: pyquil.quilbase.ClassicalComparison

The EQ comparison instruction.

op = 'EQ'
class pyquil.quilbase.ClassicalExchange(left, right)[source]

Bases: pyquil.quilbase.AbstractInstruction

The EXCHANGE instruction.

op = 'EXCHANGE'
out()[source]
class pyquil.quilbase.ClassicalExclusiveOr(left, right)[source]

Bases: pyquil.quilbase.LogicalBinaryOp

The XOR instruction.

op = 'XOR'
class pyquil.quilbase.ClassicalFalse(target)[source]

Bases: pyquil.quilbase.ClassicalMove

Deprecated class.

class pyquil.quilbase.ClassicalGreaterEqual(target, left, right)[source]

Bases: pyquil.quilbase.ClassicalComparison

The GE comparison instruction.

op = 'GE'
class pyquil.quilbase.ClassicalGreaterThan(target, left, right)[source]

Bases: pyquil.quilbase.ClassicalComparison

The GT comparison instruction.

op = 'GT'
class pyquil.quilbase.ClassicalInclusiveOr(left, right)[source]

Bases: pyquil.quilbase.LogicalBinaryOp

The IOR instruction.

op = 'IOR'
class pyquil.quilbase.ClassicalLessEqual(target, left, right)[source]

Bases: pyquil.quilbase.ClassicalComparison

The LE comparison instruction.

op = 'LE'
class pyquil.quilbase.ClassicalLessThan(target, left, right)[source]

Bases: pyquil.quilbase.ClassicalComparison

The LT comparison instruction.

op = 'LT'
class pyquil.quilbase.ClassicalLoad(target, left, right)[source]

Bases: pyquil.quilbase.AbstractInstruction

The LOAD instruction.

op = 'LOAD'
out()[source]
class pyquil.quilbase.ClassicalMove(left, right)[source]

Bases: pyquil.quilbase.AbstractInstruction

The MOVE instruction.

WARNING: In pyQuil 2.0, the order of operands is as MOVE <target> <source>.
In pyQuil 1.9, the order of operands was MOVE <source> <target>. These have reversed.
op = 'MOVE'
out()[source]
class pyquil.quilbase.ClassicalMul(left, right)[source]

Bases: pyquil.quilbase.ArithmeticBinaryOp

The MUL instruction.

op = 'MUL'
class pyquil.quilbase.ClassicalNeg(target)[source]

Bases: pyquil.quilbase.UnaryClassicalInstruction

The NEG instruction.

op = 'NEG'
class pyquil.quilbase.ClassicalNot(target)[source]

Bases: pyquil.quilbase.UnaryClassicalInstruction

The NOT instruction.

op = 'NOT'
class pyquil.quilbase.ClassicalOr(left, right)[source]

Bases: pyquil.quilbase.ClassicalInclusiveOr

Deprecated class.

class pyquil.quilbase.ClassicalStore(target, left, right)[source]

Bases: pyquil.quilbase.AbstractInstruction

The STORE instruction.

op = 'STORE'
out()[source]
class pyquil.quilbase.ClassicalSub(left, right)[source]

Bases: pyquil.quilbase.ArithmeticBinaryOp

The SUB instruction.

op = 'SUB'
class pyquil.quilbase.ClassicalTrue(target)[source]

Bases: pyquil.quilbase.ClassicalMove

Deprecated class.

class pyquil.quilbase.Declare(name, memory_type, memory_size=1, shared_region=None, offsets=None)[source]

Bases: pyquil.quilbase.AbstractInstruction

A DECLARE directive.

This is printed in Quil as:

DECLARE <name> <memory-type> (SHARING <other-name> (OFFSET <amount> <type>)* )?
asdict()[source]
out()[source]
class pyquil.quilbase.DefGate(name, matrix, parameters=None)[source]

Bases: pyquil.quilbase.AbstractInstruction

A DEFGATE directive.

Parameters:
  • name (string) – The name of the newly defined gate.
  • matrix (array-like) – {list, nparray, np.matrix} The matrix defining this gate.
  • parameters (list) – list of parameters that are used in this gate
get_constructor()[source]
Returns:A function that constructs this gate on variable qubit indices. E.g. mygate.get_constructor()(1) applies the gate to qubit 1.
num_args()[source]
Returns:The number of qubit arguments the gate takes.
Return type:int
out()[source]

Prints a readable Quil string representation of this gate.

Returns:String representation of a gate
Return type:string
class pyquil.quilbase.Gate(name, params, qubits)[source]

Bases: pyquil.quilbase.AbstractInstruction

This is the pyQuil object for a quantum gate instruction.

get_qubits(indices=True)[source]
out()[source]
class pyquil.quilbase.Halt[source]

Bases: pyquil.quilbase.SimpleInstruction

The HALT instruction.

op = 'HALT'
class pyquil.quilbase.Jump(target)[source]

Bases: pyquil.quilbase.AbstractInstruction

Representation of an unconditional jump instruction (JUMP).

out()[source]
class pyquil.quilbase.JumpConditional(target, condition)[source]

Bases: pyquil.quilbase.AbstractInstruction

Abstract representation of an conditional jump instruction.

out()[source]
class pyquil.quilbase.JumpTarget(label)[source]

Bases: pyquil.quilbase.AbstractInstruction

Representation of a target that can be jumped to.

out()[source]
class pyquil.quilbase.JumpUnless(target, condition)[source]

Bases: pyquil.quilbase.JumpConditional

The JUMP-UNLESS instruction.

op = 'JUMP-UNLESS'
class pyquil.quilbase.JumpWhen(target, condition)[source]

Bases: pyquil.quilbase.JumpConditional

The JUMP-WHEN instruction.

op = 'JUMP-WHEN'
class pyquil.quilbase.LogicalBinaryOp(left, right)[source]

Bases: pyquil.quilbase.AbstractInstruction

The abstract class for binary logical classical instructions.

out()[source]
class pyquil.quilbase.Measurement(qubit, classical_reg=None)[source]

Bases: pyquil.quilbase.AbstractInstruction

This is the pyQuil object for a Quil measurement instruction.

get_qubits(indices=True)[source]
out()[source]
class pyquil.quilbase.Nop[source]

Bases: pyquil.quilbase.SimpleInstruction

The NOP instruction.

op = 'NOP'
class pyquil.quilbase.Pragma(command, args=(), freeform_string='')[source]

Bases: pyquil.quilbase.AbstractInstruction

A PRAGMA instruction.

This is printed in QUIL as:

PRAGMA <command> <arg1> <arg2> ... <argn> "<freeform_string>"
out()[source]
class pyquil.quilbase.RawInstr(instr_str)[source]

Bases: pyquil.quilbase.AbstractInstruction

A raw instruction represented as a string.

out()[source]
class pyquil.quilbase.Reset[source]

Bases: pyquil.quilbase.SimpleInstruction

The RESET instruction.

op = 'RESET'
class pyquil.quilbase.ResetQubit(qubit)[source]

Bases: pyquil.quilbase.AbstractInstruction

This is the pyQuil object for a Quil targeted reset instruction.

get_qubits(indices=True)[source]
out()[source]
class pyquil.quilbase.SimpleInstruction[source]

Bases: pyquil.quilbase.AbstractInstruction

Abstract class for simple instructions with no arguments.

out()[source]
class pyquil.quilbase.UnaryClassicalInstruction(target)[source]

Bases: pyquil.quilbase.AbstractInstruction

The abstract class for unary classical instructions.

out()[source]
class pyquil.quilbase.Wait[source]

Bases: pyquil.quilbase.SimpleInstruction

The WAIT instruction.

op = 'WAIT'

pyquil.wavefunction

Module containing the Wavefunction object and methods for working with wavefunctions.

class pyquil.wavefunction.Wavefunction(amplitude_vector)[source]

Bases: object

Encapsulate a wavefunction representing a quantum state as returned by the QVM.

Note

The elements of the wavefunction are ordered by bitstring. E.g., for two qubits the order is 00, 01, 10, 11, where the the bits are ordered in reverse by the qubit index, i.e., for qubits 0 and 1 the bitstring 01 indicates that qubit 0 is in the state 1. See also the related documentation section in the QVM Overview.

Initializes a wavefunction

Parameters:amplitude_vector – A numpy array of complex amplitudes
static from_bit_packed_string(coef_string)[source]

From a bit packed string, unpacks to get the wavefunction :param bytes coef_string: :return:

get_outcome_probs()[source]

Parses a wavefunction (array of complex amplitudes) and returns a dictionary of outcomes and associated probabilities.

Returns:A dict with outcomes as keys and probabilities as values.
Return type:dict
static ground(qubit_num)[source]
plot(qubit_subset=None)[source]

Plots a bar chart with bitstring on the x axis and probability on the y axis.

Parameters:qubit_subset (list) – Optional parameter used for plotting a subset of the Hilbert space.
pretty_print(decimal_digits=2)[source]

Returns a string repr of the wavefunction, ignoring all outcomes with approximately zero amplitude (up to a certain number of decimal digits) and rounding the amplitudes to decimal_digits.

Parameters:decimal_digits (int) – The number of digits to truncate to.
Returns:A dict with outcomes as keys and complex amplitudes as values.
Return type:str
pretty_print_probabilities(decimal_digits=2)[source]

Prints outcome probabilities, ignoring all outcomes with approximately zero probabilities (up to a certain number of decimal digits) and rounding the probabilities to decimal_digits.

Parameters:decimal_digits (int) – The number of digits to truncate to.
Returns:A dict with outcomes as keys and probabilities as values.
Return type:dict
probabilities()[source]

Returns an array of probabilities in lexicographical order

sample_bitstrings(n_samples)[source]

Sample bitstrings from the distribution defined by the wavefunction.

Parameters:n_samples – The number of bitstrings to sample
Returns:An array of shape (n_samples, n_qubits)
static zeros(qubit_num)[source]

Constructs the groundstate wavefunction for a given number of qubits.

Parameters:qubit_num (int) –
Returns:A Wavefunction in the ground state
Return type:Wavefunction
pyquil.wavefunction.get_bitstring_from_index(index, qubit_num)[source]

Returns the bitstring in lexical order that corresponds to the given index in 0 to 2^(qubit_num) :param int index: :param int qubit_num: :return: the bitstring :rtype: str

Changelog

v2.0.0 (November 1, 2018)

PyQuil 2.0 is a major release of pyQuil, Rigetti’s toolkit for constructing and running quantum programs. This release contains many major changes including:

  1. The introduction of Quantum Cloud Services. Access Rigetti’s QPUs from co-located classical compute resources for minimal latency. The web API for running QVM and QPU jobs has been deprecated and cannot be accessed with pyQuil 2.0
  2. Advances in classical control systems and compilation allowing the pre-compilation of parametric binary executables for rapid hybrid algorithm iteration.
  3. Changes to Quil—our quantum instruction language—to provide easier ways of interacting with classical memory.

The new QCS access model and features will allow you to execute hybrid quantum algorithms several orders of magnitude (!) faster than the previous web endpoint. However, to fully exploit these speed increases you must update your programs to use the latest pyQuil features and APIs. Please read Forest 2.0: Migration Guide for a comprehensive migration guide.

An incomplete list of significant changes:

  • Python 2 is no longer supported. Please use Python 3.6+
  • Parametric gates are now normal functions. You can no longer write RX(pi/2)(0) to get a Quil RX(pi/2) 0 instruction. Just use RX(pi/2, 0).
  • Gates support keyword arguments, so you can write RX(angle=pi/2, qubit=0).
  • All async methods have been removed from QVMConnection and QVMConnection is deprecated. QPUConnection has been removed in accordance with the QCS access model. Use pyquil.get_qc() as the primary means of interacting with the QVM or QPU.
  • WavefunctionSimulator allows unfettered access to wavefunction properties and routines. These methods and properties previously lived on QVMConnection and have been deprecated there.
  • Classical memory in Quil must be declared with a name and type. Please read Forest 2.0: Migration Guide for more.
  • Compilation has changed. There are now different Compiler objects that target either the QPU or QVM. You must explicitly compile your programs to run on a QPU or a realistic QVM.

Version 2.0.1 was released on November 9, 2018 and includes documentation changes only. This release is only available as a git tag. We have not pushed a new package to PyPI.

v1.9 (June 6, 2018)

We’re happy to announce the release of pyQuil 1.9. PyQuil is Rigetti’s toolkit for constructing and running quantum programs. This release is the latest in our series of regular releases, and it’s filled with convenience features, enhancements, bug fixes, and documentation improvements.

Special thanks to community members sethuiyer, vtomole, rht, akarazeev, ejdanderson, markf94, playadust, and kadora626 for contributing to this release!

Qubit placeholders

One of the focuses of this release is a re-worked concept of “Qubit Placeholders”. These are logical qubits that can be used to construct programs. Now, a program containing qubit placeholders must be “addressed” prior to running on a QPU or QVM. The addressing stage involves mapping each qubit placeholder to a physical qubit (represented as an integer). For example, if you have a 3 qubit circuit that you want to run on different sections of the Agave chip, you now can prepare one Program and address it to many different subgraphs of the chip topology. Check out the QubitPlaceholder example notebook for more.

To support this idea, we’ve refactored parts of Pyquil to remove the assumption that qubits can be “sorted”. While true for integer qubit labels, this probably isn’t true in general. A notable change can be found in the construction of a PauliSum: now terms will stay in the order they were constructed.

  • PauliTerm now remembers the order of its operations. sX(1)*sZ(2) will compile to different Quil code than sZ(2)*sX(1), although the terms will still be equal according to the __eq__ method. During PauliSum combination of like terms, a warning will be emitted if two terms are combined that have different orders of operation.
  • PauliTerm.id() takes an optional argument sort_ops which defaults to True for backwards compatibility. However, this function should not be used for comparing term-type like it has been used previously. Use PauliTerm.operations_as_set() instead. In the future, sort_ops will default to False and will eventually be removed.
  • Program.alloc() has been deprecated. Please instantiate QubitPlaceholder() directly or request a “register” (list) of n placeholders by using the class constructor QubitPlaceholder.register(n)().
  • Programs must contain either (1) all instantiated qubits with integer indexes or (2) all placeholder qubits of type QubitPlaceholder. We have found that most users use (1) but (2) will become useful with larger and more diverse devices.
  • Programs that contain qubit placeholders must be explicitly addressed prior to execution. Previously, qubits would be assigned “under the hood” to integers 0…N. Now, you must use address_qubits() which returns a new program with all qubits indexed depending on the qubit_mapping argument. The original program is unaffected and can be “readdressed” multiple times.
  • PauliTerm can now accept QubitPlaceholder in addition to integers.
  • QubitPlaceholder is no longer a subclass of Qubit. LabelPlaceholder is no longer a subclass of Label.
  • QuilAtom subclasses’ hash functions have changed.
Randomized benchmarking sequence generation

Pyquil now includes support for performing a simple benchmarking routine - randomized benchmarking. There is a new method in the CompilerConnection that will return sequences of pyquil programs, corresponding to elements of the Clifford group. These programs are uniformly randomly sampled, and have the property that they compose to the identity. When concatenated and run as one program, these programs can be used in a procedure called randomized benchmarking to gain insight about the fidelity of operations on a QPU.

In addition, the CompilerConnection has another new method, apply_clifford_to_pauli() which conjugates PauliTerms by Program that are composed of Clifford gates. That is to say, given a circuit C, that contains only gates corresponding to elements of the Clifford group, and a tensor product of elements P, from the Pauli group, this method will compute $PCP^{dagger}$. Such a procedure can be used in various ways. An example is predicting the effect a Clifford circuit will have on an input state modeled as a density matrix, which can be written as a sum of Pauli matrices.

Ease of Use

This release includes some quality-of-life improvements such as the ability to initialize programs with generator expressions, sensible defaults for Program.measure_all(), and sensible defaults for classical_addresses in run() methods.

  • Program can be initiated with a generator expression.
  • Program.measure_all() (with no arguments) will measure all qubits in a program.
  • classical_addresses is now optional in QVM and QPU run() methods. By default, any classical addresses targeted by MEASURE will be returned.
  • QVMConnection.pauli_expectation() accepts PauliSum as arguments. This offers a more sensible API compared to QVMConnection.expectation().
  • pyQuil will now retry jobs every 10 seconds if the QPU is re-tuning.
  • CompilerConnection.compile() now takes an optional argument isa that allows per-compilation specification of the target ISA.
  • An empty program will trigger an exception if you try to run it.
Supported versions of Python

We strongly support using Python 3 with Pyquil. Although this release works with Python 2, we are dropping official support for this legacy language and moving to community support for Python 2. The next major release of Pyquil will introduce Python 3.5+ only features and will no longer work without modification for Python 2.

Bug fixes
  • shift_quantum_gates has been removed. Users who relied on this functionality should use QubitPlaceholder and address_qubits() to achieve the same result. Users should also double-check data resulting from use of this function as there were several edge cases which would cause the shift to be applied incorrectly resulting in badly-addressed qubits.
  • Slightly perturbed angles when performing RX gates under a Kraus noise model could result in incorrect behavior.
  • The quantum die example returned incorrect values when n = 2^m.

Introduction to Quantum Computing

With every breakthrough in science there is the potential for new technology. For over twenty years, researchers have done inspiring work in quantum mechanics, transforming it from a theory for understanding nature into a fundamentally new way to engineer computing technology. This field, quantum computing, is beautifully interdisciplinary, and impactful in two major ways:

  1. It reorients the relationship between physics and computer science. Physics does not just place restrictions on what computers we can design, it also grants new power and inspiration.
  2. It can simulate nature at its most fundamental level, allowing us to solve deep problems in quantum chemistry, materials discovery, and more.

Quantum computing has come a long way, and in the next few years there will be significant breakthroughs in the field. To get here, however, we have needed to change our intuition for computation in many ways. As with other paradigms — such as object-oriented programming, functional programming, distributed programming, or any of the other marvelous ways of thinking that have been expressed in code over the years — even the basic tenants of quantum computing opens up vast new potential for computation.

However, unlike other paradigms, quantum computing goes further. It requires an extension of classical probability theory. This extension, and the core of quantum computing, can be formulated in terms of linear algebra. Therefore, we begin our investigation into quantum computing with linear algebra and probability.

From Bit to Qubit

Probabilistic Bits as Vector Spaces

From an operational perspective, a bit is described by the results of measurements performed on it. Let the possible results of measuring a bit (0 or 1) be represented by orthonormal basis vectors \(\vec{0}\) and \(\vec{1}\). We will call these vectors outcomes. These outcomes span a two-dimensional vector space that represents a probabilistic bit. A probabilistic bit can be represented as a vector

\[\vec{v} = a\,\vec{0} + b\,\vec{1},\]

where \(a\) represents the probability of the bit being 0 and \(b\) represents the probability of the bit being 1. This clearly also requires that \(a+b=1\). In this picture the system (the probabilistic bit) is a two-dimensional real vector space and a state of a system is a particular vector in that vector space.

import numpy as np
import matplotlib.pyplot as plt

outcome_0 = np.array([1.0, 0.0])
outcome_1 = np.array([0.0, 1.0])
a = 0.75
b = 0.25

prob_bit = a*outcome_0 + b*outcome_1

X,Y = prob_bit
plt.figure()
ax = plt.gca()
ax.quiver(X,Y,angles='xy',scale_units='xy',scale=1)
ax.set_xlim([0,1])
ax.set_ylim([0,1])
plt.draw()
plt.show()
_images/output_4_1.png

Given some state vector, like the one plotted above, we can find the probabilities associated with each outcome by projecting the vector onto the basis outcomes. This gives us the following rule:

\[\begin{split}\operatorname{Pr}(0) = \vec{v}^T.\vec{0} = a \\ \operatorname{Pr}(1) = \vec{v}^T.\vec{1} = b,\end{split}\]

where Pr(0) and Pr(1) are the probabilities of the 0 and 1 outcomes respectively.

Dirac Notation

Physicists have introduced a convenient notation for the vector transposes and dot products we used in the previous example. This notation, called Dirac notation in honor of the great theoretical physicist Paul Dirac, allows us to define

\[\begin{split}\vec{v} = |\,v\rangle \\ \vec{v}^T = \langle v\,| \\ \vec{u}^T.\vec{v} = \langle u\,|\,v \rangle.\end{split}\]

Thus, we can rewrite our “measurement rule” in this notation as

\[\begin{split}Pr(0) = \langle v\,|\,0 \rangle = a \\ Pr(1) = \langle v\,|\,1 \rangle = b.\end{split}\]

We will use this notation throughout the rest of this introduction.

Multiple Probabilistic Bits

This vector space interpretation of a single probabilistic bit can be straightforwardly extended to multiple bits. Let us take two coins as an example (labelled 0 and 1 instead of H and T since we are programmers). Their states can be represented as

\[\begin{split} |\,u\rangle = \frac{1}{2}|\,0_u\rangle + \frac{1}{2}|\,1_u\rangle \\ |\,v\rangle = \frac{1}{2}|\,0_v\rangle + \frac{1}{2}|\,1_v\rangle,\end{split}\]

where \(1_u\) represents the 1 outcome on coin \(u\). The combined system of the two coins has four possible outcomes \(\{ 0_u0_v,\;0_u1_v,\;1_u0_v,\;1_u1_v \}\) that are the basis states of a larger four-dimensional vector space. The rule for constructing a combined state is to take the tensor product of individual states, e.g.

\[|\,u\rangle\otimes|\,v\rangle = \frac{1}{4}|\,0_u0_v\rangle+\frac{1}{4}|\,0_u1_v\rangle+\frac{1}{4}|\,1_u0_v\rangle+\frac{1}{4}|\,1_u1_v\rangle.\]

Then, the combined space is simply the space spanned by the tensor products of all pairs of basis vectors of the two smaller spaces.

We will talk more about these larger spaces in the quantum case, but it is important to note that not all composite states can be written as tensor products of sub-states. (Consider the state \(\frac{1}{2}|\,0_u0_v\rangle + \frac{1}{2}|\,1_u1_v\rangle\).) In general, the combined state for \(n\) probabilistic bits is a vector of size \(2^n\) and is given by \(\bigotimes_{i=0}^{n-1}|\,v_i\rangle\).

Qubits

Quantum mechanics rewrites these rules to some extent. A quantum bit, called a qubit, is the quantum analog of a bit in that it has two outcomes when it is measured. Similar to the previous section, a qubit can also be represented in a vector space, but with complex coefficients instead of real ones. A qubit system is a two-dimensional complex vector space, and the state of a qubit is a complex vector in that space. Again we will define a basis of outcomes \(\{|\,0\rangle, |\,1\rangle\}\) and let a generic qubit state be written as

\[\alpha |\,0\rangle + \beta |\,1\rangle.\]

Since these coefficients can be imaginary, they cannot be simply interpreted as probabilities of their associated outcomes. Instead we rewrite the rule for outcomes in the following manner:

\[\begin{split}\operatorname{Pr}(0) = |\langle v\,|\,0 \rangle|^2 = |\alpha|^2 \\ \operatorname{Pr}(1) = |\langle v\,|\,1 \rangle|^2 = |\beta|^2,\end{split}\]

and as long as \(|\alpha|^2 + |\beta|^2 = 1\) we are able to recover acceptable probabilities for outcomes based on our new complex vector.

This switch to complex vectors means that rather than representing a state vector in a plane, we instead to represent the vector on a sphere (called the Bloch sphere in quantum mechanics literature). From this perspective the quantum state corresponding to an outcome of 0 is represented by:

_images/bloch_1.png

Notice that the two axes in the horizontal plane have been labeled \(x\) and \(y\), implying that \(z\) is the vertical axis (not labeled). Physicists use the convention that a qubit’s \(\{|\,0\rangle, |\,1\rangle\}\) states are the positive and negative unit vectors along the z axis, respectively. These axes will be useful later in this document.

Multiple qubits are represented in precisely the same way, but taking tensor products of the spaces and states. Thus \(n\) qubits have \(2^n\) possible states.

An Important Distinction

An important distinction between the probabilistic case described above and the quantum case is that probabilistic states may just mask out ignorance. For example a coin is physically only 0 or 1 and the probabilistic view merely represents our ignorance about which it actually is. This is not the case in quantum mechanics. Assuming events cannot instantaneously influence one another, the quantum states — as far as we know — cannot mask any underlying state. This is what people mean when they say that there is no local hidden variable theory for quantum mechanics. These probabilistic quantum states are as real as it gets: they don’t describe our knowledge of the quantum system, they describe the physical reality of the system.

Some Code

Let us take a look at some code in pyQuil to see how these quantum states play out. We will dive deeper into quantum operations and pyQuil in the following sections. Note that in order to run these examples you will need to install pyQuil and set up a connection to the Forest API. Each of the code snippets below will be immediately followed by its output.

# Imports for pyQuil (ignore for now)
import numpy as np
from pyquil.quil import Program
from pyquil.api import QVMConnection
quantum_simulator = QVMConnection()

# pyQuil is based around operations (or gates) so we will start with the most
# basic one: the identity operation, called I. I takes one argument, the index
# of the qubit that it should be applied to.
from pyquil.gates import I

# Make a quantum program that allocates one qubit (qubit #0) and does nothing to it
p = Program(I(0))

# Quantum states are called wavefunctions for historical reasons.
# We can run this basic program on our connection to the simulator.
# This call will return the state of our qubits after we run program p.
# This api call returns a tuple, but we'll ignore the second value for now.
wavefunction = quantum_simulator.wavefunction(p)

# wavefunction is a Wavefunction object that stores a quantum state as a list of amplitudes
alpha, beta = wavefunction

print("Our qubit is in the state alpha={} and beta={}".format(alpha, beta))
print("The probability of measuring the qubit in outcome 0 is {}".format(abs(alpha)**2))
print("The probability of measuring the qubit in outcome 1 is {}".format(abs(beta)**2))

Our qubit is in the state alpha=(1+0j) and beta=0j
The probability of measuring the qubit in outcome 0 is 1.0
The probability of measuring the qubit in outcome 1 is 0.0

Applying an operation to our qubit affects the probability of each outcome.

# We can import the qubit "flip" operation, called X, and see what it does.
# We will learn more about this operation in the next section.
from pyquil.gates import X

p = Program(X(0))

wavefunc = quantum_simulator.wavefunction(p)
alpha, beta = wavefunc

print("Our qubit is in the state alpha={} and beta={}".format(alpha, beta))
print("The probability of measuring the qubit in outcome 0 is {}".format(abs(alpha)**2))
print("The probability of measuring the qubit in outcome 1 is {}".format(abs(beta)**2))

Our qubit is in the state alpha=0j and beta=(1+0j)
The probability of measuring the qubit in outcome 0 is 0.0
The probability of measuring the qubit in outcome 1 is 1.0

In this case we have flipped the probability of outcome 0 into the probability of outcome 1 for our qubit. We can also investigate what happens to the state of multiple qubits. We’d expect the state of multiple qubits to grow exponentially in size, as their vectors are tensored together.

# Multiple qubits also produce the expected scaling of the state.
p = Program(I(0), I(1))
wavefunction = quantum_simulator.wavefunction(p)
print("The quantum state is of dimension:", len(wavefunction.amplitudes))

p = Program(I(0), I(1), I(2), I(3))
wavefunction = quantum_simulator.wavefunction(p)
print("The quantum state is of dimension:", len(wavefunction.amplitudes))

p = Program()
for x in range(10):
    p += I(x)
wavefunction = quantum_simulator.wavefunction(p)
print("The quantum state is of dimension:", len(wavefunction.amplitudes)  )

The quantum state is of dimension: 4
The quantum state is of dimension: 16
The quantum state is of dimension: 1024

Let’s look at the actual value for the state of two qubits combined. The resulting dictionary of this method contains outcomes as keys and the probabilities of those outcomes as values.

# wavefunction(Program) returns a coefficient array that corresponds to outcomes in the following order
wavefunction = quantum_simulator.wavefunction(Program(I(0), I(1)))
print(wavefunction.get_outcome_probs())

{'00': 1.0, '01': 0.0, '10': 0.0, '11': 0.0}

Qubit Operations

In the previous section we introduced our first two operations: the I (or identity) operation and the X operation. In this section we will get into some more details on what these operations are.

Quantum states are complex vectors on the Bloch sphere, and quantum operations are matrices with two properties:

  1. They are reversible.
  2. When applied to a state vector on the Bloch sphere, the resulting vector is also on the Bloch sphere.

Matrices that satisfy these two properties are called unitary matrices. Applying an operation to a quantum state is the same as multiplying a vector by one of these matrices. Such an operation is called a gate.

Since individual qubits are two-dimensional vectors, operations on individual qubits are 2x2 matrices. The identity matrix leaves the state vector unchanged:

\[\begin{split}I = \left(\begin{matrix} 1 & 0\\ 0 & 1 \end{matrix}\right)\end{split}\]

so the program that applies this operation to the zero state is just

\[\begin{split} I\,|\,0\rangle = \left(\begin{matrix} 1 & 0\\ 0 & 1 \end{matrix}\right)\left(\begin{matrix} 1 \\ 0 \end{matrix}\right) = \left(\begin{matrix} 1 \\ 0 \end{matrix}\right) = |\,0\rangle\end{split}\]
p = Program(I(0))
print(quantum_simulator.wavefunction(p))

(1+0j)|0>
Pauli Operators

Let’s revisit the X gate introduced above. It is one of three important single-qubit gates, called the Pauli operators:

\[\begin{split}X = \left(\begin{matrix} 0 & 1\\ 1 & 0 \end{matrix}\right) \qquad Y = \left(\begin{matrix} 0 & -i\\ i & 0 \end{matrix}\right) \qquad Z = \left(\begin{matrix} 1 & 0\\ 0 & -1 \end{matrix}\right)\end{split}\]
from pyquil.gates import X, Y, Z

p = Program(X(0))
wavefunction = quantum_simulator.wavefunction(p)
print("X|0> = ", wavefunction)
print("The outcome probabilities are", wavefunction.get_outcome_probs())
print("This looks like a bit flip.\n")

p = Program(Y(0))
wavefunction = quantum_simulator.wavefunction(p)
print("Y|0> = ", wavefunction)
print("The outcome probabilities are", wavefunction.get_outcome_probs())
print("This also looks like a bit flip.\n")

p = Program(Z(0))
wavefunction = quantum_simulator.wavefunction(p)
print("Z|0> = ", wavefunction)
print("The outcome probabilities are", wavefunction.get_outcome_probs())
print("This state looks unchanged.")

X|0> =  (1+0j)|1>
The outcome probabilities are {'0': 0.0, '1': 1.0}
This looks like a bit flip.

Y|0> =  1j|1>
The outcome probabilities are {'0': 0.0, '1': 1.0}
This also looks like a bit flip.

Z|0> =  (1+0j)|0>
The outcome probabilities are {'0': 1.0, '1': 0.0}
This state looks unchanged.

The Pauli matrices have a visual interpretation: they perform 180-degree rotations of qubit state vectors on the Bloch sphere. They operate about their respective axes as shown in the Bloch sphere depicted above. For example, the X gate performs a 180-degree rotation about the \(x\) axis. This explains the results of our code above: for a state vector initially in the +\(z\) direction, both X and Y gates will rotate it to -\(z\), and the Z gate will leave it unchanged.

However, notice that while the X and Y gates produce the same outcome probabilities, they actually produce different states. These states are not distinguished if they are measured immediately, but they produce different results in larger programs.

Quantum programs are built by applying successive gate operations:

# Composing qubit operations is the same as multiplying matrices sequentially
p = Program(X(0), Y(0), Z(0))
wavefunction = quantum_simulator.wavefunction(p)

print("ZYX|0> = ", wavefunction)
print("With outcome probabilities\n", wavefunction.get_outcome_probs())

ZYX|0> =  [ 0.-1.j  0.+0.j]
With outcome probabilities
{'0': 1.0, '1': 0.0}
Multi-Qubit Operations

Operations can also be applied to composite states of multiple qubits. One common example is the controlled-NOT or CNOT gate that works on two qubits. Its matrix form is:

\[\begin{split}CNOT = \left(\begin{matrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 \\ \end{matrix}\right)\end{split}\]

Let’s take a look at how we could use a CNOT gate in pyQuil.

from pyquil.gates import CNOT

p = Program(CNOT(0, 1))
wavefunction = quantum_simulator.wavefunction(p)
print("CNOT|00> = ", wavefunction)
print("With outcome probabilities\n", wavefunction.get_outcome_probs())

p = Program(X(0), CNOT(0, 1))
wavefunction = quantum_simulator.wavefunction(p)
print("CNOT|01> = ", wavefunction)
print("With outcome probabilities\n", wavefunction.get_outcome_probs())

p = Program(X(1), CNOT(0, 1))
wavefunction = quantum_simulator.wavefunction(p)
print("CNOT|10> = ", wavefunction)
print("With outcome probabilities\n", wavefunction.get_outcome_probs())

p = Program(X(0), X(1), CNOT(0, 1))
wavefunction = quantum_simulator.wavefunction(p)
print("CNOT|11> = ", wavefunction)
print("With outcome probabilities\n", wavefunction.get_outcome_probs())

CNOT|00> =  (1+0j)|00>
With outcome probabilities
 {'00': 1.0, '01': 0.0, '10': 0.0, '11': 0.0}

CNOT|01> =  (1+0j)|11>
With outcome probabilities
 {'00': 0.0, '01': 0.0, '10': 0.0, '11': 1.0}

CNOT|10> =  (1+0j)|10>
With outcome probabilities
 {'00': 0.0, '01': 0.0, '10': 1.0, '11': 0.0}

CNOT|11> =  (1+0j)|01>
With outcome probabilities
 {'00': 0.0, '01': 1.0, '10': 0.0, '11': 0.0}

The CNOT gate does what its name implies: the state of the second qubit is flipped (negated) if and only if the state of the first qubit is 1 (true).

Another two-qubit gate example is the SWAP gate, which swaps the \( |01\rangle \) and \(|10\rangle \) states:

\[\begin{split}SWAP = \left(\begin{matrix} 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ \end{matrix}\right)\end{split}\]
from pyquil.gates import SWAP
p = Program(X(0), SWAP(0,1))
wavefunction = quantum_simulator.wavefunction(p)

print("SWAP|01> = ", wavefunction)
print("With outcome probabilities\n", wavefunction.get_outcome_probs())

SWAP|01> =  (1+0j)|10>
With outcome probabilities
 {'00': 0.0, '01': 0.0, '10': 1.0, '11': 0.0}

In summary, quantum computing operations are composed of a series of complex matrices applied to complex vectors. These matrices must be unitary (meaning that their complex conjugate transpose is equal to their inverse) because the overall probability of all outcomes must always sum to one.

The Quantum Abstract Machine

We now have enough background to introduce the programming model that underlies Quil. This is a hybrid quantum-classical model in which \(N\) qubits interact with \(M\) classical bits:

_images/qam.png

These qubits and classical bits come with a defined gate set, e.g. which gate operations can be applied to which qubits. Different kinds of quantum computing hardware place different limitations on what gates can be applied, and the fixed gate set represents these limitations.

Full details on the Quantum Abstract Machine and Quil can be found in the Quil whitepaper.

The next section on measurements will describe the interaction between the classical and quantum parts of a Quantum Abstract Machine (QAM).

Qubit Measurements

Measurements have two effects:

  1. They project the state vector onto one of the basic outcomes
  2. (optional) They store the outcome of the measurement in a classical bit.

Here’s a simple example:

# Create a program that stores the outcome of measuring qubit #0 into classical register [0]
classical_register_index = 0
p = Program(I(0)).measure(0, classical_register_index)

Up until this point we have used the quantum simulator to cheat a little bit — we have actually looked at the wavefunction that comes back. However, on real quantum hardware, we are unable to directly look at the wavefunction. Instead we only have access to the classical bits that are affected by measurements. This functionality is emulated by the run command.

# Choose which classical registers to look in at the end of the computation
classical_regs = [0, 1]
print(quantum_simulator.run(p, classical_regs))

[[0, 0]]

We see that both registers are zero. However, if we had flipped the qubit before measurement then we obtain:

classical_register_index = 0
p = Program(X(0)) # Flip the qubit
p.measure(0, classical_register_index) # Measure the qubit

classical_regs = [0, 1]
print(quantum_simulator.run(p, classical_regs))

[[1, 0]]

These measurements are deterministic, e.g. if we make them multiple times then we always get the same outcome:

classical_register_index = 0
p = Program(X(0)) # Flip the qubit
p.measure(0, classical_register_index) # Measure the qubit

classical_regs = [0]
trials = 10
print(quantum_simulator.run(p, classical_regs, trials))

[[1], [1], [1], [1], [1], [1], [1], [1], [1], [1]]
Classical/Quantum Interaction

However this is not the case in general — measurements can affect the quantum state as well. In fact, measurements act like projections onto the outcome basis states. To show how this works, we first introduce a new single-qubit gate, the Hadamard gate. The matrix form of the Hadamard gate is:

\[\begin{split}H = \frac{1}{\sqrt{2}}\left(\begin{matrix} 1 & 1\\ 1 & -1 \end{matrix}\right)\end{split}\]

The following pyQuil code shows how we can use the Hadamard gate:

from pyquil.gates import H

# The Hadamard produces what is called a superposition state
coin_program = Program(H(0))
wavefunction = quantum_simulator.wavefunction(coin_program)

print("H|0> = ", wavefunction)
print("With outcome probabilities\n", wavefunction.get_outcome_probs())

H|0> =  (0.7071067812+0j)|0> + (0.7071067812+0j)|1>
With outcome probabilities
{'0': 0.49999999999999989, '1': 0.49999999999999989}

A qubit in this state will be measured half of the time in the \( |0\rangle \) state, and half of the time in the \( |1\rangle \) state. In a sense, this qubit truly is a random variable representing a coin. In fact, there are many wavefunctions that will give this same operational outcome. There is a continuous family of states of the form

\[\frac{1}{\sqrt{2}}\left(|\,0\rangle + e^{i\theta}|\,1\rangle\right)\]

that represent the outcomes of an unbiased coin. Being able to work with all of these different new states is part of what gives quantum computing extra power over regular bits.

# Introduce measurement
classical_reg = 0
coin_program = Program(H(0)).measure(0, classical_reg)
trials = 10

# We see probabilistic results of about half 1's and half 0's
print(quantum_simulator.run(coin_program, [0], trials))

[[0], [1], [1], [0], [1], [0], [0], [1], [0], [0]]

pyQuil allows us to look at the wavefunction after a measurement as well:

classical_reg = 0
coin_program = Program(H(0))
print("Before measurement: H|0> = ", quantum_simulator.wavefunction(coin_program))

coin_program.measure(0, classical_reg)
for x in range(5):
    print("After measurement: ", quantum_simulator.wavefunction(coin_program))

Before measurement: H|0> =  [ 0.70710678+0.j  0.70710678+0.j]

After measurement:  (1+0j)|1>
After measurement:  (1+0j)|0>
After measurement:  (1+0j)|0>
After measurement:  (1+0j)|1>
After measurement:  (1+0j)|1>

We can clearly see that measurement has an effect on the quantum state independent of what is stored classically. We begin in a state that has a 50-50 probability of being \( |0\rangle \) or \( |1\rangle \). After measurement, the state changes into being entirely in \( |0\rangle \) or entirely in \( |1\rangle \) according to which outcome was obtained. This is the phenomenon referred to as the collapse of the wavefunction. Mathematically, the wavefunction is being projected onto the vector of the obtained outcome and subsequently rescaled to unit norm.

# This happens with bigger systems too
classical_reg = 0

# This program prepares something called a Bell state (a special kind of "entangled state")
bell_program = Program(H(0), CNOT(0, 1))
wavefunction = quantum_simulator.wavefunction(bell_program)
print("Before measurement: Bell state = ", wavefunction)

bell_program.measure(0, classical_reg)
for x in range(5):
    wavefunction = quantum_simulator.wavefunction(bell_program)
    print("After measurement: ", wavefunction.get_outcome_probs())

Before measurement: Bell state =  (0.7071067812+0j)|00> + (0.7071067812+0j)|11>

After measurement:  {'00': 1.0, '01': 0.0, '10': 0.0, '11': 0.0}
After measurement:  {'00': 0.0, '01': 0.0, '10': 0.0, '11': 1.0}
After measurement:  {'00': 1.0, '01': 0.0, '10': 0.0, '11': 0.0}
After measurement:  {'00': 1.0, '01': 0.0, '10': 0.0, '11': 0.0}
After measurement:  {'00': 0.0, '01': 0.0, '10': 0.0, '11': 1.0}

The above program prepares entanglement because, even though there are random outcomes, after every measurement both qubits are in the same state. They are either both \( |0\rangle \) or both \( |1\rangle \). This special kind of correlation is part of what makes quantum mechanics so unique and powerful.

Classical Control

There are also ways of introducing classical control of quantum programs. For example, we can use the state of classical bits to determine what quantum operations to run.

true_branch = Program(X(7)) # if branch
false_branch = Program(I(7)) # else branch

# Branch on classical reg [1]
p = Program(X(0)).measure(0, 1).if_then(1, true_branch, false_branch)

# Measure qubit #7 into classical register [7]
p.measure(7, 7)

# Run and check register [7]
print(quantum_simulator.run(p, [7]))

[[1]]

A [1] here means that qubit 7 was indeed flipped.

_images/branch.png
Example: The Probabilistic Halting Problem

A fun example is to create a program that has an exponentially increasing chance of halting, but that may run forever!

inside_loop = Program(H(0)).measure(0, 1)

p = Program().inst(X(0)).while_do(1, inside_loop)

# Run and check register [1]
print(quantum_simulator.run(p, [1]))

[[0]]
_images/loop.png

Next Steps

We hope that you have enjoyed your whirlwind tour of quantum computing. You are now ready to check out the Installation and Getting Started guide!

If you would like to learn more, Nielsen and Chuang’s Quantum Computation and Quantum Information is a particularly excellent resource for newcomers to the field.

If you’re interested in learning about the software behind quantum computing, take a look at our blog posts on The Quantum Software Challenge.

Program

class pyquil.quil.Program(*instructions)[source]

A list of pyQuil instructions that comprise a quantum program.

>>> from pyquil import Program
>>> from pyquil.gates import *
>>> p = Program()
>>> p += H(0)
>>> p += CNOT(0, 1)

Attributes

instructions Fill in any placeholders and return a list of quil AbstractInstructions.
defined_gates A list of defined gates on the program.
out() Serializes the Quil program to a string suitable for submitting to the QVM or QPU.
get_qubits([indices]) Returns all of the qubit indices used in this program, including gate applications and allocated qubits.
is_protoquil() Protoquil programs may only contain gates, Pragmas, and an initial global RESET.

Program Construction

__iadd__(other) Concatenate two programs together using +=, returning a new one.
__add__(other) Concatenate two programs together, returning a new one.
inst(*instructions) Mutates the Program object by appending new instructions.
gate(name, params, qubits) Add a gate to the program.
defgate(name, matrix[, parameters]) Define a new static gate.
define_noisy_gate(name, qubit_indices, kraus_ops) Overload a static ideal gate with a noisy one defined in terms of a Kraus map.
define_noisy_readout(qubit, p00, p11) For this program define a classical bit flip readout error channel parametrized by p00 and p11.
no_noise() Prevent a noisy gate definition from being applied to the immediately following Gate instruction.
measure(qubit_index[, classical_reg]) Measures a qubit at qubit_index and puts the result in classical_reg
reset([qubit_index]) Reset all qubits or just a specific qubit at qubit_index.
measure_all(*qubit_reg_pairs) Measures many qubits into their specified classical bits, in the order they were entered.
alloc() Get a new qubit.
declare(name[, memory_type, memory_size, …]) DECLARE a quil variable
wrap_in_numshots_loop(shots) Wraps a Quil program in a loop that re-runs the same program many times.

Control Flow

while_do(classical_reg, q_program) While a classical register at index classical_reg is 1, loop q_program
if_then(classical_reg, if_program[, …]) If the classical register at index classical reg is 1, run if_program, else run else_program.

Utility Methods

copy() Perform a shallow copy of this program.
pop() Pops off the last instruction.
dagger([inv_dict, suffix]) Creates the conjugate transpose of the Quil program.
__getitem__(index) Allows indexing into the program to get an action.

Utility Functions

pyquil.quil.get_default_qubit_mapping(program)[source]

Takes a program which contains qubit placeholders and provides a mapping to the integers 0 through N-1.

The output of this function is suitable for input to address_qubits().

Parameters:program – A program containing qubit placeholders
Returns:A dictionary mapping qubit placeholder to an addressed qubit from 0 through N-1.
pyquil.quil.address_qubits(program, qubit_mapping=None)[source]

Takes a program which contains placeholders and assigns them all defined values.

Either all qubits must be defined or all undefined. If qubits are undefined, you may provide a qubit mapping to specify how placeholders get mapped to actual qubits. If a mapping is not provided, integers 0 through N are used.

This function will also instantiate any label placeholders.

Parameters:
  • program – The program.
  • qubit_mapping – A dictionary-like object that maps from QubitPlaceholder to Qubit or int (but not both).
Returns:

A new Program with all qubit and label placeholders assigned to real qubits and labels.
pyquil.quil.instantiate_labels(instructions)[source]

Takes an iterable of instructions which may contain label placeholders and assigns them all defined values.

Returns:list of instructions with all label placeholders assigned to real labels.
pyquil.quil.implicitly_declare_ro(instructions)[source]

Implicitly declare a register named ro for backwards compatibility with Quil 1.

There used to be one un-named hunk of classical memory. Now there are variables with declarations. Instead of:

MEASURE 0 [0]

You must now measure into a named register, idiomatically:

MEASURE 0 ro[0]

The MEASURE instruction will emit this (with a deprecation warning) if you’re still using bare integers for classical addresses. However, you must also declare memory in the new scheme:

DECLARE ro BIT[8]
MEASURE 0 ro[0]

This method will determine if you are in “backwards compatibility mode” and will declare a read-out ro register for you. If you program contains any DECLARE commands or if it does not have any MEASURE x ro[x], this will not do anything.

This behavior is included for backwards compatibility and will be removed in future releases of PyQuil. Please DECLARE all memory including ro.

pyquil.quil.merge_with_pauli_noise(prog_list, probabilities, qubits)[source]

Insert pauli noise channels between each item in the list of programs. This noise channel is implemented as a single noisy identity gate acting on the provided qubits. This method does not rely on merge_programs and so avoids the inclusion of redundant Kraus Pragmas that would occur if merge_programs was called directly on programs with distinct noisy gate definitions.

Parameters:
  • prog_list (Iterable[+T_co]) – an iterable such as a program or a list of programs. If a program is provided, a single noise gate will be applied after each gate in the program. If a list of programs is provided, the noise gate will be applied after each program.
  • probabilities (List[~T]) – The 4^num_qubits list of probabilities specifying the desired pauli channel. There should be either 4 or 16 probabilities specified in the order I, X, Y, Z or II, IX, IY, IZ, XI, XX, XY, etc respectively.
  • qubits (List[~T]) – a list of the qubits that the noisy gate should act on.
Returns:

A single program with noisy gates inserted between each element of the program list.
Return type:

Program
pyquil.quil.merge_programs(prog_list)[source]

Merges a list of pyQuil programs into a single one by appending them in sequence. If multiple programs in the list contain the same gate and/or noisy gate definition with identical name, this definition will only be applied once. If different definitions with the same name appear multiple times in the program list, each will be applied once in the order of last occurrence.

Parameters:prog_list (list) – A list of pyquil programs
Returns:a single pyQuil program
Return type:Program
pyquil.quil.get_classical_addresses_from_program(program)[source]

Returns a sorted list of classical addresses found in the MEASURE instructions in the program.

Parameters:program (Program) – The program from which to get the classical addresses.
Return type:Dict[str, List[int]]
Returns:A mapping from memory region names to lists of offsets appearing in the program.
pyquil.quil.percolate_declares(program)[source]

Move all the DECLARE statements to the top of the program. Return a fresh obejct.

Parameters:program (Program) – Perhaps jumbled program.
Return type:Program
Returns:Program with DECLAREs all at the top and otherwise the same sorted contents.
pyquil.quil.validate_protoquil(program)[source]

Ensure that a program is valid ProtoQuil, otherwise raise a ValueError. Protoquil allows a global RESET before any gates, and MEASUREs on each qubit after any gates on that qubit. Pragmas are always allowed, and a final Halt instruction is allowed.

Parameters:program (Program) – The Quil program to validate.
Return type:None

Quantum Computer

pyquil.get_qc(name, *, as_qvm=None, noisy=None, connection=None)[source]

Get a quantum computer.

A quantum computer is an object of type QuantumComputer and can be backed either by a QVM simulator (“Quantum/Quil Virtual Machine”) or a physical Rigetti QPU (“Quantum Processing Unit”) made of superconducting qubits.

You can choose the quantum computer to target through a combination of its name and optional flags. There are multiple ways to get the same quantum computer. The following are equivalent:

>>> qc = get_qc("Aspen-0-12Q-A-noisy-qvm")
>>> qc = get_qc("Aspen-0-12Q-A", as_qvm=True, noisy=True)

and will construct a simulator of the 8q-agave chip with a noise model based on device characteristics. We also provide a means for constructing generic quantum simulators that are not related to a given piece of Rigetti hardware:

>>> qc = get_qc("9q-square-qvm")
>>> qc = get_qc("9q-square", as_qvm=True)

Finally, you can get request a QVM with “no” topology of a given number of qubits (technically, it’s a fully connected graph among the given number of qubits) with:

>>> qc = get_qc("5q-qvm") # or "6q-qvm", or "34q-qvm", ...

These less-realistic, fully-connected QVMs will also be more lenient on what types of programs they will run. Specifically, you do not need to do any compilation. For the other, realistic QVMs you must use qc.compile() or qc.compiler.native_quil_to_executable() prior to qc.run().

Redundant flags are acceptable, but conflicting flags will raise an exception:

>>> qc = get_qc("9q-square-qvm") # qc is fully specified by its name
>>> qc = get_qc("9q-square-qvm", as_qvm=True) # redundant, but ok
>>> qc = get_qc("9q-square-qvm", as_qvm=False) # Error!

Use list_quantum_computers() to retrieve a list of known qc names.

This method is provided as a convenience to quickly construct and use QVM’s and QPU’s. Power users may wish to have more control over the specification of a quantum computer (e.g. custom noise models, bespoke topologies, etc.). This is possible by constructing a QuantumComputer object by hand. Please refer to the documentation on QuantumComputer for more information.

Parameters:
  • name (str) – The name of the desired quantum computer. This should correspond to a name returned by list_quantum_computers(). Names ending in “-qvm” will return a QVM. Names ending in “-noisy-qvm” will return a QVM with a noise model. Otherwise, we will return a QPU with the given name.
  • as_qvm (Optional[bool]) – An optional flag to force construction of a QVM (instead of a QPU). If specified and set to True, a QVM-backed quantum computer will be returned regardless of the name’s suffix
  • noisy (Optional[bool]) – An optional flag to force inclusion of a noise model. If specified and set to True, a quantum computer with a noise model will be returned regardless of the name’s suffix. The noise model for QVM’s based on a real QPU is an empirically parameterized model based on real device noise characteristics. The generic QVM noise model is simple T1 and T2 noise plus readout error. See decoherance_noise_with_asymmetric_ro().
  • connection (Optional[ForestConnection]) – An optional :py:class:ForestConnection` object. If not specified, the default values for URL endpoints, ping time, and status time will be used. Your user id and API key will be read from ~/.pyquil_config. If you deign to change any of these parameters, pass your own ForestConnection object.
Return type:

QuantumComputer
Returns:

A pre-configured QuantumComputer
pyquil.list_quantum_computers(connection=None, qpus=True, qvms=True)[source]

List the names of available quantum computers

Parameters:
  • connection (Optional[ForestConnection]) – An optional :py:class:ForestConnection` object. If not specified, the default values for URL endpoints will be used, and your API key will be read from ~/.pyquil_config. If you deign to change any of these parameters, pass your own ForestConnection object.
  • qpus (bool) – Whether to include QPU’s in the list.
  • qvms (bool) – Whether to include QVM’s in the list.
Return type:

List[str]
class pyquil.api.QuantumComputer(*, name, qam, device, compiler, symmetrize_readout=False)[source]

A quantum computer for running quantum programs.

A quantum computer has various characteristics like supported gates, qubits, qubit topologies, gate fidelities, and more. A quantum computer also has the ability to run quantum programs.

A quantum computer can be a real Rigetti QPU that uses superconducting transmon qubits to run quantum programs, or it can be an emulator like the Rigetti QVM with noise models and mimicked topologies.

Parameters:
  • name (str) – A string identifying this particular quantum computer.
  • qam (QAM) – A quantum abstract machine which handles executing quantum programs. This dispatches to a QVM or QPU.
  • device (AbstractDevice) – A collection of connected qubits and associated specs and topology.
  • symmetrize_readout (bool) – Whether to apply readout error symmetrization. See run_symmetrized_readout() for a complete description.

Methods

run(executable) Run a quil executable.
run_and_measure(program, trials) Run the provided state preparation program and measure all qubits.
run_symmetrized_readout(program, trials) Run a quil program in such a way that the readout error is made collectively symmetric
qubits() Return a sorted list of this QuantumComputer’s device’s qubits
qubit_topology() Return a NetworkX graph representation of this QuantumComputer’s device’s qubit connectivity.
get_isa([oneq_type, twoq_type]) Return a target ISA for this QuantumComputer’s device.
compile(program[, to_native_gates, optimize]) A high-level interface to program compilation.

Compilers

An appropriate compiler is automatically created when using get_qc() and it is stored on the QuantumComputer object as the compiler attribute.

The exact process for compilation depends on whether you’re targeting a QPU or a QVM, and you can conceive of other compilation strategies than those included with pyQuil by default. Therefore, we define an abstract interface that all compilers must follow. See AbstractCompiler for more, or use one of the listed compilers below.

_qac.AbstractCompiler The abstract interface for a compiler.
QVMCompiler(endpoint, device) Client to communicate with the Compiler Server.
LocalQVMCompiler(endpoint, device) Client to communicate with a locally executing quilc instance.
QPUCompiler(endpoint, device) Client to communicate with the Compiler Server.

QAMs

An appropriate QAM is automatically created when using get_qc() and it is stored on the QuantumComputer object as the qam attribute.

The Quantum Abstract Machine (QAM) provides an abstract interface for running hybrid quantum/classical quil programs on either a Quantum Virtual Machine (QVM, a classical simulator) or a Quantum Processor Unit (QPU, a real quantum device).

_qam.QAM() The platonic ideal of this class is as a generic interface describing how a classical computer interacts with a live quantum computer.
QPU(endpoint[, user]) A connection to the QPU.
QVM(connection[, noise_model, gate_noise, …]) A virtual machine that classically emulates the execution of Quil programs.

Devices

An appropriate Device is automatically created when using get_qc() and it is stored on the QuantumComputer object as the device attribute.

There are properties of real quantum computers that go beyond the quantum abstract machine (QAM) abstraction. Real devices have performance specs, limited ISAs, and restricted topologies. AbstractDevice provides an abstract interface for accessing properties of a real quantum device or for mocking out relevant properties for a more realistic QVM.

AbstractDevice
Device(name, raw) A device (quantum chip) that can accept programs.
NxDevice(topology) A shim over the AbstractDevice API backed by a NetworkX graph.

The data structures used are documented here

ISA Basic Instruction Set Architecture specification.
Specs Basic specifications for the device, such as gate fidelities and coherence times.

Utility functions

pyquil.device.isa_from_graph(graph, oneq_type='Xhalves', twoq_type='CZ')[source]

Generate an ISA object from a NetworkX graph.

Parameters:
  • graph (Graph) – The graph
  • oneq_type – The type of 1-qubit gate. Currently ‘Xhalves’
  • twoq_type – The type of 2-qubit gate. One of ‘CZ’ or ‘CPHASE’.
Return type:

ISA
pyquil.device.specs_from_graph(graph)[source]

Generate a Specs object from a NetworkX graph with placeholder values for the actual specs.

Parameters:graph (Graph) – The graph
pyquil.device.isa_to_graph(isa)[source]

Construct a NetworkX qubit topology from an ISA object.

This discards information about supported gates.

Parameters:isa (ISA) – The ISA.
Return type:Graph
pyquil.device.gates_in_isa(isa)[source]

Generate the full gateset associated with an ISA.

Parameters:isa (ISA) – The instruction set architecture for a QPU.
Returns:A sequence of Gate objects encapsulating all gates compatible with the ISA.
Return type:Sequence[Gate]

QVM Man Page

NAME

qvm - a quantum virtual machine for executing Quil

SYNOPSIS

qvm <options> -e # Execute Mode

qvm <options> -S # Server Mode

qvm <options> --benchmark [n] # Benchmark Mode

DESCRIPTION

The Rigetti QVM is a high-performance, classical implementation of a quantum abstract machine. Specifically, it is capable of executing Quil in a variety of ways. The QVM has three main modes of operation: Execute Mode, Server Mode, and Benchmark Mode.

Execute Mode runs the QVM on a single Quil file, printing out information about the execution, as well as a textual representation of the wavefunction if --verbose is provided. (If one would like full access to the wavefunction as an efficient representation, one should use --shared.)

Server Mode runs the QVM as an HTTP server, taking simple POST requests with JSON payloads which are known to the companion library PyQuil. The server is useful even to a single user wishing to run a variety of computations.

Benchmark Mode is used for stress testing the QVM on a computer.

The QVM implements both a Quil interpreter and a just-in-time (JIT) compiler. (Note that “compilation” here does not refer to the sense related to translation of gate sets, but rather translation to machine code for rapid execution.) Interpreted mode is enabled by default and works for all modes of operation. JIT compilation is enabled by supplying the --compile option.

The QVM does not have explicit options for running programs with noise models. Instead, the Quil program itself specifies PRAGMAs for defining Kraus operators and readout POVMs.

OPTIONS

-e, --execute
(Execute Mode) Run the QVM in Execute Mode. Execute the Quil program supplied from stdin and print some information about the course of evaluation.
-S, --server
(Server Mode) Run the QVM in Server Mode. This starts an HTTP server.
--benchmark [<n>]
(Benchmark Mode) Run the QVM benchmark on <n> qubits. The default is 26.
-p <port>, --port <port>
(Server Mode) Run the QVM server on port <port>. The default is 5000.
--memory-limit <num-octets>
Limit the amount of classical memory to <num-octets> octets usable by an individual Quil program. The default is 65536.
-w <n>, --num-workers <n>
Force the number of parallel workers to be <n>. By default, this is the number of logical cores of the host machine.
--time-limit <limit-ms>
(Server Mode) Limit the amount of time for a single request to approximately <limit-ms> milliseconds. By default, this value is 0, which indicates an unlimited amount of time is allowed.
--qubit-limit <n>
(Server Mode) Limit requests to consuming <n> qubits.
--benchmark-type <name>
(Benchmark Mode) Run the benchmark named <name>. Benchmarks include “bell”, “qft”, “hadamard”.
-h, --help
Show the help message.
-v, --version
Show the version.
--verbose

Print every instruction transition of the QVM with information about timing and allocation.

(Execute Mode) Print each basis state and corresponding amplitude of the evolved wavefunction.

-c, --compile
JIT compile the Quil programs to make them run faster.
--safe-include-directory <dir>
Only allow <dir> to be included from with the Quil INCLUDE directive.
--shared <name>
(Server Mode) Run the QVM in shared memory mode. This allocates the wavefunction in POSIX shared memory named <name>. If <name> is an empty string, then a name will be generated. The –qubits argument must be specified.

EXAMPLES

qvm -e < file.quil
Run a Quil file on the QVM.
printf "H 0\nCNOT 0 1\nCNOT 1 2" | qvm --verbose -e
Create a 3-qubit Bell state, printing information about the execution along the way.
qvm -S -p 1234
Start a QVM server for use with PyQuil on port 1234.
qvm -c --benchmark 25 --benchmark-type qft
Benchmark a 25-qubit quantum Fourier transform in compiled mode.

BUGS

Shared memory mode does not work with QVMs executing noisy programs (i.e., ones where Kraus operators or POVMs are specified).

The WAIT instruction does nothing.

SUPPORT

Contact <support@rigetti.com>.

SEE ALSO

quilc(1)

version 0.16.0 (qvm: 0.20.0) [1b48c69] 21 September 2018 QVM(1)

QUILC Man Page

NAME

quilc - an optimizing, architecture-independent Quil compiler

SYNOPSIS

quilc <options>

DESCRIPTION

The Rigetti Quil compiler, quilc, is an optimizing compiler for Quil. It takes a general Quil program along with a qubit architecture, called an ISA, and produces another Quil program that is executable on that architecture. The compiler will also attempt to optimize the program by producing fewer gates and shorter gate depths. The compiler may either be run as a server which takes requests (as is used with PyQuil), or it may be run as a batch program reading from standard input.

Server Mode runs the compiler as an HTTP server, taking simple POST requests with JSON payloads which are known to the companion library pyQuil.

OPTIONS

-S, --server
(Server Mode) Run the compiler in Server Mode. This starts an HTTP server.
-?, -h, --help
Show the help message.
-v, --version
Show the version.
--verbose
Print what the compiler is thinking. (Warning: It thinks a lot.)
--isa <string>
Compile for the qubit architecture defined by <string>, which can be either 8Q, 20Q, 16QMUX, or a path to a QPU description file.
-p, --protoquil
Prescribe that the input and output must be ProtoQuil, which is Quil that is comprised of gates and measurements, with no control flow.
--port <port>
(Server Mode) Run quilc in server mode on port <port>.
-d, --compute-gate-depth
Print a calculated gate depth for the provided circuit as an appended Quil comment. (Requires -p.)
-2, --compute-2Q-gate-depth
Print a calculated multiqubit gate depth for the provided circuit as an appended Quil comment. (Requires -p. Ignores the blacklist and whitelist.)
--compute-gate-volume
Print a calculated gate volume for the provided circuit as an appended Quil comment. (Requires -p.)
-r, --compute-runtime
Print a calculated estimated runtime for the provided circuit as an appended Quil comment. (Requires -p.)
-f, --compute-fidelity
Print a calculated estimated compiled circuit fidelity for the provided circuit as an appended Quil comment. (Requires -p.)
-u, --compute-unused-qubits
Print a list of unused qubits as an appended Quil comment. (Requires -p.)
-t, --show-topological-overhead
Print the number of SWAP gates incurred for topological reasons for the provided circuit as an appended Quil comment. (Requires -p.)
--gate-blacklist <gate-list>
When calculating statistics, ignore the gates present in the comma-separated list of names of <gate-list>.
--gate-whitelist <gate-list>
When calculating statistics, consider only the gates present in the comma-separated list of names of <gate-list>.
--time-limit <limit-ms>
(Server Mode) Limit the amount of time for a single request to approximately <limit-ms> milliseconds. By default, this value is 0, which indicates an unlimited amount of time is allowed.
--without-pretty-printing
Disable pretty printing of numerical quantities (e.g., multiples of pi) in compiled output.
--prefer-gate-ladders
Use gate ladders, instead of the SWAP gate, to implement long-ranged gates, when possible.
-j, --json-serialize
Serialize the output of compilation as a JSON object.
-s, --print-logical-schedule
Include the logically parallelized schedule in JSON output. (Requires -p.)
-m, --compute-matrix-reps
Print the matrix representation of a compiled ProtoQuil program. Additionally, verify that this matrix matches the matrix representation of the input program. (Requires -p. Note that this is a very expensive operation.)
--enable-state-prep-reductions
Perform program optimizations by assuming that the quantum state starts in the zero state.

EXAMPLES

quilc --isa "8Q" < file.quil
Compile a Quil file (printing the result to stdout) for an eight qubit ring.

SUPPORT

Contact <support@rigetti.com>.

SEE ALSO

qvm(1)

0.13.0 (cl-quil: 0.19.0) [e9b41e3] 24 September 2018 QUILC(1)

Indices and Tables