# Source code for pyquil.numpy_simulator

```
##############################################################################
# Copyright 2018 Rigetti Computing
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
##############################################################################
from typing import List, Union, Sequence
import numpy as np
from numpy.random.mtrand import RandomState
from pyquil import Program
from pyquil.gate_matrices import QUANTUM_GATES
from pyquil.paulis import PauliTerm, PauliSum
from pyquil.quilbase import Gate
from pyquil.reference_simulator import AbstractQuantumSimulator
# The following function is lovingly copied from the Cirq project
# https://github.com/quantumlib/Cirq
#
# With the original copyright disclaimer:
# Copyright 2018 The Cirq Developers
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# https://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
from pyquil.unitary_tools import all_bitstrings
[docs]def targeted_einsum(gate: np.ndarray,
wf: np.ndarray,
wf_target_inds: List[int]
) -> np.ndarray:
"""Left-multiplies the given axes of the wf tensor by the given gate matrix.
Note that the matrix must have a compatible tensor structure.
For example, if you have an 6-qubit state vector ``wf`` with shape
(2, 2, 2, 2, 2, 2), and a 2-qubit unitary operation ``op`` with shape
(2, 2, 2, 2), and you want to apply ``op`` to the 5th and 3rd qubits
within ``input_state``, then the output state vector is computed as follows::
output_state = targeted_einsum(op, input_state, [5, 3])
This method also works when the right hand side is a matrix instead of a
vector. If a unitary circuit's matrix is ``old_effect``, and you append
a CNOT(q1, q4) operation onto the circuit, where the control q1 is the qubit
at offset 1 and the target q4 is the qubit at offset 4, then the appended
circuit's unitary matrix is computed as follows::
new_effect = targeted_left_multiply(CNOT.reshape((2, 2, 2, 2)), old_effect, [1, 4])
:param gate: What to left-multiply the target tensor by.
:param wf: A tensor to carefully broadcast a left-multiply over.
:param wf_target_inds: Which axes of the target are being operated on.
:returns: The output tensor.
"""
k = len(wf_target_inds)
d = len(wf.shape)
work_indices = tuple(range(k))
data_indices = tuple(range(k, k + d))
used_data_indices = tuple(data_indices[q] for q in wf_target_inds)
input_indices = work_indices + used_data_indices
output_indices = list(data_indices)
for w, t in zip(work_indices, wf_target_inds):
output_indices[t] = w
# TODO: `out` does not work if input matrices share memory with outputs, as is usually
# TODO: the case when propogating a wavefunction. This might be fixed in numpy 1.15
# https://github.com/numpy/numpy/pull/11286
# It might be worth re-investigating memory savings with `out` when numpy 1.15 becomes
# commonplace.
return np.einsum(gate, input_indices,
wf, data_indices,
output_indices)
[docs]def targeted_tensordot(gate: np.ndarray,
wf: np.ndarray,
wf_target_inds: Sequence[int]
) -> np.ndarray:
"""Left-multiplies the given axes of the wf tensor by the given gate matrix.
Compare with :py:func:`targeted_einsum`. The semantics of these two functions should be
identical, except this uses ``np.tensordot`` instead of ``np.einsum``.
:param gate: What to left-multiply the target tensor by.
:param wf: A tensor to carefully broadcast a left-multiply over.
:param wf_target_inds: Which axes of the target are being operated on.
:returns: The output tensor.
"""
gate_n_qubits = gate.ndim // 2
n_qubits = wf.ndim
# the indices we want to sum over are the final half
gate_inds = np.arange(gate_n_qubits, 2 * gate_n_qubits)
assert len(wf_target_inds) == len(gate_inds), (wf_target_inds, gate_inds)
wf = np.tensordot(gate, wf, (gate_inds, wf_target_inds))
# tensordot dumps "output" indices into 0, 1, .. gate_n_qubits
# we need to move them to the right place.
# First create a list of all the "unaffected" indices which is everything but the
# first `gate_n_qubits`
axes_ordering = list(range(gate_n_qubits, n_qubits))
# We want to "insert" the affected indices into the right place. This means
# we have to be extra careful about calling list.insert in the correct order.
# Namely, we have to insert low target indices first.
where_td_put_them = np.arange(gate_n_qubits)
sorty = np.argsort(wf_target_inds)
where_td_put_them = where_td_put_them[sorty]
sorted_targets = np.asarray(wf_target_inds)[sorty]
# now that everything is sorted, we can do the insertion.
for target_ind, from_ind in zip(sorted_targets, where_td_put_them):
axes_ordering.insert(target_ind, from_ind)
# A quick call to transpose gives us the right thing.
return wf.transpose(axes_ordering)
def get_measure_probabilities(wf, qubit):
"""
Get the probabilities of measuring a qubit.
:param wf: The statevector with a dimension for each qubit
:param qubit: The qubit to measure. We will sum over every axis except this one.
:return: A vector of classical probabilities.
"""
n_qubits = len(wf.shape)
all_inds = list(range(n_qubits))
return np.einsum(np.conj(wf), all_inds,
wf, all_inds,
[int(qubit)])
def _get_gate_tensor_and_qubits(gate: Gate):
"""Given a gate ``Instruction``, turn it into a matrix and extract qubit indices.
:param gate: the instruction
:return: tensor, qubit_inds.
"""
if len(gate.params) > 0:
matrix = QUANTUM_GATES[gate.name](*gate.params)
else:
matrix = QUANTUM_GATES[gate.name]
qubit_inds = [q.index for q in gate.qubits]
# e.g. 2-qubit matrix is 4x4; turns into (2,2,2,2) tensor.
tensor = np.reshape(matrix, (2,) * len(qubit_inds) * 2)
return tensor, qubit_inds
def _term_expectation(wf, term: PauliTerm):
# Computes <psi|XYZ..XXZ|psi>
wf2 = wf
for qubit_i, op_str in term._ops.items():
# Re-use QUANTUM_GATES since it has X, Y, Z
op_mat = QUANTUM_GATES[op_str]
wf2 = targeted_tensordot(gate=op_mat, wf=wf2, wf_target_inds=[qubit_i])
# `wf2` is XYZ..XXZ|psi>
# hit it with a <psi| i.e. `wf.dag`
return term.coefficient * np.tensordot(wf.conj(), wf2, axes=len(wf.shape))
[docs]class NumpyWavefunctionSimulator(AbstractQuantumSimulator):
def __init__(self, n_qubits, rs: RandomState = None):
"""
A wavefunction simulator that uses numpy's tensordot or einsum to update a state vector
Please consider using
:py:class:`PyQVM(..., quantum_simulator_type=NumpyWavefunctionSimulator)` rather
than using this class directly.
This class uses a n_qubit-dim ndarray to store wavefunction
amplitudes. The array is indexed into with a tuple of n_qubits 1's and 0's, with
qubit 0 as the leftmost bit. This is the opposite convention of the Rigetti Lisp QVM.
:param n_qubits: Number of qubits to simulate.
:param rs: a RandomState (should be shared with the owning :py:class:`PyQVM`) for
doing anything stochastic. A value of ``None`` disallows doing anything stochastic.
"""
self.n_qubits = n_qubits
self.rs = rs
self.wf = np.zeros((2,) * n_qubits, dtype=np.complex128)
self.wf[(0,) * n_qubits] = complex(1.0, 0)
def sample_bitstrings(self, n_samples):
"""
Sample bitstrings from the distribution defined by the wavefunction.
Qubit 0 is at ``out[:, 0]``.
:param n_samples: The number of bitstrings to sample
:return: An array of shape (n_samples, n_qubits)
"""
if self.rs is None:
raise ValueError("You have tried to perform a stochastic operation without setting the "
"random state of the simulator. Might I suggest using a PyQVM object?")
# note on reshape: it puts bitstrings in lexicographical order.
# would you look at that .. _all_bitstrings returns things in lexicographical order!
# reminder: qubit 0 is on the left in einsum simulator.
probabilities = np.abs(self.wf.reshape(-1)) ** 2
possible_bitstrings = all_bitstrings(self.n_qubits)
inds = self.rs.choice(2 ** self.n_qubits, n_samples, p=probabilities)
return possible_bitstrings[inds, :]
def do_measurement(self, qubit: int) -> int:
"""
Measure a qubit, collapse the wavefunction, and return the measurement result.
:param qubit: Index of the qubit to measure.
:return: measured bit
"""
if self.rs is None:
raise ValueError("You have tried to perform a stochastic operation without setting the "
"random state of the simulator. Might I suggest using a PyQVM object?")
# Get probabilities
measurement_probs = get_measure_probabilities(self.wf, qubit)
# Flip a coin and record the result
measured_bit = int(np.argmax(self.rs.uniform() < np.cumsum(measurement_probs)))
# Zero out amplitudes corresponding to non-measured bistrings
other_bit = (measured_bit + 1) % 2
other_bit_indices = (slice(None),) * qubit + \
(other_bit,) + \
(slice(None),) * (self.n_qubits - qubit - 1)
self.wf[other_bit_indices] = 0
# Re-normalize amplitudes corresponding to measured bistrings
meas_bit_indices = (slice(None),) * qubit + \
(measured_bit,) + \
(slice(None),) * (self.n_qubits - qubit - 1)
self.wf[meas_bit_indices] /= np.sqrt(measurement_probs[measured_bit])
return measured_bit
def do_gate(self, gate: Gate):
"""
Perform a gate.
:return: ``self`` to support method chaining.
"""
gate_matrix, qubit_inds = _get_gate_tensor_and_qubits(gate=gate)
# Note to developers: you can use either einsum- or tensordot- based functions.
# tensordot seems a little faster, but feel free to experiment.
# self.wf = targeted_einsum(gate=gate_matrix, wf=self.wf, wf_target_inds=qubit_inds)
self.wf = targeted_tensordot(gate=gate_matrix, wf=self.wf, wf_target_inds=qubit_inds)
return self
def do_gate_matrix(self, matrix: np.ndarray,
qubits: Sequence[int]) -> 'AbstractQuantumSimulator':
"""
Apply an arbitrary unitary; not necessarily a named gate.
:param matrix: The unitary matrix to apply. No checks are done
:param qubits: A list of qubits to apply the unitary to.
:return: ``self`` to support method chaining.
"""
# e.g. 2-qubit matrix is 4x4; turns into (2,2,2,2) tensor.
tensor = np.reshape(matrix, (2,) * len(qubits) * 2)
# Note to developers: you can use either einsum- or tensordot- based functions.
# tensordot seems a little faster, but feel free to experiment.
# self.wf = targeted_einsum(gate=gate_matrix, wf=self.wf, wf_target_inds=qubits)
self.wf = targeted_tensordot(gate=tensor, wf=self.wf, wf_target_inds=qubits)
return self
def expectation(self, operator: Union[PauliTerm, PauliSum]):
"""
Compute the expectation of an operator.
:param operator: The operator
:return: The operator's expectation value
"""
if not isinstance(operator, PauliSum):
operator = PauliSum([operator])
return sum(_term_expectation(self.wf, term) for term in operator)
def reset(self):
"""
Reset the wavefunction to the |000...00> state.
:return: ``self`` to support method chaining.
"""
self.wf.fill(0)
self.wf[(0,) * self.n_qubits] = complex(1.0, 0)
return self
def do_post_gate_noise(self, noise_type: str, noise_prob: float,
qubits: List[int]) -> 'AbstractQuantumSimulator':
raise NotImplementedError("The numpy simulator cannot handle noise")
```