Implement Display and conj. transpose for Unitary

- Add method to multiply unitary by another unitary

- Add tests to make sure multiplying a unitary, u, against its conjugate
  transpose yields the identity

- Fix bug in conj_t method that didn't account for 2D to 1D mapping

CHANGELOG.md

@@ -2,6 +2,10 @@
 
 ## 0.5.0
 
+### (Minor) Breaking changes
+*Since Spinoza is pre `v1.0.0`, and under heavy ...*
+
+
 ### New features
 
 - Add the controlled Hadamard gate.

spinoza/src/circuit.rs

@@ -201,6 +201,7 @@ impl QuantumCircuit {
         &self.state
     }
 
+    // TODO(saveliy): look into fusing the two loops
     /// Invert this circuit
     pub fn inverse(&mut self) {
         self.transformations.reverse();

spinoza/src/gates.rs

@@ -84,7 +84,10 @@ impl Gate {
             Self::RZ(theta) => Self::RZ(-theta),
             Self::U(theta, phi, lambda) => Self::U(-theta, -lambda, -phi),
             Self::M | Self::BitFlipNoise(_) => unimplemented!(),
-            Self::Unitary(_) => todo!(),
+            Self::Unitary(mut unitary) => {
+                unitary.conj_t();
+                Self::Unitary(unitary)
+            }
         }
     }
 
@@ -1300,6 +1303,7 @@ fn unitary_col_vec_mul(
     let mut reals = Vec::with_capacity(vec_reals.len());
     let mut imags = Vec::with_capacity(vec_imags.len());
 
+    // TODO(saveliy): parallelize this
     unitary
         .reals
         .chunks_exact(chunk_size)
@@ -1945,6 +1949,9 @@ mod tests {
         assert_float_closeness(identity[3].im, 0.0, 0.001);
     }
 
+    #[test]
+    fn unitary_inverse() {}
+
     #[test]
     #[should_panic]
     fn m_inverse() {

spinoza/src/unitaries.rs

@@ -1,6 +1,9 @@
 //! Functionality for applying large 2^n * 2^n matrices to the state
 //! Ideally, this should be a last resort
 use crate::{core::State, gates::Gate, math::Float};
+use rayon::prelude::*;
+use std::fmt;
+use std::fmt::Formatter;
 
 /// A representation of a Unitary Matrix
 #[derive(Clone)]
@@ -13,7 +16,23 @@ pub struct Unitary {
     pub width: usize,
 }
 
+impl fmt::Display for Unitary {
+    fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
+        self.reals
+            .chunks_exact(self.width)
+            .zip(self.imags.chunks_exact(self.width))
+            .for_each(|(re, im)| {
+                re.iter().zip(im.iter()).for_each(|(z_re, z_im)| {
+                    write!(f, "{z_re}+i{z_im} ").unwrap();
+                });
+                writeln!(f).unwrap();
+            });
+        Ok(())
+    }
+}
+
 impl Unitary {
+    // TODO(saveliy): look into using strategy 2 here
     /// Construct a Unitary from a single qubit gate
     pub fn from_single_qubit_gate(state: &State, gate: Gate, target: usize) -> Self {
         let g = gate.to_matrix();
@@ -47,6 +66,65 @@ impl Unitary {
             width,
         }
     }
+
+    /// Return the conjugate transpose of this `Unitary`, in-place
+    pub fn conj_t(&mut self) {
+        self.imags.par_iter_mut().for_each(|z_im| {
+            *z_im = -(*z_im);
+        });
+
+        for i in 0..self.height {
+            for j in i + 1..self.width {
+                self.reals.swap(i * self.width + j, j * self.width + i);
+                self.imags.swap(i * self.width + j, j * self.width + i);
+            }
+        }
+    }
+
+    /// Multiply this unitary matrix by another unitary matrix
+    pub fn multiply(&self, other: &Unitary) -> Self {
+        assert_eq!(self.width, other.height);
+
+        let mut result_reals = vec![0.0; self.height * other.width];
+        let mut result_imags = vec![0.0; self.height * other.width];
+
+        // self.reals
+        //     .chunks_exact(self.width)
+        //     .zip(self.imags.chunks_exact(self.width))
+        //     .enumerate()
+        //     .map(|(i, (z_re, z_im))| {
+        //         for k in 0..other.height {
+        //
+        //         }
+        //     });
+
+        for i in 0..self.height {
+            for j in 0..other.width {
+                let mut sum_real = 0.0;
+                let mut sum_imag = 0.0;
+
+                for k in 0..self.width {
+                    let a_real = self.reals[i * self.width + k];
+                    let a_imag = self.imags[i * self.width + k];
+                    let b_real = other.reals[k * other.width + j];
+                    let b_imag = other.imags[k * other.width + j];
+
+                    sum_real += a_real * b_real - a_imag * b_imag;
+                    sum_imag += a_real * b_imag + a_imag * b_real;
+                }
+
+                result_reals[i * other.width + j] = sum_real;
+                result_imags[i * other.width + j] = sum_imag;
+            }
+        }
+
+        Unitary {
+            reals: result_reals,
+            imags: result_imags,
+            height: self.height,
+            width: other.width,
+        }
+    }
 }
 
 /// Applies a unitary matrix to the Quantum State Vector
@@ -85,13 +163,13 @@ pub fn apply_unitary(state: &State, unitary: &Unitary) -> State {
 
 #[cfg(test)]
 mod tests {
+    use super::*;
+    use crate::utils::{assert_float_closeness, gen_random_state};
     use crate::{
         gates::{apply, Gate},
         math::SQRT_ONE_HALF,
     };
 
-    use super::*;
-
     #[test]
     fn test_hxi_from_single_qubit_gate() {
         let mut reals = vec![
@@ -259,4 +337,38 @@ mod tests {
         assert_eq!(s.reals, s1.reals);
         assert_eq!(s.imags, s1.imags);
     }
+
+    #[test]
+    fn display() {
+        const N: usize = 1;
+        let state = gen_random_state(N);
+        let u = Unitary::from_single_qubit_gate(&state, Gate::X, 0);
+
+        let x_gate_as_str = "0+i0 1+i0 \n1+i0 0+i0 \n".to_string();
+        let u_as_str = u.to_string();
+        assert_eq!(u_as_str, x_gate_as_str);
+    }
+
+    #[test]
+    fn conjugate_transpose() {
+        const N: usize = 2;
+        let state = State::new(N);
+        let u = Unitary::from_single_qubit_gate(&state, Gate::H, 0);
+        let mut u_ct = Unitary::from_single_qubit_gate(&state, Gate::H, 0);
+        u_ct.conj_t();
+
+        let identity = u_ct.multiply(&u);
+
+        for i in 0..u.height {
+            for j in 0..u.width {
+                if i == j {
+                    assert_float_closeness(identity.reals[i * identity.width + j], 1.0, 0.00001);
+                    assert_float_closeness(identity.imags[i * identity.width + j], 0.0, 0.00001);
+                } else {
+                    assert_float_closeness(identity.reals[i * identity.width + j], 0.0, 0.00001);
+                    assert_float_closeness(identity.imags[i * identity.width + j], 0.0, 0.00001);
+                }
+            }
+        }
+    }
 }