from itertools import product, chain
from operator import mul
from functools import reduce
import numpy as np
from ..operations import (
Gate,
Measurement,
expand_operator,
)
from qutip import basis, ket2dm, Qobj, tensor
import warnings
__all__ = ["CircuitSimulator", "CircuitResult"]
def _flatten(lst):
"""
Helper to flatten lists.
"""
return [item for sublist in lst for item in sublist]
def _mult_sublists(tensor_list, overall_inds, U, inds):
"""
Calculate the revised indices and tensor list by multiplying a new unitary
U applied to inds.
Parameters
----------
tensor_list : list of Qobj
List of gates (unitaries) acting on disjoint qubits.
overall_inds : list of list of int
List of qubit indices corresponding to each gate in tensor_list.
U: Qobj
Unitary to be multiplied with the the unitary specified by tensor_list.
inds: list of int
List of qubit indices corresponding to U.
Returns
-------
tensor_list_revised: list of Qobj
List of gates (unitaries) acting on disjoint qubits incorporating U.
overall_inds_revised: list of list of int
List of qubit indices corresponding to each gate in tensor_list_revised.
Examples
--------
First, we get some imports out of the way,
>>> from qutip_qip.operations.gates import _mult_sublists
>>> from qutip_qip.operations.gates import x_gate, y_gate, toffoli, z_gate
Suppose we have a unitary list of already processed gates,
X, Y, Z applied on qubit indices 0, 1, 2 respectively and
encounter a new TOFFOLI gate on qubit indices (0, 1, 3).
>>> tensor_list = [x_gate(), y_gate(), z_gate()]
>>> overall_inds = [[0], [1], [2]]
>>> U = toffoli()
>>> U_inds = [0, 1, 3]
Then, we can use _mult_sublists to produce a new list of unitaries by
multiplying TOFFOLI (and expanding) only on the qubit indices involving
TOFFOLI gate (and any multiplied gates).
>>> U_list, overall_inds = _mult_sublists(tensor_list, overall_inds, U, U_inds)
>>> np.testing.assert_allclose(U_list[0]) == z_gate())
>>> toffoli_xy = toffoli() * tensor(x_gate(), y_gate(), identity(2))
>>> np.testing.assert_allclose(U_list[1]), toffoli_xy)
>>> overall_inds = [[2], [0, 1, 3]]
"""
tensor_sublist = []
inds_sublist = []
tensor_list_revised = []
overall_inds_revised = []
for sub_inds, sub_U in zip(overall_inds, tensor_list):
if len(set(sub_inds).intersection(inds)) > 0:
tensor_sublist.append(sub_U)
inds_sublist.append(sub_inds)
else:
overall_inds_revised.append(sub_inds)
tensor_list_revised.append(sub_U)
inds_sublist = _flatten(inds_sublist)
U_sublist = tensor(tensor_sublist)
revised_inds = list(set(inds_sublist).union(set(inds)))
N = len(revised_inds)
sorted_positions = sorted(range(N), key=lambda key: revised_inds[key])
ind_map = {ind: pos for ind, pos in zip(revised_inds, sorted_positions)}
U_sublist = expand_operator(
U_sublist, dims=[2] * N, targets=[ind_map[ind] for ind in inds_sublist]
)
U = expand_operator(
U, dims=[2] * N, targets=[ind_map[ind] for ind in inds]
)
U_sublist = U * U_sublist
inds_sublist = revised_inds
overall_inds_revised.append(inds_sublist)
tensor_list_revised.append(U_sublist)
return tensor_list_revised, overall_inds_revised
def _expand_overall(tensor_list, overall_inds):
"""
Tensor unitaries in tensor list and then use expand_operator to rearrange
them appropriately according to the indices in overall_inds.
"""
U_overall = tensor(tensor_list)
overall_inds = _flatten(overall_inds)
U_overall = expand_operator(
U_overall, dims=[2] * len(overall_inds), targets=overall_inds
)
overall_inds = sorted(overall_inds)
return U_overall, overall_inds
def _gate_sequence_product(U_list, ind_list):
"""
Calculate the overall unitary matrix for a given list of unitary operations
that are still of original dimension.
Parameters
----------
U_list : list of Qobj
List of gates(unitaries) implementing the quantum circuit.
ind_list : list of list of int
List of qubit indices corresponding to each gate in tensor_list.
Returns
-------
U_overall : qobj
Unitary matrix corresponding to U_list.
overall_inds : list of int
List of qubit indices on which U_overall applies.
Examples
--------
First, we get some imports out of the way,
>>> from qutip_qip.operations.gates import _gate_sequence_product
>>> from qutip_qip.operations.gates import x_gate, y_gate, toffoli, z_gate
Suppose we have a circuit with gates X, Y, Z, TOFFOLI
applied on qubit indices 0, 1, 2 and [0, 1, 3] respectively.
>>> tensor_lst = [x_gate(), y_gate(), z_gate(), toffoli()]
>>> overall_inds = [[0], [1], [2], [0, 1, 3]]
Then, we can use _gate_sequence_product to produce a single unitary
obtained by multiplying unitaries in the list using heuristic methods
to reduce the size of matrices being multiplied.
>>> U_list, overall_inds = _gate_sequence_product(tensor_lst, overall_inds)
"""
num_qubits = len(set(chain(*ind_list)))
sorted_inds = sorted(set(_flatten(ind_list)))
ind_list = [[sorted_inds.index(ind) for ind in inds] for inds in ind_list]
U_overall = 1
overall_inds = []
for i, (U, inds) in enumerate(zip(U_list, ind_list)):
# when the tensor_list covers the full dimension of the circuit, we
# expand the tensor_list to a unitary and call _gate_sequence_product
# recursively on the rest of the U_list.
if len(overall_inds) == 1 and len(overall_inds[0]) == num_qubits:
U_overall, overall_inds = _expand_overall(
tensor_list, overall_inds
)
U_left, rem_inds = _gate_sequence_product(U_list[i:], ind_list[i:])
U_left = expand_operator(
U_left, dims=[2] * num_qubits, targets=rem_inds
)
return U_left * U_overall, [
sorted_inds[ind] for ind in overall_inds
]
# special case for first unitary in the list
if U_overall == 1:
U_overall = U_overall * U
overall_inds = [ind_list[0]]
tensor_list = [U_overall]
continue
# case where the next unitary interacts on some subset of qubits
# with the unitaries already in tensor_list.
elif len(set(_flatten(overall_inds)).intersection(set(inds))) > 0:
tensor_list, overall_inds = _mult_sublists(
tensor_list, overall_inds, U, inds
)
# case where the next unitary does not interact with any unitary in
# tensor_list
else:
overall_inds.append(inds)
tensor_list.append(U)
U_overall, overall_inds = _expand_overall(tensor_list, overall_inds)
return U_overall, [sorted_inds[ind] for ind in overall_inds]
def _gate_sequence_product_with_expansion(U_list, left_to_right=True):
"""
Calculate the overall unitary matrix for a given list of unitary
operations, assuming that all operations have the same dimension.
This is only for backward compatibility.
Parameters
----------
U_list : list
List of gates(unitaries) implementing the quantum circuit.
left_to_right : Boolean
Check if multiplication is to be done from left to right.
Returns
-------
U_overall : qobj
Unitary matrix corresponding to U_list.
"""
U_overall = 1
for U in U_list:
if left_to_right:
U_overall = U * U_overall
else:
U_overall = U_overall * U
return U_overall
[docs]
class CircuitSimulator:
"""
Operator based circuit simulator.
"""
def __init__(
self,
qc,
mode="state_vector_simulator",
precompute_unitary=False,
):
"""
Simulate state evolution for Quantum Circuits.
Parameters
----------
qc : :class:`.QubitCircuit`
Quantum Circuit to be simulated.
mode: string, optional
Specify if input state (and therefore computation) is in
state-vector mode or in density matrix mode.
In state_vector_simulator mode, the input must be a ket
and with each measurement, one of the collapsed
states is the new state (when using run()).
In density_matrix_simulator mode, the input can be a ket or a
density matrix and after measurement, the new state is the
mixed ensemble state obtained after the measurement.
If in density_matrix_simulator mode and given
a state vector input, the output must be assumed to
be a density matrix.
"""
self._qc = qc
self.dims = qc.dims
self.mode = mode
if precompute_unitary:
warnings.warn(
"Precomputing the full unitary is no longer supported. Switching to normal simulation mode."
)
@property
def qc(self):
return self._qc
[docs]
def initialize(self, state=None, cbits=None, measure_results=None):
"""
Reset Simulator state variables to start a new run.
Parameters
----------
state: ket or oper
ket or density matrix
cbits: list of int, optional
initial value of classical bits
measure_results : tuple of ints, optional
optional specification of each measurement result to enable
post-selection. If specified, the measurement results are
set to the tuple of bits (sequentially) instead of being
chosen at random.
"""
# Initializing the unitary operators.
if cbits and len(cbits) == self.qc.num_cbits:
self.cbits = cbits
elif self.qc.num_cbits > 0:
self.cbits = [0] * self.qc.num_cbits
else:
self.cbits = None
# Parameters that will be updated during the simulation.
# self._state keeps track of the current state of the evolution.
# It is not guaranteed to be a Qobj and could be reshaped.
# Use self.state to return the Qobj representation.
if state is not None:
if self.mode == "density_matrix_simulator" and state.isket:
self._state = ket2dm(state)
else:
self._state = state
else:
# Just computing the full unitary, no state
self._state = None
self._state_dims = (
state.dims.copy()
) # Record the dimension of the state.
self._probability = 1
self._op_index = 0
self._measure_results = measure_results
self._measure_ind = 0
if self.mode == "state_vector_simulator":
self._tensor_dims = self._state_dims[0].copy()
if state.type == "oper":
# apply the gate to a unitary, add an ancillary axis.
self._state_mat_shape = [
reduce(mul, self._state_dims[0], 1)
] * 2
self._tensor_dims += [reduce(mul, self._state_dims[0], 1)]
else:
self._state_mat_shape = [
reduce(mul, self._state_dims[0], 1),
1,
]
self._tensor_dims = tuple(self._tensor_dims)
self._state_mat_shape = tuple(self._state_mat_shape)
@property
def state(self):
"""
The current state of the simulator as a `qutip.Qobj`
Returns:
`qutip.Qobj`: The current state of the simulator.
"""
if not isinstance(self._state, Qobj) and self._state is not None:
self._state = self._state.reshape(self._state_mat_shape)
return Qobj(self._state, dims=self._state_dims)
else:
return self._state
[docs]
def run(self, state, cbits=None, measure_results=None):
"""
Calculate the result of one instance of circuit run.
Parameters
----------
state : ket or oper
state vector or density matrix input.
cbits : List of ints, optional
initialization of the classical bits.
measure_results : tuple of ints, optional
optional specification of each measurement result to enable
post-selection. If specified, the measurement results are
set to the tuple of bits (sequentially) instead of being
chosen at random.
Returns
-------
result: CircuitResult
Return a CircuitResult object containing
output state and probability.
"""
self.initialize(state, cbits, measure_results)
for _ in range(len(self._qc.gates)):
self.step()
if self._state is None:
# TODO This only happens if there is predefined post-selection on the measurement results and the measurement results is exactly 0. This needs to be improved.
break
return CircuitResult(self.state, self._probability, self.cbits)
[docs]
def run_statistics(self, state, cbits=None):
"""
Calculate all the possible outputs of a circuit
(varied by measurement gates).
Parameters
----------
state : ket
state to be observed on specified by density matrix.
cbits : List of ints, optional
initialization of the classical bits.
Returns
-------
result: CircuitResult
Return a CircuitResult object containing
output states and and their probabilities.
"""
probabilities = []
states = []
cbits_results = []
num_measurements = len(
list(filter(lambda x: isinstance(x, Measurement), self._qc.gates))
)
for results in product("01", repeat=num_measurements):
run_result = self.run(state, cbits=cbits, measure_results=results)
final_state = run_result.get_final_states(0)
probability = run_result.get_probabilities(0)
states.append(final_state)
probabilities.append(probability)
cbits_results.append(self.cbits)
return CircuitResult(states, probabilities, cbits_results)
[docs]
def step(self):
"""
Return state after one step of circuit evolution
(gate or measurement).
Returns
-------
state : ket or oper
state after one evolution step.
"""
def _decimal_to_binary(decimal, length):
binary = [int(s) for s in "{0:#b}".format(decimal)[2:]]
return [0] * (length - len(binary)) + binary
def _check_classical_control_value(operation, cbits):
"""Check if the gate should be executed, depending on the current value of classical bits."""
matched = np.empty(len(operation.classical_controls), dtype=bool)
cbits_conditions = _decimal_to_binary(
operation.classical_control_value,
len(operation.classical_controls),
)
for i in range(len(operation.classical_controls)):
cbit_index = operation.classical_controls[i]
control_value = cbits_conditions[i]
matched[i] = cbits[cbit_index] == control_value
return all(matched)
op = self._qc.gates[self._op_index]
self._op_index += 1
current_state = self._state
if isinstance(op, Measurement):
state = self._apply_measurement(op, current_state)
elif isinstance(op, Gate):
if op.classical_controls is not None:
apply_gate = _check_classical_control_value(op, self.cbits)
else:
apply_gate = True
if not apply_gate:
return current_state
if self.mode == "state_vector_simulator":
state = self._evolve_state_einsum(op, current_state)
else:
state = self._evolve_state(op, current_state)
self._state = state
def _evolve_state(self, operation, state):
"""
Applies unitary to state.
Parameters
----------
U: Qobj
unitary to be applied.
"""
if operation.name == "GLOBALPHASE":
# This is just a complex number.
U = np.exp(1.0j * operation.arg_value)
else:
# We need to use the circuit because the custom gates
# are still saved in circuit instance.
# This should be changed once that is deprecated.
U = self.qc._get_gate_unitary(operation)
U = expand_operator(
U,
dims=self.dims,
targets=operation.get_all_qubits(),
)
if self.mode == "state_vector_simulator":
state = U * state
elif self.mode == "density_matrix_simulator":
state = U * state * U.dag()
else:
raise NotImplementedError(
"mode {} is not available.".format(self.mode)
)
return state
def _evolve_state_einsum(self, gate, state):
if gate.name == "GLOBALPHASE":
# This is just a complex number.
return np.exp(1.0j * gate.arg_value) * state
# Prepare the state tensor.
targets_indices = gate.get_all_qubits()
if isinstance(state, Qobj):
# If it is a Qobj, transform it to the array representation.
state = state.full()
# Transform the gate and state array to the corresponding
# tensor form.
state = state.reshape(self._tensor_dims)
# Prepare the gate tensor.
gate = self.qc._get_gate_unitary(gate)
gate_array = gate.full().reshape(gate.dims[0] + gate.dims[1])
# Compute the tensor indices and call einsum.
num_site = len(state.shape)
ancillary_indices = range(num_site, num_site + len(targets_indices))
index_list = range(num_site)
new_index_list = list(index_list)
for j, k in enumerate(targets_indices):
new_index_list[k] = j + num_site
state = np.einsum(
gate_array,
list(ancillary_indices) + list(targets_indices),
state,
index_list,
new_index_list,
)
return state
def _apply_measurement(self, operation, state):
"""
Applies measurement gate specified by operation to current state.
Parameters
----------
operation: :class:`.Measurement`
Measurement gate in a circuit object.
"""
states, probabilities = operation.measurement_comp_basis(self.state)
if self.mode == "state_vector_simulator":
if self._measure_results:
i = int(self._measure_results[self._measure_ind])
self._measure_ind += 1
else:
probabilities = [p / sum(probabilities) for p in probabilities]
i = np.random.choice([0, 1], p=probabilities)
self._probability *= probabilities[i]
state = states[i]
if operation.classical_store is not None:
self.cbits[operation.classical_store] = i
elif self.mode == "density_matrix_simulator":
states = list(filter(lambda x: x is not None, states))
probabilities = list(filter(lambda x: x != 0, probabilities))
state = sum(p * s for s, p in zip(states, probabilities))
else:
raise NotImplementedError(
"mode {} is not available.".format(self.mode)
)
return state
[docs]
class CircuitResult:
"""
Result of a quantum circuit simulation.
"""
def __init__(self, final_states, probabilities, cbits=None):
"""
Store result of CircuitSimulator.
Parameters
----------
final_states: list of Qobj.
List of output kets or density matrices.
probabilities: list of float.
List of probabilities of obtaining each output state.
cbits: list of list of int, optional
List of cbits for each output.
"""
if isinstance(final_states, Qobj) or final_states is None:
self.final_states = [final_states]
self.probabilities = [probabilities]
if cbits:
self.cbits = [cbits]
else:
inds = list(
filter(
lambda x: final_states[x] is not None,
range(len(final_states)),
)
)
self.final_states = [final_states[i] for i in inds]
self.probabilities = [probabilities[i] for i in inds]
if cbits:
self.cbits = [cbits[i] for i in inds]
[docs]
def get_final_states(self, index=None):
"""
Return list of output states.
Parameters
----------
index: int
Indicates i-th state to be returned.
Returns
-------
final_states: Qobj or list of Qobj.
List of output kets or density matrices.
"""
if index is not None:
return self.final_states[index]
return self.final_states
[docs]
def get_probabilities(self, index=None):
"""
Return list of probabilities corresponding to the output states.
Parameters
----------
index: int
Indicates i-th probability to be returned.
Returns
-------
probabilities: float or list of float
Probabilities associated with each output state.
"""
if index is not None:
return self.probabilities[index]
return self.probabilities
[docs]
def get_cbits(self, index=None):
"""
Return list of classical bit outputs corresponding to the results.
Parameters
----------
index: int
Indicates i-th output, probability pair to be returned.
Returns
-------
cbits: list of int or list of list of int
list of classical bit outputs
"""
if index is not None:
return self.cbits[index]
return self.cbits