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):