Merge pull request #177 from UCL-CCS/clifford_sim

Added circuit simulator functionality - efficient if one uses exclusi…

symmer/evolution/__init__.py

@@ -3,4 +3,5 @@
 from .gate_library import *
 from .utils import get_CNOT_connectivity_graph, topology_match_score
 from .decomposition import qasm_to_PauliwordOp, PauliwordOp_to_QuantumCircuit
+from .circuit_symmerlator import CircuitSymmerlator
 from .variational_optimization import VQE_Driver, ADAPT_VQE

symmer/evolution/circuit_symmerlator.py

@@ -0,0 +1,198 @@
+from typing import List
+from symmer import PauliwordOp, QuantumState
+from qiskit import QuantumCircuit
+import numpy as np
+
+class CircuitSymmerlator:
+    """
+    Symmer circuit simulator
+
+    Each Clifford gate is expanded as a sequence of pi/2 rotations
+    (correct up to some global phase that cancels in the inner product calculation).
+
+    Non-Clifford gates also included, although note the use of these is inefficient.
+    
+    """
+    def __init__(self, n_qubits: int) -> None:
+        """
+        """
+        self.n_qubits = n_qubits
+        self.sequence = []
+        self.gate_map = {
+            'x':self.X,'y':self.Y,'z':self.Z,
+            'rx':self.RX,'ry':self.RY,'rz':self.RZ,
+            'sx':self.sqrtX,'sy':self.sqrtY,'sz':self.sqrtZ,
+            'cx':self.CX,'cy':self.CY,'cz':self.CZ,
+            'h':self.H,'s':self.S,'sdg':self.Sdag,
+            '':self.R,'t':self.T,'ccx':self.Toffoli,'swap':self.SWAP
+        }
+
+    def get_rotation_string(self, pauli: str, indices: List[int]):
+        """ Given a Pauli string and list of corresponding qubit indices, return the PauliwordOp
+        """
+        pauli = list(pauli)
+        assert len(pauli)==len(indices), 'Number of Paulis and indices do not match'
+        assert set(pauli).issubset({'I','X','Y','Z'}), 'Pauli operators are either I, X, Y or Z.'
+        R = ['I']*self.n_qubits
+        for i,P in zip(indices, pauli):
+            R[i] = P
+        return PauliwordOp.from_list([''.join(R)])
+
+    def pi_2_multiple(self, multiple: int):
+        """ 
+        For multiple % 4 = 0,1,2,3 and rotation operator R we have the following behaviour:
+            0: +I
+            1: +R
+            2: -I
+            3: -R
+        multiplied by the anticommuting component of the target operator.
+        """
+        return np.pi/2*multiple
+
+    #################################################################
+    ######################  CLIFFORD GATES  #########################
+    #################################################################
+
+    def X(self, index: int) -> None:
+        """ Pauli iX gate """
+        self.sequence.append( (self.get_rotation_string('X', [index]), self.pi_2_multiple(+2)) )
+
+    def Y(self, index: int) -> None:
+        """ Pauli iY gate """
+        self.sequence.append( (self.get_rotation_string('Y', [index]), self.pi_2_multiple(+2)) )
+
+    def Z(self, index: int) -> None:
+        """ Pauli iZ gate """
+        self.sequence.append( (self.get_rotation_string('Z', [index]), self.pi_2_multiple(+2)) )
+
+    def H(self, index: int) -> None:
+        """ iH Hadamard gate """
+        self.sequence.append( (self.get_rotation_string('Z', [index]), self.pi_2_multiple(+2)) )
+        self.sequence.append( (self.get_rotation_string('Y', [index]), self.pi_2_multiple(+1)) )
+
+    def S(self, index: int) -> None:
+        """ √(+i)S gate """
+        self.sequence.append( (self.get_rotation_string('Z', [index]), self.pi_2_multiple(+1)) )
+
+    def Sdag(self, index: int) -> None:
+        """ -√(-i)S† gate """
+        self.sequence.append( (self.get_rotation_string('Z', [index]), self.pi_2_multiple(+3)) )
+
+    def sqrtX(self, index: int) -> None:
+        """ √(iX) gate """
+        self.sequence.append( (self.get_rotation_string('X', [index]), self.pi_2_multiple(+1)) )
+
+    def sqrtY(self, index: int) -> None:
+        """ √(iY) gate """
+        self.sequence.append( (self.get_rotation_string('Y', [index]), self.pi_2_multiple(+1)) )
+
+    def sqrtZ(self, index: int) -> None:
+        """ √(iZ) gate """
+        self.sequence.append( (self.get_rotation_string('Z', [index]), self.pi_2_multiple(+1)) )
+
+    def CX(self, control: int, target: int) -> None:
+        """ √(-i)CX gate """
+        self.sequence.append( (self.get_rotation_string('ZX', [control, target]), self.pi_2_multiple(+1)) )
+        self.sequence.append( (self.get_rotation_string('ZI', [control, target]), self.pi_2_multiple(+3)) )
+        self.sequence.append( (self.get_rotation_string('IX', [control, target]), self.pi_2_multiple(+3)) )
+
+    def CY(self, control: int, target: int) -> None:
+        """ √(-i)CY gate """
+        self.sequence.append( (self.get_rotation_string('ZY', [control, target]), self.pi_2_multiple(+1)) )
+        self.sequence.append( (self.get_rotation_string('ZI', [control, target]), self.pi_2_multiple(+3)) )
+        self.sequence.append( (self.get_rotation_string('IY', [control, target]), self.pi_2_multiple(+3)) )
+
+    def CZ(self, control: int, target: int) -> None:
+        """ √(-i)CZ gate """
+        self.sequence.append( (self.get_rotation_string('ZZ', [control, target]), self.pi_2_multiple(+1)) )
+        self.sequence.append( (self.get_rotation_string('ZI', [control, target]), self.pi_2_multiple(+3)) )
+        self.sequence.append( (self.get_rotation_string('IZ', [control, target]), self.pi_2_multiple(+3)) )
+
+    def SWAP(self, qubit_1: int, qubit_2: int) -> None:
+        """ Swap qubits 1 and 2 """
+        self.CX(qubit_1,qubit_2)
+        self.CX(qubit_2,qubit_1)
+        self.CX(qubit_1,qubit_2)
+
+    #################################################################
+    ####################  NON-CLIFFORD GATES  #######################
+    #################### WARNING: INEFFICIENT #######################
+    #################################################################
+
+    def R(self, pauli: str, indices: List[int], angle: float) -> None:
+        """ Arbitrary rotation gate """
+        self.sequence.append( (self.get_rotation_string(pauli, indices), -angle) )
+
+    def RX(self, index: int, angle: float) -> None:
+        """ Pauli X rotation """
+        self.R('X', [index], angle)
+
+    def RY(self, index: int, angle: float) -> None:
+        """ Pauli Y rotation """
+        self.R('Y', [index], angle)
+
+    def RZ(self, index: int, angle: float) -> None:
+        """ Pauli Z rotation """
+        self.R('Z', [index], angle)
+
+    def T(self, index: int, angle: float) -> None:
+        """ T gate """
+        raise NotImplementedError()
+
+    def Toffoli(self, control_1: int, control_2: int, target: int) -> None:
+        """ Doubly-controlled X gate """
+        raise NotImplementedError()
+
+    #################################################################
+    ######################  GATE EXECUTION  #########################
+    #################################################################
+
+    def apply_sequence(self, operator: PauliwordOp) -> PauliwordOp:
+        """ Apply the stored sequence of rotations on the input operator
+        """
+        assert operator.n_qubits == self.n_qubits, 'The operator is defined over a different number of qubits'
+        return operator.perform_rotations(self.sequence[::-1])
+
+    def evaluate(self, operator: PauliwordOp) -> float:
+        """ Evaluate the stored rotations on the input operator
+        """
+        rotated_op = self.apply_sequence(operator)
+        expval = 0
+        for rotated_str, coeff in rotated_op.to_dictionary.items():
+            if set(rotated_str).issubset({'I','Z'}):
+                expval += coeff
+        return expval
+
+    @classmethod
+    def from_qasm(cls, qasm: str, angle_factor: int=1) -> "CircuitSymmerlator":
+        """ Initialize the simulator from a QASM circuit
+        """
+        instructions = qasm.split(';\n')[:-1]
+        qasm_version = instructions.pop(0)
+        inclusions   = instructions.pop(0)
+        registers    = instructions.pop(0)
+        n_qubits     = int(registers.split(' ')[1][2:-1])
+
+        self = cls(n_qubits)
+        pi = np.pi # for evaluating strings like '3*pi/2'
+        for step in instructions:
+            gate, qubits  = step.split(' ')
+            qubits = [int(q[2:-1]) for q in qubits.split(',')]
+            extract_angle = gate.split('(')
+            if len(extract_angle) == 1:
+                gate  = extract_angle[0]
+                angle = None
+            else:
+                gate, angle = extract_angle
+                angle = eval(angle[:-1])
+            if angle is not None:
+                self.gate_map[gate](*qubits, angle=angle_factor*angle)
+            else:
+                self.gate_map[gate](*qubits)
+        return self
+
+    @classmethod
+    def from_qiskit(cls, circuit: QuantumCircuit) -> "CircuitSymmerlator":
+        """ Initialize the simulator from a Qiskit QuantumCircuit
+        """
+        return cls.from_qasm(circuit.qasm(), angle_factor=-1)

