Added Pauli noise support to sparse simulator (#1971)

Added support for parametric Pauli noise to sparce simulator. (px, py,
pz) can be provided. Sparse simulator will apply X, Y or Z gates after
each gate and before each measurement with corresponding probability.
- This allows for simulation with parametric Pauli noise: bit flip,
phase flip, depolarizing, etc.
- For two-qubit gates noise is applied to each qubit.
- Allocation and service functions do not apply noise. Idle (Identity)
gate is also noiseless.
- Reset applies noise.
- Added tests.
- Released qubits are no longer checked to be in zero state by consumers
of the backend via calling qubit_is_zero. Instead, qubit_release returns
Boolean. True indicates, that the release was valid - the qubit was in a
zero state in a noiseless simulation. It will always be true in noisy
simulations when the release of a qubit is always valid.
- The functionality is not yet exposed outside Rust code. Will come in
subsequent PRs.

---------

Co-authored-by: Dmitry Vasilevsky <[email protected]>

Author: DmitryVasilevsky

Diff for: compiler/qsc_circuit/src/builder.rs

@@ -174,8 +174,9 @@ impl Backend for Builder {
         self.remapper.qubit_allocate()
     }
 
-    fn qubit_release(&mut self, q: usize) {
+    fn qubit_release(&mut self, q: usize) -> bool {
         self.remapper.qubit_release(q);
+        true
     }
 
     fn qubit_swap_id(&mut self, q0: usize, q1: usize) {
@@ -187,8 +188,8 @@ impl Backend for Builder {
     }
 
     fn qubit_is_zero(&mut self, _q: usize) -> bool {
-        // Because `qubit_is_zero` is called on every qubit release, this must return
-        // true to avoid a panic.
+        // We don't simulate quantum execution here. So we don't know if the qubit
+        // is zero or not. Returning true avoids potential panics.
         true
     }
 

Diff for: compiler/qsc_eval/src/backend.rs

@@ -1,12 +1,15 @@
 // Copyright (c) Microsoft Corporation.
 // Licensed under the MIT License.
 
+use crate::noise::PauliNoise;
+use crate::val::Value;
 use num_bigint::BigUint;
 use num_complex::Complex;
 use quantum_sparse_sim::QuantumSim;
-use rand::RngCore;
+use rand::{rngs::StdRng, Rng, RngCore, SeedableRng};
 
-use crate::val::Value;
+#[cfg(test)]
+mod noise_tests;
 
 /// The trait that must be implemented by a quantum backend, whose functions will be invoked when
 /// quantum intrinsics are called.
@@ -82,7 +85,11 @@ pub trait Backend {
     fn qubit_allocate(&mut self) -> usize {
         unimplemented!("qubit_allocate operation");
     }
-    fn qubit_release(&mut self, _q: usize) {
+    /// `false` indicates that the qubit was in a non-zero state before the release,
+    /// but should have been in the zero state.
+    /// `true` otherwise. This includes the case when the qubit was in
+    /// a non-zero state during a noisy simulation, which is allowed.
+    fn qubit_release(&mut self, _q: usize) -> bool {
         unimplemented!("qubit_release operation");
     }
     fn qubit_swap_id(&mut self, _q0: usize, _q1: usize) {
@@ -102,7 +109,14 @@ pub trait Backend {
 
 /// Default backend used when targeting sparse simulation.
 pub struct SparseSim {
+    /// Noiseless Sparse simulator to be used by this instance.
     pub sim: QuantumSim,
+    /// Pauli noise that is applied after a gate or before a measurement is executed.
+    /// Service functions aren't subject to noise.
+    pub noise: PauliNoise,
+    /// Random number generator to sample Pauli noise.
+    /// Noise is not applied when rng is None.
+    pub rng: Option<StdRng>,
 }
 
 impl Default for SparseSim {
@@ -116,132 +130,213 @@ impl SparseSim {
     pub fn new() -> Self {
         Self {
             sim: QuantumSim::new(None),
+            noise: PauliNoise::default(),
+            rng: None,
+        }
+    }
+
+    #[must_use]
+    pub fn new_with_noise(noise: &PauliNoise) -> Self {
+        Self {
+            sim: QuantumSim::new(None),
+            noise: *noise,
+            rng: if noise.is_noiseless() {
+                None
+            } else {
+                Some(StdRng::from_entropy())
+            },
         }
     }
+
+    #[must_use]
+    fn is_noiseless(&self) -> bool {
+        self.rng.is_none()
+    }
+
+    fn apply_noise(&mut self, q: usize) {
+        if let Some(rng) = &mut self.rng {
+            let p = rng.gen_range(0.0..1.0);
+            if p >= self.noise.distribution[2] {
+                // In the most common case we don't apply noise
+            } else if p < self.noise.distribution[0] {
+                self.sim.x(q);
+            } else if p < self.noise.distribution[1] {
+                self.sim.y(q);
+            } else {
+                self.sim.z(q);
+            }
+        }
+        // No noise applied if rng is None.
+    }
 }
 
 impl Backend for SparseSim {
     type ResultType = bool;
 
     fn ccx(&mut self, ctl0: usize, ctl1: usize, q: usize) {
         self.sim.mcx(&[ctl0, ctl1], q);
+        self.apply_noise(ctl0);
+        self.apply_noise(ctl1);
+        self.apply_noise(q);
     }
 
     fn cx(&mut self, ctl: usize, q: usize) {
         self.sim.mcx(&[ctl], q);
+        self.apply_noise(ctl);
+        self.apply_noise(q);
     }
 
     fn cy(&mut self, ctl: usize, q: usize) {
         self.sim.mcy(&[ctl], q);
+        self.apply_noise(ctl);
+        self.apply_noise(q);
     }
 
     fn cz(&mut self, ctl: usize, q: usize) {
         self.sim.mcz(&[ctl], q);
+        self.apply_noise(ctl);
+        self.apply_noise(q);
     }
 
     fn h(&mut self, q: usize) {
         self.sim.h(q);
+        self.apply_noise(q);
     }
 
     fn m(&mut self, q: usize) -> Self::ResultType {
+        self.apply_noise(q);
         self.sim.measure(q)
     }
 
     fn mresetz(&mut self, q: usize) -> Self::ResultType {
+        self.apply_noise(q); // Applying noise before measurement
         let res = self.sim.measure(q);
         if res {
             self.sim.x(q);
         }
+        self.apply_noise(q); // Applying noise after reset
         res
     }
 
     fn reset(&mut self, q: usize) {
         self.mresetz(q);
+        // Noise applied in mresetz.
     }
 
     fn rx(&mut self, theta: f64, q: usize) {
         self.sim.rx(theta, q);
+        self.apply_noise(q);
     }
 
     fn rxx(&mut self, theta: f64, q0: usize, q1: usize) {
-        self.h(q0);
-        self.h(q1);
-        self.rzz(theta, q0, q1);
-        self.h(q1);
-        self.h(q0);
+        self.sim.h(q0);
+        self.sim.h(q1);
+        self.sim.mcx(&[q1], q0);
+        self.sim.rz(theta, q0);
+        self.sim.mcx(&[q1], q0);
+        self.sim.h(q1);
+        self.sim.h(q0);
+        self.apply_noise(q0);
+        self.apply_noise(q1);
     }
 
     fn ry(&mut self, theta: f64, q: usize) {
         self.sim.ry(theta, q);
+        self.apply_noise(q);
     }
 
     fn ryy(&mut self, theta: f64, q0: usize, q1: usize) {
-        self.h(q0);
-        self.s(q0);
-        self.h(q0);
-        self.h(q1);
-        self.s(q1);
-        self.h(q1);
-        self.rzz(theta, q0, q1);
-        self.h(q1);
-        self.sadj(q1);
-        self.h(q1);
-        self.h(q0);
-        self.sadj(q0);
-        self.h(q0);
+        self.sim.h(q0);
+        self.sim.s(q0);
+        self.sim.h(q0);
+        self.sim.h(q1);
+        self.sim.s(q1);
+        self.sim.h(q1);
+        self.sim.mcx(&[q1], q0);
+        self.sim.rz(theta, q0);
+        self.sim.mcx(&[q1], q0);
+        self.sim.h(q1);
+        self.sim.sadj(q1);
+        self.sim.h(q1);
+        self.sim.h(q0);
+        self.sim.sadj(q0);
+        self.sim.h(q0);
+        self.apply_noise(q0);
+        self.apply_noise(q1);
     }
 
     fn rz(&mut self, theta: f64, q: usize) {
         self.sim.rz(theta, q);
+        self.apply_noise(q);
     }
 
     fn rzz(&mut self, theta: f64, q0: usize, q1: usize) {
-        self.cx(q1, q0);
-        self.rz(theta, q0);
-        self.cx(q1, q0);
+        self.sim.mcx(&[q1], q0);
+        self.sim.rz(theta, q0);
+        self.sim.mcx(&[q1], q0);
+        self.apply_noise(q0);
+        self.apply_noise(q1);
     }
 
     fn sadj(&mut self, q: usize) {
         self.sim.sadj(q);
+        self.apply_noise(q);
     }
 
     fn s(&mut self, q: usize) {
         self.sim.s(q);
+        self.apply_noise(q);
     }
 
     fn swap(&mut self, q0: usize, q1: usize) {
         self.sim.swap_qubit_ids(q0, q1);
+        self.apply_noise(q0);
+        self.apply_noise(q1);
     }
 
     fn tadj(&mut self, q: usize) {
         self.sim.tadj(q);
+        self.apply_noise(q);
     }
 
     fn t(&mut self, q: usize) {
         self.sim.t(q);
+        self.apply_noise(q);
     }
 
     fn x(&mut self, q: usize) {
         self.sim.x(q);
+        self.apply_noise(q);
     }
 
     fn y(&mut self, q: usize) {
         self.sim.y(q);
+        self.apply_noise(q);
     }
 
     fn z(&mut self, q: usize) {
         self.sim.z(q);
+        self.apply_noise(q);
     }
 
     fn qubit_allocate(&mut self) -> usize {
+        // Fresh qubit start in ground state even with noise.
         self.sim.allocate()
     }
 
-    fn qubit_release(&mut self, q: usize) {
-        self.sim.release(q);
+    fn qubit_release(&mut self, q: usize) -> bool {
+        if self.is_noiseless() {
+            let was_zero = self.sim.qubit_is_zero(q);
+            self.sim.release(q);
+            was_zero
+        } else {
+            self.sim.release(q);
+            true
+        }
     }
 
     fn qubit_swap_id(&mut self, q0: usize, q1: usize) {
+        // This is a service function rather than a gate so it doesn't incur noise.
         self.sim.swap_qubit_ids(q0, q1);
     }
 
@@ -266,10 +361,12 @@ impl Backend for SparseSim {
     }
 
     fn qubit_is_zero(&mut self, q: usize) -> bool {
+        // This is a service function rather than a measurement so it doesn't incur noise.
         self.sim.qubit_is_zero(q)
     }
 
     fn custom_intrinsic(&mut self, name: &str, arg: Value) -> Option<Result<Value, String>> {
+        // These intrinsics aren't subject to noise.
         match name {
             "GlobalPhase" => {
                 // Apply a global phase to the simulation by doing an Rz to a fresh qubit.
@@ -301,9 +398,16 @@ impl Backend for SparseSim {
     }
 
     fn set_seed(&mut self, seed: Option<u64>) {
-        match seed {
-            Some(seed) => self.sim.set_rng_seed(seed),
-            None => self.sim.set_rng_seed(rand::thread_rng().next_u64()),
+        if let Some(seed) = seed {
+            if !self.is_noiseless() {
+                self.rng = Some(StdRng::seed_from_u64(seed));
+            }
+            self.sim.set_rng_seed(seed);
+        } else {
+            if !self.is_noiseless() {
+                self.rng = Some(StdRng::from_entropy());
+            }
+            self.sim.set_rng_seed(rand::thread_rng().next_u64());
         }
     }
 }
@@ -458,9 +562,9 @@ where
         self.main.qubit_allocate()
     }
 
-    fn qubit_release(&mut self, q: usize) {
-        self.chained.qubit_release(q);
-        self.main.qubit_release(q);
+    fn qubit_release(&mut self, q: usize) -> bool {
+        let _ = self.chained.qubit_release(q);
+        self.main.qubit_release(q)
     }
 
     fn qubit_swap_id(&mut self, q0: usize, q1: usize) {