diff --git a/src/Simulation/qdk_sim_rs/Cargo.toml b/src/Simulation/qdk_sim_rs/Cargo.toml index 3855ffe6a13..b6937572a9b 100644 --- a/src/Simulation/qdk_sim_rs/Cargo.toml +++ b/src/Simulation/qdk_sim_rs/Cargo.toml @@ -2,16 +2,19 @@ # Licensed under the MIT License. [package] -name = "qdk_sim_experimental" +name = "qdk_sim_rs" version = "0.0.1-alpha" authors = ["Microsoft"] edition = "2018" license = "MIT" -description = "Experimental simulators for use with the Quantum Development Kit." +description = "Rust-based simulators for use with the Quantum Development Kit." homepage = "https://github.com/microsoft/qsharp-runtime" repository = "https://github.com/microsoft/qsharp-runtime" readme = "README.md" +# Verified with cargo-msrv. +rust-version = "1.51.0" + exclude = [ # Exclude files specific to QDK build pipelines. "*.template", "*.csx", "*.ps1", "NuGet.Config", "drop", @@ -21,72 +24,85 @@ exclude = [ "*.egg-info", "qdk_sim_experimental", "setup.py", "*.whl", "pyproject.toml" ] +# Enable LaTeX on docs.rs. +# See https://stackoverflow.com/a/54573800/267841 and +# https://github.com/rust-num/num/pull/226/files for why this works. +[package.metadata.docs.rs] +rustdoc-args = [ "--html-in-header", "docs-includes/header.html", "--html-after-content", "docs-includes/after.html" ] +features = ["document-features"] + [lib] +crate-type = ["rlib", "staticlib", "cdylib"] name = "qdk_sim" path = "src/lib.rs" -crate-type = ["rlib", "staticlib", "cdylib"] # Optional build-time features: we use this to create Python and WASM bindings. [features] default = [] -wasm = ["web-sys"] + # When Python bindings are enabled, we also need the pyo3 dependency. +## Enables Python bindings for this crate. python = ["pyo3", "numpy"] -# Enable LaTeX on docs.rs. -# See https://stackoverflow.com/a/54573800/267841 and -# https://github.com/rust-num/num/pull/226/files for why this works. -[package.metadata.docs.rs] -rustdoc-args = [ "--html-in-header", "docs-includes/header.html", "--html-after-content", "docs-includes/after.html" ] +## Ensures that the crate is compatible with usage from WebAssembly. +wasm = ["web-sys"] [profile.release] -opt-level = 3 codegen-units = 1 # Reduce number of codegen units to increase optimizations. +opt-level = 3 panic = 'unwind' [dependencies] -ndarray = { version = "0.15.2", features = ["serde"] } -num-complex = { version = "0.4", features = ["serde"] } -num-traits = "0.2" -derive_more = "0.99.10" -itertools = "0.9.0" -rand = { version = "0.7.3", features = ["alloc"] } -serde = { version = "1.0", features = ["derive"] } -serde_json = "1.0" -lazy_static = "1.4.0" +anyhow = "1.0.56" +# We use built to expose compile-time metadata about how this crate +# was built to C and Rust callers. +built = "0.5.0" +cauchy = "0.4.0" cfg-if = "1.0.0" -num_enum = "0.5.1" +derive_more = "0.99.10" +# We use document-features to automatically generate feature documentation from +# Cargo.toml, following the example at +# https://docs.rs/document-features/latest/document_features#examples. +document-features = { version = "0.2", optional = true } # See https://github.com/rust-random/rand/issues/990 # and https://docs.rs/getrandom/0.1.15/getrandom/index.html#support-for-webassembly-and-asmjs # for why this is needed. # NB: We depend on 0.1.15, since that's what gets brought in transitively # by rand and rand_core. getrandom = { version = "0.1.15", features = ["wasm-bindgen"] } - -# We only need web-sys when compiling with the wasm feature. -web-sys = { version = "0.3.4", features = ['console'], optional = true } - +itertools = "0.9.0" +lazy_static = "1.4.0" +miette = "4.3.0" +ndarray = { version = "0.15.2", features = ["serde"] } +num-complex = { version = "0.4", features = ["serde"] } +num-traits = "0.2" +num_enum = "0.5.1" +numpy = { version = "0.13.1", optional = true } # We only need PyO3 when generating Python bindings. pyo3 = { version = "0.13.2", features = ["extension-module"], optional = true } -numpy = {version = "0.13.1", optional = true} +rand = { version = "0.7.3", features = ["alloc"] } +serde = { version = "1.0", features = ["derive"] } +serde_json = "1.0" +thiserror = "1.0.30" -# We use built to expose compile-time metadata about how this crate -# was built to C and Rust callers. -built = "0.5.0" +# We only need web-sys when compiling with the wasm feature. +web-sys = { version = "0.3.4", features = ['console'], optional = true } [build-dependencies] built = "0.5.0" [dev-dependencies] +approx = { version = "0.5.1", features = ["num-complex"] } assert-json-diff = "2.0.1" criterion = { version = "0.3", features = ['html_reports', 'csv_output'] } +ndarray = { version = "0.15.4", features = ["approx"] } [[bench]] -name = "c_api_benchmark" harness = false +name = "c_api_benchmark" [[bench]] -name = "microbenchmark" harness = false +name = "microbenchmark" diff --git a/src/Simulation/qdk_sim_rs/README.md b/src/Simulation/qdk_sim_rs/README.md index 0d88dd46c5a..2a2f7630100 100644 --- a/src/Simulation/qdk_sim_rs/README.md +++ b/src/Simulation/qdk_sim_rs/README.md @@ -4,7 +4,7 @@ To generate and view the documentation for this crate locally, please run: - $ cargo +nightly doc --features python --open + $ cargo doc --features python,document-features --open --> # Quantum Development Kit Preview Simulators @@ -28,11 +28,6 @@ The [`c_api`] module allows for using the simulation functionality in this crate Similarly, the [`python`] module allows exposing data structures in this crate to Python programs. -## Cargo Features - -- **`python`**: Enables Python bindings for this crate. -- **`wasm`**: Ensures that the crate is compatible with usage from WebAssembly. - ## Representing quantum systems This crate provides several different data structures for representing quantum systems in a variety of different conventions: @@ -148,3 +143,6 @@ TODO - Stabilizer simulation not yet exposed via C API. - Test and microbenchmark coverage still incomplete. - Too many APIs `panic!` or `unwrap`, and need replaced with `Result` returns instead. + +# Crate features + diff --git a/src/Simulation/qdk_sim_rs/benches/microbenchmark.rs b/src/Simulation/qdk_sim_rs/benches/microbenchmark.rs index 89ccd07b809..efecba8d185 100644 --- a/src/Simulation/qdk_sim_rs/benches/microbenchmark.rs +++ b/src/Simulation/qdk_sim_rs/benches/microbenchmark.rs @@ -8,11 +8,13 @@ //! In particular, optimizing these benchmarks may not translate into improved //! performance in user code. +use cauchy::c64; use criterion::{criterion_group, criterion_main, Criterion}; +use ndarray::array; use qdk_sim::{ common_matrices, common_matrices::nq_eye, - linalg::{extend_one_to_n, extend_two_to_n, Tensor}, + linalg::{extend_one_to_n, extend_two_to_n, Inv, Tensor}, }; fn linalg(c: &mut Criterion) { @@ -79,6 +81,68 @@ fn linalg(c: &mut Criterion) { let _result = cnot.tensor(&x); }) }); + group.bench_function("inv 4x4 f64", |b| { + let x = array![ + [ + 0.23935896217435304, + 0.34333031120985236, + 0.8201953415286973, + 0.8074588350909441 + ], + [ + 0.11957583380425751, + 0.16906445210054732, + 0.21728173861409317, + 0.7120594445167554 + ], + [ + 0.04023516190513021, + 0.9635112441739464, + 0.9209190516642924, + 0.114251355434274 + ], + [ + 0.8749507948480983, + 0.2661348079904513, + 0.17485566324545554, + 0.2934138616881069 + ], + ]; + b.iter(|| { + let _inv = x.inv(); + }); + }); + group.bench_function("inv 4x4 c64", |b| { + let x = array![ + [ + c64::new(0.30874277550419704, 0.8167808814398533), + c64::new(0.9303782008146939, 0.8925538040143673), + c64::new(0.11573522743286513, 0.6357551264716991), + c64::new(0.7869240102858357, 0.28376515716360073) + ], + [ + c64::new(0.9638410081049803, 0.4460520369459663), + c64::new(0.043516097874141346, 0.1652124014187376), + c64::new(0.05938096491956191, 0.7696366269843138), + c64::new(0.9636976605227736, 0.8125701401805293) + ], + [ + c64::new(0.9548859426476123, 0.7825350828251003), + c64::new(0.4223649577868721, 0.4522018603906839), + c64::new(0.36001119456757835, 0.22138920205104395), + c64::new(0.0044511389785256705, 0.14148562531641973) + ], + [ + c64::new(0.6066852151129799, 0.5140547256960247), + c64::new(0.23439110687939924, 0.3064074735518828), + c64::new(0.9963759728056912, 0.401859040365666), + c64::new(0.7176238235495955, 0.3948597214947167) + ], + ]; + b.iter(|| { + let _inv = x.inv(); + }); + }); group.finish(); } diff --git a/src/Simulation/qdk_sim_rs/src/c_api.rs b/src/Simulation/qdk_sim_rs/src/c_api.rs index caab540f118..70703cbf4b7 100644 --- a/src/Simulation/qdk_sim_rs/src/c_api.rs +++ b/src/Simulation/qdk_sim_rs/src/c_api.rs @@ -1,13 +1,9 @@ // Copyright (c) Microsoft Corporation. // Licensed under the MIT License. -// The following two attributes include the README.md for this module when -// building docs (requires +nightly). -// See https://github.com/rust-lang/rust/issues/82768#issuecomment-803935643 -// for discussion. -#![cfg_attr(doc, feature(extended_key_value_attributes))] -#![cfg_attr(doc, cfg_attr(doc, doc = include_str!("../docs/c-api.md")))] +#![cfg_attr(all(), doc = include_str!("../docs/c-api.md"))] +use crate::error::{QdkSimError, QdkSimError::*}; use crate::{built_info, NoiseModel, Process, State}; use lazy_static::lazy_static; use serde_json::json; @@ -32,12 +28,12 @@ lazy_static! { /// Exposes a result to C callers by setting LAST_ERROR in the Error /// case, and generating an appropriate error code. -fn as_capi_err Result<(), String>>(result_fn: F) -> i64 { +fn as_capi_err Result<(), QdkSimError>>(result_fn: F) -> i64 { let result = result_fn(); match result { Ok(_) => 0, - Err(msg) => { - *LAST_ERROR.lock().unwrap() = Some(msg); + Err(err) => { + *LAST_ERROR.lock().unwrap() = Some(err.to_string()); -1 } } @@ -47,7 +43,7 @@ fn apply &Process>( sim_id: usize, idxs: &[usize], channel_fn: F, -) -> Result<(), String> { +) -> Result<(), QdkSimError> { let state = &mut *STATE.lock().unwrap(); if let Some(sim_state) = state.get_mut(&sim_id) { let channel = channel_fn(&sim_state.noise_model); @@ -59,7 +55,10 @@ fn apply &Process>( Err(err) => Err(err), } } else { - return Err(format!("No simulator with id {}.", sim_id)); + Err(NoSuchSimulator { + invalid_id: sim_id, + expected: state.keys().into_iter().cloned().collect(), + }) } } @@ -121,11 +120,15 @@ pub unsafe extern "C" fn init( ) -> i64 { as_capi_err(|| { if representation.is_null() { - return Err("init called with null pointer for representation".to_string()); + return Err(NullPointer("representation".to_string())); } - let representation = CStr::from_ptr(representation) - .to_str() - .map_err(|e| format!("UTF-8 error decoding representation argument: {}", e))?; + let representation = + CStr::from_ptr(representation) + .to_str() + .map_err(|e| InvalidUtf8InArgument { + arg_name: "representation".to_string(), + source: e, + })?; let state = &mut *STATE.lock().unwrap(); let id = 1 + state.keys().fold(std::usize::MIN, |a, b| a.max(*b)); @@ -136,12 +139,7 @@ pub unsafe extern "C" fn init( "mixed" => State::new_mixed(initial_capacity), "pure" => State::new_pure(initial_capacity), "stabilizer" => State::new_stabilizer(initial_capacity), - _ => { - return Err(format!( - "Unknown initial state representation {}.", - representation - )) - } + _ => return Err(InvalidRepresentation(representation.to_string())), }, noise_model: NoiseModel::ideal(), }, @@ -161,7 +159,10 @@ pub extern "C" fn destroy(sim_id: usize) -> i64 { state.remove(&sim_id); Ok(()) } else { - Err(format!("No simulator with id {} exists.", sim_id)) + Err(NoSuchSimulator { + invalid_id: sim_id, + expected: state.keys().into_iter().cloned().collect(), + }) } }) } @@ -249,7 +250,10 @@ pub unsafe extern "C" fn m(sim_id: usize, idx: usize, result_out: *mut usize) -> *result_out = result; Ok(()) } else { - Err(format!("No simulator with id {} exists.", sim_id)) + Err(NoSuchSimulator { + invalid_id: sim_id, + expected: state.keys().into_iter().cloned().collect(), + }) } }) } @@ -288,7 +292,10 @@ pub extern "C" fn get_noise_model_by_name( as_capi_err(|| { let name = unsafe { CStr::from_ptr(name) } .to_str() - .map_err(|e| format!("UTF-8 error decoding representation argument: {}", e))?; + .map_err(|e| InvalidUtf8InArgument { + arg_name: "name".to_string(), + source: e, + })?; let noise_model = NoiseModel::get_by_name(name)?; let noise_model = CString::new(noise_model.as_json()).unwrap(); unsafe { @@ -320,20 +327,28 @@ pub extern "C" fn get_noise_model_by_name( #[no_mangle] pub extern "C" fn get_noise_model(sim_id: usize, noise_model_json: *mut *const c_char) -> i64 { as_capi_err(|| { - let state = &*STATE + let state = STATE .lock() - .map_err(|e| format!("Lock poisoning error: {}", e))?; + .map_err(|_| { + // Note that as per https://github.com/dtolnay/anyhow/issues/81#issuecomment-609247231, + // common practice is for poison errors to indicate that the containing thread + // has been irrevocably corrupted and must panic. + panic!("The lock on shared state for the C API has been poisoned."); + }) + .unwrap(); if let Some(sim_state) = state.get(&sim_id) { - let c_str = CString::new(sim_state.noise_model.as_json().as_str()).map_err(|e| { - format!("Null error while converting noise model to C string: {}", e) - })?; + let c_str = CString::new(sim_state.noise_model.as_json().as_str()) + .map_err(|e| UnanticipatedCApiError(anyhow::Error::new(e)))?; unsafe { *noise_model_json = c_str.into_raw(); - } + }; + Ok(()) } else { - return Err(format!("No simulator with id {} exists.", sim_id)); + Err(NoSuchSimulator { + invalid_id: sim_id, + expected: state.keys().into_iter().cloned().collect(), + }) } - Ok(()) }) } @@ -354,7 +369,7 @@ pub extern "C" fn get_noise_model(sim_id: usize, noise_model_json: *mut *const c pub unsafe extern "C" fn set_noise_model(sim_id: usize, new_model: *const c_char) -> i64 { as_capi_err(|| { if new_model.is_null() { - return Err("set_noise_model called with null pointer".to_string()); + return Err(NullPointer("new_model".to_string())); } let c_str = CStr::from_ptr(new_model); @@ -366,25 +381,18 @@ pub unsafe extern "C" fn set_noise_model(sim_id: usize, new_model: *const c_char sim_state.noise_model = noise_model; Ok(()) } else { - Err(format!("No simulator with id {} exists.", sim_id)) + Err(NoSuchSimulator { + invalid_id: sim_id, + expected: state.keys().into_iter().cloned().collect(), + }) } } - Err(serialization_error) => Err(format!( - "{} error deserializing noise model at line {}, column {}.", - match serialization_error.classify() { - serde_json::error::Category::Data => "Data / schema", - serde_json::error::Category::Eof => "End-of-file", - serde_json::error::Category::Io => "I/O", - serde_json::error::Category::Syntax => "Syntax", - }, - serialization_error.line(), - serialization_error.column() - )), + Err(err) => Err(JsonDeserializationError(err)), }, - Err(msg) => Err(format!( - "UTF-8 error decoding serialized noise model; was valid until byte {}.", - msg.valid_up_to() - )), + Err(msg) => Err(InvalidUtf8InArgument { + arg_name: "new_model".to_string(), + source: msg, + }), } }) } @@ -406,19 +414,25 @@ pub unsafe extern "C" fn set_noise_model(sim_id: usize, new_model: *const c_char pub unsafe extern "C" fn set_noise_model_by_name(sim_id: usize, name: *const c_char) -> i64 { as_capi_err(|| { if name.is_null() { - return Err("set_noise_model_by_name called with null pointer".to_string()); + return Err(NullPointer("name".to_string())); } let name = CStr::from_ptr(name) .to_str() - .map_err(|e| format!("UTF-8 error decoding name: {}", e))?; + .map_err(|e| InvalidUtf8InArgument { + arg_name: "name".to_string(), + source: e, + })?; let noise_model = NoiseModel::get_by_name(name)?; let state = &mut *STATE.lock().unwrap(); if let Some(sim_state) = state.get_mut(&sim_id) { sim_state.noise_model = noise_model; Ok(()) } else { - Err(format!("No simulator with id {} exists.", sim_id)) + Err(NoSuchSimulator { + invalid_id: sim_id, + expected: state.keys().into_iter().cloned().collect(), + }) } }) } diff --git a/src/Simulation/qdk_sim_rs/src/error.rs b/src/Simulation/qdk_sim_rs/src/error.rs new file mode 100644 index 00000000000..ec6ccf65ca4 --- /dev/null +++ b/src/Simulation/qdk_sim_rs/src/error.rs @@ -0,0 +1,116 @@ +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. + +//! Module defining common errors that can occur during quantum simulations. + +use std::str::Utf8Error; + +use miette::Diagnostic; +use thiserror::Error; + +/// Represents errors that can occur during linear algebra operations. +#[derive(Debug, Diagnostic, Error)] +pub enum QdkSimError { + // NB: As a design note, please consider if a more specific error is better + // suited for your usecase before returning `MiscError`. + /// Raised on miscellaneous errors. + #[error("{0}")] + #[diagnostic(code(qdk_sim::other))] + MiscError(String), + + /// Raised when functionality that has not yet been implemented is called. + #[error("Not yet implemented: {0}")] + #[diagnostic(code(qdk_sim::not_yet_implemented))] + NotYetImplemented(String), + + /// Raised when the wrong number of qubits is provided for a quantum + /// process. + #[error("Channel acts on {expected} qubits, but was applied to an {actual}-qubit state.")] + #[diagnostic(code(qdk_sim::process::wrong_n_qubits))] + WrongNumberOfQubits { + /// The number of qubits that was expected, as given by the size of the + /// channel to be applied. + expected: usize, + /// The actual number of qubits for the given state. + actual: usize, + }, + + /// Raised when a channel cannot be applied to a given state due to a + /// mismatch between channel and state kinds. + #[error("Unsupported quantum process variant {channel_variant} for applying to state variant {state_variant}.")] + #[diagnostic(code(qdk_sim::process::unsupported_apply))] + UnsupportedApply { + /// The enum variant of the channel to be applied. + channel_variant: &'static str, + /// The enum variant of the state that the channel is to be applied to. + state_variant: &'static str, + }, + + /// Raised when a matrix is singular, and thus does not have an inverse. + #[error("expected invertible matrix, but got a singular or very poorly conditioned matrix (det = {det})")] + #[diagnostic(code(qdk_sim::linalg::singular))] + Singular { + /// Actual determinant of the matrix which caused this error. + det: f64, + }, + + /// Raised when an algorithm requires a matrix to be square, but a + /// rectangular matrix was passed instead. + #[error("expected square matrix, but got shape `{0}` × `{1}")] + #[diagnostic(code(qdk_sim::linalg::not_square))] + NotSquare(usize, usize), + + /// Raised when an algorithm needs to convert an element between two + /// different scalar types, but no such conversion exists for those types. + #[error("could not convert value of type `{0}` into element type `{1}`")] + #[diagnostic(code(qdk_sim::linalg::cannot_convert_element))] + CannotConvertElement(String, String), + + /// Raised when no noise model exists for a given name. + #[error("{0} is not the name of any valid noise model")] + #[diagnostic(code(qdk_sim::noise_model::invalid_repr))] + InvalidNoiseModel(String), + + /// Raised when an initial state representation is invalid. + #[error("C API error: {0} is not a valid initial state representation")] + #[diagnostic(code(qdk_sim::c_api::invalid_repr))] + InvalidRepresentation(String), + + /// Raised when a null pointer is passed through the C API. + #[error("C API error: {0} was null")] + #[diagnostic(code(qdk_sim::c_api::nullptr))] + NullPointer(String), + + /// Raised when an invalid simulator ID is passed to the C API. + #[error("C API error: No simulator with ID {invalid_id} exists. Expected: {expected:?}.")] + #[diagnostic(code(qdk_sim::c_api::invalid_sim))] + NoSuchSimulator { + /// The invalid simulator id which caused this error. + invalid_id: usize, + /// A list of valid simulator ids at the point when this error occured. + expected: Vec, + }, + + /// Raised when a string passed to the C API contains could not be decoded + /// as a UTF-8 string. + #[error("C API error: UTF-8 error decoding {arg_name} argument: {source}")] + #[diagnostic(code(qdk_sim::c_api::utf8))] + InvalidUtf8InArgument { + /// The name of the argument containing invalid UTF-8 data. + arg_name: String, + + /// The underlying UTF-8 error that caused this error. + #[source] + source: Utf8Error, + }, + + /// Raised when a JSON serialization error occurs during a C API call. + #[error(transparent)] + #[diagnostic(code(qdk_sim::c_api::json_deser))] + JsonDeserializationError(#[from] serde_json::Error), + + /// Raised when an unanticipated error occurs during a C API call. + #[error(transparent)] + #[diagnostic(code(qdk_sim::c_api::unanticipated))] + UnanticipatedCApiError(#[from] anyhow::Error), +} diff --git a/src/Simulation/qdk_sim_rs/src/lib.rs b/src/Simulation/qdk_sim_rs/src/lib.rs index e451a567a3b..ae6e721dcc9 100644 --- a/src/Simulation/qdk_sim_rs/src/lib.rs +++ b/src/Simulation/qdk_sim_rs/src/lib.rs @@ -1,12 +1,11 @@ // Copyright (c) Microsoft Corporation. // Licensed under the MIT License. -// The following two attributes include the README.md for this crate when -// building docs (requires +nightly). -// See https://github.com/rust-lang/rust/issues/82768#issuecomment-803935643 -// for discussion. -#![cfg_attr(doc, feature(extended_key_value_attributes))] -#![cfg_attr(doc, cfg_attr(doc, doc = include_str!("../README.md")))] +#![cfg_attr(all(), doc = include_str!("../README.md"))] +#![cfg_attr( + feature = "document-features", + cfg_attr(doc, doc = ::document_features::document_features!()) +)] // Set linting rules for documentation. We will stop the build on missing docs, // or docs with known broken links. We only enable this when all relevant // features are enabled, otherwise the docs build will fail on links to @@ -17,7 +16,7 @@ // that are missing an `# Example` section. Currently, that raises a lot of // warnings when building docs, but ideally we should make sure to address // warnings going forward by adding relevant examples. -#![cfg_attr(doc, warn(missing_doc_code_examples))] +#![cfg_attr(doc, warn(rustdoc::missing_doc_code_examples))] #[macro_use(array, s)] extern crate ndarray; @@ -30,6 +29,7 @@ use std::usize; pub mod c_api; mod chp_decompositions; pub mod common_matrices; +pub mod error; mod instrument; pub mod linalg; mod noise_model; diff --git a/src/Simulation/qdk_sim_rs/src/linalg/array_ext.rs b/src/Simulation/qdk_sim_rs/src/linalg/array_ext.rs new file mode 100644 index 00000000000..cf2aa45fefb --- /dev/null +++ b/src/Simulation/qdk_sim_rs/src/linalg/array_ext.rs @@ -0,0 +1,129 @@ +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. + +//! Private extension traits and methods for working with arrays. + +use ndarray::{ + Array, ArrayBase, Axis, Data, DataMut, Dimension, Ix2, OwnedRepr, RawData, RemoveAxis, +}; +use num_traits::Zero; + +use crate::error::QdkSimError; + +pub(crate) trait ArrayBaseShapeExt { + fn require_square(&self) -> Result; +} + +impl ArrayBaseShapeExt for ArrayBase +where + S: RawData, + D: Dimension, +{ + fn require_square(&self) -> Result { + let dim = self.raw_dim(); + if dim[0] != dim[1] { + Err(QdkSimError::NotSquare(dim[0], dim[1])) + } else { + Ok(dim[0]) + } + } +} + +pub trait ArrayBaseMatrixExt +where + Self: Sized, +{ + type Output; + fn lower_triangular(&self) -> Self::Output; + fn upper_triangular(&self) -> Self::Output; +} +impl ArrayBaseMatrixExt for ArrayBase +where + S: Data, + A: Zero + Clone, +{ + type Output = ArrayBase, Ix2>; + + /// Returns a copy of the matrix with all elements above the diagonal set to zero. + fn lower_triangular(&self) -> Self::Output { + Array::from_shape_fn(self.raw_dim(), |(i, j)| { + if i >= j { + self[(i, j)].clone() + } else { + A::zero() + } + }) + } + + fn upper_triangular(&self) -> Self::Output { + Array::from_shape_fn(self.raw_dim(), |(i, j)| { + if i < j { + self[(i, j)].clone() + } else { + A::zero() + } + }) + } +} + +pub(crate) trait ArrayBaseRemoveAxisExt { + fn swap_index_axis(&mut self, axis: Axis, idx_source: usize, idx_dest: usize); +} +impl ArrayBaseRemoveAxisExt for ArrayBase +where + S: Data + DataMut, + ::Elem: Clone, + D: Dimension + RemoveAxis, +{ + fn swap_index_axis(&mut self, axis: Axis, idx_source: usize, idx_dest: usize) { + // TODO[perf]: Avoid copying both; this is needed currently to avoid mutably + // borrowing something which has an outstanding immutable + // borrow. + let source = self.index_axis(axis, idx_source).clone().to_owned(); + let dest = self.index_axis(axis, idx_dest).clone().to_owned(); + self.index_axis_mut(axis, idx_source).assign(&dest); + self.index_axis_mut(axis, idx_dest).assign(&source); + } +} + +#[cfg(test)] +mod tests { + use super::{ArrayBaseMatrixExt, ArrayBaseRemoveAxisExt}; + use ndarray::{array, Array2}; + + #[test] + fn upper_triangular_works() { + let mtx: Array2 = array![[6.0, 18.0, 3.0], [2.0, 12.0, 1.0], [4.0, 15.0, 3.0]]; + let upper = mtx.upper_triangular(); + assert_eq!( + upper, + array![[0.0, 18.0, 3.0], [0.0, 0.0, 1.0], [0.0, 0.0, 0.0]] + ); + } + + #[test] + fn lower_triangular_works() { + let mtx: Array2 = array![[6.0, 18.0, 3.0], [2.0, 12.0, 1.0], [4.0, 15.0, 3.0]]; + let lower = mtx.lower_triangular(); + assert_eq!( + lower, + array![[6.0, 0.0, 0.0], [2.0, 12.0, 0.0], [4.0, 15.0, 3.0]] + ); + } + + #[test] + fn swap_index_axis_swaps_rows_correctly() { + let mut mtx = array![[6.0, 18.0, 3.0], [2.0, 12.0, 1.0], [4.0, 15.0, 3.0]]; + mtx.swap_index_axis(ndarray::Axis(0), 1, 2); + let expected = array![[6.0, 18.0, 3.0], [4.0, 15.0, 3.0], [2.0, 12.0, 1.0]]; + assert_eq!(mtx, expected); + } + + #[test] + fn swap_index_axis_swaps_columns_correctly() { + let mut mtx = array![[6.0, 18.0, 3.0], [2.0, 12.0, 1.0], [4.0, 15.0, 3.0]]; + mtx.swap_index_axis(ndarray::Axis(1), 1, 2); + let expected = array![[6.0, 3.0, 18.0], [2.0, 1.0, 12.0], [4.0, 3.0, 15.0]]; + assert_eq!(mtx, expected); + } +} diff --git a/src/Simulation/qdk_sim_rs/src/linalg/decompositions/lu.rs b/src/Simulation/qdk_sim_rs/src/linalg/decompositions/lu.rs new file mode 100644 index 00000000000..95589e9419f --- /dev/null +++ b/src/Simulation/qdk_sim_rs/src/linalg/decompositions/lu.rs @@ -0,0 +1,293 @@ +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. + +use std::any::type_name; + +use cauchy::Scalar; +use ndarray::{Array1, ArrayBase, Axis, Data, Ix1, Ix2, OwnedRepr, RawData}; +use num_traits::FromPrimitive; + +use crate::{ + error::QdkSimError, + linalg::array_ext::{ArrayBaseRemoveAxisExt, ArrayBaseShapeExt}, +}; + +/// Represents the output of an LU decomposition acting on a matrix. +#[derive(Debug)] +pub struct LU +where + S: Data, + A: Scalar, +{ + pub(self) factors: ArrayBase, + pub(self) pivots: Array1, +} + +/// Represents types that represent matrices that are decomposable into lower- +/// and upper-triangular parts. +pub trait LUDecomposable { + /// The element type for the output matrices. + type Elem: Scalar; + + /// Type of owned data in the output matrices. + type OwnedRepr: Data + RawData; + + /// The output type for LU decompositions on this type. + type Output: LUDecomposition; + + /// The type used to represent errors in the decomposition. + type Error: std::error::Error; + + /// Performs an LU decomposition on the given type. + fn lu(&self) -> Result; +} + +impl LUDecomposable for ArrayBase +where + S: Data, + A: Scalar, +{ + type Elem = A; + type OwnedRepr = OwnedRepr; + type Error = QdkSimError; + type Output = LU>; + + // We follow the approach of + // https://github.com/mathnet/mathnet-numerics/blob/18c99a57562f181aab979b42ba206d3152f86cf2/src/Numerics/LinearAlgebra/Complex/Factorization/UserLU.cs. + fn lu(&self) -> Result { + let n_rows = self.require_square()?; + let threshold = A::Real::from_f64(1e-10).ok_or_else(|| { + QdkSimError::CannotConvertElement("f64".to_string(), type_name::().to_string()) + })?; + + let factors = &mut (*self).to_owned(); + let mut pivots: Array1<_> = (0..n_rows).collect::>().into(); + + for j in 0..n_rows { + let mut pivot_col = factors.index_axis(Axis(1), j).to_owned(); + + // We do this in a for loop since the slicing to get the lower-triangular + // part comes out harder to read. + for i in 0..n_rows { + let mut s = A::zero(); + for k in 0..std::cmp::min(i, j) { + s += factors[(i, k)] * pivot_col[k]; + } + pivot_col[i] -= s; + factors[(i, j)] = pivot_col[i]; + } + + // Find the pivot and exchange if needed. + let mut idx_pivot = j; + for i in j + 1..n_rows { + if pivot_col[i].abs() > pivot_col[idx_pivot].abs() { + idx_pivot = i; + } + } + + if idx_pivot != j { + factors.swap_index_axis(Axis(0), idx_pivot, j); + + // Record the exchange in the pivot permutation. + pivots[j] = idx_pivot; + } + + // Compute multipliers. + if j < n_rows && factors[(j, j)].abs() >= threshold { + for i in j + 1..n_rows { + let elem = factors[(j, j)]; + factors[(i, j)] /= elem; + } + } + } + + Ok(LU { + factors: factors.clone(), + pivots, + }) + } +} + +/// Represents the result of decomposing a square matrix $A$ into lower- and +/// upper-triangular components $L$ and $U$, respectively. +pub trait LUDecomposition +where + S: Data, + A: Scalar, +{ + /// The type used to represent errors in using this decomposition to solve + /// matrix equations. + type Error: std::error::Error; + + /// The type resulting from using this decomposition to solve problems + /// of the form $A\vec{x} = \vec{y}$. + type VectorSolution; + + /// The type resulting from using this decomposition to solve problems + /// of the form $AX = Y$. + type MatrixSolution; + + /// The size $n$ of the matrix whose decomposition is represented. In + /// particular, $L\Pi U$ is taken to be an $n \times n$ matrix, where $L$ + /// and $U$ are the lower- and upper-triangular factors, and where $\Pi$ + /// is a permutation matrix. + fn order(&self) -> usize; + + // TODO: Change signature to be in-place. + /// Uses this decomposition to solve an equation of the form + /// $A\vec{x} = \vec{y}$ for $\vec{x}$ given $\vec{y}$, where $A$ is + /// implicitly represented by this decomposition. + fn solve_vector(&self, rhs: &ArrayBase) -> Result; + + // TODO: Change signature to be in-place. + /// Uses this decomposition to solve an equation of the form + /// $AX = Y$ for $X$ given $Y$, where $A$ is + /// implicitly represented by this decomposition. + fn solve_matrix(&self, rhs: &ArrayBase) -> Result; +} + +impl LUDecomposition for LU +where + S: Data, + A: Scalar, +{ + type Error = QdkSimError; + type VectorSolution = ArrayBase, Ix1>; + type MatrixSolution = ArrayBase, Ix2>; + + fn solve_vector(&self, rhs: &ArrayBase) -> Result { + todo!() + } + + fn solve_matrix(&self, rhs: &ArrayBase) -> Result { + let n_rows = rhs.require_square()?; + let mut result = rhs.to_owned(); + for (idx_col, idx_pivot) in self.pivots.iter().enumerate() { + if idx_col != *idx_pivot { + result.swap_index_axis(Axis(0), idx_col, *idx_pivot); + } + } + + // NB: We use explicit for loops here as using vectorized indexing + // for triangular slices is somewhat challenging to do robustly. + // + // Note that by convention, we have stored both the lower- and + // upper-triangular parts of the LU decomposition in the same + // array, so restricting how we index the decomposition array gives + // us each part of the decomposition. + + // Start by solving by using the lower-triangular part. + for k in 0..n_rows { + for i in k + 1..n_rows { + for j in 0..n_rows { + let temp = result[(k, j)] * self.factors[(i, k)]; + result[(i, j)] -= temp; + } + } + } + + // Solve the upper-triangular part. + for k in (0..n_rows).rev() { + for j in 0..n_rows { + result[(k, j)] /= self.factors[(k, k)]; + } + for i in 0..k { + for j in 0..n_rows { + let temp = result[(k, j)] * self.factors[(i, k)]; + result[(i, j)] -= temp; + } + } + } + + Ok(result) + } + + fn order(&self) -> usize { + self.factors.shape()[0] + } +} + +#[cfg(test)] +mod tests { + use approx::assert_abs_diff_eq; + use cauchy::c64; + use ndarray::{array, Array2, OwnedRepr}; + + use crate::{ + c64, + error::QdkSimError, + linalg::decompositions::{LUDecomposable, LU}, + }; + + #[test] + fn lu_decomposition_works_f64() -> Result<(), QdkSimError> { + let mtx = array![[6.0, 18.0, 3.0], [2.0, 12.0, 1.0], [4.0, 15.0, 3.0]]; + let lu: LU> = mtx.lu()?; + + // NB: This tests the internal structure of the LU decomposition, and + // may validly fail if the algorithm above is modified. + let expected_factors = array![ + [6.0, 18.0, 3.0], + [0.3333333333333333, 6.0, 0.0], + [0.6666666666666666, 0.5, 1.0], + ]; + for (actual, expected) in lu.factors.iter().zip(expected_factors.iter()) { + assert_abs_diff_eq!(actual, expected, epsilon = 1e-6); + } + + let expected_pivots = vec![0, 1, 2]; + assert_eq!(lu.pivots.to_vec(), expected_pivots); + Ok(()) + } + + #[test] + fn lu_decomposition_works_c64() -> Result<(), QdkSimError> { + // In [1]: import scipy.linalg as la + // In [2]: la.lu([ + // ...: [-1, 1j, -2], + // ...: [3, 0, -4j], + // ...: [-1, 5, -1] + // ...: ]) + // Out[2]: (array([[0., 0., 1.], + // [1., 0., 0.], + // [0., 1., 0.]]), + // array([[ 1. +0.j , 0. +0.j , 0. +0.j ], + // [-0.33333333+0.j , 1. +0.j , 0. +0.j ], + // [-0.33333333+0.j , 0. +0.2j, 1. +0.j ]]), + // array([[ 3. +0.j , 0. +0.j , + // -0. -4.j ], + // [ 0. +0.j , 5. +0.j , + // -1. -1.33333333j], + // [ 0. +0.j , 0. +0.j , + // -2.26666667-1.13333333j]])) + let mtx: Array2 = array![ + [c64!(-1.0), c64!(1.0 i), c64!(-2.0)], + [c64!(3.0), c64!(0.0), c64!(-4.0 i)], + [c64!(-1.0), c64!(5.0), c64!(-1.0)] + ]; + let lu: LU> = mtx.lu()?; + + // NB: This tests the internal structure of the LU decomposition, and + // may validly fail if the algorithm above is modified. + let expected_factors = array![ + [c64::new(3.0, 0.0), c64::new(0.0, 0.0), c64::new(0.0, -4.0)], + [ + c64::new(-0.3333333333333333, 0.0), + c64::new(5.0, 0.0), + c64::new(-1.0, -1.3333333333333333) + ], + [ + c64::new(-0.3333333333333333, 0.0), + c64::new(0.0, 0.2), + c64::new(-2.26666667, -1.13333333) + ], + ]; + for (actual, expected) in lu.factors.iter().zip(expected_factors.iter()) { + assert_abs_diff_eq!(actual, expected, epsilon = 1e-6); + } + + let expected_pivots = vec![1, 2, 2]; + assert_eq!(lu.pivots.to_vec(), expected_pivots); + Ok(()) + } +} diff --git a/src/Simulation/qdk_sim_rs/src/linalg/decompositions/mod.rs b/src/Simulation/qdk_sim_rs/src/linalg/decompositions/mod.rs new file mode 100644 index 00000000000..6b95aa272ec --- /dev/null +++ b/src/Simulation/qdk_sim_rs/src/linalg/decompositions/mod.rs @@ -0,0 +1,4 @@ +//! Support for various decompositions used in linear algebra. + +mod lu; +pub use lu::*; diff --git a/src/Simulation/qdk_sim_rs/src/linalg/inv.rs b/src/Simulation/qdk_sim_rs/src/linalg/inv.rs new file mode 100644 index 00000000000..d6c89a4beb4 --- /dev/null +++ b/src/Simulation/qdk_sim_rs/src/linalg/inv.rs @@ -0,0 +1,784 @@ +use cauchy::Scalar; +use ndarray::{Array1, Array2, OwnedRepr}; + +use crate::linalg::decompositions::{LUDecomposable, LUDecomposition}; + +/// Types that support the matrix inverse $A^{-1}$. +pub trait Inv { + /// Errors that can result from the [`inv`] method. + /// + /// [`inv`]: Inv::inv + type Error; + /// Type used to represent $A^{-1}$ when calling [`inv`]. + /// + /// [`inv`]: Inv::inv + type Output; + + /// Computes the inverse of a given matrix. + fn inv(&self) -> Result; +} + +impl Inv for M +where + M: LUDecomposable, Error = E>, + A: Scalar, + // TODO: Allow decomposables and decompositions to have different errors. + ::Output: LUDecomposition, Error = E>, +{ + type Error = M::Error; + type Output = + <::Output as LUDecomposition>>::MatrixSolution; + + fn inv(&self) -> Result { + // TODO: Check determinant here. + let lu = self.lu()?; + // TODO: Change to use Identity trait. + let id_mtx = Array2::from_diag(&Array1::ones(lu.order())); + lu.solve_matrix(&id_mtx) + } +} + +#[cfg(test)] +mod tests { + use approx::assert_abs_diff_eq; + use cauchy::c64; + use ndarray::array; + use num_traits::Zero; + + use crate::{c64, error::QdkSimError, linalg::Inv}; + + #[test] + fn inv_works_f64() -> Result<(), QdkSimError> { + let mtx = array![[6.0, 18.0, 3.0], [2.0, 12.0, 1.0], [4.0, 15.0, 3.0]]; + // TODO: Actually write the test! + let inv = mtx.inv()?; + println!("{:?}", inv); + Ok(()) + } + + #[test] + fn inv_works_c64() -> Result<(), QdkSimError> { + let mtx = array![ + [c64::zero(), c64::new(0.0, 1.0), c64::new(2.0, 0.0)], + [c64::new(0.0, 3.0), c64::new(4.0, 0.0), c64::new(0.0, 5.0)], + [c64::new(6.0, 0.0), c64::new(0.0, 7.0), c64::new(8.0, 0.0)] + ]; + let inv = mtx.inv()?; + let expected_inv = array![ + [c64!(-0.69791667), c64!(-0.0625 i), c64!(0.13541667),], + [c64!(-0.0625 i), c64!(0.125), c64!(-0.0625 i)], + [c64!(0.46875), c64!(-0.0625 i), c64!(-0.03125)], + ]; + + for (actual, expected) in inv.iter().zip(expected_inv.iter()) { + assert_abs_diff_eq!(actual, expected, epsilon = 1e-6); + } + Ok(()) + } + + #[test] + fn inv_works_c64_5x5() -> Result<(), QdkSimError> { + // Python code to generate test case: + // import numpy as np + // import scipy.linalg as la + // X = np.random.random((5, 5)) + np.random.random((5, 5)) * 1j + // def print_arr_as_rs(arr): + // print("array![") + // for row in arr: + // print(" [", end="") + // print(", ".join( + // f"c64::new({element.real}, {element.imag})" + // for element in row + // ), end="") + // print("],") + // print("]") + // print_arr_as_rs(X) + // print_arr_as_rs(la.inv(X)) + let mtx = array![ + [ + c64::new(0.25759564633928533, 0.9499197824003289), + c64::new(0.36288815045785816, 0.7542794140806953), + c64::new(0.3575366182927219, 0.15718447037949934), + c64::new(0.6258618223229145, 0.606098740015765), + c64::new(0.4262277886124155, 0.952615526372946) + ], + [ + c64::new(0.7574275030894558, 0.5980955544457704), + c64::new(0.8519089404515027, 0.37274304468840325), + c64::new(0.6484719497482538, 0.29885541821603523), + c64::new(0.047156400699850165, 0.7003281373877304), + c64::new(0.7573593614243241, 0.02022619925054081) + ], + [ + c64::new(0.06565657022052329, 0.8964247633381773), + c64::new(0.43971180912111907, 0.0936752171206513), + c64::new(0.8565111965854413, 0.716416438691636), + c64::new(0.4620853548253515, 0.08859493859607226), + c64::new(0.6761502941167822, 0.7629177012545356) + ], + [ + c64::new(0.28592498922650456, 0.11505276541596976), + c64::new(0.8638866172586472, 0.3860930727399199), + c64::new(0.06899006463605606, 0.8630942739465155), + c64::new(0.07967078190474364, 0.6271701884478902), + c64::new(0.9889623172686699, 0.6831913961093329) + ], + [ + c64::new(0.9198323331490859, 0.9418968234232035), + c64::new(0.7956810140209064, 0.7897894731423162), + c64::new(0.9390139934709593, 0.9133215782294428), + c64::new(0.4983734769859526, 0.0994474246120769), + c64::new(0.2619209371310002, 0.004972757766929403) + ], + ]; + let inv = mtx.inv()?; + let expected_inv = array![ + [ + c64::new(-0.7909619701222462, 0.4675463791528054), + c64::new(-0.09950302744992989, -1.5959803664920833), + c64::new(0.681218090153333, -1.0192017365612318), + c64::new(0.4942583087733383, 1.4655410303969532), + c64::new(0.3918368731842419, 0.555047115830301) + ], + [ + c64::new(0.14432483967352205, -1.2656338362747568), + c64::new(1.3951911164644337, 1.1214902968878036), + c64::new(-0.21697568748595358, 1.7970639543704376), + c64::new(-1.0823513631779629, -0.5997321814949002), + c64::new(-0.4824879562743729, -0.7262272675611384) + ], + [ + c64::new(0.4217193550149313, 0.8408538401399359), + c64::new(-0.7486966253548832, 0.4967589619972676), + c64::new(-0.31853978279704065, -0.535332891638263), + c64::new(0.4293214150441409, -0.680836516539264), + c64::new(0.4312617193077752, -0.528250661249761) + ], + [ + c64::new(0.7120409190490832, -0.057771642957427746), + c64::new(-0.8431484241081729, -0.4785306364974677), + c64::new(-0.23128119608361064, 0.041150277748116), + c64::new(0.17415737286549077, -0.3041063301615835), + c64::new(0.30229406523570523, 0.20878202324721712) + ], + [ + c64::new(0.24696145093874475, -0.038144733375202966), + c64::new(-0.3996802757158412, 0.41676978509463203), + c64::new(-0.26054584996722907, -0.6309725521979278), + c64::new(0.6465176542968839, -0.39453347662338617), + c64::new(-0.16680206233005246, 0.18386759857791737) + ], + ]; + + for (actual, expected) in inv.iter().zip(expected_inv.iter()) { + assert_abs_diff_eq!(actual, expected, epsilon = 1e-6); + } + Ok(()) + } + + #[test] + fn inv_works_c64_16x16() -> Result<(), QdkSimError> { + // Python code to generate test case: + // import numpy as np + // import scipy.linalg as la + // X = np.random.random((16, 16)) + np.random.random((16, 16)) * 1j + // def print_arr_as_rs(arr): + // print("array![") + // for row in arr: + // print(" [", end="") + // print(", ".join( + // f"c64::new({element.real}, {element.imag})" + // for element in row + // ), end="") + // print("],") + // print("]") + // print_arr_as_rs(X) + // print_arr_as_rs(la.inv(X)) + let mtx = array![ + [ + c64::new(0.24343550539460002, 0.7096397032845722), + c64::new(0.629546686432831, 0.5301006031714603), + c64::new(0.4199349954349497, 0.358375134500281), + c64::new(0.400595688139045, 0.9645759115525548), + c64::new(0.08639871723999515, 0.7655041461418464), + c64::new(0.5272421501371002, 0.49548242232242845), + c64::new(0.201862871149052, 0.691752157021674), + c64::new(0.566606570806358, 0.9090447061650156), + c64::new(0.5224821650372808, 0.5254829260378098), + c64::new(0.5264656718847738, 0.37651456343277556), + c64::new(0.5561885182402057, 0.9143124378366423), + c64::new(0.49710125441166253, 0.030484298760797723), + c64::new(0.1631863300384664, 0.8113573847147066), + c64::new(0.8400531600269201, 0.43703798604805766), + c64::new(0.1730150235896153, 0.4988400829932095), + c64::new(0.07911501549474143, 0.4510850896893345) + ], + [ + c64::new(0.4182284703927561, 0.17768533289509136), + c64::new(0.5436803677118057, 0.28670911836907875), + c64::new(0.8660479956675148, 0.6718934426880514), + c64::new(0.8224549142195182, 0.09964377556322512), + c64::new(0.34496329518768976, 0.9801731197806263), + c64::new(0.5849863010766326, 0.4490005338446601), + c64::new(0.07836771055042646, 0.5420110965117988), + c64::new(0.4669545625108412, 0.1407565997515524), + c64::new(0.42484012971040686, 0.8759094071016321), + c64::new(0.5205341023964746, 0.9923357553551803), + c64::new(0.7014879416155709, 0.37485464528319457), + c64::new(0.5471165052558605, 0.20796849150802799), + c64::new(0.6199629816663412, 0.7519415024152558), + c64::new(0.9037790360470818, 0.666480231673042), + c64::new(0.2086744777910019, 0.07827526780188987), + c64::new(0.4670844989586508, 0.9537412452187163) + ], + [ + c64::new(0.4140469741259334, 0.4560336908973007), + c64::new(0.2546277825179136, 0.08143353851401769), + c64::new(0.6537603463090524, 0.6570297683362842), + c64::new(0.4009562738772193, 0.8805488266979159), + c64::new(0.02207598382742204, 0.5920605348013798), + c64::new(0.05517423266462318, 0.793992529375867), + c64::new(0.5632193027908459, 0.5977005530132938), + c64::new(0.9017864875558482, 0.5339011375670782), + c64::new(0.3828142959921751, 0.339242061558683), + c64::new(0.3385209625213713, 0.8324049401182788), + c64::new(0.30970070432395513, 0.005197774068218308), + c64::new(0.20695380217097692, 0.472246800325997), + c64::new(0.20221796644122614, 0.6857628071658661), + c64::new(0.17174699519185188, 0.6929524277609821), + c64::new(0.8041808362434856, 0.3371802625347794), + c64::new(0.9362896880145364, 0.38138437444507334) + ], + [ + c64::new(0.4696066202823814, 0.03198447459798637), + c64::new(0.6323638366681098, 0.9258645924825338), + c64::new(0.48515708183133144, 0.11024424103997121), + c64::new(0.609317530558061, 0.2972302001516389), + c64::new(0.30060843150464833, 0.8100153459051707), + c64::new(0.7986668256449622, 0.3922734957863946), + c64::new(0.2729653457634389, 0.4953984142309359), + c64::new(0.06098942703216337, 0.2004705039464001), + c64::new(0.4106442591731715, 0.47107371888912153), + c64::new(0.8932416419338051, 0.35348623017807734), + c64::new(0.9437579522965495, 0.3028603115275772), + c64::new(0.8997707657703546, 0.7534540438416928), + c64::new(0.8732369875427841, 0.6374743043364531), + c64::new(0.381827893689629, 0.3927530492921386), + c64::new(0.16011119980494248, 0.1146828362583655), + c64::new(0.4209728872731916, 0.9462487108084525) + ], + [ + c64::new(0.9179230088359127, 0.5038538240712301), + c64::new(0.2889594017105449, 0.787940771523522), + c64::new(0.4598053875867901, 0.41623746452910215), + c64::new(0.8095036581534074, 0.21551533218788865), + c64::new(0.4268228153778082, 0.9889060962179577), + c64::new(0.4304494380785887, 0.8040452167212958), + c64::new(0.223311800504025, 0.4970151956113278), + c64::new(0.007997571773524781, 0.6317969612629845), + c64::new(0.7469084890757939, 0.3363308133241615), + c64::new(0.7207259005808849, 0.6053258370249136), + c64::new(0.22517978068201328, 0.76379017292531), + c64::new(0.012292832735091963, 0.35381918017361436), + c64::new(0.8538576693766322, 0.464070671170861), + c64::new(0.429655967580578, 0.19559316787868042), + c64::new(0.935606881298734, 0.7033910463384877), + c64::new(0.4549095375495651, 0.2058092116780722) + ], + [ + c64::new(0.3865520894048742, 0.002013314615649353), + c64::new(0.22015446091828872, 0.38043374130648455), + c64::new(0.5048931185357252, 0.4194081288354866), + c64::new(0.6158067241441327, 0.7767694835770059), + c64::new(0.7584962310602148, 0.9546886708614651), + c64::new(0.5822353268156216, 0.537889230679487), + c64::new(0.8662438490163971, 0.7545107340070697), + c64::new(0.9083377609146424, 0.6927881140843738), + c64::new(0.9898651658147615, 0.11855995178629275), + c64::new(0.5420545051674228, 0.7449456850548131), + c64::new(0.6300620336271961, 0.3901867594702042), + c64::new(0.34711633901876726, 0.018969229982951585), + c64::new(0.2583485690095335, 0.3627546968287746), + c64::new(0.6679984633215221, 0.8079392177620631), + c64::new(0.456269593518148, 0.672929467478175), + c64::new(0.7703161905109068, 0.5312644373073533) + ], + [ + c64::new(0.19545257016136797, 0.7486267918001699), + c64::new(0.45956607922957327, 0.02495271845680236), + c64::new(0.23801167172055393, 0.5234015133852364), + c64::new(0.8151487989582498, 0.9925066835945822), + c64::new(0.4173531176676154, 0.9983961026357865), + c64::new(0.08850701855742737, 0.6651580674279115), + c64::new(0.14664093713712056, 0.9616446116985102), + c64::new(0.06556846579674958, 0.4332711601853556), + c64::new(0.00999669313077256, 0.9503289570965522), + c64::new(0.8492820744902253, 0.19441515542717458), + c64::new(0.1196791277954703, 0.3135973114594467), + c64::new(0.8974320738101824, 0.5088961076620361), + c64::new(0.3077492417211729, 0.0587507288460859), + c64::new(0.6718542468314663, 0.8784564572108206), + c64::new(0.8058709416507477, 0.5693768630795327), + c64::new(0.5608670604971727, 0.5311021334681032) + ], + [ + c64::new(0.34560208388480096, 0.16739020607281618), + c64::new(0.23331826916254472, 0.6130815663685395), + c64::new(0.9397427753841658, 0.6655565440854049), + c64::new(0.8821156688165074, 0.7944792018803679), + c64::new(0.42739938631087837, 0.2878194026469886), + c64::new(0.9233635811530989, 0.09216121753894313), + c64::new(0.885269983678575, 0.10561646132268732), + c64::new(0.7306049360915047, 0.28794390992481966), + c64::new(0.7180647224566964, 0.38575181855530294), + c64::new(0.6028287170510034, 0.1452395558557621), + c64::new(0.5858041396491549, 0.3765295350341745), + c64::new(0.5234838270893098, 0.9291711699124525), + c64::new(0.7220257096923388, 0.681906888190176), + c64::new(0.44671808608453434, 0.25263171502629855), + c64::new(0.15387870406855686, 0.01859599479014662), + c64::new(0.7108714026423096, 0.9921637786908944) + ], + [ + c64::new(0.5522540902713062, 0.5289418043806777), + c64::new(0.17955033836453782, 0.6641512242937374), + c64::new(0.2870982224237145, 0.8253475828830054), + c64::new(0.6953135064981454, 0.5389399858548773), + c64::new(0.9044154166667651, 0.7291375690536017), + c64::new(0.014969049881469187, 0.5823645756346624), + c64::new(0.4106110274919824, 0.9772441717239208), + c64::new(0.7899475513249096, 0.33718757706048264), + c64::new(0.12564982302572625, 0.08164867163164302), + c64::new(0.8277448644423485, 0.6684800005333317), + c64::new(0.8063722001578879, 0.293139957786377), + c64::new(0.7616443403756338, 0.41248468565508944), + c64::new(0.48370112983282043, 0.5237165655219297), + c64::new(0.9852738590227202, 0.7814762792425599), + c64::new(0.3150523000233134, 0.25232529923596003), + c64::new(0.7055164700072192, 0.4364210816486238) + ], + [ + c64::new(0.47480263870579453, 0.5668225452197742), + c64::new(0.8861112169309894, 0.7839593615514723), + c64::new(0.47984998067122964, 0.6940783683811662), + c64::new(0.6365911151977862, 0.8508530098815925), + c64::new(0.9221748970317584, 0.37814148012826587), + c64::new(0.41247888626124285, 0.6036385621285446), + c64::new(0.2077081234948941, 0.6237160782342986), + c64::new(0.3014578023119979, 0.4601857335517109), + c64::new(0.2594737250836574, 0.6075684486769833), + c64::new(0.8128732699081882, 0.41023156765679436), + c64::new(0.36232588883648964, 0.7338270593218202), + c64::new(0.5175615608791986, 0.04730733357986405), + c64::new(0.16569218047093304, 0.20214314530371735), + c64::new(0.4906398311215432, 0.6566575126029983), + c64::new(0.6122744130095827, 0.5855960128693636), + c64::new(0.8069280585173866, 0.6446311686762304) + ], + [ + c64::new(0.8873728806853252, 0.5005611313637354), + c64::new(0.9770223512150816, 0.6343132757665193), + c64::new(0.4061148961549872, 0.6729021869161408), + c64::new(0.457189125698762, 0.06439923211784304), + c64::new(0.8080764133185422, 0.8427563110745654), + c64::new(0.14742320294056432, 0.7461451243453329), + c64::new(0.6621587725591924, 0.07592063914586178), + c64::new(0.9119688176628091, 0.5791055651876627), + c64::new(0.05277334761942254, 0.7149032903155043), + c64::new(0.9617330181619003, 0.3459175480617631), + c64::new(0.8273379741828939, 0.10046702842390531), + c64::new(0.0053714122905957895, 0.46193759316661875), + c64::new(0.38494207909599165, 0.9376983173170867), + c64::new(0.7355626689153292, 0.744918797249697), + c64::new(0.3565603523759998, 0.6328091577590916), + c64::new(0.48732908110439943, 0.9961922465946875) + ], + [ + c64::new(0.8535294237098889, 0.4555765978417987), + c64::new(0.7986881125668084, 0.51780181927999), + c64::new(0.7401082539212396, 0.4260186037296265), + c64::new(0.5005532719840523, 0.8663437529938097), + c64::new(0.7666829073963346, 0.053633375518612136), + c64::new(0.09618596396671064, 0.14132877650296938), + c64::new(0.41325974785966146, 0.5188530456833758), + c64::new(0.4951335463859594, 0.7681452338571635), + c64::new(0.13337087025353667, 0.4534550761637307), + c64::new(0.264802459443372, 0.6960503802635601), + c64::new(0.10620477984934751, 0.16952615428155038), + c64::new(0.41326087740330664, 0.11093454026305205), + c64::new(0.11885855965831427, 0.8643900810326837), + c64::new(0.9531576858961567, 0.6405661730490178), + c64::new(0.2256965951584058, 0.6929381202000422), + c64::new(0.8775026501055405, 0.06819532335632383) + ], + [ + c64::new(0.05366306018663092, 0.6786935353784738), + c64::new(0.7049618674230367, 0.5343466515014068), + c64::new(0.1971345668641724, 0.9945961106878076), + c64::new(0.21850330503480253, 0.9315673412181285), + c64::new(0.46085694564476165, 0.41750086472796577), + c64::new(0.7797969253742099, 0.37698595587188244), + c64::new(0.22352377642177046, 0.6202647253113842), + c64::new(0.3947773505412384, 0.30392363084842966), + c64::new(0.79796785522656, 0.3573353874612716), + c64::new(0.23931143583970305, 0.712883482793276), + c64::new(0.14995225507944243, 0.6591874658307774), + c64::new(0.15535778984760662, 0.6646611213320233), + c64::new(0.7203947478912175, 0.2850133023325271), + c64::new(0.6466176117241828, 0.7808674648656392), + c64::new(0.5041129975238132, 0.36145725474426016), + c64::new(0.5302947663343138, 0.7083319689929208) + ], + [ + c64::new(0.7320943239049572, 0.908797021842243), + c64::new(0.7427211309634909, 0.3356960467687501), + c64::new(0.22286966092460725, 0.27543792599706995), + c64::new(0.05098355106161934, 0.12880818021612772), + c64::new(0.42913896501867776, 0.48225626640930075), + c64::new(0.7247114618619713, 0.6201582971419118), + c64::new(0.2990209549883077, 0.6065353578264744), + c64::new(0.8803034486159311, 0.46632710634376084), + c64::new(0.3900337102002004, 0.8603415217544679), + c64::new(0.03161503771138252, 0.2020000054053933), + c64::new(0.03285549876078819, 0.7197691844048204), + c64::new(0.8861914311445265, 0.4524597317254516), + c64::new(0.919139071989545, 0.9903753036352839), + c64::new(0.26413815520887607, 0.6125517285372406), + c64::new(0.9557958183246573, 0.5678539269316171), + c64::new(0.40410119324680294, 0.27980193142477106) + ], + [ + c64::new(0.27285675504845885, 0.14225223021480682), + c64::new(0.15871896780671346, 0.30806072726128164), + c64::new(0.28647798096884347, 0.37659944477098006), + c64::new(0.9539651068234646, 0.42730784410934874), + c64::new(0.9918958737388014, 0.5885319671531076), + c64::new(0.13253799282053358, 0.6799795473579588), + c64::new(0.5697443682755279, 0.9692621115480067), + c64::new(0.9408520393874307, 0.635189510953812), + c64::new(0.6378103987865228, 0.94962613572105), + c64::new(0.8980087424095724, 0.2840861303932122), + c64::new(0.5177631962541133, 0.5298970606217572), + c64::new(0.9173844985722672, 0.5722269153975366), + c64::new(0.10624346021142306, 0.6437026744926997), + c64::new(0.8234387192061203, 0.6901214486979598), + c64::new(0.6394729661348116, 0.6883637984832083), + c64::new(0.7781863797181444, 0.9197816550850407) + ], + [ + c64::new(0.8221647342763111, 0.3724858860937563), + c64::new(0.8384533817923859, 0.12061755516837247), + c64::new(0.14283137274059887, 0.2633143962286675), + c64::new(0.8479111175389686, 0.4143606483497789), + c64::new(0.8586200226659483, 0.04144888283592396), + c64::new(0.5924050968149841, 0.2936642057302724), + c64::new(0.26548598304261006, 0.5566420247616555), + c64::new(0.8270951953157784, 0.43319207911455204), + c64::new(0.6449007431399179, 0.16135807058447116), + c64::new(0.8489376152853729, 0.2423209662369905), + c64::new(0.9080821967375533, 0.1456211521082592), + c64::new(0.6165589096032493, 0.33794567982637114), + c64::new(0.39510580490087865, 0.4541868882549819), + c64::new(0.6767219080169666, 0.015218396907329734), + c64::new(0.5013860268093684, 0.21134538382323254), + c64::new(0.8745499679879252, 0.553499342738368) + ], + ]; + let inv = mtx.inv()?; + let expected_inv = array![ + [ + c64::new(-0.33341646484875576, 1.1052853101713702), + c64::new(-0.9524670156810291, -1.4260083616552919), + c64::new(-0.28008723834080096, 2.370239467388982), + c64::new(1.9975641630467906, 2.807918119574881), + c64::new(-2.145372013203242, 0.372298143790908), + c64::new(1.9668787833875745, -1.4178060687148868), + c64::new(1.628648632004574, -3.7986548511472975), + c64::new(-0.3841234618774047, -1.6266964977599572), + c64::new(0.5422313413841024, 1.2455849319689138), + c64::new(0.903361443678693, 1.0898735025274768), + c64::new(1.1133856248111635, -0.4701055609971547), + c64::new(-0.24962622846943747, -1.7771615655750181), + c64::new(0.5505173511853978, 1.158698308900222), + c64::new(0.5203113052796122, 0.2383987562681391), + c64::new(-3.692746889492212, 2.0725054682902377), + c64::new(0.3038728322291176, -2.214101409455073) + ], + [ + c64::new(0.6983382732954325, -1.5158076344322589), + c64::new(-0.20357552550341684, 1.5929538750963994), + c64::new(1.2274983134845052, -2.3218771130111637), + c64::new(0.10693506334445249, -2.6282640973556752), + c64::new(1.4343432776702896, 0.24146199302301508), + c64::new(-2.4976548307366, 1.3294740940569563), + c64::new(-2.802338558934232, 2.1485292671458254), + c64::new(-1.2647768773084678, 1.4373996801440505), + c64::new(-0.7052630223190702, -0.9023615614347303), + c64::new(-0.09089496809290643, -0.6402466613665145), + c64::new(-0.40136582493590267, 0.062206713362174204), + c64::new(0.22716356730332987, 1.4219236621098785), + c64::new(0.9769380171392674, -1.7439870854541086), + c64::new(-0.6518850947761342, -0.09063262579860987), + c64::new(3.8490253998520467, -0.6602384272549691), + c64::new(-1.266366905888641, 1.2948650691721733) + ], + [ + c64::new(-2.1150334917169067, -0.7524716105762088), + c64::new(2.865705855106625, -1.7919232385720578), + c64::new(-4.6372368665191495, 0.17211604185024365), + c64::new(-4.2548274092020835, 4.405149160028291), + c64::new(-1.273424038031718, -4.607022884361642), + c64::new(3.538122864770967, 3.7768700142222427), + c64::new(7.1939705238336025, 2.432924405835742), + c64::new(2.9121489577395616, -0.6869034106342744), + c64::new(-2.2868064317950814, 1.7817946421397064), + c64::new(-1.1327618966999562, 2.1697296568784603), + c64::new(0.8446373052778617, 1.2833732909034243), + c64::new(3.9395250724162416, -1.3807357928110635), + c64::new(-2.3853552783633574, 1.0480719154218856), + c64::new(-0.8893681746460103, 1.466190836254347), + c64::new(-4.414495537980002, -6.84540935784314), + c64::new(2.6889373223511144, -0.2690773578058888) + ], + [ + c64::new(-0.9411498074258554, -0.21883261987691854), + c64::new(0.734218533330576, -0.25138850289824616), + c64::new(-1.413778097867341, 0.13946356232182566), + c64::new(-0.9230913277697869, 1.1939182811997437), + c64::new(-0.500623877942419, -0.9632715774539262), + c64::new(1.2147048376728817, 0.7973868770253624), + c64::new(2.124908888672486, 0.0620633188311287), + c64::new(0.8530892013630627, -0.47480774043129137), + c64::new(-0.41954658192217403, 0.6143444446321021), + c64::new(-0.16212043344560645, 0.455921823842463), + c64::new(0.1452607123772317, 0.5140894248791609), + c64::new(0.879848490746153, -0.835664894813614), + c64::new(-0.40773386118785726, 0.595894534716184), + c64::new(-0.3107685669271478, 0.632628292161437), + c64::new(-1.6066598484534516, -1.4202036674117309), + c64::new(1.0679617005904698, -0.21927353401787894) + ], + [ + c64::new(-1.2006792624267157, -3.0785366301386503), + c64::new(2.9998484902996343, 2.0506177173361895), + c64::new(-2.786643678486076, -5.270377319461184), + c64::new(-6.332372870212514, -2.7133499614347927), + c64::new(3.7738891982330425, -3.402493113787302), + c64::new(-1.7781003206660229, 5.718899556657502), + c64::new(0.8339112283930712, 8.762743527504561), + c64::new(2.0569636038790033, 3.146955624252856), + c64::new(-2.7031084237762477, -1.200977331949071), + c64::new(-2.4430054443688007, -0.4059754102625929), + c64::new(-1.185157131906472, 1.4890565683775576), + c64::new(3.590908753088401, 3.0941791756116297), + c64::new(-2.641272411857855, -2.0351254037035407), + c64::new(-1.5648879937245017, -0.011039956994125966), + c64::new(4.899209965296281, -8.154920908424456), + c64::new(1.7041853352045995, 3.986228161505659) + ], + [ + c64::new(0.7098321944918885, 1.8899960715632942), + c64::new(-1.2876995397035196, -0.4974779624183936), + c64::new(1.5979867371106224, 2.18302312402332), + c64::new(2.9317303907773473, 0.1464447495528882), + c64::new(-1.1309352560779358, 2.3978248444784915), + c64::new(0.833438793503704, -3.7067752030654693), + c64::new(-1.6585616430787469, -3.8407526888470445), + c64::new(-0.8982175686022119, -1.0025945939580605), + c64::new(1.5575446743219135, -0.38889435787413484), + c64::new(1.1772769305079582, -0.6778245090880448), + c64::new(-0.15792138663043898, -0.9309931883669093), + c64::new(-2.66343386833843, -0.5715520673153744), + c64::new(0.9557077119414871, 0.8656183410600455), + c64::new(1.441553858828804, -0.3040776192978809), + c64::new(-1.7175922339016105, 5.159005168251366), + c64::new(-0.43154772399617425, -1.4184624078914763) + ], + [ + c64::new(0.7067369639484306, -0.11677774934896479), + c64::new(-1.683473630140353, 1.3019416191833875), + c64::new(2.6652317489455895, -0.7333921335041848), + c64::new(2.355056156923994, -3.1632176945121766), + c64::new(0.8497205963164971, 2.5732497308004616), + c64::new(-2.0154697095526872, -1.3128840122359084), + c64::new(-4.4422946390839915, -0.5506581875054117), + c64::new(-1.479170402306385, 1.173540079477825), + c64::new(1.0644798429183335, -1.3106697173813893), + c64::new(0.5830433337637937, -1.4688418402948535), + c64::new(-0.45785525478165967, -0.4195897948227327), + c64::new(-1.8983086728039444, 0.8913322959702731), + c64::new(1.2560059050813184, -1.6095836772765104), + c64::new(0.6008114921264225, -0.6653443016353698), + c64::new(3.2493010280194836, 2.997293768718692), + c64::new(-2.4054246583174637, 0.9739433780581721) + ], + [ + c64::new(1.4239159811426512, 0.8492138766545766), + c64::new(-1.3991192281292535, -0.43373158118485333), + c64::new(2.174929687527398, 1.3879134381263984), + c64::new(2.0535605031600834, 0.30290545784992706), + c64::new(-0.9507775086955514, 1.5330343629000758), + c64::new(-0.5995075478759522, -2.5158922223417863), + c64::new(-1.961409868797519, -3.0719263946349766), + c64::new(-1.1532845082224947, -1.3535437102253036), + c64::new(1.2536363004520334, -0.39574740266623964), + c64::new(0.878499113182353, -0.015575249533619606), + c64::new(0.14506197677793486, -0.6095788748412565), + c64::new(-2.0006000760113993, -0.36270160074571994), + c64::new(1.5014401345887063, 0.9124332027445035), + c64::new(0.6224797643970621, -0.41935426116973706), + c64::new(-0.08903298168132473, 4.263757880031095), + c64::new(-0.9093231367777892, -1.4242001367849448) + ], + [ + c64::new(-1.252615813160683, -3.0046048352701398), + c64::new(3.6937074952945474, 0.8629851151639598), + c64::new(-4.424158892742261, -4.819197382549424), + c64::new(-7.548435879349295, 0.4123445749516162), + c64::new(2.989181436478559, -5.518161108503455), + c64::new(-0.3971571217791885, 7.354687327960869), + c64::new(4.5312026486551575, 8.878293476735566), + c64::new(2.5972543971580055, 1.7949272103216405), + c64::new(-4.724266957695443, -0.17086789993622736), + c64::new(-3.7587317704878878, 1.225753184207096), + c64::new(-0.1388339158026093, 2.338843606827039), + c64::new(5.43195196867949, 2.4055348833010304), + c64::new(-2.3261483013555857, -0.9883963113909155), + c64::new(-2.709252611717642, 0.34074394535010155), + c64::new(3.0393053848361538, -11.037706160752986), + c64::new(3.3744887199210103, 2.3998987532495195) + ], + [ + c64::new(-0.6653756801965123, -0.1036593934409793), + c64::new(1.5205374386230774, -0.8841741070156326), + c64::new(-2.394043208594155, -0.8045471482882556), + c64::new(-2.814975456215832, 1.5932935660927487), + c64::new(0.11248961579460955, -2.369649518617432), + c64::new(0.993651193365689, 2.0349442465753826), + c64::new(3.4164791453229926, 2.358359283330307), + c64::new(1.5385342334885286, 0.1632682172084824), + c64::new(-1.7279134524806181, 0.7690062252967498), + c64::new(-0.8684676203130114, 1.1327657637345285), + c64::new(0.587069210347443, 1.0816539677189212), + c64::new(1.953176918767738, -0.17502354955401722), + c64::new(-1.4319425641758028, 0.39723015900232195), + c64::new(-0.7953530169702814, 0.3018693314086077), + c64::new(-1.1095546459915613, -3.91769613866538), + c64::new(2.0871995699347408, -0.34994087477466085) + ], + [ + c64::new(1.0428659661743744, -1.627000038293237), + c64::new(0.8935286888615769, 2.9905846717173947), + c64::new(1.0755083030676815, -4.0177728554390955), + c64::new(-3.0359003674327036, -5.027961108809421), + c64::new(4.673494661011926, 0.21701930072847653), + c64::new(-3.927807792882151, 2.0369036075192652), + c64::new(-4.030989423360226, 6.59990309935656), + c64::new(0.26335030205361054, 2.5483491799122), + c64::new(-0.49435775657652165, -3.1792386277930893), + c64::new(-2.0204286938145746, -2.9045120248583745), + c64::new(-1.8049059670607743, 0.49938863667012034), + c64::new(-0.22188049558232192, 4.281258983711006), + c64::new(-1.2928614778042757, -2.1001096045646483), + c64::new(-0.693613251969597, -1.5908493889838762), + c64::new(7.277239343072571, -2.0979995176164756), + c64::new(-0.12787942570107225, 3.8509947575701413) + ], + [ + c64::new(-0.9030035567896482, -1.285148731288444), + c64::new(2.226551915191729, 0.9198169968078413), + c64::new(-2.570296540921288, -2.542050486376615), + c64::new(-4.000007296285765, -0.3495888233487804), + c64::new(1.5322785903488643, -2.8364745006616494), + c64::new(0.05457640225206961, 3.8086232939559075), + c64::new(2.4132022722967004, 5.021382620580456), + c64::new(1.575482666377672, 1.2949366880751572), + c64::new(-2.2798286177934255, -0.23202522027280154), + c64::new(-1.9453429360673578, 0.9086022093574165), + c64::new(-0.36782592552388826, 1.1282298110284499), + c64::new(2.7909851000471724, 1.1901581765442475), + c64::new(-2.0449147956434595, -1.1027820508657102), + c64::new(-0.957402145821622, 0.2219440862791633), + c64::new(1.5311134129121113, -6.106774871371683), + c64::new(1.9750485314567656, 1.6036470887259389) + ], + [ + c64::new(-1.647229157806129, 0.788861290474623), + c64::new(0.6148722774868443, -1.7793373228544274), + c64::new(-2.3174465361956824, 2.1040057464401194), + c64::new(-0.4420474653129551, 2.7791610362888512), + c64::new(-2.069639274156649, -1.155550866182406), + c64::new(3.1026708695195566, 0.207718180381943), + c64::new(4.081353810713267, -1.5985484153287735), + c64::new(1.7387648401193387, -1.1811673575712405), + c64::new(0.7078217807594732, 1.7697384542515109), + c64::new(0.3896853684770497, 1.072939257271109), + c64::new(0.7616221623406603, 0.18279301165029538), + c64::new(0.7572481411523853, -2.699298614275696), + c64::new(-0.9170378278185156, 1.6906134084938873), + c64::new(0.5337827371636874, 0.8724411387783653), + c64::new(-5.276156593076657, -1.156578578366272), + c64::new(1.1874546730575855, -0.8966724612603467) + ], + [ + c64::new(0.9313726674656437, 3.0718691038815065), + c64::new(-3.0215888657993535, -2.0907788661535767), + c64::new(2.3356696345973313, 5.596976766853978), + c64::new(6.609499575230608, 2.6858335287924033), + c64::new(-4.114103323857318, 3.5840334133467096), + c64::new(2.1436259729553644, -6.129539492479053), + c64::new(-0.7790731721951514, -9.373889045466344), + c64::new(-2.2086446458204954, -2.9265263723081896), + c64::new(3.37294887937093, 1.4327653829016842), + c64::new(2.4774112974245375, 0.4697634315107313), + c64::new(1.204678119876694, -1.9976531671605633), + c64::new(-3.4061059135734686, -3.4243155383617183), + c64::new(2.724539266310142, 1.9301312701773405), + c64::new(1.6583588462898349, 0.20482421031998577), + c64::new(-5.3022656451878865, 8.359424149418086), + c64::new(-2.0536748485066942, -3.4518838281996613) + ], + [ + c64::new(-0.3400038949681029, 0.45128888584303883), + c64::new(0.302834399057287, -0.9502246062639901), + c64::new(-0.6068476034774675, 1.6660485163780343), + c64::new(0.19115500499373536, 2.173061763191289), + c64::new(-0.9538667784471226, -0.9977729990492099), + c64::new(1.2388604002033339, -0.8371530287267638), + c64::new(2.553109099137643, -1.2809339536966684), + c64::new(0.23292377813052328, -0.8646085600315921), + c64::new(-0.5896410383024675, 1.033670974129205), + c64::new(-0.022540559446400155, 1.0698412350894706), + c64::new(0.4127485917223416, 0.11538759781118935), + c64::new(0.697367920179159, -0.952172194056835), + c64::new(-0.27389037619737855, 1.0604722509205344), + c64::new(0.02919900182455054, 0.13786814316808738), + c64::new(-2.259132868623872, -0.5163004615257191), + c64::new(0.5043165094609598, -1.1438622215101546) + ], + [ + c64::new(3.212745888447596, 2.913982534731071), + c64::new(-5.707956217016549, 0.561560814070595), + c64::new(8.47650395832462, 3.578564034019665), + c64::new(10.124896409465327, -4.184969047203229), + c64::new(-1.170061311708083, 8.686707346447756), + c64::new(-3.313821036277578, -8.896990083424384), + c64::new(-10.7526969515756, -9.519756126873094), + c64::new(-4.9127993200515085, -1.1059785971618572), + c64::new(5.749895822184144, -1.1833266940320586), + c64::new(4.508761401024904, -3.0900236338081783), + c64::new(-0.6673814135567939, -3.3461695674261254), + c64::new(-7.889103766073837, -0.4075117015034545), + c64::new(4.399990035115406, -0.34774164990939216), + c64::new(2.6333247662163712, -1.4466247972859447), + c64::new(2.014420530325137, 15.024581881327247), + c64::new(-5.8578217429440445, -1.7952789388416184) + ], + ]; + + for (actual, expected) in inv.iter().zip(expected_inv.iter()) { + assert_abs_diff_eq!(actual, expected, epsilon = 1e-6); + } + Ok(()) + } +} diff --git a/src/Simulation/qdk_sim_rs/src/linalg.rs b/src/Simulation/qdk_sim_rs/src/linalg/mod.rs similarity index 73% rename from src/Simulation/qdk_sim_rs/src/linalg.rs rename to src/Simulation/qdk_sim_rs/src/linalg/mod.rs index 98e41a02f8c..d324fd020ad 100644 --- a/src/Simulation/qdk_sim_rs/src/linalg.rs +++ b/src/Simulation/qdk_sim_rs/src/linalg/mod.rs @@ -3,12 +3,29 @@ //! Provides common linear algebra functions and traits. -use ndarray::{Array, Array1, Array2, ArrayView1, ArrayView2}; +use ndarray::{Array, Array2, ArrayView2}; use num_traits::Zero; -use std::{convert::TryInto, ops::Mul}; +use std::convert::TryInto; use crate::{common_matrices::nq_eye, C64}; +// Define the public API surface for submodules. + +mod tensor; +pub use tensor::*; + +mod inv; +pub use inv::*; + +mod pow; +pub use pow::*; + +pub mod decompositions; + +// Define private modules as well. + +mod array_ext; + /// Represents types that have hermitian conjugates (e.g.: $A^\dagger$ for /// a matrix $A$ is defined as the complex conjugate transpose of $A$, /// $(A^\dagger)\_{ij} = A\_{ji}^*$). @@ -53,73 +70,6 @@ impl ConjBy for Array2 { } } -/// The tensor product operator ($\otimes$). -pub trait Tensor { - /// The resulting type after applying the tensor product. - type Output; - - /// Performs the tensor product. - /// - /// # Example - /// ``` - /// // TODO - /// ``` - fn tensor(self, rhs: Rhs) -> Self::Output; -} - -impl, T: Copy + Mul> Tensor for ArrayView1<'_, T> { - type Output = Array1; - fn tensor(self, other: Other) -> Self::Output { - let other: Self = other.into(); - let unflat = Array::from_shape_fn((self.shape()[0], other.shape()[0]), |(i, j)| { - self[(i)] * other[(j)] - }); - unflat - .into_shape(self.shape()[0] * other.shape()[0]) - .unwrap() - } -} - -impl, T: Copy + Mul> Tensor for &Array1 { - type Output = Array1; - - fn tensor(self, other: Other) -> Self::Output { - let other: Self = other.into(); - self.view().tensor(other) - } -} - -impl, T: Copy + Mul> Tensor for ArrayView2<'_, T> { - type Output = Array2; - fn tensor(self, other: Other) -> Self::Output { - let other: Self = other.into(); - let unflat = Array::from_shape_fn( - ( - self.shape()[0], - other.shape()[0], - self.shape()[1], - other.shape()[1], - ), - |(i, j, k, l)| self[(i, k)] * other[(j, l)], - ); - unflat - .into_shape(( - self.shape()[0] * other.shape()[0], - self.shape()[1] * other.shape()[1], - )) - .unwrap() - } -} - -impl, T: Copy + Mul> Tensor for &Array2 { - type Output = Array2; - - fn tensor(self, other: Other) -> Self::Output { - let other: Self = other.into(); - self.view().tensor(other).to_owned() - } -} - /// Represents types for which the trace can be computed. pub trait Trace { /// The type returned by the trace. diff --git a/src/Simulation/qdk_sim_rs/src/linalg/pow.rs b/src/Simulation/qdk_sim_rs/src/linalg/pow.rs new file mode 100644 index 00000000000..78df3335db6 --- /dev/null +++ b/src/Simulation/qdk_sim_rs/src/linalg/pow.rs @@ -0,0 +1,179 @@ +use cauchy::Scalar; +use ndarray::{linalg::Dot, Array, Array2, Ix2}; +use num_traits::Zero; + +use crate::linalg::Inv; + +/// Types for which there is a meaningful matrix identity value. +pub trait Identity: Sized { + /// Returns a matrix identity of a given size. + fn eye(size: usize) -> Self; + + /// Returns a matrix identity the same size as a given input. + fn eye_like(&self) -> Self; +} + +impl Identity for Array +where + A: Scalar + Zero + Clone, +{ + fn eye(size: usize) -> Self { + Self::eye(size) + } + + fn eye_like(&self) -> Self { + as Identity>::eye(std::cmp::min(self.shape()[0], self.shape()[1])) + } +} + +/// Types that support matrix powers $A^{n}$. +pub trait MatrixPower: Sized { + /// Errors that can result from the [`matrix_power`] method. + /// + /// [`matrix_power`]: MatrixPower::matrix_power + type Error; + + /// Returns an integer power of the given matrix. + fn matrix_power(&self, pow: i64) -> Result; +} + +impl MatrixPower for M +where + M: Inv + Dot + Identity, +{ + type Error = M::Error; + + fn matrix_power(&self, pow: i64) -> Result { + if pow < 0 { + return self.inv()?.matrix_power(-pow); + } + + let mut result = self.eye_like(); + + // TODO(perf): use binary exponentiation; the simple for loop is just to + // quickly get something that can be tested. Once we have + // good unit test coverage, we can use that to make sure + // binary exponentiation works right. + for _ in 0..pow { + result = result.dot(self); + } + + Ok(result) + } +} + +#[cfg(test)] +mod tests { + use approx::assert_abs_diff_eq; + use cauchy::c64; + use ndarray::{array, Array2}; + + use crate::{c64, error::QdkSimError, linalg::MatrixPower}; + + #[test] + fn pauli_x_squares_to_ident() -> Result<(), QdkSimError> { + let x: Array2 = array![[0.0, 1.0], [1.0, 0.0]]; + let actual = x.matrix_power(2)?; + let expected = array![[1.0, 0.0], [0.0, 1.0]]; + for (actual, expected) in actual.iter().zip(expected.iter()) { + assert_abs_diff_eq!(actual, expected, epsilon = 1e-6); + } + Ok(()) + } + + #[test] + fn sqrtx_to_the_fourth_is_ident() -> Result<(), QdkSimError> { + let sqrtx: Array2 = c64!(0.5) + * array![ + [c64!(1.0 + 1.0 i), c64!(1.0 - 1.0 i)], + [c64!(1.0 - 1.0 i), c64!(1.0 + 1.0 i)] + ]; + let actual = sqrtx.matrix_power(4)?; + let expected = array![[c64!(1.0), c64!(0.0)], [c64!(0.0), c64!(1.0)]]; + for (actual, expected) in actual.iter().zip(expected.iter()) { + assert_abs_diff_eq!(actual, expected, epsilon = 1e-6); + } + Ok(()) + } + + #[test] + fn seventh_root_of_x_to_the_seventh_is_x() -> Result<(), QdkSimError> { + let seventh_root_x: Array2 = array![ + [ + c64::new(0.9504844339512094, 0.21694186955877903), + c64::new(0.04951556604879037, -0.21694186955877903) + ], + [ + c64::new(0.04951556604879037, -0.21694186955877903), + c64::new(0.9504844339512094, 0.21694186955877903) + ], + ]; + let actual = seventh_root_x.matrix_power(7)?; + let expected = array![[c64!(0.0), c64!(1.0)], [c64!(1.0), c64!(0.0)]]; + for (actual, expected) in actual.iter().zip(expected.iter()) { + assert_abs_diff_eq!(actual, expected, epsilon = 1e-6); + } + Ok(()) + } + + #[test] + fn rand_4x4_c64_pow() -> Result<(), QdkSimError> { + let arr: Array2 = array![ + [ + c64::new(0.4197614535218097, 0.6393638574841801), + c64::new(0.021249529076351803, 0.4916053336160465), + c64::new(0.7598776596987344, 0.06591304637672646), + c64::new(0.5633311068416154, 0.031336286212855446) + ], + [ + c64::new(0.0269250407737347, 0.20379518145743047), + c64::new(0.30555508065274883, 0.15808290828587157), + c64::new(0.5104720129393309, 0.6948521818460662), + c64::new(0.19865943993696022, 0.9287971169482773) + ], + [ + c64::new(0.2905701510819163, 0.09833520319062727), + c64::new(0.23402610846324, 0.8393135547393529), + c64::new(0.5401748112924776, 0.9821393064000913), + c64::new(0.0927979325237448, 0.2477695883264862) + ], + [ + c64::new(0.7717353949875485, 0.42140532067756675), + c64::new(0.5528553121605853, 0.9640588921494402), + c64::new(0.7349192227517806, 0.5792000992940656), + c64::new(0.12338346370340081, 0.07289997402380066) + ], + ]; + let actual = arr.matrix_power(13)?; + let expected = array![ + [ + c64::new(7238.141066326559, -9911.274649764955), + c64::new(17648.478137710346, -9418.948356812798), + c64::new(18297.504198477975, -18154.35549609193), + c64::new(10159.63982578618, -5485.709652173136) + ], + [ + c64::new(12959.04402150731, -4661.525869682588), + c64::new(22325.641696782357, 2341.2319188120773), + c64::new(28420.99091079291, -5370.043643156111), + c64::new(12892.072719114653, 1288.7133317204193) + ], + [ + c64::new(12882.930101604788, -2963.4083140513117), + c64::new(21010.506984145948, 4779.901376069587), + c64::new(27700.037816647793, -1873.6300238965507), + c64::new(12139.8605952774, 2699.9073450820306) + ], + [ + c64::new(12071.951917001059, -9465.554479542623), + c64::new(24427.99605457002, -5339.098219764661), + c64::new(28174.33364730465, -15627.076130106801), + c64::new(14083.953080905429, -3149.97121499003) + ], + ]; + for (actual, expected) in actual.iter().zip(expected.iter()) { + assert_abs_diff_eq!(actual, expected, epsilon = 1e-6); + } + Ok(()) + } +} diff --git a/src/Simulation/qdk_sim_rs/src/linalg/tensor.rs b/src/Simulation/qdk_sim_rs/src/linalg/tensor.rs new file mode 100644 index 00000000000..ef5d9ff7574 --- /dev/null +++ b/src/Simulation/qdk_sim_rs/src/linalg/tensor.rs @@ -0,0 +1,73 @@ +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. + +use std::ops::Mul; + +use ndarray::{Array, Array1, Array2, ArrayView1, ArrayView2}; + +/// The tensor product operator ($\otimes$). +pub trait Tensor { + /// The resulting type after applying the tensor product. + type Output; + + /// Performs the tensor product. + /// + /// # Example + /// ``` + /// // TODO + /// ``` + fn tensor(self, rhs: Rhs) -> Self::Output; +} + +impl, T: Copy + Mul> Tensor for ArrayView1<'_, T> { + type Output = Array1; + fn tensor(self, other: Other) -> Self::Output { + let other: Self = other.into(); + let unflat = Array::from_shape_fn((self.shape()[0], other.shape()[0]), |(i, j)| { + self[(i)] * other[(j)] + }); + unflat + .into_shape(self.shape()[0] * other.shape()[0]) + .unwrap() + } +} + +impl, T: Copy + Mul> Tensor for &Array1 { + type Output = Array1; + + fn tensor(self, other: Other) -> Self::Output { + let other: Self = other.into(); + self.view().tensor(other) + } +} + +impl, T: Copy + Mul> Tensor for ArrayView2<'_, T> { + type Output = Array2; + fn tensor(self, other: Other) -> Self::Output { + let other: Self = other.into(); + let unflat = Array::from_shape_fn( + ( + self.shape()[0], + other.shape()[0], + self.shape()[1], + other.shape()[1], + ), + |(i, j, k, l)| self[(i, k)] * other[(j, l)], + ); + unflat + .into_shape(( + self.shape()[0] * other.shape()[0], + self.shape()[1] * other.shape()[1], + )) + .unwrap() + } +} + +impl, T: Copy + Mul> Tensor for &Array2 { + type Output = Array2; + + fn tensor(self, other: Other) -> Self::Output { + let other: Self = other.into(); + self.view().tensor(other).to_owned() + } +} diff --git a/src/Simulation/qdk_sim_rs/src/noise_model.rs b/src/Simulation/qdk_sim_rs/src/noise_model.rs index 50fdd79e12a..80c153781b7 100644 --- a/src/Simulation/qdk_sim_rs/src/noise_model.rs +++ b/src/Simulation/qdk_sim_rs/src/noise_model.rs @@ -1,5 +1,6 @@ use crate::chp_decompositions::ChpOperation; use crate::common_matrices; +use crate::error::QdkSimError; use crate::instrument::Instrument; use crate::linalg::HasDagger; use crate::processes::Process; @@ -90,11 +91,11 @@ impl NoiseModel { /// # use qdk_sim::NoiseModel; /// let noise_model = NoiseModel::get_by_name("ideal"); /// ``` - pub fn get_by_name(name: &str) -> Result { + pub fn get_by_name(name: &str) -> Result { match name { "ideal" => Ok(NoiseModel::ideal()), "ideal_stabilizer" => Ok(NoiseModel::ideal_stabilizer()), - _ => Err(format!("Unrecognized noise model name {}.", name)), + _ => Err(QdkSimError::InvalidNoiseModel(name.to_string())), } } diff --git a/src/Simulation/qdk_sim_rs/src/processes/apply.rs b/src/Simulation/qdk_sim_rs/src/processes/apply.rs index cd5cca76b5a..97cbcc843c6 100644 --- a/src/Simulation/qdk_sim_rs/src/processes/apply.rs +++ b/src/Simulation/qdk_sim_rs/src/processes/apply.rs @@ -9,8 +9,8 @@ use ndarray::{Array, Array2, Array3, ArrayView2, Axis}; use rand::{distributions::WeightedIndex, prelude::Distribution, thread_rng}; use crate::{ - chp_decompositions::ChpOperation, linalg::ConjBy, log, log_as_err, Pauli, Process, - ProcessData::*, State, StateData::*, Tableau, C64, + chp_decompositions::ChpOperation, error::QdkSimError, linalg::ConjBy, log, log_as_err, Pauli, + Process, ProcessData::*, State, StateData::*, Tableau, VariantName, C64, }; use super::promote_pauli_channel; @@ -18,12 +18,12 @@ use super::promote_pauli_channel; impl Process { /// Applies this process to a quantum register with a given /// state, returning the new state of that register. - pub fn apply(&self, state: &State) -> Result { + pub fn apply(&self, state: &State) -> Result { if state.n_qubits != self.n_qubits { - return Err(format!( - "Channel acts on {} qubits, but was applied to {}-qubit state.", - self.n_qubits, state.n_qubits - )); + return Err(QdkSimError::WrongNumberOfQubits { + expected: self.n_qubits, + actual: state.n_qubits, + }); } match &self.data { @@ -39,13 +39,16 @@ impl Process { Ok(acc_state) } ChpDecomposition(_operations) => todo!(), - Unsupported => Err("Unsupported quantum process.".to_string()), + Unsupported => Err(QdkSimError::UnsupportedApply { + channel_variant: self.variant_name(), + state_variant: state.variant_name(), + }), } } /// Applies this process to the given qubits in a register with a given /// state, returning the new state of that register. - pub fn apply_to(&self, idx_qubits: &[usize], state: &State) -> Result { + pub fn apply_to(&self, idx_qubits: &[usize], state: &State) -> Result { // If we have a sequence, we can apply each in turn and exit early. if let Sequence(channels) = &self.data { // TODO[perf]: eliminate the extraneous clone here. @@ -137,7 +140,7 @@ fn apply_chp_decomposition_to( n_qubits: usize, idx_qubits: &[usize], tableau: &Tableau, -) -> Result { +) -> Result { let mut new_tableau = tableau.clone(); for operation in operations { match *operation { @@ -155,22 +158,25 @@ fn apply_chp_decomposition_to( }) } -pub(crate) fn apply_unitary(u: &Array2, state: &State) -> Result { +pub(crate) fn apply_unitary(u: &Array2, state: &State) -> Result { Ok(State { n_qubits: state.n_qubits, data: match &state.data { Pure(psi) => Pure(u.dot(psi)), Mixed(rho) => Mixed(rho.conjugate_by(&u.into())), Stabilizer(_tableau) => { - return Err( + return Err(QdkSimError::NotYetImplemented( "TODO: Promote stabilizer state to state vector and recurse.".to_string(), - ) + )) } }, }) } -pub(crate) fn apply_kraus_decomposition(ks: &Array3, state: &State) -> Result { +pub(crate) fn apply_kraus_decomposition( + ks: &Array3, + state: &State, +) -> Result { Ok(State { n_qubits: state.n_qubits, data: match &state.data { @@ -195,9 +201,9 @@ pub(crate) fn apply_kraus_decomposition(ks: &Array3, state: &State) -> Resu sum }), Stabilizer(_tableau) => { - return Err( + return Err(QdkSimError::NotYetImplemented( "TODO: Promote stabilizer state to state vector and recurse.".to_string(), - ) + )) } }, }) @@ -206,7 +212,7 @@ pub(crate) fn apply_kraus_decomposition(ks: &Array3, state: &State) -> Resu pub(crate) fn apply_pauli_channel( paulis: &[(f64, Vec)], state: &State, -) -> Result { +) -> Result { Ok(State { n_qubits: state.n_qubits, data: match &state.data { diff --git a/src/Simulation/qdk_sim_rs/src/processes/mod.rs b/src/Simulation/qdk_sim_rs/src/processes/mod.rs index 4aacc049d96..073ab37d52b 100644 --- a/src/Simulation/qdk_sim_rs/src/processes/mod.rs +++ b/src/Simulation/qdk_sim_rs/src/processes/mod.rs @@ -6,14 +6,15 @@ mod apply; use crate::chp_decompositions::ChpOperation; use crate::linalg::{extend_one_to_n, extend_two_to_n, zeros_like}; use crate::processes::ProcessData::{KrausDecomposition, MixedPauli, Unitary}; -use crate::NoiseModel; use crate::QubitSized; use crate::C64; use crate::{AsUnitary, Pauli}; +use crate::{NoiseModel, VariantName}; use itertools::Itertools; use ndarray::{Array, Array2, Array3, Axis, NewAxis}; use num_complex::Complex; use num_traits::{One, Zero}; + use serde::{Deserialize, Serialize}; use std::convert::TryInto; use std::ops::Add; @@ -55,6 +56,19 @@ pub enum ProcessData { Unsupported, } +impl VariantName for Process { + fn variant_name(&self) -> &'static str { + match &self.data { + MixedPauli(_) => "MixedPauli", + Unitary(_) => "Unitary", + KrausDecomposition(_) => "KrausDecomposition", + ProcessData::Sequence(_) => "Sequence", + ProcessData::ChpDecomposition(_) => "ChpDecomposition", + ProcessData::Unsupported => "Unsupported", + } + } +} + impl Process { /// Returns a new Pauli channel, given a mixture of Pauli operators. pub fn new_pauli_channel(data: T) -> Self { diff --git a/src/Simulation/qdk_sim_rs/src/states.rs b/src/Simulation/qdk_sim_rs/src/states.rs index 71f36e51b6b..cdb6f44adf9 100644 --- a/src/Simulation/qdk_sim_rs/src/states.rs +++ b/src/Simulation/qdk_sim_rs/src/states.rs @@ -7,6 +7,7 @@ use crate::linalg::Trace; use crate::states::StateData::{Mixed, Pure, Stabilizer}; use crate::tableau::Tableau; use crate::QubitSized; +use crate::VariantName; use crate::C64; use core::fmt::Display; use ndarray::{Array1, Array2, Axis}; @@ -49,6 +50,16 @@ impl Display for State { } } +impl VariantName for State { + fn variant_name(&self) -> &'static str { + match self.data { + StateData::Pure(_) => "Pure", + StateData::Mixed(_) => "Mixed", + StateData::Stabilizer(_) => "Stabilizer", + } + } +} + impl Display for StateData { fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::result::Result<(), std::fmt::Error> { match self { diff --git a/src/Simulation/qdk_sim_rs/src/tableau.rs b/src/Simulation/qdk_sim_rs/src/tableau.rs index ccbf5ddd98f..0a65828e456 100644 --- a/src/Simulation/qdk_sim_rs/src/tableau.rs +++ b/src/Simulation/qdk_sim_rs/src/tableau.rs @@ -3,7 +3,10 @@ use std::fmt::Display; -use crate::utils::{set_row_to_row_sum, set_vec_to_row_sum, swap_columns}; +use crate::{ + error::QdkSimError, + utils::{set_row_to_row_sum, set_vec_to_row_sum, swap_columns}, +}; use ndarray::{s, Array, Array1, Array2}; use serde::{Deserialize, Serialize}; @@ -91,16 +94,18 @@ impl Tableau { /// /// If the assertion would pass, `Ok(())` is returned, otherwise an [`Err`] /// describing the assertion failure is returned. - pub fn assert_meas(&self, idx_target: usize, expected: bool) -> Result<(), String> { - let actual = self.determinstic_result(idx_target).ok_or(format!( - "Expected {}, but measurement result would be random.", - expected - ))?; + pub fn assert_meas(&self, idx_target: usize, expected: bool) -> Result<(), QdkSimError> { + let actual = self + .determinstic_result(idx_target) + .ok_or(QdkSimError::MiscError(format!( + "Expected {}, but measurement result would be random.", + expected + )))?; if actual != expected { - Err(format!( + Err(QdkSimError::MiscError(format!( "Expected {}, but measurement result would actually be {}.", expected, actual - )) + ))) } else { Ok(()) } diff --git a/src/Simulation/qdk_sim_rs/src/utils.rs b/src/Simulation/qdk_sim_rs/src/utils.rs index 753ea0f5d46..ea4ddaf7b31 100644 --- a/src/Simulation/qdk_sim_rs/src/utils.rs +++ b/src/Simulation/qdk_sim_rs/src/utils.rs @@ -4,6 +4,14 @@ use ndarray::{s, Array1, Array2, ArrayView1}; use num_complex::Complex; +use crate::error::QdkSimError; + +/// Used for enums whose variant names are meaningful. +/// In the future, this should be derived and not hard-coded. +pub(crate) trait VariantName { + fn variant_name(&self) -> &'static str; +} + /// A complex number with two 64-bit floating point fields. /// That is, the analogy of [`f64`] to complex values. pub type C64 = Complex; @@ -31,9 +39,9 @@ fn log_message(msg: &str) { } /// Prints a message as an error, and returns it as a [`Result`]. -pub fn log_as_err(msg: String) -> Result { +pub fn log_as_err(msg: String) -> Result { log(&msg); - Err(msg) + Err(QdkSimError::MiscError(msg)) } /// Prints a message as an error. @@ -115,3 +123,34 @@ pub fn phase_product(row1: &ArrayView1, row2: &ArrayView1) -> bool { ((if r1 { 2 } else { 0 }) + (if r2 { 2 } else { 0 }) + acc) % 4 == 2 } + +/// Macro for quickly declaring complex numbers of type [`cauchy::c64`]. +/// +/// # Example +/// ``` +/// # use qdk_sim::c64; +/// # use cauchy::c64; +/// let re_unit = c64!(1.0); +/// assert_eq!(re_unit, c64::new(1.0, 0.0)); +/// +/// let im_unit = c64!(1.0 i); +/// assert_eq!(im_unit, c64::new(0.0, 1.0)); +/// +/// let z = c64!(1.0 + 2.0 i); +/// assert_eq!(z, c64::new(1.0, 2.0)); +/// ``` +#[macro_export] +macro_rules! c64 { + ($re:literal) => { + c64::new($re, 0.0) + }; + ($re:literal + $im:literal i) => { + c64::new($re, $im) + }; + ($re:literal - $im:literal i) => { + c64::new($re, -$im) + }; + ($im:literal i) => { + c64::new(0.0, $im) + }; +} diff --git a/src/Simulation/qdk_sim_rs/tests/chp_simulation_tests.rs b/src/Simulation/qdk_sim_rs/tests/chp_simulation_tests.rs index c709d664cd3..8bdd2342f26 100644 --- a/src/Simulation/qdk_sim_rs/tests/chp_simulation_tests.rs +++ b/src/Simulation/qdk_sim_rs/tests/chp_simulation_tests.rs @@ -1,10 +1,10 @@ // Copyright (c) Microsoft Corporation. // Licensed under the MIT License. -use qdk_sim::{Pauli, Process, State, Tableau}; +use qdk_sim::{error::QdkSimError, Pauli, Process, State, Tableau}; #[test] -fn pauli_channel_applies_correctly() -> Result<(), String> { +fn pauli_channel_applies_correctly() -> Result<(), QdkSimError> { let x = Process::new_pauli_channel(Pauli::X); let state = State::new_stabilizer(1); let output_state = x.apply(&state)?;