Merge pull request #303 from quantumlib/qsim-expval-clean

Support aggregation of expectation values in qsimcirq

pybind_interface/pybind_main.cpp

@@ -482,10 +482,33 @@ class SimulatorHelper {
                             unsigned>>& opsums_and_qubit_counts,
       bool is_noisy, const StateType& input_state) {
     auto helper = SimulatorHelper(options, is_noisy);
-    if (!helper.is_valid || !helper.simulate(input_state)) {
+    if (!helper.is_valid) {
       return {};
     }
-    return helper.get_expectation_value(opsums_and_qubit_counts);
+    if (!is_noisy) {
+      if (!helper.simulate(input_state)) {
+        return {};
+      }
+      return helper.get_expectation_value(opsums_and_qubit_counts);
+    }
+
+    // Aggregate expectation values for noisy circuits.
+    std::vector<std::complex<double>> results(
+      opsums_and_qubit_counts.size(), 0);
+    for (unsigned rep = 0; rep < helper.noisy_reps; ++rep) {
+      if (!helper.simulate(input_state)) {
+        return {};
+      }
+      auto evs = helper.get_expectation_value(opsums_and_qubit_counts);
+      for (unsigned i = 0; i < evs.size(); ++i) {
+        results[i] += evs[i];
+      }
+    }
+    double inverse_num_reps = 1.0 / helper.noisy_reps;
+    for (unsigned i = 0; i < results.size(); ++i) {
+      results[i] *= inverse_num_reps;
+    }
+    return results;
   }
 
  private:
@@ -498,6 +521,7 @@ class SimulatorHelper {
       if (is_noisy) {
         ncircuit = getNoisyCircuit(options);
         num_qubits = parseOptions<unsigned>(options, "n\0");
+        noisy_reps = parseOptions<unsigned>(options, "r\0");
       } else {
         circuit = getCircuit(options);
         num_qubits = circuit.num_qubits;
@@ -551,12 +575,16 @@ class SimulatorHelper {
   template <typename StateType>
   bool simulate(const StateType& input_state) {
     init_state(input_state);
+    bool result = false;
     if (is_noisy) {
       std::vector<uint64_t> stat;
-      return NoisyRunner::Run(
+      result = NoisyRunner::Run(
         get_noisy_params(), num_qubits, ncircuit, seed, scratch, state, stat);
+    } else {
+      result = Runner::Run(get_params(), circuit, state);
     }
-    return Runner::Run(get_params(), circuit, state);
+    seed += 1;
+    return result;
   }
 
   py::array_t<float> release_state_to_python() {
@@ -601,6 +629,7 @@ class SimulatorHelper {
 
   unsigned num_qubits;
   unsigned num_threads;
+  unsigned noisy_reps;
   unsigned max_fused_size;
   unsigned verbosity;
   unsigned seed;

qsimcirq/qsim_simulator.py

@@ -80,9 +80,13 @@ def __init__(self, qsim_options: dict = {},
             applied to all circuits run using this simulator. Accepted keys and
             their behavior are as follows:
                 - 'f': int (> 0). Maximum size of fused gates. Default: 2.
+                - 'r': int (> 0). Noisy repetitions (see below). Default: 1.
                 - 't': int (> 0). Number of threads to run on. Default: 1.
                 - 'v': int (>= 0). Log verbosity. Default: 0.
             See qsim/docs/usage.md for more details on these options.
+            "Noisy repetitions" specifies how many repetitions to aggregate
+            over when calculating expectation values for a noisy circuit.
+            Note that this does not apply to other simulation types.
         seed: A random state or seed object, as defined in cirq.value.
     
     Raises:
@@ -94,7 +98,7 @@ def __init__(self, qsim_options: dict = {},
           'used in QSimCircuit instantiation.'
       )
     self._prng = value.parse_random_state(seed)
-    self.qsim_options = {'t': 1, 'f': 2, 'v': 0}
+    self.qsim_options = {'t': 1, 'f': 2, 'v': 0, 'r': 1}
     self.qsim_options.update(qsim_options)
 
   def get_seed(self):

qsimcirq_tests/qsimcirq_test.py

@@ -778,6 +778,33 @@ def test_multi_qubit_noise(cx_qubits, noise_type):
   assert all(result_count > 0 for result_count in result_hist)
 
 
+def test_noise_aggregation():
+  q0 = cirq.LineQubit(0)
+  # damp_prob is set high to minimize test variance.
+  # Even with this setting, estimation of states and expectation values from
+  # noisy circuits is highly variable, so this test uses wide tolerances.
+  damp_prob = 0.4
+  circuit = cirq.Circuit(
+    cirq.X(q0), cirq.amplitude_damp(gamma=damp_prob).on(q0),
+  )
+  psum1 = cirq.Z(q0)
+  psum2 = cirq.X(q0)
+
+  # Test expectation value aggregation over repetitions of a noisy circuit.
+  # Repetitions are handled in C++, so overhead costs are minimal.
+  qsim_simulator = qsimcirq.QSimSimulator(qsim_options={"r": 10000}, seed=1)
+  qsim_evs = qsim_simulator.simulate_expectation_values_sweep(
+    circuit, [psum1, psum2], params={}, permit_terminal_measurements=True)
+  assert len(qsim_evs) == 1
+  assert len(qsim_evs[0]) == 2
+
+  # <Z> = (-1) * (probability of |1>) + 1 * (probability of |0>)
+  # For damp_prob = 0.4, <Z> == -0.2
+  damped_zval = damp_prob - (1 - damp_prob)
+  expected_evs = [[damped_zval, 0]]
+  assert cirq.approx_eq(qsim_evs, expected_evs, atol=0.05)
+
+
 def test_multi_qubit_fusion():
   q0, q1, q2, q3 = cirq.LineQubit.range(4)
   qubits = [q0, q1, q2, q3]