Density matrix simulator EVs (#3979)

Fixes #3964. CC @zaqqwerty

This is a near-exact copy of the behavior and tests for `sparse_simulator`. If we see the need in the future, this code could probably be folded into a `SimulatesEVsAndFullState` interface as the default `simulate_expectation_values_sweep` behavior, but that doesn't seem particularly urgent at the moment.

Relevant changes from the `sparse_simulator` code:
- `obs.expectation_from_density_matrix(result.final_density_matrix, qmap)` instead of state-vector equivalent
- added `test_simulate_noisy_expectation_values` test
- that's it

Author: 95-martin-orion

cirq-core/cirq/sim/density_matrix_simulator.py

@@ -31,6 +31,7 @@ class DensityMatrixSimulator(
         'DensityMatrixSimulatorState',
         act_on_density_matrix_args.ActOnDensityMatrixArgs,
     ],
+    simulator.SimulatesExpectationValues,
 ):
     """A simulator for density matrices and noisy quantum circuits.
 
@@ -220,6 +221,40 @@ def _create_simulator_trial_result(
             params=params, measurements=measurements, final_simulator_state=final_simulator_state
         )
 
+    # TODO(#4209): Deduplicate with identical code in sparse_simulator.
+    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_density_matrix(result.final_density_matrix, qmap)
+                    for obs in pslist
+                ]
+            )
+        return swept_evs
+
 
 class DensityMatrixStepResult(simulator.StepResult['DensityMatrixSimulatorState']):
     """A single step in the simulation of the DensityMatrixSimulator.

cirq-core/cirq/sim/density_matrix_simulator_test.py

@@ -808,6 +808,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()