symmer/evolution/gate_library.py

@@ -111,6 +111,22 @@ def CX(n_qubits:int, control:int, target:int) -> PauliwordOp:
     _Had = Had(n_qubits, target)
     return _Had * CZ(n_qubits, control, target) * _Had
 
+def CY(n_qubits:int, control:int, target:int) -> PauliwordOp:
+    """ 
+    Controlled Y gate
+
+    Args:
+        n_qubits (int): Number of qubits.
+        control (int): Qubit index at which will act as a control qubit.
+        target (int): Qubit index at which the operation 'X' has to be applied if the control qubit is in |1> state.
+
+    Returns:
+        PauliwordOp representing the Controlled-X (CX) gate applied to the specified control and target qubits in a system of 'n_qubits' qubits.
+    """
+    _Had = Had(n_qubits, target)
+    _S   = S(n_qubits, target)
+    return _S * _Had * CZ(n_qubits, control, target) * _Had * _S.dagger
+
 def RX(n_qubits:int, index:int, angle:float) -> PauliwordOp:
     """ 
     Rotation-X gate

symmer/operators/base.py

@@ -829,7 +829,7 @@ def expval(self, psi: "QuantumState") -> complex:
         Returns:
             complex: The expectation value.
         """
-        if self.n_terms > psi.n_terms:
+        if self.n_terms > psi.n_terms and psi.n_terms > 10:
             return (psi.dagger * self * psi).real
         else:
             if self.n_terms > 1:

