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)