Support ApplyUnitary in simulation (#2051)

This adds a new library function `ApplyUnitary` that given a matrix of
complex numbers will apply that matrix to the given qubits. This
requires a new version of the simulator that exposes the functionality
publicly.

Starting as draft to wait for
qir-alliance/qir-runner#208

Cargo.toml

@@ -58,6 +58,7 @@ log = "0.4"
 miette = { version = "7.2", features = ["fancy-no-syscall"] }
 thiserror = "1.0"
 nalgebra = { version = "0.33" }
+ndarray = "0.15.4"
 num-bigint = "0.4"
 num-complex = "0.4"
 num-traits = "0.2"
@@ -77,7 +78,7 @@ wasm-bindgen-futures = "0.4"
 rand = "0.8"
 serde_json = "1.0"
 pyo3 = "0.22"
-quantum-sparse-sim = { git = "https://github.com/qir-alliance/qir-runner", tag = "v0.7.4" }
+quantum-sparse-sim = { git = "https://github.com/qir-alliance/qir-runner", rev = "562e2c11ad685dd01bfc1ae975e00d4133615995" }
 async-trait = "0.1"
 tokio = { version = "1.35", features = ["macros", "rt"] }
 

compiler/qsc_eval/Cargo.toml

@@ -10,6 +10,7 @@ license.workspace = true
 
 [dependencies]
 miette = { workspace = true }
+ndarray = { workspace = true }
 num-bigint = { workspace = true }
 num-complex = { workspace = true }
 num-traits = { workspace = true }

compiler/qsc_eval/src/backend.rs

@@ -1,8 +1,9 @@
 // Copyright (c) Microsoft Corporation.
 // Licensed under the MIT License.
 
-use crate::noise::PauliNoise;
 use crate::val::Value;
+use crate::{noise::PauliNoise, val::unwrap_tuple};
+use ndarray::Array2;
 use num_bigint::BigUint;
 use num_complex::Complex;
 use quantum_sparse_sim::QuantumSim;
@@ -419,6 +420,29 @@ impl Backend for SparseSim {
                 self.apply_noise(q);
                 Some(Ok(Value::unit()))
             }
+            "Apply" => {
+                let [matrix, qubits] = unwrap_tuple(arg);
+                let qubits = qubits
+                    .unwrap_array()
+                    .iter()
+                    .filter_map(|q| q.clone().unwrap_qubit().try_deref().map(|q| q.0))
+                    .collect::<Vec<_>>();
+                let matrix = unwrap_matrix_as_array2(matrix, &qubits);
+
+                // Confirm the matrix is unitary by checking if multiplying it by its adjoint gives the identity matrix (up to numerical precision).
+                let adj = matrix.t().map(Complex::<f64>::conj);
+                if (matrix.dot(&adj) - Array2::<Complex<f64>>::eye(1 << qubits.len()))
+                    .map(|x| x.norm())
+                    .sum()
+                    > 1e-9
+                {
+                    return Some(Err("matrix is not unitary".to_string()));
+                }
+
+                self.sim.apply(&matrix, &qubits, None);
+
+                Some(Ok(Value::unit()))
+            }
             _ => None,
         }
     }
@@ -438,6 +462,27 @@ impl Backend for SparseSim {
     }
 }
 
+fn unwrap_matrix_as_array2(matrix: Value, qubits: &[usize]) -> Array2<Complex<f64>> {
+    let matrix: Vec<Vec<Complex<f64>>> = matrix
+        .unwrap_array()
+        .iter()
+        .map(|row| {
+            row.clone()
+                .unwrap_array()
+                .iter()
+                .map(|elem| {
+                    let [re, im] = unwrap_tuple(elem.clone());
+                    Complex::<f64>::new(re.unwrap_double(), im.unwrap_double())
+                })
+                .collect::<Vec<_>>()
+        })
+        .collect::<Vec<_>>();
+
+    Array2::from_shape_fn((1 << qubits.len(), 1 << qubits.len()), |(i, j)| {
+        matrix[i][j]
+    })
+}
+
 /// Simple struct that chains two backends together so that the chained
 /// backend is called before the main backend.
 /// For any intrinsics that return a value,

compiler/qsc_eval/src/intrinsic.rs

