Introduction to Cirq

Cirq is a framework for writing quantum algorithms for noisy intermediate scale quantum (NISQ) devices. Roughly speaking, NISQ devices are those with O(100) qubits that can enact O(1000) gates. Because the resources for NISQ devices are so constrained, we believe that a framework for writing programs on these devices needs to be aware of all of the architectural properties of the device on which the algorithm is written. This is in contrast to other frameworks where there is a clean separation between the abstract model being used and the details of the device.

In this tutorial we will walk through the basics of writing quantum algorithms in Cirq. Our final goal will be to write a variational ansatz for use in an optimization algorithm.

Installing Cirq

To use Cirq one first needs to install Cirq. Installation instructions are available at https://quantumai.google/cirq/start/install. For the purpose of this tutorial, we run pip install cirq as shown in the following code cell to install the latest release of Cirq.

Different notebook execution systems exist, but for most part they have "run" button on a cell which you can click, or "shift + enter" is often the shortcut to run the cell.

try:
    import cirq
except ImportError:
    print("installing cirq...")
    !pip install --quiet cirq
    print("installed cirq.")
    import cirq

import matplotlib.pyplot as plt
import numpy as np

Let's check that Cirq has been successfully installed by importing Cirq and printing out a diagram of Google's Sycamore device.

Google's Sycamore chip

"""Test successful installation by printing out the Sycamore device."""
import cirq_google

print(cirq_google.Sycamore)

(0, 5)───(0, 6)
                                             │        │
                                             │        │
                                    (1, 4)───(1, 5)───(1, 6)───(1, 7)
                                    │        │        │        │
                                    │        │        │        │
                           (2, 3)───(2, 4)───(2, 5)───(2, 6)───(2, 7)───(2, 8)
                           │        │        │        │        │        │
                           │        │        │        │        │        │
                  (3, 2)───(3, 3)───(3, 4)───(3, 5)───(3, 6)───(3, 7)───(3, 8)───(3, 9)
                  │        │        │        │        │        │        │        │
                  │        │        │        │        │        │        │        │
         (4, 1)───(4, 2)───(4, 3)───(4, 4)───(4, 5)───(4, 6)───(4, 7)───(4, 8)───(4, 9)
         │        │        │        │        │        │        │        │
         │        │        │        │        │        │        │        │
(5, 0)───(5, 1)───(5, 2)───(5, 3)───(5, 4)───(5, 5)───(5, 6)───(5, 7)───(5, 8)
         │        │        │        │        │        │        │
         │        │        │        │        │        │        │
         (6, 1)───(6, 2)───(6, 3)───(6, 4)───(6, 5)───(6, 6)───(6, 7)
                  │        │        │        │        │
                  │        │        │        │        │
                  (7, 2)───(7, 3)───(7, 4)───(7, 5)───(7, 6)
                           │        │        │
                           │        │        │
                           (8, 3)───(8, 4)───(8, 5)
                                    │
                                    │
                                    (9, 4)

This cell should run successfully, and the output should in fact be the grid of qubits for the Sycamore device. If so, the install worked!

Qubits, Operations, Moments and Circuits

In Cirq, circuits are represented by a Circuit object. Conceptually:

A Circuit is a collection of Moments.
A Moment is a collection of Operations that all act during the same abstract time slice.
An Operation is an effect that operates on a specific subset of Qubits.
- The most common type of Operation is a Gate applied to several qubits (a "GateOperation").
The Qubits of a circuit are implicitly defined by the operations - you can't allocate qubits to a Circuit.

These ideas are illustrated by the above diagram, where time goes from left to right, horizontal lines are different qubits, operations acting on qubits are the boxes (sometimes spanning multiple qubits) and a moment is a group of operations that "happen at the same time".

Create a `Circuit`

A typical way to create a Circuit is shown below.

"""Creating a circuit."""
# Define three qubits.
a = cirq.NamedQubit("a")
b = cirq.NamedQubit("b")
c = cirq.NamedQubit("c")

# Define a list of operations.
ops = [cirq.H(a), cirq.H(b), cirq.CNOT(b, c), cirq.H(b)]

# Create a circuit from the list of operations.
circuit = cirq.Circuit(ops)
print("Circuit:\n")
print(circuit)

Circuit:

a: ───H───────────

b: ───H───@───H───
          │
c: ───────X───────

Exercise: Create a Circuit

Write a program to create the following circuit.

          ┌──┐
0: ───H─────@────────
            │
1: ───H────@┼────H───
           ││
2: ────────X┼────────
            │
3: ─────────X────────
          └──┘

Note that the circuit has a total of 3 different moments.

# @title Attempt the solution here

# Define 4 qubits.

# Define a list of operations.

# Create a circuit from the list of operations.

# @title Expand to view the solution
"""Creating a circuit."""
# Define four line qubits. NamedQubits would also work, however this demonstrates a more succint syntax.
q = cirq.LineQubit.range(4)

# Define a list of operations.
ops = [cirq.H(q[0]), cirq.H(q[1]), cirq.CNOT(q[1], q[2]), cirq.CNOT(q[0], q[3]), cirq.H(q[1])]

# Create a circuit from the list of operations.
print("Circuit:\n")
print(cirq.Circuit(ops))

Circuit:

          ┌──┐
0: ───H─────@────────
            │
1: ───H────@┼────H───
           ││
2: ────────X┼────────
            │
3: ─────────X────────
          └──┘

Unpacking the circuit

We can unpack this a bit and see all of the components for the circuit.

The first thing we do is pick some qubits to use. There are many different types of qubits in Cirq, and you can define your own by inheriting from the cirq.Qid class. There's nothing inherently special or magical about these quantum id types such as cirq.NamedQubit. They simply identify what you wish to operate on, which is relevant when you are targeting a specific device. For example, if we were creating a circuit for the Sycamore device and wanted to refer to the qubit in the left-most position, we would use cirq.GridQubit(5, 0). (See the first diagram of the Sycamore device we printed out.) For simplicity, in the previous cell we defined cirq.NamedQubits which are simply qubits that can be identified by a name.

Next, we encounter the object cirq.H which is a Hadamard gate with unitary

\[ H = {1 \over \sqrt{2} } \left[ \begin{array}[cc] & 1 & 1 \\ 1 & -1 \end{array}\right] . \]

In Cirq, cirq.H is an instance of the cirq.HPowGate class, which itself is a subclass of Gate (along with other classes). We can use Cirq to see the unitary matrix of Gate objects as follows.

"""Get the unitary of a gate, here the Hadamard gate."""
cirq.unitary(cirq.H)

array([[ 0.70710678+0.j,  0.70710678+0.j],
       [ 0.70710678+0.j, -0.70710678+0.j]])

We see that this agrees with the unitary for the Hadamard gate above.

Gate objects have the ability to be applied "on" one or more qubits. There are two ways to do this for gates, either using the on method or by directly calling the gate on the qubits as if the gate were a function and the qubits were arguments. For example to apply the H on qubit a we can say cirq.H.on(a) or cirq.H(a).

The result of those expressions is typically a GateOperation object, which is a type of Operation.

Note: In Cirq, there is a strong distinction between Operations and Gates. An Operation is associated with specific qubits and can be put in Circuits. A Gate has unspecified qubits, and will produce an operation when acting on qubits.

Once you have a collection of operations, you can construct a Circuit by passing the operations into the constructor for a Circuit:

ops = [list of operations]
circuit = cirq.Circuit(ops)

The last thing we did in the example code was use the (surprisingly useful) ability to print the circuit as a text diagram.

The diagram is visually helpful, but it doesn't really get into the internal details of how the Circuit is represented. As mentioned, a Circuit is made up of a sequence of Moment objects, and each Moment object is a list of non-overlapping Operations. To see this internal structure, we can iterate over the Moments in the Circuit and print them out.

"""Print out the moments in a circuit."""
print("Circuit:\n")
print(circuit)

# Inspecting individual moments.
print("\nMoments in the circuit:\n")
for i, moment in enumerate(circuit):
    print(f'Moment {i}: \n{moment}')

Circuit:

a: ───H───────────

b: ───H───@───H───
          │
c: ───────X───────

Moments in the circuit:

Moment 0: 
  ╷ None
╶─┼──────
a │ H
  │
b │ H
  │
Moment 1: 
  ╷ None
╶─┼──────
b │ @
  │ │
c │ X
  │
Moment 2: 
  ╷ None
╶─┼──────
b │ H
  │

We see that this circuit consists of three moments. For even more on the underlying structure of a circuit, we can print the circuit's repr. This returns a more detailed (and usually less readable) expression.

"""Print the repr of a circuit."""
print(repr(circuit))

cirq.Circuit([
    cirq.Moment(
        cirq.H(cirq.NamedQubit('a')),
        cirq.H(cirq.NamedQubit('b')),
    ),
    cirq.Moment(
        cirq.CNOT(cirq.NamedQubit('b'), cirq.NamedQubit('c')),
    ),
    cirq.Moment(
        cirq.H(cirq.NamedQubit('b')),
    ),
])

Although it is less readable, the usefulness of printing the repr is that it includes all the gory details which can be useful when debugging. The repr is also a valid python expression that evaluates to the circuit. For example, if we notice that a circuit generated in some complicated way triggers a bug in a simulator, copy-pasting the generated circuit's repr into a test, and then working from there, is a simple way to decouple the reproduction of the bug from the circuit generation code.

More ways to create `Circuit`s

Above we created a Circuit by passing in a list of operations to its constructor. In Cirq, there are many ways to construct and modify circuits, and each of these is useful in different contexts. Here are a few examples:

Circuit(...): This is the simplest way to make a circuit. Give this method some operations, and out pops a circuit.
append: Circuits are mutable. You can start with an empty circuit = cirq.Circuit() and simply circuit.append(operations) to add on more and more operations .
insert: Instead of appending, you can insert before a particular moment location (labeled by an integer index).

One interesting, and extremely convenient, fact about Circuit(...), append, and insert is that they "auto flatten" whatever you give them. You can give them a list of operations, but you can also give them

a list of lists of operations,
a generator function that sometimes yields tuples of operations and other times yields individual operations,
or just a single operation (without a list around it).

If it can recursively iterated into individual operations, these three methods will take it.

The above idea uses a concept we call an OP_TREE in Cirq. An OP_TREE is not a class, but a contract. The basic idea is that if the input can be iteratively flattened into a list of operations, then the input is an OP_TREE.

The main place where auto-flattening is useful is when you are building a circuit's operations using generators.

