Add then method for composing two Clifford tableaux

cirq-core/cirq/qis/clifford_tableau.py

@@ -17,7 +17,6 @@
 
 import cirq
 from cirq import protocols
-from cirq.ops.dense_pauli_string import DensePauliString
 from cirq.value import big_endian_int_to_digits
 
 
@@ -84,6 +83,10 @@ def rs(self, new_rs: np.array) -> None:
         assert np.shape(new_rs) == (2 * self.n,)
         self._rs[:-1] = np.array(new_rs).astype(bool)
 
+    def matrix(self) -> np.array:
+        """Returns a 2n * 2n matrix of Clifford table."""
+        return np.concatenate([self.xs, self.zs], axis=1)
+
     def _json_dict_(self) -> Dict[str, Any]:
         return protocols.obj_to_dict_helper(self, ['n', 'rs', 'xs', 'zs'])
 
@@ -176,6 +179,48 @@ def _str_full_(self) -> str:
 
         return string
 
+    def then(self, second: 'CliffordTableau') -> 'CliffordTableau':
+        """Returns a composed CliffordTableau of this tableau and the second tableau.
+
+        Then composed tableau is equal to (up to global phase) the composed
+        unitary operation of the two tableaux, i.e. equivalent to applying the unitary operation
+        this CliffordTableau then applying the second one.
+        """
+        m1 = self.matrix().astype(int)
+        m2 = second.matrix().astype(int)
+
+        p1 = self.rs.astype(int)
+        p2 = second.rs.astype(int)
+
+        merged_m = np.mod(m1.dot(m2), 2)
+        phase = np.mod(p1 + m1.dot(p2), 2)
+
+        # we need more phase correction for expanding Y to XZ and swapping Z_iX_i order.
+        for k in range(2 * self.n):
+            swap_phase = 0  # value betwen 0 and 3 representing [1, i, -1, -i] respectively.
+            prev_row_sum = np.zeros([2 * self.n])
+            swap_phase += np.sum(m1[k, : self.n] * m1[k, self.n :])  # Y gate => iXZ
+            for i, v in enumerate(m1[k]):
+                if v == 0:
+                    continue
+                swap_phase += np.sum(m2[i, : self.n] * m2[i, self.n :])  # Y gate => iXZ
+                swap_phase += 2 * np.sum(m2[i, : self.n] * prev_row_sum[self.n :])  # swap Z_i X_i
+                prev_row_sum += m2[i]
+            prev_row_sum = np.mod(prev_row_sum, 2)
+            swap_phase -= np.sum(prev_row_sum[: self.n] * prev_row_sum[self.n :])  # XZ => -iY
+            phase[k] = (phase[k] + (swap_phase % 4) / 2) % 2
+
+        merged_tableau = CliffordTableau(num_qubits=self.n)
+        merged_tableau.xs = merged_m[:, : self.n]
+        merged_tableau.zs = merged_m[:, self.n :]
+        merged_tableau.rs = phase
+        return merged_tableau
+
+    def __matmul__(self, second: 'CliffordTableau'):
+        if not isinstance(second, CliffordTableau):
+            return NotImplemented
+        return second.then(self)
+
     def _rowsum(self, q1, q2):
         """Implements the "rowsum" routine defined by
         Aaronson and Gottesman.
@@ -202,7 +247,7 @@ def g(x1, z1, x2, z2):
         self._xs[q1, :] ^= self._xs[q2, :]
         self._zs[q1, :] ^= self._zs[q2, :]
 
-    def _row_to_dense_pauli(self, i: int) -> DensePauliString:
+    def _row_to_dense_pauli(self, i: int) -> 'cirq.DensePauliString':
         """
         Args:
             i: index of the row in the tableau.
@@ -226,12 +271,12 @@ def _row_to_dense_pauli(self, i: int) -> DensePauliString:
                 pauli_mask += "I"
         return cirq.DensePauliString(pauli_mask, coefficient=coefficient)
 
-    def stabilizers(self) -> List[DensePauliString]:
+    def stabilizers(self) -> List['cirq.DensePauliString']:
         """Returns the stabilizer generators of the state. These
         are n operators {S_1,S_2,...,S_n} such that S_i |psi> = |psi>"""
         return [self._row_to_dense_pauli(i) for i in range(self.n, 2 * self.n)]
 
-    def destabilizers(self) -> List[DensePauliString]:
+    def destabilizers(self) -> List['cirq.DensePauliString']:
         """Returns the destabilizer generators of the state. These
         are n operators {S_1,S_2,...,S_n} such that along with the stabilizer
         generators above generate the full Pauli group on n qubits."""

cirq-core/cirq/qis/clifford_tableau_test.py

@@ -322,3 +322,115 @@ def test_deprecated_clifford_location():
         from cirq.sim import CliffordTableau
 
         CliffordTableau(num_qubits=1)
