Add cirq.final_density_matrix utility method (#2487)

Fixes #2421

cirq/__init__.py

@@ -298,6 +298,7 @@
 from cirq.sim import (
     bloch_vector_from_state_vector,
     StabilizerStateChForm,
+    CIRCUIT_LIKE,
     CliffordSimulator,
     CliffordState,
     CliffordSimulatorStepResult,
@@ -311,6 +312,7 @@
     dirac_notation,
     measure_density_matrix,
     measure_state_vector,
+    final_density_matrix,
     final_wavefunction,
     sample,
     sample_density_matrix,

cirq/protocols/json_serialization_test.py

@@ -166,6 +166,7 @@ def test_fail_to_resolve():
     'SupportsUnitary',
 
     # mypy types:
+    'CIRCUIT_LIKE',
     'DURATION_LIKE',
     'NOISE_MODEL_LIKE',
     'OP_TREE',

cirq/sim/__init__.py

@@ -29,6 +29,8 @@
 )
 
 from cirq.sim.mux import (
+    CIRCUIT_LIKE,
+    final_density_matrix,
     final_wavefunction,
     sample,
     sample_sweep,

cirq/sim/density_matrix_simulator.py

@@ -16,15 +16,14 @@
 
 import collections
 
-from typing import Dict, Iterator, List, Type, Union, TYPE_CHECKING
+from typing import Dict, Iterator, List, Type, Union
 
 import numpy as np
 
 from cirq import circuits, ops, protocols, study, value, devices
 from cirq.sim import density_matrix_utils, simulator
 
-if TYPE_CHECKING:
-    import cirq
+import cirq
 
 
 class _StateAndBuffers:
@@ -120,14 +119,36 @@ def __init__(self,
                  *,
                  dtype: Type[np.number] = np.complex64,
                  noise: 'cirq.NOISE_MODEL_LIKE' = None,
-                 seed: value.RANDOM_STATE_LIKE = None):
+                 seed: value.RANDOM_STATE_LIKE = None,
+                 ignore_measurement_results: bool = False):
         """Density matrix simulator.
 
          Args:
             dtype: The `numpy.dtype` used by the simulation. One of
                 `numpy.complex64` or `numpy.complex128`
             noise: A noise model to apply while simulating.
             seed: The random seed to use for this simulator.
+            ignore_measurement_results: if True, then the simulation
+                will treat measurement as dephasing instead of collapsing
+                process.
+
+                Example:
+                >>> (q0,) = cirq.LineQubit.range(1)
+                >>> circuit = cirq.Circuit(cirq.H(q0), cirq.measure(q0))
+
+                Default case (ignore_measurement_results = False):
+                >>> simulator = cirq.DensityMatrixSimulator()
+                >>> result = simulator.run(circuit)
+
+                The measurement result will be strictly one of 0 or 1.
+
+                In the other case:
+                >>> simulator = cirq.DensityMatrixSimulator(
+                ...     ignore_measurement_results = True)
+                >>> result = simulator.run(circuit)
+
+                The measurement result will be the maximally mixed state
+                with equal probability for 0 and 1.
         """
         if dtype not in {np.complex64, np.complex128}:
             raise ValueError(
@@ -136,6 +157,7 @@ def __init__(self,
         self._dtype = dtype
         self._prng = value.parse_random_state(seed)
         self.noise = devices.NoiseModel.from_noise_model_like(noise)
+        self._ignore_measurement_results = (ignore_measurement_results)
 
     def _run(self, circuit: circuits.Circuit,
              param_resolver: study.ParamResolver,
@@ -277,20 +299,27 @@ def keep(potential_op: ops.Operation) -> bool:
                 if isinstance(op.gate, ops.MeasurementGate):
                     meas = op.gate
                     if perform_measurements:
-                        invert_mask = meas.full_invert_mask()
-                        # Measure updates inline.
-                        bits, _ = density_matrix_utils.measure_density_matrix(
-                            state.tensor,
-                            indices,
-                            qid_shape=qid_shape,
-                            out=state.tensor,
-                            seed=self._prng)
-                        corrected = [
-                            bit ^ (bit < 2 and mask)
-                            for bit, mask in zip(bits, invert_mask)
-                        ]
-                        key = protocols.measurement_key(meas)
-                        measurements[key].extend(corrected)
+                        if self._ignore_measurement_results:
+                            for i, q in enumerate(op.qubits):
+                                self._apply_op_channel(
+                                    cirq.phase_damp(1).on(q), state,
+                                    [indices[i]])
+                        else:
+                            invert_mask = meas.full_invert_mask()
+                            # Measure updates inline.
+                            bits, _ = (
+                                density_matrix_utils.measure_density_matrix(
+                                    state.tensor,
+                                    indices,
+                                    qid_shape=qid_shape,
+                                    out=state.tensor,
+                                    seed=self._prng))
+                            corrected = [
+                                bit ^ (bit < 2 and mask)
+                                for bit, mask in zip(bits, invert_mask)
+                            ]
+                            key = protocols.measurement_key(meas)
+                            measurements[key].extend(corrected)
                 else:
                     # TODO: Use apply_channel similar to apply_unitary.
                     self._apply_op_channel(op, state, indices)

cirq/sim/density_matrix_simulator_test.py

@@ -97,6 +97,22 @@ def test_run_bit_flips(dtype):
                                     {'0': [[b0]], '1': [[b1]]})
 
 
+@pytest.mark.parametrize('dtype', [np.complex64, np.complex128])
+def test_run_bit_flips_with_dephasing(dtype):
+    q0, q1 = cirq.LineQubit.range(2)
+    simulator = cirq.DensityMatrixSimulator(dtype=dtype,
+                                            ignore_measurement_results=True)
+    for b0 in [0, 1]:
+        for b1 in [0, 1]:
+            circuit = cirq.Circuit((cirq.X**b0)(q0), (cirq.X**b1)(q1),
+                                   cirq.measure(q0), cirq.measure(q1))
+            result = simulator.run(circuit)
+            np.testing.assert_equal(result.measurements, {
+                '0': [[b0]],
+                '1': [[b1]]
+            })
+
+
 @pytest.mark.parametrize('dtype', [np.complex64, np.complex128])
 def test_run_qudit_increments(dtype):
     q0, q1 = cirq.LineQid.for_qid_shape((3, 4))
@@ -848,22 +864,44 @@ def test_density_matrix_trial_result_repr():
                 "qubit_map={cirq.LineQubit(0): 0}))""")
 
 
+class XAsOp(cirq.Operation):
+
+    def __init__(self, q):
+        # coverage: ignore
+        self.q = q
+
+    @property
+    def qubits(self):
+        # coverage: ignore
+        return self.q,
+
+    def with_qubits(self, *new_qubits):
+        # coverage: ignore
+        return XAsOp(new_qubits[0])
+
+    def _channel_(self):
+        # coverage: ignore
+        return cirq.channel(cirq.X)
+
+
 def test_works_on_operation():
 
     class XAsOp(cirq.Operation):
 
         def __init__(self, q):
+            # coverage: ignore
             self.q = q
 
         @property
         def qubits(self):
+            # coverage: ignore
             return self.q,
 
         def with_qubits(self, *new_qubits):
-            # coverage: ignore
-            return XAsOp(new_qubits[0])
+            raise NotImplementedError()
 
         def _channel_(self):
+            # coverage: ignore
             return cirq.channel(cirq.X)
 
     s = cirq.DensityMatrixSimulator()
@@ -873,6 +911,30 @@ def _channel_(self):
                                atol=1e-8)
 
 
+def test_works_on_operation_dephased():
+
+    class HAsOp(cirq.Operation):
+
+        def __init__(self, q):
+            self.q = q
+
+        @property
+        def qubits(self):
+            return self.q,
+
+        def with_qubits(self, *new_qubits):
+            raise NotImplementedError()
+
+        def _channel_(self):
+            return cirq.channel(cirq.H)
+
+    s = cirq.DensityMatrixSimulator(ignore_measurement_results=True)
+    c = cirq.Circuit(HAsOp(cirq.LineQubit(0)))
+    np.testing.assert_allclose(s.simulate(c).final_density_matrix,
+                               [[0.5 + 0.j, 0.5 + 0.j], [0.5 + 0.j, 0.5 + 0.j]],
+                               atol=1e-8)
+
+
 def test_works_on_pauli_string_phasor():
     a, b = cirq.LineQubit.range(2)
     c = cirq.Circuit(np.exp(1j * np.pi * cirq.X(a) * cirq.X(b)))
@@ -882,7 +944,6 @@ def test_works_on_pauli_string_phasor():
                                np.diag([0, 0, 0, 1]),
                                atol=1e-8)
 
-
 def test_works_on_pauli_string():
     a, b = cirq.LineQubit.range(2)
     c = cirq.Circuit(cirq.X(a) * cirq.X(b))