Recall that, in Python, functions that have a yield statement are generators. Generators are functions that act as iterators.

In this context, auto-flattening means that generators producing operations for a circuit can simply yield sub-generators (instead of iterating over them and yielding their items). We show an example of this below.

"""Creating a circuit from generator functions."""


def xor_swap(a, b):
    """Swaps two qubits with three CNOTs."""
    yield cirq.CNOT(a, b)  # |a> |b> --> |a> |a ^ b>
    yield cirq.CNOT(b, a)  # |a> |a ^ b> --> |a ^ a ^ b> | a ^ b> = |b>|a^b>
    yield cirq.CNOT(a, b)  # |b> |a ^ b> --> |b>|a ^ b ^ b> = |b> |a>


cirq.Circuit(xor_swap(a, b))

Exercise: Create a circuit to left rotate 5 qubits.

Now that we've learnt how to build a circuit to swap the state of two qubits, use this to build a circuit which left rotates 5 qubits i.e. on applying the circuit on 5 qubits (0 - 4), the state of qubits should change as follows:

0 --> 4 
4 --> 3
3 --> 2
2 --> 1
1 --> 0

# @title Attempt the solution here

# @title Expand to view the solution
def left_rotate(qubits):
    """Rotates qubits to the left."""
    for i in range(len(qubits) - 1):
        a, b = qubits[i : i + 2]
        yield xor_swap(a, b)


# Get five qubits on a line.
line = cirq.LineQubit.range(5)

# Create a circuit which rotates the qubits to the left.
# uncomment the next line to compare with your circuit
# print(cirq.Circuit(left_rotate(line)))

One can see how this method of creating circuits is quite powerful.

Note that cirq.SWAP is a pre-defined gate in Cirq. We used three cirq.CNOTs instead of cirq.SWAP in the above example to demonstrate auto-flattening with generators.

Insert strategies

You may have noticed that there is a hole in what we've explained so far. We have been passing a one-dimensional sequence of operations, but the output is a two-dimensional circuit (a list-of-lists-of-operations). There is a degree of freedom that hasn't been account for. Specifically, how does Cirq choose the moment that each operation will be placed within?

The answer is the concept of a cirq.InsertStrategy. An InsertStrategy defines how Operations are placed in a Circuit when requested to be inserted at a given location. Here a location is identified by the index of the Moment in the Circuit that operations should be placed before.

Note: In the case of Circuit.append this means inserting at the index len(circuit) which is one more than the largest moment index and so represents the end of the circuit.

There are currently four insertion strategies in Cirq:

InsertStrategy.EARLIEST (the default),
InsertStrategy.NEW,
InsertStrategy.INLINE,
InsertStrategy.NEW_THEN_INLINE.

The strategy InsertStrategy.EARLIEST is defined as follows:

InsertStrategy.EARLIEST: Scans backward from the insert location until a moment with operations touching qubits affected by the operation to insert is found. The operation is added into the moment just after that location.

For example, if we first create an Operation in a single moment, and then use InsertStrategy.EARLIEST the Operation can slide back to this first Moment if there is space.

"""Appending operations with InsertStrategy.EARLIEST."""
# Create an empty circuit.
circuit = cirq.Circuit()

# Append an operation.
# Note: InsertStrategy.EARLIEST is used by default if not otherwise specified.
circuit.append([cirq.CZ(a, b)])

# Append more operations.
# Note: InsertStrategy.EARLIEST is used by default if not otherwise specified.
circuit.append([cirq.H(a), cirq.H(b), cirq.H(c)])

# Display the circuit.
print("Circuit:\n")
print(circuit)

Circuit:

a: ───@───H───
      │
b: ───@───H───

c: ───H───────

After creating the first moment with a CZ gate, the second append uses the InsertStrategy.EARLIEST strategy. The H on a and b cannot slide back, while the H on c can and so ends up in the first Moment.

While InsertStrategy.EARLIEST is the default strategy, the second most important strategy is InsertStrategy.NEW_THEN_INLINE, defined as follows:

InsertStrategy.NEW_THEN_INLINE: For the first operation, add it to a new Moment the insertion point. Attempts to add the operation after the first operation to insert into the moment just before the desired insert location. But, if there's already an existing operation affecting any of the qubits touched by the operation to insert, a new moment is created instead and this Moment is the one that is subsequently used for insertions.

To see an example of this strategy, we create a circuit with the same operations but inserting them with a different strategy.

"""Appending operations with InsertStrategy.NEW_THEN_INLINE."""
# Create an empty circuit.
circuit = cirq.Circuit()

# Append an operation.
circuit.append([cirq.CZ(a, b)], strategy=cirq.InsertStrategy.NEW_THEN_INLINE)

# Append more operations.
circuit.append([cirq.H(a), cirq.H(b), cirq.H(c)], strategy=cirq.InsertStrategy.NEW_THEN_INLINE)

# Display the circuit.
print("Circuit:\n")
print(circuit)

Circuit:

a: ───@───H───
      │
b: ───@───H───

c: ───────H───

In contrast to the previous codeblock using InsertStrategy.EARLIEST, we see that the three cirq.H gates appended after the cirq.CZ gate appear in the same moment when we use InsertStrategy.NEW_THEN_INLINE.

Exercise: Create the given circuit using least number of appends

Now that you've learned about InsertStrategys, here is an exercise to validate your understanding. Create, using the least number of appends, the following circuit:

a: ───@───H───────────H───H───
      │
b: ───@───────H───@───H───────
                  │
c: ───H───────────@───────────

Here imagine that you want exactly the moments indicated by the spacing of the circuit so that there are six moments in this circuit.

# @title Attempt the solution here

# @title Expand to view the solution
# Define three qubits.
a = cirq.NamedQubit('a')
b = cirq.NamedQubit('b')
c = cirq.NamedQubit('c')

# Get an empty circuit.
circuit = cirq.Circuit()

# Append these gates using cirq.InsertStrategy.EARLIEST (the default strategy).
circuit.append([cirq.CZ(a, b), cirq.H(c), cirq.H(a)])

# Append these gates using cirq.InsertStrategy.NEW_THEN_INLINE.
circuit.append(
    [cirq.H(b), cirq.CZ(b, c), cirq.H(b), cirq.H(a), cirq.H(a)],
    strategy=cirq.InsertStrategy.NEW_THEN_INLINE,
)

# Display the circuit.
print("Circuit:\n")
print(circuit)

Circuit:

a: ───@───H───────────H───H───
      │
b: ───@───────H───@───H───────
                  │
c: ───H───────────@───────────

Simulations of a Circuit

Now that we know how to construct Circuits in Cirq, let's see how to execute them on a simulator. First we create a simple circuit to simulate in the following cell.

"""Get a circuit to simulate."""


def basic_circuit(measure=True):
    """Returns a simple circuit with some one- and two-qubit gates,
    as well as (optionally) measurements.
    """
    # Gates we will use in the circuit.
    sqrt_x = cirq.X**0.5
    cz = cirq.CZ

    # Yield the operations.
    yield sqrt_x(a), sqrt_x(b)
    yield cz(a, b)
    yield sqrt_x(a), sqrt_x(b)
    if measure:
        yield cirq.measure(a, b)


# Create a circuit including measurements.
circuit = cirq.Circuit(basic_circuit())
print(circuit)

a: ───X^0.5───@───X^0.5───M───
              │           │
b: ───X^0.5───@───X^0.5───M───

The main simulator in Cirq is the cirq.Simulator. The general pattern of simulation is to instantiate this simulator, then pass in a circuit to either the run or simulate methods (more on this below).

"""Example of simulating a circuit in Cirq."""
# Get a simulator.
simulator = cirq.Simulator()

# Pass the circuit to the simulator.run method.
result = simulator.run(circuit, repetitions=1)
print("Measurement results:")
print(result)

Measurement results:
a,b=1, 1

Running this multiple times should result in different measurement results, since the circuit produces a superposition over all computational basis states.

Above we used the run method of the simulator. In Cirq, run methods mimic the actual hardware in that they don't give one access to unphysical objects like the wavefunction. The repetitions argument is how many times to sample from the circuit.

If one wants to get the wavefunction, the simulate methods can be used as shown below.

"""Simulating a circuit with the `simulate` method."""
# Get a circuit without measurements.
circuit = cirq.Circuit(basic_circuit(measure=False))

# Simulate the circuit.
result = simulator.simulate(circuit, qubit_order=[a, b])

# Print the final state vector (wavefunction).
print("State vector:")
print(np.around(result.final_state_vector, 3))

# Print the state vector in Dirac notation.
print("\nDirac notation:")
print(result.dirac_notation())

State vector:
[0.5+0.j  0. +0.5j 0. +0.5j 0.5+0.j ]

Dirac notation:
0.5|00⟩ + 0.5j|01⟩ + 0.5j|10⟩ + 0.5|11⟩

Notice that we passed a qubit_order into the simulate method. This order helps define the order of the kronecker (tensor) product used in the resulting final_state_vector.

Note: The qubit_order argument is optional. When it is omitted, qubits are sorted ascending according to the ordering methods defined by their Python class (for example cirq.NamedQubit sorts lexicographically by name). If there are multiple types of qubits in one circuit, the name of the type is used as a tie breaker.

The simplest qubit_order value you can provide is a list of the qubits in the desired order. Any qubits from the circuit that are not in the list will be ordered using the default __str__ ordering, but come after qubits that are in the list.

Note: Be aware that all qubits in the list are included in the simulation, even if they are not operated on by the circuit.

The mapping from the order of the qubits to the order of the amplitudes in the wave function can be tricky to understand. Basically, it is the same as the ordering used by numpy.kron.

If the state vector is the array

(0.1, 0.2, 0.3, 0.4),

then this is

0.1|00⟩ + 0.2|01⟩ + 0.3|10⟩ + 0.4|11⟩

in Dirac notation. If

qubit order = [a, b]

then |00> means qubit a is in 0 and qubit b is in 0, |01> means qubit a is 0 and qubit b is 1, etc.

Another way to think about the qubit-to-amplitude ordering is as "for loop ordering":

for a in [0, 1]:
    for b in [0, 1]:
        print(a, b)

The first index (the outermost loop) is the slowest to vary.

Repetitions and histograms

As mentioned, the simulator run methods also take an option for repeating the circuit, namely, the repetitions argument. If the measurements in the circuit are terminal and all other operations are unitary, this simulator is optimized to not recompute the state vector before sampling from the circuit.

"""Simulate a circuit using 1000 repetitions."""
# Get a circuit with terminal measurements to simulate.
circuit = cirq.Circuit(basic_circuit())

