Merge pull request tensorflow#381 from quantumlib/qsim-noise-model

Support noise model in simulator construction

qsimcirq/qsim_simulator.py

@@ -18,7 +18,7 @@
 
 from cirq import (
     circuits,
-    linalg,
+    devices,
     ops,
     protocols,
     sim,
@@ -191,6 +191,7 @@ def __init__(
         self,
         qsim_options: Union[None, Dict, QSimOptions] = None,
         seed: value.RANDOM_STATE_OR_SEED_LIKE = None,
+        noise: "devices.NOISE_MODEL_LIKE" = None,
         circuit_memoization_size: int = 0,
     ):
         """Creates a new QSimSimulator using the given options and seed.
@@ -200,6 +201,8 @@ def __init__(
                 to use for all circuits run using this simulator. See the
                 QSimOptions class for details.
             seed: A random state or seed object, as defined in cirq.value.
+            noise: A cirq.NoiseModel to apply to all circuits simulated with
+                this simulator.
             circuit_memoization_size: The number of last translated circuits
                 to be memoized from simulation executions, to eliminate
                 translation overhead. Every simulation will perform a linear
@@ -225,6 +228,7 @@ def __init__(
         self._prng = value.parse_random_state(seed)
         self.qsim_options = QSimOptions().as_dict()
         self.qsim_options.update(qsim_options)
+        self.noise = devices.NoiseModel.from_noise_model_like(noise)
 
         # module to use for simulation
         if self.qsim_options["g"]:
@@ -303,11 +307,18 @@ def _sample_measure_results(
             ValueError: If there are multiple MeasurementGates with the same key,
                 or if repetitions is negative.
         """
-        if not isinstance(program, qsimc.QSimCircuit):
-            program = qsimc.QSimCircuit(program, device=program.device)
+
+        # Add noise to the circuit if a noise model was provided.
+        all_qubits = program.all_qubits()
+        program = qsimc.QSimCircuit(
+            self.noise.noisy_moments(program, sorted(all_qubits))
+            if self.noise is not devices.NO_NOISE
+            else program,
+            device=program.device,
+        )
 
         # Compute indices of measured qubits
-        ordered_qubits = ops.QubitOrder.DEFAULT.order_for(program.all_qubits())
+        ordered_qubits = ops.QubitOrder.DEFAULT.order_for(all_qubits)
         num_qubits = len(ordered_qubits)
 
         qubit_map = {qubit: index for index, qubit in enumerate(ordered_qubits)}
@@ -424,13 +435,18 @@ def compute_amplitudes_sweep(
         Returns:
             List of amplitudes.
         """
-        if not isinstance(program, qsimc.QSimCircuit):
-            program = qsimc.QSimCircuit(program, device=program.device)
 
-        # qsim numbers qubits in reverse order from cirq
-        cirq_order = ops.QubitOrder.as_qubit_order(qubit_order).order_for(
-            program.all_qubits()
+        # Add noise to the circuit if a noise model was provided.
+        all_qubits = program.all_qubits()
+        program = qsimc.QSimCircuit(
+            self.noise.noisy_moments(program, sorted(all_qubits))
+            if self.noise is not devices.NO_NOISE
+            else program,
+            device=program.device,
         )
+
+        # qsim numbers qubits in reverse order from cirq
+        cirq_order = ops.QubitOrder.as_qubit_order(qubit_order).order_for(all_qubits)
         num_qubits = len(cirq_order)
         bitstrings = [
             format(bitstring, "b").zfill(num_qubits)[::-1] for bitstring in bitstrings
@@ -502,17 +518,22 @@ def simulate_sweep(
             initial_state = 0
         if not isinstance(initial_state, (int, np.ndarray)):
             raise TypeError("initial_state must be an int or state vector.")
-        if not isinstance(program, qsimc.QSimCircuit):
-            program = qsimc.QSimCircuit(program, device=program.device)
+
+        # Add noise to the circuit if a noise model was provided.
+        all_qubits = program.all_qubits()
+        program = qsimc.QSimCircuit(
+            self.noise.noisy_moments(program, sorted(all_qubits))
+            if self.noise is not devices.NO_NOISE
+            else program,
+            device=program.device,
+        )
 
         options = {}
         options.update(self.qsim_options)
 
         param_resolvers = study.to_resolvers(params)
         # qsim numbers qubits in reverse order from cirq
-        cirq_order = ops.QubitOrder.as_qubit_order(qubit_order).order_for(
-            program.all_qubits()
-        )
+        cirq_order = ops.QubitOrder.as_qubit_order(qubit_order).order_for(all_qubits)
         qsim_order = list(reversed(cirq_order))
         num_qubits = len(qsim_order)
         if isinstance(initial_state, np.ndarray):
@@ -610,9 +631,8 @@ def simulate_expectation_values_sweep(
             observables = [observables]
         psumlist = [ops.PauliSum.wrap(pslike) for pslike in observables]
 
-        cirq_order = ops.QubitOrder.as_qubit_order(qubit_order).order_for(
-            program.all_qubits()
-        )
+        all_qubits = program.all_qubits()
+        cirq_order = ops.QubitOrder.as_qubit_order(qubit_order).order_for(all_qubits)
         qsim_order = list(reversed(cirq_order))
         num_qubits = len(qsim_order)
         qubit_map = {qubit: index for index, qubit in enumerate(qsim_order)}
@@ -635,8 +655,14 @@ def simulate_expectation_values_sweep(
             initial_state = 0
         if not isinstance(initial_state, (int, np.ndarray)):
             raise TypeError("initial_state must be an int or state vector.")
-        if not isinstance(program, qsimc.QSimCircuit):
-            program = qsimc.QSimCircuit(program, device=program.device)
+
+        # Add noise to the circuit if a noise model was provided.
+        program = qsimc.QSimCircuit(
+            self.noise.noisy_moments(program, sorted(all_qubits))
+            if self.noise is not devices.NO_NOISE
+            else program,
+            device=program.device,
+        )
 
         options = {}
         options.update(self.qsim_options)

qsimcirq_tests/qsimcirq_test.py

@@ -988,6 +988,37 @@ def test_noise_aggregation():
     assert cirq.approx_eq(qsim_evs, expected_evs, atol=0.05)
 
 
+def test_noise_model():
+    q0, q1 = cirq.LineQubit.range(2)
+
+    circuit = cirq.Circuit(cirq.X(q0), cirq.CNOT(q0, q1), cirq.measure(q0, q1, key="m"))
+    quiet_sim = qsimcirq.QSimSimulator()
+    quiet_results = quiet_sim.run(circuit, repetitions=100)
+    assert quiet_results.histogram(key="m")[0b11] == 100
+
+    class ReadoutError(cirq.NoiseModel):
+        def noisy_operation(self, operation: "cirq.Operation") -> "cirq.OP_TREE":
+            if isinstance(operation.gate, cirq.MeasurementGate):
+                return [cirq.X.on_each(*operation.qubits), operation]
+            return [operation]
+
+    noisy_sim = qsimcirq.QSimSimulator(noise=ReadoutError())
+    noisy_results = noisy_sim.run(circuit, repetitions=100)
+    # ReadoutError will flip both qubits.
+    assert noisy_results.histogram(key="m")[0b00] == 100
+
+    noisy_state = noisy_sim.simulate(circuit)
+    assert cirq.approx_eq(noisy_state.state_vector(), [1, 0, 0, 0])
+
+    obs = cirq.Z(q0) + cirq.Z(q1)
+    noisy_evs = noisy_sim.simulate_expectation_values(
+        circuit,
+        observables=obs,
+        permit_terminal_measurements=True,
+    )
+    assert noisy_evs == [2]
+
+
 def test_multi_qubit_fusion():
     q0, q1, q2, q3 = cirq.LineQubit.range(4)
     qubits = [q0, q1, q2, q3]