diff --git a/Cargo.lock b/Cargo.lock index 996cf6a64819..c09000cabd2a 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -30,6 +30,16 @@ dependencies = [ "log", ] +[[package]] +name = "approx" +version = "0.5.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "cab112f0a86d568ea0e627cc1d6be74a1e9cd55214684db5561995f6dad897c6" +dependencies = [ + "num-complex", + "num-traits", +] + [[package]] name = "ariadne" version = "0.3.0" @@ -793,6 +803,7 @@ version = "0.15.6" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "adb12d4e967ec485a5f71c6311fe28158e9d6f4bc4a447b474184d0f91a8fa32" dependencies = [ + "approx", "matrixmultiply", "num-complex", "num-integer", @@ -1078,6 +1089,7 @@ dependencies = [ "pyo3-build-config", "pyo3-ffi", "pyo3-macros", + "smallvec 1.13.1", "unindent", ] @@ -1149,6 +1161,7 @@ name = "qiskit_accelerate" version = "1.1.0" dependencies = [ "ahash", + "approx", "faer", "faer-core", "hashbrown 0.14.3", diff --git a/crates/accelerate/Cargo.toml b/crates/accelerate/Cargo.toml index 5c4ba81496d7..fa855aede219 100644 --- a/crates/accelerate/Cargo.toml +++ b/crates/accelerate/Cargo.toml @@ -33,11 +33,15 @@ features = ["union"] [dependencies.pyo3] workspace = true -features = ["hashbrown", "indexmap", "num-complex", "num-bigint"] +features = ["hashbrown", "indexmap", "num-complex", "num-bigint", "smallvec"] [dependencies.ndarray] version = "^0.15.6" -features = ["rayon"] +features = ["rayon", "approx-0_5"] + +[dependencies.approx] +version = "0.5" +features = ["num-complex"] [dependencies.hashbrown] workspace = true diff --git a/crates/accelerate/src/euler_one_qubit_decomposer.rs b/crates/accelerate/src/euler_one_qubit_decomposer.rs index e74f70006696..6af60400a910 100644 --- a/crates/accelerate/src/euler_one_qubit_decomposer.rs +++ b/crates/accelerate/src/euler_one_qubit_decomposer.rs @@ -11,9 +11,11 @@ // that they have been altered from the originals. #![allow(clippy::too_many_arguments)] +#![allow(clippy::upper_case_acronyms)] use hashbrown::HashMap; use num_complex::{Complex64, ComplexFloat}; +use smallvec::{smallvec, SmallVec}; use std::cmp::Ordering; use std::f64::consts::PI; @@ -27,7 +29,7 @@ use numpy::PyReadonlyArray2; use crate::utils::SliceOrInt; -const DEFAULT_ATOL: f64 = 1e-12; +pub const ANGLE_ZERO_EPSILON: f64 = 1e-12; #[pyclass(module = "qiskit._accelerate.euler_one_qubit_decomposer")] pub struct OneQubitGateErrorMap { @@ -61,11 +63,13 @@ impl OneQubitGateErrorMap { #[pyclass(sequence)] pub struct OneQubitGateSequence { - gates: Vec<(String, Vec)>, + pub gates: Vec<(String, SmallVec<[f64; 3]>)>, #[pyo3(get)] - global_phase: f64, + pub global_phase: f64, } +type OneQubitGateSequenceState = (Vec<(String, SmallVec<[f64; 3]>)>, f64); + #[pymethods] impl OneQubitGateSequence { #[new] @@ -75,11 +79,11 @@ impl OneQubitGateSequence { global_phase: 0., } } - fn __getstate__(&self) -> (Vec<(String, Vec)>, f64) { + fn __getstate__(&self) -> OneQubitGateSequenceState { (self.gates.clone(), self.global_phase) } - fn __setstate__(&mut self, state: (Vec<(String, Vec)>, f64)) { + fn __setstate__(&mut self, state: OneQubitGateSequenceState) { self.gates = state.0; self.global_phase = state.1; } @@ -93,7 +97,7 @@ impl OneQubitGateSequence { SliceOrInt::Slice(slc) => { let len = self.gates.len().try_into().unwrap(); let indices = slc.indices(len)?; - let mut out_vec: Vec<(String, Vec)> = Vec::new(); + let mut out_vec: Vec<(String, SmallVec<[f64; 3]>)> = Vec::new(); // Start and stop will always be positive the slice api converts // negatives to the index for example: // list(range(5))[-1:-3:-1] @@ -145,10 +149,10 @@ fn circuit_kak( let mut lam = lam; let mut theta = theta; let mut phi = phi; - let mut circuit: Vec<(String, Vec)> = Vec::with_capacity(3); + let mut circuit: Vec<(String, SmallVec<[f64; 3]>)> = Vec::with_capacity(3); let mut atol = match atol { Some(atol) => atol, - None => DEFAULT_ATOL, + None => ANGLE_ZERO_EPSILON, }; if !simplify { atol = -1.0; @@ -161,7 +165,7 @@ fn circuit_kak( // slippage coming from _mod_2pi injecting multiples of 2pi. lam = mod_2pi(lam, atol); if lam.abs() > atol { - circuit.push((String::from(k_gate), vec![lam])); + circuit.push((String::from(k_gate), smallvec![lam])); global_phase += lam / 2.; } return OneQubitGateSequence { @@ -182,13 +186,13 @@ fn circuit_kak( lam = mod_2pi(lam, atol); if lam.abs() > atol { global_phase += lam / 2.; - circuit.push((String::from(k_gate), vec![lam])); + circuit.push((String::from(k_gate), smallvec![lam])); } - circuit.push((String::from(a_gate), vec![theta])); + circuit.push((String::from(a_gate), smallvec![theta])); phi = mod_2pi(phi, atol); if phi.abs() > atol { global_phase += phi / 2.; - circuit.push((String::from(k_gate), vec![phi])); + circuit.push((String::from(k_gate), smallvec![phi])); } OneQubitGateSequence { gates: circuit, @@ -207,12 +211,12 @@ fn circuit_u3( let mut circuit = Vec::new(); let atol = match atol { Some(atol) => atol, - None => DEFAULT_ATOL, + None => ANGLE_ZERO_EPSILON, }; let phi = mod_2pi(phi, atol); let lam = mod_2pi(lam, atol); if !simplify || theta.abs() > atol || phi.abs() > atol || lam.abs() > atol { - circuit.push((String::from("u3"), vec![theta, phi, lam])); + circuit.push((String::from("u3"), smallvec![theta, phi, lam])); } OneQubitGateSequence { gates: circuit, @@ -231,7 +235,7 @@ fn circuit_u321( let mut circuit = Vec::new(); let mut atol = match atol { Some(atol) => atol, - None => DEFAULT_ATOL, + None => ANGLE_ZERO_EPSILON, }; if !simplify { atol = -1.0; @@ -239,17 +243,17 @@ fn circuit_u321( if theta.abs() < atol { let tot = mod_2pi(phi + lam, atol); if tot.abs() > atol { - circuit.push((String::from("u1"), vec![tot])); + circuit.push((String::from("u1"), smallvec![tot])); } } else if (theta - PI / 2.).abs() < atol { circuit.push(( String::from("u2"), - vec![mod_2pi(phi, atol), mod_2pi(lam, atol)], + smallvec![mod_2pi(phi, atol), mod_2pi(lam, atol)], )); } else { circuit.push(( String::from("u3"), - vec![theta, mod_2pi(phi, atol), mod_2pi(lam, atol)], + smallvec![theta, mod_2pi(phi, atol), mod_2pi(lam, atol)], )); } OneQubitGateSequence { @@ -269,7 +273,7 @@ fn circuit_u( let mut circuit = Vec::new(); let mut atol = match atol { Some(atol) => atol, - None => DEFAULT_ATOL, + None => ANGLE_ZERO_EPSILON, }; if !simplify { atol = -1.0; @@ -277,7 +281,7 @@ fn circuit_u( let phi = mod_2pi(phi, atol); let lam = mod_2pi(lam, atol); if theta.abs() > atol || phi.abs() > atol || lam.abs() > atol { - circuit.push((String::from("u"), vec![theta, phi, lam])); + circuit.push((String::from("u"), smallvec![theta, phi, lam])); } OneQubitGateSequence { gates: circuit, @@ -310,7 +314,7 @@ where }; let mut atol = match atol { Some(atol) => atol, - None => DEFAULT_ATOL, + None => ANGLE_ZERO_EPSILON, }; if !simplify { atol = -1.0; @@ -370,7 +374,7 @@ fn circuit_rr( let mut circuit = Vec::new(); let mut atol = match atol { Some(atol) => atol, - None => DEFAULT_ATOL, + None => ANGLE_ZERO_EPSILON, }; if !simplify { atol = -1.0; @@ -379,19 +383,22 @@ fn circuit_rr( if mod_2pi((phi + lam) / 2., atol).abs() < atol { // This can be expressed as a single R gate if theta.abs() > atol { - circuit.push((String::from("r"), vec![theta, mod_2pi(PI / 2. + phi, atol)])); + circuit.push(( + String::from("r"), + smallvec![theta, mod_2pi(PI / 2. + phi, atol)], + )); } } else { // General case: use two R gates if (theta - PI).abs() > atol { circuit.push(( String::from("r"), - vec![theta - PI, mod_2pi(PI / 2. - lam, atol)], + smallvec![theta - PI, mod_2pi(PI / 2. - lam, atol)], )); } circuit.push(( String::from("r"), - vec![PI, mod_2pi(0.5 * (phi - lam + PI), atol)], + smallvec![PI, mod_2pi(0.5 * (phi - lam + PI), atol)], )); } @@ -403,7 +410,7 @@ fn circuit_rr( #[pyfunction] pub fn generate_circuit( - target_basis: &str, + target_basis: &EulerBasis, theta: f64, phi: f64, lam: f64, @@ -412,17 +419,17 @@ pub fn generate_circuit( atol: Option, ) -> PyResult { let res = match target_basis { - "ZYZ" => circuit_kak(theta, phi, lam, phase, "rz", "ry", simplify, atol), - "ZXZ" => circuit_kak(theta, phi, lam, phase, "rz", "rx", simplify, atol), - "XZX" => circuit_kak(theta, phi, lam, phase, "rx", "rz", simplify, atol), - "XYX" => circuit_kak(theta, phi, lam, phase, "rx", "ry", simplify, atol), - "U3" => circuit_u3(theta, phi, lam, phase, simplify, atol), - "U321" => circuit_u321(theta, phi, lam, phase, simplify, atol), - "U" => circuit_u(theta, phi, lam, phase, simplify, atol), - "PSX" => { + EulerBasis::ZYZ => circuit_kak(theta, phi, lam, phase, "rz", "ry", simplify, atol), + EulerBasis::ZXZ => circuit_kak(theta, phi, lam, phase, "rz", "rx", simplify, atol), + EulerBasis::XZX => circuit_kak(theta, phi, lam, phase, "rx", "rz", simplify, atol), + EulerBasis::XYX => circuit_kak(theta, phi, lam, phase, "rx", "ry", simplify, atol), + EulerBasis::U3 => circuit_u3(theta, phi, lam, phase, simplify, atol), + EulerBasis::U321 => circuit_u321(theta, phi, lam, phase, simplify, atol), + EulerBasis::U => circuit_u(theta, phi, lam, phase, simplify, atol), + EulerBasis::PSX => { let mut inner_atol = match atol { Some(atol) => atol, - None => DEFAULT_ATOL, + None => ANGLE_ZERO_EPSILON, }; if !simplify { inner_atol = -1.0; @@ -430,11 +437,11 @@ pub fn generate_circuit( let fnz = |circuit: &mut OneQubitGateSequence, phi: f64| { let phi = mod_2pi(phi, inner_atol); if phi.abs() > inner_atol { - circuit.gates.push((String::from("p"), vec![phi])); + circuit.gates.push((String::from("p"), smallvec![phi])); } }; let fnx = |circuit: &mut OneQubitGateSequence| { - circuit.gates.push((String::from("sx"), Vec::new())); + circuit.gates.push((String::from("sx"), SmallVec::new())); }; circuit_psx_gen( @@ -449,10 +456,10 @@ pub fn generate_circuit( None::>, ) } - "ZSX" => { + EulerBasis::ZSX => { let mut inner_atol = match atol { Some(atol) => atol, - None => DEFAULT_ATOL, + None => ANGLE_ZERO_EPSILON, }; if !simplify { inner_atol = -1.0; @@ -460,12 +467,12 @@ pub fn generate_circuit( let fnz = |circuit: &mut OneQubitGateSequence, phi: f64| { let phi = mod_2pi(phi, inner_atol); if phi.abs() > inner_atol { - circuit.gates.push((String::from("rz"), vec![phi])); + circuit.gates.push((String::from("rz"), smallvec![phi])); circuit.global_phase += phi / 2.; } }; let fnx = |circuit: &mut OneQubitGateSequence| { - circuit.gates.push((String::from("sx"), Vec::new())); + circuit.gates.push((String::from("sx"), SmallVec::new())); }; circuit_psx_gen( theta, @@ -479,10 +486,10 @@ pub fn generate_circuit( None::>, ) } - "U1X" => { + EulerBasis::U1X => { let mut inner_atol = match atol { Some(atol) => atol, - None => DEFAULT_ATOL, + None => ANGLE_ZERO_EPSILON, }; if !simplify { inner_atol = -1.0; @@ -490,12 +497,12 @@ pub fn generate_circuit( let fnz = |circuit: &mut OneQubitGateSequence, phi: f64| { let phi = mod_2pi(phi, inner_atol); if phi.abs() > inner_atol { - circuit.gates.push((String::from("u1"), vec![phi])); + circuit.gates.push((String::from("u1"), smallvec![phi])); } }; let fnx = |circuit: &mut OneQubitGateSequence| { circuit.global_phase += PI / 4.; - circuit.gates.push((String::from("rx"), vec![PI / 2.])); + circuit.gates.push((String::from("rx"), smallvec![PI / 2.])); }; circuit_psx_gen( theta, @@ -509,10 +516,10 @@ pub fn generate_circuit( None::>, ) } - "ZSXX" => { + EulerBasis::ZSXX => { let mut inner_atol = match atol { Some(atol) => atol, - None => DEFAULT_ATOL, + None => ANGLE_ZERO_EPSILON, }; if !simplify { inner_atol = -1.0; @@ -520,15 +527,15 @@ pub fn generate_circuit( let fnz = |circuit: &mut OneQubitGateSequence, phi: f64| { let phi = mod_2pi(phi, inner_atol); if phi.abs() > inner_atol { - circuit.gates.push((String::from("rz"), vec![phi])); + circuit.gates.push((String::from("rz"), smallvec![phi])); circuit.global_phase += phi / 2.; } }; let fnx = |circuit: &mut OneQubitGateSequence| { - circuit.gates.push((String::from("sx"), Vec::new())); + circuit.gates.push((String::from("sx"), SmallVec::new())); }; let fnxpi = |circuit: &mut OneQubitGateSequence| { - circuit.gates.push((String::from("x"), Vec::new())); + circuit.gates.push((String::from("x"), SmallVec::new())); }; circuit_psx_gen( theta, @@ -542,32 +549,88 @@ pub fn generate_circuit( Some(fnxpi), ) } - "RR" => circuit_rr(theta, phi, lam, phase, simplify, atol), - other => { - return Err(PyTypeError::new_err(format!( - "Invalid target basis: {other}" - ))) - } + EulerBasis::RR => circuit_rr(theta, phi, lam, phase, simplify, atol), }; Ok(res) } +#[derive(Clone, Debug, Copy)] +#[pyclass(module = "qiskit._accelerate.euler_one_qubit_decomposer")] +pub enum EulerBasis { + U321, + U3, + U, + PSX, + ZSX, + ZSXX, + U1X, + RR, + ZYZ, + ZXZ, + XYX, + XZX, +} + +#[pymethods] +impl EulerBasis { + #![allow(clippy::wrong_self_convention)] + pub fn to_str(&self) -> String { + match self { + Self::U321 => "U321".to_string(), + Self::U3 => "U3".to_string(), + Self::U => "U".to_string(), + Self::PSX => "PSX".to_string(), + Self::ZSX => "ZSX".to_string(), + Self::ZSXX => "ZSXX".to_string(), + Self::U1X => "U1X".to_string(), + Self::RR => "RR".to_string(), + Self::ZYZ => "ZYZ".to_string(), + Self::ZXZ => "ZXZ".to_string(), + Self::XYX => "XYX".to_string(), + Self::XZX => "XZX".to_string(), + } + } + + #[staticmethod] + pub fn from_string(input: &str) -> PyResult { + let res = match input { + "U321" => EulerBasis::U321, + "U3" => EulerBasis::U3, + "U" => EulerBasis::U, + "PSX" => EulerBasis::PSX, + "ZSX" => EulerBasis::ZSX, + "ZSXX" => EulerBasis::ZSXX, + "U1X" => EulerBasis::U1X, + "RR" => EulerBasis::RR, + "ZYZ" => EulerBasis::ZYZ, + "ZXZ" => EulerBasis::ZXZ, + "XYX" => EulerBasis::XYX, + "XZX" => EulerBasis::XZX, + basis => { + return Err(PyTypeError::new_err(format!( + "Invalid target basis {basis}" + ))); + } + }; + Ok(res) + } +} + #[inline] -fn angles_from_unitary(unitary: ArrayView2, target_basis: &str) -> [f64; 4] { +pub fn angles_from_unitary(unitary: ArrayView2, target_basis: EulerBasis) -> [f64; 4] { match target_basis { - "U321" => params_u3_inner(unitary), - "U3" => params_u3_inner(unitary), - "U" => params_u3_inner(unitary), - "PSX" => params_u1x_inner(unitary), - "ZSX" => params_u1x_inner(unitary), - "ZSXX" => params_u1x_inner(unitary), - "U1X" => params_u1x_inner(unitary), - "RR" => params_zyz_inner(unitary), - "ZYZ" => params_zyz_inner(unitary), - "ZXZ" => params_zxz_inner(unitary), - "XYX" => params_xyx_inner(unitary), - "XZX" => params_xzx_inner(unitary), - &_ => unreachable!(), + EulerBasis::U321 => params_u3_inner(unitary), + EulerBasis::U3 => params_u3_inner(unitary), + EulerBasis::U => params_u3_inner(unitary), + EulerBasis::PSX => params_u1x_inner(unitary), + EulerBasis::ZSX => params_u1x_inner(unitary), + EulerBasis::ZSXX => params_u1x_inner(unitary), + EulerBasis::U1X => params_u1x_inner(unitary), + EulerBasis::RR => params_zyz_inner(unitary), + EulerBasis::ZYZ => params_zyz_inner(unitary), + EulerBasis::ZXZ => params_zxz_inner(unitary), + EulerBasis::XYX => params_xyx_inner(unitary), + EulerBasis::XZX => params_xzx_inner(unitary), } } @@ -592,7 +655,7 @@ fn compare_error_fn( } fn compute_error( - gates: &[(String, Vec)], + gates: &[(String, SmallVec<[f64; 3]>)], error_map: Option<&OneQubitGateErrorMap>, qubit: usize, ) -> (f64, usize) { @@ -620,7 +683,7 @@ pub fn compute_error_one_qubit_sequence( #[pyfunction] pub fn compute_error_list( - circuit: Vec<(String, Vec)>, + circuit: Vec<(String, SmallVec<[f64; 3]>)>, qubit: usize, error_map: Option<&OneQubitGateErrorMap>, ) -> (f64, usize) { @@ -637,33 +700,46 @@ pub fn unitary_to_gate_sequence( simplify: bool, atol: Option, ) -> PyResult> { - const VALID_BASES: [&str; 12] = [ - "U321", "U3", "U", "PSX", "ZSX", "ZSXX", "U1X", "RR", "ZYZ", "ZXZ", "XYX", "XZX", - ]; - for basis in &target_basis_list { - if !VALID_BASES.contains(basis) { - return Err(PyTypeError::new_err(format!( - "Invalid target basis {basis}" - ))); - } + let mut target_basis_vec: Vec = Vec::with_capacity(target_basis_list.len()); + for basis in target_basis_list { + let basis_enum = EulerBasis::from_string(basis)?; + target_basis_vec.push(basis_enum) } let unitary_mat = unitary.as_array(); - let best_result = target_basis_list + Ok(unitary_to_gate_sequence_inner( + unitary_mat, + &target_basis_vec, + qubit, + error_map, + simplify, + atol, + )) +} + +#[inline] +pub fn unitary_to_gate_sequence_inner( + unitary_mat: ArrayView2, + target_basis_list: &[EulerBasis], + qubit: usize, + error_map: Option<&OneQubitGateErrorMap>, + simplify: bool, + atol: Option, +) -> Option { + target_basis_list .iter() .map(|target_basis| { - let [theta, phi, lam, phase] = angles_from_unitary(unitary_mat, target_basis); + let [theta, phi, lam, phase] = angles_from_unitary(unitary_mat, *target_basis); generate_circuit(target_basis, theta, phi, lam, phase, simplify, atol).unwrap() }) .min_by(|a, b| { let error_a = compare_error_fn(a, &error_map, qubit); let error_b = compare_error_fn(b, &error_map, qubit); error_a.partial_cmp(&error_b).unwrap_or(Ordering::Equal) - }); - Ok(best_result) + }) } #[inline] -fn det_one_qubit(mat: ArrayView2) -> Complex64 { +pub fn det_one_qubit(mat: ArrayView2) -> Complex64 { mat[[0, 0]] * mat[[1, 1]] - mat[[0, 1]] * mat[[1, 0]] } diff --git a/crates/accelerate/src/two_qubit_decompose.rs b/crates/accelerate/src/two_qubit_decompose.rs index d9f451a43682..cb3c1308829e 100644 --- a/crates/accelerate/src/two_qubit_decompose.rs +++ b/crates/accelerate/src/two_qubit_decompose.rs @@ -18,59 +18,127 @@ // of real components and one of imaginary components. // In order to avoid copying we want to use `MatRef` or `MatMut`. -use num_complex::Complex; +use approx::abs_diff_eq; +use num_complex::{Complex, Complex64, ComplexFloat}; +use pyo3::exceptions::PyIndexError; +use pyo3::import_exception; use pyo3::prelude::*; use pyo3::wrap_pyfunction; use pyo3::Python; +use smallvec::{smallvec, SmallVec}; use std::f64::consts::PI; use faer::IntoFaerComplex; -use faer::{mat, prelude::*, scale, Mat, MatRef}; +use faer::IntoNdarray; +use faer::IntoNdarrayComplex; +use faer::Side::Lower; +use faer::{prelude::*, scale, Mat, MatRef}; use faer_core::{c64, ComplexField}; +use ndarray::linalg::kron; +use ndarray::prelude::*; +use ndarray::Zip; use numpy::PyReadonlyArray2; +use numpy::ToPyArray; +use crate::euler_one_qubit_decomposer::{ + angles_from_unitary, det_one_qubit, unitary_to_gate_sequence_inner, EulerBasis, + ANGLE_ZERO_EPSILON, +}; use crate::utils; +use rand::prelude::*; +use rand_distr::StandardNormal; +use rand_pcg::Pcg64Mcg; + const PI2: f64 = PI / 2.0; const PI4: f64 = PI / 4.0; const PI32: f64 = 3.0 * PI2; -const C0: c64 = c64 { re: 0.0, im: 0.0 }; +const TWO_PI: f64 = 2.0 * PI; + const C1: c64 = c64 { re: 1.0, im: 0.0 }; -const C1_NEG: c64 = c64 { re: -1.0, im: 0.0 }; -const C1IM: c64 = c64 { re: 0.0, im: 1.0 }; -const C1_NEG_IM: c64 = c64 { re: 0.0, im: -1.0 }; -const C0_5: c64 = c64 { re: 0.5, im: 0.0 }; -const C0_5_NEG: c64 = c64 { re: -0.5, im: 0.0 }; -const C0_5_IM: c64 = c64 { re: 0.0, im: 0.5 }; -const C0_5_IM_NEG: c64 = c64 { re: 0.0, im: -0.5 }; - -const B_NON_NORMALIZED: [c64; 16] = [ - C1, C1IM, C0, C0, C0, C0, C1IM, C1, C0, C0, C1IM, C1_NEG, C1, C1_NEG_IM, C0, C0, + +const ONE_QUBIT_IDENTITY: [[Complex64; 2]; 2] = [ + [Complex64::new(1., 0.), Complex64::new(0., 0.)], + [Complex64::new(0., 0.), Complex64::new(1., 0.)], ]; -const B_NON_NORMALIZED_DAGGER: [c64; 16] = [ - C0_5, - C0, - C0, - C0_5, - C0_5_IM_NEG, - C0, - C0, - C0_5_IM, - C0, - C0_5_IM_NEG, - C0_5_IM_NEG, - C0, - C0, - C0_5, - C0_5_NEG, - C0, +const B_NON_NORMALIZED: [[Complex64; 4]; 4] = [ + [ + Complex64::new(1.0, 0.), + Complex64::new(0., 1.), + Complex64::new(0., 0.), + Complex64::new(0., 0.), + ], + [ + Complex64::new(0., 0.), + Complex64::new(0., 0.), + Complex64::new(0., 1.), + Complex64::new(1.0, 0.0), + ], + [ + Complex64::new(0., 0.), + Complex64::new(0., 0.), + Complex64::new(0., 1.), + Complex64::new(-1., 0.), + ], + [ + Complex64::new(1., 0.), + Complex64::new(0., -1.), + Complex64::new(0., 0.), + Complex64::new(0., 0.), + ], ]; -fn transform_from_magic_basis(unitary: Mat) -> Mat { - let _b_nonnormalized = mat::from_row_major_slice::(&B_NON_NORMALIZED, 4, 4); - let _b_nonnormalized_dagger = mat::from_row_major_slice::(&B_NON_NORMALIZED_DAGGER, 4, 4); - _b_nonnormalized_dagger * unitary * _b_nonnormalized +const B_NON_NORMALIZED_DAGGER: [[Complex64; 4]; 4] = [ + [ + Complex64::new(0.5, 0.), + Complex64::new(0., 0.), + Complex64::new(0., 0.), + Complex64::new(0.5, 0.0), + ], + [ + Complex64::new(0., -0.5), + Complex64::new(0., 0.), + Complex64::new(0., 0.), + Complex64::new(0., 0.5), + ], + [ + Complex64::new(0., 0.), + Complex64::new(0., -0.5), + Complex64::new(0., -0.5), + Complex64::new(0., 0.), + ], + [ + Complex64::new(0., 0.), + Complex64::new(0.5, 0.), + Complex64::new(-0.5, 0.), + Complex64::new(0., 0.), + ], +]; + +enum MagicBasisTransform { + Into, + OutOf, +} + +fn magic_basis_transform( + unitary: ArrayView2, + direction: MagicBasisTransform, +) -> Array2 { + let _b_nonnormalized = aview2(&B_NON_NORMALIZED); + let _b_nonnormalized_dagger = aview2(&B_NON_NORMALIZED_DAGGER); + match direction { + MagicBasisTransform::OutOf => _b_nonnormalized_dagger.dot(&unitary).dot(&_b_nonnormalized), + MagicBasisTransform::Into => _b_nonnormalized.dot(&unitary).dot(&_b_nonnormalized_dagger), + } +} + +fn transform_from_magic_basis(u: Mat) -> Mat { + let unitary: ArrayView2 = u.as_ref().into_ndarray_complex(); + magic_basis_transform(unitary, MagicBasisTransform::OutOf) + .view() + .into_faer_complex() + .to_owned() } // faer::c64 and num_complex::Complex are both structs @@ -98,6 +166,57 @@ impl Arg for c64 { } } +pub trait TraceToFidelity { + /// Average gate fidelity is :math:`Fbar = (d + |Tr (Utarget \\cdot U^dag)|^2) / d(d+1)` + /// M. Horodecki, P. Horodecki and R. Horodecki, PRA 60, 1888 (1999) + fn trace_to_fid(self) -> f64; +} + +impl TraceToFidelity for Complex64 { + fn trace_to_fid(self) -> f64 { + (4.0 + self.abs().powi(2)) / 20.0 + } +} + +impl TraceToFidelity for c64 { + fn trace_to_fid(self) -> f64 { + (4.0 + self.faer_abs2()) / 20.0 + } +} + +fn decompose_two_qubit_product_gate( + special_unitary: ArrayView2, +) -> PyResult<(Array2, Array2, f64)> { + let mut r: Array2 = special_unitary.slice(s![..2, ..2]).to_owned(); + let mut det_r = det_one_qubit(r.view()); + if det_r.abs() < 0.1 { + r = special_unitary.slice(s![2.., ..2]).to_owned(); + det_r = det_one_qubit(r.view()); + } + if det_r.abs() < 0.1 { + import_exception!(qiskit, QiskitError); + return Err(QiskitError::new_err( + "decompose_two_qubit_product_gate: unable to decompose: detR < 0.1", + )); + } + r.mapv_inplace(|x| x / det_r.sqrt()); + let r_t_conj: Array2 = r.t().mapv(|x| x.conj()); + let eye = aview2(&ONE_QUBIT_IDENTITY); + let mut temp = kron(&eye, &r_t_conj); + temp = special_unitary.dot(&temp); + let mut l = temp.slice(s![..;2, ..;2]).to_owned(); + let det_l = det_one_qubit(l.view()); + if det_l.abs() < 0.9 { + import_exception!(qiskit, QiskitError); + return Err(QiskitError::new_err( + "decompose_two_qubit_product_gate: unable to decompose: detL < 0.9", + )); + } + l.mapv_inplace(|x| x / det_l.sqrt()); + let phase = det_l.arg() / 2.; + Ok((l, r, phase)) +} + fn __weyl_coordinates(unitary: MatRef) -> [f64; 3] { let uscaled = scale(C1 / unitary.determinant().powf(0.25)) * unitary; let uup = transform_from_magic_basis(uscaled); @@ -180,20 +299,926 @@ fn __num_basis_gates(basis_b: f64, basis_fidelity: f64, unitary: MatRef) -> traces .into_iter() .enumerate() - .map(|(idx, trace)| (idx, trace_to_fid(&trace) * basis_fidelity.powi(idx as i32))) + .map(|(idx, trace)| (idx, trace.trace_to_fid() * basis_fidelity.powi(idx as i32))) .min_by(|(_idx1, fid1), (_idx2, fid2)| fid2.partial_cmp(fid1).unwrap()) .unwrap() .0 } -/// Average gate fidelity is :math:`Fbar = (d + |Tr (Utarget \\cdot U^dag)|^2) / d(d+1)` -/// M. Horodecki, P. Horodecki and R. Horodecki, PRA 60, 1888 (1999) -fn trace_to_fid(trace: &c64) -> f64 { - (4.0 + trace.faer_abs2()) / 20.0 +/// A good approximation to the best value x to get the minimum +/// trace distance for :math:`U_d(x, x, x)` from :math:`U_d(a, b, c)`. +fn closest_partial_swap(a: f64, b: f64, c: f64) -> f64 { + let m = (a + b + c) / 3.; + let [am, bm, cm] = [a - m, b - m, c - m]; + let [ab, bc, ca] = [a - b, b - c, c - a]; + m + am * bm * cm * (6. + ab * ab + bc * bc + ca * ca) / 18. +} + +fn rx_matrix(theta: f64) -> Array2 { + let half_theta = theta / 2.; + let cos = Complex64::new(half_theta.cos(), 0.); + let isin = Complex64::new(0., -half_theta.sin()); + array![[cos, isin], [isin, cos]] +} + +fn ry_matrix(theta: f64) -> Array2 { + let half_theta = theta / 2.; + let cos = Complex64::new(half_theta.cos(), 0.); + let sin = Complex64::new(half_theta.sin(), 0.); + array![[cos, -sin], [sin, cos]] +} + +fn rz_matrix(theta: f64) -> Array2 { + let ilam2 = Complex64::new(0., 0.5 * theta); + array![ + [(-ilam2).exp(), Complex64::new(0., 0.)], + [Complex64::new(0., 0.), ilam2.exp()] + ] +} + +const DEFAULT_FIDELITY: f64 = 1.0 - 1.0e-9; +const C1_IM: Complex64 = Complex64::new(0.0, 1.0); + +#[derive(Clone, Debug, Copy)] +#[pyclass(module = "qiskit._accelerate.two_qubit_decompose")] +enum Specialization { + General, + IdEquiv, + SWAPEquiv, + PartialSWAPEquiv, + PartialSWAPFlipEquiv, + ControlledEquiv, + MirrorControlledEquiv, + // These next 3 gates use the definition of fSim from eq (1) in: + // https://arxiv.org/pdf/2001.08343.pdf + #[allow(non_camel_case_types)] + fSimaabEquiv, + #[allow(non_camel_case_types)] + fSimabbEquiv, + #[allow(non_camel_case_types)] + fSimabmbEquiv, +} + +impl Specialization { + fn to_u8(self) -> u8 { + match self { + Specialization::General => 0, + Specialization::IdEquiv => 1, + Specialization::SWAPEquiv => 2, + Specialization::PartialSWAPEquiv => 3, + Specialization::PartialSWAPFlipEquiv => 4, + Specialization::ControlledEquiv => 5, + Specialization::MirrorControlledEquiv => 6, + Specialization::fSimaabEquiv => 7, + Specialization::fSimabbEquiv => 8, + Specialization::fSimabmbEquiv => 9, + } + } + + fn from_u8(val: u8) -> Self { + match val { + 0 => Specialization::General, + 1 => Specialization::IdEquiv, + 2 => Specialization::SWAPEquiv, + 3 => Specialization::PartialSWAPEquiv, + 4 => Specialization::PartialSWAPFlipEquiv, + 5 => Specialization::ControlledEquiv, + 6 => Specialization::MirrorControlledEquiv, + 7 => Specialization::fSimaabEquiv, + 8 => Specialization::fSimabbEquiv, + 9 => Specialization::fSimabmbEquiv, + _ => unreachable!("Invalid specialization value"), + } + } +} + +#[derive(Clone, Debug)] +#[allow(non_snake_case)] +#[pyclass(module = "qiskit._accelerate.two_qubit_decompose", subclass)] +pub struct TwoQubitWeylDecomposition { + #[pyo3(get)] + a: f64, + #[pyo3(get)] + b: f64, + #[pyo3(get)] + c: f64, + #[pyo3(get)] + global_phase: f64, + K1l: Array2, + K2l: Array2, + K1r: Array2, + K2r: Array2, + #[pyo3(get)] + specialization: Specialization, + default_euler_basis: EulerBasis, + #[pyo3(get)] + requested_fidelity: Option, + #[pyo3(get)] + calculated_fidelity: f64, + unitary_matrix: Array2, +} + +impl TwoQubitWeylDecomposition { + fn weyl_gate( + &self, + simplify: bool, + sequence: &mut TwoQubitSequenceVec, + atol: f64, + global_phase: &mut f64, + ) { + match self.specialization { + Specialization::MirrorControlledEquiv => { + sequence.push(("swap".to_string(), SmallVec::new(), smallvec![0, 1])); + sequence.push(( + "rzz".to_string(), + smallvec![(PI4 - self.c) * 2.], + smallvec![0, 1], + )); + *global_phase += PI4 + } + Specialization::SWAPEquiv => { + sequence.push(("swap".to_string(), SmallVec::new(), smallvec![0, 1])); + *global_phase -= 3. * PI / 4. + } + _ => { + if !simplify || self.a.abs() > atol { + sequence.push(("rxx".to_string(), smallvec![-self.a * 2.], smallvec![0, 1])); + } + if !simplify || self.b.abs() > atol { + sequence.push(("ryy".to_string(), smallvec![-self.b * 2.], smallvec![0, 1])); + } + if !simplify || self.c.abs() > atol { + sequence.push(("rzz".to_string(), smallvec![-self.c * 2.], smallvec![0, 1])); + } + } + } + } +} + +const IPZ: [[Complex64; 2]; 2] = [ + [C1_IM, Complex64::new(0., 0.)], + [Complex64::new(0., 0.), Complex64::new(0., -1.)], +]; +const IPY: [[Complex64; 2]; 2] = [ + [Complex64::new(0., 0.), Complex64::new(1., 0.)], + [Complex64::new(-1., 0.), Complex64::new(0., 0.)], +]; +const IPX: [[Complex64; 2]; 2] = [ + [Complex64::new(0., 0.), C1_IM], + [C1_IM, Complex64::new(0., 0.)], +]; + +#[pymethods] +impl TwoQubitWeylDecomposition { + fn __getstate__(&self, py: Python) -> ([f64; 5], [PyObject; 5], u8, String, Option) { + let specialization = self.specialization.to_u8(); + ( + [ + self.a, + self.b, + self.c, + self.global_phase, + self.calculated_fidelity, + ], + [ + self.K1l.to_pyarray(py).into(), + self.K1r.to_pyarray(py).into(), + self.K2l.to_pyarray(py).into(), + self.K2r.to_pyarray(py).into(), + self.unitary_matrix.to_pyarray(py).into(), + ], + specialization, + self.default_euler_basis.to_str(), + self.requested_fidelity, + ) + } + + fn __setstate__( + &mut self, + state: ( + [f64; 5], + [PyReadonlyArray2; 5], + u8, + String, + Option, + ), + ) { + self.a = state.0[0]; + self.b = state.0[1]; + self.c = state.0[2]; + self.global_phase = state.0[3]; + self.calculated_fidelity = state.0[4]; + self.K1l = state.1[0].as_array().to_owned(); + self.K1r = state.1[1].as_array().to_owned(); + self.K2l = state.1[2].as_array().to_owned(); + self.K2r = state.1[3].as_array().to_owned(); + self.unitary_matrix = state.1[4].as_array().to_owned(); + self.default_euler_basis = EulerBasis::from_string(&state.3).unwrap(); + self.requested_fidelity = state.4; + self.specialization = Specialization::from_u8(state.2); + } + + fn __getnewargs__(&self, py: Python) -> (PyObject, Option, Option, bool) { + ( + self.unitary_matrix.to_pyarray(py).into(), + self.requested_fidelity, + None, + true, + ) + } + + #[new] + #[pyo3(signature=(unitary_matrix, fidelity=DEFAULT_FIDELITY, _specialization=None, _pickle_context=false))] + fn new( + unitary_matrix: PyReadonlyArray2, + fidelity: Option, + _specialization: Option, + _pickle_context: bool, + ) -> PyResult { + // If we're in a pickle context just make the closest to an empty + // object as we can with minimal allocations and effort. All the + // data will be filled in during deserialization from __setstate__. + if _pickle_context { + return Ok(TwoQubitWeylDecomposition { + a: 0., + b: 0., + c: 0., + global_phase: 0., + K1l: Array2::zeros((0, 0)), + K2l: Array2::zeros((0, 0)), + K1r: Array2::zeros((0, 0)), + K2r: Array2::zeros((0, 0)), + specialization: Specialization::General, + default_euler_basis: EulerBasis::ZYZ, + requested_fidelity: fidelity, + calculated_fidelity: 0., + unitary_matrix: Array2::zeros((0, 0)), + }); + } + let ipz: ArrayView2 = aview2(&IPZ); + let ipy: ArrayView2 = aview2(&IPY); + let ipx: ArrayView2 = aview2(&IPX); + + let mut u = unitary_matrix.as_array().to_owned(); + let unitary_matrix = unitary_matrix.as_array().to_owned(); + let det_u = u.view().into_faer_complex().determinant().to_num_complex(); + let det_pow = det_u.powf(-0.25); + u.mapv_inplace(|x| x * det_pow); + let mut global_phase = det_u.arg() / 4.; + let u_p = magic_basis_transform(u.view(), MagicBasisTransform::OutOf); + let m2 = u_p.t().dot(&u_p); + let default_euler_basis = EulerBasis::ZYZ; + + // M2 is a symmetric complex matrix. We need to decompose it as M2 = P D P^T where + // P ∈ SO(4), D is diagonal with unit-magnitude elements. + // + // We can't use raw `eig` directly because it isn't guaranteed to give us real or othogonal + // eigenvectors. Instead, since `M2` is complex-symmetric, + // M2 = A + iB + // for real-symmetric `A` and `B`, and as + // M2^+ @ M2 = A^2 + B^2 + i [A, B] = 1 + // we must have `A` and `B` commute, and consequently they are simultaneously diagonalizable. + // Mixing them together _should_ account for any degeneracy problems, but it's not + // guaranteed, so we repeat it a little bit. The fixed seed is to make failures + // deterministic; the value is not important. + let mut state = Pcg64Mcg::seed_from_u64(2023); + let mut found = false; + let mut d: Array1 = Array1::zeros(0); + let mut p: Array2 = Array2::zeros((0, 0)); + for i in 0..100 { + let rand_a: f64; + let rand_b: f64; + // For debugging the algorithm use the same RNG values from the + // previous Python implementation for the first random trial. + // In most cases this loop only executes a single iteration and + // using the same rng values rules out possible RNG differences + // as the root cause of a test failure + if i == 0 { + rand_a = 1.2602066112249388; + rand_b = 0.22317849046722027; + } else { + rand_a = state.sample(StandardNormal); + rand_b = state.sample(StandardNormal); + } + let m2_real = m2.mapv(|val| rand_a * val.re + rand_b * val.im); + let p_inner = m2_real + .view() + .into_faer() + .selfadjoint_eigendecomposition(Lower) + .u() + .into_ndarray() + .mapv(Complex64::from); + let d_inner = p_inner.t().dot(&m2).dot(&p_inner).diag().to_owned(); + let mut diag_d: Array2 = Array2::zeros((4, 4)); + diag_d + .diag_mut() + .iter_mut() + .enumerate() + .for_each(|(index, x)| *x = d_inner[index]); + + let compare = p_inner.dot(&diag_d).dot(&p_inner.t()); + found = abs_diff_eq!(compare.view(), m2, epsilon = 1.0e-13); + if found { + p = p_inner; + d = d_inner; + break; + } + } + if !found { + import_exception!(qiskit, QiskitError); + return Err(QiskitError::new_err(format!( + "TwoQubitWeylDecomposition: failed to diagonalize M2. Please report this at https://github.com/Qiskit/qiskit-terra/issues/4159. Input: {:?}", unitary_matrix + ))); + } + let mut d = -d.map(|x| x.arg() / 2.); + d[3] = -d[0] - d[1] - d[2]; + let mut cs: SmallVec<[f64; 3]> = (0..3) + .map(|i| ((d[i] + d[3]) / 2.0).rem_euclid(TWO_PI)) + .collect(); + let cstemp: SmallVec<[f64; 3]> = cs + .iter() + .map(|x| x.rem_euclid(PI2)) + .map(|x| x.min(PI2 - x)) + .collect(); + let mut order = utils::arg_sort(&cstemp); + (order[0], order[1], order[2]) = (order[1], order[2], order[0]); + (cs[0], cs[1], cs[2]) = (cs[order[0]], cs[order[1]], cs[order[2]]); + (d[0], d[1], d[2]) = (d[order[0]], d[order[1]], d[order[2]]); + let mut p_orig = p.clone(); + for (i, item) in order.iter().enumerate().take(3) { + let slice_a = p.slice_mut(s![.., i]); + let slice_b = p_orig.slice_mut(s![.., *item]); + Zip::from(slice_a).and(slice_b).for_each(::std::mem::swap); + } + if p.view().into_faer_complex().determinant().re < 0. { + p.slice_mut(s![.., -1]).mapv_inplace(|x| -x); + } + let mut temp: Array2 = Array2::zeros((4, 4)); + temp.diag_mut() + .iter_mut() + .enumerate() + .for_each(|(index, x)| *x = (C1_IM * d[index]).exp()); + let k1 = magic_basis_transform(u_p.dot(&p).dot(&temp).view(), MagicBasisTransform::Into); + let k2 = magic_basis_transform(p.t(), MagicBasisTransform::Into); + + #[allow(non_snake_case)] + let (mut K1l, mut K1r, phase_l) = decompose_two_qubit_product_gate(k1.view())?; + #[allow(non_snake_case)] + let (K2l, mut K2r, phase_r) = decompose_two_qubit_product_gate(k2.view())?; + global_phase += phase_l + phase_r; + + // Flip into Weyl chamber + if cs[0] > PI2 { + cs[0] -= PI32; + K1l = K1l.dot(&ipy); + K1r = K1r.dot(&ipy); + global_phase += PI2; + } + if cs[1] > PI2 { + cs[1] -= PI32; + K1l = K1l.dot(&ipx); + K1r = K1r.dot(&ipx); + global_phase += PI2; + } + let mut conjs = 0; + if cs[0] > PI4 { + cs[0] = PI2 - cs[0]; + K1l = K1l.dot(&ipy); + K2r = ipy.dot(&K2r); + conjs += 1; + global_phase -= PI2; + } + if cs[1] > PI4 { + cs[1] = PI2 - cs[1]; + K1l = K1l.dot(&ipx); + K2r = ipx.dot(&K2r); + conjs += 1; + global_phase += PI2; + if conjs == 1 { + global_phase -= PI; + } + } + if cs[2] > PI2 { + cs[2] -= PI32; + K1l = K1l.dot(&ipz); + K1r = K1r.dot(&ipz); + global_phase += PI2; + if conjs == 1 { + global_phase -= PI; + } + } + if conjs == 1 { + cs[2] = PI2 - cs[2]; + K1l = K1l.dot(&ipz); + K2r = ipz.dot(&K2r); + global_phase += PI2; + } + if cs[2] > PI4 { + cs[2] -= PI2; + K1l = K1l.dot(&ipz); + K1r = K1r.dot(&ipz); + global_phase -= PI2; + } + let [a, b, c] = [cs[1], cs[0], cs[2]]; + let is_close = |ap: f64, bp: f64, cp: f64| -> bool { + let [da, db, dc] = [a - ap, b - bp, c - cp]; + let tr = 4. + * Complex64::new( + da.cos() * db.cos() * dc.cos(), + da.sin() * db.sin() * dc.sin(), + ); + match fidelity { + Some(fid) => tr.trace_to_fid() >= fid, + // Set to false here to default to general specialization in the absence of a + // fidelity and provided specialization. + None => false, + } + }; + + let closest_abc = closest_partial_swap(a, b, c); + let closest_ab_minus_c = closest_partial_swap(a, b, -c); + let mut flipped_from_original = false; + let specialization = match _specialization { + Some(specialization) => specialization, + None => { + if is_close(0., 0., 0.) { + Specialization::IdEquiv + } else if is_close(PI4, PI4, PI4) || is_close(PI4, PI4, -PI4) { + Specialization::SWAPEquiv + } else if is_close(closest_abc, closest_abc, closest_abc) { + Specialization::PartialSWAPEquiv + } else if is_close(closest_ab_minus_c, closest_ab_minus_c, -closest_ab_minus_c) { + Specialization::PartialSWAPFlipEquiv + } else if is_close(a, 0., 0.) { + Specialization::ControlledEquiv + } else if is_close(PI4, PI4, c) { + Specialization::MirrorControlledEquiv + } else if is_close((a + b) / 2., (a + b) / 2., c) { + Specialization::fSimaabEquiv + } else if is_close(a, (b + c) / 2., (b + c) / 2.) { + Specialization::fSimabbEquiv + } else if is_close(a, (b - c) / 2., (c - b) / 2.) { + Specialization::fSimabmbEquiv + } else { + Specialization::General + } + } + }; + let general = TwoQubitWeylDecomposition { + a, + b, + c, + global_phase, + K1l, + K1r, + K2l, + K2r, + specialization: Specialization::General, + default_euler_basis, + requested_fidelity: fidelity, + calculated_fidelity: -1.0, + unitary_matrix, + }; + let mut specialized: TwoQubitWeylDecomposition = match specialization { + // :math:`U \sim U_d(0,0,0) \sim Id` + // + // This gate binds 0 parameters, we make it canonical by setting + // :math:`K2_l = Id` , :math:`K2_r = Id`. + Specialization::IdEquiv => TwoQubitWeylDecomposition { + specialization, + a: 0., + b: 0., + c: 0., + K1l: general.K1l.dot(&general.K2l), + K1r: general.K1r.dot(&general.K2r), + K2l: Array2::eye(2), + K2r: Array2::eye(2), + ..general + }, + // :math:`U \sim U_d(\pi/4, \pi/4, \pi/4) \sim U(\pi/4, \pi/4, -\pi/4) \sim \text{SWAP}` + // + // This gate binds 0 parameters, we make it canonical by setting + // :math:`K2_l = Id` , :math:`K2_r = Id`. + Specialization::SWAPEquiv => { + if c > 0. { + TwoQubitWeylDecomposition { + specialization, + a: PI4, + b: PI4, + c: PI4, + K1l: general.K1l.dot(&general.K2r), + K1r: general.K1r.dot(&general.K2l), + K2l: Array2::eye(2), + K2r: Array2::eye(2), + ..general + } + } else { + flipped_from_original = true; + TwoQubitWeylDecomposition { + specialization, + a: PI4, + b: PI4, + c: PI4, + global_phase: global_phase + PI2, + K1l: general.K1l.dot(&ipz).dot(&general.K2r), + K1r: general.K1r.dot(&ipz).dot(&general.K2l), + K2l: Array2::eye(2), + K2r: Array2::eye(2), + ..general + } + } + } + // :math:`U \sim U_d(\alpha\pi/4, \alpha\pi/4, \alpha\pi/4) \sim \text{SWAP}^\alpha` + // + // This gate binds 3 parameters, we make it canonical by setting: + // + // :math:`K2_l = Id`. + Specialization::PartialSWAPEquiv => { + let closest = closest_partial_swap(a, b, c); + let mut k2l_dag = general.K2l.t().to_owned(); + k2l_dag.view_mut().mapv_inplace(|x| x.conj()); + TwoQubitWeylDecomposition { + specialization, + a: closest, + b: closest, + c: closest, + K1l: general.K1l.dot(&general.K2l), + K1r: general.K1r.dot(&general.K2l), + K2r: k2l_dag.dot(&general.K2r), + K2l: Array2::eye(2), + ..general + } + } + // :math:`U \sim U_d(\alpha\pi/4, \alpha\pi/4, -\alpha\pi/4) \sim \text{SWAP}^\alpha` + // + // (a non-equivalent root of SWAP from the TwoQubitWeylPartialSWAPEquiv + // similar to how :math:`x = (\pm \sqrt(x))^2`) + // + // This gate binds 3 parameters, we make it canonical by setting: + // + // :math:`K2_l = Id` + Specialization::PartialSWAPFlipEquiv => { + let closest = closest_partial_swap(a, b, -c); + let mut k2l_dag = general.K2l.t().to_owned(); + k2l_dag.mapv_inplace(|x| x.conj()); + TwoQubitWeylDecomposition { + specialization, + a: closest, + b: closest, + c: -closest, + K1l: general.K1l.dot(&general.K2l), + K1r: general.K1r.dot(&ipz).dot(&general.K2l).dot(&ipz), + K2r: ipz.dot(&k2l_dag).dot(&ipz).dot(&general.K2r), + K2l: Array2::eye(2), + ..general + } + } + // :math:`U \sim U_d(\alpha, 0, 0) \sim \text{Ctrl-U}` + // + // This gate binds 4 parameters, we make it canonical by setting: + // + // :math:`K2_l = Ry(\theta_l) Rx(\lambda_l)` , + // :math:`K2_r = Ry(\theta_r) Rx(\lambda_r)` . + Specialization::ControlledEquiv => { + let euler_basis = EulerBasis::XYX; + let [k2ltheta, k2lphi, k2llambda, k2lphase] = + angles_from_unitary(general.K2l.view(), euler_basis); + let [k2rtheta, k2rphi, k2rlambda, k2rphase] = + angles_from_unitary(general.K2r.view(), euler_basis); + TwoQubitWeylDecomposition { + specialization, + a, + b: 0., + c: 0., + global_phase: global_phase + k2lphase + k2rphase, + K1l: general.K1l.dot(&rx_matrix(k2lphi)), + K1r: general.K1r.dot(&rx_matrix(k2rphi)), + K2l: ry_matrix(k2ltheta).dot(&rx_matrix(k2llambda)), + K2r: ry_matrix(k2rtheta).dot(&rx_matrix(k2rlambda)), + default_euler_basis: euler_basis, + ..general + } + } + // :math:`U \sim U_d(\pi/4, \pi/4, \alpha) \sim \text{SWAP} \cdot \text{Ctrl-U}` + // + // This gate binds 4 parameters, we make it canonical by setting: + // + // :math:`K2_l = Ry(\theta_l)\cdot Rz(\lambda_l)` , :math:`K2_r = Ry(\theta_r)\cdot Rz(\lambda_r)` + Specialization::MirrorControlledEquiv => { + let [k2ltheta, k2lphi, k2llambda, k2lphase] = + angles_from_unitary(general.K2l.view(), EulerBasis::ZYZ); + let [k2rtheta, k2rphi, k2rlambda, k2rphase] = + angles_from_unitary(general.K2r.view(), EulerBasis::ZYZ); + TwoQubitWeylDecomposition { + specialization, + a: PI4, + b: PI4, + c, + global_phase: global_phase + k2lphase + k2rphase, + K1l: general.K1l.dot(&rz_matrix(k2rphi)), + K1r: general.K1r.dot(&rz_matrix(k2lphi)), + K2l: ry_matrix(k2ltheta).dot(&rz_matrix(k2llambda)), + K2r: ry_matrix(k2rtheta).dot(&rz_matrix(k2rlambda)), + ..general + } + } + // :math:`U \sim U_d(\alpha, \alpha, \beta), \alpha \geq |\beta|` + // + // This gate binds 5 parameters, we make it canonical by setting: + // + // :math:`K2_l = Ry(\theta_l)\cdot Rz(\lambda_l)`. + Specialization::fSimaabEquiv => { + let [k2ltheta, k2lphi, k2llambda, k2lphase] = + angles_from_unitary(general.K2l.view(), EulerBasis::ZYZ); + TwoQubitWeylDecomposition { + specialization, + a: (a + b) / 2., + b: (a + b) / 2., + c, + global_phase: global_phase + k2lphase, + K1r: general.K1r.dot(&rz_matrix(k2lphi)), + K1l: general.K1l.dot(&rz_matrix(k2lphi)), + K2l: ry_matrix(k2ltheta).dot(&rz_matrix(k2llambda)), + K2r: rz_matrix(-k2lphi).dot(&general.K2r), + ..general + } + } + // :math:`U \sim U_d(\alpha, \beta, -\beta), \alpha \geq \beta \geq 0` + // + // This gate binds 5 parameters, we make it canonical by setting: + // + // :math:`K2_l = Ry(\theta_l)Rx(\lambda_l)` + Specialization::fSimabbEquiv => { + let euler_basis = EulerBasis::XYX; + let [k2ltheta, k2lphi, k2llambda, k2lphase] = + angles_from_unitary(general.K2l.view(), euler_basis); + TwoQubitWeylDecomposition { + specialization, + a, + b: (b + c) / 2., + c: (b + c) / 2., + global_phase: global_phase + k2lphase, + K1r: general.K1r.dot(&rx_matrix(k2lphi)), + K1l: general.K1l.dot(&rx_matrix(k2lphi)), + K2l: ry_matrix(k2ltheta).dot(&rx_matrix(k2llambda)), + K2r: rx_matrix(-k2lphi).dot(&general.K2r), + default_euler_basis: euler_basis, + ..general + } + } + // :math:`U \sim U_d(\alpha, \beta, -\beta), \alpha \geq \beta \geq 0` + // + // This gate binds 5 parameters, we make it canonical by setting: + // + // :math:`K2_l = Ry(\theta_l)Rx(\lambda_l)` + Specialization::fSimabmbEquiv => { + let euler_basis = EulerBasis::XYX; + let [k2ltheta, k2lphi, k2llambda, k2lphase] = + angles_from_unitary(general.K2l.view(), euler_basis); + TwoQubitWeylDecomposition { + specialization, + a, + b: (b - c) / 2., + c: -((b - c) / 2.), + global_phase: global_phase + k2lphase, + K1l: general.K1l.dot(&rx_matrix(k2lphi)), + K1r: general.K1r.dot(&ipz).dot(&rx_matrix(k2lphi)).dot(&ipz), + K2l: ry_matrix(k2ltheta).dot(&rx_matrix(k2llambda)), + K2r: ipz.dot(&rx_matrix(-k2lphi)).dot(&ipz).dot(&general.K2r), + default_euler_basis: euler_basis, + ..general + } + } + // U has no special symmetry. + // + // This gate binds all 6 possible parameters, so there is no need to make the single-qubit + // pre-/post-gates canonical. + Specialization::General => general, + }; + + let tr = if flipped_from_original { + let [da, db, dc] = [ + PI2 - a - specialized.a, + b - specialized.b, + -c - specialized.c, + ]; + 4. * Complex64::new( + da.cos() * db.cos() * dc.cos(), + da.sin() * db.sin() * dc.sin(), + ) + } else { + let [da, db, dc] = [a - specialized.a, b - specialized.b, c - specialized.c]; + 4. * Complex64::new( + da.cos() * db.cos() * dc.cos(), + da.sin() * db.sin() * dc.sin(), + ) + }; + specialized.calculated_fidelity = tr.trace_to_fid(); + if let Some(fid) = specialized.requested_fidelity { + if specialized.calculated_fidelity + 1.0e-13 < fid { + import_exception!(qiskit, QiskitError); + return Err(QiskitError::new_err(format!( + "Specialization: {:?} calculated fidelity: {} is worse than requested fidelity: {}", + specialized.specialization, + specialized.calculated_fidelity, + fid + ))); + } + } + specialized.global_phase += tr.arg(); + Ok(specialized) + } + + #[allow(non_snake_case)] + #[getter] + fn K1l(&self, py: Python) -> PyObject { + self.K1l.to_pyarray(py).into() + } + + #[allow(non_snake_case)] + #[getter] + fn K1r(&self, py: Python) -> PyObject { + self.K1r.to_pyarray(py).into() + } + + #[allow(non_snake_case)] + #[getter] + fn K2l(&self, py: Python) -> PyObject { + self.K2l.to_pyarray(py).into() + } + + #[allow(non_snake_case)] + #[getter] + fn K2r(&self, py: Python) -> PyObject { + self.K2r.to_pyarray(py).into() + } + + #[getter] + fn unitary_matrix(&self, py: Python) -> PyObject { + self.unitary_matrix.to_pyarray(py).into() + } + + #[pyo3(signature = (euler_basis=None, simplify=false, atol=None))] + fn circuit( + &self, + euler_basis: Option<&str>, + simplify: bool, + atol: Option, + ) -> PyResult { + let euler_basis: EulerBasis = match euler_basis { + Some(basis) => EulerBasis::from_string(basis)?, + None => self.default_euler_basis, + }; + let target_1q_basis_list: Vec = vec![euler_basis]; + + let mut gate_sequence = Vec::new(); + let mut global_phase: f64 = self.global_phase; + + let c2r = unitary_to_gate_sequence_inner( + self.K2r.view(), + &target_1q_basis_list, + 0, + None, + simplify, + atol, + ) + .unwrap(); + for gate in c2r.gates { + gate_sequence.push((gate.0, gate.1, smallvec![0])) + } + global_phase += c2r.global_phase; + let c2l = unitary_to_gate_sequence_inner( + self.K2l.view(), + &target_1q_basis_list, + 1, + None, + simplify, + atol, + ) + .unwrap(); + for gate in c2l.gates { + gate_sequence.push((gate.0, gate.1, smallvec![1])) + } + global_phase += c2l.global_phase; + self.weyl_gate( + simplify, + &mut gate_sequence, + atol.unwrap_or(ANGLE_ZERO_EPSILON), + &mut global_phase, + ); + let c1r = unitary_to_gate_sequence_inner( + self.K1r.view(), + &target_1q_basis_list, + 0, + None, + simplify, + atol, + ) + .unwrap(); + for gate in c1r.gates { + gate_sequence.push((gate.0, gate.1, smallvec![0])) + } + global_phase += c2r.global_phase; + let c1l = unitary_to_gate_sequence_inner( + self.K1l.view(), + &target_1q_basis_list, + 1, + None, + simplify, + atol, + ) + .unwrap(); + for gate in c1l.gates { + gate_sequence.push((gate.0, gate.1, smallvec![1])) + } + Ok(TwoQubitGateSequence { + gates: gate_sequence, + global_phase, + }) + } +} + +type TwoQubitSequenceVec = Vec<(String, SmallVec<[f64; 3]>, SmallVec<[u8; 2]>)>; + +#[pyclass(sequence)] +pub struct TwoQubitGateSequence { + gates: TwoQubitSequenceVec, + #[pyo3(get)] + global_phase: f64, +} + +#[pymethods] +impl TwoQubitGateSequence { + #[new] + fn new() -> Self { + TwoQubitGateSequence { + gates: Vec::new(), + global_phase: 0., + } + } + + fn __getstate__(&self) -> (TwoQubitSequenceVec, f64) { + (self.gates.clone(), self.global_phase) + } + + fn __setstate__(&mut self, state: (TwoQubitSequenceVec, f64)) { + self.gates = state.0; + self.global_phase = state.1; + } + + fn __len__(&self) -> PyResult { + Ok(self.gates.len()) + } + + fn __getitem__(&self, py: Python, idx: utils::SliceOrInt) -> PyResult { + match idx { + utils::SliceOrInt::Slice(slc) => { + let len = self.gates.len().try_into().unwrap(); + let indices = slc.indices(len)?; + let mut out_vec: TwoQubitSequenceVec = Vec::new(); + // Start and stop will always be positive the slice api converts + // negatives to the index for example: + // list(range(5))[-1:-3:-1] + // will return start=4, stop=2, and step=- + let mut pos: isize = indices.start; + let mut cond = if indices.step < 0 { + pos > indices.stop + } else { + pos < indices.stop + }; + while cond { + if pos < len as isize { + out_vec.push(self.gates[pos as usize].clone()); + } + pos += indices.step; + if indices.step < 0 { + cond = pos > indices.stop; + } else { + cond = pos < indices.stop; + } + } + Ok(out_vec.into_py(py)) + } + utils::SliceOrInt::Int(idx) => { + let len = self.gates.len() as isize; + if idx >= len || idx < -len { + Err(PyIndexError::new_err(format!("Invalid index, {idx}"))) + } else if idx < 0 { + let len = self.gates.len(); + Ok(self.gates[len - idx.unsigned_abs()].to_object(py)) + } else { + Ok(self.gates[idx as usize].to_object(py)) + } + } + } + } } #[pymodule] pub fn two_qubit_decompose(_py: Python, m: &PyModule) -> PyResult<()> { m.add_wrapped(wrap_pyfunction!(_num_basis_gates))?; + m.add_class::()?; + m.add_class::()?; + m.add_class::()?; Ok(()) } diff --git a/qiskit/__init__.py b/qiskit/__init__.py index 782d8aaabfbd..4d91f6d6e86a 100644 --- a/qiskit/__init__.py +++ b/qiskit/__init__.py @@ -70,6 +70,7 @@ sys.modules["qiskit._accelerate.convert_2q_block_matrix"] = ( qiskit._accelerate.convert_2q_block_matrix ) +sys.modules["qiskit._accelerate.two_qubit_decompose"] = qiskit._accelerate.two_qubit_decompose # qiskit errors operator from qiskit.exceptions import QiskitError, MissingOptionalLibraryError diff --git a/qiskit/synthesis/two_qubit/two_qubit_decompose.py b/qiskit/synthesis/two_qubit/two_qubit_decompose.py index 9985d7cf50c6..d8d6b0d8033d 100644 --- a/qiskit/synthesis/two_qubit/two_qubit_decompose.py +++ b/qiskit/synthesis/two_qubit/two_qubit_decompose.py @@ -29,21 +29,21 @@ import io import base64 import warnings -from typing import ClassVar, Optional, Type +from typing import Optional, Type import logging import numpy as np from qiskit.circuit import QuantumRegister, QuantumCircuit, Gate -from qiskit.circuit.library.standard_gates import CXGate, RXGate, RYGate, RZGate +from qiskit.circuit.library.standard_gates import CXGate, RXGate from qiskit.exceptions import QiskitError from qiskit.quantum_info.operators import Operator -from qiskit.synthesis.two_qubit.weyl import transform_to_magic_basis from qiskit.synthesis.one_qubit.one_qubit_decompose import ( OneQubitEulerDecomposer, DEFAULT_ATOL, ) +from qiskit.utils.deprecation import deprecate_func from qiskit._accelerate import two_qubit_decompose logger = logging.getLogger(__name__) @@ -151,237 +151,30 @@ class TwoQubitWeylDecomposition: requested_fidelity: Optional[float] # None means no automatic specialization calculated_fidelity: float # Fidelity after specialization - _original_decomposition: "TwoQubitWeylDecomposition" - _is_flipped_from_original: bool # The approx is closest to a Weyl reflection of the original? - - _default_1q_basis: ClassVar[str] = "ZYZ" # Default one qubit basis (explicit parameterization) - - def __init_subclass__(cls, **kwargs): - """Subclasses should be concrete, not factories. - - Make explicitly-instantiated subclass __new__ call base __new__ with fidelity=None""" - super().__init_subclass__(**kwargs) - cls.__new__ = lambda cls, *a, fidelity=None, **k: TwoQubitWeylDecomposition.__new__( - cls, *a, fidelity=None, **k - ) - - @staticmethod - def __new__(cls, unitary_matrix, *, fidelity=(1.0 - 1.0e-9), _unpickling=False): - """Perform the Weyl chamber decomposition, and optionally choose a specialized subclass.""" - - # The flip into the Weyl Chamber is described in B. Kraus and J. I. Cirac, Phys. Rev. A 63, - # 062309 (2001). - # - # FIXME: There's a cleaner-seeming method based on choosing branch cuts carefully, in Andrew - # M. Childs, Henry L. Haselgrove, and Michael A. Nielsen, Phys. Rev. A 68, 052311, but I - # wasn't able to get that to work. - # - # The overall decomposition scheme is taken from Drury and Love, arXiv:0806.4015 [quant-ph]. - - if _unpickling: - return super().__new__(cls) - - pi = np.pi - pi2 = np.pi / 2 - pi4 = np.pi / 4 - - # Make U be in SU(4) - U = np.array(unitary_matrix, dtype=complex, copy=True) - detU = np.linalg.det(U) - U *= detU ** (-0.25) - global_phase = cmath.phase(detU) / 4 - - Up = transform_to_magic_basis(U, reverse=True) - M2 = Up.T.dot(Up) - - # M2 is a symmetric complex matrix. We need to decompose it as M2 = P D P^T where - # P ∈ SO(4), D is diagonal with unit-magnitude elements. - # - # We can't use raw `eig` directly because it isn't guaranteed to give us real or othogonal - # eigenvectors. Instead, since `M2` is complex-symmetric, - # M2 = A + iB - # for real-symmetric `A` and `B`, and as - # M2^+ @ M2 = A^2 + B^2 + i [A, B] = 1 - # we must have `A` and `B` commute, and consequently they are simultaneously diagonalizable. - # Mixing them together _should_ account for any degeneracy problems, but it's not - # guaranteed, so we repeat it a little bit. The fixed seed is to make failures - # deterministic; the value is not important. - state = np.random.default_rng(2020) - for _ in range(100): # FIXME: this randomized algorithm is horrendous - M2real = state.normal() * M2.real + state.normal() * M2.imag - _, P = np.linalg.eigh(M2real) - D = P.T.dot(M2).dot(P).diagonal() - if np.allclose(P.dot(np.diag(D)).dot(P.T), M2, rtol=0, atol=1.0e-13): - break - else: - raise QiskitError( - "TwoQubitWeylDecomposition: failed to diagonalize M2." - " Please report this at https://github.com/Qiskit/qiskit-terra/issues/4159." - f" Input: {U.tolist()}" - ) - - d = -np.angle(D) / 2 - d[3] = -d[0] - d[1] - d[2] - cs = np.mod((d[:3] + d[3]) / 2, 2 * np.pi) - - # Reorder the eigenvalues to get in the Weyl chamber - cstemp = np.mod(cs, pi2) - np.minimum(cstemp, pi2 - cstemp, cstemp) - order = np.argsort(cstemp)[[1, 2, 0]] - cs = cs[order] - d[:3] = d[order] - P[:, :3] = P[:, order] - - # Fix the sign of P to be in SO(4) - if np.real(np.linalg.det(P)) < 0: - P[:, -1] = -P[:, -1] - - # Find K1, K2 so that U = K1.A.K2, with K being product of single-qubit unitaries - K1 = transform_to_magic_basis(Up @ P @ np.diag(np.exp(1j * d))) - K2 = transform_to_magic_basis(P.T) - - K1l, K1r, phase_l = decompose_two_qubit_product_gate(K1) - K2l, K2r, phase_r = decompose_two_qubit_product_gate(K2) - global_phase += phase_l + phase_r - - K1l = K1l.copy() - - # Flip into Weyl chamber - if cs[0] > pi2: - cs[0] -= 3 * pi2 - K1l = K1l.dot(_ipy) - K1r = K1r.dot(_ipy) - global_phase += pi2 - if cs[1] > pi2: - cs[1] -= 3 * pi2 - K1l = K1l.dot(_ipx) - K1r = K1r.dot(_ipx) - global_phase += pi2 - conjs = 0 - if cs[0] > pi4: - cs[0] = pi2 - cs[0] - K1l = K1l.dot(_ipy) - K2r = _ipy.dot(K2r) - conjs += 1 - global_phase -= pi2 - if cs[1] > pi4: - cs[1] = pi2 - cs[1] - K1l = K1l.dot(_ipx) - K2r = _ipx.dot(K2r) - conjs += 1 - global_phase += pi2 - if conjs == 1: - global_phase -= pi - if cs[2] > pi2: - cs[2] -= 3 * pi2 - K1l = K1l.dot(_ipz) - K1r = K1r.dot(_ipz) - global_phase += pi2 - if conjs == 1: - global_phase -= pi - if conjs == 1: - cs[2] = pi2 - cs[2] - K1l = K1l.dot(_ipz) - K2r = _ipz.dot(K2r) - global_phase += pi2 - if cs[2] > pi4: - cs[2] -= pi2 - K1l = K1l.dot(_ipz) - K1r = K1r.dot(_ipz) - global_phase -= pi2 - - a, b, c = cs[1], cs[0], cs[2] - - # Save the non-specialized decomposition for later comparison - od = super().__new__(TwoQubitWeylDecomposition) - od.a = a - od.b = b - od.c = c - od.K1l = K1l - od.K1r = K1r - od.K2l = K2l - od.K2r = K2r - od.global_phase = global_phase - od.requested_fidelity = fidelity - od.calculated_fidelity = 1.0 - od.unitary_matrix = np.array(unitary_matrix, dtype=complex, copy=True) - od.unitary_matrix.setflags(write=False) - od._original_decomposition = None - od._is_flipped_from_original = False - - def is_close(ap, bp, cp): - da, db, dc = a - ap, b - bp, c - cp - tr = 4 * complex( - math.cos(da) * math.cos(db) * math.cos(dc), - math.sin(da) * math.sin(db) * math.sin(dc), - ) - fid = trace_to_fid(tr) - return fid >= fidelity - - if fidelity is None: # Don't specialize if None - instance = super().__new__( - TwoQubitWeylGeneral if cls is TwoQubitWeylDecomposition else cls - ) - elif is_close(0, 0, 0): - instance = super().__new__(TwoQubitWeylIdEquiv) - elif is_close(pi4, pi4, pi4) or is_close(pi4, pi4, -pi4): - instance = super().__new__(TwoQubitWeylSWAPEquiv) - elif (lambda x: is_close(x, x, x))(_closest_partial_swap(a, b, c)): - instance = super().__new__(TwoQubitWeylPartialSWAPEquiv) - elif (lambda x: is_close(x, x, -x))(_closest_partial_swap(a, b, -c)): - instance = super().__new__(TwoQubitWeylPartialSWAPFlipEquiv) - elif is_close(a, 0, 0): - instance = super().__new__(TwoQubitWeylControlledEquiv) - elif is_close(pi4, pi4, c): - instance = super().__new__(TwoQubitWeylMirrorControlledEquiv) - elif is_close((a + b) / 2, (a + b) / 2, c): - instance = super().__new__(TwoQubitWeylfSimaabEquiv) - elif is_close(a, (b + c) / 2, (b + c) / 2): - instance = super().__new__(TwoQubitWeylfSimabbEquiv) - elif is_close(a, (b - c) / 2, (c - b) / 2): - instance = super().__new__(TwoQubitWeylfSimabmbEquiv) - else: - instance = super().__new__(TwoQubitWeylGeneral) - - instance._original_decomposition = od - return instance + _specializations = two_qubit_decompose.Specialization def __init__( self, - unitary_matrix: list[list[complex]] | np.ndarray[complex], - fidelity: float | None = None, + unitary_matrix: np.ndarray, + fidelity: float | None = 1.0 - 1.0e-9, + *, + _specialization: two_qubit_decompose.Specialization | None = None, ): - """ - Args: - unitary_matrix: The unitary to decompose. - fidelity: The target fidelity of the decomposed operation. - """ - del unitary_matrix # unused in __init__ (used in new) - od = self._original_decomposition - self.a, self.b, self.c = od.a, od.b, od.c - self.K1l, self.K1r = od.K1l, od.K1r - self.K2l, self.K2r = od.K2l, od.K2r - self.global_phase = od.global_phase - self.unitary_matrix = od.unitary_matrix + unitary_matrix = np.asarray(unitary_matrix, dtype=complex) + self._inner_decomposition = two_qubit_decompose.TwoQubitWeylDecomposition( + unitary_matrix, fidelity=fidelity, _specialization=_specialization + ) + self.a = self._inner_decomposition.a + self.b = self._inner_decomposition.b + self.c = self._inner_decomposition.c + self.global_phase = self._inner_decomposition.global_phase + self.K1l = self._inner_decomposition.K1l + self.K1r = self._inner_decomposition.K1r + self.K2l = self._inner_decomposition.K2l + self.K2r = self._inner_decomposition.K2r + self.unitary_matrix = unitary_matrix self.requested_fidelity = fidelity - self._is_flipped_from_original = False - self.specialize() - - # Update the phase after specialization: - if self._is_flipped_from_original: - da, db, dc = (np.pi / 2 - od.a) - self.a, od.b - self.b, -od.c - self.c - tr = 4 * complex( - math.cos(da) * math.cos(db) * math.cos(dc), - math.sin(da) * math.sin(db) * math.sin(dc), - ) - else: - da, db, dc = od.a - self.a, od.b - self.b, od.c - self.c - tr = 4 * complex( - math.cos(da) * math.cos(db) * math.cos(dc), - math.sin(da) * math.sin(db) * math.sin(dc), - ) - self.global_phase += cmath.phase(tr) - self.calculated_fidelity = trace_to_fid(tr) + self.calculated_fidelity = self._inner_decomposition.calculated_fidelity if logger.isEnabledFor(logging.DEBUG): actual_fidelity = self.actual_fidelity() logger.debug( @@ -395,62 +188,34 @@ def __init__( "Requested fidelity different from actual by %s", self.calculated_fidelity - actual_fidelity, ) - if self.requested_fidelity and self.calculated_fidelity + 1.0e-13 < self.requested_fidelity: - raise QiskitError( - f"{self.__class__.__name__}: " - f"calculated fidelity: {self.calculated_fidelity} " - f"is worse than requested fidelity: {self.requested_fidelity}." - ) + @deprecate_func(since="1.1.0", removal_timeline="in the 2.0.0 release") def specialize(self): - """Make changes to the decomposition to comply with any specialization.""" + """Make changes to the decomposition to comply with any specializations. - # Do update a, b, c, k1l, k1r, k2l, k2r, _is_flipped_from_original to round to the - # specialization. Do not update the global phase, since this gets done in generic - # __init__() + This method will always raise a ``NotImplementedError`` because + there are no specializations to comply with in the current implementation. + """ raise NotImplementedError def circuit( self, *, euler_basis: str | None = None, simplify: bool = False, atol: float = DEFAULT_ATOL ) -> QuantumCircuit: """Returns Weyl decomposition in circuit form.""" - - # simplify, atol arguments are passed to OneQubitEulerDecomposer - if euler_basis is None: - euler_basis = self._default_1q_basis - oneq_decompose = OneQubitEulerDecomposer(euler_basis) - c1l, c1r, c2l, c2r = ( - oneq_decompose(k, simplify=simplify, atol=atol) - for k in (self.K1l, self.K1r, self.K2l, self.K2r) + circuit_sequence = self._inner_decomposition.circuit( + euler_basis=euler_basis, simplify=simplify, atol=atol ) - circ = QuantumCircuit(2, global_phase=self.global_phase) - circ.compose(c2r, [0], inplace=True) - circ.compose(c2l, [1], inplace=True) - self._weyl_gate(simplify, circ, atol) - circ.compose(c1r, [0], inplace=True) - circ.compose(c1l, [1], inplace=True) + circ = QuantumCircuit(2, global_phase=circuit_sequence.global_phase) + for name, params, qubits in circuit_sequence: + getattr(circ, name)(*params, *qubits) return circ - def _weyl_gate(self, simplify, circ: QuantumCircuit, atol): - """Appends U_d(a, b, c) to the circuit. - - Can be overridden in subclasses for special cases""" - if not simplify or abs(self.a) > atol: - circ.rxx(-self.a * 2, 0, 1) - if not simplify or abs(self.b) > atol: - circ.ryy(-self.b * 2, 0, 1) - if not simplify or abs(self.c) > atol: - circ.rzz(-self.c * 2, 0, 1) - def actual_fidelity(self, **kwargs) -> float: """Calculates the actual fidelity of the decomposed circuit to the input unitary.""" circ = self.circuit(**kwargs) trace = np.trace(Operator(circ).data.T.conj() @ self.unitary_matrix) return trace_to_fid(trace) - def __getnewargs_ex__(self): - return (self.unitary_matrix,), {"_unpickling": True} - def __repr__(self): """Represent with enough precision to allow copy-paste debugging of all corner cases""" prefix = f"{type(self).__qualname__}.from_bytes(" @@ -461,12 +226,15 @@ def __repr__(self): b64ascii[-1] += "," pretty = [f"# {x.rstrip()}" for x in str(self).splitlines()] indent = "\n" + " " * 4 + specialization_variant = str(self._inner_decomposition.specialization).split(".")[1] + specialization_repr = f"{type(self).__qualname__}._specializations.{specialization_variant}" lines = ( [prefix] + pretty + b64ascii + [ f"requested_fidelity={self.requested_fidelity},", + f"_specialization={specialization_repr},", f"calculated_fidelity={self.calculated_fidelity},", f"actual_fidelity={self.actual_fidelity()},", f"abc={(self.a, self.b, self.c)})", @@ -476,7 +244,12 @@ def __repr__(self): @classmethod def from_bytes( - cls, bytes_in: bytes, *, requested_fidelity: float, **kwargs + cls, + bytes_in: bytes, + *, + requested_fidelity: float, + _specialization: two_qubit_decompose.Specialization | None, + **kwargs, ) -> "TwoQubitWeylDecomposition": """Decode bytes into :class:`.TwoQubitWeylDecomposition`.""" # Used by __repr__ @@ -484,121 +257,15 @@ def from_bytes( b64 = base64.decodebytes(bytes_in) with io.BytesIO(b64) as f: arr = np.load(f, allow_pickle=False) - return cls(arr, fidelity=requested_fidelity) + return cls(arr, fidelity=requested_fidelity, _specialization=_specialization) def __str__(self): - pre = f"{self.__class__.__name__}(\n\t" + specialization = str(self._inner_decomposition.specialization).split(".")[1] + pre = f"{self.__class__.__name__} [specialization={specialization}] (\n\t" circ_indent = "\n\t".join(self.circuit(simplify=True).draw("text").lines(-1)) return f"{pre}{circ_indent}\n)" -class TwoQubitWeylIdEquiv(TwoQubitWeylDecomposition): - r""":math:`U \sim U_d(0,0,0) \sim Id` - - This gate binds 0 parameters, we make it canonical by setting - :math:`K2_l = Id` , :math:`K2_r = Id`. - """ - - def specialize(self): - self.a = self.b = self.c = 0.0 - self.K1l = self.K1l @ self.K2l - self.K1r = self.K1r @ self.K2r - self.K2l = _id.copy() - self.K2r = _id.copy() - - -class TwoQubitWeylSWAPEquiv(TwoQubitWeylDecomposition): - r""":math:`U \sim U_d(\pi/4, \pi/4, \pi/4) \sim U(\pi/4, \pi/4, -\pi/4) \sim \text{SWAP}` - - This gate binds 0 parameters, we make it canonical by setting - :math:`K2_l = Id` , :math:`K2_r = Id`. - """ - - def specialize(self): - if self.c > 0: - self.K1l = self.K1l @ self.K2r - self.K1r = self.K1r @ self.K2l - else: - self._is_flipped_from_original = True - self.K1l = self.K1l @ _ipz @ self.K2r - self.K1r = self.K1r @ _ipz @ self.K2l - self.global_phase = self.global_phase + np.pi / 2 - self.a = self.b = self.c = np.pi / 4 - self.K2l = _id.copy() - self.K2r = _id.copy() - - def _weyl_gate(self, simplify, circ: QuantumCircuit, atol): - del self, simplify, atol # unused - circ.swap(0, 1) - circ.global_phase -= 3 * np.pi / 4 - - -def _closest_partial_swap(a, b, c) -> float: - r"""A good approximation to the best value x to get the minimum - trace distance for :math:`U_d(x, x, x)` from :math:`U_d(a, b, c)`. - """ - m = (a + b + c) / 3 - am, bm, cm = a - m, b - m, c - m - ab, bc, ca = a - b, b - c, c - a - - return m + am * bm * cm * (6 + ab * ab + bc * bc + ca * ca) / 18 - - -class TwoQubitWeylPartialSWAPEquiv(TwoQubitWeylDecomposition): - r""":math:`U \sim U_d(\alpha\pi/4, \alpha\pi/4, \alpha\pi/4) \sim \text{SWAP}^\alpha` - This gate binds 3 parameters, we make it canonical by setting: - :math:`K2_l = Id`. - """ - - def specialize(self): - self.a = self.b = self.c = _closest_partial_swap(self.a, self.b, self.c) - self.K1l = self.K1l @ self.K2l - self.K1r = self.K1r @ self.K2l - self.K2r = self.K2l.T.conj() @ self.K2r - self.K2l = _id.copy() - - -class TwoQubitWeylPartialSWAPFlipEquiv(TwoQubitWeylDecomposition): - r""":math:`U \sim U_d(\alpha\pi/4, \alpha\pi/4, -\alpha\pi/4) \sim \text{SWAP}^\alpha` - (a non-equivalent root of SWAP from the TwoQubitWeylPartialSWAPEquiv - similar to how :math:`x = (\pm \sqrt(x))^2`) - This gate binds 3 parameters, we make it canonical by setting: - :math:`K2_l = Id`. - """ - - def specialize(self): - self.a = self.b = _closest_partial_swap(self.a, self.b, -self.c) - self.c = -self.a - self.K1l = self.K1l @ self.K2l - self.K1r = self.K1r @ _ipz @ self.K2l @ _ipz - self.K2r = _ipz @ self.K2l.T.conj() @ _ipz @ self.K2r - self.K2l = _id.copy() - - -_oneq_xyx = OneQubitEulerDecomposer("XYX") -_oneq_zyz = OneQubitEulerDecomposer("ZYZ") - - -class TwoQubitWeylControlledEquiv(TwoQubitWeylDecomposition): - r""":math:`U \sim U_d(\alpha, 0, 0) \sim \text{Ctrl-U}` - This gate binds 4 parameters, we make it canonical by setting: - :math:`K2_l = Ry(\theta_l) Rx(\lambda_l)` , - :math:`K2_r = Ry(\theta_r) Rx(\lambda_r)` . - """ - - _default_1q_basis = "XYX" - - def specialize(self): - self.b = self.c = 0 - k2ltheta, k2lphi, k2llambda, k2lphase = _oneq_xyx.angles_and_phase(self.K2l) - k2rtheta, k2rphi, k2rlambda, k2rphase = _oneq_xyx.angles_and_phase(self.K2r) - self.global_phase += k2lphase + k2rphase - self.K1l = self.K1l @ np.asarray(RXGate(k2lphi)) - self.K1r = self.K1r @ np.asarray(RXGate(k2rphi)) - self.K2l = np.asarray(RYGate(k2ltheta)) @ np.asarray(RXGate(k2llambda)) - self.K2r = np.asarray(RYGate(k2rtheta)) @ np.asarray(RXGate(k2rlambda)) - - class TwoQubitControlledUDecomposer: r"""Decompose two-qubit unitary in terms of a desired :math:`U \sim U_d(\alpha, 0, 0) \sim \text{Ctrl-U}` @@ -627,18 +294,23 @@ def __init__(self, rxx_equivalent_gate: Type[Gate]): circ = QuantumCircuit(2) circ.rxx(test_angle, 0, 1) - decomposer_rxx = TwoQubitWeylControlledEquiv(Operator(circ).data) + decomposer_rxx = TwoQubitWeylDecomposition( + Operator(circ).data, + fidelity=None, + _specialization=two_qubit_decompose.Specialization.ControlledEquiv, + ) circ = QuantumCircuit(2) circ.append(rxx_equivalent_gate(test_angle), qargs=[0, 1]) - decomposer_equiv = TwoQubitWeylControlledEquiv(Operator(circ).data) + decomposer_equiv = TwoQubitWeylDecomposition( + Operator(circ).data, + fidelity=None, + _specialization=two_qubit_decompose.Specialization.ControlledEquiv, + ) scale = decomposer_rxx.a / decomposer_equiv.a - if ( - not isinstance(decomp, TwoQubitWeylControlledEquiv) - or abs(decomp.a * 2 - test_angle / scale) > atol - ): + if abs(decomp.a * 2 - test_angle / scale) > atol: raise QiskitError( f"{rxx_equivalent_gate.__name__} is not equivalent to an RXXGate." ) @@ -706,7 +378,7 @@ def _to_rxx_gate(self, angle: float) -> QuantumCircuit: circ = QuantumCircuit(2) circ.append(self.rxx_equivalent_gate(self.scale * angle), qargs=[0, 1]) - decomposer_inv = TwoQubitWeylControlledEquiv(Operator(circ).data) + decomposer_inv = TwoQubitWeylDecomposition(Operator(circ).data) oneq_decompose = OneQubitEulerDecomposer("ZYZ") @@ -763,97 +435,6 @@ def _weyl_gate(self, circ: QuantumCircuit, atol=1.0e-13): return circ -class TwoQubitWeylMirrorControlledEquiv(TwoQubitWeylDecomposition): - r""":math:`U \sim U_d(\pi/4, \pi/4, \alpha) \sim \text{SWAP} \cdot \text{Ctrl-U}` - - This gate binds 4 parameters, we make it canonical by setting: - :math:`K2_l = Ry(\theta_l)\cdot Rz(\lambda_l)` , :math:`K2_r = Ry(\theta_r)\cdot Rz(\lambda_r)`. - """ - - def specialize(self): - self.a = self.b = np.pi / 4 - k2ltheta, k2lphi, k2llambda, k2lphase = _oneq_zyz.angles_and_phase(self.K2l) - k2rtheta, k2rphi, k2rlambda, k2rphase = _oneq_zyz.angles_and_phase(self.K2r) - self.global_phase += k2lphase + k2rphase - self.K1r = self.K1r @ np.asarray(RZGate(k2lphi)) - self.K1l = self.K1l @ np.asarray(RZGate(k2rphi)) - self.K2l = np.asarray(RYGate(k2ltheta)) @ np.asarray(RZGate(k2llambda)) - self.K2r = np.asarray(RYGate(k2rtheta)) @ np.asarray(RZGate(k2rlambda)) - - def _weyl_gate(self, simplify, circ: QuantumCircuit, atol): - circ.swap(0, 1) - circ.rzz((np.pi / 4 - self.c) * 2, 0, 1) - circ.global_phase += np.pi / 4 - - -# These next 3 gates use the definition of fSim from https://arxiv.org/pdf/2001.08343.pdf eq (1) -class TwoQubitWeylfSimaabEquiv(TwoQubitWeylDecomposition): - r""":math:`U \sim U_d(\alpha, \alpha, \beta), \alpha \geq |\beta|` - - This gate binds 5 parameters, we make it canonical by setting: - :math:`K2_l = Ry(\theta_l)\cdot Rz(\lambda_l)`. - """ - - def specialize(self): - self.a = self.b = (self.a + self.b) / 2 - k2ltheta, k2lphi, k2llambda, k2lphase = _oneq_zyz.angles_and_phase(self.K2l) - self.global_phase += k2lphase - self.K1r = self.K1r @ np.asarray(RZGate(k2lphi)) - self.K1l = self.K1l @ np.asarray(RZGate(k2lphi)) - self.K2l = np.asarray(RYGate(k2ltheta)) @ np.asarray(RZGate(k2llambda)) - self.K2r = np.asarray(RZGate(-k2lphi)) @ self.K2r - - -class TwoQubitWeylfSimabbEquiv(TwoQubitWeylDecomposition): - r""":math:`U \sim U_d(\alpha, \beta, -\beta), \alpha \geq \beta \geq 0` - - This gate binds 5 parameters, we make it canonical by setting: - :math:`K2_l = Ry(\theta_l)Rx(\lambda_l)`. - """ - - _default_1q_basis = "XYX" - - def specialize(self): - self.b = self.c = (self.b + self.c) / 2 - k2ltheta, k2lphi, k2llambda, k2lphase = _oneq_xyx.angles_and_phase(self.K2l) - self.global_phase += k2lphase - self.K1r = self.K1r @ np.asarray(RXGate(k2lphi)) - self.K1l = self.K1l @ np.asarray(RXGate(k2lphi)) - self.K2l = np.asarray(RYGate(k2ltheta)) @ np.asarray(RXGate(k2llambda)) - self.K2r = np.asarray(RXGate(-k2lphi)) @ self.K2r - - -class TwoQubitWeylfSimabmbEquiv(TwoQubitWeylDecomposition): - r""":math:`U \sim U_d(\alpha, \beta, -\beta), \alpha \geq \beta \geq 0` - - This gate binds 5 parameters, we make it canonical by setting: - :math:`K2_l = Ry(\theta_l)Rx(\lambda_l)`. - """ - - _default_1q_basis = "XYX" - - def specialize(self): - self.b = (self.b - self.c) / 2 - self.c = -self.b - k2ltheta, k2lphi, k2llambda, k2lphase = _oneq_xyx.angles_and_phase(self.K2l) - self.global_phase += k2lphase - self.K1r = self.K1r @ _ipz @ np.asarray(RXGate(k2lphi)) @ _ipz - self.K1l = self.K1l @ np.asarray(RXGate(k2lphi)) - self.K2l = np.asarray(RYGate(k2ltheta)) @ np.asarray(RXGate(k2llambda)) - self.K2r = _ipz @ np.asarray(RXGate(-k2lphi)) @ _ipz @ self.K2r - - -class TwoQubitWeylGeneral(TwoQubitWeylDecomposition): - """U has no special symmetry. - - This gate binds all 6 possible parameters, so there is no need to make the single-qubit - pre-/post-gates canonical. - """ - - def specialize(self): - pass # Nothing to do - - def Ud(a, b, c): r"""Generates the array :math:`e^{(i a XX + i b YY + i c ZZ)}`""" return np.array( diff --git a/releasenotes/notes/rust-two-qubit-weyl-ec551f3f9c812124.yaml b/releasenotes/notes/rust-two-qubit-weyl-ec551f3f9c812124.yaml new file mode 100644 index 000000000000..12d5c27d9af8 --- /dev/null +++ b/releasenotes/notes/rust-two-qubit-weyl-ec551f3f9c812124.yaml @@ -0,0 +1,27 @@ +--- +features_synthesis: + - | + The :class:`.TwoQubitWeylDecomposition` synthesis class has been rewritten + in Rust for better performance. +upgrade_synthesis: + - | + The :class:`.TwoQubitWeylDecomposition` no longer will self-specialize into + a subclass on creation. This was an internal detail of the :class:`.TwoQubitWeylDecomposition` + previously, and was not a documented public behavior as all the subclasses behaved + the same and were only used for internal dispatch. However, as it was discoverable + behavior this release note is to document that this will no longer occur and all + instances of :class:`.TwoQubitWeylDecomposition` will be of the same type. There is no + change in behavior for public methods of the class. + +deprecations_synthesis: + - | + The :meth:`.TwoQubitWeylDecomposition.specialize` method is now deprecated + and will be removed in the Qiskit 2.0.0 release. This method never had + a public purpose and was unsafe for an end user to call as it would + mutate the calculated decomposition in the object and produce invalid + fields in the object. It was only used internally to construct a new + :class:`.TwoQubitWeylDecomposition` object. Despite this it was still a + documented part of the public API for the class and is now being + deprecated without any potential replacement. This release it always will + raise a ``NotImplementedError`` when called because the specialization + subclassing has been removed as part of the Rust rewrite of the class. diff --git a/test/python/synthesis/test_synthesis.py b/test/python/synthesis/test_synthesis.py index a5b360deb6cf..d3f8560cb9b2 100644 --- a/test/python/synthesis/test_synthesis.py +++ b/test/python/synthesis/test_synthesis.py @@ -56,16 +56,6 @@ from qiskit.synthesis.one_qubit.one_qubit_decompose import OneQubitEulerDecomposer from qiskit.synthesis.two_qubit.two_qubit_decompose import ( TwoQubitWeylDecomposition, - TwoQubitWeylIdEquiv, - TwoQubitWeylSWAPEquiv, - TwoQubitWeylPartialSWAPEquiv, - TwoQubitWeylPartialSWAPFlipEquiv, - TwoQubitWeylfSimaabEquiv, - TwoQubitWeylfSimabbEquiv, - TwoQubitWeylfSimabmbEquiv, - TwoQubitWeylControlledEquiv, - TwoQubitWeylMirrorControlledEquiv, - TwoQubitWeylGeneral, two_qubit_cnot_decompose, TwoQubitBasisDecomposer, TwoQubitControlledUDecomposer, @@ -73,6 +63,7 @@ decompose_two_qubit_product_gate, TwoQubitDecomposeUpToDiagonal, ) +from qiskit._accelerate.two_qubit_decompose import Specialization from qiskit.synthesis.unitary import qsd from test import combine # pylint: disable=wrong-import-order from test import QiskitTestCase # pylint: disable=wrong-import-order @@ -229,8 +220,16 @@ def check_two_qubit_weyl_specialization( # Loop to check both for implicit and explicity specialization for decomposer in (TwoQubitWeylDecomposition, expected_specialization): - with self.assertDebugOnly(): - decomp = decomposer(target_unitary, fidelity=fidelity) + if isinstance(decomposer, TwoQubitWeylDecomposition): + with self.assertDebugOnly(): + decomp = decomposer(target_unitary, fidelity=fidelity) + decomp_name = decomp.specialization + else: + with self.assertDebugOnly(): + decomp = TwoQubitWeylDecomposition( + target_unitary, fidelity=None, _specialization=expected_specialization + ) + decomp_name = expected_specialization self.assertRoundTrip(decomp) self.assertRoundTripPickle(decomp) self.assertEqual( @@ -238,27 +237,33 @@ def check_two_qubit_weyl_specialization( 0, "Incorrect saved unitary in the decomposition.", ) - self.assertIsInstance(decomp, expected_specialization, "Incorrect Weyl specialization.") + self.assertEqual( + decomp._inner_decomposition.specialization, + expected_specialization, + "Incorrect Weyl specialization.", + ) circ = decomp.circuit(simplify=True) self.assertDictEqual( - dict(circ.count_ops()), expected_gates, f"Gate counts of {decomposer.__name__}" + dict(circ.count_ops()), expected_gates, f"Gate counts of {decomp_name}" ) actual_fid = decomp.actual_fidelity() self.assertAlmostEqual(decomp.calculated_fidelity, actual_fid, places=13) - self.assertGreaterEqual(actual_fid, fidelity, f"fidelity of {decomposer.__name__}") + self.assertGreaterEqual(actual_fid, fidelity, f"fidelity of {decomp_name}") actual_unitary = Operator(circ).data trace = np.trace(actual_unitary.T.conj() @ target_unitary) - self.assertAlmostEqual( - trace.imag, 0, places=13, msg=f"Real trace for {decomposer.__name__}" - ) + self.assertAlmostEqual(trace.imag, 0, places=13, msg=f"Real trace for {decomp_name}") with self.assertDebugOnly(): - decomp2 = expected_specialization(target_unitary, fidelity=None) # Shouldn't raise + decomp2 = TwoQubitWeylDecomposition( + target_unitary, fidelity=None, _specialization=expected_specialization + ) # Shouldn't raise self.assertRoundTrip(decomp2) self.assertRoundTripPickle(decomp2) - if expected_specialization is not TwoQubitWeylGeneral: + if expected_specialization != Specialization.General: with self.assertRaises(QiskitError) as exc: - _ = expected_specialization(target_unitary, fidelity=1.0) - self.assertIn("worse than requested", exc.exception.message) + _ = TwoQubitWeylDecomposition( + target_unitary, fidelity=1.0, _specialization=expected_specialization + ) + self.assertIn("worse than requested", str(exc.exception)) def check_exact_decomposition( self, target_unitary, decomposer, tolerance=1.0e-12, num_basis_uses=None @@ -828,7 +833,7 @@ def test_weyl_specialize_id(self): self.check_two_qubit_weyl_specialization( k1 @ Ud(a + da, b + db, c + dc) @ k2, 0.999, - TwoQubitWeylIdEquiv, + Specialization.IdEquiv, {"rz": 4, "ry": 2}, ) @@ -842,7 +847,7 @@ def test_weyl_specialize_swap(self): self.check_two_qubit_weyl_specialization( k1 @ Ud(a + da, b + db, c + dc) @ k2, 0.999, - TwoQubitWeylSWAPEquiv, + Specialization.SWAPEquiv, {"rz": 4, "ry": 2, "swap": 1}, ) @@ -856,7 +861,7 @@ def test_weyl_specialize_flip_swap(self): self.check_two_qubit_weyl_specialization( k1 @ Ud(a + da, b + db, c + dc) @ k2, 0.999, - TwoQubitWeylSWAPEquiv, + Specialization.SWAPEquiv, {"rz": 4, "ry": 2, "swap": 1}, ) @@ -870,7 +875,7 @@ def test_weyl_specialize_pswap(self, theta=0.123): self.check_two_qubit_weyl_specialization( k1 @ Ud(a + da, b + db, c + dc) @ k2, 0.999, - TwoQubitWeylPartialSWAPEquiv, + Specialization.PartialSWAPEquiv, {"rz": 6, "ry": 3, "rxx": 1, "ryy": 1, "rzz": 1}, ) @@ -884,7 +889,7 @@ def test_weyl_specialize_flip_pswap(self, theta=0.123): self.check_two_qubit_weyl_specialization( k1 @ Ud(a + da, b + db, c + dc) @ k2, 0.999, - TwoQubitWeylPartialSWAPFlipEquiv, + Specialization.PartialSWAPFlipEquiv, {"rz": 6, "ry": 3, "rxx": 1, "ryy": 1, "rzz": 1}, ) @@ -898,7 +903,7 @@ def test_weyl_specialize_fsim_aab(self, aaa=0.456, bbb=0.132): self.check_two_qubit_weyl_specialization( k1 @ Ud(a + da, b + db, c + dc) @ k2, 0.999, - TwoQubitWeylfSimaabEquiv, + Specialization.fSimaabEquiv, {"rz": 7, "ry": 4, "rxx": 1, "ryy": 1, "rzz": 1}, ) @@ -912,7 +917,7 @@ def test_weyl_specialize_fsim_abb(self, aaa=0.456, bbb=0.132): self.check_two_qubit_weyl_specialization( k1 @ Ud(a + da, b + db, c + dc) @ k2, 0.999, - TwoQubitWeylfSimabbEquiv, + Specialization.fSimabbEquiv, {"rx": 7, "ry": 4, "rxx": 1, "ryy": 1, "rzz": 1}, ) @@ -926,7 +931,7 @@ def test_weyl_specialize_fsim_abmb(self, aaa=0.456, bbb=0.132): self.check_two_qubit_weyl_specialization( k1 @ Ud(a + da, b + db, c + dc) @ k2, 0.999, - TwoQubitWeylfSimabmbEquiv, + Specialization.fSimabmbEquiv, {"rx": 7, "ry": 4, "rxx": 1, "ryy": 1, "rzz": 1}, ) @@ -940,7 +945,7 @@ def test_weyl_specialize_ctrl(self, aaa=0.456): self.check_two_qubit_weyl_specialization( k1 @ Ud(a + da, b + db, c + dc) @ k2, 0.999, - TwoQubitWeylControlledEquiv, + Specialization.ControlledEquiv, {"rx": 6, "ry": 4, "rxx": 1}, ) @@ -954,7 +959,7 @@ def test_weyl_specialize_mirror_ctrl(self, aaa=-0.456): self.check_two_qubit_weyl_specialization( k1 @ Ud(a + da, b + db, c + dc) @ k2, 0.999, - TwoQubitWeylMirrorControlledEquiv, + Specialization.MirrorControlledEquiv, {"rz": 6, "ry": 4, "rzz": 1, "swap": 1}, ) @@ -968,7 +973,7 @@ def test_weyl_specialize_general(self, aaa=0.456, bbb=0.345, ccc=0.123): self.check_two_qubit_weyl_specialization( k1 @ Ud(a + da, b + db, c + dc) @ k2, 0.999, - TwoQubitWeylGeneral, + Specialization.General, {"rz": 8, "ry": 4, "rxx": 1, "ryy": 1, "rzz": 1}, )