Allow expectation values for density matrix, MPS, Clifford

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."""

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)

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.

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."""

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()

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)

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)
+    # <Z> = (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)
+    # <X> = (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()

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)

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

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)