Stim clifford simulator integration

.github/workflows/continuous_integration.yml

@@ -44,6 +44,7 @@ jobs:
         pip install cirq
         pip install projectq
         pip install pennylane
+        pip install stim
       if: always()
 
     - name: Install Microsoft qsharp/qdk

tangelo/helpers/math.py

@@ -53,3 +53,29 @@ def bool_col_echelon(bool_array):
                 pivot -= 1
 
     return bool_array
+
+
+def arrays_almost_equal_up_to_global_phase(array1, array2, atol=1e-6):
+    """
+    Checks if two arrays are almost equal up to a global phase.
+
+    Args:
+        array1 (array): Self-explanatory.
+        array2 (array): Self-explanatory.
+        atol (float) : Optional, absolute tolerance
+
+    Returns:
+        bool : True if arrays are almost equal up to a global phase, False otherwise.
+    """
+    if len(array1) != len(array2):
+        return False
+
+    array1 = np.asarray(array1)
+    array2 = np.asarray(array2)
+
+    if np.allclose(array1, array2, atol=atol):
+        return True
+
+    # Check for global phase difference
+    phase_diff = np.angle(array1[0] / array2[0])
+    return np.allclose(array1, array2 * np.exp(1j * phase_diff), atol=atol)

tangelo/helpers/utils.py

@@ -73,11 +73,12 @@ def new_func(*args, **kwargs):
 
 
 # List all built-in backends supported
-all_backends = {"qulacs", "qiskit", "cirq", "braket", "projectq", "qdk", "pennylane", "sympy"}
+all_backends = {"qulacs", "qiskit", "cirq", "braket", "projectq", "qdk", "pennylane", "sympy", "stim"}
 all_backends_simulator = {"qulacs", "qiskit", "cirq", "qdk"}
 sv_backends_simulator = {"qulacs", "qiskit", "cirq"}
 symbolic_backends = {"sympy"}
 chem_backends = {"pyscf", "psi4"}
+clifford_backends_simulator = {"stim"}
 
 # Dictionary mapping package names to their identifier in this module
 packages = {p: p for p in all_backends}
@@ -88,6 +89,7 @@ def new_func(*args, **kwargs):
 installed_simulator = installed_backends & all_backends_simulator
 installed_sv_simulator = installed_backends & sv_backends_simulator
 installed_chem_backends = {p_id for p_id in chem_backends if is_package_installed(p_id)}
+installed_clifford_simulators = {p_id for p_id in clifford_backends_simulator if is_package_installed(p_id)}
 
 # Check if qulacs installed (better performance for tests). If not, defaults to cirq (always installed with openfermion)
 default_simulator = "qulacs" if "qulacs" in installed_sv_simulator else "cirq"

tangelo/linq/gate.py

@@ -16,7 +16,7 @@
 gate operation, without tying it to a particular backend or an underlying
 mathematical operation.
 """
-from math import pi
+from math import pi, isclose
 from typing import Union
 
 from numpy import integer, ndarray, floating
@@ -36,6 +36,8 @@
                     "CNOT", "CX", "CY", "CZ", "CRX", "CRY", "CRZ", "CPHASE", "XX", "SWAP",
                     "CSWAP"}
 
+CLIFFORD_GATES = {"H", "S", "X", "Z", "Y", "SDAG", "CNOT", "CX", "CY", "CZ", "SWAP", "CXSWAP"}
+
 
 class Gate(dict):
     """An abstract gate class that exposes all the gate information such as gate
@@ -180,6 +182,24 @@ def inverse(self):
             raise AttributeError(f"{self.name} is not an invertible gate when parameter is {self.parameter}")
         return Gate(self.name, self.target, self.control, new_parameter, self.is_variational)
 