@@ -10,13 +10,12 @@ use crate::{
     backend::Backend,
     error::PackageSpan,
     output::Receiver,
-    val::{self, Value},
+    val::{self, unwrap_tuple, Value},
     Error, Rc,
 };
 use num_bigint::BigInt;
 use rand::{rngs::StdRng, Rng};
 use rustc_hash::{FxHashMap, FxHashSet};
-use std::array;
 use std::convert::TryFrom;
 
 #[allow(clippy::too_many_lines)]
@@ -426,8 +425,3 @@ pub fn qubit_relabel(
 
     Ok(Value::unit())
 }
-
-fn unwrap_tuple<const N: usize>(value: Value) -> [Value; N] {
-    let values = value.unwrap_tuple();
-    array::from_fn(|i| values[i].clone())
-}

compiler/qsc_eval/src/val.rs

@@ -5,6 +5,7 @@ use num_bigint::BigInt;
 use qsc_data_structures::{display::join, functors::FunctorApp};
 use qsc_fir::fir::{Functor, Pauli, StoreItemId};
 use std::{
+    array,
     fmt::{self, Display, Formatter},
     rc::{Rc, Weak},
 };
@@ -578,3 +579,9 @@ pub fn update_functor_app(functor: Functor, app: FunctorApp) -> FunctorApp {
         },
     }
 }
+
+#[must_use]
+pub fn unwrap_tuple<const N: usize>(value: Value) -> [Value; N] {
+    let values = value.unwrap_tuple();
+    array::from_fn(|i| values[i].clone())
+}

library/src/tests.rs

@@ -16,7 +16,7 @@ mod table_lookup;
 
 use indoc::indoc;
 use qsc::{
-    interpret::{GenericReceiver, Interpreter, Result, Value},
+    interpret::{self, GenericReceiver, Interpreter, Result, Value},
     target::Profile,
     Backend, LanguageFeatures, PackageType, SourceMap, SparseSim,
 };
@@ -30,6 +30,15 @@ pub fn test_expression(expr: &str, expected: &Value) -> String {
     test_expression_with_lib(expr, "", expected)
 }
 
+pub fn test_expression_fails(expr: &str) -> String {
+    test_expression_fails_with_lib_and_profile_and_sim(
+        expr,
+        "",
+        Profile::Unrestricted,
+        &mut SparseSim::default(),
+    )
+}
+
 pub fn test_expression_with_lib(expr: &str, lib: &str, expected: &Value) -> String {
     test_expression_with_lib_and_profile(expr, lib, Profile::Unrestricted, expected)
 }
@@ -93,6 +102,45 @@ pub fn test_expression_with_lib_and_profile_and_sim(
     String::from_utf8(stdout).expect("stdout should be valid utf8")
 }
 
+pub fn test_expression_fails_with_lib_and_profile_and_sim(
+    expr: &str,
+    lib: &str,
+    profile: Profile,
+    sim: &mut impl Backend<ResultType = impl Into<Result>>,
+) -> String {
+    let mut stdout = vec![];
+    let mut out = GenericReceiver::new(&mut stdout);
+
+    let sources = SourceMap::new([("test".into(), lib.into())], Some(expr.into()));
+
+    let (std_id, store) = qsc::compile::package_store_with_stdlib(profile.into());
+
+    let mut interpreter = Interpreter::new(
+        sources,
+        PackageType::Exe,
+        profile.into(),
+        LanguageFeatures::default(),
+        store,
+        &[(std_id, None)],
+    )
+    .expect("test should compile");
+
+    let result = interpreter
+        .eval_entry_with_sim(sim, &mut out)
+        .expect_err("test should run successfully");
+
+    assert!(
+        result.len() == 1,
+        "Expected a single error, got {:?}",
+        result.len()
+    );
+    let interpret::Error::Eval(err) = &result[0] else {
+        panic!("Expected an Eval error, got {:?}", result[0]);
+    };
+
+    err.error().error().to_string()
+}
+
 /// # Panics
 ///
 /// Will panic if f64 values are significantly different.

library/src/tests/intrinsic.rs

@@ -7,7 +7,7 @@ use expect_test::expect;
 use indoc::indoc;
 use qsc::{interpret::Value, target::Profile, SparseSim};
 
-use super::{test_expression, test_expression_with_lib_and_profile_and_sim};
+use super::{test_expression, test_expression_fails, test_expression_with_lib_and_profile_and_sim};
 
 // These tests verify multi-controlled decomposition logic for gate operations. Each test
 // manually allocates 2N qubits, performs the decomposed operation from the library on the first N,
