diff --git a/cirq/contrib/quimb/mps_simulator.py b/cirq/contrib/quimb/mps_simulator.py index aec48cbc02a..b34ccfdccdc 100644 --- a/cirq/contrib/quimb/mps_simulator.py +++ b/cirq/contrib/quimb/mps_simulator.py @@ -225,6 +225,9 @@ def __str__(self) -> str: final = self._final_simulator_state return f'measurements: {samples}\noutput state: {final}' + def expectation_from_state(self, obs: 'cirq.PauliSum'): + return obs.expectation_from_state_vector(self.final_state.state_vector(), self.qubit_map) + class MPSSimulatorStepResult(simulator.StepResult['MPSState']): """A `StepResult` that can perform measurements.""" diff --git a/cirq/contrib/quimb/mps_simulator_test.py b/cirq/contrib/quimb/mps_simulator_test.py index ade54e10f98..3e270f9c0a1 100644 --- a/cirq/contrib/quimb/mps_simulator_test.py +++ b/cirq/contrib/quimb/mps_simulator_test.py @@ -507,3 +507,76 @@ def test_state_act_on_args_initializer(): ) assert s.axes == (2,) assert s.log_of_measurement_results == {'test': 4} + + +def test_simulate_expectation_values(): + # Compare with test_expectation_from_state_vector_two_qubit_states + # in file: cirq/ops/linear_combinations_test.py + q0, q1 = cirq.LineQubit.range(2) + psum1 = cirq.Z(q0) + 3.2 * cirq.Z(q1) + psum2 = -1 * cirq.X(q0) + 2 * cirq.X(q1) + c1 = cirq.Circuit(cirq.I(q0), cirq.X(q1)) + simulator = ccq.mps_simulator.MPSSimulator() + result = simulator.simulate_expectation_values(c1, [psum1, psum2]) + assert cirq.approx_eq(result[0], -2.2, atol=1e-6) + assert cirq.approx_eq(result[1], 0, atol=1e-6) + + c2 = cirq.Circuit(cirq.H(q0), cirq.H(q1)) + result = simulator.simulate_expectation_values(c2, [psum1, psum2]) + assert cirq.approx_eq(result[0], 0, atol=1e-6) + assert cirq.approx_eq(result[1], 1, atol=1e-6) + + psum3 = cirq.Z(q0) + cirq.X(q1) + c3 = cirq.Circuit(cirq.I(q0), cirq.H(q1)) + result = simulator.simulate_expectation_values(c3, psum3) + assert cirq.approx_eq(result[0], 2, atol=1e-6) + + +def test_simulate_expectation_values_terminal_measure(): + q0 = cirq.LineQubit(0) + circuit = cirq.Circuit(cirq.H(q0), cirq.measure(q0)) + obs = cirq.Z(q0) + simulator = ccq.mps_simulator.MPSSimulator() + with pytest.raises(ValueError): + _ = simulator.simulate_expectation_values(circuit, obs) + + results = {-1: 0, 1: 0} + for _ in range(100): + result = simulator.simulate_expectation_values( + circuit, obs, permit_terminal_measurements=True + ) + if cirq.approx_eq(result[0], -1, atol=1e-6): + results[-1] += 1 + if cirq.approx_eq(result[0], 1, atol=1e-6): + results[1] += 1 + + # With a measurement after H, the Z-observable expects a specific state. + assert results[-1] > 0 + assert results[1] > 0 + assert results[-1] + results[1] == 100 + + circuit = cirq.Circuit(cirq.H(q0)) + results = {0: 0} + for _ in range(100): + result = simulator.simulate_expectation_values( + circuit, obs, permit_terminal_measurements=True + ) + if cirq.approx_eq(result[0], 0, atol=1e-6): + results[0] += 1 + + # Without measurement after H, the Z-observable is indeterminate. + assert results[0] == 100 + + +def test_simulate_expectation_values_qubit_order(): + q0, q1, q2 = cirq.LineQubit.range(3) + circuit = cirq.Circuit(cirq.H(q0), cirq.H(q1), cirq.X(q2)) + obs = cirq.X(q0) + cirq.X(q1) - cirq.Z(q2) + simulator = ccq.mps_simulator.MPSSimulator() + + result = simulator.simulate_expectation_values(circuit, obs) + assert cirq.approx_eq(result[0], 3, atol=1e-6) + + # Adjusting the qubit order has no effect on the observables. + result_flipped = simulator.simulate_expectation_values(circuit, obs, qubit_order=[q1, q2, q0]) + assert cirq.approx_eq(result_flipped[0], 3, atol=1e-6) diff --git a/cirq/ops/linear_combinations.py b/cirq/ops/linear_combinations.py index a8433d414e6..1ad58a8b693 100644 --- a/cirq/ops/linear_combinations.py +++ b/cirq/ops/linear_combinations.py @@ -452,7 +452,7 @@ def expectation_from_state_vector( See `PauliString.expectation_from_state_vector`. Args: - state: An array representing a valid state vector. + state_vector: An array representing a valid state vector. qubit_map: A map from all qubits used in this PauliSum to the indices of the qubits that `state_vector` is defined over. atol: Absolute numerical tolerance. diff --git a/cirq/sim/clifford/clifford_simulator.py b/cirq/sim/clifford/clifford_simulator.py index 03d31e33e79..a08c368c29d 100644 --- a/cirq/sim/clifford/clifford_simulator.py +++ b/cirq/sim/clifford/clifford_simulator.py @@ -169,6 +169,9 @@ def __str__(self) -> str: final = self._final_simulator_state return f'measurements: {samples}\noutput state: {final}' + def expectation_from_state(self, obs: 'cirq.PauliSum'): + return obs.expectation_from_state_vector(self.final_state.state_vector(), self.qubit_map) + class CliffordSimulatorStepResult(simulator.StepResult['CliffordState']): """A `StepResult` that includes `StateVectorMixin` methods.""" diff --git a/cirq/sim/clifford/clifford_simulator_test.py b/cirq/sim/clifford/clifford_simulator_test.py index ebeebdb4f84..51f7ac631f0 100644 --- a/cirq/sim/clifford/clifford_simulator_test.py +++ b/cirq/sim/clifford/clifford_simulator_test.py @@ -324,6 +324,42 @@ def test_clifford_circuit(): ) +def test_clifford_circuit_expectation_value(): + q0, q1 = cirq.LineQubit.range(2) + circuit = cirq.Circuit() + + np.random.seed(0) + + for _ in range(100): + x = np.random.randint(7) + + if x == 0: + circuit.append(cirq.X(np.random.choice((q0, q1)))) + elif x == 1: + circuit.append(cirq.Z(np.random.choice((q0, q1)))) + elif x == 2: + circuit.append(cirq.Y(np.random.choice((q0, q1)))) + elif x == 3: + circuit.append(cirq.S(np.random.choice((q0, q1)))) + elif x == 4: + circuit.append(cirq.H(np.random.choice((q0, q1)))) + elif x == 5: + circuit.append(cirq.CNOT(q0, q1)) + elif x == 6: + circuit.append(cirq.CZ(q0, q1)) + + clifford_simulator = cirq.CliffordSimulator() + state_vector_simulator = cirq.Simulator() + + psum1 = cirq.Z(q0) + 3.2 * cirq.Z(q1) + psum2 = -1 * cirq.X(q0) + 2 * cirq.X(q1) + np.testing.assert_almost_equal( + clifford_simulator.simulate_expectation_values(circuit, [psum1, psum2]), + state_vector_simulator.simulate_expectation_values(circuit, [psum1, psum2]), + decimal=5, + ) + + @pytest.mark.parametrize("qubits", [cirq.LineQubit.range(2), cirq.LineQubit.range(4)]) def test_clifford_circuit_2(qubits): circuit = cirq.Circuit() diff --git a/cirq/sim/density_matrix_simulator.py b/cirq/sim/density_matrix_simulator.py index 6c1b23bc580..36e640fc18d 100644 --- a/cirq/sim/density_matrix_simulator.py +++ b/cirq/sim/density_matrix_simulator.py @@ -515,3 +515,6 @@ def __repr__(self) -> str: f'params={self.params!r}, measurements={self.measurements!r}, ' f'final_simulator_state={self._final_simulator_state!r})' ) + + def expectation_from_state(self, obs: 'cirq.PauliSum'): + return obs.expectation_from_density_matrix(self.final_density_matrix, self.qubit_map) diff --git a/cirq/sim/density_matrix_simulator_test.py b/cirq/sim/density_matrix_simulator_test.py index 6e652f8e325..0674a416a6c 100644 --- a/cirq/sim/density_matrix_simulator_test.py +++ b/cirq/sim/density_matrix_simulator_test.py @@ -791,6 +791,106 @@ def test_simulate_moment_steps_intermediate_measurement(dtype): np.testing.assert_almost_equal(step.density_matrix(), expected) +@pytest.mark.parametrize('dtype', [np.complex64, np.complex128]) +def test_simulate_expectation_values(dtype): + # Compare with test_expectation_from_state_vector_two_qubit_states + # in file: cirq/ops/linear_combinations_test.py + q0, q1 = cirq.LineQubit.range(2) + psum1 = cirq.Z(q0) + 3.2 * cirq.Z(q1) + psum2 = -1 * cirq.X(q0) + 2 * cirq.X(q1) + c1 = cirq.Circuit(cirq.I(q0), cirq.X(q1)) + simulator = cirq.DensityMatrixSimulator(dtype=dtype) + result = simulator.simulate_expectation_values(c1, [psum1, psum2]) + assert cirq.approx_eq(result[0], -2.2, atol=1e-6) + assert cirq.approx_eq(result[1], 0, atol=1e-6) + + c2 = cirq.Circuit(cirq.H(q0), cirq.H(q1)) + result = simulator.simulate_expectation_values(c2, [psum1, psum2]) + assert cirq.approx_eq(result[0], 0, atol=1e-6) + assert cirq.approx_eq(result[1], 1, atol=1e-6) + + psum3 = cirq.Z(q0) + cirq.X(q1) + c3 = cirq.Circuit(cirq.I(q0), cirq.H(q1)) + result = simulator.simulate_expectation_values(c3, psum3) + assert cirq.approx_eq(result[0], 2, atol=1e-6) + + +@pytest.mark.parametrize('dtype', [np.complex64, np.complex128]) +def test_simulate_noisy_expectation_values(dtype): + q0 = cirq.LineQubit(0) + psums = [cirq.Z(q0), cirq.X(q0)] + c1 = cirq.Circuit( + cirq.X(q0), + cirq.amplitude_damp(gamma=0.1).on(q0), + ) + simulator = cirq.DensityMatrixSimulator(dtype=dtype) + result = simulator.simulate_expectation_values(c1, psums) + # = (gamma - 1) + gamma = -0.8 + assert cirq.approx_eq(result[0], -0.8, atol=1e-6) + assert cirq.approx_eq(result[1], 0, atol=1e-6) + + c2 = cirq.Circuit( + cirq.H(q0), + cirq.depolarize(p=0.3).on(q0), + ) + result = simulator.simulate_expectation_values(c2, psums) + assert cirq.approx_eq(result[0], 0, atol=1e-6) + # = (1 - p) + (-p / 3) = 0.6 + assert cirq.approx_eq(result[1], 0.6, atol=1e-6) + + +@pytest.mark.parametrize('dtype', [np.complex64, np.complex128]) +def test_simulate_expectation_values_terminal_measure(dtype): + q0 = cirq.LineQubit(0) + circuit = cirq.Circuit(cirq.H(q0), cirq.measure(q0)) + obs = cirq.Z(q0) + simulator = cirq.DensityMatrixSimulator(dtype=dtype) + with pytest.raises(ValueError): + _ = simulator.simulate_expectation_values(circuit, obs) + + results = {-1: 0, 1: 0} + for _ in range(100): + result = simulator.simulate_expectation_values( + circuit, obs, permit_terminal_measurements=True + ) + if cirq.approx_eq(result[0], -1, atol=1e-6): + results[-1] += 1 + if cirq.approx_eq(result[0], 1, atol=1e-6): + results[1] += 1 + + # With a measurement after H, the Z-observable expects a specific state. + assert results[-1] > 0 + assert results[1] > 0 + assert results[-1] + results[1] == 100 + + circuit = cirq.Circuit(cirq.H(q0)) + results = {0: 0} + for _ in range(100): + result = simulator.simulate_expectation_values( + circuit, obs, permit_terminal_measurements=True + ) + if cirq.approx_eq(result[0], 0, atol=1e-6): + results[0] += 1 + + # Without measurement after H, the Z-observable is indeterminate. + assert results[0] == 100 + + +@pytest.mark.parametrize('dtype', [np.complex64, np.complex128]) +def test_simulate_expectation_values_qubit_order(dtype): + q0, q1, q2 = cirq.LineQubit.range(3) + circuit = cirq.Circuit(cirq.H(q0), cirq.H(q1), cirq.X(q2)) + obs = cirq.X(q0) + cirq.X(q1) - cirq.Z(q2) + simulator = cirq.DensityMatrixSimulator(dtype=dtype) + + result = simulator.simulate_expectation_values(circuit, obs) + assert cirq.approx_eq(result[0], 3, atol=1e-6) + + # Adjusting the qubit order has no effect on the observables. + result_flipped = simulator.simulate_expectation_values(circuit, obs, qubit_order=[q1, q2, q0]) + assert cirq.approx_eq(result_flipped[0], 3, atol=1e-6) + + def test_density_matrix_simulator_state_eq(): q0, q1 = cirq.LineQubit.range(2) eq = cirq.testing.EqualsTester() diff --git a/cirq/sim/simulator.py b/cirq/sim/simulator.py index 3b5e3b5d48b..3fa35b3fad3 100644 --- a/cirq/sim/simulator.py +++ b/cirq/sim/simulator.py @@ -302,7 +302,9 @@ def simulate_expectation_values_sweep( """ -class SimulatesFinalState(Generic[TSimulationTrialResult], metaclass=abc.ABCMeta): +class SimulatesFinalState( + Generic[TSimulationTrialResult], SimulatesExpectationValues, metaclass=abc.ABCMeta +): """Simulator that allows access to the simulator's final state. Implementors of this interface should implement the simulate_sweep @@ -372,6 +374,33 @@ def simulate_sweep( """ raise NotImplementedError() + def simulate_expectation_values_sweep( + self, + program: 'cirq.Circuit', + observables: Union['cirq.PauliSumLike', List['cirq.PauliSumLike']], + params: 'study.Sweepable', + qubit_order: ops.QubitOrderOrList = ops.QubitOrder.DEFAULT, + initial_state: Any = None, + permit_terminal_measurements: bool = False, + ) -> List[List[float]]: + if not permit_terminal_measurements and program.are_any_measurements_terminal(): + raise ValueError( + 'Provided circuit has terminal measurements, which may ' + 'skew expectation values. If this is intentional, set ' + 'permit_terminal_measurements=True.' + ) + swept_evs = [] + qubit_order = ops.QubitOrder.as_qubit_order(qubit_order) + if not isinstance(observables, List): + observables = [observables] + pslist = [ops.PauliSum.wrap(pslike) for pslike in observables] + for param_resolver in study.to_resolvers(params): + result = self.simulate( + program, param_resolver, qubit_order=qubit_order, initial_state=initial_state + ) + swept_evs.append([result.expectation_from_state(obs) for obs in pslist]) + return swept_evs + class SimulatesIntermediateState( Generic[TStepResult, TSimulationTrialResult, TSimulatorState], @@ -736,6 +765,12 @@ def qubit_map(self) -> Dict[ops.Qid, int]: def _qid_shape_(self) -> Tuple[int, ...]: return _qubit_map_to_shape(self.qubit_map) + def expectation_from_state(self, obs: 'cirq.PauliSum'): + """Subclasses can implement this to allow expectation values derived + from the final state. + """ + raise NotImplementedError() + def _qubit_map_to_shape(qubit_map: Dict[ops.Qid, int]) -> Tuple[int, ...]: qid_shape: List[int] = [-1] * len(qubit_map) diff --git a/cirq/sim/sparse_simulator.py b/cirq/sim/sparse_simulator.py index 51a35438eee..d48f6b5a61f 100644 --- a/cirq/sim/sparse_simulator.py +++ b/cirq/sim/sparse_simulator.py @@ -16,14 +16,12 @@ import collections from typing import ( - Any, Dict, Iterator, List, Type, TYPE_CHECKING, DefaultDict, - Union, cast, ) @@ -46,7 +44,6 @@ class Simulator( simulator.SimulatesSamples, state_vector_simulator.SimulatesIntermediateStateVector['SparseSimulatorStep'], - simulator.SimulatesExpectationValues, ): """A sparse matrix state vector simulator that uses numpy. @@ -268,39 +265,6 @@ def _base_iterator( ) sim_state.log_of_measurement_results.clear() - def simulate_expectation_values_sweep( - self, - program: 'cirq.Circuit', - observables: Union['cirq.PauliSumLike', List['cirq.PauliSumLike']], - params: 'study.Sweepable', - qubit_order: ops.QubitOrderOrList = ops.QubitOrder.DEFAULT, - initial_state: Any = None, - permit_terminal_measurements: bool = False, - ) -> List[List[float]]: - if not permit_terminal_measurements and program.are_any_measurements_terminal(): - raise ValueError( - 'Provided circuit has terminal measurements, which may ' - 'skew expectation values. If this is intentional, set ' - 'permit_terminal_measurements=True.' - ) - swept_evs = [] - qubit_order = ops.QubitOrder.as_qubit_order(qubit_order) - qmap = {q: i for i, q in enumerate(qubit_order.order_for(program.all_qubits()))} - if not isinstance(observables, List): - observables = [observables] - pslist = [ops.PauliSum.wrap(pslike) for pslike in observables] - for param_resolver in study.to_resolvers(params): - result = self.simulate( - program, param_resolver, qubit_order=qubit_order, initial_state=initial_state - ) - swept_evs.append( - [ - obs.expectation_from_state_vector(result.final_state_vector, qmap) - for obs in pslist - ] - ) - return swept_evs - class SparseSimulatorStep( state_vector.StateVectorMixin, state_vector_simulator.StateVectorStepResult diff --git a/cirq/sim/state_vector_simulator.py b/cirq/sim/state_vector_simulator.py index 4e6e013833e..42217b75eee 100644 --- a/cirq/sim/state_vector_simulator.py +++ b/cirq/sim/state_vector_simulator.py @@ -193,3 +193,6 @@ def __repr__(self) -> str: f'measurements={self.measurements!r}, ' f'final_simulator_state={self._final_simulator_state!r})' ) + + def expectation_from_state(self, obs: 'cirq.PauliSum'): + return obs.expectation_from_state_vector(self.final_state_vector, self.qubit_map)