Update: improve gate application efficiency

Cargo.toml

@@ -12,6 +12,9 @@ path = "src/main.rs"
 name = "simquil"
 path = "src/lib.rs"
 
+[profile.release]
+debug = true
+
 [dependencies]
 structopt = "0.3.13"
 quil = { git = "https://github.com/rigetti/quil-rust.git", branch = "derive-error-for-evaluation-error" }

src/gates/mod.rs

@@ -7,8 +7,8 @@ use crate::gates::standard::swap;
 
 #[derive(Default, Clone, Debug)]
 pub struct QGate {
-    matrix: Array2<Complex64>,
-    qubits: Vec<usize>,
+    pub matrix: Array2<Complex64>,
+    pub qubits: Vec<usize>,
 }
 
 impl QGate {

src/wavefunction.rs

@@ -1,5 +1,7 @@
-use num::{complex::Complex64, pow};
 use std::fmt;
+
+use ndarray::Array;
+use num::{complex::Complex64, pow};
 use thiserror::Error;
 
 use crate::{
@@ -12,8 +14,8 @@ pub enum WavefunctionError {}
 
 #[derive(PartialEq, Default)]
 pub struct Wavefunction {
-    wfn: ndarray::Array1<Complex64>,
-    n_qubits: usize,
+    pub wfn: ndarray::Array1<Complex64>,
+    pub n_qubits: usize,
 }
 
 impl Wavefunction {
@@ -69,8 +71,7 @@ impl Wavefunction {
     /// # Arguments
     ///
     /// * `qubit` - The measured qubit
-    /// * `excited_probability` - The probability that `qubit` is in an excited
-    ///   state in the wavefunction
+    /// * `excited_probability` - The probability that `qubit` is in an excited state
     /// * `measured` - The result of the measurement (`0` or `1`)
     pub fn collapse_wavefunction(
         &mut self,
@@ -93,12 +94,75 @@ impl Wavefunction {
         }
     }
 
-    /// Apply a gate matrix to the wavefunction
+    /// Apply a gate to the wavefunction
     pub fn apply(&mut self, gate: &QGate) {
-        self.wfn = gate.lift(self.n_qubits).dot(&self.wfn)
+        let subwfns_indices = get_indices_for_qubits(gate.qubits.clone(), self.n_qubits);
+        for subwfn_indices in subwfns_indices {
+            let subwfn_amplitudes = Array::from_vec(
+                subwfn_indices
+                    .iter()
+                    .map(|i| self.wfn[*i])
+                    .collect::<Vec<_>>(),
+            );
+            let subwfn_new_amplitudes = gate.matrix.dot(&subwfn_amplitudes);
+            subwfn_indices
+                .iter()
+                .zip(subwfn_new_amplitudes)
+                .for_each(|(i, a)| self.wfn[*i] = a);
+        }
     }
 }
 
+/// Generate all bit strings (as integers up to 2^n with each bitstring having the restriction that for each qubit index
+/// in `qubits` the corresponding value in the bitstring is specified in `restrictions`.
+///
+/// For example, `bitstrings_with_restrictions(vec![2, 0], vec![1, 0], 3)` produces all 2^3 bitstrings (as integers)
+/// where the third bit is always 1 and the first bit is always 0, i.e [0b100, 0b110].
+fn bitstrings_with_restrictions(
+    qubits: Vec<usize>,
+    restrictions: Vec<usize>,
+    n_qubits: usize,
+) -> Vec<usize> {
+    (0..(2u32.pow(n_qubits as u32)))
+        .filter(|i| {
+            qubits
+                .iter()
+                .zip(&restrictions)
+                .all(|(q, v)| ((i >> *q) & 0b1) == (*v as u32))
+        })
+        .map(|i| i as usize)
+        .collect()
+}
+
+/// Generate all bitstrings up to 2^n with each bitstring represented by a vector of bits.
+fn bitstrings(n: usize) -> Vec<Vec<usize>> {
+    (0..(2u32.pow(n as u32)))
+        .map(|v| (0..n).rev().map(|i| ((v >> i) & 1) as usize).collect())
+        .collect()
+}
+
+fn get_indices_for_qubits(qubits: Vec<usize>, n_qubits: usize) -> Vec<Vec<usize>> {
+    let b = bitstrings(qubits.len())
+        .iter()
+        .map(|r| bitstrings_with_restrictions(qubits.clone(), r.clone(), n_qubits))
+        .collect::<Vec<_>>();
+
+    let mut results: Vec<Vec<usize>> = vec![];
+    let rows = b.len();
+    let cols = b[0].len();
+    // Essentially a matrix transpose
+    for j in 0..cols {
+        let mut new_row = vec![];
+        #[allow(clippy::needless_range_loop)]
+        for i in 0..rows {
+            new_row.push(b[i][j])
+        }
+        results.push(new_row);
+    }
+
+    results
+}
+
 impl fmt::Debug for Wavefunction {
     fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
         self.wfn.iter().enumerate().for_each(|(i, amplitude)| {
@@ -123,12 +187,13 @@ impl fmt::Debug for Wavefunction {
 mod tests {
     use ndarray::arr1;
 
-    use super::Wavefunction;
     use crate::{
         gates,
         matrix::{C0, C1},
     };
 
+    use super::Wavefunction;
+
     #[test]
     fn apply_1q_gates() {
         let mut wfn = Wavefunction::ground_state_wavefunction(1);