@@ -3008,3 +3008,142 @@ fn test_exp() {
     "#]]
     .assert_eq(&dump);
 }
+
+#[test]
+fn test_apply_unitary_with_h_matrix() {
+    let dump = test_expression(
+        indoc! {"
+        {
+            open Std.Math;
+            open Std.Diagnostics;
+            use q = Qubit();
+            let one_sqrt_2 = new Complex { Real = 1.0 / Sqrt(2.0), Imag = 0.0 };
+            ApplyUnitary(
+                [
+                    [one_sqrt_2, one_sqrt_2],
+                    [one_sqrt_2, NegationC(one_sqrt_2)]
+                ],
+                [q]
+            );
+            DumpMachine();
+            Reset(q);
+        }
+        "},
+        &Value::unit(),
+    );
+
+    expect![[r#"
+        STATE:
+        |0⟩: 0.7071+0.0000𝑖
+        |1⟩: 0.7071+0.0000𝑖
+    "#]]
+    .assert_eq(&dump);
+}
+
+#[test]
+fn test_apply_unitary_with_swap_matrix() {
+    let dump = test_expression(
+        indoc! {"
+        {
+            open Std.Math;
+            open Std.Diagnostics;
+            use qs = Qubit[2];
+            H(qs[0]);
+            DumpMachine();
+            let one = new Complex { Real = 1.0, Imag = 0.0 };
+            let zero = new Complex { Real = 0.0, Imag = 0.0 };
+            ApplyUnitary(
+                [
+                    [one, zero, zero, zero],
+                    [zero, zero, one, zero],
+                    [zero, one, zero, zero],
+                    [zero, zero, zero, one]
+                ],
+                qs
+            );
+            DumpMachine();
+            ResetAll(qs);
+        }
+        "},
+        &Value::unit(),
+    );
+
+    expect![[r#"
+        STATE:
+        |00⟩: 0.7071+0.0000𝑖
+        |10⟩: 0.7071+0.0000𝑖
+        STATE:
+        |00⟩: 0.7071+0.0000𝑖
+        |01⟩: 0.7071+0.0000𝑖
+    "#]]
+    .assert_eq(&dump);
+}
+
+#[test]
+fn test_apply_unitary_fails_when_matrix_not_square() {
+    let err = test_expression_fails(indoc! {"
+        {
+            open Std.Math;
+            open Std.Diagnostics;
+            use q = Qubit();
+            ApplyUnitary(
+                [
+                    [new Complex { Real = 1.0, Imag = 0.0 }],
+                    [new Complex { Real = 0.0, Imag = 0.0 }]
+                ],
+                [q]
+            );
+            DumpMachine();
+            Reset(q);
+        }
+        "});
+
+    expect!["program failed: matrix passed to ApplyUnitary must be square."].assert_eq(&err);
+}
+
+#[test]
+fn test_apply_unitary_fails_when_matrix_wrong_size() {
+    let err = test_expression_fails(indoc! {"
+        {
+            open Std.Math;
+            open Std.Diagnostics;
+            use qs = Qubit[2];
+            let one_sqrt_2 = new Complex { Real = 1.0 / Sqrt(2.0), Imag = 0.0 };
+            ApplyUnitary(
+                [
+                    [one_sqrt_2, one_sqrt_2],
+                    [one_sqrt_2, NegationC(one_sqrt_2)]
+                ],
+                qs
+            );
+            DumpMachine();
+            ResetAll(qs);
+        }
+        "});
+
+    expect!["program failed: matrix passed to ApplyUnitary must have dimensions 2^Length(qubits)."]
+        .assert_eq(&err);
+}
+
+#[test]
+fn test_apply_unitary_fails_when_matrix_not_unitary() {
+    let err = test_expression_fails(indoc! {"
+        {
+            open Std.Math;
+            open Std.Diagnostics;
+            use q = Qubit();
+            let zero = new Complex { Real = 0.0, Imag = 0.0 };
+            ApplyUnitary(
+                [
+                    [zero, zero],
+                    [zero, zero]
+                ],
+                [q]
+            );
+            DumpMachine();
+            Reset(q);
+        }
+        "});
+
+    expect!["intrinsic callable `Apply` failed: matrix is not unitary"].assert_eq(&err);
+}