feat: Replace QermitPauli internals with UnitaryTableau (#201)

* Use only apply_gate

* Remove commute_coeff

* Use qubit pauli tensor in is_measureable

* Add from qubit_pauli_tensor

* Use qubit pauli tensors within reduce_qubits

* Use qubit pauli tensor in equality check

* Remove is_identity

* Use qubit pauli tensor in dagger

* Remove random pauli

* Use qubit pauli tensor in __hash__

* Remove pre_multiply

* Use qubit pauli tensor in from_pauli_list

* Use qubit pauli tensor in get_control_circuit

* Use qubit pauli tensor in circuit

* Use Op rather than OpType in apply_gate

* Correct wrong phase

* Use UnitaryTableau

* Remove unnecessary imports

* Repair high compute test

* Add qulacs dependency

* Move qulacs dependency

* Add projectq

* Reduce projectq version dep

* Revert projectq dependency change

* Do not use init

* Initialise with qubit pauli tensor

* Remove from pauli list

* Complete documentation

* Do not call dagger

* Remove notebook

qermit/coherent_pauli_checks/pauli_sampler.py

@@ -8,7 +8,7 @@
 from pytket import Circuit
 from pytket._tket.unit_id import UnitID
 from pytket.circuit import Bit, CircBox, Command, OpType, Qubit
-from pytket.pauli import Pauli, QubitPauliString
+from pytket.pauli import Pauli, QubitPauliString, QubitPauliTensor
 
 from qermit.noise_model.noise_model import NoiseModel
 from qermit.noise_model.qermit_pauli import QermitPauli
@@ -100,7 +100,7 @@ def add_pauli_checks_to_circbox(
                     pauli_check_circuit.add_barrier(command.args + [control_qubit])
 
                 end_stabiliser_list = [
-                    start_stabiliser.dagger()
+                    start_stabiliser.dagger
                     for start_stabiliser in start_stabiliser_list
                 ]
 
@@ -124,9 +124,8 @@ def add_pauli_checks_to_circbox(
                         # is not Clifford. It could be worth raising a clearer
                         # error.
                         end_stabiliser.apply_gate(
-                            clifford_command.op.type,
-                            clifford_command.qubits,
-                            params=clifford_command.op.params,
+                            op=clifford_command.op,
+                            qubits=clifford_command.qubits,
                         )
 
                     pauli_check_circuit.add_barrier(command.args + [control_qubit])
@@ -190,9 +189,12 @@ def sample(self, circ: Circuit) -> List[QermitPauli]:
         """
         return [
             QermitPauli(
-                Z_list=[1] * circ.n_qubits,
-                X_list=[0] * circ.n_qubits,
-                qubit_list=circ.qubits,
+                QubitPauliTensor(
+                    string=QubitPauliString(
+                        map={qubit: Pauli.Z for qubit in circ.qubits}
+                    ),
+                    coeff=1,
+                )
             )
         ]
 
@@ -208,9 +210,12 @@ def sample(self, circ: Circuit) -> List[QermitPauli]:
         """
         return [
             QermitPauli(
-                Z_list=[0] * circ.n_qubits,
-                X_list=[1] * circ.n_qubits,
-                qubit_list=circ.qubits,
+                QubitPauliTensor(
+                    string=QubitPauliString(
+                        map={qubit: Pauli.X for qubit in circ.qubits}
+                    ),
+                    coeff=1,
+                )
             )
         ]
 
@@ -239,21 +244,21 @@ def sample(
         :return: Random Pauli of length equal to the size of the circuit.
         """
         # TODO: Make sure sampling is done without replacement
-
         stabiliser_list: List[QermitPauli] = []
         while len(stabiliser_list) < self.n_checks:
-            Z_list = [self.rng.integers(2) for _ in circ.qubits]
-            X_list = [self.rng.integers(2) for _ in circ.qubits]
-
-            # Avoids using the identity string as it commutes with all errors
-            if any(Z == 1 for Z in Z_list) or any(X == 1 for X in X_list):
-                stabiliser_list.append(
-                    QermitPauli(
-                        Z_list=Z_list,
-                        X_list=X_list,
-                        qubit_list=circ.qubits,
-                    )
-                )
+            qpt = QubitPauliTensor(
+                string=QubitPauliString(
+                    map={
+                        qubit: self.rng.choice(
+                            numpy.array([Pauli.X, Pauli.Y, Pauli.Z, Pauli.I])
+                        )
+                        for qubit in circ.qubits
+                    }
+                ),
+                coeff=1,
+            )
+            if qpt != QubitPauliTensor():
+                stabiliser_list.append(QermitPauli(qpt))
 
         return stabiliser_list
 
@@ -347,6 +352,6 @@ def sample(
         #     )
 
         return [
-            QermitPauli.from_qubit_pauli_string(smallest_commute_prob_pauli)
+            QermitPauli(QubitPauliTensor(string=smallest_commute_prob_pauli, coeff=1))
             for smallest_commute_prob_pauli in smallest_commute_prob_pauli_list
         ]

qermit/noise_model/noise_model.py

@@ -12,7 +12,7 @@
 from numpy.typing import NDArray
 from pytket import Circuit, Qubit
 from pytket.circuit import OpType
-from pytket.pauli import Pauli, QubitPauliString
+from pytket.pauli import Pauli, QubitPauliString, QubitPauliTensor
 from scipy.linalg import fractional_matrix_power  # type: ignore
 
 from .qermit_pauli import QermitPauli
@@ -99,9 +99,14 @@ def to_ptm(self) -> Tuple[NDArray, Dict[Tuple[Pauli, ...], int]]:
         # PTM entry as a sum pf error weights multiplied by +/-1
         # Depending on commutation relations.
         for pauli_tuple, index in pauli_index.items():
-            pauli = QermitPauli.from_pauli_iterable(
-                pauli_iterable=pauli_tuple,
-                qubit_list=[Qubit(i) for i in range(self.n_qubits)],
+            pauli = QermitPauli(
+                QubitPauliTensor(
+                    string=QubitPauliString(
+                        paulis=list(pauli_tuple),
+                        qubits=[Qubit(i) for i in range(self.n_qubits)],
+                    ),
+                    coeff=1,
+                )
             )
 
             # Can add the identity error rate.
@@ -111,14 +116,23 @@ def to_ptm(self) -> Tuple[NDArray, Dict[Tuple[Pauli, ...], int]]:
             ptm[index][index] += self.identity_error_rate
 
             for error, error_rate in self.distribution.items():
-                error_pauli = QermitPauli.from_pauli_iterable(
-                    pauli_iterable=error,
-                    qubit_list=[Qubit(i) for i in range(self.n_qubits)],
+                error_pauli = QermitPauli(
+                    QubitPauliTensor(
+                        string=QubitPauliString(
+                            paulis=list(error),
+                            qubits=[Qubit(i) for i in range(self.n_qubits)],
+                        ),
+                        coeff=1,
+                    )
                 )
-
-                ptm[index][index] += error_rate * QermitPauli.commute_coeff(
-                    pauli_one=pauli, pauli_two=error_pauli
+                commute_coeff = (
+                    1
+                    if pauli.qubit_pauli_tensor.commutes_with(
+                        error_pauli.qubit_pauli_tensor
+                    )
+                    else -1
                 )
+                ptm[index][index] += error_rate * commute_coeff
 
         # Some checks that the form of the PTM is correct.
         identity = tuple(Pauli.I for _ in range(self.n_qubits))
@@ -181,17 +195,32 @@ def from_ptm(
         # is the matrix of commutation values.
         commutation_matrix = np.zeros(ptm.shape)
         for pauli_one_tuple, index_one in pauli_index.items():
-            pauli_one = QermitPauli.from_pauli_iterable(
-                pauli_iterable=pauli_one_tuple,
-                qubit_list=[Qubit(i) for i in range(len(pauli_one_tuple))],
+            pauli_one = QermitPauli(
+                QubitPauliTensor(
+                    string=QubitPauliString(
+                        paulis=list(pauli_one_tuple),
+                        qubits=[Qubit(i) for i in range(len(pauli_one_tuple))],
+                    ),
+                    coeff=1,
+                )
             )
             for pauli_two_tuple, index_two in pauli_index.items():
-                pauli_two = QermitPauli.from_pauli_iterable(
-                    pauli_iterable=pauli_two_tuple,
-                    qubit_list=[Qubit(i) for i in range(len(pauli_two_tuple))],
+                pauli_two = QermitPauli(
+                    QubitPauliTensor(
+                        string=QubitPauliString(
+                            paulis=list(pauli_two_tuple),
+                            qubits=[Qubit(i) for i in range(len(pauli_two_tuple))],
+                        ),
+                        coeff=1,
+                    )
                 )
-                commutation_matrix[index_one][index_two] = QermitPauli.commute_coeff(
-                    pauli_one=pauli_one, pauli_two=pauli_two
+
+                commutation_matrix[index_one][index_two] = (
+                    1
+                    if pauli_one.qubit_pauli_tensor.commutes_with(
+                        pauli_two.qubit_pauli_tensor
+                    )
+                    else -1
                 )
 
         error_rate_list = np.matmul(ptm.diagonal(), np.linalg.inv(commutation_matrix))
@@ -607,7 +636,8 @@ def counter_propagate(
         for _ in range(n_counts):
             pauli_error = self.random_propagate(cliff_circ, **kwargs)
 
-            if not pauli_error.is_identity:
+            # Check if the error is the identity.
+            if pauli_error.qubit_pauli_tensor.string != QubitPauliString():
                 error_counter.update([pauli_error])
 
         return error_counter
@@ -627,13 +657,13 @@ def random_propagate(
         :raises Exception: Raised if direction is invalid.
         :return: Resulting logical error.
         """
-
-        # Create identity error.
-        qubit_list = cliff_circ.qubits
         pauli_error = QermitPauli(
-            Z_list=[0] * len(qubit_list),
-            X_list=[0] * len(qubit_list),
-            qubit_list=qubit_list,
+            QubitPauliTensor(
+                string=QubitPauliString(
+                    map={qubit: Pauli.I for qubit in cliff_circ.qubits}
+                ),
+                coeff=1,
+            )
         )
 
         # Commands are ordered in reverse or original order depending on which
@@ -656,9 +686,8 @@ def random_propagate(
             if direction == Direction.forward:
                 # Apply gate to total error.
                 pauli_error.apply_gate(
-                    op_type=command.op.type,
+                    op=command.op,
                     qubits=cast(List[Qubit], command.args),
-                    params=command.op.params,
                 )
             # Add noise operation if appropriate.
             if command.op.type in self.noisy_gates:
@@ -685,9 +714,8 @@ def random_propagate(
                 # which has the same effect on the pauli as pushing through the
                 # dagger.
                 pauli_error.apply_gate(
-                    op_type=command.op.dagger.type,
+                    op=command.op.dagger,
                     qubits=cast(List[Qubit], command.args),
-                    params=command.op.dagger.params,
                 )
 
         return pauli_error