Merge pull request #225 from BoxiLi/simulator_refactor.py

Refactor the gate level `CircuitSimulator`
- Avoid expanding all propagators and then applying them one by one. Only expand when applying the unitary. This saves memory and also make it easier for further improvement, i.e., exploring other ways to perform the matrix multiplication (e.g. `einsum`).
- Switch some public members to private. Those members that track the execution of the circuit evaluation should not be public.
- Small changes to improve the clarity.

Author: BoxiLi

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
 =========================

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):