diff --git a/doc/source/qip-simulator.rst b/doc/source/qip-simulator.rst index 553652a3..c0f1c767 100644 --- a/doc/source/qip-simulator.rst +++ b/doc/source/qip-simulator.rst @@ -169,7 +169,8 @@ The :class:`.CircuitSimulator` class also enables stepping through the circuit: .. testcode:: - print(sim.step()) + sim.step() + print(sim.state) **Output**: @@ -191,44 +192,6 @@ This only executes one gate in the circuit and allows for a better understanding of how the state evolution takes place. The method steps through both the gates and the measurements. -Precomputing the unitary -======================== - -By default, the :class:`.CircuitSimulator` class is initialized such that -the circuit evolution is conducted by applying each unitary to the state interactively. -However, by setting the argument ``precompute_unitary=True``, :class:`.CircuitSimulator` -precomputes the product of the unitaries (in between the measurements): - -.. testcode:: - - sim = CircuitSimulator(qc, precompute_unitary=True) - - print(sim.ops) - -.. testoutput:: - :options: +NORMALIZE_WHITESPACE - - [Quantum object: dims = [[2, 2, 2], [2, 2, 2]], shape = (8, 8), type = oper, isherm = False - Qobj data = - [[ 0. 0.57735 0. -0.57735 0. 0.40825 0. -0.40825] - [ 0.57735 0. -0.57735 0. 0.40825 0. -0.40825 0. ] - [ 0.57735 0. 0.57735 0. 0.40825 0. 0.40825 0. ] - [ 0. 0.57735 0. 0.57735 0. 0.40825 0. 0.40825] - [ 0.57735 0. 0. 0. -0.8165 0. 0. 0. ] - [ 0. 0.57735 0. 0. 0. -0.8165 0. 0. ] - [ 0. 0. 0.57735 0. 0. 0. -0.8165 0. ] - [ 0. 0. 0. 0.57735 0. 0. 0. -0.8165 ]], - Measurement(M0, target=[0], classical_store=0), - Measurement(M1, target=[1], classical_store=1), - Measurement(M2, target=[2], classical_store=2)] - - -Here, ``sim.ops`` stores all the circuit operations that are going to be applied during -state evolution. As observed above, all the unitaries of the circuit are compressed into -a single unitary product with the precompute optimization enabled. -This is more efficient if one runs the same circuit one multiple initial states. -However, as the number of qubits increases, this will consume more and more memory -and become unfeasible. Density Matrix Simulation ========================= diff --git a/src/qutip_qip/circuit/circuit.py b/src/qutip_qip/circuit/circuit.py index d74ad665..c3f8d52e 100644 --- a/src/qutip_qip/circuit/circuit.py +++ b/src/qutip_qip/circuit/circuit.py @@ -18,13 +18,12 @@ Measurement, expand_operator, GATE_CLASS_MAP, - gate_sequence_product, ) from .circuitsimulator import ( CircuitSimulator, CircuitResult, ) -from qutip import basis, Qobj +from qutip import Qobj, qeye try: @@ -481,10 +480,6 @@ def run( post-selection. If specified, the measurement results are set to the tuple of bits (sequentially) instead of being chosen at random. - 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. Returns ------- @@ -499,7 +494,6 @@ def run( raise TypeError("State is not a ket or a density matrix.") sim = CircuitSimulator( self, - U_list, mode, precompute_unitary, ) @@ -520,10 +514,6 @@ def run_statistics( initialization of the classical bits. U_list: list of Qobj, optional list of predefined unitaries corresponding to circuit. - 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. Returns ------- @@ -537,12 +527,7 @@ def run_statistics( mode = "density_matrix_simulator" else: raise TypeError("State is not a ket or a density matrix.") - sim = CircuitSimulator( - self, - U_list, - mode, - precompute_unitary, - ) + sim = CircuitSimulator(self, mode, precompute_unitary) return sim.run_statistics(state, cbits) def resolve_gates(self, basis=["CNOT", "RX", "RY", "RZ"]): @@ -892,48 +877,46 @@ def propagators(self, expand=True, ignore_measurement=False): "Cannot compute the propagator of a measurement operator." "Please set ignore_measurement=True." ) - for gate in gates: if gate.name == "GLOBALPHASE": qobj = gate.get_qobj(self.N) - U_list.append(qobj) - continue - - if gate.name in self.user_gates: - if gate.controls is not None: - raise ValueError( - "A user defined gate {} takes only " - "`targets` variable.".format(gate.name) - ) - func_or_oper = self.user_gates[gate.name] - if inspect.isfunction(func_or_oper): - func = func_or_oper - para_num = len(inspect.getfullargspec(func)[0]) - if para_num == 0: - qobj = func() - elif para_num == 1: - qobj = func(gate.arg_value) - else: - raise ValueError( - "gate function takes at most one parameters." - ) - elif isinstance(func_or_oper, Qobj): - qobj = func_or_oper - else: - raise ValueError("gate is neither function nor operator") + else: + qobj = self._get_gate_unitary(gate) if expand: all_targets = gate.get_all_qubits() qobj = expand_operator( qobj, dims=self.dims, targets=all_targets ) - else: - if expand: - qobj = gate.get_qobj(self.N, self.dims) - else: - qobj = gate.get_compact_qobj() U_list.append(qobj) return U_list + def _get_gate_unitary(self, gate): + if gate.name in self.user_gates: + if gate.controls is not None: + raise ValueError( + "A user defined gate {} takes only " + "`targets` variable.".format(gate.name) + ) + func_or_oper = self.user_gates[gate.name] + if inspect.isfunction(func_or_oper): + func = func_or_oper + para_num = len(inspect.getfullargspec(func)[0]) + if para_num == 0: + qobj = func() + elif para_num == 1: + qobj = func(gate.arg_value) + else: + raise ValueError( + "gate function takes at most one parameters." + ) + elif isinstance(func_or_oper, Qobj): + qobj = func_or_oper + else: + raise ValueError("gate is neither function nor operator") + else: + qobj = gate.get_compact_qobj() + return qobj + def compute_unitary(self): """Evaluates the matrix of all the gates in a quantum circuit. @@ -942,8 +925,9 @@ def compute_unitary(self): circuit_unitary : :class:`qutip.Qobj` Product of all gate arrays in the quantum circuit. """ - gate_list = self.propagators() - circuit_unitary = gate_sequence_product(gate_list) + sim = CircuitSimulator(self) + result = sim.run(qeye(self.dims)) + circuit_unitary = result.get_final_states()[0] return circuit_unitary def latex_code(self): diff --git a/src/qutip_qip/circuit/circuitsimulator.py b/src/qutip_qip/circuit/circuitsimulator.py index ddcb9aee..1e73d343 100644 --- a/src/qutip_qip/circuit/circuitsimulator.py +++ b/src/qutip_qip/circuit/circuitsimulator.py @@ -1,14 +1,12 @@ from itertools import product, chain -import os - +from operator import mul +from functools import reduce 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 import warnings @@ -225,7 +223,9 @@ def _gate_sequence_product(U_list, ind_list): def _gate_sequence_product_with_expansion(U_list, left_to_right=True): """ - Calculate the overall unitary matrix for a given list of unitary operations. + 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 ---------- @@ -259,12 +259,8 @@ class CircuitSimulator: def __init__( self, qc, - U_list=None, mode="state_vector_simulator", precompute_unitary=False, - state=None, - cbits=None, - measure_results=None, ): """ Simulate state evolution for Quantum Circuits. @@ -274,9 +270,6 @@ def __init__( qc : :class:`.QubitCircuit` Quantum Circuit to be simulated. - U_list: list of Qobj, optional - list of predefined unitaries corresponding to circuit. - mode: string, optional Specify if input state (and therefore computation) is in state-vector mode or in density matrix mode. @@ -289,123 +282,19 @@ def __init__( 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._qc = qc self.dims = qc.dims 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() - - if any(p is not None for p in (state, cbits, measure_results)): warnings.warn( - "Initializing the quantum state, cbits and measure_results " - "when initializing the simulator is deprecated. " - "The inputs are ignored. " - "They should, instead, be provided when running the simulation." + "Precomputing the full unitary is no longer supported. Switching to normal simulation mode." ) - 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 + @property + def qc(self): + return self._qc def initialize(self, state=None, cbits=None, measure_results=None): """ @@ -419,16 +308,13 @@ def initialize(self, state=None, cbits=None, measure_results=None): 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. """ - + # Initializing the unitary operators. if cbits and len(cbits) == self.qc.num_cbits: self.cbits = cbits elif self.qc.num_cbits > 0: @@ -436,47 +322,54 @@ def initialize(self, state=None, cbits=None, measure_results=None): else: self.cbits = None - self.state = 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) + self._state = ket2dm(state) else: - self.state = state - - self.probability = 1 - self.op_index = 0 - self.measure_results = measure_results - self.measure_ind = 0 + 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) - def _compute_unitary(self, U_list, inds_list): + @property + def state(self): """ - Compute unitary corresponding to a product of unitaries in U_list - and expand it to size of circuit. + The current state of the simulator as a `qutip.Qobj` - 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 + Returns: + `qutip.Qobj`: The current state of the simulator. """ - - 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, dims=self.qc.dims, targets=overall_inds - ) - return U_overall + 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 def run(self, state, cbits=None, measure_results=None): """ @@ -500,12 +393,13 @@ def run(self, state, cbits=None, measure_results=None): 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: + 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) + return CircuitResult(self.state, self._probability, self.cbits) def run_statistics(self, state, cbits=None): """ @@ -531,7 +425,7 @@ def run_statistics(self, state, cbits=None): cbits_results = [] num_measurements = len( - list(filter(lambda x: isinstance(x, Measurement), self.qc.gates)) + list(filter(lambda x: isinstance(x, Measurement), self._qc.gates)) ) for results in product("01", repeat=num_measurements): @@ -559,38 +453,40 @@ def _decimal_to_binary(decimal, length): binary = [int(s) for s in "{0:#b}".format(decimal)[2:]] return [0] * (length - len(binary)) + binary - 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[cbit_index] == control_value - for cbit_index, control_value in zip( - operation.classical_controls, - _decimal_to_binary( - operation.classical_control_value, - len(operation.classical_controls), - ), - ) - ] + 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), ) - if apply_gate: - if self.precompute_unitary: - U = expand_operator( - U, - dims=self.qc.dims, - targets=operation.get_all_qubits(), - ) - self._evolve_state(U) - else: - self._evolve_state(op) + 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) - self.op_index += 1 - return self.state + op = self._qc.gates[self._op_index] + self._op_index += 1 - def _evolve_state(self, U): + 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. @@ -599,41 +495,61 @@ def _evolve_state(self, U): 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": - self._evolve_ket(U) + state = U * state elif self.mode == "density_matrix_simulator": - self._evolve_dm(U) + 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 _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): + def _apply_measurement(self, operation, state): """ Applies measurement gate specified by operation to current state. @@ -642,29 +558,29 @@ def _apply_measurement(self, operation): 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 + 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] + 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)) - self.state = sum(p * s for s, p in zip(states, probabilities)) + state = sum(p * s for s, p in zip(states, probabilities)) else: raise NotImplementedError( "mode {} is not available.".format(self.mode) ) + return state class CircuitResult: diff --git a/tests/test_circuit.py b/tests/test_circuit.py index 98032bb2..e3ea1218 100644 --- a/tests/test_circuit.py +++ b/tests/test_circuit.py @@ -61,24 +61,6 @@ def _measurement_circuit(): return qc -def _simulators_sv(qc): - - sim_sv_precompute = CircuitSimulator(qc, mode="state_vector_simulator", - precompute_unitary=True) - sim_sv = CircuitSimulator(qc, mode="state_vector_simulator") - - return [sim_sv_precompute, sim_sv] - - -def _simulators_dm(qc): - - sim_dm_precompute = CircuitSimulator(qc, mode="density_matrix_simulator", - precompute_unitary=True) - sim_dm = CircuitSimulator(qc, mode="density_matrix_simulator") - - return [sim_dm_precompute, sim_dm] - - class TestQubitCircuit: """ A test class for the QuTiP functions for Circuit resolution. @@ -586,17 +568,16 @@ def test_runstatistics_teleportation(self): def test_measurement_circuit(self): qc = _measurement_circuit() - simulators = _simulators_sv(qc) + simulator = CircuitSimulator(qc) labels = ["00", "01", "10", "11"] for label in labels: state = bell_state(label) - for i, simulator in enumerate(simulators): - simulator.run(state) - if label[0] == "0": - assert simulator.cbits[0] == simulator.cbits[1] - else: - assert simulator.cbits[0] != simulator.cbits[1] + simulator.run(state) + if label[0] == "0": + assert simulator.cbits[0] == simulator.cbits[1] + else: + assert simulator.cbits[0] != simulator.cbits[1] def test_circuit_with_selected_measurement_result(self): qc = QubitCircuit(N=1, num_cbits=1) @@ -654,17 +635,17 @@ def test_wstate(self): _, probs_initial = fourth.measurement_comp_basis(state) - simulators = _simulators_sv(qc) + simulator = CircuitSimulator(qc) + + result = simulator.run_statistics(state) + final_states = result.get_final_states() + result_cbits = result.get_cbits() - for simulator in simulators: - result = simulator.run_statistics(state) - final_states = result.get_final_states() - result_cbits = result.get_cbits() + for i, final_state in enumerate(final_states): + _, probs_final = fourth.measurement_comp_basis(final_state) + np.testing.assert_allclose(probs_initial, probs_final) + assert sum(result_cbits[i]) == 1 - for i, final_state in enumerate(final_states): - _, probs_final = fourth.measurement_comp_basis(final_state) - np.testing.assert_allclose(probs_initial, probs_final) - assert sum(result_cbits[i]) == 1 _latex_template = r""" \documentclass[border=3pt]{standalone}