From b5af9da79ff98c650260e4efc2de27f96eb464da Mon Sep 17 00:00:00 2001 From: Dax Fohl Date: Wed, 19 Jan 2022 09:30:07 -0800 Subject: [PATCH] Fix act-on specialization for all the Clifford simulators (#4748) Fixes #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) --- cirq-core/cirq/__init__.py | 1 + cirq-core/cirq/ops/common_channels.py | 10 - cirq-core/cirq/ops/common_gates.py | 247 ------------------ cirq-core/cirq/ops/global_phase_op.py | 14 - cirq-core/cirq/ops/random_gate_channel.py | 15 -- cirq-core/cirq/ops/swap_gates.py | 21 +- cirq-core/cirq/sim/__init__.py | 1 + cirq-core/cirq/sim/clifford/__init__.py | 4 + .../clifford/act_on_clifford_tableau_args.py | 165 ++++++++---- .../sim/clifford/act_on_stabilizer_args.py | 165 ++++++++++++ .../act_on_stabilizer_ch_form_args.py | 172 +++++++----- 11 files changed, 392 insertions(+), 423 deletions(-) create mode 100644 cirq-core/cirq/sim/clifford/act_on_stabilizer_args.py diff --git a/cirq-core/cirq/__init__.py b/cirq-core/cirq/__init__.py index 2e9d2f54cf0..d35ee90a00a 100644 --- a/cirq-core/cirq/__init__.py +++ b/cirq-core/cirq/__init__.py @@ -414,6 +414,7 @@ ActOnCliffordTableauArgs, ActOnDensityMatrixArgs, ActOnStabilizerCHFormArgs, + ActOnStabilizerArgs, ActOnStateVectorArgs, StabilizerStateChForm, CIRCUIT_LIKE, diff --git a/cirq-core/cirq/ops/common_channels.py b/cirq-core/cirq/ops/common_channels.py index 936e19f4faf..753db55dc78 100644 --- a/cirq-core/cirq/ops/common_channels.py +++ b/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: diff --git a/cirq-core/cirq/ops/common_gates.py b/cirq-core/cirq/ops/common_gates.py index 4bef5936e90..75938f9a385 100644 --- a/cirq-core/cirq/ops/common_gates.py +++ b/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 diff --git a/cirq-core/cirq/ops/global_phase_op.py b/cirq-core/cirq/ops/global_phase_op.py index 5588172e12e..139f51ab177 100644 --- a/cirq-core/cirq/ops/global_phase_op.py +++ b/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) diff --git a/cirq-core/cirq/ops/random_gate_channel.py b/cirq-core/cirq/ops/random_gate_channel.py index fbb0a6d66ae..21689b297dc 100644 --- a/cirq-core/cirq/ops/random_gate_channel.py +++ b/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']) diff --git a/cirq-core/cirq/ops/swap_gates.py b/cirq-core/cirq/ops/swap_gates.py index fd464d6cd6b..d0b556a23ca 100644 --- a/cirq-core/cirq/ops/swap_gates.py +++ b/cirq-core/cirq/ops/swap_gates.py @@ -24,7 +24,7 @@ EigenGate. """ -from typing import Optional, Tuple, TYPE_CHECKING, List, Sequence +from typing import Optional, Tuple, TYPE_CHECKING, List import numpy as np import sympy @@ -94,25 +94,6 @@ def _has_stabilizer_effect_(self) -> Optional[bool]: return None return self.exponent % 1 == 0 - def _act_on_(self, args: 'cirq.ActOnArgs', qubits: Sequence['cirq.Qid']): - from cirq import ops, sim, protocols - - if isinstance(args, (sim.ActOnStabilizerCHFormArgs, sim.ActOnCliffordTableauArgs)): - if not self._has_stabilizer_effect_(): - return NotImplemented - if isinstance(args, sim.ActOnStabilizerCHFormArgs): - args.state.omega *= 1j ** (2 * self.global_shift * self._exponent) - - if self._exponent % 2 == 1: - protocols.act_on(ops.CNOT, args, qubits) - protocols.act_on(ops.CNOT, args, tuple(reversed(qubits))) - protocols.act_on(ops.CNOT, args, qubits) - - # An even exponent does not change anything except the global phase above. - return True - - return NotImplemented - def _apply_unitary_(self, args: 'protocols.ApplyUnitaryArgs') -> Optional[np.ndarray]: if self._exponent != 1: return NotImplemented diff --git a/cirq-core/cirq/sim/__init__.py b/cirq-core/cirq/sim/__init__.py index e09f43eae59..2394f785b18 100644 --- a/cirq-core/cirq/sim/__init__.py +++ b/cirq-core/cirq/sim/__init__.py @@ -90,6 +90,7 @@ from cirq.sim.clifford import ( ActOnCliffordTableauArgs, ActOnStabilizerCHFormArgs, + ActOnStabilizerArgs, StabilizerSampler, StabilizerStateChForm, CliffordSimulator, diff --git a/cirq-core/cirq/sim/clifford/__init__.py b/cirq-core/cirq/sim/clifford/__init__.py index 2f8a2e5a310..42ee7bf7ad8 100644 --- a/cirq-core/cirq/sim/clifford/__init__.py +++ b/cirq-core/cirq/sim/clifford/__init__.py @@ -7,6 +7,10 @@ ActOnStabilizerCHFormArgs, ) +from cirq.sim.clifford.act_on_stabilizer_args import ( + ActOnStabilizerArgs, +) + from cirq.sim.clifford.stabilizer_state_ch_form import ( StabilizerStateChForm, ) diff --git a/cirq-core/cirq/sim/clifford/act_on_clifford_tableau_args.py b/cirq-core/cirq/sim/clifford/act_on_clifford_tableau_args.py index 9655f17411f..9629e7a7a9b 100644 --- a/cirq-core/cirq/sim/clifford/act_on_clifford_tableau_args.py +++ b/cirq-core/cirq/sim/clifford/act_on_clifford_tableau_args.py @@ -14,27 +14,20 @@ """A protocol for implementing high performance clifford tableau evolutions for Clifford Simulator.""" -from typing import Any, Dict, TYPE_CHECKING, List, Sequence, Union +from typing import Any, Dict, TYPE_CHECKING, List, Sequence import numpy as np from cirq.ops import common_gates -from cirq.ops import pauli_gates -from cirq.ops.clifford_gate import SingleQubitCliffordGate -from cirq.protocols import has_unitary, num_qubits, unitary -from cirq.sim.act_on_args import ActOnArgs -from cirq.type_workarounds import NotImplementedType +from cirq.ops import global_phase_op +from cirq.sim.clifford.act_on_stabilizer_args import ActOnStabilizerArgs if TYPE_CHECKING: import cirq -class ActOnCliffordTableauArgs(ActOnArgs): - """State and context for an operation acting on a clifford tableau. - - To act on this object, directly edit the `tableau` property, which is - storing the density matrix of the quantum system with one axis per qubit. - """ +class ActOnCliffordTableauArgs(ActOnStabilizerArgs): + """State and context for an operation acting on a clifford tableau.""" def __init__( self, @@ -59,25 +52,6 @@ def __init__( super().__init__(prng, qubits, log_of_measurement_results) self.tableau = tableau - def _act_on_fallback_( - self, - action: Union['cirq.Operation', 'cirq.Gate'], - qubits: Sequence['cirq.Qid'], - allow_decompose: bool = True, - ) -> Union[bool, NotImplementedType]: - strats = [] - if allow_decompose: - strats.append(_strat_act_on_clifford_tableau_from_single_qubit_decompose) - for strat in strats: - result = strat(action, self, qubits) - if result is False: - break # coverage: ignore - if result is True: - return True - assert result is NotImplemented, str(result) - - return NotImplemented - def _perform_measurement(self, qubits: Sequence['cirq.Qid']) -> List[int]: """Returns the measurement from the tableau.""" return [self.tableau._measure(self.qubit_map[q], self.prng) for q in qubits] @@ -94,24 +68,111 @@ def sample( # Unnecessary for now but can be added later if there is a use case. raise NotImplementedError() - -def _strat_act_on_clifford_tableau_from_single_qubit_decompose( - val: Any, args: 'cirq.ActOnCliffordTableauArgs', qubits: Sequence['cirq.Qid'] -) -> bool: - if num_qubits(val) == 1: - if not has_unitary(val): - return NotImplemented - u = unitary(val) - clifford_gate = SingleQubitCliffordGate.from_unitary(u) - if clifford_gate is not None: - for axis, quarter_turns in clifford_gate.decompose_rotation(): - if axis == pauli_gates.X: - common_gates.XPowGate(exponent=quarter_turns / 2)._act_on_(args, qubits) - elif axis == pauli_gates.Y: - common_gates.YPowGate(exponent=quarter_turns / 2)._act_on_(args, qubits) - else: - assert axis == pauli_gates.Z - common_gates.ZPowGate(exponent=quarter_turns / 2)._act_on_(args, qubits) - return True - - return NotImplemented + def _x(self, g: common_gates.XPowGate, axis: int): + exponent = g.exponent + if exponent % 2 == 0: + return + if exponent % 0.5 != 0.0: + raise ValueError('X exponent must be half integer') # coverage: ignore + tableau = self.tableau + effective_exponent = exponent % 2 + if effective_exponent == 0.5: + tableau.xs[:, axis] ^= tableau.zs[:, axis] + tableau.rs[:] ^= tableau.xs[:, axis] & tableau.zs[:, axis] + elif effective_exponent == 1: + tableau.rs[:] ^= tableau.zs[:, axis] + elif effective_exponent == 1.5: + tableau.rs[:] ^= tableau.xs[:, axis] & tableau.zs[:, axis] + tableau.xs[:, axis] ^= tableau.zs[:, axis] + + def _y(self, g: common_gates.YPowGate, axis: int): + exponent = g.exponent + if exponent % 2 == 0: + return + if exponent % 0.5 != 0.0: + raise ValueError('Y exponent must be half integer') # coverage: ignore + tableau = self.tableau + effective_exponent = exponent % 2 + if effective_exponent == 0.5: + tableau.rs[:] ^= tableau.xs[:, axis] & (~tableau.zs[:, axis]) + (tableau.xs[:, axis], tableau.zs[:, axis]) = ( + tableau.zs[:, axis].copy(), + tableau.xs[:, axis].copy(), + ) + elif effective_exponent == 1: + tableau.rs[:] ^= tableau.xs[:, axis] ^ tableau.zs[:, axis] + elif effective_exponent == 1.5: + tableau.rs[:] ^= ~(tableau.xs[:, axis]) & tableau.zs[:, axis] + (tableau.xs[:, axis], tableau.zs[:, axis]) = ( + tableau.zs[:, axis].copy(), + tableau.xs[:, axis].copy(), + ) + + def _z(self, g: common_gates.ZPowGate, axis: int): + exponent = g.exponent + if exponent % 2 == 0: + return + if exponent % 0.5 != 0.0: + raise ValueError('Z exponent must be half integer') # coverage: ignore + tableau = self.tableau + effective_exponent = exponent % 2 + if effective_exponent == 0.5: + tableau.rs[:] ^= tableau.xs[:, axis] & tableau.zs[:, axis] + tableau.zs[:, axis] ^= tableau.xs[:, axis] + elif effective_exponent == 1: + tableau.rs[:] ^= tableau.xs[:, axis] + elif effective_exponent == 1.5: + tableau.rs[:] ^= tableau.xs[:, axis] & (~tableau.zs[:, axis]) + tableau.zs[:, axis] ^= tableau.xs[:, axis] + + def _h(self, g: common_gates.HPowGate, axis: int): + exponent = g.exponent + if exponent % 2 == 0: + return + if exponent % 1 != 0: + raise ValueError('H exponent must be integer') # coverage: ignore + self._y(common_gates.YPowGate(exponent=0.5), axis) + self._x(common_gates.XPowGate(), axis) + + def _cz(self, g: common_gates.CZPowGate, control_axis: int, target_axis: int): + exponent = g.exponent + if exponent % 2 == 0: + return + if exponent % 1 != 0: + raise ValueError('CZ exponent must be integer') # coverage: ignore + tableau = self.tableau + (tableau.xs[:, target_axis], tableau.zs[:, target_axis]) = ( + tableau.zs[:, target_axis].copy(), + tableau.xs[:, target_axis].copy(), + ) + tableau.rs[:] ^= tableau.xs[:, target_axis] & tableau.zs[:, target_axis] + tableau.rs[:] ^= ( + tableau.xs[:, control_axis] + & tableau.zs[:, target_axis] + & (~(tableau.xs[:, target_axis] ^ tableau.zs[:, control_axis])) + ) + tableau.xs[:, target_axis] ^= tableau.xs[:, control_axis] + tableau.zs[:, control_axis] ^= tableau.zs[:, target_axis] + (tableau.xs[:, target_axis], tableau.zs[:, target_axis]) = ( + tableau.zs[:, target_axis].copy(), + tableau.xs[:, target_axis].copy(), + ) + tableau.rs[:] ^= tableau.xs[:, target_axis] & tableau.zs[:, target_axis] + + def _cx(self, g: common_gates.CXPowGate, control_axis: int, target_axis: int): + exponent = g.exponent + if exponent % 2 == 0: + return + if exponent % 1 != 0: + raise ValueError('CX exponent must be integer') # coverage: ignore + tableau = self.tableau + tableau.rs[:] ^= ( + tableau.xs[:, control_axis] + & tableau.zs[:, target_axis] + & (~(tableau.xs[:, target_axis] ^ tableau.zs[:, control_axis])) + ) + tableau.xs[:, target_axis] ^= tableau.xs[:, control_axis] + tableau.zs[:, control_axis] ^= tableau.zs[:, target_axis] + + def _global_phase(self, g: global_phase_op.GlobalPhaseGate): + pass diff --git a/cirq-core/cirq/sim/clifford/act_on_stabilizer_args.py b/cirq-core/cirq/sim/clifford/act_on_stabilizer_args.py new file mode 100644 index 00000000000..24be68544e4 --- /dev/null +++ b/cirq-core/cirq/sim/clifford/act_on_stabilizer_args.py @@ -0,0 +1,165 @@ +# Copyright 2021 The Cirq Developers +# +# 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 +# +# https://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 abc +from typing import Any, Sequence, TYPE_CHECKING, Union + +import numpy as np + +from cirq import ops, protocols, linalg +from cirq.ops import common_gates, global_phase_op, matrix_gates, swap_gates +from cirq.ops.clifford_gate import SingleQubitCliffordGate +from cirq.protocols import has_unitary, num_qubits, unitary +from cirq.sim.act_on_args import ActOnArgs +from cirq.type_workarounds import NotImplementedType + +if TYPE_CHECKING: + import cirq + + +class ActOnStabilizerArgs(ActOnArgs, metaclass=abc.ABCMeta): + """Abstract wrapper around a stabilizer state for the act_on protocol.""" + + def _act_on_fallback_( + self, + action: Union['cirq.Operation', 'cirq.Gate'], + qubits: Sequence['cirq.Qid'], + allow_decompose: bool = True, + ) -> Union[bool, NotImplementedType]: + strats = [ + self._strat_apply_gate, + self._strat_apply_mixture, + ] + if allow_decompose: + strats.append(self._strat_decompose) + strats.append(self._strat_act_from_single_qubit_decompose) + for strat in strats: + result = strat(action, qubits) # type: ignore + if result is True: + return True + assert result is NotImplemented, str(result) + + return NotImplemented + + @abc.abstractmethod + def _x(self, g: common_gates.XPowGate, axis: int): + """Apply an X gate""" + + @abc.abstractmethod + def _y(self, g: common_gates.YPowGate, axis: int): + """Apply a Y gate""" + + @abc.abstractmethod + def _z(self, g: common_gates.ZPowGate, axis: int): + """Apply a Z gate""" + + @abc.abstractmethod + def _h(self, g: common_gates.HPowGate, axis: int): + """Apply an H gate""" + + @abc.abstractmethod + def _cz(self, g: common_gates.CZPowGate, control_axis: int, target_axis: int): + """Apply a CZ gate""" + + @abc.abstractmethod + def _cx(self, g: common_gates.CXPowGate, control_axis: int, target_axis: int): + """Apply a CX gate""" + + @abc.abstractmethod + def _global_phase(self, g: global_phase_op.GlobalPhaseGate): + """Apply global phase""" + + def _swap(self, g: swap_gates.SwapPowGate, control_axis: int, target_axis: int): + """Apply a SWAP gate""" + if g.exponent % 1 != 0: + raise ValueError('Swap exponent must be integer') # coverage: ignore + self._cx(common_gates.CX, control_axis, target_axis) + self._cx( + common_gates.CXPowGate(exponent=g.exponent, global_shift=g.global_shift), + target_axis, + control_axis, + ) + self._cx(common_gates.CX, control_axis, target_axis) + + def _strat_apply_gate(self, val: Any, qubits: Sequence['cirq.Qid']) -> bool: + if not protocols.has_stabilizer_effect(val): + return NotImplemented + gate = val.gate if isinstance(val, ops.Operation) else val + axes = self.get_axes(qubits) + if isinstance(gate, common_gates.XPowGate): + self._x(gate, axes[0]) + elif isinstance(gate, common_gates.YPowGate): + self._y(gate, axes[0]) + elif isinstance(gate, common_gates.ZPowGate): + self._z(gate, axes[0]) + elif isinstance(gate, common_gates.HPowGate): + self._h(gate, axes[0]) + elif isinstance(gate, common_gates.CXPowGate): + self._cx(gate, axes[0], axes[1]) + elif isinstance(gate, common_gates.CZPowGate): + self._cz(gate, axes[0], axes[1]) + elif isinstance(gate, global_phase_op.GlobalPhaseGate): + self._global_phase(gate) + elif isinstance(gate, swap_gates.SwapPowGate): + self._swap(gate, axes[0], axes[1]) + else: + return NotImplemented + return True + + def _strat_apply_mixture(self, val: Any, qubits: Sequence['cirq.Qid']) -> bool: + mixture = protocols.mixture(val, None) + if mixture is None: + return NotImplemented + if not all(linalg.is_unitary(m) for _, m in mixture): + return NotImplemented + probabilities, unitaries = zip(*mixture) + index = self.prng.choice(len(unitaries), p=probabilities) + return self._strat_act_from_single_qubit_decompose( + matrix_gates.MatrixGate(unitaries[index]), qubits + ) + + def _strat_act_from_single_qubit_decompose( + self, val: Any, qubits: Sequence['cirq.Qid'] + ) -> bool: + if num_qubits(val) == 1: + if not has_unitary(val): + return NotImplemented + u = unitary(val) + clifford_gate = SingleQubitCliffordGate.from_unitary(u) + if clifford_gate is not None: + # Gather the effective unitary applied so as to correct for the + # global phase later. + final_unitary = np.eye(2) + for axis, quarter_turns in clifford_gate.decompose_rotation(): + gate = axis ** (quarter_turns / 2) + self._strat_apply_gate(gate, qubits) + final_unitary = np.matmul(unitary(gate), final_unitary) + + # Find the entry with the largest magnitude in the input unitary. + k = max(np.ndindex(*u.shape), key=lambda t: abs(u[t])) + # Correct the global phase that wasn't conserved in the above + # decomposition. + self._global_phase(global_phase_op.GlobalPhaseGate(u[k] / final_unitary[k])) + return True + + return NotImplemented + + def _strat_decompose(self, val: Any, qubits: Sequence['cirq.Qid']) -> bool: + gate = val.gate if isinstance(val, ops.Operation) else val + operations = protocols.decompose_once_with_qubits(gate, qubits, None) + if operations is None or not all(protocols.has_stabilizer_effect(op) for op in operations): + return NotImplemented + for op in operations: + protocols.act_on(op, self) + return True diff --git a/cirq-core/cirq/sim/clifford/act_on_stabilizer_ch_form_args.py b/cirq-core/cirq/sim/clifford/act_on_stabilizer_ch_form_args.py index 0525e407bcd..9a479ad8944 100644 --- a/cirq-core/cirq/sim/clifford/act_on_stabilizer_ch_form_args.py +++ b/cirq-core/cirq/sim/clifford/act_on_stabilizer_ch_form_args.py @@ -12,28 +12,21 @@ # See the License for the specific language governing permissions and # limitations under the License. -from typing import Any, Dict, TYPE_CHECKING, List, Sequence, Union +from typing import Any, Dict, TYPE_CHECKING, List, Sequence import numpy as np from cirq import value, ops, protocols from cirq.ops import common_gates, pauli_gates -from cirq.ops.clifford_gate import SingleQubitCliffordGate -from cirq.protocols import has_unitary, num_qubits, unitary -from cirq.sim.act_on_args import ActOnArgs -from cirq.type_workarounds import NotImplementedType +from cirq.ops import global_phase_op +from cirq.sim.clifford.act_on_stabilizer_args import ActOnStabilizerArgs if TYPE_CHECKING: import cirq - from typing import Optional -class ActOnStabilizerCHFormArgs(ActOnArgs): - """Wrapper around a stabilizer state in CH form for the act_on protocol. - - To act on this object, directly edit the `state` property, which is - storing the stabilizer state of the quantum system with one axis per qubit. - """ +class ActOnStabilizerCHFormArgs(ActOnStabilizerArgs): + """Wrapper around a stabilizer state in CH form for the act_on protocol.""" def __init__( self, @@ -43,6 +36,7 @@ def __init__( qubits: Sequence['cirq.Qid'] = None, ): """Initializes with the given state and the axes for the operation. + Args: state: The StabilizerStateChForm to act on. Operations are expected to perform inplace edits of this object. @@ -57,23 +51,6 @@ def __init__( super().__init__(prng, qubits, log_of_measurement_results) self.state = state - def _act_on_fallback_( - self, - action: Union['cirq.Operation', 'cirq.Gate'], - qubits: Sequence['cirq.Qid'], - allow_decompose: bool = True, - ) -> Union[bool, NotImplementedType]: - strats = [] - if allow_decompose: - strats.append(_strat_act_on_stabilizer_ch_form_from_single_qubit_decompose) - for strat in strats: - result = strat(action, self, qubits) - if result is True: - return True - assert result is NotImplemented, str(result) - - return NotImplemented - def _perform_measurement(self, qubits: Sequence['cirq.Qid']) -> List[int]: """Returns the measurement from the stabilizer state form.""" return [self.state._measure(self.qubit_map[q], self.prng) for q in qubits] @@ -96,39 +73,104 @@ def sample( protocols.act_on(op, ch_form_args) return np.array(list(measurements.values()), dtype=bool) - -def _strat_act_on_stabilizer_ch_form_from_single_qubit_decompose( - val: Any, args: 'cirq.ActOnStabilizerCHFormArgs', qubits: Sequence['cirq.Qid'] -) -> bool: - if num_qubits(val) == 1: - if not has_unitary(val): - return NotImplemented - u = unitary(val) - clifford_gate = SingleQubitCliffordGate.from_unitary(u) - if clifford_gate is not None: - # Gather the effective unitary applied so as to correct for the - # global phase later. - final_unitary = np.eye(2) - for axis, quarter_turns in clifford_gate.decompose_rotation(): - gate: Optional[cirq.Gate] = None - if axis == pauli_gates.X: - gate = common_gates.XPowGate(exponent=quarter_turns / 2) - assert gate._act_on_(args, qubits) - elif axis == pauli_gates.Y: - gate = common_gates.YPowGate(exponent=quarter_turns / 2) - assert gate._act_on_(args, qubits) - else: - assert axis == pauli_gates.Z - gate = common_gates.ZPowGate(exponent=quarter_turns / 2) - assert gate._act_on_(args, qubits) - - final_unitary = np.matmul(unitary(gate), final_unitary) - - # Find the entry with the largest magnitude in the input unitary. - k = max(np.ndindex(*u.shape), key=lambda t: abs(u[t])) - # Correct the global phase that wasn't conserved in the above - # decomposition. - args.state.omega *= u[k] / final_unitary[k] - return True - - return NotImplemented + def _x(self, g: common_gates.XPowGate, axis: int): + exponent = g.exponent + if exponent % 2 != 0: + if exponent % 0.5 != 0.0: + raise ValueError('Y exponent must be half integer') # coverage: ignore + self._h(common_gates.H, axis) + self._z(common_gates.ZPowGate(exponent=exponent), axis) + self._h(common_gates.H, axis) + self.state.omega *= _phase(g) + + def _y(self, g: common_gates.YPowGate, axis: int): + exponent = g.exponent + if exponent % 0.5 != 0.0: + raise ValueError('Y exponent must be half integer') # coverage: ignore + if exponent % 2 == 0: + self.state.omega *= _phase(g) + elif exponent % 2 == 0.5: + self._z(pauli_gates.Z, axis) + self._h(common_gates.H, axis) + self.state.omega *= _phase(g) * (1 + 1j) / (2 ** 0.5) + elif exponent % 2 == 1: + self._z(pauli_gates.Z, axis) + self._h(common_gates.H, axis) + self._z(pauli_gates.Z, axis) + self._h(common_gates.H, axis) + self.state.omega *= _phase(g) * 1j + elif exponent % 2 == 1.5: + self._h(common_gates.H, axis) + self._z(pauli_gates.Z, axis) + self.state.omega *= _phase(g) * (1 - 1j) / (2 ** 0.5) + + def _z(self, g: common_gates.ZPowGate, axis: int): + exponent = g.exponent + state = self.state + if exponent % 2 != 0: + if exponent % 0.5 != 0.0: + raise ValueError('Z exponent must be half integer') # coverage: ignore + effective_exponent = exponent % 2 + for _ in range(int(effective_exponent * 2)): + # Prescription for S left multiplication. + # Reference: https://arxiv.org/abs/1808.00128 Proposition 4 end + state.M[axis, :] ^= state.G[axis, :] + state.gamma[axis] = (state.gamma[axis] - 1) % 4 + state.omega *= _phase(g) + + def _h(self, g: common_gates.HPowGate, axis: int): + exponent = g.exponent + state = self.state + if exponent % 2 != 0: + if exponent % 1 != 0: + raise ValueError('H exponent must be integer') # coverage: ignore + # Prescription for H left multiplication + # Reference: https://arxiv.org/abs/1808.00128 + # Equations 48, 49 and Proposition 4 + t = state.s ^ (state.G[axis, :] & state.v) + u = state.s ^ (state.F[axis, :] & (~state.v)) ^ (state.M[axis, :] & state.v) + alpha = sum(state.G[axis, :] & (~state.v) & state.s) % 2 + beta = sum(state.M[axis, :] & (~state.v) & state.s) + beta += sum(state.F[axis, :] & state.v & state.M[axis, :]) + beta += sum(state.F[axis, :] & state.v & state.s) + beta %= 2 + delta = (state.gamma[axis] + 2 * (alpha + beta)) % 4 + state.update_sum(t, u, delta=delta, alpha=alpha) + state.omega *= _phase(g) + + def _cz(self, g: common_gates.CZPowGate, control_axis: int, target_axis: int): + exponent = g.exponent + state = self.state + if exponent % 2 != 0: + if exponent % 1 != 0: + raise ValueError('CZ exponent must be integer') # coverage: ignore + # Prescription for CZ left multiplication. + # Reference: https://arxiv.org/abs/1808.00128 Proposition 4 end + state.M[control_axis, :] ^= state.G[target_axis, :] + state.M[target_axis, :] ^= state.G[control_axis, :] + state.omega *= _phase(g) + + def _cx(self, g: common_gates.CXPowGate, control_axis: int, target_axis: int): + exponent = g.exponent + state = self.state + if exponent % 2 != 0: + if exponent % 1 != 0: + raise ValueError('CX exponent must be integer') # coverage: ignore + # Prescription for CX left multiplication. + # Reference: https://arxiv.org/abs/1808.00128 Proposition 4 end + state.gamma[control_axis] = ( + state.gamma[control_axis] + + state.gamma[target_axis] + + 2 * (sum(state.M[control_axis, :] & state.F[target_axis, :]) % 2) + ) % 4 + state.G[target_axis, :] ^= state.G[control_axis, :] + state.F[control_axis, :] ^= state.F[target_axis, :] + state.M[control_axis, :] ^= state.M[target_axis, :] + state.omega *= _phase(g) + + def _global_phase(self, g: global_phase_op.GlobalPhaseGate): + self.state.omega *= g.coefficient + + +def _phase(gate): + return np.exp(1j * np.pi * gate.global_shift * gate.exponent)