+    def is_clifford(self, abs_tol=1e-4):
+        """
+        Check if a quantum gate is a Clifford gate.
+
+        Args:
+            abs_tol (float) : Optional, Absolute tolerance for the difference between gate parameter clifford parameter values
+
+        Returns:
+            bool: True if the gate is in Clifford gate list or has a parameter corresponding to Clifford points,
+                  False otherwise.
+        """
+        if self.name in CLIFFORD_GATES:
+            return True
+        elif self.name in {"RX", "RY", "RZ", "PHASE"}:
+            return isclose(self.parameter % (pi / 2), 0, abs_tol=abs_tol)
+        else:
+            return False
+
     def serialize(self):
         return {"type": "Gate",
                 "params": {"name": self.name, "target": self.target,

tangelo/linq/helpers/circuits/clifford_circuits.py

@@ -0,0 +1,89 @@
+# Copyright 2023 Good Chemistry Company.
+#
+# 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 math import pi, isclose
+
+from tangelo.linq import Gate
+
+
+def decompose_gate_to_cliffords(gate, abs_tol=1e-4):
+    """
+    Decomposes a single qubit parameterized gate into Clifford gates.
+
+    Args:
+        gate (Gate): The gate to be decomposed.
+        abs_tol (float): Optional, absolute tolerance for value comparison (default: 1e-4).
+
+    Returns:
+        list: A list of Clifford gates representing the decomposition.
+
+    Raises:
+        ValueError: If parameterized gate cannot be decomposed into Clifford gates.
+
+    """
+    # Return an error if the gate is not in the Clifford gate set
+    if not gate.is_clifford(abs_tol):
+        raise ValueError(f"Error. The following gate cannot be decomposed into Clifford gates:\n {gate}")
+    # If gate is not in the parameterized gate set, it decomposes to itself
+    elif gate.name not in {"RX", "RY", "RZ", "PHASE"}:
+        return gate
+    # If gate parameter is close to 0, the gate operation is the identity
+    elif isclose(gate.parameter, 0, abs_tol=abs_tol):
+        return []
+
+    # Find which Clifford parameter gate parameter corresponds to.
+    clifford_values = [0, pi, pi / 2, -pi / 2]
+    clifford_parameter = next((value for value in clifford_values if
+                               isclose(gate.parameter % (2 * pi), value % (2 * pi), abs_tol=abs_tol)), None)
+
+    if clifford_parameter is None:
+        raise ValueError(f"Error: Parameterized gate {gate} cannot be decomposed into Clifford gates")
+
+    gate_list = []
+    if gate.name == "RY":
+        if clifford_parameter == -pi / 2:
+            gate_list = [Gate("H", gate.target), Gate("Z", gate.target)]
+        elif clifford_parameter == pi / 2:
+            gate_list = [Gate("Z", gate.target), Gate("H", gate.target)]
+
+        elif clifford_parameter == pi:
+            gate_list = [Gate("SDAG", gate.target), Gate("Y", gate.target), Gate("SDAG", gate.target)]
+
+    elif gate.name == "RX":
+        if clifford_parameter == -pi / 2:
+            gate_list = [Gate("S", gate.target), Gate("H", gate.target), Gate("S", gate.target)]
+        elif clifford_parameter == pi / 2:
+            gate_list = [Gate("SDAG", gate.target), Gate("H", gate.target), Gate("SDAG", gate.target)]
+        elif clifford_parameter == pi:
+            gate_list = [Gate("SDAG", gate.target), Gate("X", gate.target), Gate("SDAG", gate.target)]
+
+    elif gate.name == "RZ":
+        if clifford_parameter == -pi / 2:
+            gate_list = [Gate("H", gate.target), Gate("S", gate.target), Gate("H", gate.target), Gate("S", gate.target),
+                         Gate("H", gate.target)]
+        elif clifford_parameter == pi / 2:
+            gate_list = [Gate("H", gate.target), Gate("SDAG", gate.target), Gate("H", gate.target),
+                         Gate("SDAG", gate.target), Gate("H", gate.target)]
+        elif clifford_parameter == pi:
+            gate_list = [Gate("H", gate.target), Gate("SDAG", gate.target), Gate("X", gate.target), Gate("SDAG", gate.target), Gate("H", gate.target)]
+
+    elif gate.name == "PHASE":
+        if clifford_parameter == -pi / 2:
+            gate_list = [Gate("SDAG", gate.target)]
+        elif clifford_parameter == pi / 2:
+            gate_list = [Gate("S", gate.target)]
+        elif clifford_parameter == pi:
+            gate_list = [Gate("Z", gate.target)]
+
+    return gate_list

tangelo/linq/helpers/circuits/tests/test_clifford_circuits.py

@@ -0,0 +1,43 @@
+# Copyright 2023 Good Chemistry Company.
+#
+# 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.
+
+import unittest
+
+import numpy as np
+
+from tangelo.linq import Gate, Circuit
+from tangelo.linq.helpers.circuits.clifford_circuits import decompose_gate_to_cliffords
+from tangelo.linq.translator.translate_cirq import translate_c_to_cirq
+from tangelo.helpers.math import arrays_almost_equal_up_to_global_phase
+
+clifford_angles = [0, np.pi/2, np.pi, -np.pi/2, 3*np.pi/2, 7*np.pi, -5*np.pi]
+non_clifford_gate = Gate("RY", 0, parameter=np.pi/3)
+
+
+class CliffordCircuitTest(unittest.TestCase):
+
+    def test_decompose_gate_to_cliffords(self):
+        """Test if gate decomposition returns correct sequence. Uses Cirq to validate that the unitaries returned
+        are equal up to a global phase"""
+        for gate in ["RY", "RX", "RZ", "PHASE"]:
+            for angle in clifford_angles:
+                ref_gate = Gate(gate, 0, parameter=angle)
+                u_ref = translate_c_to_cirq(Circuit([ref_gate]))
+                u_decompose = translate_c_to_cirq(Circuit(decompose_gate_to_cliffords(ref_gate)))
+                arrays_almost_equal_up_to_global_phase(u_ref.unitary(), u_decompose.unitary())
+
+    def test_non_clifford_gates(self):
+        """Test if non-clifford gate raises value error"""
+
+        self.assertRaises(ValueError, decompose_gate_to_cliffords, non_clifford_gate)

tangelo/linq/target/__init__.py

@@ -18,10 +18,11 @@
 from .target_qulacs import QulacsSimulator
 from .target_qdk import QDKSimulator
 from .target_sympy import SympySimulator
-from tangelo.helpers.utils import all_backends_simulator
+from .target_stim import StimSimulator
+from tangelo.helpers.utils import all_backends_simulator, clifford_backends_simulator
 
 
-target_dict = {"qiskit": QiskitSimulator, "cirq": CirqSimulator, "qdk": QDKSimulator, "qulacs": QulacsSimulator, "sympy": SympySimulator}
+target_dict = {"qiskit": QiskitSimulator, "cirq": CirqSimulator, "qdk": QDKSimulator, "qulacs": QulacsSimulator, "sympy": SympySimulator, "stim": StimSimulator}
 
 # Generate backend info dictionary
 backend_info = {sim_id: target_dict[sim_id].backend_info() for sim_id in all_backends_simulator}

tangelo/linq/target/backend.py

@@ -456,6 +456,7 @@ def _get_expectation_value_from_statevector(self, qubit_operator, state_prep_cir
             complex: The expectation value of this operator with regards to the
                 state preparation.
         """
+
         n_qubits = state_prep_circuit.width
 
         expectation_value = 0.