View on QuantumAI | Run in Google Colab | View source on GitHub | Download notebook |

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

Cirq comes with built-in Python simulators for testing
small circuits. The two main types of simulations that Cirq
supports are pure state and mixed state. The pure state simulations
are supported by `cirq.Simulator`

and the mixed state
simulators are supported by `cirq.DensityMatrixSimulator`

.

The names *pure state simulator* and *mixed state
simulators* refer to the fact that these simulations are
for quantum circuits; including unitary, measurements, and noise
that keeps the evolution in a pure state or a mixed state.
In other words, there are some noisy evolutions
that are supported by the pure state simulator as long as they
preserve the purity of the state.

Some external high-performance simulators also provide an interface to Cirq. These can often provide results faster than Cirq's built-in simulators, especially when working with larger circuits. For details on these tools, see the external simulators section.

## Introduction to pure state simulation

Here is a simple circuit:

```
import cirq
q0 = cirq.GridQubit(0, 0)
q1 = cirq.GridQubit(1, 0)
def basic_circuit(meas=True):
sqrt_x = cirq.X**0.5
yield sqrt_x(q0), sqrt_x(q1)
yield cirq.CZ(q0, q1)
yield sqrt_x(q0), sqrt_x(q1)
if meas:
yield cirq.measure(q0, key='q0'), cirq.measure(q1, key='q1')
circuit = cirq.Circuit()
circuit.append(basic_circuit())
print(circuit)
```

(0, 0): ───X^0.5───@───X^0.5───M('q0')─── │ (1, 0): ───X^0.5───@───X^0.5───M('q1')───

We can simulate this by creating a `<a href="https://quantumai.google/reference/python/cirq/sim/Simulator"><code>cirq.Simulator</code></a>`

and
passing the circuit into its `run`

method:

```
from cirq import Simulator
simulator = Simulator()
result = simulator.run(circuit)
print(result)
```

q0=0 q1=1

The method `run()`

returns a
`Result`

.
As you can see, the object `result`

contains the result of any
measurements for the simulation run.

The actual measurement results depend on the seeding from
`random`

seed generator (`numpy`

). You can set this using
`numpy.random_seed`

.

Another run can result in a different set of measurement results:

```
result = simulator.run(circuit)
print(result)
```

q0=1 q1=0

The simulator is designed to mimic what running a program on a quantum computer is actually like.

### Accessing the state vector

In particular, the `run()`

methods (`run()`

and `run_sweep()`

)
on the simulator do not give access to the state vector of the
quantum computer (since one doesn't have access to this on the
actual quantum hardware). Instead, the `simulate()`

methods
(`simulate()`

, `simulate_sweep()`

, `simulate_moment_steps()`

)
should be used if one wants to debug the circuit and get access
to the full state vector:

```
import numpy as np
circuit = cirq.Circuit()
circuit.append(basic_circuit(False))
result = simulator.simulate(circuit, qubit_order=[q0, q1])
print(np.around(result.final_state_vector, 3))
```

[0.5+0.j 0. +0.5j 0. +0.5j 0.5+0.j ]

`simulate()`

returns a `SimulationTrialResult`

containing
the final state, as seen above. The built-in Cirq simulator returns a
`StateVectorTrialResult`

,
which includes a number of utilities for analyzing the final state
vector.

Note that the simulator uses numpy's `float32`

precision
(which is `complex64`

for complex numbers) by default,
but that the simulator can take in a `dtype`

of `np.complex128`

if higher precision is needed.

### Expectation values

For applications that measure expectation values of observables,
the `simulate_expectation_values()`

method provides a simple
interface for returning just the desired expectation values.
This can be more efficient than returning the entire state vector,
particularly when handling multiple results at once, or when using
an external simulator.

```
XX_obs = cirq.X(q0) * cirq.X(q1)
ZZ_obs = cirq.Z(q0) * cirq.Z(q1)
ev_list = simulator.simulate_expectation_values(
cirq.Circuit(basic_circuit(False)),
observables=[XX_obs, ZZ_obs],
)
print(ev_list)
```

[(1+0j), 0j]

`simulate_expectation_values()`

returns a list of
expectation values, one for each observable provided.

## Qubit and Amplitude Ordering

The `qubit_order`

argument to the simulator's `run()`

method
determines the ordering of some results, such as the
amplitudes in the final wave function. The `qubit_order`

argument is optional: when it is omitted, qubits are ordered
ascending by their name (i.e., what their `__str__`

method returns).

The simplest `qubit_order`

