Fix act-on specialization for all the Clifford simulators (quantumlib#4748)

Fixes quantumlib#4639, ref [tinyurl.com/clifford-simulators-refactor](https://tinyurl.com/clifford-simulators-refactor)

This change makes the Clifford simulators more consistent with other simulators. Rather than the gate having all the logic of how to update the simulator state, which seemed out-of-place and limited the ability to reuse in any new kinds of Clifford simulators, this PR moves all the logic to the simulators themselves.

It also creates an ABC for Clifford simulator states, allowing reuse of the evolution logic in the base class, such that the subclasses just have to implement the specific gate applicators.

Note I *tried* doing this via a new Clifford protocol as well ([link](https://github.com/quantumlib/Cirq/compare/master...daxfohl:clifford3?expand=1)) as we originally discussed, and it works but I didn't like the result. The protocol ended up just duplicating all the information of the gate, in a weird structure, and felt like an extra, unnecessary, hacky wrapper around just using the gate itself. So I prefer the approach in this PR that just uses the gate instances directly, though I'm open to suggestion.

@tanujkhattar @viathor (cc @ybc1991)

Author: daxfohl

cirq-core/cirq/__init__.py

@@ -414,6 +414,7 @@
     ActOnCliffordTableauArgs,
     ActOnDensityMatrixArgs,
     ActOnStabilizerCHFormArgs,
+    ActOnStabilizerArgs,
     ActOnStateVectorArgs,
     StabilizerStateChForm,
     CIRCUIT_LIKE,

cirq-core/cirq/ops/common_channels.py

@@ -316,16 +316,6 @@ def __str__(self) -> str:
             return f"depolarize(p={self._p})"
         return f"depolarize(p={self._p},n_qubits={self._n_qubits})"
 
-    def _act_on_(self, args: 'cirq.ActOnArgs', qubits: Sequence['cirq.Qid']) -> bool:
-        from cirq.sim import clifford
-
-        if isinstance(args, clifford.ActOnCliffordTableauArgs):
-            if args.prng.random() < self._p:
-                gate = args.prng.choice([pauli_gates.X, pauli_gates.Y, pauli_gates.Z])
-                protocols.act_on(gate, args, qubits)
-            return True
-        return NotImplemented
-
     def _circuit_diagram_info_(self, args: 'protocols.CircuitDiagramInfoArgs') -> Tuple[str, ...]:
         result: Tuple[str, ...]
         if args.precision is not None:

cirq-core/cirq/ops/common_gates.py

@@ -62,12 +62,6 @@
 """
 
 
-def _act_with_gates(args, qubits, *gates: 'cirq.SupportsActOnQubits') -> None:
-    """Act on the given args with the given gates in order."""
-    for gate in gates:
-        assert gate._act_on_(args, qubits)
-
-
 def _pi(rads):
     return sympy.pi if protocols.is_parameterized(rads) else np.pi
 
@@ -108,35 +102,6 @@ def _apply_unitary_(self, args: 'protocols.ApplyUnitaryArgs') -> Optional[np.nda
             args.available_buffer *= p
         return args.available_buffer
 
-    def _act_on_(self, args: 'cirq.ActOnArgs', qubits: Sequence['cirq.Qid']):
-        from cirq.sim import clifford
-
-        if isinstance(args, clifford.ActOnCliffordTableauArgs):
-            if not protocols.has_stabilizer_effect(self):
-                return NotImplemented
-            tableau = args.tableau
-            q = args.qubit_map[qubits[0]]
-            effective_exponent = self._exponent % 2
-            if effective_exponent == 0.5:
-                tableau.xs[:, q] ^= tableau.zs[:, q]
-                tableau.rs[:] ^= tableau.xs[:, q] & tableau.zs[:, q]
-            elif effective_exponent == 1:
-                tableau.rs[:] ^= tableau.zs[:, q]
-            elif effective_exponent == 1.5:
-                tableau.rs[:] ^= tableau.xs[:, q] & tableau.zs[:, q]
-                tableau.xs[:, q] ^= tableau.zs[:, q]
-            return True
-
-        if isinstance(args, clifford.ActOnStabilizerCHFormArgs):
-            if not protocols.has_stabilizer_effect(self):
-                return NotImplemented
-            _act_with_gates(args, qubits, H, ZPowGate(exponent=self._exponent), H)
-            # Adjust the global phase based on the global_shift parameter.
-            args.state.omega *= np.exp(1j * np.pi * self.global_shift * self.exponent)
-            return True
-
-        return NotImplemented
-
     def in_su2(self) -> 'Rx':
         """Returns an equal-up-global-phase gate from the group SU2."""
         return Rx(rads=self._exponent * _pi(self._exponent))
@@ -362,51 +327,6 @@ def _apply_unitary_(self, args: 'protocols.ApplyUnitaryArgs') -> Optional[np.nda
             args.available_buffer *= p
         return args.available_buffer
 
-    def _act_on_(self, args: 'cirq.ActOnArgs', qubits: Sequence['cirq.Qid']):
-        from cirq.sim import clifford
-
-        if isinstance(args, clifford.ActOnCliffordTableauArgs):
-            if not protocols.has_stabilizer_effect(self):
-                return NotImplemented
-            tableau = args.tableau
-            q = args.qubit_map[qubits[0]]
-            effective_exponent = self._exponent % 2
-            if effective_exponent == 0.5:
-                tableau.rs[:] ^= tableau.xs[:, q] & (~tableau.zs[:, q])
-                (tableau.xs[:, q], tableau.zs[:, q]) = (
-                    tableau.zs[:, q].copy(),
-                    tableau.xs[:, q].copy(),
-                )
-            elif effective_exponent == 1:
-                tableau.rs[:] ^= tableau.xs[:, q] ^ tableau.zs[:, q]
-            elif effective_exponent == 1.5:
-                tableau.rs[:] ^= ~(tableau.xs[:, q]) & tableau.zs[:, q]
-                (tableau.xs[:, q], tableau.zs[:, q]) = (
-                    tableau.zs[:, q].copy(),
-                    tableau.xs[:, q].copy(),
-                )
-            return True
-
-        if isinstance(args, clifford.ActOnStabilizerCHFormArgs):
-            if not protocols.has_stabilizer_effect(self):
-                return NotImplemented
-            effective_exponent = self._exponent % 2
-            state = args.state
-            Z = ZPowGate()
-            if effective_exponent == 0.5:
-                _act_with_gates(args, qubits, Z, H)
-                state.omega *= (1 + 1j) / (2 ** 0.5)
-            elif effective_exponent == 1:
-                _act_with_gates(args, qubits, Z, H, Z, H)
-                state.omega *= 1j
-            elif effective_exponent == 1.5:
-                _act_with_gates(args, qubits, H, Z)
-                state.omega *= (1 - 1j) / (2 ** 0.5)
-            # Adjust the global phase based on the global_shift parameter.
-            args.state.omega *= np.exp(1j * np.pi * self.global_shift * self.exponent)
-            return True
-        return NotImplemented
-
     def in_su2(self) -> 'Ry':
         """Returns an equal-up-global-phase gate from the group SU2."""
         return Ry(rads=self._exponent * _pi(self._exponent))
@@ -580,42 +500,6 @@ def _apply_unitary_(self, args: 'protocols.ApplyUnitaryArgs') -> Optional[np.nda
             args.target_tensor *= p
         return args.target_tensor
 
-    def _act_on_(self, args: 'cirq.ActOnArgs', qubits: Sequence['cirq.Qid']):
-        from cirq.sim import clifford
-
-        if isinstance(args, clifford.ActOnCliffordTableauArgs):
-            if not protocols.has_stabilizer_effect(self):
-                return NotImplemented
-            tableau = args.tableau
-            q = args.qubit_map[qubits[0]]
-            effective_exponent = self._exponent % 2
-            if effective_exponent == 0.5:
-                tableau.rs[:] ^= tableau.xs[:, q] & tableau.zs[:, q]
-                tableau.zs[:, q] ^= tableau.xs[:, q]
-            elif effective_exponent == 1:
-                tableau.rs[:] ^= tableau.xs[:, q]
-            elif effective_exponent == 1.5:
-                tableau.rs[:] ^= tableau.xs[:, q] & (~tableau.zs[:, q])
-                tableau.zs[:, q] ^= tableau.xs[:, q]
-            return True
-
-        if isinstance(args, clifford.ActOnStabilizerCHFormArgs):
-            if not protocols.has_stabilizer_effect(self):
-                return NotImplemented
-            q = args.qubit_map[qubits[0]]
-            effective_exponent = self._exponent % 2
-            state = args.state
-            for _ in range(int(effective_exponent * 2)):
-                # Prescription for S left multiplication.
-                # Reference: https://arxiv.org/abs/1808.00128 Proposition 4 end
-                state.M[q, :] ^= state.G[q, :]
-                state.gamma[q] = (state.gamma[q] - 1) % 4
-            # Adjust the global phase based on the global_shift parameter.
-            args.state.omega *= np.exp(1j * np.pi * self.global_shift * self.exponent)
-            return True
-
-        return NotImplemented
-
     def _decompose_into_clifford_with_qubits_(self, qubits):
         from cirq.ops.clifford_gate import SingleQubitCliffordGate
 
@@ -899,49 +783,6 @@ def _apply_unitary_(self, args: 'protocols.ApplyUnitaryArgs') -> Optional[np.nda
         args.target_tensor *= np.sqrt(2) * p
         return args.target_tensor
 
-    def _act_on_(self, args: 'cirq.ActOnArgs', qubits: Sequence['cirq.Qid']):
-        from cirq.sim import clifford
-
-        if isinstance(args, clifford.ActOnCliffordTableauArgs):
-            if not protocols.has_stabilizer_effect(self):
-                return NotImplemented
-            tableau = args.tableau
-            q = args.qubit_map[qubits[0]]
-            if self._exponent % 2 == 1:
-                (tableau.xs[:, q], tableau.zs[:, q]) = (
-                    tableau.zs[:, q].copy(),
-                    tableau.xs[:, q].copy(),
-                )
-                tableau.rs[:] ^= tableau.xs[:, q] & tableau.zs[:, q]
-            return True
-
-        if isinstance(args, clifford.ActOnStabilizerCHFormArgs):
-            if not protocols.has_stabilizer_effect(self):
-                return NotImplemented
-            q = args.qubit_map[qubits[0]]
-            state = args.state
-            if self._exponent % 2 == 1:
-                # Prescription for H left multiplication
-                # Reference: https://arxiv.org/abs/1808.00128
-                # Equations 48, 49 and Proposition 4
-                t = state.s ^ (state.G[q, :] & state.v)
-                u = state.s ^ (state.F[q, :] & (~state.v)) ^ (state.M[q, :] & state.v)
-
-                alpha = sum(state.G[q, :] & (~state.v) & state.s) % 2
-                beta = sum(state.M[q, :] & (~state.v) & state.s)
-                beta += sum(state.F[q, :] & state.v & state.M[q, :])
-                beta += sum(state.F[q, :] & state.v & state.s)
-                beta %= 2
-
-                delta = (state.gamma[q] + 2 * (alpha + beta)) % 4
-
-                state.update_sum(t, u, delta=delta, alpha=alpha)
-            # Adjust the global phase based on the global_shift parameter.
-            args.state.omega *= np.exp(1j * np.pi * self.global_shift * self.exponent)
-            return True
-
-        return NotImplemented
-
     def _decompose_(self, qubits):
         q = qubits[0]
 
@@ -1063,52 +904,6 @@ def _apply_unitary_(
             args.target_tensor *= p
         return args.target_tensor
 
-    def _act_on_(self, args: 'cirq.ActOnArgs', qubits: Sequence['cirq.Qid']):
-        from cirq.sim import clifford
-
-        if isinstance(args, clifford.ActOnCliffordTableauArgs):
-            if not protocols.has_stabilizer_effect(self):
-                return NotImplemented
-            tableau = args.tableau
-            q1 = args.qubit_map[qubits[0]]
-            q2 = args.qubit_map[qubits[1]]
-            if self._exponent % 2 == 1:
-                (tableau.xs[:, q2], tableau.zs[:, q2]) = (
-                    tableau.zs[:, q2].copy(),
-                    tableau.xs[:, q2].copy(),
-                )
-                tableau.rs[:] ^= tableau.xs[:, q2] & tableau.zs[:, q2]
-                tableau.rs[:] ^= (
-                    tableau.xs[:, q1]
-                    & tableau.zs[:, q2]
-                    & (~(tableau.xs[:, q2] ^ tableau.zs[:, q1]))
-                )
-                tableau.xs[:, q2] ^= tableau.xs[:, q1]
-                tableau.zs[:, q1] ^= tableau.zs[:, q2]
-                (tableau.xs[:, q2], tableau.zs[:, q2]) = (
-                    tableau.zs[:, q2].copy(),
-                    tableau.xs[:, q2].copy(),
-                )
-                tableau.rs[:] ^= tableau.xs[:, q2] & tableau.zs[:, q2]
-            return True
-
-        if isinstance(args, clifford.ActOnStabilizerCHFormArgs):
-            if not protocols.has_stabilizer_effect(self):
-                return NotImplemented
-            q1 = args.qubit_map[qubits[0]]
-            q2 = args.qubit_map[qubits[1]]
-            state = args.state
-            if self._exponent % 2 == 1:
-                # Prescription for CZ left multiplication.
-                # Reference: https://arxiv.org/abs/1808.00128 Proposition 4 end
-                state.M[q1, :] ^= state.G[q2, :]
-                state.M[q2, :] ^= state.G[q1, :]
-            # Adjust the global phase based on the global_shift parameter.
-            args.state.omega *= np.exp(1j * np.pi * self.global_shift * self.exponent)
-            return True
-
-        return NotImplemented
-
     def _pauli_expansion_(self) -> value.LinearDict[str]:
         if protocols.is_parameterized(self):
             return NotImplemented
@@ -1291,48 +1086,6 @@ def _apply_unitary_(self, args: 'protocols.ApplyUnitaryArgs') -> Optional[np.nda
             args.target_tensor *= p
         return args.target_tensor
 
-    def _act_on_(self, args: 'cirq.ActOnArgs', qubits: Sequence['cirq.Qid']):
-        from cirq.sim import clifford
-
-        if isinstance(args, clifford.ActOnCliffordTableauArgs):
-            if not protocols.has_stabilizer_effect(self):
-                return NotImplemented
-            tableau = args.tableau
-            q1 = args.qubit_map[qubits[0]]
-            q2 = args.qubit_map[qubits[1]]
-            if self._exponent % 2 == 1:
-                tableau.rs[:] ^= (
-                    tableau.xs[:, q1]
-                    & tableau.zs[:, q2]
-                    & (~(tableau.xs[:, q2] ^ tableau.zs[:, q1]))
-                )
-                tableau.xs[:, q2] ^= tableau.xs[:, q1]
-                tableau.zs[:, q1] ^= tableau.zs[:, q2]
-            return True
-
-        if isinstance(args, clifford.ActOnStabilizerCHFormArgs):
-            if not protocols.has_stabilizer_effect(self):
-                return NotImplemented
-            q1 = args.qubit_map[qubits[0]]
-            q2 = args.qubit_map[qubits[1]]
-            state = args.state
-            if self._exponent % 2 == 1:
-                # Prescription for CX left multiplication.
-                # Reference: https://arxiv.org/abs/1808.00128 Proposition 4 end
-                state.gamma[q1] = (
-                    state.gamma[q1]
-                    + state.gamma[q2]
-                    + 2 * (sum(state.M[q1, :] & state.F[q2, :]) % 2)
-                ) % 4
-                state.G[q2, :] ^= state.G[q1, :]
-                state.F[q1, :] ^= state.F[q2, :]
-                state.M[q1, :] ^= state.M[q2, :]
-            # Adjust the global phase based on the global_shift parameter.
-            args.state.omega *= np.exp(1j * np.pi * self.global_shift * self.exponent)
-            return True
-
-        return NotImplemented
-
     def _pauli_expansion_(self) -> value.LinearDict[str]:
         if protocols.is_parameterized(self):
             return NotImplemented

cirq-core/cirq/ops/global_phase_op.py

@@ -87,20 +87,6 @@ def _apply_unitary_(self, args) -> np.ndarray:
     def _has_stabilizer_effect_(self) -> bool:
         return True
 
-    def _act_on_(self, args: 'cirq.ActOnArgs', qubits):
-        from cirq.sim import clifford
-
-        if isinstance(args, clifford.ActOnCliffordTableauArgs):
-            # Since CliffordTableau does not keep track of the global phase,
-            # it's safe to just ignore it here.
-            return True
-
-        if isinstance(args, clifford.ActOnStabilizerCHFormArgs):
-            args.state.omega *= self.coefficient
-            return True
-
-        return NotImplemented
-
     def __str__(self) -> str:
         return str(self.coefficient)
 

cirq-core/cirq/ops/random_gate_channel.py

@@ -22,7 +22,6 @@
     cast,
     SupportsFloat,
     Optional,
-    Sequence,
 )
 
 import numpy as np
@@ -123,20 +122,6 @@ def _trace_distance_bound_(self) -> float:
             result *= float(self.probability)
         return result
 
-    def _act_on_(self, args: 'cirq.ActOnArgs', qubits: Sequence['cirq.Qid']):
-        from cirq.sim import clifford
-
-        if self._is_parameterized_():
-            return NotImplemented
-        if isinstance(args, clifford.ActOnCliffordTableauArgs):
-            if args.prng.random() < self.probability:
-                # Note: because we're doing this probabilistically, it's not
-                # safe to fallback to other strategies if act_on fails. Those
-                # strategies could double-count the probability.
-                protocols.act_on(self.sub_gate, args, qubits)
-            return True
-        return NotImplemented
-
     def _json_dict_(self) -> Dict[str, Any]:
         return protocols.obj_to_dict_helper(self, ['sub_gate', 'probability'])