From dae02aa02d37c9203660c45222263b8f7eae397f Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Sun, 25 Feb 2024 09:22:38 -0500 Subject: [PATCH] Oxidize TwoQubitWeylDecomposition This commit is part 1 of migrating the default 2q unitary synthesis to leverage parallel rust described in #8774, the eventual goal is to be able to run unitary synthesis in parallel for all the unitary matrices in a circuit in a single call from the `UnitarySynthesis` pass. This commit lays the initial groundwork for doing this by starting with the largest piece of the default 2q unitary synthesis code, the TwoQubitWeylDecomposition class. It migrates the entire class to be a pyclass in rust. There is still a Python subclass for it that handles the actual QuantumCircuit generation and also the string representations which are dependent on circuit generation. The goal of this is to keep the same basic algorithm in place but re-implement as-is in Rust as a common starting point for eventual improvements to the underlying algorithm as well as parallelizing the synthesis of multiple 2q unitary matrices. --- Cargo.lock | 12 + crates/accelerate/Cargo.toml | 6 +- .../src/euler_one_qubit_decomposer.rs | 39 +- crates/accelerate/src/two_qubit_decompose.rs | 929 +++++++++++++++++- crates/accelerate/src/utils.rs | 21 +- .../two_qubit/two_qubit_decompose.py | 511 +--------- test/python/synthesis/test_synthesis.py | 200 +--- 7 files changed, 982 insertions(+), 736 deletions(-) diff --git a/Cargo.lock b/Cargo.lock index 1fb29d80b45a..3bc457d6f65c 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", @@ -1149,6 +1160,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..10e2dee6f031 100644 --- a/crates/accelerate/Cargo.toml +++ b/crates/accelerate/Cargo.toml @@ -37,7 +37,11 @@ features = ["hashbrown", "indexmap", "num-complex", "num-bigint"] [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..dbfa0535da4e 100644 --- a/crates/accelerate/src/euler_one_qubit_decomposer.rs +++ b/crates/accelerate/src/euler_one_qubit_decomposer.rs @@ -27,7 +27,7 @@ use numpy::PyReadonlyArray2; use crate::utils::SliceOrInt; -const DEFAULT_ATOL: f64 = 1e-12; +pub const DEFAULT_ATOL: f64 = 1e-12; #[pyclass(module = "qiskit._accelerate.euler_one_qubit_decomposer")] pub struct OneQubitGateErrorMap { @@ -61,9 +61,9 @@ impl OneQubitGateErrorMap { #[pyclass(sequence)] pub struct OneQubitGateSequence { - gates: Vec<(String, Vec)>, + pub gates: Vec<(String, Vec)>, #[pyo3(get)] - global_phase: f64, + pub global_phase: f64, } #[pymethods] @@ -553,7 +553,7 @@ pub fn generate_circuit( } #[inline] -fn angles_from_unitary(unitary: ArrayView2, target_basis: &str) -> [f64; 4] { +pub fn angles_from_unitary(unitary: ArrayView2, target_basis: &str) -> [f64; 4] { match target_basis { "U321" => params_u3_inner(unitary), "U3" => params_u3_inner(unitary), @@ -609,6 +609,7 @@ fn compute_error( } } +#[inline] #[pyfunction] pub fn compute_error_one_qubit_sequence( circuit: &OneQubitGateSequence, @@ -618,6 +619,7 @@ pub fn compute_error_one_qubit_sequence( compute_error(&circuit.gates, error_map, qubit) } +#[inline] #[pyfunction] pub fn compute_error_list( circuit: Vec<(String, Vec)>, @@ -648,7 +650,26 @@ pub fn unitary_to_gate_sequence( } } let unitary_mat = unitary.as_array(); - let best_result = target_basis_list + Ok(unitary_to_gate_sequence_inner( + unitary_mat, + &target_basis_list, + qubit, + error_map, + simplify, + atol, + )) +} + +#[inline] +pub fn unitary_to_gate_sequence_inner( + unitary_mat: ArrayView2, + target_basis_list: &[&str], + 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); @@ -658,12 +679,11 @@ pub fn unitary_to_gate_sequence( 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]] } @@ -697,6 +717,7 @@ fn params_zxz_inner(mat: ArrayView2) -> [f64; 4] { [theta, phi + PI / 2., lam - PI / 2., phase] } +#[inline] #[pyfunction] pub fn params_zyz(unitary: PyReadonlyArray2) -> [f64; 4] { let mat = unitary.as_array(); @@ -712,6 +733,7 @@ fn params_u3_inner(mat: ArrayView2) -> [f64; 4] { [theta, phi, lam, phase - 0.5 * (phi + lam)] } +#[inline] #[pyfunction] pub fn params_u3(unitary: PyReadonlyArray2) -> [f64; 4] { let mat = unitary.as_array(); @@ -726,6 +748,7 @@ fn params_u1x_inner(mat: ArrayView2) -> [f64; 4] { [theta, phi, lam, phase - 0.5 * (theta + phi + lam)] } +#[inline] #[pyfunction] pub fn params_u1x(unitary: PyReadonlyArray2) -> [f64; 4] { let mat = unitary.as_array(); diff --git a/crates/accelerate/src/two_qubit_decompose.rs b/crates/accelerate/src/two_qubit_decompose.rs index d9f451a43682..5d34a9b6ad93 100644 --- a/crates/accelerate/src/two_qubit_decompose.rs +++ b/crates/accelerate/src/two_qubit_decompose.rs @@ -18,59 +18,118 @@ // 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, PyTypeError, PyValueError}; use pyo3::prelude::*; use pyo3::wrap_pyfunction; use pyo3::Python; 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, DEFAULT_ATOL, +}; 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 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.), + ], ]; -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_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.), + ], ]; -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 +fn transform_from_magic_basis(unitary: ArrayView2, reverse: bool) -> Array2 { + let _b_nonnormalized = aview2(&B_NON_NORMALIZED); + let _b_nonnormalized_dagger = aview2(&B_NON_NORMALIZED_DAGGER); + if reverse { + _b_nonnormalized_dagger.dot(&unitary).dot(&_b_nonnormalized) + } else { + _b_nonnormalized.dot(&unitary).dot(&_b_nonnormalized_dagger) + } +} + +fn transform_from_magic_basis_faer(u: Mat, reverse: bool) -> Mat { + let unitary: ArrayView2 = u.as_ref().into_ndarray_complex(); + let _b_nonnormalized = aview2(&B_NON_NORMALIZED); + let _b_nonnormalized_dagger = aview2(&B_NON_NORMALIZED_DAGGER); + if reverse { + _b_nonnormalized_dagger.dot(&unitary).dot(&_b_nonnormalized) + } else { + _b_nonnormalized.dot(&unitary).dot(&_b_nonnormalized_dagger) + } + .view() + .into_faer_complex() + .to_owned() } // faer::c64 and num_complex::Complex are both structs @@ -98,9 +157,33 @@ impl Arg for c64 { } } +fn decompose_two_qubit_product_gate( + unitary: ArrayView2, +) -> (Array2, Array2, f64) { + let mut r: Array2 = unitary.slice(s![..2, ..2]).to_owned(); + let mut det_r = det_one_qubit(r.view()); + if det_r.abs() < 0.1 { + r = unitary.slice(s![2.., ..2]).to_owned(); + det_r = det_one_qubit(r.view()); + } + r.mapv_inplace(|x| x / det_r.sqrt()); + let r_t_conj: Array2 = r.t().mapv(|x| x.conj()).to_owned(); + let eye = array![ + [Complex64::new(1., 0.), Complex64::new(0., 0.)], + [Complex64::new(0., 0.), Complex64::new(1., 0.)], + ]; + let mut temp = kron(&eye, &r_t_conj); + temp = unitary.dot(&temp); + let mut l = temp.slice_mut(s![..;2, ..;2]).to_owned(); + let det_l = det_one_qubit(l.view()); + l.mapv_inplace(|x| x / det_l.sqrt()); + let phase = det_l.arg() / 2.; + (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); + let uup = transform_from_magic_basis_faer(uscaled, true); let mut darg: Vec<_> = (uup.transpose() * &uup) .complex_eigenvalues() .into_iter() @@ -180,7 +263,12 @@ 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_to_fid_c64(&trace) * basis_fidelity.powi(idx as i32), + ) + }) .min_by(|(_idx1, fid1), (_idx2, fid2)| fid2.partial_cmp(fid1).unwrap()) .unwrap() .0 @@ -188,12 +276,779 @@ fn __num_basis_gates(basis_b: f64, basis_fidelity: f64, unitary: MatRef) -> /// 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 { +fn trace_to_fid_c64(trace: &c64) -> f64 { (4.0 + trace.faer_abs2()) / 20.0 } +fn trace_to_fid(trace: Complex64) -> f64 { + (4.0 + trace.abs().powi(2)) / 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 isin = Complex64::new(0., half_theta.sin()); + array![[cos, -isin], [isin, 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)] +#[pyclass(module = "qiskit._accelerate.two_qubit_decompose")] +enum Specializations { + 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 + SimaabEquiv, + SimabbEquiv, + SimabmbEquiv, +} + +#[derive(Clone)] +#[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, + specialization: Specializations, + default_euler_basis: String, + #[pyo3(get)] + requested_fidelity: Option, + #[pyo3(get)] + calculated_fidelity: f64, + unitary_matrix: Array2, +} + +impl TwoQubitWeylDecomposition { + fn weyl_gate( + &mut self, + simplify: bool, + sequence: &mut Vec<(String, Vec, [u8; 2])>, + atol: f64, + ) { + match self.specialization { + Specializations::MirrorControlledEquiv => { + sequence.push(("swap".to_string(), Vec::new(), [0, 1])); + sequence.push(("rzz".to_string(), vec![PI4 - self.c], [0, 1])); + self.global_phase += PI4; + } + Specializations::SWAPEquiv => { + sequence.push(("swap".to_string(), Vec::new(), [0, 1])); + self.global_phase -= 3. * PI / 4.; + } + _ => { + if !simplify || self.a.abs() > atol { + sequence.push(("rxx".to_string(), vec![-self.a * 2.], [0, 1])); + } + if !simplify || self.b.abs() > atol { + sequence.push(("ryy".to_string(), vec![-self.b * 2.], [0, 1])); + } + if !simplify || self.c.abs() > atol { + sequence.push(("rzz".to_string(), vec![-self.c * 2.], [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 { + #[new] + #[pyo3(signature=(unitary_matrix, fidelity=DEFAULT_FIDELITY, specialization=None))] + fn new( + unitary_matrix: PyReadonlyArray2, + fidelity: Option, + specialization: Option, + ) -> PyResult { + 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(); + u.mapv_inplace(|x| x * det_u.powf(-0.25)); + let mut global_phase = det_u.arg() / 4.; + let u_p = transform_from_magic_basis(u.view(), true); + // Use ndarray here because matmul precision in faer is lower, it seems + // to more aggressively round to zero which causes different behaviour + // during the eigen decomposition below. + let m2 = u_p.t().dot(&u_p); + let default_euler_basis = "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 _ in 0..100 { + let rand_a: f64 = state.sample(StandardNormal); + let rand_b: f64 = state.sample(StandardNormal); + let m2_real = m2.mapv(|val| rand_a * val.re + rand_b * val.im).to_owned(); + let p_inner = m2_real + .view() + .into_faer() + .selfadjoint_eigendecomposition(Lower) + .u() + .into_ndarray() + .mapv(Complex64::from) + .to_owned(); + // Uncomment this to use numpy for eigh instead of faer (useful if needed to compare) + + // let m2_real = m2.as_ref().into_ndarray().map(|val| rand_a * val.re + rand_b * val.im).to_owned(); + // let numpy_linalg = PyModule::import(py, "numpy.linalg")?; + // let eigh = numpy_linalg.getattr("eigh")?; + // let m2_real_arr = m2_real.to_pyarray(py); + // let result = eigh.call1((m2_real_arr,))?.downcast::()?; + // let p_raw = result.get_item(1)?; + // let p_arr = p_raw.extract::>()?.as_array().to_owned(); + // let p_inner = p_arr.mapv(Complex64::from).to_owned(); + + // let m2_arr: ArrayView2 = m2.as_ref().into_ndarray_complex(); + 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()).to_owned(); + found = abs_diff_eq!(compare.view(), m2, epsilon = 1.0e-13); + if found { + p = p_inner; + d = d_inner; + break; + } + } + if !found { + return Err(PyTypeError::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: Array1 = (0..3) + .map(|i| ((d[i] + d[3]) / 2.0).rem_euclid(TWO_PI)) + .collect(); + let cstemp: Vec = 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 = transform_from_magic_basis(u_p.dot(&p).dot(&temp).view(), false); + let k2 = transform_from_magic_basis(p.t(), false); + + #[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]; + conjs += 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) => trace_to_fid(tr) >= 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.) { + Specializations::IdEquiv + } else if is_close(PI4, PI4, PI4) || is_close(PI4, PI4, -PI4) { + Specializations::SWAPEquiv + } else if is_close(closest_abc, closest_abc, closest_abc) { + Specializations::PartialSWAPEquiv + } else if is_close(closest_ab_minus_c, closest_ab_minus_c, -closest_ab_minus_c) { + Specializations::PartialSWAPFlipEquiv + } else if is_close(a, 0., 0.) { + Specializations::ControlledEquiv + } else if is_close(PI4, PI4, c) { + Specializations::MirrorControlledEquiv + } else if is_close((a + b) / 2., (a + b) / 2., c) { + Specializations::SimaabEquiv + } else if is_close(a, (b + c) / 2., (b + c) / 2.) { + Specializations::SimabbEquiv + } else if is_close(a, (b - c) / 2., (c - b) / 2.) { + Specializations::SimabmbEquiv + } else { + Specializations::General + } + } + }; + + let mut specialized: TwoQubitWeylDecomposition = match specialization { + Specializations::IdEquiv => TwoQubitWeylDecomposition { + a: 0., + b: 0., + c: 0., + global_phase, + K1l: K1l.dot(&K2l), + K1r: K1r.dot(&K2r), + K2l: Array2::eye(2), + K2r: Array2::eye(2), + specialization: Specializations::IdEquiv, + default_euler_basis: default_euler_basis.to_string(), + requested_fidelity: fidelity, + calculated_fidelity: 1.0, + unitary_matrix, + }, + Specializations::SWAPEquiv => { + if c > 0. { + TwoQubitWeylDecomposition { + a: PI4, + b: PI4, + c: PI4, + global_phase, + K1l: K1l.dot(&K2r), + K1r: K1r.dot(&K2l), + K2l: Array2::eye(2), + K2r: Array2::eye(2), + specialization: Specializations::SWAPEquiv, + default_euler_basis: default_euler_basis.to_string(), + requested_fidelity: fidelity, + calculated_fidelity: 1.0, + unitary_matrix, + } + } else { + flipped_from_original = true; + TwoQubitWeylDecomposition { + a: PI4, + b: PI4, + c: PI4, + global_phase: global_phase + PI2, + K1l: K1l.dot(&ipz).dot(&K2r), + K1r: K1r.dot(&ipz).dot(&K2l), + K2l: Array2::eye(2), + K2r: Array2::eye(2), + specialization: Specializations::SWAPEquiv, + default_euler_basis: default_euler_basis.to_string(), + requested_fidelity: fidelity, + calculated_fidelity: 1.0, + unitary_matrix, + } + } + } + Specializations::PartialSWAPEquiv => { + let closest = closest_partial_swap(a, b, c); + let mut k2r_temp = K2l.t().to_owned(); + k2r_temp.view_mut().mapv_inplace(|x| x.conj()); + TwoQubitWeylDecomposition { + a: closest, + b: closest, + c: closest, + global_phase, + K1l: K1l.dot(&K2l), + K1r: K1r.dot(&K2l), + K2r: k2r_temp.dot(&K2r), + K2l: Array2::eye(2), + specialization: Specializations::PartialSWAPEquiv, + default_euler_basis: default_euler_basis.to_string(), + requested_fidelity: fidelity, + calculated_fidelity: 1.0, + unitary_matrix, + } + } + Specializations::PartialSWAPFlipEquiv => { + let closest = closest_partial_swap(a, b, -c); + let mut k2_temp = K2l.t().to_owned(); + k2_temp.mapv_inplace(|x| x.conj()); + TwoQubitWeylDecomposition { + a: closest, + b: closest, + c: -a, + global_phase, + K1l: K1l.dot(&K2l), + K1r: K1r.dot(&ipz).dot(&K2l).dot(&ipz), + K2r: ipz.dot(&k2_temp).dot(&ipz).dot(&K2r), + K2l: Array2::eye(2), + specialization: Specializations::PartialSWAPFlipEquiv, + default_euler_basis: default_euler_basis.to_string(), + requested_fidelity: fidelity, + calculated_fidelity: 1.0, + unitary_matrix, + } + } + Specializations::ControlledEquiv => { + let default_euler_basis = "XYX"; + let [k2ltheta, k2lphi, k2llambda, k2lphase] = + angles_from_unitary(K2l.view(), "XYX"); + let [k2rtheta, k2rphi, k2rlambda, k2rphase] = + angles_from_unitary(K2r.view(), "XYX"); + TwoQubitWeylDecomposition { + a, + b: 0., + c: 0., + global_phase: global_phase + k2lphase + k2rphase, + K1l: K1l.dot(&rx_matrix(k2lphi)), + K1r: K1r.dot(&rx_matrix(k2rphi)), + K2l: ry_matrix(k2ltheta).dot(&rx_matrix(k2llambda)), + K2r: ry_matrix(k2rtheta).dot(&rx_matrix(k2rlambda)), + specialization: Specializations::ControlledEquiv, + default_euler_basis: default_euler_basis.to_string(), + requested_fidelity: fidelity, + calculated_fidelity: 1.0, + unitary_matrix, + } + } + Specializations::MirrorControlledEquiv => { + let [k2ltheta, k2lphi, k2llambda, k2lphase] = + angles_from_unitary(K2l.view(), "ZYZ"); + let [k2rtheta, k2rphi, k2rlambda, k2rphase] = + angles_from_unitary(K2r.view(), "ZYZ"); + TwoQubitWeylDecomposition { + a: PI4, + b: PI4, + c, + global_phase: global_phase + k2lphase + k2rphase, + K1l: K1l.dot(&rz_matrix(k2lphi)), + K1r: K1r.dot(&rz_matrix(k2rphi)), + K2l: ry_matrix(k2ltheta).dot(&rz_matrix(k2llambda)), + K2r: ry_matrix(k2rtheta).dot(&rz_matrix(k2rlambda)), + specialization: Specializations::MirrorControlledEquiv, + default_euler_basis: default_euler_basis.to_string(), + requested_fidelity: fidelity, + calculated_fidelity: 1.0, + unitary_matrix, + } + } + Specializations::SimaabEquiv => { + let [k2ltheta, k2lphi, k2llambda, k2lphase] = + angles_from_unitary(K2l.view(), "ZYZ"); + TwoQubitWeylDecomposition { + a: (a + b) / 2., + b: (a + b) / 2., + c, + global_phase: global_phase + k2lphase, + K1r: K1r.dot(&rz_matrix(k2lphi)), + K1l: K1l.dot(&rz_matrix(k2lphi)), + K2l: ry_matrix(k2ltheta).dot(&rz_matrix(k2llambda)), + K2r: rz_matrix(-k2lphi).dot(&K2r), + specialization: Specializations::SimaabEquiv, + default_euler_basis: default_euler_basis.to_string(), + requested_fidelity: fidelity, + calculated_fidelity: 1.0, + unitary_matrix, + } + } + Specializations::SimabbEquiv => { + let default_euler_basis = "XYX"; + let [k2ltheta, k2lphi, k2llambda, k2lphase] = + angles_from_unitary(K2l.view(), "XYX"); + TwoQubitWeylDecomposition { + a, + b: (b + c) / 2., + c: (b + c) / 2., + global_phase: global_phase + k2lphase, + K1r: K1r.dot(&rx_matrix(k2lphi)), + K1l: K1l.dot(&rx_matrix(k2lphi)), + K2l: ry_matrix(k2ltheta).dot(&rz_matrix(k2llambda)), + K2r: ry_matrix(-k2lphi).dot(&K2r), + specialization: Specializations::SimabbEquiv, + default_euler_basis: default_euler_basis.to_string(), + requested_fidelity: fidelity, + calculated_fidelity: 1.0, + unitary_matrix, + } + } + Specializations::SimabmbEquiv => { + // TwoQubitWeylfSimabmbEquiv + let default_euler_basis = "XYX"; + let [k2ltheta, k2lphi, k2llambda, k2lphase] = + angles_from_unitary(K2l.view(), "XYX"); + TwoQubitWeylDecomposition { + a, + b: (b - c) / 2., + c: (b - c) / 2., + global_phase: global_phase + k2lphase, + K1l: K1l.dot(&rz_matrix(k2lphi)), + K1r: 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(&K2r), + specialization: Specializations::SimabmbEquiv, + default_euler_basis: default_euler_basis.to_string(), + requested_fidelity: fidelity, + calculated_fidelity: 1.0, + unitary_matrix, + } + } + Specializations::General => TwoQubitWeylDecomposition { + a, + b, + c, + global_phase, + K1l, + K2l, + K1r, + K2r, + specialization: Specializations::General, + default_euler_basis: default_euler_basis.to_string(), + requested_fidelity: fidelity, + calculated_fidelity: 1.0, + unitary_matrix, + }, + }; + + 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 = trace_to_fid(tr); + if let Some(fid) = specialized.requested_fidelity { + if specialized.calculated_fidelity + 1.0e-13 < fid { + return Err(PyValueError::new_err("Uh oh")); + } + } + 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( + &mut self, + euler_basis: Option<&str>, + simplify: bool, + atol: Option, + ) -> TwoQubitGateSequence { + let binding = self.default_euler_basis.clone(); + let euler_basis: &str = euler_basis.unwrap_or(&binding); + let target_1q_basis_list: Vec<&str> = 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, [0, 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, [1, 1])) + } + global_phase += c2l.global_phase; + self.weyl_gate(simplify, &mut gate_sequence, atol.unwrap_or(DEFAULT_ATOL)); + 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, [0, 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, [1, 1])) + } + TwoQubitGateSequence { + gates: gate_sequence, + global_phase, + } + } +} + +type TwoQubitSequenceVec = Vec<(String, Vec, [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: Vec<(String, Vec, [u8; 2])> = 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/crates/accelerate/src/utils.rs b/crates/accelerate/src/utils.rs index bc8a66568f9b..90d74a9bd516 100644 --- a/crates/accelerate/src/utils.rs +++ b/crates/accelerate/src/utils.rs @@ -15,7 +15,9 @@ use pyo3::types::PySlice; use faer::prelude::*; use faer::IntoFaerComplex; -use num_complex::Complex; +use faer::IntoNdarray; +use faer::Side::Lower; +use num_complex::{Complex, Complex64}; use numpy::{IntoPyArray, PyReadonlyArray2}; /// A private enumeration type used to extract arguments to pymethod @@ -53,8 +55,25 @@ pub fn eigenvalues(py: Python, unitary: PyReadonlyArray2>) -> PyObj .into() } +/// Return the eigenvalues of `unitary` as a one-dimensional `numpy.ndarray` +/// with `dtype(complex128)`. +#[pyfunction] +#[pyo3(text_signature = "(unitary, /")] +pub fn eigh(py: Python, unitary: PyReadonlyArray2) -> PyObject { + unitary + .as_array() + .into_faer() + .selfadjoint_eigendecomposition(Lower) + .u() + .into_ndarray() + .mapv(Complex64::from) + .into_pyarray(py) + .into() +} + #[pymodule] pub fn utils(_py: Python, m: &PyModule) -> PyResult<()> { m.add_wrapped(wrap_pyfunction!(eigenvalues))?; + m.add_wrapped(wrap_pyfunction!(eigh))?; Ok(()) } diff --git a/qiskit/synthesis/two_qubit/two_qubit_decompose.py b/qiskit/synthesis/two_qubit/two_qubit_decompose.py index c8af4002fb27..b34603ebd90e 100644 --- a/qiskit/synthesis/two_qubit/two_qubit_decompose.py +++ b/qiskit/synthesis/two_qubit/two_qubit_decompose.py @@ -96,7 +96,7 @@ def decompose_two_qubit_product_gate(special_unitary_matrix: np.ndarray): _id = np.array([[1, 0], [0, 1]], dtype=complex) -class TwoQubitWeylDecomposition: +class TwoQubitWeylDecomposition(two_qubit_decompose.TwoQubitWeylDecomposition): r"""Two-qubit Weyl decomposition. Decompose two-qubit unitary @@ -151,306 +151,26 @@ 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 - - def __init__( - self, - unitary_matrix: list[list[complex]] | np.ndarray[complex], - fidelity: float | 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 - 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) - if logger.isEnabledFor(logging.DEBUG): - actual_fidelity = self.actual_fidelity() - logger.debug( - "Requested fidelity: %s calculated fidelity: %s actual fidelity %s", - self.requested_fidelity, - self.calculated_fidelity, - actual_fidelity, - ) - if abs(self.calculated_fidelity - actual_fidelity) > 1.0e-12: - logger.warning( - "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}." - ) - - def specialize(self): - """Make changes to the decomposition to comply with any specialization.""" - - # 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__() - 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) - ) - 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) + circuit_sequence = super().circuit() + circ = QuantumCircuit(2, global_phase=circuit_sequence.global_phase) + for name, params, qubits in circuit_sequence: + if qubits[0] == qubits[1]: + qargs = (qubits[0],) + else: + qargs = tuple(qubits) + getattr(circ, name)(*params, *qargs) 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(" @@ -492,113 +212,6 @@ def __str__(self): 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}` @@ -623,22 +236,19 @@ def __init__(self, rxx_equivalent_gate: Type[Gate]): rxx_equivalent_gate(test_angle, label="foo") except TypeError as _: raise QiskitError("Equivalent gate needs to take exactly 1 angle parameter.") from _ - decomp = TwoQubitWeylDecomposition(rxx_equivalent_gate(test_angle)) + decomp = TwoQubitWeylDecomposition(rxx_equivalent_gate(test_angle).to_matrix()) 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.Specializations.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.Specializations.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." ) @@ -663,7 +273,7 @@ def __call__(self, unitary, *, atol=DEFAULT_ATOL) -> QuantumCircuit: """ # pylint: disable=attribute-defined-outside-init - self.decomposer = TwoQubitWeylDecomposition(unitary) + self.decomposer = TwoQubitWeylDecomposition(np.asarray(unitary, dtype=complex)) oneq_decompose = OneQubitEulerDecomposer("ZYZ") c1l, c1r, c2l, c2r = ( @@ -706,7 +316,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 +373,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/test/python/synthesis/test_synthesis.py b/test/python/synthesis/test_synthesis.py index a5b360deb6cf..be8b1f5f203f 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, @@ -141,24 +131,10 @@ def check_one_qubit_euler_angles(self, operator, basis="U3", tolerance=1e-14, si np.abs(maxdist) < tolerance, f"Operator {operator}: Worst distance {maxdist}" ) - @contextlib.contextmanager - def assertDebugOnly(self): # FIXME: when at python 3.10+ replace with assertNoLogs - """Context manager, asserts log is emitted at level DEBUG but no higher""" - with self.assertLogs("qiskit.synthesis", "DEBUG") as ctx: - yield - for i in range(len(ctx.records)): - self.assertLessEqual( - ctx.records[i].levelno, - logging.DEBUG, - msg=f"Unexpected logging entry: {ctx.output[i]}", - ) - self.assertIn("Requested fidelity:", ctx.records[i].getMessage()) - def assertRoundTrip(self, weyl1: TwoQubitWeylDecomposition): """Fail if eval(repr(weyl1)) not equal to weyl1""" repr1 = repr(weyl1) - with self.assertDebugOnly(): - weyl2: TwoQubitWeylDecomposition = eval(repr1) # pylint: disable=eval-used + weyl2: TwoQubitWeylDecomposition = eval(repr1) # pylint: disable=eval-used msg_base = f"weyl1:\n{repr1}\nweyl2:\n{repr(weyl2)}" self.assertEqual(type(weyl1), type(weyl2), msg_base) maxdiff = np.max(abs(weyl1.unitary_matrix - weyl2.unitary_matrix)) @@ -201,8 +177,7 @@ def assertRoundTripPickle(self, weyl1: TwoQubitWeylDecomposition): def check_two_qubit_weyl_decomposition(self, target_unitary, tolerance=1.0e-12): """Check TwoQubitWeylDecomposition() works for a given operator""" # pylint: disable=invalid-name - with self.assertDebugOnly(): - decomp = TwoQubitWeylDecomposition(target_unitary, fidelity=None) + decomp = TwoQubitWeylDecomposition(target_unitary, fidelity=None) # self.assertRoundTrip(decomp) # Too slow op = np.exp(1j * decomp.global_phase) * Operator(np.eye(4)) for u, qs in ( @@ -229,8 +204,7 @@ 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) + decomp = decomposer(target_unitary, fidelity=fidelity) self.assertRoundTrip(decomp) self.assertRoundTripPickle(decomp) self.assertEqual( @@ -251,8 +225,7 @@ def check_two_qubit_weyl_specialization( self.assertAlmostEqual( trace.imag, 0, places=13, msg=f"Real trace for {decomposer.__name__}" ) - with self.assertDebugOnly(): - decomp2 = expected_specialization(target_unitary, fidelity=None) # Shouldn't raise + decomp2 = expected_specialization(target_unitary, fidelity=None) # Shouldn't raise self.assertRoundTrip(decomp2) self.assertRoundTripPickle(decomp2) if expected_specialization is not TwoQubitWeylGeneral: @@ -635,13 +608,13 @@ class TestTwoQubitWeylDecomposition(CheckDecompositions): def test_TwoQubitWeylDecomposition_repr(self, seed=42): """Check that eval(__repr__) is exact round trip""" target = random_unitary(4, seed=seed) - weyl1 = TwoQubitWeylDecomposition(target, fidelity=0.99) + weyl1 = TwoQubitWeylDecomposition(target.data, fidelity=0.99) self.assertRoundTrip(weyl1) def test_TwoQubitWeylDecomposition_pickle(self, seed=42): """Check that loads(dumps()) is exact round trip""" target = random_unitary(4, seed=seed) - weyl1 = TwoQubitWeylDecomposition(target, fidelity=0.99) + weyl1 = TwoQubitWeylDecomposition(target.data, fidelity=0.99) self.assertRoundTripPickle(weyl1) def test_two_qubit_weyl_decomposition_cnot(self): @@ -815,164 +788,6 @@ def test_two_qubit_weyl_decomposition_abc(self, smallest=1e-18, factor=9.8, step ] -class TestTwoQubitWeylDecompositionSpecialization(CheckDecompositions): - """Check TwoQubitWeylDecomposition specialized subclasses""" - - def test_weyl_specialize_id(self): - """Weyl specialization for Id gate""" - a, b, c = 0.0, 0.0, 0.0 - for da, db, dc in DELTAS: - for k1l, k1r, k2l, k2r in K1K2SB: - k1 = np.kron(k1l.data, k1r.data) - k2 = np.kron(k2l.data, k2r.data) - self.check_two_qubit_weyl_specialization( - k1 @ Ud(a + da, b + db, c + dc) @ k2, - 0.999, - TwoQubitWeylIdEquiv, - {"rz": 4, "ry": 2}, - ) - - def test_weyl_specialize_swap(self): - """Weyl specialization for swap gate""" - a, b, c = np.pi / 4, np.pi / 4, np.pi / 4 - for da, db, dc in DELTAS: - for k1l, k1r, k2l, k2r in K1K2SB: - k1 = np.kron(k1l.data, k1r.data) - k2 = np.kron(k2l.data, k2r.data) - self.check_two_qubit_weyl_specialization( - k1 @ Ud(a + da, b + db, c + dc) @ k2, - 0.999, - TwoQubitWeylSWAPEquiv, - {"rz": 4, "ry": 2, "swap": 1}, - ) - - def test_weyl_specialize_flip_swap(self): - """Weyl specialization for flip swap gate""" - a, b, c = np.pi / 4, np.pi / 4, -np.pi / 4 - for da, db, dc in DELTAS: - for k1l, k1r, k2l, k2r in K1K2SB: - k1 = np.kron(k1l.data, k1r.data) - k2 = np.kron(k2l.data, k2r.data) - self.check_two_qubit_weyl_specialization( - k1 @ Ud(a + da, b + db, c + dc) @ k2, - 0.999, - TwoQubitWeylSWAPEquiv, - {"rz": 4, "ry": 2, "swap": 1}, - ) - - def test_weyl_specialize_pswap(self, theta=0.123): - """Weyl specialization for partial swap gate""" - a, b, c = theta, theta, theta - for da, db, dc in DELTAS: - for k1l, k1r, k2l, k2r in K1K2SB: - k1 = np.kron(k1l.data, k1r.data) - k2 = np.kron(k2l.data, k2r.data) - self.check_two_qubit_weyl_specialization( - k1 @ Ud(a + da, b + db, c + dc) @ k2, - 0.999, - TwoQubitWeylPartialSWAPEquiv, - {"rz": 6, "ry": 3, "rxx": 1, "ryy": 1, "rzz": 1}, - ) - - def test_weyl_specialize_flip_pswap(self, theta=0.123): - """Weyl specialization for flipped partial swap gate""" - a, b, c = theta, theta, -theta - for da, db, dc in DELTAS: - for k1l, k1r, k2l, k2r in K1K2SB: - k1 = np.kron(k1l.data, k1r.data) - k2 = np.kron(k2l.data, k2r.data) - self.check_two_qubit_weyl_specialization( - k1 @ Ud(a + da, b + db, c + dc) @ k2, - 0.999, - TwoQubitWeylPartialSWAPFlipEquiv, - {"rz": 6, "ry": 3, "rxx": 1, "ryy": 1, "rzz": 1}, - ) - - def test_weyl_specialize_fsim_aab(self, aaa=0.456, bbb=0.132): - """Weyl specialization for partial swap gate""" - a, b, c = aaa, aaa, bbb - for da, db, dc in DELTAS: - for k1l, k1r, k2l, k2r in K1K2SB: - k1 = np.kron(k1l.data, k1r.data) - k2 = np.kron(k2l.data, k2r.data) - self.check_two_qubit_weyl_specialization( - k1 @ Ud(a + da, b + db, c + dc) @ k2, - 0.999, - TwoQubitWeylfSimaabEquiv, - {"rz": 7, "ry": 4, "rxx": 1, "ryy": 1, "rzz": 1}, - ) - - def test_weyl_specialize_fsim_abb(self, aaa=0.456, bbb=0.132): - """Weyl specialization for partial swap gate""" - a, b, c = aaa, bbb, bbb - for da, db, dc in DELTAS: - for k1l, k1r, k2l, k2r in K1K2SB: - k1 = np.kron(k1l.data, k1r.data) - k2 = np.kron(k2l.data, k2r.data) - self.check_two_qubit_weyl_specialization( - k1 @ Ud(a + da, b + db, c + dc) @ k2, - 0.999, - TwoQubitWeylfSimabbEquiv, - {"rx": 7, "ry": 4, "rxx": 1, "ryy": 1, "rzz": 1}, - ) - - def test_weyl_specialize_fsim_abmb(self, aaa=0.456, bbb=0.132): - """Weyl specialization for partial swap gate""" - a, b, c = aaa, bbb, -bbb - for da, db, dc in DELTAS: - for k1l, k1r, k2l, k2r in K1K2SB: - k1 = np.kron(k1l.data, k1r.data) - k2 = np.kron(k2l.data, k2r.data) - self.check_two_qubit_weyl_specialization( - k1 @ Ud(a + da, b + db, c + dc) @ k2, - 0.999, - TwoQubitWeylfSimabmbEquiv, - {"rx": 7, "ry": 4, "rxx": 1, "ryy": 1, "rzz": 1}, - ) - - def test_weyl_specialize_ctrl(self, aaa=0.456): - """Weyl specialization for partial swap gate""" - a, b, c = aaa, 0.0, 0.0 - for da, db, dc in DELTAS: - for k1l, k1r, k2l, k2r in K1K2SB: - k1 = np.kron(k1l.data, k1r.data) - k2 = np.kron(k2l.data, k2r.data) - self.check_two_qubit_weyl_specialization( - k1 @ Ud(a + da, b + db, c + dc) @ k2, - 0.999, - TwoQubitWeylControlledEquiv, - {"rx": 6, "ry": 4, "rxx": 1}, - ) - - def test_weyl_specialize_mirror_ctrl(self, aaa=-0.456): - """Weyl specialization for partial swap gate""" - a, b, c = np.pi / 4, np.pi / 4, aaa - for da, db, dc in DELTAS: - for k1l, k1r, k2l, k2r in K1K2SB: - k1 = np.kron(k1l.data, k1r.data) - k2 = np.kron(k2l.data, k2r.data) - self.check_two_qubit_weyl_specialization( - k1 @ Ud(a + da, b + db, c + dc) @ k2, - 0.999, - TwoQubitWeylMirrorControlledEquiv, - {"rz": 6, "ry": 4, "rzz": 1, "swap": 1}, - ) - - def test_weyl_specialize_general(self, aaa=0.456, bbb=0.345, ccc=0.123): - """Weyl specialization for partial swap gate""" - a, b, c = aaa, bbb, ccc - for da, db, dc in DELTAS: - for k1l, k1r, k2l, k2r in K1K2SB: - k1 = np.kron(k1l.data, k1r.data) - k2 = np.kron(k2l.data, k2r.data) - self.check_two_qubit_weyl_specialization( - k1 @ Ud(a + da, b + db, c + dc) @ k2, - 0.999, - TwoQubitWeylGeneral, - {"rz": 8, "ry": 4, "rxx": 1, "ryy": 1, "rzz": 1}, - ) - - @ddt class TestTwoQubitDecompose(CheckDecompositions): """Test TwoQubitBasisDecomposer() for exact/approx decompositions""" @@ -1079,8 +894,7 @@ def test_approx_supercontrolled_decompose_random(self, seed): tgt = random_unitary(4, seed=state).data tgt *= np.exp(1j * tgt_phase) - with self.assertDebugOnly(): - traces_pred = decomposer.traces(TwoQubitWeylDecomposition(tgt)) + traces_pred = decomposer.traces(TwoQubitWeylDecomposition(tgt)) for i in range(4): with self.subTest(i=i):