# Sample from the circuit 1000 times.
result = simulator.run(circuit, repetitions=1000)

# Get a histogram of measurement results.
print(result.histogram(key="a,b"))

# Plot a state histogram of the result.
cirq.plot_state_histogram(result)

Counter({0: 284, 2: 259, 3: 233, 1: 224})
<Axes: title={'center': 'Result State Histogram'}, xlabel='qubit state', ylabel='result count'>

png

Here we have also demonstrated the use of the histogram method on the result which sums over all the different results for all of the different repetitions.

The histogram method can also be given a fold_func argument, in order to group measurement results under some key before counting them up. For example, we can group by whether or not the two measurement results agreed:

print(
    result.histogram(
        key="a,b", fold_func=lambda bits: "agree" if bits[0] == bits[1] else "disagree"
    )
)

Counter({'agree': 517, 'disagree': 483})

The Deutsch-Jozsa Algorithm

The very first indication that quantum computers could be more powerful than classical computers was provided by David Deutsch in his 1985 paper

David Deutsch, "Quantum Theory, the Church-Turing Principle and the Universal Quantum Computer" Proc. R. Soc. Lond. A 400 97–117. http://doi.org/10.1098/rspa.1985.0070

This algorithm was extended by Deutsch and Richard Jozsa to a more convincing algorithmic separation and what is now called the Deutsch-Jozsa algorithm. In this section we will show how to write circuits for the Deutsch algorithm and then as an exercise in using Cirq for algorithms for a small version of the Deutsch-Jozsa algorithm.

Let's begin with the Deutsch algorithm. In Deutsch's algorithm you are given access to a box which computes a one bit boolean function. That is it is a box which takes in a bit and outputs a bit. If we want to be a mathematician or theoretical computer scientist we write the function \(f\) as \(f: \{0, 1\} \rightarrow \{0, 1\}\). There are exactly four such boolean functions which we can write out in a table

\(x\)	\(f_0\)	\(f_1\)	\(f_x\)	\(f_{\bar{x} }\)
0	0	1	0	1
1	0	1	1	0

The first two of these are constant functions, \(f_0\) and \(f_1\). That is they always output a constant value (independent of the input). The other two \(f_x\) and \(f_\bar{x}\) are balanced. Over their inputs \(0\) and \(1\), they have an equal number of \(0\)s and \(1\)s in their truth table.

We can now state Deutsch's problem:

Given access to a one bit input one bit output boolean function, determine by querying the function whether the function is balanced or constant.

It shouldn't take you much to convince yourself that in order to solve this problem classically you need to call the function on both possible input values. The easiest way to see this is just to consider what happens if you query the function on one particular input and notice that, for either input, learning the value of the function does not separate the constant from balanced functions. In summary:

Classically one must query the binary function twice to distinguish the constant function from the balanced function.

Now lets turn to the quantum approach to this problem. There is one bit of book keeping we need to take care of. Above we have described a classical function on bits that is not reversible. That is, knowing the values of the output does not allow us to determine uniquely the value of the input. In order to run this on a quantum computer, however we need to make this computation reversible. A trick for taking a classical non-reversible function and making it "quantum happy" is to compute the value in an extra register and store the input. Suppose we have an \(n\) bit input \(x\) and we are computing a (potentially non-reverisble) boolean function \(f(x)\). Then we can implement this via a Unitary \(U_f\) that acts like on \(n + 1\) qubits

\[ U_f |x\rangle |y\rangle = |x\rangle | y \oplus f(x)\rangle . \]

Here \(\oplus\) is addition modulo \(2\) (XOR) and we have identified how \(U_f\) acts by its action on all computational basis states \(|x\rangle\) (\(n\) input qubits) and \(|y\rangle\) (\(1\) output qubit). To see that this is reversible one can note that applying the transformation twice returns the state to its original form.

Let's see how to implement these functions in Cirq.

\(f_0\) enacts the transform

\[ \begin{eqnarray} |00\rangle &\rightarrow& |00\rangle \\ |01\rangle &\rightarrow& |01\rangle \\ |10\rangle &\rightarrow& |10\rangle \\ |11\rangle &\rightarrow& |11\rangle \\ \end{eqnarray} \]

Well this is just the identity transform, i.e. an empty circuit.

\(f_1\) enacts the transform

\[ \begin{eqnarray} |00\rangle &\rightarrow& |01\rangle \\ |01\rangle &\rightarrow& |00\rangle \\ |10\rangle &\rightarrow& |11\rangle \\ |11\rangle &\rightarrow& |10\rangle \\ \end{eqnarray} \]

This is the cirq.X bit flip gate on the second qubit.

\(f_x\) enacts the transform

\[ \begin{eqnarray} |00\rangle &\rightarrow& |00\rangle \\ |01\rangle &\rightarrow& |01\rangle \\ |10\rangle &\rightarrow& |11\rangle \\ |11\rangle &\rightarrow& |10\rangle \\ \end{eqnarray} \]

This is nothing more than a cirq.CNOT from the first bit to the second bit.

Finally \(f_\bar{x}\) enacts the transform

\[ \begin{eqnarray} |00\rangle &\rightarrow& |01\rangle \\ |01\rangle &\rightarrow& |00\rangle \\ |10\rangle &\rightarrow& |10\rangle \\ |11\rangle &\rightarrow& |11\rangle \\ \end{eqnarray} \]

which is a cirq.CNOT from the first bit to the second bit followed by a cirq.X on the second bit.

We can encapulate these functions into a dictionary from a oracle name to the operations in the circuit needed to enact this function.

"""Store the operations to query each function in a dictionary."""
# Get qubits for the operations to act on.
q0, q1 = cirq.LineQubit.range(2)

# Define the dictionary of operations. The key of each dictionary entry
# is the subscript of the function f in the above explanatory text.
oracles = {
    '0': [],
    '1': [cirq.X(q1)],
    'x': [cirq.CNOT(q0, q1)],
    'notx': [cirq.CNOT(q0, q1), cirq.X(q1)],
}

We now turn to Deutsch's algorithm. Suppose we are given access to the reversible oracle functions we have defined above. By a similar argument for our irreversible classical functions you can show that you cannot distinguish the balanced from the constant functions by using this oracle only once. But now we can ask the question: what if we are allowed to query this box in superposition, i.e. what if we can use the power of quantum computing?

Deutsch was able to show that you could solve this problem now, with quantum computers, using only a single query. To see how this works we need two simple insights.

Suppose that we prepare the second qubit in the superposition state \(|-\rangle=\frac{1}{\sqrt{2} }(|0\rangle-|1\rangle)\) and apply the oracle. Then we can check that

\[ U_f |x\rangle |-\rangle = U_f|x\rangle \frac{1}{\sqrt{2} }(|0\rangle -|1\rangle ) = |x\rangle \frac{1}{\sqrt{2} }(|f(x)\rangle -|f(x) \oplus 1\rangle ) = (-1)^{f(x)} |x\rangle |-\rangle . \]

This is the so called "phase kickback trick". By applying \(U_f\) onto a target which is in superposition, the value of the function ends up showing up in the global phase.

How can we leverage this to distinguish between the constant and balanced functions? Note that for the constant functions the phase that is applied is the same for all inputs \(|x\rangle\), whereas for the balanced functions the phase is different for each value of \(x\). In other words, if we use the phase kickback trick then for each of the oracles we apply the following transform on the first qubit:

\[ \begin{eqnarray} f_0 \rightarrow I, && f_1 \rightarrow -I, && f_x \rightarrow Z, && f_\bar{x} \rightarrow -Z && \end{eqnarray} \]

Now we only need, on the first qubit, to distinguish between the identity gate and the \(Z\) gate. But we can do this by recalling the identity

\[ H Z H = X \]

where \(H\) is the Hamadard gate.

This means that we can turn a phase flip into a bit flip by applying Hadamards before and after the phase flip. If we look at the constant and balanced functions we see that this means that the constant functions will be proportional to \(I\) and the balanced functions will be proportional to \(X\). If we feed in \(|0\rangle\) to this register, then in the first cases we will only see \(|0\rangle\) and in the second case we will only see \(|1\rangle\). In other words we will be able to distinguish constant from balanced using a single query of the oracle.

Let's code this up.

"""Creating the circuit used in Deutsch's algorithm."""


def deutsch_algorithm(oracle):
    """Returns the circuit for Deutsch's algorithm given an input
    oracle, i.e., a sequence of operations to query a particular function.
    """
    yield cirq.X(q1)
    yield cirq.H(q0), cirq.H(q1)
    yield oracle
    yield cirq.H(q0)
    yield cirq.measure(q0)


for key, oracle in oracles.items():
    print(f"Circuit for f_{key}:")
    print(cirq.Circuit(deutsch_algorithm(oracle)), end="\n\n")

Circuit for f_0:
0: ───H───H───M───

1: ───X───H───────

Circuit for f_1:
0: ───H───H───M───

1: ───X───H───X───

Circuit for f_x:
0: ───H───────@───H───M───
              │
1: ───X───H───X───────────

Circuit for f_notx:
0: ───H───────@───H───M───
              │
1: ───X───H───X───X───────

Lets run these circuits a bunch of times to see that the measurement result ends up correctly distinguishing constant from balanced.

"""Simulate each of the circuits."""
simulator = cirq.Simulator()
for key, oracle in oracles.items():
    result = simulator.run(cirq.Circuit(deutsch_algorithm(oracle)), repetitions=10)
    print(f'oracle: f_{key:<4} results: {result}')

oracle: f_0    results: q(0)=0000000000
oracle: f_1    results: q(0)=0000000000
oracle: f_x    results: q(0)=1111111111
oracle: f_notx results: q(0)=1111111111

We interpret the simulation results as follows:

For the first two functions \(f_0\) and \(f_1\), we always measure \(0\). Therefore, we know that these functions are constant.
For the second two functions \(f_x\) and \(f_{\bar{x} }\), we always measure \(1\). Therefore, we know that these functions are balanced.

Exercise: Two Bit Deutsch-Jozsa Algorithm

All boolean functions for one input bit are either constant or balanced. For boolean functions from two input bits not all functions are constant or balanced. There are two constant functions, \(f(x_0, x_1) = 0\) and \(f(x_0, x_1)=1\), while there are \({4 \choose 2} = 6\) balanced functions. The following code gives you the operations for these functions where we take two input qubits and compute the function in the third qubit.

"""Operations to query all possible functions on two bits.
Two of these functions are constant, and six of these functions are balanced.
"""
# Define three qubits to use.
q0, q1, q2 = cirq.LineQubit.range(3)

