Source code for qutip_qip.circuit.circuitsimulator

from itertools import product, chain
import os

import numpy as np

from . import circuit_latex as _latex
from ..operations import (
    Gate,
    Measurement,
    expand_operator,
    gate_sequence_product,
)
from qutip import basis, ket2dm, Qobj, tensor


__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, N, [ind_map[ind] for ind in inds_sublist]
    )
    U = expand_operator(U, N, [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, len(overall_inds), 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, num_qubits, 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.

    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, state=None, cbits=None, U_list=None, measure_results=None, mode="state_vector_simulator", precompute_unitary=False, ): """ Simulate state evolution for Quantum Circuits. Parameters ---------- qc : :class:`.QubitCircuit` Quantum Circuit to be simulated. state: ket or oper ket or density matrix cbits: list of int, optional initial value of classical bits U_list: list of Qobj, optional list of predefined unitaries corresponding to circuit. 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. 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. precompute_unitary: Boolean, optional Specify if computation is done by pre-computing and aggregating gate unitaries. Possibly a faster method in the case of large number of repeat runs with different state inputs. """ self.qc = qc self.mode = mode self.precompute_unitary = precompute_unitary if U_list: self.U_list = U_list elif precompute_unitary: self.U_list = qc.propagators(expand=False, ignore_measurement=True) else: self.U_list = qc.propagators(ignore_measurement=True) self.ops = [] self.inds_list = [] if precompute_unitary: self._process_ops_precompute() else: self._process_ops() self.initialize(state, cbits, measure_results) def _process_ops(self): """ Process list of gates (including measurements), and stores them in self.ops (as unitaries) for further computation. """ U_list_index = 0 for operation in self.qc.gates: if isinstance(operation, Measurement): self.ops.append(operation) elif isinstance(operation, Gate): if operation.classical_controls: self.ops.append((operation, self.U_list[U_list_index])) else: self.ops.append(self.U_list[U_list_index]) U_list_index += 1 def _process_ops_precompute(self): """ Process list of gates (including measurements), aggregate gate unitaries (by multiplying) and store them in self.ops for further computation. The gate multiplication is carried out only for groups of matrices in between classically controlled gates and measurement gates. Examples -------- If we have a circuit that looks like: ----|X|-----|Y|----|M0|-----|X|---- then self.ops = [YX, M0, X] """ prev_index = 0 U_list_index = 0 for gate in self.qc.gates: if isinstance(gate, Measurement): continue else: self.inds_list.append(gate.get_all_qubits()) for operation in self.qc.gates: if isinstance(operation, Measurement): if U_list_index > prev_index: self.ops.append( self._compute_unitary( self.U_list[prev_index:U_list_index], self.inds_list[prev_index:U_list_index], ) ) prev_index = U_list_index self.ops.append(operation) elif isinstance(operation, Gate): if operation.classical_controls: if U_list_index > prev_index: self.ops.append( self._compute_unitary( self.U_list[prev_index:U_list_index], self.inds_list[prev_index:U_list_index], ) ) prev_index = U_list_index self.ops.append((operation, self.U_list[prev_index])) prev_index += 1 U_list_index += 1 else: U_list_index += 1 if U_list_index > prev_index: self.ops.append( self._compute_unitary( self.U_list[prev_index:U_list_index], self.inds_list[prev_index:U_list_index], ) ) prev_index = U_list_index + 1 U_list_index = prev_index
[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 U_list: list of Qobj, optional list of predefined unitaries corresponding to circuit. 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. """ 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 self.state = None if state is not None: if self.mode == "density_matrix_simulator" and state.isket: self.state = ket2dm(state) else: self.state = state self.probability = 1 self.op_index = 0 self.measure_results = measure_results self.measure_ind = 0
def _compute_unitary(self, U_list, inds_list): """ Compute unitary corresponding to a product of unitaries in U_list and expand it to size of circuit. Parameters ---------- U_list: list of Qobj list of predefined unitaries. inds_list: list of list of int list of qubit indices corresponding to each unitary in U_list Returns ------- U: Qobj resultant unitary """ U_overall, overall_inds = gate_sequence_product( U_list, inds_list=inds_list, expand=True ) if len(overall_inds) != self.qc.N: U_overall = expand_operator( U_overall, N=self.qc.N, targets=overall_inds ) return U_overall
[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.ops)): if self.step() is None: 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. """ op = self.ops[self.op_index] if isinstance(op, Measurement): self._apply_measurement(op) elif isinstance(op, tuple): operation, U = op apply_gate = all( [self.cbits[i] for i in operation.classical_controls] ) if apply_gate: if self.precompute_unitary: U = expand_operator( U, self.qc.N, operation.get_all_qubits() ) self._evolve_state(U) else: self._evolve_state(op) self.op_index += 1 return self.state
def _evolve_state(self, U): """ Applies unitary to state. Parameters ---------- U: Qobj unitary to be applied. """ if self.mode == "state_vector_simulator": self._evolve_ket(U) elif self.mode == "density_matrix_simulator": self._evolve_dm(U) else: raise NotImplementedError( "mode {} is not available.".format(self.mode) ) def _evolve_ket(self, U): """ Applies unitary to ket state. Parameters ---------- U: Qobj unitary to be applied. """ self.state = U * self.state def _evolve_dm(self, U): """ Applies unitary to density matrix state. Parameters ---------- U: Qobj unitary to be applied. """ self.state = U * self.state * U.dag() def _apply_measurement(self, operation): """ 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] self.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)) self.state = sum(p * s for s, p in zip(states, probabilities)) else: raise NotImplementedError( "mode {} is not available.".format(self.mode) )
[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