tests/test_evolution/test_circuit_symmerlator.py

@@ -0,0 +1,142 @@
+import pytest
+import numpy as np
+from symmer import PauliwordOp
+from symmer.evolution.gate_library import *
+from symmer.evolution.circuit_symmerlator import CircuitSymmerlator
+from qiskit.quantum_info import Statevector, random_clifford
+
+def test_random_cliffords():
+    """ Test the simulator on lots of random 5-qubit Clifford circuits
+    """
+    for _ in range(50):
+        n_q = 5
+        # generate random Clifford circuit
+        qc = random_clifford(n_q).to_circuit()
+        sv = Statevector(qc).data.reshape(-1,1)
+        # ... and a random observable
+        observable = PauliwordOp.random(n_q, 100, complex_coeffs=False)
+        observable_matrix = observable.reindex(list(range(0,n_q))[::-1]).to_sparse_matrix.toarray() # reverse order because qiskit
+        # Than, intialize the circuit simulator and check it matches direct statevector calculation.
+        CS = CircuitSymmerlator.from_qiskit(qc)
+        assert np.isclose(
+            (sv.conj().T @ observable_matrix @ sv)[0,0],
+            CS.evaluate(observable)
+        )
+
+#### INDIVIDUAL GATE TESTS BELOW ####
+
+one_qubit_paulis = [
+    'I', 'X', 'Y', 'Z'
+]
+two_qubit_paulis = [
+    'IZ','ZI','ZZ','IX','XI','XX','IY','YI','YY',
+    'XY','XZ','YX','YZ','ZX','ZY','II'
+]
+
+@pytest.mark.parametrize("P", one_qubit_paulis)
+def test_X_gate(P):
+    CS = CircuitSymmerlator(1)
+    CS.X(0)
+    P = PauliwordOp.from_list([P])
+    assert CS.apply_sequence(P) == X(1,0)*P*X(1,0)
+
+@pytest.mark.parametrize("P", one_qubit_paulis)
+def test_Y_gate(P):
+    CS = CircuitSymmerlator(1)
+    CS.Y(0)
+    P = PauliwordOp.from_list([P])
+    assert CS.apply_sequence(P) == Y(1,0)*P*Y(1,0)
+
+@pytest.mark.parametrize("P", one_qubit_paulis)
+def test_Z_gate(P):
+    CS = CircuitSymmerlator(1)
+    CS.Z(0)
+    P = PauliwordOp.from_list([P])
+    assert CS.apply_sequence(P) == Z(1,0)*P*Z(1,0)
+
+@pytest.mark.parametrize("P", one_qubit_paulis)
+def test_RX_gate(P):
+    CS = CircuitSymmerlator(1)
+    angle = np.random.random()
+    CS.RX(0, angle)
+    P = PauliwordOp.from_list([P])
+    assert CS.apply_sequence(P) == RX(1,0,-angle)*P*RX(1,0,angle)
+
+@pytest.mark.parametrize("P", one_qubit_paulis)
+def test_RY_gate(P):
+    CS = CircuitSymmerlator(1)
+    angle = np.random.random()
+    CS.RY(0, angle)
+    P = PauliwordOp.from_list([P])
+    assert CS.apply_sequence(P) == RY(1,0,-angle)*P*RY(1,0,angle)
+
+@pytest.mark.parametrize("P", one_qubit_paulis)
+def test_RZ_gate(P):
+    CS = CircuitSymmerlator(1)
+    angle = np.random.random()
+    CS.RZ(0, angle)
+    P = PauliwordOp.from_list([P])
+    assert CS.apply_sequence(P) == RZ(1,0,-angle)*P*RZ(1,0,angle)
+
+@pytest.mark.parametrize("P", one_qubit_paulis)
+def test_H_gate(P):
+    CS = CircuitSymmerlator(1)
+    CS.H(0)
+    P = PauliwordOp.from_list([P])
+    assert CS.apply_sequence(P) == Had(1,0)*P*Had(1,0)
+
+@pytest.mark.parametrize("P", one_qubit_paulis)
+def test_S_gate(P):
+    CS = CircuitSymmerlator(1)
+    CS.S(0)
+    P = PauliwordOp.from_list([P])
+    assert CS.apply_sequence(P) == S(1,0).dagger * P * S(1,0)
+
+@pytest.mark.parametrize("P", one_qubit_paulis)
+def test_Sdag_gate(P):
+    CS = CircuitSymmerlator(1)
+    CS.Sdag(0)
+    P = PauliwordOp.from_list([P])
+    assert CS.apply_sequence(P) == S(1,0) * P * S(1,0).dagger
+
+@pytest.mark.parametrize("P", one_qubit_paulis)
+def test_sqrtX_gate(P):
+    CS = CircuitSymmerlator(1)
+    CS.sqrtX(0)
+    P = PauliwordOp.from_list([P])
+    assert CS.apply_sequence(P) == RX(1,0,+np.pi/2)*P*RX(1,0,-np.pi/2)
+
+@pytest.mark.parametrize("P", one_qubit_paulis)
+def test_sqrtY_gate(P):
+    CS = CircuitSymmerlator(1)
+    CS.sqrtY(0)
+    P = PauliwordOp.from_list([P])
+    assert CS.apply_sequence(P) == RY(1,0,+np.pi/2)*P*RY(1,0,-np.pi/2)
+
+@pytest.mark.parametrize("P", one_qubit_paulis)
+def test_sqrtZ_gate(P):
+    CS = CircuitSymmerlator(1)
+    CS.sqrtZ(0)
+    P = PauliwordOp.from_list([P])
+    assert CS.apply_sequence(P) == RZ(1,0,+np.pi/2)*P*RZ(1,0,-np.pi/2)
+
+@pytest.mark.parametrize("P", two_qubit_paulis)
+def test_CX_gate(P):
+    CS = CircuitSymmerlator(2)
+    CS.CX(0,1)
+    P = PauliwordOp.from_list([P])
+    assert CS.apply_sequence(P) == CX(2,0,1)*P*CX(2,0,1)
+
+@pytest.mark.parametrize("P", two_qubit_paulis)
+def test_CY_gate(P):
+    CS = CircuitSymmerlator(2)
+    CS.CY(0,1)
+    P = PauliwordOp.from_list([P])
+    assert CS.apply_sequence(P) == CY(2,0,1)*P*CY(2,0,1)
+
+@pytest.mark.parametrize("P", two_qubit_paulis)
+def test_CZ_gate(P):
+    CS = CircuitSymmerlator(2)
+    CS.CZ(0,1)
+    P = PauliwordOp.from_list([P])
+    assert CS.apply_sequence(P) == CZ(2,0,1)*P*CZ(2,0,1)