value you can provide is a list of
the qubits in the desired ordered. 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. 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`

:

```
outside = [1, 10]
inside = [1, 2]
print(np.kron(outside, inside))
```

[ 1 2 10 20]

More concretely, `k`

'th amplitude in the wave function
will correspond to the `k`

'th case that would be encountered
when nesting loops over the possible values of each qubit.

The first qubit's computational basis values are looped over in the outermost loop, the last qubit's computational basis values are looped over in the inner-most loop, etc.:

```
i = 0
for first in [0, 1]:
for second in [0, 1]:
print('amps[{}] is for first={}, second={}'.format(i, first, second))
i += 1
```

amps[0] is for first=0, second=0 amps[1] is for first=0, second=1 amps[2] is for first=1, second=0 amps[3] is for first=1, second=1

We can check that this is in fact the ordering with a circuit that flips one qubit out of two:

```
q_stay = cirq.NamedQubit('q_stay')
q_flip = cirq.NamedQubit('q_flip')
c = cirq.Circuit(cirq.X(q_flip))
# first qubit in order flipped
result = simulator.simulate(c, qubit_order=[q_flip, q_stay])
print(abs(result.final_state_vector).round(3))
```

[0. 0. 1. 0.]

```
# second qubit in order flipped
result = simulator.simulate(c, qubit_order=[q_stay, q_flip])
print(abs(result.final_state_vector).round(3))
```

[0. 1. 0. 0.]

## Stepping through a circuit

When debugging, it is useful to not just see the end result of a circuit, but to inspect, or even modify, the state of the system at different steps in the circuit.

To support this, Cirq provides a method to return an iterator
over a `Moment`

by `Moment`

simulation. This method is named `simulate_moment_steps`

:

```
circuit = cirq.Circuit()
circuit.append(basic_circuit())
for i, step in enumerate(simulator.simulate_moment_steps(circuit)):
print('state at step %d: %s' % (i, np.around(step.state_vector(), 3)))
```

state at step 0: [0. +0.5j 0.5+0.j 0.5+0.j 0. -0.5j] state at step 1: [0. +0.5j 0.5+0.j 0.5+0.j 0. +0.5j] state at step 2: [0.5+0.j 0. +0.5j 0. +0.5j 0.5+0.j ] state at step 3: [1.+0.j 0.+0.j 0.+0.j 0.+0.j]

The object returned by the `moment_steps`

iterator is a
`StepResult`

. This object has the state along with any
measurements that occurred **during** that step (so does
not include measurement results from previous `Moments`

).

In addition, the `StepResult`

contains `set_state_vector()`

,
which can be used to set the `state`

. One can pass a valid
full state to this method by passing a numpy array. Or,
alternatively, one can pass an integer, and then the state
will be set to lie entirely in the computation basis state
for the binary expansion of the passed integer.

### Alternate stepping behavior

For simulators that do not support `simulate_moment_steps`

,
it is possible to replicate this behavior by splitting
the circuit into "chunks" and passing the results of each
chunk as the initial state for the next chunk:

```
chunks = [cirq.Circuit(moment) for moment in basic_circuit()]
next_state = 0 # represents the all-zero state
for i, chunk in enumerate(chunks):
result = simulator.simulate(chunk, initial_state=next_state)
next_state = result.final_state_vector
print(f'state at step {i}: {np.around(next_state, 3)}')
```

state at step 0: [0. +0.5j 0.5+0.j 0.5+0.j 0. -0.5j] state at step 1: [0. +0.5j 0.5+0.j 0.5+0.j 0. +0.5j] state at step 2: [0.5+0.j 0. +0.5j 0. +0.5j 0.5+0.j ] state at step 3: [1.+0.j 0.+0.j 0.+0.j 0.+0.j]

The added cost of passing state vectors around like this
is nontrivial; for this reason, this workaround should
only be used with simulators that do not support
`simulate_moment_steps`

.

## Parameterized values and studies

In addition to circuit gates with fixed values, Cirq also
supports gates which can have `Symbol`

value (see
Gates). These are values that can be resolved
at *run-time*.

For simulators, these values are resolved by
providing a `cirq.ParamResolver`

. A `cirq.ParamResolver`

provides
a map from the `Symbol`

's name to its assigned value.

```
import sympy
rot_w_gate = cirq.X**sympy.Symbol('x')
circuit = cirq.Circuit()
circuit.append([rot_w_gate(q0), rot_w_gate(q1)])
for y in range(5):
resolver = cirq.ParamResolver({'x': y / 4.0})
result = simulator.simulate(circuit, resolver)
print(np.round(result.final_state_vector, 2))
```

[1.+0.j 0.+0.j 0.+0.j 0.+0.j] [ 0.6 +0.6j 0.25-0.25j 0.25-0.25j -0.1 -0.1j ] [0. +0.5j 0.5+0.j 0.5+0.j 0. -0.5j] [-0.1 +0.1j 0.25+0.25j 0.25+0.25j 0.6 -0.6j ] [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
`sweep`

. A `sweep`

is a sequence of trials, where each
trial is a run with a particular set of parameter values.

Running a `sweep`

returns a `Result`

for each set of fixed parameter values and repetitions.

For instance:

```
resolvers = [cirq.ParamResolver({'x': y / 2.0}) for y in range(3)]
circuit = cirq.Circuit()
circuit.append([rot_w_gate(q0), rot_w_gate(q1)])
circuit.append([cirq.measure(q0, key='q0'), cirq.measure(q1, key='q1')])
results = simulator.run_sweep(program=circuit,
params=resolvers,
repetitions=2)
for result in results:
print(result)
```

q0=00 q1=00 q0=10 q1=10 q0=11 q1=11

Above we see that assigning different values to gate parameters yields
different results for each trial in the sweep, and that each trial is repeated
`repetitions`

times.

## Mixed state simulations

In addition to pure state simulation, Cirq also supports simulation of mixed states.

Even though this simulator is not as efficient as the pure state simulators, they allow for a larger class of noisy circuits to be run as well as keeping track of the simulation's density matrix. This fact can allow for more exact simulations (for example, the pure state simulator's Monte Carlo simulation only allows sampling from the density matrix, not explicitly giving the entries of the density matrix like the mixed state simulator can do).

Mixed state simulation is supported by the
`cirq.DensityMatrixSimulator`

class.

Here is a simple example of simulating a channel using the mixed state simulator:

```
q = cirq.NamedQubit('a')
circuit = cirq.Circuit(cirq.H(q), cirq.amplitude_damp(0.2)(q), cirq.measure(q))
simulator = cirq.DensityMatrixSimulator()
result = simulator.run(circuit, repetitions=100)
print(result.histogram(key='a'))
```

Counter({0: 59, 1: 41})

We create a state in an equal superposition of 0 and 1, then apply amplitude damping which takes 1 to 0 with something like a probability of 0.2.

We see that instead of about 50 percent of the timing being in 0, about 20 percent of the 1 has been converted into 0, so we end up with total around 60 percent in the 0 state.

Like the pure state simulators, the mixed state simulator
supports `run()`

and `run_sweeps()`

methods.

The `cirq.DensityMatrixSimulator`

also supports getting access
to the density matrix of the circuit at the end of simulating
the circuit, or when stepping through the circuit. These are
done by the `simulate()`

and `simulate_sweep()`

methods, or,
for stepping through the circuit, via the `simulate_moment_steps`

method. For example, we can simulate creating an equal
superposition followed by an amplitude damping channel with a
gamma of 0.2 by:

```
q = cirq.NamedQubit('a')
circuit = cirq.Circuit(cirq.H(q), cirq.amplitude_damp(0.2)(q))
simulator = cirq.DensityMatrixSimulator()
result = simulator.simulate(circuit)
print(np.around(result.final_density_matrix, 3))
```

[[0.6 +0.j 0.447+0.j] [0.447+0.j 0.4 +0.j]]

We see that we have access to the density matrix at the
end of the simulation via `final_density_matrix`

.

## External simulators

There are a few high-performance circuit simulators which
provide an interface for simulating Cirq `Circuit`

s.
These projects are listed below, along with their PyPI package
name and a description of simulator methods that they support.

Project name | PyPI package | Description |
---|---|---|

qsim | qsimcirq | Implements `cirq.SimulatesAmplitudes` and `cirq.SimulatesFinalState` . Recommended for deep circuits with up to 30 qubits (consumes 8GB RAM). Larger circuits are possible, but RAM usage doubles with each additional qubit. |

qsimh | qsimcirq | Implements `cirq.SimulatesAmplitudes` . Intended for heavy parallelization across several computers; Cirq users should generally prefer qsim. |

qFlex | qflexcirq | Implements `cirq.SimulatesAmplitudes` . Recommended for shallow / low-entanglement circuits with more than 30 qubits. RAM usage is highly dependent on the number of two-qubit gates in the circuit. |

quimb | quimb | Cirq-to-quimb translation layer provided in `contrib/quimb` . In addition to circuit simulation, this allows the use of quimb circuit-analysis tools on Cirq circuits. |