From 645243b71f419594442b41360d2d9aa838b06a52 Mon Sep 17 00:00:00 2001 From: "Stefan J. Wernli" Date: Thu, 5 Dec 2024 10:57:43 -0800 Subject: [PATCH] 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 https://github.com/qir-alliance/qir-runner/pull/208 --- Cargo.lock | 5 +- Cargo.toml | 3 +- compiler/qsc_eval/Cargo.toml | 1 + compiler/qsc_eval/src/backend.rs | 47 +++++++++- compiler/qsc_eval/src/intrinsic.rs | 8 +- compiler/qsc_eval/src/val.rs | 7 ++ library/src/tests.rs | 50 +++++++++- library/src/tests/intrinsic.rs | 141 ++++++++++++++++++++++++++++- library/std/src/Std/Intrinsic.qs | 51 ++++++++++- 9 files changed, 299 insertions(+), 14 deletions(-) diff --git a/Cargo.lock b/Cargo.lock index 9fc4a8b0a9..493fc42ba7 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -1341,6 +1341,7 @@ dependencies = [ "expect-test", "indoc", "miette", + "ndarray", "num-bigint", "num-complex", "num-traits", @@ -1625,8 +1626,8 @@ dependencies = [ [[package]] name = "quantum-sparse-sim" -version = "0.7.4" -source = "git+https://github.com/qir-alliance/qir-runner?tag=v0.7.4#7b41f9313609f8309317ce91afccf0a2c7fe5a6f" +version = "0.7.5" +source = "git+https://github.com/qir-alliance/qir-runner?rev=562e2c11ad685dd01bfc1ae975e00d4133615995#562e2c11ad685dd01bfc1ae975e00d4133615995" dependencies = [ "ndarray", "num-bigint", diff --git a/Cargo.toml b/Cargo.toml index 0925ce8a6d..57aafe9faa 100644 --- a/Cargo.toml +++ b/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"] } diff --git a/compiler/qsc_eval/Cargo.toml b/compiler/qsc_eval/Cargo.toml index 366a1ea295..2dc458b16a 100644 --- a/compiler/qsc_eval/Cargo.toml +++ b/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 } diff --git a/compiler/qsc_eval/src/backend.rs b/compiler/qsc_eval/src/backend.rs index 6939e3b18a..62cd9021fe 100644 --- a/compiler/qsc_eval/src/backend.rs +++ b/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::>(); + 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::::conj); + if (matrix.dot(&adj) - Array2::>::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> { + let matrix: Vec>> = matrix + .unwrap_array() + .iter() + .map(|row| { + row.clone() + .unwrap_array() + .iter() + .map(|elem| { + let [re, im] = unwrap_tuple(elem.clone()); + Complex::::new(re.unwrap_double(), im.unwrap_double()) + }) + .collect::>() + }) + .collect::>(); + + 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, diff --git a/compiler/qsc_eval/src/intrinsic.rs b/compiler/qsc_eval/src/intrinsic.rs index 3ea85a14be..711f34de19 100644 --- a/compiler/qsc_eval/src/intrinsic.rs +++ b/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(value: Value) -> [Value; N] { - let values = value.unwrap_tuple(); - array::from_fn(|i| values[i].clone()) -} diff --git a/compiler/qsc_eval/src/val.rs b/compiler/qsc_eval/src/val.rs index af0000f311..f3d17f11ed 100644 --- a/compiler/qsc_eval/src/val.rs +++ b/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(value: Value) -> [Value; N] { + let values = value.unwrap_tuple(); + array::from_fn(|i| values[i].clone()) +} diff --git a/library/src/tests.rs b/library/src/tests.rs index 0d801d5c78..ec5aaf503a 100644 --- a/library/src/tests.rs +++ b/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>, +) -> 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. diff --git a/library/src/tests/intrinsic.rs b/library/src/tests/intrinsic.rs index c9ff5effa5..f58223c883 100644 --- a/library/src/tests/intrinsic.rs +++ b/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); +} diff --git a/library/std/src/Std/Intrinsic.qs b/library/std/src/Std/Intrinsic.qs index f6d6fa332d..67465050f9 100644 --- a/library/std/src/Std/Intrinsic.qs +++ b/library/std/src/Std/Intrinsic.qs @@ -1140,6 +1140,55 @@ operation Z(qubit : Qubit) : Unit is Adj + Ctl { adjoint self; } +/// # Summary +/// Applies the given unitary matrix to the given qubits. The matrix is checked at runtime to ensure it's shape is square and that the matrix dimensions are `2 ^ Length(qubits)`. +/// This operation is simulator-only and is not supported on hardware. +/// +/// # Input +/// ## matrix +/// The unitary matrix to apply. +/// ## qubits +/// The qubits to which the unitary matrix should be applied. +/// +/// # Example +/// This performs a two qubit CNOT using the unitary matrix representation: +/// ```qsharp +/// import Std.Math.Complex; +/// use qs = Qubit[2]; +/// 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, one, zero, zero], +/// [zero, zero, zero, one], +/// [zero, zero, one, zero] +/// ], +/// qs +/// ); +/// ``` +@Config(Unrestricted) +operation ApplyUnitary(matrix : Complex[][], qubits : Qubit[]) : Unit { + let num_rows = Length(matrix); + for col in matrix { + if Length(col) != num_rows { + fail "matrix passed to ApplyUnitary must be square."; + } + } + + let num_qubits = Length(qubits); + if num_rows != 1 <<< num_qubits { + fail "matrix passed to ApplyUnitary must have dimensions 2^Length(qubits)."; + } + + Apply(matrix, qubits); +} + +@Config(Unrestricted) +operation Apply(matrix : Complex[][], qubits : Qubit[]) : Unit { + body intrinsic; +} + /// # Summary /// Logs a message. /// @@ -1155,4 +1204,4 @@ function Message(msg : String) : Unit { body intrinsic; } -export AND, CCNOT, CNOT, Exp, H, I, M, Measure, R, R1, R1Frac, Reset, ResetAll, RFrac, Rx, Rxx, Ry, Ryy, Rz, Rzz, S, SWAP, T, X, Y, Z, Message; +export AND, CCNOT, CNOT, Exp, H, I, M, Measure, R, R1, R1Frac, Reset, ResetAll, RFrac, Rx, Rxx, Ry, Ryy, Rz, Rzz, S, SWAP, T, X, Y, Z, ApplyUnitary, Message;