# Define the operations to query each of the two constant functions.
constant = ([], [cirq.X(q2)])

# Define the operations to query each of the six balanced functions.
balanced = (
    [cirq.CNOT(q0, q2)],
    [cirq.CNOT(q1, q2)],
    [cirq.CNOT(q0, q2), cirq.CNOT(q1, q2)],
    [cirq.CNOT(q0, q2), cirq.X(q2)],
    [cirq.CNOT(q1, q2), cirq.X(q2)],
    [cirq.CNOT(q0, q2), cirq.CNOT(q1, q2), cirq.X(q2)],
)

An extension of Deutsch's orginal algorithm is the Deutsch-Jozsa algorithm, which can distinguish constant from balanced functions like these using a single query to the oracle. The goal of this exercise is to write a quantum circuit that can distinguish these.

# @title Attempt the solution here
"""Exercise: Write a quantum circuit that can distinguish 
constant from balanced functions on two bits.
"""


def your_circuit(oracle):
    # Your code here!
    yield oracle
    # Your code here!
    yield cirq.measure(q2)

You can check your circuit by running the follow cell which simulates the circuit for all oracles.

"""Check your answer by running this cell."""
simulator = cirq.Simulator()

print("\nYour result on constant functions:")
for oracle in constant:
    result = simulator.run(cirq.Circuit(your_circuit(oracle)), repetitions=10)
    print(result)

print("\nYour result on balanced functions:")
for oracle in balanced:
    result = simulator.run(cirq.Circuit(your_circuit(oracle)), repetitions=10)
    print(result)

Your result on constant functions:
q(2)=0000000000
q(2)=1111111111

Your result on balanced functions:
q(2)=0000000000
q(2)=0000000000
q(2)=0000000000
q(2)=1111111111
q(2)=1111111111
q(2)=1111111111

# @title Expand to view the solution
def dj_circuit(oracle):
    # Phase kickback trick.
    yield cirq.X(q2), cirq.H(q2)

    # Get an equal superposition over input bits.
    yield cirq.H(q0), cirq.H(q1)

    # Query the function.
    yield oracle

    # Use interference to get result, put last qubit into |1>.
    yield cirq.H(q0), cirq.H(q1), cirq.H(q2)

    # Use a final OR gate to put result in final qubit.
    yield cirq.X(q0), cirq.X(q1), cirq.CCX(q0, q1, q2)
    yield cirq.measure(q2)

As above, we can check the solution by running the circuit with each of the oracles.

"""Simulate the Deutsch-Jozsa circuit and check the results."""
print("Result on constant functions:")
for oracle in constant:
    result = simulator.run(cirq.Circuit(dj_circuit(oracle)), repetitions=10)
    print(result)

print("\nResult on balanced functions:")
for oracle in balanced:
    result = simulator.run(cirq.Circuit(dj_circuit(oracle)), repetitions=10)
    print(result)

Result on constant functions:
q(2)=0000000000
q(2)=0000000000

Result on balanced functions:
q(2)=1111111111
q(2)=1111111111
q(2)=1111111111
q(2)=1111111111
q(2)=1111111111
q(2)=1111111111