+
+
+def _three_identical_table(num_qubits):
+    t1 = cirq.CliffordTableau(num_qubits)
+    t2 = cirq.CliffordTableau(num_qubits)
+    t3 = cirq.CliffordTableau(num_qubits)
+    return t1, t2, t3
+
+
+def test_tableau_then():
+
+    t1, t2, expected_t = _three_identical_table(1)
+    assert expected_t == t1.then(t2)
+
+    t1, t2, expected_t = _three_identical_table(1)
+    _ = [_H(t, 0) for t in (t1, expected_t)]
+    _ = [_H(t, 0) for t in (t2, expected_t)]
+    assert expected_t == t1.then(t2)
+
+    t1, t2, expected_t = _three_identical_table(1)
+    _ = [_X(t, 0) for t in (t1, expected_t)]
+    _ = [_S(t, 0) for t in (t2, expected_t)]
+    assert expected_t == t1.then(t2)
+    assert expected_t != t2.then(t1)
+
+    t1, t2, expected_t = _three_identical_table(1)
+    _ = [_X(t, 0) for t in (t1, expected_t)]
+    _ = [_H(t, 0) for t in (t1, expected_t)]
+    _ = [_Z(t, 0) for t in (t1, expected_t)]
+    _ = [_S(t, 0) for t in (t2, expected_t)]
+    _ = [_H(t, 0) for t in (t2, expected_t)]
+    assert expected_t == t1.then(t2)
+    assert expected_t != t2.then(t1)
+
+    t1, t2, expected_t = _three_identical_table(2)
+    _ = [_H(t, 0) for t in (t1, expected_t)]
+    _ = [_H(t, 1) for t in (t1, expected_t)]
+    _ = [_H(t, 0) for t in (t2, expected_t)]
+    _ = [_H(t, 1) for t in (t2, expected_t)]
+    assert expected_t == t1.then(t2)
+
+    t1, t2, expected_t = _three_identical_table(2)
+    _ = [_H(t, 0) for t in (t1, expected_t)]
+    _ = [_CNOT(t, 0, 1) for t in (t1, expected_t)]
+    _ = [_S(t, 0) for t in (t2, expected_t)]
+    _ = [_X(t, 1) for t in (t2, expected_t)]
+    assert expected_t == t1.then(t2)
+    assert expected_t != t2.then(t1)
+
+    t1, t2, expected_t = _three_identical_table(2)
+    _ = [_H(t, 0) for t in (t1, expected_t)]
+    _ = [_CNOT(t, 0, 1) for t in (t1, expected_t)]
+    _ = [_S(t, 1) for t in (t2, expected_t)]
+    _ = [_CNOT(t, 1, 0) for t in (t2, expected_t)]
+    assert expected_t == t1.then(t2)
+    assert expected_t != t2.then(t1)
+
+    def random_circuit(num_ops, num_qubits, seed=12345):
+        prng = np.random.RandomState(seed)
+        candidate_op = [_H, _S, _X, _Z]
+        if num_qubits > 1:
+            candidate_op = [_H, _S, _X, _Z, _CNOT]
+
+        seq_op = []
+        for _ in range(num_ops):
+            op = prng.randint(len(candidate_op))
+            if op != 4:
+                args = (prng.randint(num_qubits),)
+            else:
+                args = prng.choice(num_qubits, 2, replace=False)
+            seq_op.append((candidate_op[op], args))
+        return seq_op
+
+    # Do small random circuits test 100 times.
+    for seed in range(100):
+        t1, t2, expected_t = _three_identical_table(8)
+        seq_op = random_circuit(num_ops=20, num_qubits=8, seed=seed)
+        for i, (op, args) in enumerate(seq_op):
+            if i < 7:
+                _ = [op(t, *args) for t in (t1, expected_t)]
+            else:
+                _ = [op(t, *args) for t in (t2, expected_t)]
+        assert expected_t == t1.then(t2)
+
+    # Since merging Clifford Tableau operation is O(n^3),
+    # running 100 qubits case is still fast.
+    t1, t2, expected_t = _three_identical_table(100)
+    seq_op = random_circuit(num_ops=1000, num_qubits=100)
+    for i, (op, args) in enumerate(seq_op):
+        if i < 350:
+            _ = [op(t, *args) for t in (t1, expected_t)]
+        else:
+            _ = [op(t, *args) for t in (t2, expected_t)]
+    assert expected_t == t1.then(t2)
+
+
+def test_tableau_matmul():
+    t1, t2, expected_t = _three_identical_table(1)
+    _ = [_H(t, 0) for t in (t1, expected_t)]
+    _ = [_H(t, 0) for t in (t2, expected_t)]
+    assert expected_t == t2 @ t1
+
+    t1, t2, expected_t = _three_identical_table(1)
+    _ = [_X(t, 0) for t in (t1, expected_t)]
+    _ = [_S(t, 0) for t in (t2, expected_t)]
+    assert expected_t == t2 @ t1
+    assert expected_t != t1 @ t2
+
+    with pytest.raises(TypeError):
+        # pylint: disable=pointless-statement
+        t1 @ 21
+        # pylint: enable=pointless-statement