diff --git a/cirq/__init__.py b/cirq/__init__.py index 206a94611ef..9f27af266c7 100644 --- a/cirq/__init__.py +++ b/cirq/__init__.py @@ -414,6 +414,7 @@ QuantumState, quantum_state, STATE_VECTOR_LIKE, + StabilizerState, superoperator_to_choi, superoperator_to_kraus, to_valid_density_matrix, diff --git a/cirq/qis/__init__.py b/cirq/qis/__init__.py index d2ca57d32f9..9134f0b8c33 100644 --- a/cirq/qis/__init__.py +++ b/cirq/qis/__init__.py @@ -25,7 +25,7 @@ superoperator_to_kraus, ) -from cirq.qis.clifford_tableau import CliffordTableau +from cirq.qis.clifford_tableau import CliffordTableau, StabilizerState from cirq.qis.measures import ( entanglement_fidelity, diff --git a/cirq/qis/clifford_tableau.py b/cirq/qis/clifford_tableau.py index d231c58268b..017082961f0 100644 --- a/cirq/qis/clifford_tableau.py +++ b/cirq/qis/clifford_tableau.py @@ -12,17 +12,118 @@ # See the License for the specific language governing permissions and # limitations under the License. +import abc from typing import Any, Dict, List, TYPE_CHECKING import numpy as np from cirq import protocols -from cirq.value import big_endian_int_to_digits +from cirq.value import big_endian_int_to_digits, linear_dict if TYPE_CHECKING: import cirq -class CliffordTableau: +class StabilizerState(metaclass=abc.ABCMeta): + """Interface for quantum stabilizer state representations. + + This interface is used for CliffordTableau and StabilizerChForm quantum + state representations, allowing simulators to act on them abstractly. + """ + + @abc.abstractmethod + def apply_x(self, axis: int, exponent: float = 1, global_shift: float = 0): + """Apply an X operation to the state. + + Args: + axis: The axis to which the operation should be applied. + exponent: The exponent of the X operation, must be a half-integer. + global_shift: The global phase shift of the raw operation, prior to + exponentiation. Typically the value in `gate.global_shift`. + Raises: + ValueError: If the exponent is not half-integer. + """ + + @abc.abstractmethod + def apply_y(self, axis: int, exponent: float = 1, global_shift: float = 0): + """Apply an Y operation to the state. + + Args: + axis: The axis to which the operation should be applied. + exponent: The exponent of the Y operation, must be a half-integer. + global_shift: The global phase shift of the raw operation, prior to + exponentiation. Typically the value in `gate.global_shift`. + Raises: + ValueError: If the exponent is not half-integer. + """ + + @abc.abstractmethod + def apply_z(self, axis: int, exponent: float = 1, global_shift: float = 0): + """Apply a Z operation to the state. + + Args: + axis: The axis to which the operation should be applied. + exponent: The exponent of the Z operation, must be a half-integer. + global_shift: The global phase shift of the raw operation, prior to + exponentiation. Typically the value in `gate.global_shift`. + Raises: + ValueError: If the exponent is not half-integer. + """ + + @abc.abstractmethod + def apply_h(self, axis: int, exponent: float = 1, global_shift: float = 0): + """Apply an H operation to the state. + + Args: + axis: The axis to which the operation should be applied. + exponent: The exponent of the H operation, must be an integer. + global_shift: The global phase shift of the raw operation, prior to + exponentiation. Typically the value in `gate.global_shift`. + Raises: + ValueError: If the exponent is not an integer. + """ + + @abc.abstractmethod + def apply_cz( + self, control_axis: int, target_axis: int, exponent: float = 1, global_shift: float = 0 + ): + """Apply a CZ operation to the state. + + Args: + control_axis: The control axis of the operation. + target_axis: The axis to which the operation should be applied. + exponent: The exponent of the CZ operation, must be an integer. + global_shift: The global phase shift of the raw operation, prior to + exponentiation. Typically the value in `gate.global_shift`. + Raises: + ValueError: If the exponent is not an integer. + """ + + @abc.abstractmethod + def apply_cx( + self, control_axis: int, target_axis: int, exponent: float = 1, global_shift: float = 0 + ): + """Apply a CX operation to the state. + + Args: + control_axis: The control axis of the operation. + target_axis: The axis to which the operation should be applied. + exponent: The exponent of the CX operation, must be an integer. + global_shift: The global phase shift of the raw operation, prior to + exponentiation. Typically the value in `gate.global_shift`. + Raises: + ValueError: If the exponent is not an integer. + """ + + @abc.abstractmethod + def apply_global_phase(self, coefficient: linear_dict.Scalar): + """Apply a global phase to the state. + + Args: + coefficient: The global phase to apply. + """ + + +class CliffordTableau(StabilizerState): """Tableau representation of a stabilizer state (based on Aaronson and Gottesman 2006). @@ -375,3 +476,105 @@ def _measure(self, q, prng: np.random.RandomState) -> int: self.rs[p] = bool(prng.randint(2)) return int(self.rs[p]) + + def apply_x(self, axis: int, exponent: float = 1, global_shift: float = 0): + if exponent % 2 == 0: + return + if exponent % 0.5 != 0.0: + raise ValueError('X exponent must be half integer') # coverage: ignore + effective_exponent = exponent % 2 + if effective_exponent == 0.5: + self.xs[:, axis] ^= self.zs[:, axis] + self.rs[:] ^= self.xs[:, axis] & self.zs[:, axis] + elif effective_exponent == 1: + self.rs[:] ^= self.zs[:, axis] + elif effective_exponent == 1.5: + self.rs[:] ^= self.xs[:, axis] & self.zs[:, axis] + self.xs[:, axis] ^= self.zs[:, axis] + + def apply_y(self, axis: int, exponent: float = 1, global_shift: float = 0): + if exponent % 2 == 0: + return + if exponent % 0.5 != 0.0: + raise ValueError('Y exponent must be half integer') # coverage: ignore + effective_exponent = exponent % 2 + if effective_exponent == 0.5: + self.rs[:] ^= self.xs[:, axis] & (~self.zs[:, axis]) + (self.xs[:, axis], self.zs[:, axis]) = ( + self.zs[:, axis].copy(), + self.xs[:, axis].copy(), + ) + elif effective_exponent == 1: + self.rs[:] ^= self.xs[:, axis] ^ self.zs[:, axis] + elif effective_exponent == 1.5: + self.rs[:] ^= ~(self.xs[:, axis]) & self.zs[:, axis] + (self.xs[:, axis], self.zs[:, axis]) = ( + self.zs[:, axis].copy(), + self.xs[:, axis].copy(), + ) + + def apply_z(self, axis: int, exponent: float = 1, global_shift: float = 0): + if exponent % 2 == 0: + return + if exponent % 0.5 != 0.0: + raise ValueError('Z exponent must be half integer') # coverage: ignore + effective_exponent = exponent % 2 + if effective_exponent == 0.5: + self.rs[:] ^= self.xs[:, axis] & self.zs[:, axis] + self.zs[:, axis] ^= self.xs[:, axis] + elif effective_exponent == 1: + self.rs[:] ^= self.xs[:, axis] + elif effective_exponent == 1.5: + self.rs[:] ^= self.xs[:, axis] & (~self.zs[:, axis]) + self.zs[:, axis] ^= self.xs[:, axis] + + def apply_h(self, axis: int, exponent: float = 1, global_shift: float = 0): + if exponent % 2 == 0: + return + if exponent % 1 != 0: + raise ValueError('H exponent must be integer') # coverage: ignore + self.apply_y(axis, 0.5) + self.apply_x(axis) + + def apply_cz( + self, control_axis: int, target_axis: int, exponent: float = 1, global_shift: float = 0 + ): + if exponent % 2 == 0: + return + if exponent % 1 != 0: + raise ValueError('CZ exponent must be integer') # coverage: ignore + (self.xs[:, target_axis], self.zs[:, target_axis]) = ( + self.zs[:, target_axis].copy(), + self.xs[:, target_axis].copy(), + ) + self.rs[:] ^= self.xs[:, target_axis] & self.zs[:, target_axis] + self.rs[:] ^= ( + self.xs[:, control_axis] + & self.zs[:, target_axis] + & (~(self.xs[:, target_axis] ^ self.zs[:, control_axis])) + ) + self.xs[:, target_axis] ^= self.xs[:, control_axis] + self.zs[:, control_axis] ^= self.zs[:, target_axis] + (self.xs[:, target_axis], self.zs[:, target_axis]) = ( + self.zs[:, target_axis].copy(), + self.xs[:, target_axis].copy(), + ) + self.rs[:] ^= self.xs[:, target_axis] & self.zs[:, target_axis] + + def apply_cx( + self, control_axis: int, target_axis: int, exponent: float = 1, global_shift: float = 0 + ): + if exponent % 2 == 0: + return + if exponent % 1 != 0: + raise ValueError('CX exponent must be integer') # coverage: ignore + self.rs[:] ^= ( + self.xs[:, control_axis] + & self.zs[:, target_axis] + & (~(self.xs[:, target_axis] ^ self.zs[:, control_axis])) + ) + self.xs[:, target_axis] ^= self.xs[:, control_axis] + self.zs[:, control_axis] ^= self.zs[:, target_axis] + + def apply_global_phase(self, coefficient: linear_dict.Scalar): + pass diff --git a/cirq/sim/clifford/act_on_clifford_tableau_args.py b/cirq/sim/clifford/act_on_clifford_tableau_args.py index f2be23883ad..e8608c6b93c 100644 --- a/cirq/sim/clifford/act_on_clifford_tableau_args.py +++ b/cirq/sim/clifford/act_on_clifford_tableau_args.py @@ -18,15 +18,14 @@ import numpy as np -from cirq.ops import common_gates -from cirq.ops import global_phase_op +from cirq.qis import clifford_tableau from cirq.sim.clifford.act_on_stabilizer_args import ActOnStabilizerArgs if TYPE_CHECKING: import cirq -class ActOnCliffordTableauArgs(ActOnStabilizerArgs): +class ActOnCliffordTableauArgs(ActOnStabilizerArgs[clifford_tableau.CliffordTableau]): """State and context for an operation acting on a clifford tableau.""" def __init__( @@ -53,19 +52,23 @@ def __init__( simulation. """ super().__init__( + state=tableau, prng=prng, qubits=qubits, log_of_measurement_results=log_of_measurement_results, classical_data=classical_data, ) - self.tableau = tableau + + @property + def tableau(self) -> 'cirq.CliffordTableau': + return self.state 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] + return [self.state._measure(self.qubit_map[q], self.prng) for q in qubits] def _on_copy(self, target: 'ActOnCliffordTableauArgs', deep_copy_buffers: bool = True): - target.tableau = self.tableau.copy() + target._state = self.state.copy() def sample( self, @@ -75,112 +78,3 @@ def sample( ) -> np.ndarray: # Unnecessary for now but can be added later if there is a use case. raise NotImplementedError() - - 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/sim/clifford/act_on_stabilizer_args.py b/cirq/sim/clifford/act_on_stabilizer_args.py index 24be68544e4..d9051f21147 100644 --- a/cirq/sim/clifford/act_on_stabilizer_args.py +++ b/cirq/sim/clifford/act_on_stabilizer_args.py @@ -13,11 +13,11 @@ # limitations under the License. import abc -from typing import Any, Sequence, TYPE_CHECKING, Union +from typing import Any, Dict, Generic, Optional, Sequence, TYPE_CHECKING, TypeVar, Union import numpy as np -from cirq import ops, protocols, linalg +from cirq import linalg, ops, protocols 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 @@ -28,9 +28,47 @@ import cirq -class ActOnStabilizerArgs(ActOnArgs, metaclass=abc.ABCMeta): +TStabilizerState = TypeVar('TStabilizerState', bound='cirq.StabilizerState') + + +class ActOnStabilizerArgs(ActOnArgs, Generic[TStabilizerState], metaclass=abc.ABCMeta): """Abstract wrapper around a stabilizer state for the act_on protocol.""" + def __init__( + self, + state: TStabilizerState, + prng: Optional[np.random.RandomState] = None, + log_of_measurement_results: Optional[Dict[str, Any]] = None, + qubits: Optional[Sequence['cirq.Qid']] = None, + classical_data: Optional['cirq.ClassicalDataStore'] = None, + ): + """Initializes the ActOnStabilizerArgs. + + Args: + state: The quantum stabilizer state to use in the simulation or + act_on invocation. + prng: The pseudo random number generator to use for probabilistic + effects. + qubits: Determines the canonical ordering of the qubits. This + is often used in specifying the initial state, i.e. the + ordering of the computational basis states. + log_of_measurement_results: A mutable object that measurements are + being recorded into. + classical_data: The shared classical data container for this + simulation. + """ + super().__init__( + prng=prng, + qubits=qubits, + log_of_measurement_results=log_of_measurement_results, + classical_data=classical_data, + ) + self._state = state + + @property + def state(self) -> TStabilizerState: + return self._state + def _act_on_fallback_( self, action: Union['cirq.Operation', 'cirq.Gate'], @@ -52,45 +90,15 @@ def _act_on_fallback_( 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): + def _swap( + self, control_axis: int, target_axis: int, exponent: float = 1, global_shift: float = 0 + ): """Apply a SWAP gate""" - if g.exponent % 1 != 0: + if 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) + self._state.apply_cx(control_axis, target_axis) + self._state.apply_cx(target_axis, control_axis, exponent, global_shift) + self._state.apply_cx(control_axis, target_axis) def _strat_apply_gate(self, val: Any, qubits: Sequence['cirq.Qid']) -> bool: if not protocols.has_stabilizer_effect(val): @@ -98,21 +106,21 @@ def _strat_apply_gate(self, val: Any, qubits: Sequence['cirq.Qid']) -> bool: 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]) + self._state.apply_x(axes[0], gate.exponent, gate.global_shift) elif isinstance(gate, common_gates.YPowGate): - self._y(gate, axes[0]) + self._state.apply_y(axes[0], gate.exponent, gate.global_shift) elif isinstance(gate, common_gates.ZPowGate): - self._z(gate, axes[0]) + self._state.apply_z(axes[0], gate.exponent, gate.global_shift) elif isinstance(gate, common_gates.HPowGate): - self._h(gate, axes[0]) + self._state.apply_h(axes[0], gate.exponent, gate.global_shift) elif isinstance(gate, common_gates.CXPowGate): - self._cx(gate, axes[0], axes[1]) + self._state.apply_cx(axes[0], axes[1], gate.exponent, gate.global_shift) elif isinstance(gate, common_gates.CZPowGate): - self._cz(gate, axes[0], axes[1]) + self._state.apply_cz(axes[0], axes[1], gate.exponent, gate.global_shift) elif isinstance(gate, global_phase_op.GlobalPhaseGate): - self._global_phase(gate) + self._state.apply_global_phase(gate.coefficient) elif isinstance(gate, swap_gates.SwapPowGate): - self._swap(gate, axes[0], axes[1]) + self._swap(axes[0], axes[1], gate.exponent, gate.global_shift) else: return NotImplemented return True @@ -150,7 +158,7 @@ def _strat_act_from_single_qubit_decompose( 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])) + self._state.apply_global_phase(u[k] / final_unitary[k]) return True return NotImplemented diff --git a/cirq/sim/clifford/act_on_stabilizer_ch_form_args.py b/cirq/sim/clifford/act_on_stabilizer_ch_form_args.py index c84f0020406..6775fafa2ce 100644 --- a/cirq/sim/clifford/act_on_stabilizer_ch_form_args.py +++ b/cirq/sim/clifford/act_on_stabilizer_ch_form_args.py @@ -17,16 +17,16 @@ import numpy as np from cirq import _compat, value, ops, protocols -from cirq.ops import common_gates, pauli_gates -from cirq.ops import global_phase_op -from cirq.sim.clifford import clifford_simulator +from cirq.sim.clifford import stabilizer_state_ch_form from cirq.sim.clifford.act_on_stabilizer_args import ActOnStabilizerArgs if TYPE_CHECKING: import cirq -class ActOnStabilizerCHFormArgs(ActOnStabilizerArgs): +class ActOnStabilizerCHFormArgs( + ActOnStabilizerArgs[stabilizer_state_ch_form.StabilizerStateChForm] +): """Wrapper around a stabilizer state in CH form for the act_on protocol.""" @_compat.deprecated_parameter( @@ -58,41 +58,48 @@ def __init__( being recorded into. initial_state: The initial state for the simulation. This can be a full CH form passed by reference which will be modified inplace, - or a big-endian int in the computational basis. + or a big-endian int in the computational basis. If the state is + an integer, qubits must be provided in order to determine + array sizes. classical_data: The shared classical data container for this simulation. + + Raises: + ValueError: If initial state is an integer but qubits are not + provided. """ + initial_state = state or initial_state + if isinstance(initial_state, int): + if qubits is None: + raise ValueError('Must specify qubits if initial state is integer') + initial_state = stabilizer_state_ch_form.StabilizerStateChForm( + len(qubits), initial_state + ) super().__init__( + state=initial_state, prng=prng, qubits=qubits, log_of_measurement_results=log_of_measurement_results, classical_data=classical_data, ) - initial_state = state or initial_state - if isinstance(initial_state, int): - qubit_map = {q: i for i, q in enumerate(self.qubits)} - initial_state = clifford_simulator.CliffordState( - qubit_map, initial_state=initial_state - ).ch_form - self.state = initial_state 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] def _on_copy(self, target: 'ActOnStabilizerCHFormArgs', deep_copy_buffers: bool = True): - target.state = self.state.copy() + target._state = self.state.copy() def _on_kronecker_product( self, other: 'cirq.ActOnStabilizerCHFormArgs', target: 'cirq.ActOnStabilizerCHFormArgs' ): - target.state = self.state.kron(other.state) + target._state = self.state.kron(other.state) def _on_transpose_to_qubit_order( self, qubits: Sequence['cirq.Qid'], target: 'cirq.ActOnStabilizerCHFormArgs' ): axes = [self.qubit_map[q] for q in qubits] - target.state = self.state.reindex(axes) + target._state = self.state.reindex(axes) def sample( self, @@ -113,105 +120,3 @@ def sample( ) protocols.act_on(op, ch_form_args) return np.array(list(measurements.measurements.values()), dtype=bool) - - 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) diff --git a/cirq/sim/clifford/act_on_stabilizer_ch_form_args_test.py b/cirq/sim/clifford/act_on_stabilizer_ch_form_args_test.py index 0cfbbc16344..11c93536a54 100644 --- a/cirq/sim/clifford/act_on_stabilizer_ch_form_args_test.py +++ b/cirq/sim/clifford/act_on_stabilizer_ch_form_args_test.py @@ -24,6 +24,8 @@ def test_init_state(): initial_state=1, ) np.testing.assert_allclose(args.state.state_vector(), [0, 1]) + with pytest.raises(ValueError, match='Must specify qubits'): + _ = cirq.ActOnStabilizerCHFormArgs(initial_state=1) def test_deprecated_warning(): diff --git a/cirq/sim/clifford/stabilizer_state_ch_form.py b/cirq/sim/clifford/stabilizer_state_ch_form.py index baf28623f48..8a53ef5b3e2 100644 --- a/cirq/sim/clifford/stabilizer_state_ch_form.py +++ b/cirq/sim/clifford/stabilizer_state_ch_form.py @@ -16,14 +16,12 @@ import numpy as np import cirq -from cirq import protocols, value -from cirq.ops import pauli_gates -from cirq.sim import clifford +from cirq import protocols, qis, value from cirq.value import big_endian_int_to_digits @value.value_equality -class StabilizerStateChForm: +class StabilizerStateChForm(qis.StabilizerState): r"""A representation of stabilizer states using the CH form, $|\psi> = \omega U_C U_H |s>$ @@ -56,22 +54,11 @@ def __init__(self, num_qubits: int, initial_state: int = 0) -> None: self.omega: complex = 1 # Apply X for every non-zero element of initial_state - qubits = cirq.LineQubit.range(num_qubits) - args = clifford.ActOnStabilizerCHFormArgs( - prng=np.random.RandomState(), - log_of_measurement_results={}, - qubits=qubits, - initial_state=self, - ) for (i, val) in enumerate( big_endian_int_to_digits(initial_state, digit_count=num_qubits, base=2) ): if val: - protocols.act_on( - pauli_gates.X, - args, - [qubits[i]], - ) + self.apply_x(i) def _json_dict_(self) -> Dict[str, Any]: return protocols.obj_to_dict_helper(self, ['n', 'G', 'F', 'M', 'gamma', 'v', 's', 'omega']) @@ -304,3 +291,100 @@ def reindex(self, axes: Sequence[int]) -> 'cirq.StabilizerStateChForm': copy.s = self.s[axes] copy.omega = self.omega return copy + + def apply_x(self, axis: int, exponent: float = 1, global_shift: float = 0): + if exponent % 2 != 0: + if exponent % 0.5 != 0.0: + raise ValueError('X exponent must be half integer') # coverage: ignore + self.apply_h(axis) + self.apply_z(axis, exponent) + self.apply_h(axis) + self.omega *= _phase(exponent, global_shift) + + def apply_y(self, axis: int, exponent: float = 1, global_shift: float = 0): + if exponent % 0.5 != 0.0: + raise ValueError('Y exponent must be half integer') # coverage: ignore + shift = _phase(exponent, global_shift) + if exponent % 2 == 0: + self.omega *= shift + elif exponent % 2 == 0.5: + self.apply_z(axis) + self.apply_h(axis) + self.omega *= shift * (1 + 1j) / (2 ** 0.5) + elif exponent % 2 == 1: + self.apply_z(axis) + self.apply_h(axis) + self.apply_z(axis) + self.apply_h(axis) + self.omega *= shift * 1j + elif exponent % 2 == 1.5: + self.apply_h(axis) + self.apply_z(axis) + self.omega *= shift * (1 - 1j) / (2 ** 0.5) + + def apply_z(self, axis: int, exponent: float = 1, global_shift: float = 0): + 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 + self.M[axis, :] ^= self.G[axis, :] + self.gamma[axis] = (self.gamma[axis] - 1) % 4 + self.omega *= _phase(exponent, global_shift) + + def apply_h(self, axis: int, exponent: float = 1, global_shift: float = 0): + 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 = self.s ^ (self.G[axis, :] & self.v) + u = self.s ^ (self.F[axis, :] & (~self.v)) ^ (self.M[axis, :] & self.v) + alpha = sum(self.G[axis, :] & (~self.v) & self.s) % 2 + beta = sum(self.M[axis, :] & (~self.v) & self.s) + beta += sum(self.F[axis, :] & self.v & self.M[axis, :]) + beta += sum(self.F[axis, :] & self.v & self.s) + beta %= 2 + delta = (self.gamma[axis] + 2 * (alpha + beta)) % 4 + self.update_sum(t, u, delta=delta, alpha=alpha) + self.omega *= _phase(exponent, global_shift) + + def apply_cz( + self, control_axis: int, target_axis: int, exponent: float = 1, global_shift: float = 0 + ): + 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 + self.M[control_axis, :] ^= self.G[target_axis, :] + self.M[target_axis, :] ^= self.G[control_axis, :] + self.omega *= _phase(exponent, global_shift) + + def apply_cx( + self, control_axis: int, target_axis: int, exponent: float = 1, global_shift: float = 0 + ): + 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 + self.gamma[control_axis] = ( + self.gamma[control_axis] + + self.gamma[target_axis] + + 2 * (sum(self.M[control_axis, :] & self.F[target_axis, :]) % 2) + ) % 4 + self.G[target_axis, :] ^= self.G[control_axis, :] + self.F[control_axis, :] ^= self.F[target_axis, :] + self.M[control_axis, :] ^= self.M[target_axis, :] + self.omega *= _phase(exponent, global_shift) + + def apply_global_phase(self, coefficient: value.Scalar): + self.omega *= coefficient + + +def _phase(exponent, global_shift): + return np.exp(1j * np.pi * global_shift * exponent)