As with the single-bit case (Deutsch's algorithm), we always measure \(0\) for constant functions and always measure \(1\) for balanced functions.

Gates

Cirq comes with a plethora of common quantum gates. Here we show a few of them.

"""Examples of common gates defined in Cirq."""
# Get some qubits.
q0, q1, q2 = cirq.LineQubit.range(3)

# Get a bunch of common gates defined in Cirq.
ops = [
    cirq.X(q0),  # Pauli-X.
    cirq.Y(q1),  # Pauli-Y.
    cirq.Z(q2),  # Pauli-Z.
    cirq.CZ(q0, q1),  # Controlled-Z gate.
    cirq.CNOT(q1, q2),  # Controlled-X gate.
    cirq.H(q0),  # Hadamard gate.
    cirq.T(q1),  # T gate.
    cirq.S(q2),  # S gate.
    cirq.CCZ(q0, q1, q2),  # Controlled CZ gate.
    cirq.SWAP(q0, q1),  # Swap gate.
    cirq.CSWAP(q0, q1, q2),  # Controlled swap gate.
    cirq.CCX(q0, q1, q2),  # Toffoli (CCNOT) gate.
    cirq.ISWAP(q0, q1),  # ISWAP gate.
    cirq.Rx(rads=0.5 * np.pi)(q0),  # Rotation about X.
    cirq.Ry(rads=0.5 * np.pi)(q1),  # Rotation about Y.
    cirq.Rz(rads=0.5 * np.pi)(q2),  # Rotation about Z.
    cirq.X(q0) ** 0.5,  # Sqrt of NOT gate.
]

# Display a circuit with all of these operations.
print(cirq.Circuit(ops))

0: ───X───@───H───────@───×───@───@───iSwap──────Rx(0.5π)───X^0.5───
          │           │   │   │   │   │
1: ───Y───@───@───T───@───×───×───@───iSwap──────Ry(0.5π)───────────
              │       │       │   │
2: ───Z───────X───S───@───────×───X───Rz(0.5π)──────────────────────

For each of these gates, you can figure out how they act on the computational basis by calling cirq.unitary on the gate. For example, to see the unitary of CNOT, we can do:

"""Get the unitary of CNOT."""
print(cirq.unitary(cirq.CNOT))

[[1.+0.j 0.+0.j 0.+0.j 0.+0.j]
 [0.+0.j 1.+0.j 0.+0.j 0.+0.j]
 [0.+0.j 0.+0.j 0.+0.j 1.+0.j]
 [0.+0.j 0.+0.j 1.+0.j 0.+0.j]]

For single qubit gates, we have named gates like cirq.H for the Hadmard gate as well as the single qubit rotation gates defined as follows.

\[ {\tt cirq.Rx(\theta)}: \exp(-i \frac{\theta}{2} X) = cos \frac{\theta}{2} I - i \sin \frac{\theta}{2} X =\left[ \begin{array} ~\cos \frac{\theta}{2} & -i \sin \frac{\theta}{2} \\ -i \sin \frac{\theta}{2} & \cos \frac{\theta}{2}\end{array} \right] \\ \]

\[ {\tt cirq.Ry(\theta)}: \exp(-i \frac{\theta}{2} Y) = cos \frac{\theta}{2} I - i \sin \frac{\theta}{2} Y =\left[ \begin{array} ~\cos \frac{\theta}{2} & -\sin \frac{\theta}{2} \\ \sin \frac{\theta}{2} & \cos \frac{\theta}{2}\end{array} \right] \\ \]

\[ {\tt cirq.Rz(\theta)}: \exp(-i \frac{\theta}{2} Z) = cos \frac{\theta}{2} I - i \sin \frac{\theta}{2} Z =\left[ \begin{array} ~e^{i \frac{\theta}{2} } & 0 \\ 0 & e^{-i \frac{\theta}{2} } \end{array} \right] \\ \]

In addition to cirq.unitary another important method (behind the scenes, anyways) is cirq.apply_unitary. This allows you to apply a unitary gate onto a state. Of course we could have applied the unitary directly to the state, using cirq.unitary. We'll see below in understanding how these methods are implemented that the cirq.apply_unitary can be used to apply the gate more directly onto the state and can save allocations of memory to store the unitary.

If we apply cirq.Rx to a state we can see how it rotates the state. To do this let us introduce a new simulate method simulate_moment_steps. This allows us to simulate the circuit Moment by Moment. At each point we can access the state. For example here we can use this to create a circuit that is a series of small cirq.Rx rotations and plot the probability of measuring the state in the \(|0\rangle\) state:

"""Plot the probability of measuring a qubit in the ground state."""
# Get a qubit.
a = cirq.NamedQubit('a')

# Get a circuit of a bunch of X rotations.
num_angles = 200
circuit = cirq.Circuit([cirq.Rx(rads=np.pi / 50.0)(a) for theta in range(num_angles)])

# List to store probabilities of the ground state.
probs = []

# Step through the simulation results.
for step in simulator.simulate_moment_steps(circuit):
    prob = np.abs(step.state_vector(copy=True)) ** 2
    probs.append(prob[0])

# Plot the probability of the ground state at each simulation step.
plt.style.use('seaborn-v0_8-whitegrid')
plt.plot(probs, 'o')
plt.xlabel("Step")
plt.ylabel("Probability of ground state");

png

Above we have given ourselves direct access to the wave function and calculated the exact probabilities. Suppose we wanted to sample from the wave function at each point instead.

"""Plot the probability of measuring a qubit in the ground state by sampling."""
# Number of times to sample.
repetitions = 100

# List to store the probability of the ground state.
sampled_probs = []

for i, step in enumerate(simulator.simulate_moment_steps(circuit)):
    samples = step.sample([a], repetitions=repetitions)
    prob = np.sum(samples, axis=0)[0] / repetitions
    sampled_probs.append(prob)


# Plot the probability of the ground state at each simulation step.
plt.style.use('seaborn-v0_8-whitegrid')
plt.plot(sampled_probs, 'o')
plt.xlabel("Step")
plt.ylabel("Probability of ground state");

png

Custom gates

Suppose there is a gate that you want Cirq to support, but it is not implemented in Cirq. How do you go about defining your new gate?

Cirq aims to be Pythonic. One way in which it does this is that it relies on Python's protocol pattern. Protocols are similar to interfaces, in that they define a collection of methods that an object must support to implement a protocol, but different in that this requirement is more informal and not a part of a class or interface declaration. An object supports a protocol if it implements the methods that the protocol defines. You're probably familiar with this if you've ever done something like defined your own Container in Python. To do this for an object you simply define the __contains__, __setitem__, and __getitem__ methods on your object, and then you can use this object anywere the Container protocol is supported.

Let's see how this works for defining a custom gate. The gate we will define is a single qubit gate that has only rational amplitudes. This is based on the famous 3, 4, 5 triangle you may remember from a long ago math class: \(3^2 + 4^2 = 5^2\). Using this observation we can construct normalized vectors and a unitary transform using the ratios of \(3\), \(4\), and \(5\):

\[ \zeta =\left[ \begin{array} ~\frac{3}{5} & \frac{4}{5} \\ -\frac{4}{5} & \frac{3}{5} \end{array} \right] \]

Below is a simple implementation of this gate in Cirq. To do this we simply define a class that inherits from cirq.Gate and implements the cirq.SupportsUnitary protocol by implementing the _unitary_(self) method. We also define an optional __str__ representation which Cirq will use when printing this gate out in a circuit diagram.

"""Example of defining a custom gate in Cirq."""


class RationalGate(cirq.Gate):
    def _num_qubits_(self) -> int:
        return 1

    def _unitary_(self):
        return np.array([[3 / 5, 4 / 5], [-4 / 5, 3 / 5]])

    def __str__(self):
        return 'ζ'

We can now use this custom gate just like any other gate in Cirq.

"""Using the custom gate in a circuit."""
a = cirq.NamedQubit('a')
rg = RationalGate()
print(cirq.Circuit(rg(a)))

a: ───ζ───

We can also get its unitary, as shown below, because the RationalGate defines a _unitary_ method.

print(cirq.unitary(rg))

[[ 0.6  0.8]
 [-0.8  0.6]]

Let's check that we can use this gate in a simulation.

"""Simulate a circuit with a custom gate."""
circuit = cirq.Circuit(rg(a))
simulator = cirq.Simulator()
result = simulator.simulate(circuit)
print(result.final_state_vector)

[ 0.6+0.j -0.8+0.j]

Note on simulating circuits with custom gates. The _unitary_ method is extremely inefficient for gates over many qubits. In most cases the method _apply_unitary_ will be used instead, if it is available. This method allows much more fine grained control on how a unitary is applied to a state, but it is harder to implement, for example because it is expected to use the pre-allocated workspace buffer that was given to it. Almost all of the basic gates we have defined in Cirq have this method implemented. If you need to get performant, custom multi-qubit gates, you should implement a custom _apply_unitary_ method for such gates.

Exercise: Custom Controlled Rx gate

Recall that the cirq.Rx gate is a rotation about the \(X\) Pauli axis:

\[ {\tt cirq.Rx(θ)}: \exp(-i \frac{\theta}{2} X) = cos \frac{\theta}{2} I - i \sin \frac{\theta}{2} X =\left[ \begin{array} ~\cos \frac{\theta}{2} & -i \sin \frac{\theta}{2} \\ -i \sin \frac{\theta}{2} & \cos \frac{\theta}{2}\end{array} \right] . \\ \]

As an exercise, create a two-qubit controlled cirq.Rx gate defined as follows:

\[ {\tt CRx(\theta)}: \left[\begin{array} ~1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & \cos \frac{\theta}{2} & -i \sin \frac{\theta}{2} \\ 0 & 0 & -i \sin \frac{\theta}{2} & \cos \frac{\theta}{2} \end{array} \right] . \]

# @title Attempt the solution here
"""Define a custom controlled cirq.rx gate here."""


class CRx(cirq.Gate):
    def __init__(self, theta):
        self.theta = theta

    def _num_qubits_(self):
        return 2

    def _unitary_(self):
        return np.array(
            [
                # Your code here!
            ]
        )


# Print out its unitary.
print(np.around(cirq.unitary(CRx(0.5 * np.pi)), 3))

[]

# @title Expand to view the solution
"""Defining a custom controlled cirq.Rx gate."""


class CRx(cirq.Gate):
    def __init__(self, theta):
        self.theta = theta

    def _num_qubits_(self):
        return 2

    def _unitary_(self):
        return np.array(
            [
                [1, 0, 0, 0],
                [0, 1, 0, 0],
                [0, 0, np.cos(self.theta / 2), -1j * np.sin(self.theta / 2)],
                [0, 0, -1j * np.sin(self.theta / 2), np.cos(self.theta / 2)],
            ]
        )

    def _circuit_diagram_info_(self, args):
        return '@', f'Rx({self.theta / np.pi}π)'


# Print out its unitary.
print(np.around(cirq.unitary(CRx(0.5 * np.pi)), 3))

[[1.   +0.j    0.   +0.j    0.   +0.j    0.   +0.j   ]
 [0.   +0.j    1.   +0.j    0.   +0.j    0.   +0.j   ]
 [0.   +0.j    0.   +0.j    0.707+0.j    0.   -0.707j]
 [0.   +0.j    0.   +0.j    0.   -0.707j 0.707+0.j   ]]

Note that we also define the _circuit_diagram_info_ method which tells Cirq how to display the gate in a circuit diagram. The first string in the tuple is the symbol for the top wire, and the second string in the tuple is the symbol for the bottom wire. We can use this in a circuit to see the diagram info as shown below.

"""Display a circuit with the custom gate."""
# Get qubits.
a = cirq.NamedQubit('a')
b = cirq.NamedQubit('b')

# Display the circuit.
print('Circuit diagram:')
print(cirq.Circuit(CRx(0.25 * np.pi)(a, b)))

Circuit diagram:
a: ───@───────────
      │
b: ───Rx(0.25π)───

Decompose Protocol

Cirq also supports creating custom "composite" gates which can be defined in terms of other "simpler" gates. The custom gate should implement the cirq.SupportsDecompose protocol by implementing _decompose_ method.

Note that implementing _decompose_ on custom gates is convenient because if a gate has its _decompose_ method defined, other protocols, like cirq.SupportsUnitary, are automatically supported by composing the values obtained from the decomposed operations.

Below we show how one can construct a custom gate that can be decomposed into two gates.

"""Example of a custom gate which supports the decompose protocol."""


class HXGate(cirq.Gate):
    def _num_qubits_(self) -> int:
        return 1

    def _decompose_(self, qubits):
        return cirq.H(*qubits), cirq.X(*qubits)

    def __str__(self):
        return 'HX'

Even though the above gate does not define a _unitary_ method, we can use it in a circuit as follows.

"""Use the gate in a circuit."""
HX = HXGate()

a = cirq.NamedQubit('a')
circuit = cirq.Circuit(HX(a))
assert cirq.has_unitary(HX(a))
print(circuit, cirq.unitary(circuit), sep="\n\n")

a: ───HX───

[[ 0.70710678+0.j -0.70710678+0.j]
 [ 0.70710678+0.j  0.70710678+0.j]]

The symbol HX is a single gate, not a product of two gates. We can decompose the HXGate using cirq.decompose as shown below.

"""Decompose the gate."""
print(cirq.Circuit(cirq.decompose(circuit)))

a: ───Y^0.5───X───X───

Note that this not only decomposed the HX gate into H and X, it also decomposed H into Y**0.5 and X. In order to decompose only once, one can use cirq.decompose_once:

"""Decompose the gate once."""
print(cirq.Circuit(cirq.decompose_once(HX(a))))

a: ───H───X───

Gatesets

In many contexts, the notion of what gate you can apply is defined either by the physical hardware you are running or perhaps by the quantum error correcting code you are working with. In quantum computing we typically talk about gate sets and work with respect to a given gate set.

Cirq provides a cirq.Gateset abstraction to create custom gatesets, which is useful for

Describing the set of allowed gates in a human readable format
Validating a given gate / optree against the set of allowed gates

For example, a gateset that supports any x/y/z single qubit rotation and CNOT can be defind as follows:

# Insert a type (eg: cirq.XPowGate) to accept all instances of that type.
# Insert an instance (eg: cirq.CNOT) to accept only one specific instance of the type.
gateset = cirq.Gateset(cirq.XPowGate, cirq.YPowGate, cirq.ZPowGate, cirq.CNOT)

# Valid gates and operations are accepted by the gateset.
assert cirq.CNOT(*cirq.LineQubit.range(2)) in gateset
assert cirq.X**0.5 in gateset

# Arbitrary powers of cirq.CXPowGate are not part of the gateset.
assert cirq.CNOT**0.5 not in gateset

Parameterized Circuits

In addition to circuit gates with fixed values, Cirq also supports parameterized gates with symbolic values via sympy. These are placeholder values, such as sympy.Symbol('x'), that will only be resolved at run-time. For simulators these values are resolved by providing a ParamResolver. A ParamResolver provides a map from the Symbol's name to its assigned value.

Plain Python dictionaries can also be used whenever a ParamResolver is needed.

"""Define a circuit with parameterized gates."""
# Import sympy for parameterized values.
import sympy as sp

# Get qubits to use in the circuit.
a = cirq.NamedQubit("a")
b = cirq.NamedQubit("b")

# Define a parameterized value.
val = sp.Symbol("s")

# Create a circuit.
circuit = cirq.Circuit(cirq.X.on(a) ** val, cirq.X.on(b) ** val)

# Display it.
print("Circuit with parameterized gates:\n")
print(circuit)

Circuit with parameterized gates:

a: ───X^s───

b: ───X^s───

When we simulate this circuit, we must provide a param_resolver as mentioned.

"""Simulate the circuit at multiple parameter values."""
simulator = cirq.Simulator()

# Simulate the circuit for several values of the parameter.
num_params = 5
for y in range(num_params):
    result = simulator.simulate(circuit, param_resolver={"s": y / 4.0})
    print(f"s={y}: {np.around(result.final_state_vector, 2)}\n")

s=0: [1.+0.j 0.+0.j 0.+0.j 0.+0.j]

s=1: [ 0.6 +0.6j   0.25-0.25j  0.25-0.25j -0.1 -0.1j ]

s=2: [0. +0.5j 0.5+0.j  0.5+0.j  0. -0.5j]

s=3: [-0.1 +0.1j   0.25+0.25j  0.25+0.25j  0.6 -0.6j ]

s=4: [0.+0.j 0.+0.j 0.+0.j 1.+0.j]

Here we see that the Symbol is used in two gates, and then the resolver provides this value at run time.

Parameterized values are most useful in defining what we call a cirq.Sweep. A cirq.Sweep is a collection of parameter resolvers. Our cirq.Simulator and cirq.Sampler (for real devices) interfaces support parameter sweeps. The run_sweep method allows the user to a parameterized circuit with each of the configurations. Within each run, the user might specify the number of repetitions as well which determines the sample size for each run. Running parameter sweeps returns a list of cirq.Results, one per set of fixed parameter values and repetitions. Example:

"""Simulate the circuit at multiple parameter values."""
# Get a list of param resolvers.
num_params = 5
resolvers = [cirq.ParamResolver({'s': y / 8.0}) for y in range(num_params)]

# Add measurements to the circuit.
circuit.append([cirq.measure(a), cirq.measure(b)])

# Simulate the circuit using run_sweep.
results = simulator.run_sweep(program=circuit, params=resolvers, repetitions=10)

for i, result in enumerate(results):
    print(f'params: {result.params.param_dict}\n{result}\n')

params: OrderedDict([('s', 0.0)])
a=0000000000
b=0000000000

params: OrderedDict([('s', 0.125)])
a=0000000000
b=0000000000

params: OrderedDict([('s', 0.25)])
a=0010000000
b=0000000010

params: OrderedDict([('s', 0.375)])
a=1000100000
b=0001010001

params: OrderedDict([('s', 0.5)])
a=0100101000
b=0001001110

A very similar method to run_sweep is the sample method. This returns a pandas DataFrame, where each column is a parameter/measurement key and each row is a repetition.

results = simulator.sample(program=circuit, params=resolvers, repetitions=10)

results.describe()

Above we passed in a list of ParamResolvers to the params parameter of run_sweep and sample. But one can also pass in a Sweepable. There are some useful methods for generating Sweepables, for example to generate an equally spaced set of ParamResolvers one can use Linspace

"""Alternative method of getting a sequence of param resolvers."""
linspace = cirq.Linspace(start=0, stop=1.0, length=11, key='x')
for p in linspace:
    print(p)

cirq.ParamResolver({'x': 0.0})
cirq.ParamResolver({'x': 0.1})
cirq.ParamResolver({'x': 0.2})
cirq.ParamResolver({'x': 0.3})
cirq.ParamResolver({'x': 0.4})
cirq.ParamResolver({'x': 0.5})
cirq.ParamResolver({'x': 0.6})
cirq.ParamResolver({'x': 0.7})
cirq.ParamResolver({'x': 0.8})
cirq.ParamResolver({'x': 0.9})
cirq.ParamResolver({'x': 1.0})

Exercise: Rotate a qubit

Let's do the equivalent of a Rabi-flop experiment. That is, let's apply a XPowGate rotating about the X axis for a linearly spaced set of values followed by a computational basis measurement. The end result should be a plot of the sampled fraction that were \(|1\rangle\) as a function of gates of \(X^t\) for \(t\) between 0 and \(1\) for 100 values of \(t\) and each result sampled 100 times.

# @title Attempt the solution here
import pandas

q = cirq.NamedQubit("q")
theta = sp.Symbol("theta")
parameterized_circuit = cirq.Circuit(
    # your code here
    cirq.measure(q)
)
# replace None with something else
param_resolvers = None
repetitions = 100
results = cirq.Simulator().sample(
    program=parameterized_circuit, params=param_resolvers, repetitions=repetitions
)

# You can test with the following plot
# pandas.crosstab(results.theta, results.q).plot()

# @title Expand to view the solution
import pandas

q = cirq.NamedQubit("q")
parameterized_circuit = cirq.Circuit([cirq.Rx(rads=sp.Symbol("theta"))(q), cirq.measure(q)])
param_resolvers = cirq.Linspace(start=0, stop=np.pi, length=100, key='theta')
results = cirq.Simulator().sample(
    program=parameterized_circuit, params=param_resolvers, repetitions=repetitions
)
pandas.crosstab(results.theta, results.q).plot()

<Axes: xlabel='theta'>

png

Noise

In addition to circuits with unitary gates, Cirq also has support for modeling noisy quantum evolutions. This is useful when modeling what will happen when running on actual hardware.

Cirq currently supports noise that fits within the context of operator sum representations of noise (a.k.a quantum operations, quantum dyanamical maps, superoperators, etc). This formalism models the evolution of a density matrix via

\[ \rho \rightarrow \sum_k A_k \rho A_k^\dagger \]

where the \(A_k\) are Kraus operators. These operators are not necessarily unitary and satisfy the property

\[ \sum_k A_k^\dagger A_k = I . \]

An example of a noise operator is the depolarizing channel on one qubit. This takes

\[ \rho \rightarrow (1-p) \rho + \frac{p}{3} (X \rho X + Y \rho Y + Z \rho Z) . \]

In Cirq we can define such a channel and use it in a quantum circuit:

"""Create a circuit with a depolarizing channel."""
circuit = cirq.Circuit(cirq.depolarize(0.2)(a), cirq.measure(a))
print(circuit)

a: ───D(0.2)───M───

Previously we saw that gates could implement that _unitary_ protocol, and by doing so they could be used to perform wave function simulation. For noise the gates implement the _kraus_ protocol. Classes that implement this protocol return the Krauss operators on their _kraus_ method. Thus

for i, kraus in enumerate(cirq.kraus(cirq.depolarize(0.2))):
    print(f"Kraus operator {i} is:", kraus, sep="\n", end="\n\n")

Kraus operator 0 is:
[[0.89442719 0.        ]
 [0.         0.89442719]]

Kraus operator 1 is:
[[0.        +0.j 0.25819889+0.j]
 [0.25819889+0.j 0.        +0.j]]

Kraus operator 2 is:
[[0.+0.j         0.-0.25819889j]
 [0.+0.25819889j 0.+0.j        ]]

Kraus operator 3 is:
[[ 0.25819889+0.j  0.        +0.j]
 [ 0.        +0.j -0.25819889+0.j]]

The Kraus operators are often more conveniently represented in a Pauli basis. We can do this in Cirq as shown below.

for i, kraus in enumerate(cirq.kraus(cirq.depolarize(0.2))):
    pauli_ex = cirq.expand_matrix_in_orthogonal_basis(kraus, cirq.PAULI_BASIS)
    print(f"Kraus operator {i} is:", pauli_ex, sep="\n", end="\n\n")

Kraus operator 0 is:
0.894*I

Kraus operator 1 is:
0.258*X

Kraus operator 2 is:
0.258*Y

Kraus operator 3 is:
0.258*Z

In addition to the wavefunction simulator, Cirq also has a density matrix simulator. Instead of keeping track of the wavefunction, this simulator keeps track of the density matrix. It has the same run and simulate type methods. For example we can use this to simulate depolarizing channel and return the final density matrix of the system.

"""Example of simulating a noisy circuit with the density matrix simulator."""
# Circuit to simulate.
circuit = cirq.Circuit(cirq.depolarize(0.2)(a))
print(f'Circuit:\n{circuit}\n')

# Get the density matrix simulator.
simulator = cirq.DensityMatrixSimulator()

# Simulate the circuit and get the final density matrix.
matrix = simulator.simulate(circuit).final_density_matrix
print(f'Final density matrix:\n{matrix}')

Circuit:
a: ───D(0.2)───

Final density matrix:
[[0.8666666 +0.j 0.        +0.j]
 [0.        +0.j 0.13333333+0.j]]

One thing to note is that the density matrix simulator simulates measurement statistically, and not as a channel where the outcome is not known. Consider the following example.

"""Simulating a circuit with measurements using the DensityMatrixSimulator."""
# Get a circuit with measurements.
circuit = cirq.Circuit(cirq.depolarize(0.5)(a), cirq.measure(a))

# Simulate with the density matrix multiple times.
dmat1 = simulator.simulate(circuit).final_density_matrix
dmat2 = simulator.simulate(circuit).final_density_matrix

print(np.allclose(dmat1, dmat2))

True

Because the final density matrix is statistical due to the measurements, the output of the above cell will change when executed multiple times.

Monte carlo simulations

Density matrix simulations are more expensive than pure state wave function simulations. However some channels allow an interpreation of randomly applying one of a fixed set of unitaries with differing probabilites. For example the depolarizing channel above can be interpretted as:

With probability \(1-p\) apply the identity to the state, and
with probability \(p\) apply one of the three Pauli matrices \(X\), \(Y\), or \(Z\) with equal probability.

Channels that can be interpretted in this form can be simulating using a wavefunction simulator: when this channel is simulated the simulation will sample a unitary with the appropriate probability.

For channels of these type, the channel can, instead of implementing the _kraus_ protocol, implement the _mixture_ protocol:

"""Use the cirq.mixture protocol on the cirq.depolarize channel."""
for p, u in cirq.mixture(cirq.depolarize(0.2)):
    print(f"prob = {p}\nunitary: \n{u}\n")

prob = 0.8
unitary: 
[[1. 0.]
 [0. 1.]]

prob = 0.06666666666666667
unitary: 
[[0.+0.j 1.+0.j]
 [1.+0.j 0.+0.j]]

prob = 0.06666666666666667
unitary: 
[[0.+0.j 0.-1.j]
 [0.+1.j 0.+0.j]]

prob = 0.06666666666666667
unitary: 
[[ 1.+0.j  0.+0.j]
 [ 0.+0.j -1.+0.j]]

In fact the depolarizing channel does not implement _kraus_. Instead it only implements _mixture_ and the cirq.kraus method notices this and derives the kraus operators from the mixture.

"""Check if cirq.depolarize has _kraus_ and _mixture_ methods."""
# Get a depolarizing channel.
d = cirq.depolarize(0.2)

# Check if it has _kraus_ implemented.
print(f"does cirq.depolarize(0.2) have _kraus_? {'yes' if getattr(d, '_kraus_', None) else 'no'}")

# Check if it has _mixture_ implemented.
print(
    f"does cirq.depolarize(0.2) have _mixture_? {'yes' if getattr(d, '_mixture_', None) else 'no'}"
)

does cirq.depolarize(0.2) have _kraus_? no
does cirq.depolarize(0.2) have _mixture_? yes

When channels implement mixture then, as we said, we can use the wavefunction simulator:

"""Use the wavefunction simulator on a channel that implements the mixture protocol."""
circuit = cirq.Circuit(cirq.depolarize(0.5).on(a), cirq.measure(a))
simulator = cirq.Simulator()
result = simulator.run(circuit, repetitions=10)
print(result)

a=1010010010

Because the unitary Kraus operators are applied stochastically, executing the above cell multiple times will produce different outputs.

Adding noise to circuits and simulations

To add noise to circuits or during simulations, we provide the notion of a NoiseModel. A NoiseModel may add noise operation by operation, or it may add noise moment by moment, or it may add noise across a list of moments.

For example we can define a noise model that add a single qubit depolarizing for every qubit in each moment.

"""Adding noise to a circuit."""
# Get a noiseless circuit.
noise = cirq.ConstantQubitNoiseModel(cirq.depolarize(0.2))
circuit = cirq.Circuit(cirq.H(a), cirq.CNOT(a, b), cirq.measure(a, b))
print(f'Circuit with no noise:\n{circuit}\n')

# Add noise to the circuit.
system_qubits = sorted(circuit.all_qubits())
noisy_circuit = cirq.Circuit()
for moment in circuit:
    noisy_circuit.append(noise.noisy_moment(moment, system_qubits))
print(f'Circuit with noise:\n{noisy_circuit}')

Circuit with no noise:
a: ───H───@───M───
          │   │
b: ───────X───M───

Circuit with noise:
a: ───H───D(0.2)[<virtual>]───@───D(0.2)[<virtual>]───M───D(0.2)[<virtual>]───
                              │                       │
b: ───────D(0.2)[<virtual>]───X───D(0.2)[<virtual>]───M───D(0.2)[<virtual>]───

We can also pass a noise model into the cirq.DensityMatrixSimulator and execute a noisy circuit in this manner.

"""Perform noisy simulation by defining a density matrix simulator with a noise model."""
# Define a noise model.
noise = cirq.ConstantQubitNoiseModel(cirq.depolarize(0.2))

# Pass this noise model into the simulator.
simulator = cirq.DensityMatrixSimulator(noise=noise)

# Get a circuit to simulate.
circuit = cirq.Circuit(cirq.H(a), cirq.CNOT(a, b), cirq.measure(a, b))

# Simulate the circuit in steps.
for i, step in enumerate(simulator.simulate_moment_steps(circuit)):
    print(f'After step {i} state was\n{step.density_matrix()}\n')

After step 0 state was
[[0.43333328+0.j 0.        +0.j 0.31777775+0.j 0.        +0.j]
 [0.        +0.j 0.06666666+0.j 0.        +0.j 0.04888888+0.j]
 [0.31777775+0.j 0.        +0.j 0.43333328+0.j 0.        +0.j]
 [0.        +0.j 0.04888888+0.j 0.        +0.j 0.06666666+0.j]]

After step 1 state was
[[0.34859255+0.j 0.        +0.j 0.        +0.j 0.17089382+0.j]
 [0.        +0.j 0.15140739+0.j 0.02629136+0.j 0.        +0.j]
 [0.        +0.j 0.02629136+0.j 0.15140739+0.j 0.        +0.j]
 [0.17089382+0.j 0.        +0.j 0.        +0.j 0.34859255+0.j]]

After step 2 state was
[[0.01777777+0.j 0.        +0.j 0.        +0.j 0.        +0.j]
 [0.        +0.j 0.11555552+0.j 0.        +0.j 0.        +0.j]
 [0.        +0.j 0.        +0.j 0.11555552+0.j 0.        +0.j]
 [0.        +0.j 0.        +0.j 0.        +0.j 0.7511109 +0.j]]

Devices

NISQ algorithms work in a regime where every gate counts. A key philosophy behind Cirq is that we believe the details of the hardware, the performance characteristics, as well as device constraints, will be key to getting the most out of NISQ algorithms. Towards this end these hardware features are contained in the Device class.

For example, here is Google's Sycamore device which we printed out at the start of this notebook.

print(cirq_google.Sycamore)

(0, 5)───(0, 6)
                                             │        │
                                             │        │
                                    (1, 4)───(1, 5)───(1, 6)───(1, 7)
                                    │        │        │        │
                                    │        │        │        │
                           (2, 3)───(2, 4)───(2, 5)───(2, 6)───(2, 7)───(2, 8)
                           │        │        │        │        │        │
                           │        │        │        │        │        │
                  (3, 2)───(3, 3)───(3, 4)───(3, 5)───(3, 6)───(3, 7)───(3, 8)───(3, 9)
                  │        │        │        │        │        │        │        │
                  │        │        │        │        │        │        │        │
         (4, 1)───(4, 2)───(4, 3)───(4, 4)───(4, 5)───(4, 6)───(4, 7)───(4, 8)───(4, 9)
         │        │        │        │        │        │        │        │
         │        │        │        │        │        │        │        │
(5, 0)───(5, 1)───(5, 2)───(5, 3)───(5, 4)───(5, 5)───(5, 6)───(5, 7)───(5, 8)
         │        │        │        │        │        │        │
         │        │        │        │        │        │        │
         (6, 1)───(6, 2)───(6, 3)───(6, 4)───(6, 5)───(6, 6)───(6, 7)
                  │        │        │        │        │
                  │        │        │        │        │
                  (7, 2)───(7, 3)───(7, 4)───(7, 5)───(7, 6)
                           │        │        │
                           │        │        │
                           (8, 3)───(8, 4)───(8, 5)
                                    │
                                    │
                                    (9, 4)

In a future version, we intend for each Device to define a noise model.

Devices also contain more information about the timing of the device. For example here we can calculate the duration of an X on the Sycamore device.

"""Get the duration of an operation."""
op = cirq.X.on(cirq.GridQubit(5, 5))
gate_durations = cirq_google.Sycamore.metadata.gate_durations
for gate_family in gate_durations:
    if op in gate_family:
        print(gate_durations[gate_family])

25 ns

Another property of devices is that they can be used to enforce constraints from the hardware, both checking that these constraints are satisfied, but also enforcing the constraints on the device. For example, on the Sycamore device, a two-qubit gate can be only be performed between adjacent qubits. So, for example, validating an operation on non-adjacent qubits will raise an error.

"""Validate operations on a device."""
# Get non-adjacent qubits on the Sycamore device.
q55 = cirq.GridQubit(5, 5)
q56 = cirq.GridQubit(5, 6)
q66 = cirq.GridQubit(6, 6)

# Operations on adjacent qubits will be validated.
cirq_google.Sycamore.validate_operation(cirq.SQRT_ISWAP(q55, q56))
cirq_google.Sycamore.validate_operation(cirq.SQRT_ISWAP(q56, q66))

# Operation on non-adjacent qubits will raise an error.
ops = [cirq.SQRT_ISWAP(q55, q66)]
circuit = cirq.Circuit(ops)
print(circuit)

try:
    cirq_google.Sycamore.validate_circuit(circuit)
except ValueError as ex:
    print(f"error, as expected: \n{ex}")

(5, 5): ───iSwap───────
           │
(6, 6): ───iSwap^0.5───
error, as expected: 
Qubit pair is not valid on device: (cirq.GridQubit(5, 5), cirq.GridQubit(6, 6)).

Exercise: Make a Device

Construct a device that acts on a square sized lattice, and only allows Hadamard, CZ, and measurement gates.

# @title Attempt the solution here

# @title Expand to view the solution
class SquareDevice(cirq.Device):
    """A Square Grid Device.

    The device that only allows
        1) Grid Qubits from (0, 0) to (grid_size - 1, grid_size - 1)
        2) H, CZ and MeasurementGate gates.
    """

    def __init__(self, grid_size):
        self.qubits = []
        for i in range(grid_size):
            for j in range(grid_size):
                self.qubits.append(cirq.GridQubit(i, j))

    def validate_operation(self, operation: 'cirq.Operation') -> None:
        if not isinstance(operation, cirq.GateOperation):
            raise ValueError(f"Unsupported operation {operation}")
        if not (
            operation.gate in [cirq.H, cirq.CZ] or isinstance(operation.gate, cirq.MeasurementGate)
        ):
            raise ValueError(f"Unsupported gate {operation.gate}")
        for qubit in operation.qubits:
            if qubit not in self.qubits:
                raise ValueError(f"Qubit {qubit} not on device")

Quantum Virtual Machine

Cirq provides the Quantum Virtual Machine as a way to run circuits using the same interface that quantum hardware uses and get results that approximate those the hardware would produce. It combines Devices, Noisy Simulation, realistic Noise Models and a Virtual Engine Interface to produce a cirq.SimulatedLocalEngine that you can treat almost interchangeably with hardware devices. See the Quantum Virtual Machine page to learn more about how it works, or jump to the QVM Creation Template notebook to just start using it.

Transformers for circuit compilation

A transformer in Cirq is any callable, that satisfies the cirq.TRANSFORMER API, and transforms an input circuit into an output circuit. A custom transformer can be created by simply decorating the method/class with @cirq.transformer decorator.

Cirq also provides useful transformer primitives which can be composed together to implement common compilation patterns. For example:

cirq.map_operations: Applies local transformations on operations, by calling map_func(op) for each op.
cirq.merge_operations: Merges connected component of operations by calling merge_func(op1, op2) on iteratively for every mergeable operation op1 and op2.

Below, we give an example of how you can create a circuit optimizer using the above primitives.

"""Example of writing a custom cirq transformer."""


@cirq.transformer
def xz_optimizer(circuit, *, context=None):
    """Replaces an X followed by a Z with a Y."""

    def merge_func(op1, op2):
        return cirq.Y(*op1.qubits) if op1.gate == cirq.X and op2.gate == cirq.Z else None

    return cirq.merge_operations(circuit, merge_func)


circuit = cirq.Circuit(cirq.X(a), cirq.Z(a), cirq.CZ(a, b), cirq.X(a))
print(f"Before optimizing:\n{circuit}\n")
circuit = xz_optimizer(circuit)
print(f"After optimizing:\n{circuit}")

Before optimizing:
a: ───X───Z───@───X───
              │
b: ───────────@───────

After optimizing:
a: ───Y───────@───X───
              │
b: ───────────@───────

Exercise: Simplify flipped CNOTs

Write a transformer, using transformer primitives, which (greedily) performs the simplification that

a: ───H───@───H───
          │
b: ───H───X───H───

is equal to

a: ───X───
      │
b: ───@───

# @title Attempt the solution here
@cirq.transformer
def simplify_flipped_cnots(circuit, *, context=None):
    # Use transformer primitives to simplify your circuit.
    pass


"""Test your optimizer on this circuit."""
circuit = cirq.Circuit(
    cirq.H.on_each(a, b, c), cirq.CNOT(a, b), cirq.H.on_each(a, b), cirq.CZ(a, b)
)
print(simplify_flipped_cnots(circuit))

None

# @title Expand to view the solution
def simplify_flipped_cnots(circuit, *, context=None):
    """Replaces a CX surrounded by Hadamards"""

    def can_merge(ops1, ops2):
        merged = cirq.Circuit(ops1, ops2)
        return (
            len(merged) <= 3
            and all(o.gate == cirq.H for o in merged[0].operations)
            and merged[1].operations[0].gate == cirq.CNOT
            and (len(merged) < 3 or all(o.gate == cirq.H for o in merged[2].operations))
        )

    merged_tag = "h_cx_h"
    circuit = cirq.merge_operations_to_circuit_op(
        circuit, can_merge, merged_circuit_op_tag=merged_tag
    )

    def map_func(op, _):
        if merged_tag not in op.tags:
            return op
        cops = [*op.untagged.mapped_circuit().all_operations()]
        return cirq.CNOT(*cops[2].qubits[::-1]) if len(cops) == 5 else cops

    return cirq.map_operations(circuit, map_func)


circuit = cirq.Circuit(
    cirq.H.on_each(a, b, c), cirq.CNOT(a, b), cirq.H.on_each(a, b), cirq.CZ(a, b)
)
print(f"Before optimizing:\n{circuit}\n")
print(f"After optimizing:\n{simplify_flipped_cnots(circuit)}")

Before optimizing:
a: ───H───@───H───@───
          │       │
b: ───H───X───H───@───

c: ───H───────────────

After optimizing:
a: ───────X───────@───
          │       │
b: ───────@───────@───

c: ───H───────────────

Compiling to Target Gatesets

Cirq's philosophy on compiling circuits for execution on a NISQ target device or simulator is that it would often require running only a handful of individual compilation passes on the input circuit, one after the other.

cirq.CompilationTargetGateset is an abstraction in Cirq to represent such compilation targets as well as the bundles of transformer passes which should be executed to compile a circuit to this target. Cirq has implementations for common target gatesets like cirq.CZTargetGateset, cirq.SqrtIswapTargetGateset etc.

Below, we give an example of compiling a circuit to a specific target gateset.

circuit = cirq.testing.random_circuit(qubits=4, n_moments=6, op_density=0.8, random_state=1234)
print(f"Original Circuit (depth {len(circuit)}):", circuit)
cz_compiled_circuit = cirq.optimize_for_target_gateset(circuit, gateset=cirq.CZTargetGateset())
cirq.testing.assert_circuits_with_terminal_measurements_are_equivalent(
    circuit, cz_compiled_circuit, atol=1e-6
)
print(f"Compiled circuit for CZ Target (depth {len(cz_compiled_circuit)}):", cz_compiled_circuit)
sqrt_iswap_compiled_circuit = cirq.optimize_for_target_gateset(
    circuit, gateset=cirq.SqrtIswapTargetGateset()
)
cirq.testing.assert_circuits_with_terminal_measurements_are_equivalent(
    circuit, sqrt_iswap_compiled_circuit, atol=1e-6
)
print(
    f"Compiled circuit for Sqrt-Iswap Target (depth {len(sqrt_iswap_compiled_circuit)}):",
    sqrt_iswap_compiled_circuit,
)

Original Circuit (depth 6):               ┌──┐                       ┌──┐
0: ───iSwap──────────────iSwap───iSwap────@─────
      │                  │       │        │
1: ───┼─────────×────T───iSwap───iSwap────┼×────
      │         │                         ││
2: ───┼────────T┼────────X────────────────@┼────
      │         │                          │
3: ───iSwap─────×──────────────────────────×────
              └──┘                       └──┘
Compiled circuit for CZ Target (depth 17):                                                                                                                                                                                                         ┌──────────────────┐   ┌───────────────────────┐
0: ───PhXZ(a=0.25,x=-0.5,z=0)─────────@───PhXZ(a=0.25,x=0.5,z=5.55e-17)───@───PhXZ(a=-0.75,x=0.5,z=0.5)────────────────────────────────────────────────────────────────PhXZ(a=0,x=0,z=1)──────────────────────────────────────────────────────────────@─────────────────────────────────────────────────────────────────
                                      │                                   │                                                                                                                                                                           │
1: ───────────────────────────────────┼───────────────────────────────────┼───────────────────────────────@───PhXZ(a=0.5,x=0.5,z=0)───@───PhXZ(a=0.5,x=-0.5,z=0)───@───PhXZ(a=-1.11e-16,x=0.25,z=0.5)────@──────────────────────PhXZ(a=0,x=-0.5,z=0)──┼────@───PhXZ(a=-2.22e-16,x=0.5,z=0)───@───PhXZ(a=0,x=0,z=-0.5)───
                                      │                                   │                               │                           │                            │                                     │                                            │    │                                 │
2: ───────────────────────────────────┼───────────────────────────────────┼───────────────────────────────┼───────────────────────────┼────────────────────────────┼───PhXZ(a=0,x=0,z=0.25)──────────────┼PhXZ(a=0,x=1,z=0)───────────────────────────@────┼─────────────────────────────────┼──────────────────────────
                                      │                                   │                               │                           │                            │                                     │                                                 │                                 │
3: ───PhXZ(a=0.25,x=0.5,z=5.55e-17)───@───PhXZ(a=0.25,x=-0.5,z=0)─────────@───PhXZ(a=0.25,x=0.5,z=0.5)────@───PhXZ(a=0.5,x=0.5,z=0)───@───PhXZ(a=0.5,x=-0.5,z=0)───@───PhXZ(a=-0.5,x=0,z=-1)─────────────@──────────────────────PhXZ(a=0.5,x=-0.5,z=0)─────@───PhXZ(a=0.5,x=0.5,z=0)─────────@───PhXZ(a=0,x=0,z=0.5)────
                                                                                                                                                                                                        └──────────────────┘   └───────────────────────┘
Compiled circuit for Sqrt-Iswap Target (depth 16):                                                                                                                                                                                                            ┌──────────────────┐                                ┌──────────────────┐
0: ───PhXZ(a=-4.44e-16,x=0,z=0.5)───iSwap───────iSwap──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────PhXZ(a=-0.5,x=0.5,z=0.5)──────────────iSwap────────PhXZ(a=-0.5,x=1,z=0)───────────────────iSwap────────PhXZ(a=-1,x=0.5,z=1.04e-08)───────────────────────────────────────────
                                    │           │                                                                                                                                                                    │                                                   │
1: ─────────────────────────────────┼───────────┼─────────────────────────────────────iSwap───────PhXZ(a=-0.83,x=0.5,z=0.83)───iSwap───────PhXZ(a=0.5,x=0.5,z=0)───iSwap───────PhXZ(a=-0.5,x=0.92,z=0.5)────iSwap────┼────────────PhXZ(a=0.17,x=0.5,z=-0.17)────iSwap────┼────────────PhXZ(a=0.5,x=0.5,z=0)─────────iSwap───────PhXZ(a=0.5,x=0.83,z=-0.5)───
                                    │           │                                     │                                        │                                   │                                        │        │                                          │        │                                          │
2: ─────────────────────────────────┼───────────┼─────────────────────────────────────┼────────────────────────────────────────┼───────────────────────────────────┼───────────PhXZ(a=-1,x=0.5,z=0)─────────┼────────iSwap^0.5──────────────────────────────────┼────────iSwap^0.5────PhXZ(a=-1,x=0.5,z=0.25)───────┼───────────────────────────────────────
                                    │           │                                     │                                        │                                   │                                        │                                                   │                                                   │
3: ─────────────────────────────────iSwap^0.5───iSwap^0.5───PhXZ(a=-1.0,x=0,z=-0.5)───iSwap^0.5───PhXZ(a=-0.83,x=0.5,z=0.83)───iSwap^0.5───PhXZ(a=0.5,x=0.5,z=0)───iSwap^0.5───PhXZ(a=0.5,x=0.83,z=0.5)─────iSwap^0.5─────────────PhXZ(a=0.17,x=0.5,z=-0.17)────iSwap^0.5─────────────PhXZ(a=0.5,x=0.5,z=0)─────────iSwap^0.5───PhXZ(a=0.5,x=0.83,z=-1.0)───
                                                                                                                                                                                                           └──────────────────┘                                └──────────────────┘

Other interesting things in Cirq

Experiments. The cirq.experiments package can perform and plot the results of some basic experiments for understanding how well a system is performing.

qubit = cirq.LineQubit(0)
result = cirq.experiments.single_qubit_state_tomography(
    sampler=cirq.Simulator(),  # In case of Google QCS or other hardware providers, sampler could point at real hardware.
    qubit=qubit,
    circuit=cirq.Circuit(cirq.Z(qubit), cirq.X(qubit)),
    repetitions=1000,
)
result.plot();

png

Testing. The cirq.testing package has useful debugging and testing methods like cirq.testing.assert_implements_consistent_protocols and cirq.testing.assert_allclose_up_to_global_phase.

class InconsistentXGate(cirq.Gate):
    def _num_qubits_(self) -> int:
        return 1

    def _decompose_(self, qubits):
        yield cirq.H(qubits[0])
        yield cirq.Z(qubits[0])
        yield cirq.H(qubits[0])

    def _unitary_(self):
        return np.array([[0, -1j], [1j, 0]])  # Oops! Y instead of X!


# cirq.testing.assert_decompose_is_consistent_with_unitary(InconsistentXGate())

Export. You can export a circuit as Qasm.

"""Export a circuit to Qasm."""
a, b, c = cirq.LineQubit.range(3)
circuit = cirq.Circuit(cirq.H(a), cirq.H(c), cirq.CNOT(a, b), cirq.CCZ(a, b, c))
print(circuit.to_qasm())

// Generated from Cirq v1.4.1

OPENQASM 2.0;
include "qelib1.inc";


// Qubits: [q(0), q(1), q(2)]
qreg q[3];


h q[0];
h q[2];
cx q[0],q[1];

// Gate: CCZ
h q[2];
ccx q[0],q[1],q[2];
h q[2];

You can also export a circuit as QUIL. For that functionality, install the cirq-rigetti package.

You can also turn a circuit into a link to the drag-and-drop web simulation Quirk (though somewhat inconveniently).

"""Export a circuit to a Quirk URL."""
from cirq.contrib.quirk.export_to_quirk import circuit_to_quirk_url

print(circuit_to_quirk_url(circuit))

http://algassert.com/quirk#circuit=%7B%22cols%22%3A%5B%5B%22H%22%2C1%2C%22H%22%5D%2C%5B%22%E2%80%A2%22%2C%22X%22%5D%2C%5B%22%E2%80%A2%22%2C%22%E2%80%A2%22%2C%22Z%22%5D%5D%7D

Introduction to Cirq

Installing Cirq

Qubits, Operations, Moments and Circuits

Create a Circuit

Exercise: Create a Circuit

Unpacking the circuit

More ways to create Circuits

Exercise: Create a circuit to left rotate 5 qubits.

Insert strategies

Exercise: Create the given circuit using least number of appends

Simulations of a Circuit

Repetitions and histograms

The Deutsch-Jozsa Algorithm

Exercise: Two Bit Deutsch-Jozsa Algorithm

Gates

Custom gates

Exercise: Custom Controlled Rx gate

Decompose Protocol

Gatesets

Parameterized Circuits

Exercise: Rotate a qubit

Noise

Monte carlo simulations

Adding noise to circuits and simulations

Devices

Exercise: Make a Device

Quantum Virtual Machine

Transformers for circuit compilation

Exercise: Simplify flipped CNOTs

Compiling to Target Gatesets

Other interesting things in Cirq

Create a `Circuit`

More ways to create `Circuit`s