From c7d675fd9832a510fed349b225e6a13c87765652 Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Sun, 25 Feb 2024 09:22:38 -0500 Subject: [PATCH 01/46] 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 9985d7cf50c6..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): From ce569292ef1023ff1fc2bea28865964d82aef88e Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Wed, 6 Mar 2024 06:11:20 -0500 Subject: [PATCH 02/46] Fix typo in closest_partial_swap() This commit fixes a typo the formula in the function. This is the same fix from #11953. Co-authored-by: Shelly Garion <46566946+ShellyGarion@users.noreply.github.com> --- crates/accelerate/src/two_qubit_decompose.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/crates/accelerate/src/two_qubit_decompose.rs b/crates/accelerate/src/two_qubit_decompose.rs index 5d34a9b6ad93..23c09d59f7df 100644 --- a/crates/accelerate/src/two_qubit_decompose.rs +++ b/crates/accelerate/src/two_qubit_decompose.rs @@ -290,7 +290,7 @@ 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. + m + am * bm * cm * (6. + ab * ab + bc * bc + ca * ca) / 18. } fn rx_matrix(theta: f64) -> Array2 { From f2323dbe3b288b40587c5edec10ffc19301ceee0 Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Wed, 6 Mar 2024 13:42:58 -0500 Subject: [PATCH 03/46] Run black and lint --- .../two_qubit/two_qubit_decompose.py | 19 +++++++++++++------ test/python/synthesis/test_synthesis.py | 9 +++------ 2 files changed, 16 insertions(+), 12 deletions(-) diff --git a/qiskit/synthesis/two_qubit/two_qubit_decompose.py b/qiskit/synthesis/two_qubit/two_qubit_decompose.py index b34603ebd90e..0ad159d4d567 100644 --- a/qiskit/synthesis/two_qubit/two_qubit_decompose.py +++ b/qiskit/synthesis/two_qubit/two_qubit_decompose.py @@ -29,17 +29,16 @@ import io import base64 import warnings -from typing import ClassVar, Optional, Type +from typing import Optional, Type import logging import numpy as np from qiskit.circuit import QuantumRegister, QuantumCircuit, Gate -from qiskit.circuit.library.standard_gates import CXGate, RXGate, RYGate, RZGate +from qiskit.circuit.library.standard_gates import CXGate, RXGate from qiskit.exceptions import QiskitError from qiskit.quantum_info.operators import Operator -from qiskit.synthesis.two_qubit.weyl import transform_to_magic_basis from qiskit.synthesis.one_qubit.one_qubit_decompose import ( OneQubitEulerDecomposer, DEFAULT_ATOL, @@ -155,7 +154,7 @@ def circuit( self, *, euler_basis: str | None = None, simplify: bool = False, atol: float = DEFAULT_ATOL ) -> QuantumCircuit: """Returns Weyl decomposition in circuit form.""" - circuit_sequence = super().circuit() + circuit_sequence = super().circuit(euler_basis=euler_basis, simplify=simplify, atol=atol) circ = QuantumCircuit(2, global_phase=circuit_sequence.global_phase) for name, params, qubits in circuit_sequence: if qubits[0] == qubits[1]: @@ -240,11 +239,19 @@ def __init__(self, rxx_equivalent_gate: Type[Gate]): circ = QuantumCircuit(2) circ.rxx(test_angle, 0, 1) - decomposer_rxx = TwoQubitWeylDecomposition(Operator(circ).data, fidelity=None, specialization=two_qubit_decompose.Specializations.ControlledEquiv) + 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 = TwoQubitWeylDecomposition(Operator(circ).data, fidelity=None, specialization=two_qubit_decompose.Specializations.ControlledEquiv) + decomposer_equiv = TwoQubitWeylDecomposition( + Operator(circ).data, + fidelity=None, + specialization=two_qubit_decompose.Specializations.ControlledEquiv, + ) scale = decomposer_rxx.a / decomposer_equiv.a diff --git a/test/python/synthesis/test_synthesis.py b/test/python/synthesis/test_synthesis.py index be8b1f5f203f..76ac5606aaac 100644 --- a/test/python/synthesis/test_synthesis.py +++ b/test/python/synthesis/test_synthesis.py @@ -14,8 +14,6 @@ import pickle import unittest -import contextlib -import logging import math import numpy as np import scipy @@ -228,10 +226,9 @@ def check_two_qubit_weyl_specialization( decomp2 = expected_specialization(target_unitary, fidelity=None) # Shouldn't raise self.assertRoundTrip(decomp2) self.assertRoundTripPickle(decomp2) - if expected_specialization is not TwoQubitWeylGeneral: - with self.assertRaises(QiskitError) as exc: - _ = expected_specialization(target_unitary, fidelity=1.0) - self.assertIn("worse than requested", exc.exception.message) + with self.assertRaises(QiskitError) as exc: + _ = expected_specialization(target_unitary, fidelity=1.0) + self.assertIn("worse than requested", exc.exception.message) def check_exact_decomposition( self, target_unitary, decomposer, tolerance=1.0e-12, num_basis_uses=None From 201c7631d02d3e61ad58fe68d7a2dbc5a0a75f39 Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Thu, 7 Mar 2024 08:46:22 -0500 Subject: [PATCH 04/46] Fix potential imaginary component sign flip in determinant --- crates/accelerate/src/two_qubit_decompose.rs | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/crates/accelerate/src/two_qubit_decompose.rs b/crates/accelerate/src/two_qubit_decompose.rs index 23c09d59f7df..8d7f86ffe53d 100644 --- a/crates/accelerate/src/two_qubit_decompose.rs +++ b/crates/accelerate/src/two_qubit_decompose.rs @@ -318,7 +318,7 @@ fn rz_matrix(theta: f64) -> Array2 { const DEFAULT_FIDELITY: f64 = 1.0 - 1.0e-9; const C1_IM: Complex64 = Complex64::new(0.0, 1.0); -#[derive(Clone)] +#[derive(Clone, Debug)] #[pyclass(module = "qiskit._accelerate.two_qubit_decompose")] enum Specializations { General, @@ -417,10 +417,15 @@ impl TwoQubitWeylDecomposition { 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)); + // Faer sometimes returns negative 0s which will throw off the signs + // after the powf we do below, normalize to 0. instead by adding a + // zero complex. + let det_u = u.view().into_faer_complex().determinant().to_num_complex() + Complex64::new(0., 0.); + let det_pow = det_u.powf(-0.25); + u.mapv_inplace(|x| x * det_pow); let mut global_phase = det_u.arg() / 4.; let u_p = transform_from_magic_basis(u.view(), true); // Use ndarray here because matmul precision in faer is lower, it seems @@ -476,6 +481,7 @@ impl TwoQubitWeylDecomposition { .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 { From f087ee899618479b754ff3e668b8c04e2d318a44 Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Thu, 7 Mar 2024 08:48:06 -0500 Subject: [PATCH 05/46] Run cargo fmt --- crates/accelerate/src/two_qubit_decompose.rs | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/crates/accelerate/src/two_qubit_decompose.rs b/crates/accelerate/src/two_qubit_decompose.rs index 8d7f86ffe53d..7024527c5b3c 100644 --- a/crates/accelerate/src/two_qubit_decompose.rs +++ b/crates/accelerate/src/two_qubit_decompose.rs @@ -423,7 +423,8 @@ impl TwoQubitWeylDecomposition { // Faer sometimes returns negative 0s which will throw off the signs // after the powf we do below, normalize to 0. instead by adding a // zero complex. - let det_u = u.view().into_faer_complex().determinant().to_num_complex() + Complex64::new(0., 0.); + let det_u = + u.view().into_faer_complex().determinant().to_num_complex() + Complex64::new(0., 0.); let det_pow = det_u.powf(-0.25); u.mapv_inplace(|x| x * det_pow); let mut global_phase = det_u.arg() / 4.; From 368c950cf61c5f229cd41ed2bf84f33ab6cb0bea Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Thu, 7 Mar 2024 09:01:59 -0500 Subject: [PATCH 06/46] Simplify using numpy eigh example comment --- crates/accelerate/src/two_qubit_decompose.rs | 20 +++++++++----------- 1 file changed, 9 insertions(+), 11 deletions(-) diff --git a/crates/accelerate/src/two_qubit_decompose.rs b/crates/accelerate/src/two_qubit_decompose.rs index 7024527c5b3c..12d2cc17bdcb 100644 --- a/crates/accelerate/src/two_qubit_decompose.rs +++ b/crates/accelerate/src/two_qubit_decompose.rs @@ -464,17 +464,15 @@ impl TwoQubitWeylDecomposition { .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 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_inner = p_raw + // .extract::>()? + // .as_array() + // .mapv(Complex64::from); let d_inner = p_inner.t().dot(&m2).dot(&p_inner).diag().to_owned(); let mut diag_d: Array2 = Array2::zeros((4, 4)); diag_d From 7302a2957ff4031f2670ab1ef57e86f98d89b03a Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Sun, 10 Mar 2024 10:18:16 +0900 Subject: [PATCH 07/46] Add missing checks to decompose_two_qubit_product_gate() --- crates/accelerate/src/two_qubit_decompose.rs | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/crates/accelerate/src/two_qubit_decompose.rs b/crates/accelerate/src/two_qubit_decompose.rs index 12d2cc17bdcb..89225bbbeb6e 100644 --- a/crates/accelerate/src/two_qubit_decompose.rs +++ b/crates/accelerate/src/two_qubit_decompose.rs @@ -166,6 +166,10 @@ fn decompose_two_qubit_product_gate( r = unitary.slice(s![2.., ..2]).to_owned(); det_r = det_one_qubit(r.view()); } + assert!( + det_r.abs() < 0.1, + "decompose_two_qubit_product_gate: unable to decompose: detR < 0 .1" + ); r.mapv_inplace(|x| x / det_r.sqrt()); let r_t_conj: Array2 = r.t().mapv(|x| x.conj()).to_owned(); let eye = array![ @@ -176,6 +180,10 @@ fn decompose_two_qubit_product_gate( temp = unitary.dot(&temp); let mut l = temp.slice_mut(s![..;2, ..;2]).to_owned(); let det_l = det_one_qubit(l.view()); + assert!( + det_l.abs() < 0.9, + "decompose_two_qubit_product_gate: unable to decompose: detL < 0.9" + ); l.mapv_inplace(|x| x / det_l.sqrt()); let phase = det_l.arg() / 2.; (l, r, phase) From 207225fe9e5ed64b4ddb4cb7ddc6a2af62bab605 Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Sun, 10 Mar 2024 10:22:31 +0900 Subject: [PATCH 08/46] Use rng first trial from Python implementation To aid in debugging and rule out rng differences causing different results this commit switches the first iteration of the randomized loop to have hard coded values that are identical to what the rng in numpy was returning. It is very unlikely that this will have any impact because the specific random numbers used shouldn't matter, this is mostly to just rule it out as a possibility in debugging the remaining test failures. --- crates/accelerate/src/two_qubit_decompose.rs | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/crates/accelerate/src/two_qubit_decompose.rs b/crates/accelerate/src/two_qubit_decompose.rs index 89225bbbeb6e..debcf4f51254 100644 --- a/crates/accelerate/src/two_qubit_decompose.rs +++ b/crates/accelerate/src/two_qubit_decompose.rs @@ -459,9 +459,21 @@ impl TwoQubitWeylDecomposition { 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); + for i in 0..100 { + let rand_a: f64; + let rand_b: f64; + // For debugging the algorithm use the same RNG values from the + // previous Python implementation for the first random trial. + // In most cases this loop only executes a single iteration and + // using the same rng values rules out possible RNG differences + // as the root cause of a test failure + if i == 0 { + rand_a = 1.2602066112249388; + rand_b = 0.22317849046722027; + } else { + rand_a = state.sample(StandardNormal); + rand_b = state.sample(StandardNormal); + } let m2_real = m2.mapv(|val| rand_a * val.re + rand_b * val.im).to_owned(); let p_inner = m2_real .view() From a47938204392affa63d29014517c53cc8c0be467 Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Sun, 10 Mar 2024 10:33:51 +0900 Subject: [PATCH 09/46] Fix whitespace in error message --- crates/accelerate/src/two_qubit_decompose.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/crates/accelerate/src/two_qubit_decompose.rs b/crates/accelerate/src/two_qubit_decompose.rs index debcf4f51254..3cbcc5b60dcf 100644 --- a/crates/accelerate/src/two_qubit_decompose.rs +++ b/crates/accelerate/src/two_qubit_decompose.rs @@ -168,7 +168,7 @@ fn decompose_two_qubit_product_gate( } assert!( det_r.abs() < 0.1, - "decompose_two_qubit_product_gate: unable to decompose: detR < 0 .1" + "decompose_two_qubit_product_gate: unable to decompose: detR < 0.1" ); r.mapv_inplace(|x| x / det_r.sqrt()); let r_t_conj: Array2 = r.t().mapv(|x| x.conj()).to_owned(); From 2eb7f86ae5414e0c04de92867ab24a86fb3049b8 Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Mon, 11 Mar 2024 14:31:18 +0900 Subject: [PATCH 10/46] Fix assert check --- crates/accelerate/src/two_qubit_decompose.rs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/crates/accelerate/src/two_qubit_decompose.rs b/crates/accelerate/src/two_qubit_decompose.rs index 3cbcc5b60dcf..55d89f2f7577 100644 --- a/crates/accelerate/src/two_qubit_decompose.rs +++ b/crates/accelerate/src/two_qubit_decompose.rs @@ -167,7 +167,7 @@ fn decompose_two_qubit_product_gate( det_r = det_one_qubit(r.view()); } assert!( - det_r.abs() < 0.1, + det_r.abs() >= 0.1, "decompose_two_qubit_product_gate: unable to decompose: detR < 0.1" ); r.mapv_inplace(|x| x / det_r.sqrt()); @@ -181,7 +181,7 @@ fn decompose_two_qubit_product_gate( let mut l = temp.slice_mut(s![..;2, ..;2]).to_owned(); let det_l = det_one_qubit(l.view()); assert!( - det_l.abs() < 0.9, + det_l.abs() >= 0.9, "decompose_two_qubit_product_gate: unable to decompose: detL < 0.9" ); l.mapv_inplace(|x| x / det_l.sqrt()); From c245f2527ef0e861bf9db4bedcbcf9ea8c9887a4 Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Mon, 11 Mar 2024 15:17:50 +0900 Subject: [PATCH 11/46] Fix missing multiplier on specialized weyl gate --- crates/accelerate/src/two_qubit_decompose.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/crates/accelerate/src/two_qubit_decompose.rs b/crates/accelerate/src/two_qubit_decompose.rs index 55d89f2f7577..4cd179060af9 100644 --- a/crates/accelerate/src/two_qubit_decompose.rs +++ b/crates/accelerate/src/two_qubit_decompose.rs @@ -378,7 +378,7 @@ impl TwoQubitWeylDecomposition { 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])); + sequence.push(("rzz".to_string(), vec![(PI4 - self.c) * 2.], [0, 1])); self.global_phase += PI4; } Specializations::SWAPEquiv => { From 0a59e79f0da332beae5ee1ea3c7a995b010695a0 Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Mon, 11 Mar 2024 15:46:55 +0900 Subject: [PATCH 12/46] Fix various mistakes This commit fixes two fundamental issues in the code. The first is the rz and ry matrix were being incorrectly constructed for a given angle. This caused the specializations that were computing the 1q matrices in the decomposition based on a product with these gates' matrices to return invalid results. The second issue is for the MirrorControlledEquiv specialization had the angles backwards for computing the matrix of the rz gates used in the products for the 1q matrices: `K1l = K1l @ Rz(K1r)` and `K1r = K1r @ Rz(K1l)` not `K1l = K1l @ Rz(K1l)` and `K1r = K1r @ Rz(K1r)` This was a typo from the original transcription. --- crates/accelerate/src/two_qubit_decompose.rs | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/crates/accelerate/src/two_qubit_decompose.rs b/crates/accelerate/src/two_qubit_decompose.rs index 4cd179060af9..fd432fc43882 100644 --- a/crates/accelerate/src/two_qubit_decompose.rs +++ b/crates/accelerate/src/two_qubit_decompose.rs @@ -311,14 +311,14 @@ fn rx_matrix(theta: f64) -> Array2 { 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]] + let sin = Complex64::new(half_theta.sin(), 0.); + array![[cos, -sin], [sin, cos]] } fn rz_matrix(theta: f64) -> Array2 { let ilam2 = Complex64::new(0., 0.5 * theta); array![ - [-ilam2.exp(), Complex64::new(0., 0.)], + [(-ilam2).exp(), Complex64::new(0., 0.)], [Complex64::new(0., 0.), ilam2.exp()] ] } @@ -343,7 +343,7 @@ enum Specializations { SimabmbEquiv, } -#[derive(Clone)] +#[derive(Clone, Debug)] #[allow(non_snake_case)] #[pyclass(module = "qiskit._accelerate.two_qubit_decompose", subclass)] pub struct TwoQubitWeylDecomposition { @@ -605,7 +605,6 @@ impl TwoQubitWeylDecomposition { 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. @@ -650,7 +649,6 @@ impl TwoQubitWeylDecomposition { } } }; - let mut specialized: TwoQubitWeylDecomposition = match specialization { Specializations::IdEquiv => TwoQubitWeylDecomposition { a: 0., @@ -746,9 +744,9 @@ impl TwoQubitWeylDecomposition { Specializations::ControlledEquiv => { let default_euler_basis = "XYX"; let [k2ltheta, k2lphi, k2llambda, k2lphase] = - angles_from_unitary(K2l.view(), "XYX"); + angles_from_unitary(K2l.view(), default_euler_basis); let [k2rtheta, k2rphi, k2rlambda, k2rphase] = - angles_from_unitary(K2r.view(), "XYX"); + angles_from_unitary(K2r.view(), default_euler_basis); TwoQubitWeylDecomposition { a, b: 0., @@ -775,8 +773,8 @@ impl TwoQubitWeylDecomposition { b: PI4, c, global_phase: global_phase + k2lphase + k2rphase, - K1l: K1l.dot(&rz_matrix(k2lphi)), - K1r: K1r.dot(&rz_matrix(k2rphi)), + K1l: K1l.dot(&rz_matrix(k2rphi)), + K1r: K1r.dot(&rz_matrix(k2lphi)), K2l: ry_matrix(k2ltheta).dot(&rz_matrix(k2llambda)), K2r: ry_matrix(k2rtheta).dot(&rz_matrix(k2rlambda)), specialization: Specializations::MirrorControlledEquiv, From 38e7c62e9f37800ba45a034c8385fea0c0fc0df4 Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Mon, 11 Mar 2024 17:26:55 +0900 Subject: [PATCH 13/46] Add pickle serialization support for TwoQubitWeylDecomposition --- crates/accelerate/src/two_qubit_decompose.rs | 105 ++++++++++++++++++- 1 file changed, 103 insertions(+), 2 deletions(-) diff --git a/crates/accelerate/src/two_qubit_decompose.rs b/crates/accelerate/src/two_qubit_decompose.rs index fd432fc43882..3ed6db560ab1 100644 --- a/crates/accelerate/src/two_qubit_decompose.rs +++ b/crates/accelerate/src/two_qubit_decompose.rs @@ -326,7 +326,7 @@ fn rz_matrix(theta: f64) -> Array2 { const DEFAULT_FIDELITY: f64 = 1.0 - 1.0e-9; const C1_IM: Complex64 = Complex64::new(0.0, 1.0); -#[derive(Clone, Debug)] +#[derive(Clone, Debug, Copy)] #[pyclass(module = "qiskit._accelerate.two_qubit_decompose")] enum Specializations { General, @@ -415,13 +415,114 @@ const IPX: [[Complex64; 2]; 2] = [ #[pymethods] impl TwoQubitWeylDecomposition { + fn __getstate__(&self, py: Python) -> ([f64; 5], [PyObject; 5], u8, String, Option) { + let specialization = match self.specialization { + Specializations::General => 0, + Specializations::IdEquiv => 1, + Specializations::SWAPEquiv => 2, + Specializations::PartialSWAPEquiv => 3, + Specializations::PartialSWAPFlipEquiv => 4, + Specializations::ControlledEquiv => 5, + Specializations::MirrorControlledEquiv => 6, + Specializations::SimaabEquiv => 7, + Specializations::SimabbEquiv => 8, + Specializations::SimabmbEquiv => 9, + }; + ( + [ + self.a, + self.b, + self.c, + self.global_phase, + self.calculated_fidelity, + ], + [ + self.K1l.to_pyarray(py).into(), + self.K1r.to_pyarray(py).into(), + self.K2l.to_pyarray(py).into(), + self.K2r.to_pyarray(py).into(), + self.unitary_matrix.to_pyarray(py).into(), + ], + specialization, + self.default_euler_basis.clone(), + self.requested_fidelity, + ) + } + + fn __setstate__( + &mut self, + state: ( + [f64; 5], + [PyReadonlyArray2; 5], + u8, + String, + Option, + ), + ) { + self.a = state.0[0]; + self.b = state.0[1]; + self.c = state.0[2]; + self.global_phase = state.0[3]; + self.calculated_fidelity = state.0[4]; + self.K1l = state.1[0].as_array().to_owned(); + self.K1r = state.1[1].as_array().to_owned(); + self.K2l = state.1[2].as_array().to_owned(); + self.K2r = state.1[3].as_array().to_owned(); + self.unitary_matrix = state.1[4].as_array().to_owned(); + self.default_euler_basis = state.3; + self.requested_fidelity = state.4; + self.specialization = match state.2 { + 0 => Specializations::General, + 1 => Specializations::IdEquiv, + 2 => Specializations::SWAPEquiv, + 3 => Specializations::PartialSWAPEquiv, + 4 => Specializations::PartialSWAPFlipEquiv, + 5 => Specializations::ControlledEquiv, + 6 => Specializations::MirrorControlledEquiv, + 7 => Specializations::SimaabEquiv, + 8 => Specializations::SimabbEquiv, + 9 => Specializations::SimabmbEquiv, + _ => unreachable!("Invalid specialization value"), + }; + } + + fn __getnewargs__(&self, py: Python) -> (PyObject, Option, Option, bool) { + ( + self.unitary_matrix.to_pyarray(py).into(), + self.requested_fidelity, + None, + true, + ) + } + #[new] - #[pyo3(signature=(unitary_matrix, fidelity=DEFAULT_FIDELITY, specialization=None))] + #[pyo3(signature=(unitary_matrix, fidelity=DEFAULT_FIDELITY, specialization=None, _pickle_context=false))] fn new( unitary_matrix: PyReadonlyArray2, fidelity: Option, specialization: Option, + _pickle_context: bool, ) -> PyResult { + // If we're in a pickle context just make the closest to an empty + // object as we can with minimal allocations and effort. All the + // data will be filled in during deserialization from __setstate__. + if _pickle_context { + return Ok(TwoQubitWeylDecomposition { + a: 0., + b: 0., + c: 0., + global_phase: 0., + K1l: Array2::zeros((0, 0)), + K2l: Array2::zeros((0, 0)), + K1r: Array2::zeros((0, 0)), + K2r: Array2::zeros((0, 0)), + specialization: Specializations::General, + default_euler_basis: "ZYZ".to_string(), + requested_fidelity: fidelity, + calculated_fidelity: 0., + unitary_matrix: Array2::zeros((0, 0)), + }); + } let ipz: ArrayView2 = aview2(&IPZ); let ipy: ArrayView2 = aview2(&IPY); let ipx: ArrayView2 = aview2(&IPX); From 9ae72880d1496f9c8eec259175bfc4e795fcc37b Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Mon, 11 Mar 2024 18:29:31 +0900 Subject: [PATCH 14/46] Remove duplicate conjs increment --- crates/accelerate/src/two_qubit_decompose.rs | 1 - 1 file changed, 1 deletion(-) diff --git a/crates/accelerate/src/two_qubit_decompose.rs b/crates/accelerate/src/two_qubit_decompose.rs index 3ed6db560ab1..a7f93a68a3f6 100644 --- a/crates/accelerate/src/two_qubit_decompose.rs +++ b/crates/accelerate/src/two_qubit_decompose.rs @@ -675,7 +675,6 @@ impl TwoQubitWeylDecomposition { } if cs[1] > PI4 { cs[1] = PI2 - cs[1]; - conjs += 1; K1l = K1l.dot(&ipx); K2r = ipx.dot(&K2r); conjs += 1; From 9b58e341feabd143e5969ff51614575b97b3a834 Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Mon, 11 Mar 2024 18:39:32 +0900 Subject: [PATCH 15/46] Fix fSim abmb equivalent specialization --- crates/accelerate/src/two_qubit_decompose.rs | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/crates/accelerate/src/two_qubit_decompose.rs b/crates/accelerate/src/two_qubit_decompose.rs index a7f93a68a3f6..7a04a01f909d 100644 --- a/crates/accelerate/src/two_qubit_decompose.rs +++ b/crates/accelerate/src/two_qubit_decompose.rs @@ -924,14 +924,13 @@ impl TwoQubitWeylDecomposition { } } Specializations::SimabmbEquiv => { - // TwoQubitWeylfSimabmbEquiv let default_euler_basis = "XYX"; let [k2ltheta, k2lphi, k2llambda, k2lphase] = - angles_from_unitary(K2l.view(), "XYX"); + angles_from_unitary(K2l.view(), default_euler_basis); TwoQubitWeylDecomposition { a, b: (b - c) / 2., - c: (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), @@ -982,7 +981,12 @@ impl TwoQubitWeylDecomposition { 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")); + return Err(PyValueError::new_err(format!( + "Specialization: {:?} calculated fidelity: {} is worse than requested fidelity: {}", + specialized.specialization, + specialized.calculated_fidelity, + fid + ))); } } specialized.global_phase += tr.arg(); From 096550bb87b9e5926838b9c011a1a9d07c97442b Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Mon, 11 Mar 2024 19:28:56 +0900 Subject: [PATCH 16/46] Add specialized test class back --- crates/accelerate/src/two_qubit_decompose.rs | 44 ++-- qiskit/__init__.py | 1 + .../two_qubit/two_qubit_decompose.py | 10 +- test/python/synthesis/test_synthesis.py | 194 +++++++++++++++++- 4 files changed, 217 insertions(+), 32 deletions(-) diff --git a/crates/accelerate/src/two_qubit_decompose.rs b/crates/accelerate/src/two_qubit_decompose.rs index 7a04a01f909d..27671b1bdf0b 100644 --- a/crates/accelerate/src/two_qubit_decompose.rs +++ b/crates/accelerate/src/two_qubit_decompose.rs @@ -359,6 +359,7 @@ pub struct TwoQubitWeylDecomposition { K2l: Array2, K1r: Array2, K2r: Array2, + #[pyo3(get)] specialization: Specializations, default_euler_basis: String, #[pyo3(get)] @@ -374,16 +375,17 @@ impl TwoQubitWeylDecomposition { simplify: bool, sequence: &mut Vec<(String, Vec, [u8; 2])>, atol: f64, + global_phase: &mut f64, ) { match self.specialization { Specializations::MirrorControlledEquiv => { sequence.push(("swap".to_string(), Vec::new(), [0, 1])); sequence.push(("rzz".to_string(), vec![(PI4 - self.c) * 2.], [0, 1])); - self.global_phase += PI4; + *global_phase += PI4 } Specializations::SWAPEquiv => { sequence.push(("swap".to_string(), Vec::new(), [0, 1])); - self.global_phase -= 3. * PI / 4.; + *global_phase -= 3. * PI / 4. } _ => { if !simplify || self.a.abs() > atol { @@ -762,7 +764,7 @@ impl TwoQubitWeylDecomposition { specialization: Specializations::IdEquiv, default_euler_basis: default_euler_basis.to_string(), requested_fidelity: fidelity, - calculated_fidelity: 1.0, + calculated_fidelity: -1.0, unitary_matrix, }, Specializations::SWAPEquiv => { @@ -779,7 +781,7 @@ impl TwoQubitWeylDecomposition { specialization: Specializations::SWAPEquiv, default_euler_basis: default_euler_basis.to_string(), requested_fidelity: fidelity, - calculated_fidelity: 1.0, + calculated_fidelity: -1.0, unitary_matrix, } } else { @@ -796,7 +798,7 @@ impl TwoQubitWeylDecomposition { specialization: Specializations::SWAPEquiv, default_euler_basis: default_euler_basis.to_string(), requested_fidelity: fidelity, - calculated_fidelity: 1.0, + calculated_fidelity: -1.0, unitary_matrix, } } @@ -817,7 +819,7 @@ impl TwoQubitWeylDecomposition { specialization: Specializations::PartialSWAPEquiv, default_euler_basis: default_euler_basis.to_string(), requested_fidelity: fidelity, - calculated_fidelity: 1.0, + calculated_fidelity: -1.0, unitary_matrix, } } @@ -828,7 +830,7 @@ impl TwoQubitWeylDecomposition { TwoQubitWeylDecomposition { a: closest, b: closest, - c: -a, + c: -closest, global_phase, K1l: K1l.dot(&K2l), K1r: K1r.dot(&ipz).dot(&K2l).dot(&ipz), @@ -837,7 +839,7 @@ impl TwoQubitWeylDecomposition { specialization: Specializations::PartialSWAPFlipEquiv, default_euler_basis: default_euler_basis.to_string(), requested_fidelity: fidelity, - calculated_fidelity: 1.0, + calculated_fidelity: -1.0, unitary_matrix, } } @@ -859,7 +861,7 @@ impl TwoQubitWeylDecomposition { specialization: Specializations::ControlledEquiv, default_euler_basis: default_euler_basis.to_string(), requested_fidelity: fidelity, - calculated_fidelity: 1.0, + calculated_fidelity: -1.0, unitary_matrix, } } @@ -880,7 +882,7 @@ impl TwoQubitWeylDecomposition { specialization: Specializations::MirrorControlledEquiv, default_euler_basis: default_euler_basis.to_string(), requested_fidelity: fidelity, - calculated_fidelity: 1.0, + calculated_fidelity: -1.0, unitary_matrix, } } @@ -899,7 +901,7 @@ impl TwoQubitWeylDecomposition { specialization: Specializations::SimaabEquiv, default_euler_basis: default_euler_basis.to_string(), requested_fidelity: fidelity, - calculated_fidelity: 1.0, + calculated_fidelity: -1.0, unitary_matrix, } } @@ -914,12 +916,12 @@ impl TwoQubitWeylDecomposition { 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), + K2l: ry_matrix(k2ltheta).dot(&rx_matrix(k2llambda)), + K2r: rx_matrix(-k2lphi).dot(&K2r), specialization: Specializations::SimabbEquiv, default_euler_basis: default_euler_basis.to_string(), requested_fidelity: fidelity, - calculated_fidelity: 1.0, + calculated_fidelity: -1.0, unitary_matrix, } } @@ -932,14 +934,14 @@ impl TwoQubitWeylDecomposition { b: (b - c) / 2., c: -((b - c) / 2.), global_phase: global_phase + k2lphase, - K1l: K1l.dot(&rz_matrix(k2lphi)), + K1l: K1l.dot(&rx_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, + calculated_fidelity: -1.0, unitary_matrix, } } @@ -955,7 +957,7 @@ impl TwoQubitWeylDecomposition { specialization: Specializations::General, default_euler_basis: default_euler_basis.to_string(), requested_fidelity: fidelity, - calculated_fidelity: 1.0, + calculated_fidelity: -1.0, unitary_matrix, }, }; @@ -977,7 +979,6 @@ impl TwoQubitWeylDecomposition { 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 { @@ -1062,7 +1063,12 @@ impl TwoQubitWeylDecomposition { 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)); + self.weyl_gate( + simplify, + &mut gate_sequence, + atol.unwrap_or(DEFAULT_ATOL), + &mut global_phase, + ); let c1r = unitary_to_gate_sequence_inner( self.K1r.view(), &target_1q_basis_list, diff --git a/qiskit/__init__.py b/qiskit/__init__.py index 782d8aaabfbd..4d91f6d6e86a 100644 --- a/qiskit/__init__.py +++ b/qiskit/__init__.py @@ -70,6 +70,7 @@ sys.modules["qiskit._accelerate.convert_2q_block_matrix"] = ( qiskit._accelerate.convert_2q_block_matrix ) +sys.modules["qiskit._accelerate.two_qubit_decompose"] = qiskit._accelerate.two_qubit_decompose # qiskit errors operator from qiskit.exceptions import QiskitError, MissingOptionalLibraryError diff --git a/qiskit/synthesis/two_qubit/two_qubit_decompose.py b/qiskit/synthesis/two_qubit/two_qubit_decompose.py index 0ad159d4d567..bc10e2b55c45 100644 --- a/qiskit/synthesis/two_qubit/two_qubit_decompose.py +++ b/qiskit/synthesis/two_qubit/two_qubit_decompose.py @@ -186,6 +186,7 @@ def __repr__(self): + b64ascii + [ f"requested_fidelity={self.requested_fidelity},", + f"specialization={self.specialization}," f"calculated_fidelity={self.calculated_fidelity},", f"actual_fidelity={self.actual_fidelity()},", f"abc={(self.a, self.b, self.c)})", @@ -195,7 +196,12 @@ def __repr__(self): @classmethod def from_bytes( - cls, bytes_in: bytes, *, requested_fidelity: float, **kwargs + cls, + bytes_in: bytes, + *, + requested_fidelity: float, + specialization: two_qubit_decompose.Specializations | None, + **kwargs, ) -> "TwoQubitWeylDecomposition": """Decode bytes into :class:`.TwoQubitWeylDecomposition`.""" # Used by __repr__ @@ -203,7 +209,7 @@ def from_bytes( b64 = base64.decodebytes(bytes_in) with io.BytesIO(b64) as f: arr = np.load(f, allow_pickle=False) - return cls(arr, fidelity=requested_fidelity) + return cls(arr, fidelity=requested_fidelity, specialization=specialization) def __str__(self): pre = f"{self.__class__.__name__}(\n\t" diff --git a/test/python/synthesis/test_synthesis.py b/test/python/synthesis/test_synthesis.py index 76ac5606aaac..e609a7f63164 100644 --- a/test/python/synthesis/test_synthesis.py +++ b/test/python/synthesis/test_synthesis.py @@ -61,6 +61,7 @@ decompose_two_qubit_product_gate, TwoQubitDecomposeUpToDiagonal, ) +from qiskit._accelerate.two_qubit_decompose import Specializations from qiskit.synthesis.unitary import qsd from test import combine # pylint: disable=wrong-import-order from test import QiskitTestCase # pylint: disable=wrong-import-order @@ -202,7 +203,14 @@ def check_two_qubit_weyl_specialization( # Loop to check both for implicit and explicity specialization for decomposer in (TwoQubitWeylDecomposition, expected_specialization): - decomp = decomposer(target_unitary, fidelity=fidelity) + if isinstance(decomposer, TwoQubitWeylDecomposition): + decomp = decomposer(target_unitary, fidelity=fidelity) + decomp_name = decomp.specialization + else: + decomp = TwoQubitWeylDecomposition( + target_unitary, fidelity=None, specialization=expected_specialization + ) + decomp_name = expected_specialization self.assertRoundTrip(decomp) self.assertRoundTripPickle(decomp) self.assertEqual( @@ -210,25 +218,31 @@ def check_two_qubit_weyl_specialization( 0, "Incorrect saved unitary in the decomposition.", ) - self.assertIsInstance(decomp, expected_specialization, "Incorrect Weyl specialization.") + self.assertEqual( + decomp.specialization, expected_specialization, "Incorrect Weyl specialization." + ) circ = decomp.circuit(simplify=True) self.assertDictEqual( - dict(circ.count_ops()), expected_gates, f"Gate counts of {decomposer.__name__}" + dict(circ.count_ops()), expected_gates, f"Gate counts of {decomp_name}" ) actual_fid = decomp.actual_fidelity() self.assertAlmostEqual(decomp.calculated_fidelity, actual_fid, places=13) - self.assertGreaterEqual(actual_fid, fidelity, f"fidelity of {decomposer.__name__}") + self.assertGreaterEqual(actual_fid, fidelity, f"fidelity of {decomp_name}") actual_unitary = Operator(circ).data trace = np.trace(actual_unitary.T.conj() @ target_unitary) - self.assertAlmostEqual( - trace.imag, 0, places=13, msg=f"Real trace for {decomposer.__name__}" - ) - decomp2 = expected_specialization(target_unitary, fidelity=None) # Shouldn't raise + self.assertAlmostEqual(trace.imag, 0, places=13, msg=f"Real trace for {decomp_name}") + decomp2 = TwoQubitWeylDecomposition( + target_unitary, fidelity=None, specialization=expected_specialization + ) # Shouldn't raise self.assertRoundTrip(decomp2) self.assertRoundTripPickle(decomp2) - with self.assertRaises(QiskitError) as exc: - _ = expected_specialization(target_unitary, fidelity=1.0) - self.assertIn("worse than requested", exc.exception.message) + if expected_specialization != Specializations.General: + # TODO: When rust code emits QiskitError change assertion to match + with self.assertRaises(ValueError) as exc: + _ = TwoQubitWeylDecomposition( + target_unitary, fidelity=1.0, specialization=expected_specialization + ) + self.assertIn("worse than requested", str(exc.exception)) def check_exact_decomposition( self, target_unitary, decomposer, tolerance=1.0e-12, num_basis_uses=None @@ -785,6 +799,164 @@ 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, + Specializations.IdEquiv, + {"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, + Specializations.SWAPEquiv, + {"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, + Specializations.SWAPEquiv, + {"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, + Specializations.PartialSWAPEquiv, + {"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, + Specializations.PartialSWAPFlipEquiv, + {"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, + Specializations.SimaabEquiv, + {"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, + Specializations.SimabbEquiv, + {"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, + Specializations.SimabmbEquiv, + {"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, + Specializations.ControlledEquiv, + {"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, + Specializations.MirrorControlledEquiv, + {"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, + Specializations.General, + {"rz": 8, "ry": 4, "rxx": 1, "ryy": 1, "rzz": 1}, + ) + + @ddt class TestTwoQubitDecompose(CheckDecompositions): """Test TwoQubitBasisDecomposer() for exact/approx decompositions""" From 17e917f3d83a9f6952d7813509d3acd1fba2e485 Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Mon, 11 Mar 2024 20:24:39 +0900 Subject: [PATCH 17/46] Use QiskitError for backwards compatibility --- crates/accelerate/src/two_qubit_decompose.rs | 9 ++++++--- test/python/synthesis/test_synthesis.py | 5 ++--- 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/crates/accelerate/src/two_qubit_decompose.rs b/crates/accelerate/src/two_qubit_decompose.rs index 27671b1bdf0b..8acda292445f 100644 --- a/crates/accelerate/src/two_qubit_decompose.rs +++ b/crates/accelerate/src/two_qubit_decompose.rs @@ -20,7 +20,8 @@ use approx::abs_diff_eq; use num_complex::{Complex, Complex64, ComplexFloat}; -use pyo3::exceptions::{PyIndexError, PyTypeError, PyValueError}; +use pyo3::exceptions::PyIndexError; +use pyo3::import_exception; use pyo3::prelude::*; use pyo3::wrap_pyfunction; use pyo3::Python; @@ -613,7 +614,8 @@ impl TwoQubitWeylDecomposition { } } if !found { - return Err(PyTypeError::new_err(format!( + import_exception!(qiskit, QiskitError); + return Err(QiskitError::new_err(format!( "TwoQubitWeylDecomposition: failed to diagonalize M2. Please report this at https://github.com/Qiskit/qiskit-terra/issues/4159. Input: {:?}", unitary_matrix ))); } @@ -982,7 +984,8 @@ impl TwoQubitWeylDecomposition { 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(format!( + import_exception!(qiskit, QiskitError); + return Err(QiskitError::new_err(format!( "Specialization: {:?} calculated fidelity: {} is worse than requested fidelity: {}", specialized.specialization, specialized.calculated_fidelity, diff --git a/test/python/synthesis/test_synthesis.py b/test/python/synthesis/test_synthesis.py index e609a7f63164..520ba6310237 100644 --- a/test/python/synthesis/test_synthesis.py +++ b/test/python/synthesis/test_synthesis.py @@ -237,11 +237,10 @@ def check_two_qubit_weyl_specialization( self.assertRoundTrip(decomp2) self.assertRoundTripPickle(decomp2) if expected_specialization != Specializations.General: - # TODO: When rust code emits QiskitError change assertion to match - with self.assertRaises(ValueError) as exc: + with self.assertRaises(QiskitError) as exc: _ = TwoQubitWeylDecomposition( target_unitary, fidelity=1.0, specialization=expected_specialization - ) + ) self.assertIn("worse than requested", str(exc.exception)) def check_exact_decomposition( From d0863fb2286da71e14bcac7bb2b69c9e4b38cc3a Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Mon, 11 Mar 2024 20:42:15 +0900 Subject: [PATCH 18/46] Add release note --- .../rust-two-qubit-weyl-ec551f3f9c812124.yaml | 14 ++++++++++++++ 1 file changed, 14 insertions(+) create mode 100644 releasenotes/notes/rust-two-qubit-weyl-ec551f3f9c812124.yaml diff --git a/releasenotes/notes/rust-two-qubit-weyl-ec551f3f9c812124.yaml b/releasenotes/notes/rust-two-qubit-weyl-ec551f3f9c812124.yaml new file mode 100644 index 000000000000..5ea9d6012948 --- /dev/null +++ b/releasenotes/notes/rust-two-qubit-weyl-ec551f3f9c812124.yaml @@ -0,0 +1,14 @@ +--- +features_synthesis: + - | + The :class:`.TwoQubitWeylDecomposition` synthesis class has been rewritten + in Rust for better performance. +upgrade_synthesis: + - | + The :class:`.TwoQubitWeylDecomposition` no longer will self specialize into + a subclass on creation. This was an internal detail of the :class:`.TwoQubitWeylDecomposition` + previously, and was not a documented public behavior as all the subclasses behaved + the same and were only used for internal dispatch. However, as it was discoverable + behavior this release note is to document that this will no longer occur and all + instances of :class:`.TwoQubitWeylDecomposition` will be of the same type. There is no + change in behavior for public methods of the class. From 6b4501fef84932745e1db5f8249f0e124ecab149 Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Mon, 11 Mar 2024 20:45:26 +0900 Subject: [PATCH 19/46] Make explicit specialization private --- crates/accelerate/src/two_qubit_decompose.rs | 6 +++--- qiskit/synthesis/two_qubit/two_qubit_decompose.py | 10 +++++----- test/python/synthesis/test_synthesis.py | 6 +++--- 3 files changed, 11 insertions(+), 11 deletions(-) diff --git a/crates/accelerate/src/two_qubit_decompose.rs b/crates/accelerate/src/two_qubit_decompose.rs index 8acda292445f..e89ff4a9fbe4 100644 --- a/crates/accelerate/src/two_qubit_decompose.rs +++ b/crates/accelerate/src/two_qubit_decompose.rs @@ -499,11 +499,11 @@ impl TwoQubitWeylDecomposition { } #[new] - #[pyo3(signature=(unitary_matrix, fidelity=DEFAULT_FIDELITY, specialization=None, _pickle_context=false))] + #[pyo3(signature=(unitary_matrix, fidelity=DEFAULT_FIDELITY, _specialization=None, _pickle_context=false))] fn new( unitary_matrix: PyReadonlyArray2, fidelity: Option, - specialization: Option, + _specialization: Option, _pickle_context: bool, ) -> PyResult { // If we're in a pickle context just make the closest to an empty @@ -727,7 +727,7 @@ impl TwoQubitWeylDecomposition { 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 { + let specialization = match _specialization { Some(specialization) => specialization, None => { if is_close(0., 0., 0.) { diff --git a/qiskit/synthesis/two_qubit/two_qubit_decompose.py b/qiskit/synthesis/two_qubit/two_qubit_decompose.py index bc10e2b55c45..bb72b0067083 100644 --- a/qiskit/synthesis/two_qubit/two_qubit_decompose.py +++ b/qiskit/synthesis/two_qubit/two_qubit_decompose.py @@ -186,7 +186,7 @@ def __repr__(self): + b64ascii + [ f"requested_fidelity={self.requested_fidelity},", - f"specialization={self.specialization}," + f"_specialization={self.specialization}," f"calculated_fidelity={self.calculated_fidelity},", f"actual_fidelity={self.actual_fidelity()},", f"abc={(self.a, self.b, self.c)})", @@ -200,7 +200,7 @@ def from_bytes( bytes_in: bytes, *, requested_fidelity: float, - specialization: two_qubit_decompose.Specializations | None, + _specialization: two_qubit_decompose.Specializations | None, **kwargs, ) -> "TwoQubitWeylDecomposition": """Decode bytes into :class:`.TwoQubitWeylDecomposition`.""" @@ -209,7 +209,7 @@ def from_bytes( b64 = base64.decodebytes(bytes_in) with io.BytesIO(b64) as f: arr = np.load(f, allow_pickle=False) - return cls(arr, fidelity=requested_fidelity, specialization=specialization) + return cls(arr, fidelity=requested_fidelity, _specialization=_specialization) def __str__(self): pre = f"{self.__class__.__name__}(\n\t" @@ -248,7 +248,7 @@ def __init__(self, rxx_equivalent_gate: Type[Gate]): decomposer_rxx = TwoQubitWeylDecomposition( Operator(circ).data, fidelity=None, - specialization=two_qubit_decompose.Specializations.ControlledEquiv, + _specialization=two_qubit_decompose.Specializations.ControlledEquiv, ) circ = QuantumCircuit(2) @@ -256,7 +256,7 @@ def __init__(self, rxx_equivalent_gate: Type[Gate]): decomposer_equiv = TwoQubitWeylDecomposition( Operator(circ).data, fidelity=None, - specialization=two_qubit_decompose.Specializations.ControlledEquiv, + _specialization=two_qubit_decompose.Specializations.ControlledEquiv, ) scale = decomposer_rxx.a / decomposer_equiv.a diff --git a/test/python/synthesis/test_synthesis.py b/test/python/synthesis/test_synthesis.py index 520ba6310237..75ff43e784ae 100644 --- a/test/python/synthesis/test_synthesis.py +++ b/test/python/synthesis/test_synthesis.py @@ -208,7 +208,7 @@ def check_two_qubit_weyl_specialization( decomp_name = decomp.specialization else: decomp = TwoQubitWeylDecomposition( - target_unitary, fidelity=None, specialization=expected_specialization + target_unitary, fidelity=None, _specialization=expected_specialization ) decomp_name = expected_specialization self.assertRoundTrip(decomp) @@ -232,14 +232,14 @@ def check_two_qubit_weyl_specialization( trace = np.trace(actual_unitary.T.conj() @ target_unitary) self.assertAlmostEqual(trace.imag, 0, places=13, msg=f"Real trace for {decomp_name}") decomp2 = TwoQubitWeylDecomposition( - target_unitary, fidelity=None, specialization=expected_specialization + target_unitary, fidelity=None, _specialization=expected_specialization ) # Shouldn't raise self.assertRoundTrip(decomp2) self.assertRoundTripPickle(decomp2) if expected_specialization != Specializations.General: with self.assertRaises(QiskitError) as exc: _ = TwoQubitWeylDecomposition( - target_unitary, fidelity=1.0, specialization=expected_specialization + target_unitary, fidelity=1.0, _specialization=expected_specialization ) self.assertIn("worse than requested", str(exc.exception)) From b273c9398a05d3e9c760354adfc0eb80799953e1 Mon Sep 17 00:00:00 2001 From: Jake Lishman Date: Mon, 11 Mar 2024 18:54:29 +0000 Subject: [PATCH 20/46] Use explicit inheritance from general decomposition (#23) --- crates/accelerate/src/two_qubit_decompose.rs | 186 ++++++++----------- 1 file changed, 78 insertions(+), 108 deletions(-) diff --git a/crates/accelerate/src/two_qubit_decompose.rs b/crates/accelerate/src/two_qubit_decompose.rs index e89ff4a9fbe4..262812ce5265 100644 --- a/crates/accelerate/src/two_qubit_decompose.rs +++ b/crates/accelerate/src/two_qubit_decompose.rs @@ -753,215 +753,185 @@ impl TwoQubitWeylDecomposition { } } }; + let general = TwoQubitWeylDecomposition { + a, + b, + c, + global_phase, + K1l, + K1r, + K2l, + K2r, + specialization: Specializations::General, + default_euler_basis: default_euler_basis.to_owned(), + requested_fidelity: fidelity, + calculated_fidelity: -1.0, + unitary_matrix, + }; let mut specialized: TwoQubitWeylDecomposition = match specialization { Specializations::IdEquiv => TwoQubitWeylDecomposition { + specialization, a: 0., b: 0., c: 0., - global_phase, - K1l: K1l.dot(&K2l), - K1r: K1r.dot(&K2r), + K1l: general.K1l.dot(&general.K2l), + K1r: general.K1r.dot(&general.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, + ..general }, Specializations::SWAPEquiv => { if c > 0. { TwoQubitWeylDecomposition { + specialization, a: PI4, b: PI4, c: PI4, - global_phase, - K1l: K1l.dot(&K2r), - K1r: K1r.dot(&K2l), + K1l: general.K1l.dot(&general.K2r), + K1r: general.K1r.dot(&general.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, + ..general } } else { flipped_from_original = true; TwoQubitWeylDecomposition { + specialization, a: PI4, b: PI4, c: PI4, global_phase: global_phase + PI2, - K1l: K1l.dot(&ipz).dot(&K2r), - K1r: K1r.dot(&ipz).dot(&K2l), + K1l: general.K1l.dot(&ipz).dot(&general.K2r), + K1r: general.K1r.dot(&ipz).dot(&general.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, + ..general } } } 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()); + let mut k2l_dag = general.K2l.t().to_owned(); + k2l_dag.view_mut().mapv_inplace(|x| x.conj()); TwoQubitWeylDecomposition { + specialization, a: closest, b: closest, c: closest, - global_phase, - K1l: K1l.dot(&K2l), - K1r: K1r.dot(&K2l), - K2r: k2r_temp.dot(&K2r), + K1l: general.K1l.dot(&general.K2l), + K1r: general.K1r.dot(&general.K2l), + K2r: k2l_dag.dot(&general.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, + ..general } } 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()); + let mut k2l_dag = general.K2l.t().to_owned(); + k2l_dag.mapv_inplace(|x| x.conj()); TwoQubitWeylDecomposition { + specialization, a: closest, b: closest, c: -closest, - global_phase, - K1l: K1l.dot(&K2l), - K1r: K1r.dot(&ipz).dot(&K2l).dot(&ipz), - K2r: ipz.dot(&k2_temp).dot(&ipz).dot(&K2r), + K1l: general.K1l.dot(&general.K2l), + K1r: general.K1r.dot(&ipz).dot(&general.K2l).dot(&ipz), + K2r: ipz.dot(&k2l_dag).dot(&ipz).dot(&general.K2r), K2l: Array2::eye(2), - specialization: Specializations::PartialSWAPFlipEquiv, - default_euler_basis: default_euler_basis.to_string(), - requested_fidelity: fidelity, - calculated_fidelity: -1.0, - unitary_matrix, + ..general } } Specializations::ControlledEquiv => { - let default_euler_basis = "XYX"; + let euler_basis = "XYX".to_owned(); let [k2ltheta, k2lphi, k2llambda, k2lphase] = - angles_from_unitary(K2l.view(), default_euler_basis); + angles_from_unitary(general.K2l.view(), &euler_basis); let [k2rtheta, k2rphi, k2rlambda, k2rphase] = - angles_from_unitary(K2r.view(), default_euler_basis); + angles_from_unitary(general.K2r.view(), &euler_basis); TwoQubitWeylDecomposition { + specialization, a, b: 0., c: 0., global_phase: global_phase + k2lphase + k2rphase, - K1l: K1l.dot(&rx_matrix(k2lphi)), - K1r: K1r.dot(&rx_matrix(k2rphi)), + K1l: general.K1l.dot(&rx_matrix(k2lphi)), + K1r: general.K1r.dot(&rx_matrix(k2rphi)), K2l: ry_matrix(k2ltheta).dot(&rx_matrix(k2llambda)), K2r: ry_matrix(k2rtheta).dot(&rx_matrix(k2rlambda)), - specialization: Specializations::ControlledEquiv, - default_euler_basis: default_euler_basis.to_string(), - requested_fidelity: fidelity, - calculated_fidelity: -1.0, - unitary_matrix, + default_euler_basis: euler_basis, + ..general } } Specializations::MirrorControlledEquiv => { let [k2ltheta, k2lphi, k2llambda, k2lphase] = - angles_from_unitary(K2l.view(), "ZYZ"); + angles_from_unitary(general.K2l.view(), "ZYZ"); let [k2rtheta, k2rphi, k2rlambda, k2rphase] = - angles_from_unitary(K2r.view(), "ZYZ"); + angles_from_unitary(general.K2r.view(), "ZYZ"); TwoQubitWeylDecomposition { + specialization, a: PI4, b: PI4, c, global_phase: global_phase + k2lphase + k2rphase, - K1l: K1l.dot(&rz_matrix(k2rphi)), - K1r: K1r.dot(&rz_matrix(k2lphi)), + K1l: general.K1l.dot(&rz_matrix(k2rphi)), + K1r: general.K1r.dot(&rz_matrix(k2lphi)), K2l: ry_matrix(k2ltheta).dot(&rz_matrix(k2llambda)), K2r: ry_matrix(k2rtheta).dot(&rz_matrix(k2rlambda)), - specialization: Specializations::MirrorControlledEquiv, - default_euler_basis: default_euler_basis.to_string(), - requested_fidelity: fidelity, - calculated_fidelity: -1.0, - unitary_matrix, + ..general } } Specializations::SimaabEquiv => { let [k2ltheta, k2lphi, k2llambda, k2lphase] = - angles_from_unitary(K2l.view(), "ZYZ"); + angles_from_unitary(general.K2l.view(), "ZYZ"); TwoQubitWeylDecomposition { + specialization, 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)), + K1r: general.K1r.dot(&rz_matrix(k2lphi)), + K1l: general.K1l.dot(&rz_matrix(k2lphi)), K2l: ry_matrix(k2ltheta).dot(&rz_matrix(k2llambda)), - K2r: rz_matrix(-k2lphi).dot(&K2r), - specialization: Specializations::SimaabEquiv, - default_euler_basis: default_euler_basis.to_string(), - requested_fidelity: fidelity, - calculated_fidelity: -1.0, - unitary_matrix, + K2r: rz_matrix(-k2lphi).dot(&general.K2r), + ..general } } Specializations::SimabbEquiv => { - let default_euler_basis = "XYX"; + let euler_basis = "XYX".to_owned(); let [k2ltheta, k2lphi, k2llambda, k2lphase] = - angles_from_unitary(K2l.view(), "XYX"); + angles_from_unitary(general.K2l.view(), &euler_basis); TwoQubitWeylDecomposition { + specialization, a, b: (b + c) / 2., c: (b + c) / 2., global_phase: global_phase + k2lphase, - K1r: K1r.dot(&rx_matrix(k2lphi)), - K1l: K1l.dot(&rx_matrix(k2lphi)), + K1r: general.K1r.dot(&rx_matrix(k2lphi)), + K1l: general.K1l.dot(&rx_matrix(k2lphi)), K2l: ry_matrix(k2ltheta).dot(&rx_matrix(k2llambda)), - K2r: rx_matrix(-k2lphi).dot(&K2r), - specialization: Specializations::SimabbEquiv, - default_euler_basis: default_euler_basis.to_string(), - requested_fidelity: fidelity, - calculated_fidelity: -1.0, - unitary_matrix, + K2r: rx_matrix(-k2lphi).dot(&general.K2r), + default_euler_basis: euler_basis, + ..general } } Specializations::SimabmbEquiv => { - let default_euler_basis = "XYX"; + let euler_basis = "XYX".to_owned(); let [k2ltheta, k2lphi, k2llambda, k2lphase] = - angles_from_unitary(K2l.view(), default_euler_basis); + angles_from_unitary(general.K2l.view(), &euler_basis); TwoQubitWeylDecomposition { + specialization, a, b: (b - c) / 2., c: -((b - c) / 2.), global_phase: global_phase + k2lphase, - K1l: K1l.dot(&rx_matrix(k2lphi)), - K1r: K1r.dot(&ipz).dot(&rx_matrix(k2lphi)).dot(&ipz), + K1l: general.K1l.dot(&rx_matrix(k2lphi)), + K1r: general.K1r.dot(&ipz).dot(&rx_matrix(k2lphi)).dot(&ipz), K2l: ry_matrix(k2ltheta).dot(&rx_matrix(k2llambda)), - K2r: ipz.dot(&rx_matrix(-k2lphi)).dot(&ipz).dot(&K2r), - specialization: Specializations::SimabmbEquiv, - default_euler_basis: default_euler_basis.to_string(), - requested_fidelity: fidelity, - calculated_fidelity: -1.0, - unitary_matrix, + K2r: ipz.dot(&rx_matrix(-k2lphi)).dot(&ipz).dot(&general.K2r), + default_euler_basis: euler_basis, + ..general } } - 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, - }, + Specializations::General => general, }; let tr = if flipped_from_original { From 3bc353ff5ecdf5278774393c4c3679c93dfcbd2b Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Tue, 12 Mar 2024 04:36:01 +0900 Subject: [PATCH 21/46] Apply suggestions from code review --- crates/accelerate/src/two_qubit_decompose.rs | 12 +----------- 1 file changed, 1 insertion(+), 11 deletions(-) diff --git a/crates/accelerate/src/two_qubit_decompose.rs b/crates/accelerate/src/two_qubit_decompose.rs index 262812ce5265..fe770512150f 100644 --- a/crates/accelerate/src/two_qubit_decompose.rs +++ b/crates/accelerate/src/two_qubit_decompose.rs @@ -172,7 +172,7 @@ fn decompose_two_qubit_product_gate( "decompose_two_qubit_product_gate: unable to decompose: detR < 0.1" ); r.mapv_inplace(|x| x / det_r.sqrt()); - let r_t_conj: Array2 = r.t().mapv(|x| x.conj()).to_owned(); + let r_t_conj: Array2 = r.t().mapv(|x| x.conj()); let eye = array![ [Complex64::new(1., 0.), Complex64::new(0., 0.)], [Complex64::new(0., 0.), Complex64::new(1., 0.)], @@ -587,16 +587,6 @@ impl TwoQubitWeylDecomposition { .into_ndarray() .mapv(Complex64::from) .to_owned(); - // Uncomment this to use numpy for eigh instead of faer (useful if needed to compare) - // 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_inner = p_raw - // .extract::>()? - // .as_array() - // .mapv(Complex64::from); let d_inner = p_inner.t().dot(&m2).dot(&p_inner).diag().to_owned(); let mut diag_d: Array2 = Array2::zeros((4, 4)); diag_d From aff637eb0d830b26d16c978a50e86d52eb06bf2c Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Tue, 12 Mar 2024 05:03:08 +0900 Subject: [PATCH 22/46] Use smallvec for circuit sequence params and qubits --- Cargo.lock | 1 + crates/accelerate/Cargo.toml | 2 +- .../src/euler_one_qubit_decomposer.rs | 62 ++++++++++--------- crates/accelerate/src/two_qubit_decompose.rs | 33 +++++----- .../two_qubit/two_qubit_decompose.py | 6 +- 5 files changed, 56 insertions(+), 48 deletions(-) diff --git a/Cargo.lock b/Cargo.lock index 3bc457d6f65c..d095cadf4650 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -1089,6 +1089,7 @@ dependencies = [ "pyo3-build-config", "pyo3-ffi", "pyo3-macros", + "smallvec 1.13.1", "unindent", ] diff --git a/crates/accelerate/Cargo.toml b/crates/accelerate/Cargo.toml index 10e2dee6f031..fa855aede219 100644 --- a/crates/accelerate/Cargo.toml +++ b/crates/accelerate/Cargo.toml @@ -33,7 +33,7 @@ features = ["union"] [dependencies.pyo3] workspace = true -features = ["hashbrown", "indexmap", "num-complex", "num-bigint"] +features = ["hashbrown", "indexmap", "num-complex", "num-bigint", "smallvec"] [dependencies.ndarray] version = "^0.15.6" diff --git a/crates/accelerate/src/euler_one_qubit_decomposer.rs b/crates/accelerate/src/euler_one_qubit_decomposer.rs index dbfa0535da4e..98aff4bc27c4 100644 --- a/crates/accelerate/src/euler_one_qubit_decomposer.rs +++ b/crates/accelerate/src/euler_one_qubit_decomposer.rs @@ -14,6 +14,7 @@ use hashbrown::HashMap; use num_complex::{Complex64, ComplexFloat}; +use smallvec::{smallvec, SmallVec}; use std::cmp::Ordering; use std::f64::consts::PI; @@ -61,11 +62,13 @@ impl OneQubitGateErrorMap { #[pyclass(sequence)] pub struct OneQubitGateSequence { - pub gates: Vec<(String, Vec)>, + pub gates: Vec<(String, SmallVec<[f64; 3]>)>, #[pyo3(get)] pub global_phase: f64, } +type OneQubitGateSequenceState = (Vec<(String, SmallVec<[f64; 3]>)>, f64); + #[pymethods] impl OneQubitGateSequence { #[new] @@ -75,11 +78,11 @@ impl OneQubitGateSequence { global_phase: 0., } } - fn __getstate__(&self) -> (Vec<(String, Vec)>, f64) { + fn __getstate__(&self) -> OneQubitGateSequenceState { (self.gates.clone(), self.global_phase) } - fn __setstate__(&mut self, state: (Vec<(String, Vec)>, f64)) { + fn __setstate__(&mut self, state: OneQubitGateSequenceState) { self.gates = state.0; self.global_phase = state.1; } @@ -93,7 +96,7 @@ impl OneQubitGateSequence { SliceOrInt::Slice(slc) => { let len = self.gates.len().try_into().unwrap(); let indices = slc.indices(len)?; - let mut out_vec: Vec<(String, Vec)> = Vec::new(); + let mut out_vec: Vec<(String, SmallVec<[f64; 3]>)> = Vec::new(); // Start and stop will always be positive the slice api converts // negatives to the index for example: // list(range(5))[-1:-3:-1] @@ -145,7 +148,7 @@ fn circuit_kak( let mut lam = lam; let mut theta = theta; let mut phi = phi; - let mut circuit: Vec<(String, Vec)> = Vec::with_capacity(3); + let mut circuit: Vec<(String, SmallVec<[f64; 3]>)> = Vec::with_capacity(3); let mut atol = match atol { Some(atol) => atol, None => DEFAULT_ATOL, @@ -161,7 +164,7 @@ fn circuit_kak( // slippage coming from _mod_2pi injecting multiples of 2pi. lam = mod_2pi(lam, atol); if lam.abs() > atol { - circuit.push((String::from(k_gate), vec![lam])); + circuit.push((String::from(k_gate), smallvec![lam])); global_phase += lam / 2.; } return OneQubitGateSequence { @@ -182,13 +185,13 @@ fn circuit_kak( lam = mod_2pi(lam, atol); if lam.abs() > atol { global_phase += lam / 2.; - circuit.push((String::from(k_gate), vec![lam])); + circuit.push((String::from(k_gate), smallvec![lam])); } - circuit.push((String::from(a_gate), vec![theta])); + circuit.push((String::from(a_gate), smallvec![theta])); phi = mod_2pi(phi, atol); if phi.abs() > atol { global_phase += phi / 2.; - circuit.push((String::from(k_gate), vec![phi])); + circuit.push((String::from(k_gate), smallvec![phi])); } OneQubitGateSequence { gates: circuit, @@ -212,7 +215,7 @@ fn circuit_u3( let phi = mod_2pi(phi, atol); let lam = mod_2pi(lam, atol); if !simplify || theta.abs() > atol || phi.abs() > atol || lam.abs() > atol { - circuit.push((String::from("u3"), vec![theta, phi, lam])); + circuit.push((String::from("u3"), smallvec![theta, phi, lam])); } OneQubitGateSequence { gates: circuit, @@ -239,17 +242,17 @@ fn circuit_u321( if theta.abs() < atol { let tot = mod_2pi(phi + lam, atol); if tot.abs() > atol { - circuit.push((String::from("u1"), vec![tot])); + circuit.push((String::from("u1"), smallvec![tot])); } } else if (theta - PI / 2.).abs() < atol { circuit.push(( String::from("u2"), - vec![mod_2pi(phi, atol), mod_2pi(lam, atol)], + smallvec![mod_2pi(phi, atol), mod_2pi(lam, atol)], )); } else { circuit.push(( String::from("u3"), - vec![theta, mod_2pi(phi, atol), mod_2pi(lam, atol)], + smallvec![theta, mod_2pi(phi, atol), mod_2pi(lam, atol)], )); } OneQubitGateSequence { @@ -277,7 +280,7 @@ fn circuit_u( let phi = mod_2pi(phi, atol); let lam = mod_2pi(lam, atol); if theta.abs() > atol || phi.abs() > atol || lam.abs() > atol { - circuit.push((String::from("u"), vec![theta, phi, lam])); + circuit.push((String::from("u"), smallvec![theta, phi, lam])); } OneQubitGateSequence { gates: circuit, @@ -379,19 +382,22 @@ fn circuit_rr( if mod_2pi((phi + lam) / 2., atol).abs() < atol { // This can be expressed as a single R gate if theta.abs() > atol { - circuit.push((String::from("r"), vec![theta, mod_2pi(PI / 2. + phi, atol)])); + circuit.push(( + String::from("r"), + smallvec![theta, mod_2pi(PI / 2. + phi, atol)], + )); } } else { // General case: use two R gates if (theta - PI).abs() > atol { circuit.push(( String::from("r"), - vec![theta - PI, mod_2pi(PI / 2. - lam, atol)], + smallvec![theta - PI, mod_2pi(PI / 2. - lam, atol)], )); } circuit.push(( String::from("r"), - vec![PI, mod_2pi(0.5 * (phi - lam + PI), atol)], + smallvec![PI, mod_2pi(0.5 * (phi - lam + PI), atol)], )); } @@ -430,11 +436,11 @@ pub fn generate_circuit( let fnz = |circuit: &mut OneQubitGateSequence, phi: f64| { let phi = mod_2pi(phi, inner_atol); if phi.abs() > inner_atol { - circuit.gates.push((String::from("p"), vec![phi])); + circuit.gates.push((String::from("p"), smallvec![phi])); } }; let fnx = |circuit: &mut OneQubitGateSequence| { - circuit.gates.push((String::from("sx"), Vec::new())); + circuit.gates.push((String::from("sx"), SmallVec::new())); }; circuit_psx_gen( @@ -460,12 +466,12 @@ pub fn generate_circuit( let fnz = |circuit: &mut OneQubitGateSequence, phi: f64| { let phi = mod_2pi(phi, inner_atol); if phi.abs() > inner_atol { - circuit.gates.push((String::from("rz"), vec![phi])); + circuit.gates.push((String::from("rz"), smallvec![phi])); circuit.global_phase += phi / 2.; } }; let fnx = |circuit: &mut OneQubitGateSequence| { - circuit.gates.push((String::from("sx"), Vec::new())); + circuit.gates.push((String::from("sx"), SmallVec::new())); }; circuit_psx_gen( theta, @@ -490,12 +496,12 @@ pub fn generate_circuit( let fnz = |circuit: &mut OneQubitGateSequence, phi: f64| { let phi = mod_2pi(phi, inner_atol); if phi.abs() > inner_atol { - circuit.gates.push((String::from("u1"), vec![phi])); + circuit.gates.push((String::from("u1"), smallvec![phi])); } }; let fnx = |circuit: &mut OneQubitGateSequence| { circuit.global_phase += PI / 4.; - circuit.gates.push((String::from("rx"), vec![PI / 2.])); + circuit.gates.push((String::from("rx"), smallvec![PI / 2.])); }; circuit_psx_gen( theta, @@ -520,15 +526,15 @@ pub fn generate_circuit( let fnz = |circuit: &mut OneQubitGateSequence, phi: f64| { let phi = mod_2pi(phi, inner_atol); if phi.abs() > inner_atol { - circuit.gates.push((String::from("rz"), vec![phi])); + circuit.gates.push((String::from("rz"), smallvec![phi])); circuit.global_phase += phi / 2.; } }; let fnx = |circuit: &mut OneQubitGateSequence| { - circuit.gates.push((String::from("sx"), Vec::new())); + circuit.gates.push((String::from("sx"), SmallVec::new())); }; let fnxpi = |circuit: &mut OneQubitGateSequence| { - circuit.gates.push((String::from("x"), Vec::new())); + circuit.gates.push((String::from("x"), SmallVec::new())); }; circuit_psx_gen( theta, @@ -592,7 +598,7 @@ fn compare_error_fn( } fn compute_error( - gates: &[(String, Vec)], + gates: &[(String, SmallVec<[f64; 3]>)], error_map: Option<&OneQubitGateErrorMap>, qubit: usize, ) -> (f64, usize) { @@ -622,7 +628,7 @@ pub fn compute_error_one_qubit_sequence( #[inline] #[pyfunction] pub fn compute_error_list( - circuit: Vec<(String, Vec)>, + circuit: Vec<(String, SmallVec<[f64; 3]>)>, qubit: usize, error_map: Option<&OneQubitGateErrorMap>, ) -> (f64, usize) { diff --git a/crates/accelerate/src/two_qubit_decompose.rs b/crates/accelerate/src/two_qubit_decompose.rs index fe770512150f..104bcaf15e6f 100644 --- a/crates/accelerate/src/two_qubit_decompose.rs +++ b/crates/accelerate/src/two_qubit_decompose.rs @@ -25,6 +25,7 @@ use pyo3::import_exception; use pyo3::prelude::*; use pyo3::wrap_pyfunction; use pyo3::Python; +use smallvec::{smallvec, SmallVec}; use std::f64::consts::PI; use faer::IntoFaerComplex; @@ -374,29 +375,33 @@ impl TwoQubitWeylDecomposition { fn weyl_gate( &mut self, simplify: bool, - sequence: &mut Vec<(String, Vec, [u8; 2])>, + sequence: &mut TwoQubitSequenceVec, atol: f64, global_phase: &mut f64, ) { match self.specialization { Specializations::MirrorControlledEquiv => { - sequence.push(("swap".to_string(), Vec::new(), [0, 1])); - sequence.push(("rzz".to_string(), vec![(PI4 - self.c) * 2.], [0, 1])); + sequence.push(("swap".to_string(), SmallVec::new(), smallvec![0, 1])); + sequence.push(( + "rzz".to_string(), + smallvec![(PI4 - self.c) * 2.], + smallvec![0, 1], + )); *global_phase += PI4 } Specializations::SWAPEquiv => { - sequence.push(("swap".to_string(), Vec::new(), [0, 1])); + sequence.push(("swap".to_string(), SmallVec::new(), smallvec![0, 1])); *global_phase -= 3. * PI / 4. } _ => { if !simplify || self.a.abs() > atol { - sequence.push(("rxx".to_string(), vec![-self.a * 2.], [0, 1])); + sequence.push(("rxx".to_string(), smallvec![-self.a * 2.], smallvec![0, 1])); } if !simplify || self.b.abs() > atol { - sequence.push(("ryy".to_string(), vec![-self.b * 2.], [0, 1])); + sequence.push(("ryy".to_string(), smallvec![-self.b * 2.], smallvec![0, 1])); } if !simplify || self.c.abs() > atol { - sequence.push(("rzz".to_string(), vec![-self.c * 2.], [0, 1])); + sequence.push(("rzz".to_string(), smallvec![-self.c * 2.], smallvec![0, 1])); } } } @@ -614,7 +619,7 @@ impl TwoQubitWeylDecomposition { let mut cs: Array1 = (0..3) .map(|i| ((d[i] + d[3]) / 2.0).rem_euclid(TWO_PI)) .collect(); - let cstemp: Vec = cs + let cstemp: SmallVec<[f64; 3]> = cs .iter() .map(|x| x.rem_euclid(PI2)) .map(|x| x.min(PI2 - x)) @@ -1010,7 +1015,7 @@ impl TwoQubitWeylDecomposition { ) .unwrap(); for gate in c2r.gates { - gate_sequence.push((gate.0, gate.1, [0, 0])) + gate_sequence.push((gate.0, gate.1, smallvec![0])) } global_phase += c2r.global_phase; let c2l = unitary_to_gate_sequence_inner( @@ -1023,7 +1028,7 @@ impl TwoQubitWeylDecomposition { ) .unwrap(); for gate in c2l.gates { - gate_sequence.push((gate.0, gate.1, [1, 1])) + gate_sequence.push((gate.0, gate.1, smallvec![1])) } global_phase += c2l.global_phase; self.weyl_gate( @@ -1042,7 +1047,7 @@ impl TwoQubitWeylDecomposition { ) .unwrap(); for gate in c1r.gates { - gate_sequence.push((gate.0, gate.1, [0, 0])) + gate_sequence.push((gate.0, gate.1, smallvec![0])) } global_phase += c2r.global_phase; let c1l = unitary_to_gate_sequence_inner( @@ -1055,7 +1060,7 @@ impl TwoQubitWeylDecomposition { ) .unwrap(); for gate in c1l.gates { - gate_sequence.push((gate.0, gate.1, [1, 1])) + gate_sequence.push((gate.0, gate.1, smallvec![1])) } TwoQubitGateSequence { gates: gate_sequence, @@ -1064,7 +1069,7 @@ impl TwoQubitWeylDecomposition { } } -type TwoQubitSequenceVec = Vec<(String, Vec, [u8; 2])>; +type TwoQubitSequenceVec = Vec<(String, SmallVec<[f64; 3]>, SmallVec<[u8; 2]>)>; #[pyclass(sequence)] pub struct TwoQubitGateSequence { @@ -1101,7 +1106,7 @@ impl TwoQubitGateSequence { 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(); + let mut out_vec: TwoQubitSequenceVec = Vec::new(); // Start and stop will always be positive the slice api converts // negatives to the index for example: // list(range(5))[-1:-3:-1] diff --git a/qiskit/synthesis/two_qubit/two_qubit_decompose.py b/qiskit/synthesis/two_qubit/two_qubit_decompose.py index bb72b0067083..cbb12bbb5c62 100644 --- a/qiskit/synthesis/two_qubit/two_qubit_decompose.py +++ b/qiskit/synthesis/two_qubit/two_qubit_decompose.py @@ -157,11 +157,7 @@ def circuit( circuit_sequence = super().circuit(euler_basis=euler_basis, simplify=simplify, atol=atol) 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) + getattr(circ, name)(*params, *qubits) return circ def actual_fidelity(self, **kwargs) -> float: From 817cf4bc3d08484a3ef704271765115738900d1b Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Tue, 12 Mar 2024 05:20:22 +0900 Subject: [PATCH 23/46] Use const 2x2 identity array where possible --- crates/accelerate/src/two_qubit_decompose.rs | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/crates/accelerate/src/two_qubit_decompose.rs b/crates/accelerate/src/two_qubit_decompose.rs index 104bcaf15e6f..b5a752621634 100644 --- a/crates/accelerate/src/two_qubit_decompose.rs +++ b/crates/accelerate/src/two_qubit_decompose.rs @@ -56,6 +56,17 @@ const TWO_PI: f64 = 2.0 * PI; const C1: c64 = c64 { re: 1.0, im: 0.0 }; +const ONE_QUBIT_IDENTITY: [[Complex64; 2]; 2] = [ + [ + Complex64::new(1., 0.), + Complex64::new(0., 0.), + ], + [ + Complex64::new(0., 0.), + Complex64::new(1., 0.), + ], +]; + const B_NON_NORMALIZED: [[Complex64; 4]; 4] = [ [ Complex64::new(1.0, 0.), @@ -174,10 +185,7 @@ fn decompose_two_qubit_product_gate( ); r.mapv_inplace(|x| x / det_r.sqrt()); let r_t_conj: Array2 = r.t().mapv(|x| x.conj()); - let eye = array![ - [Complex64::new(1., 0.), Complex64::new(0., 0.)], - [Complex64::new(0., 0.), Complex64::new(1., 0.)], - ]; + let eye = aview2(&ONE_QUBIT_IDENTITY); let mut temp = kron(&eye, &r_t_conj); temp = unitary.dot(&temp); let mut l = temp.slice_mut(s![..;2, ..;2]).to_owned(); From f2e7b6d669348a2cc86f267270ae6cf8058e13a3 Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Tue, 12 Mar 2024 05:20:53 +0900 Subject: [PATCH 24/46] circuit() and weyl_gate don't need a mutable self --- crates/accelerate/src/two_qubit_decompose.rs | 17 +++++------------ 1 file changed, 5 insertions(+), 12 deletions(-) diff --git a/crates/accelerate/src/two_qubit_decompose.rs b/crates/accelerate/src/two_qubit_decompose.rs index b5a752621634..d108cfc40629 100644 --- a/crates/accelerate/src/two_qubit_decompose.rs +++ b/crates/accelerate/src/two_qubit_decompose.rs @@ -57,14 +57,8 @@ const TWO_PI: f64 = 2.0 * PI; const C1: c64 = c64 { re: 1.0, im: 0.0 }; const ONE_QUBIT_IDENTITY: [[Complex64; 2]; 2] = [ - [ - Complex64::new(1., 0.), - Complex64::new(0., 0.), - ], - [ - Complex64::new(0., 0.), - Complex64::new(1., 0.), - ], + [Complex64::new(1., 0.), Complex64::new(0., 0.)], + [Complex64::new(0., 0.), Complex64::new(1., 0.)], ]; const B_NON_NORMALIZED: [[Complex64; 4]; 4] = [ @@ -381,7 +375,7 @@ pub struct TwoQubitWeylDecomposition { impl TwoQubitWeylDecomposition { fn weyl_gate( - &mut self, + &self, simplify: bool, sequence: &mut TwoQubitSequenceVec, atol: f64, @@ -1001,13 +995,12 @@ impl TwoQubitWeylDecomposition { #[pyo3(signature = (euler_basis=None, simplify=false, atol=None))] fn circuit( - &mut self, + &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 euler_basis: &str = euler_basis.unwrap_or(&self.default_euler_basis); let target_1q_basis_list: Vec<&str> = vec![euler_basis]; let mut gate_sequence = Vec::new(); From 2c5d9634d5a643f37539a821115ed727504d1a8f Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Tue, 12 Mar 2024 05:25:54 +0900 Subject: [PATCH 25/46] Remove unnecessary inline annotations --- crates/accelerate/src/euler_one_qubit_decomposer.rs | 5 ----- 1 file changed, 5 deletions(-) diff --git a/crates/accelerate/src/euler_one_qubit_decomposer.rs b/crates/accelerate/src/euler_one_qubit_decomposer.rs index 98aff4bc27c4..460474cea113 100644 --- a/crates/accelerate/src/euler_one_qubit_decomposer.rs +++ b/crates/accelerate/src/euler_one_qubit_decomposer.rs @@ -615,7 +615,6 @@ fn compute_error( } } -#[inline] #[pyfunction] pub fn compute_error_one_qubit_sequence( circuit: &OneQubitGateSequence, @@ -625,7 +624,6 @@ pub fn compute_error_one_qubit_sequence( compute_error(&circuit.gates, error_map, qubit) } -#[inline] #[pyfunction] pub fn compute_error_list( circuit: Vec<(String, SmallVec<[f64; 3]>)>, @@ -723,7 +721,6 @@ 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(); @@ -739,7 +736,6 @@ 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(); @@ -754,7 +750,6 @@ 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(); From 99db4bf986e6ab528ccd27562a6741a26b8a1dca Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Tue, 12 Mar 2024 05:29:34 +0900 Subject: [PATCH 26/46] Rename Specializations enum Specialization --- crates/accelerate/src/two_qubit_decompose.rs | 98 +++++++++---------- .../two_qubit/two_qubit_decompose.py | 6 +- test/python/synthesis/test_synthesis.py | 26 ++--- 3 files changed, 65 insertions(+), 65 deletions(-) diff --git a/crates/accelerate/src/two_qubit_decompose.rs b/crates/accelerate/src/two_qubit_decompose.rs index d108cfc40629..c134cab8ed6b 100644 --- a/crates/accelerate/src/two_qubit_decompose.rs +++ b/crates/accelerate/src/two_qubit_decompose.rs @@ -332,7 +332,7 @@ const C1_IM: Complex64 = Complex64::new(0.0, 1.0); #[derive(Clone, Debug, Copy)] #[pyclass(module = "qiskit._accelerate.two_qubit_decompose")] -enum Specializations { +enum Specialization { General, IdEquiv, SWAPEquiv, @@ -364,7 +364,7 @@ pub struct TwoQubitWeylDecomposition { K1r: Array2, K2r: Array2, #[pyo3(get)] - specialization: Specializations, + specialization: Specialization, default_euler_basis: String, #[pyo3(get)] requested_fidelity: Option, @@ -382,7 +382,7 @@ impl TwoQubitWeylDecomposition { global_phase: &mut f64, ) { match self.specialization { - Specializations::MirrorControlledEquiv => { + Specialization::MirrorControlledEquiv => { sequence.push(("swap".to_string(), SmallVec::new(), smallvec![0, 1])); sequence.push(( "rzz".to_string(), @@ -391,7 +391,7 @@ impl TwoQubitWeylDecomposition { )); *global_phase += PI4 } - Specializations::SWAPEquiv => { + Specialization::SWAPEquiv => { sequence.push(("swap".to_string(), SmallVec::new(), smallvec![0, 1])); *global_phase -= 3. * PI / 4. } @@ -427,16 +427,16 @@ const IPX: [[Complex64; 2]; 2] = [ impl TwoQubitWeylDecomposition { fn __getstate__(&self, py: Python) -> ([f64; 5], [PyObject; 5], u8, String, Option) { let specialization = match self.specialization { - Specializations::General => 0, - Specializations::IdEquiv => 1, - Specializations::SWAPEquiv => 2, - Specializations::PartialSWAPEquiv => 3, - Specializations::PartialSWAPFlipEquiv => 4, - Specializations::ControlledEquiv => 5, - Specializations::MirrorControlledEquiv => 6, - Specializations::SimaabEquiv => 7, - Specializations::SimabbEquiv => 8, - Specializations::SimabmbEquiv => 9, + Specialization::General => 0, + Specialization::IdEquiv => 1, + Specialization::SWAPEquiv => 2, + Specialization::PartialSWAPEquiv => 3, + Specialization::PartialSWAPFlipEquiv => 4, + Specialization::ControlledEquiv => 5, + Specialization::MirrorControlledEquiv => 6, + Specialization::SimaabEquiv => 7, + Specialization::SimabbEquiv => 8, + Specialization::SimabmbEquiv => 9, }; ( [ @@ -482,21 +482,21 @@ impl TwoQubitWeylDecomposition { self.default_euler_basis = state.3; self.requested_fidelity = state.4; self.specialization = match state.2 { - 0 => Specializations::General, - 1 => Specializations::IdEquiv, - 2 => Specializations::SWAPEquiv, - 3 => Specializations::PartialSWAPEquiv, - 4 => Specializations::PartialSWAPFlipEquiv, - 5 => Specializations::ControlledEquiv, - 6 => Specializations::MirrorControlledEquiv, - 7 => Specializations::SimaabEquiv, - 8 => Specializations::SimabbEquiv, - 9 => Specializations::SimabmbEquiv, + 0 => Specialization::General, + 1 => Specialization::IdEquiv, + 2 => Specialization::SWAPEquiv, + 3 => Specialization::PartialSWAPEquiv, + 4 => Specialization::PartialSWAPFlipEquiv, + 5 => Specialization::ControlledEquiv, + 6 => Specialization::MirrorControlledEquiv, + 7 => Specialization::SimaabEquiv, + 8 => Specialization::SimabbEquiv, + 9 => Specialization::SimabmbEquiv, _ => unreachable!("Invalid specialization value"), }; } - fn __getnewargs__(&self, py: Python) -> (PyObject, Option, Option, bool) { + fn __getnewargs__(&self, py: Python) -> (PyObject, Option, Option, bool) { ( self.unitary_matrix.to_pyarray(py).into(), self.requested_fidelity, @@ -510,7 +510,7 @@ impl TwoQubitWeylDecomposition { fn new( unitary_matrix: PyReadonlyArray2, fidelity: Option, - _specialization: Option, + _specialization: Option, _pickle_context: bool, ) -> PyResult { // If we're in a pickle context just make the closest to an empty @@ -526,7 +526,7 @@ impl TwoQubitWeylDecomposition { K2l: Array2::zeros((0, 0)), K1r: Array2::zeros((0, 0)), K2r: Array2::zeros((0, 0)), - specialization: Specializations::General, + specialization: Specialization::General, default_euler_basis: "ZYZ".to_string(), requested_fidelity: fidelity, calculated_fidelity: 0., @@ -728,25 +728,25 @@ impl TwoQubitWeylDecomposition { Some(specialization) => specialization, None => { if is_close(0., 0., 0.) { - Specializations::IdEquiv + Specialization::IdEquiv } else if is_close(PI4, PI4, PI4) || is_close(PI4, PI4, -PI4) { - Specializations::SWAPEquiv + Specialization::SWAPEquiv } else if is_close(closest_abc, closest_abc, closest_abc) { - Specializations::PartialSWAPEquiv + Specialization::PartialSWAPEquiv } else if is_close(closest_ab_minus_c, closest_ab_minus_c, -closest_ab_minus_c) { - Specializations::PartialSWAPFlipEquiv + Specialization::PartialSWAPFlipEquiv } else if is_close(a, 0., 0.) { - Specializations::ControlledEquiv + Specialization::ControlledEquiv } else if is_close(PI4, PI4, c) { - Specializations::MirrorControlledEquiv + Specialization::MirrorControlledEquiv } else if is_close((a + b) / 2., (a + b) / 2., c) { - Specializations::SimaabEquiv + Specialization::SimaabEquiv } else if is_close(a, (b + c) / 2., (b + c) / 2.) { - Specializations::SimabbEquiv + Specialization::SimabbEquiv } else if is_close(a, (b - c) / 2., (c - b) / 2.) { - Specializations::SimabmbEquiv + Specialization::SimabmbEquiv } else { - Specializations::General + Specialization::General } } }; @@ -759,14 +759,14 @@ impl TwoQubitWeylDecomposition { K1r, K2l, K2r, - specialization: Specializations::General, + specialization: Specialization::General, default_euler_basis: default_euler_basis.to_owned(), requested_fidelity: fidelity, calculated_fidelity: -1.0, unitary_matrix, }; let mut specialized: TwoQubitWeylDecomposition = match specialization { - Specializations::IdEquiv => TwoQubitWeylDecomposition { + Specialization::IdEquiv => TwoQubitWeylDecomposition { specialization, a: 0., b: 0., @@ -777,7 +777,7 @@ impl TwoQubitWeylDecomposition { K2r: Array2::eye(2), ..general }, - Specializations::SWAPEquiv => { + Specialization::SWAPEquiv => { if c > 0. { TwoQubitWeylDecomposition { specialization, @@ -806,7 +806,7 @@ impl TwoQubitWeylDecomposition { } } } - Specializations::PartialSWAPEquiv => { + Specialization::PartialSWAPEquiv => { let closest = closest_partial_swap(a, b, c); let mut k2l_dag = general.K2l.t().to_owned(); k2l_dag.view_mut().mapv_inplace(|x| x.conj()); @@ -822,7 +822,7 @@ impl TwoQubitWeylDecomposition { ..general } } - Specializations::PartialSWAPFlipEquiv => { + Specialization::PartialSWAPFlipEquiv => { let closest = closest_partial_swap(a, b, -c); let mut k2l_dag = general.K2l.t().to_owned(); k2l_dag.mapv_inplace(|x| x.conj()); @@ -838,7 +838,7 @@ impl TwoQubitWeylDecomposition { ..general } } - Specializations::ControlledEquiv => { + Specialization::ControlledEquiv => { let euler_basis = "XYX".to_owned(); let [k2ltheta, k2lphi, k2llambda, k2lphase] = angles_from_unitary(general.K2l.view(), &euler_basis); @@ -858,7 +858,7 @@ impl TwoQubitWeylDecomposition { ..general } } - Specializations::MirrorControlledEquiv => { + Specialization::MirrorControlledEquiv => { let [k2ltheta, k2lphi, k2llambda, k2lphase] = angles_from_unitary(general.K2l.view(), "ZYZ"); let [k2rtheta, k2rphi, k2rlambda, k2rphase] = @@ -876,7 +876,7 @@ impl TwoQubitWeylDecomposition { ..general } } - Specializations::SimaabEquiv => { + Specialization::SimaabEquiv => { let [k2ltheta, k2lphi, k2llambda, k2lphase] = angles_from_unitary(general.K2l.view(), "ZYZ"); TwoQubitWeylDecomposition { @@ -892,7 +892,7 @@ impl TwoQubitWeylDecomposition { ..general } } - Specializations::SimabbEquiv => { + Specialization::SimabbEquiv => { let euler_basis = "XYX".to_owned(); let [k2ltheta, k2lphi, k2llambda, k2lphase] = angles_from_unitary(general.K2l.view(), &euler_basis); @@ -910,7 +910,7 @@ impl TwoQubitWeylDecomposition { ..general } } - Specializations::SimabmbEquiv => { + Specialization::SimabmbEquiv => { let euler_basis = "XYX".to_owned(); let [k2ltheta, k2lphi, k2llambda, k2lphase] = angles_from_unitary(general.K2l.view(), &euler_basis); @@ -928,7 +928,7 @@ impl TwoQubitWeylDecomposition { ..general } } - Specializations::General => general, + Specialization::General => general, }; let tr = if flipped_from_original { @@ -1151,6 +1151,6 @@ 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::()?; + m.add_class::()?; Ok(()) } diff --git a/qiskit/synthesis/two_qubit/two_qubit_decompose.py b/qiskit/synthesis/two_qubit/two_qubit_decompose.py index cbb12bbb5c62..22b4f1e4efa3 100644 --- a/qiskit/synthesis/two_qubit/two_qubit_decompose.py +++ b/qiskit/synthesis/two_qubit/two_qubit_decompose.py @@ -196,7 +196,7 @@ def from_bytes( bytes_in: bytes, *, requested_fidelity: float, - _specialization: two_qubit_decompose.Specializations | None, + _specialization: two_qubit_decompose.Specialization | None, **kwargs, ) -> "TwoQubitWeylDecomposition": """Decode bytes into :class:`.TwoQubitWeylDecomposition`.""" @@ -244,7 +244,7 @@ def __init__(self, rxx_equivalent_gate: Type[Gate]): decomposer_rxx = TwoQubitWeylDecomposition( Operator(circ).data, fidelity=None, - _specialization=two_qubit_decompose.Specializations.ControlledEquiv, + _specialization=two_qubit_decompose.Specialization.ControlledEquiv, ) circ = QuantumCircuit(2) @@ -252,7 +252,7 @@ def __init__(self, rxx_equivalent_gate: Type[Gate]): decomposer_equiv = TwoQubitWeylDecomposition( Operator(circ).data, fidelity=None, - _specialization=two_qubit_decompose.Specializations.ControlledEquiv, + _specialization=two_qubit_decompose.Specialization.ControlledEquiv, ) scale = decomposer_rxx.a / decomposer_equiv.a diff --git a/test/python/synthesis/test_synthesis.py b/test/python/synthesis/test_synthesis.py index 75ff43e784ae..bea64f9449f2 100644 --- a/test/python/synthesis/test_synthesis.py +++ b/test/python/synthesis/test_synthesis.py @@ -61,7 +61,7 @@ decompose_two_qubit_product_gate, TwoQubitDecomposeUpToDiagonal, ) -from qiskit._accelerate.two_qubit_decompose import Specializations +from qiskit._accelerate.two_qubit_decompose import Specialization from qiskit.synthesis.unitary import qsd from test import combine # pylint: disable=wrong-import-order from test import QiskitTestCase # pylint: disable=wrong-import-order @@ -236,7 +236,7 @@ def check_two_qubit_weyl_specialization( ) # Shouldn't raise self.assertRoundTrip(decomp2) self.assertRoundTripPickle(decomp2) - if expected_specialization != Specializations.General: + if expected_specialization != Specialization.General: with self.assertRaises(QiskitError) as exc: _ = TwoQubitWeylDecomposition( target_unitary, fidelity=1.0, _specialization=expected_specialization @@ -811,7 +811,7 @@ def test_weyl_specialize_id(self): self.check_two_qubit_weyl_specialization( k1 @ Ud(a + da, b + db, c + dc) @ k2, 0.999, - Specializations.IdEquiv, + Specialization.IdEquiv, {"rz": 4, "ry": 2}, ) @@ -825,7 +825,7 @@ def test_weyl_specialize_swap(self): self.check_two_qubit_weyl_specialization( k1 @ Ud(a + da, b + db, c + dc) @ k2, 0.999, - Specializations.SWAPEquiv, + Specialization.SWAPEquiv, {"rz": 4, "ry": 2, "swap": 1}, ) @@ -839,7 +839,7 @@ def test_weyl_specialize_flip_swap(self): self.check_two_qubit_weyl_specialization( k1 @ Ud(a + da, b + db, c + dc) @ k2, 0.999, - Specializations.SWAPEquiv, + Specialization.SWAPEquiv, {"rz": 4, "ry": 2, "swap": 1}, ) @@ -853,7 +853,7 @@ def test_weyl_specialize_pswap(self, theta=0.123): self.check_two_qubit_weyl_specialization( k1 @ Ud(a + da, b + db, c + dc) @ k2, 0.999, - Specializations.PartialSWAPEquiv, + Specialization.PartialSWAPEquiv, {"rz": 6, "ry": 3, "rxx": 1, "ryy": 1, "rzz": 1}, ) @@ -867,7 +867,7 @@ def test_weyl_specialize_flip_pswap(self, theta=0.123): self.check_two_qubit_weyl_specialization( k1 @ Ud(a + da, b + db, c + dc) @ k2, 0.999, - Specializations.PartialSWAPFlipEquiv, + Specialization.PartialSWAPFlipEquiv, {"rz": 6, "ry": 3, "rxx": 1, "ryy": 1, "rzz": 1}, ) @@ -881,7 +881,7 @@ def test_weyl_specialize_fsim_aab(self, aaa=0.456, bbb=0.132): self.check_two_qubit_weyl_specialization( k1 @ Ud(a + da, b + db, c + dc) @ k2, 0.999, - Specializations.SimaabEquiv, + Specialization.SimaabEquiv, {"rz": 7, "ry": 4, "rxx": 1, "ryy": 1, "rzz": 1}, ) @@ -895,7 +895,7 @@ def test_weyl_specialize_fsim_abb(self, aaa=0.456, bbb=0.132): self.check_two_qubit_weyl_specialization( k1 @ Ud(a + da, b + db, c + dc) @ k2, 0.999, - Specializations.SimabbEquiv, + Specialization.SimabbEquiv, {"rx": 7, "ry": 4, "rxx": 1, "ryy": 1, "rzz": 1}, ) @@ -909,7 +909,7 @@ def test_weyl_specialize_fsim_abmb(self, aaa=0.456, bbb=0.132): self.check_two_qubit_weyl_specialization( k1 @ Ud(a + da, b + db, c + dc) @ k2, 0.999, - Specializations.SimabmbEquiv, + Specialization.SimabmbEquiv, {"rx": 7, "ry": 4, "rxx": 1, "ryy": 1, "rzz": 1}, ) @@ -923,7 +923,7 @@ def test_weyl_specialize_ctrl(self, aaa=0.456): self.check_two_qubit_weyl_specialization( k1 @ Ud(a + da, b + db, c + dc) @ k2, 0.999, - Specializations.ControlledEquiv, + Specialization.ControlledEquiv, {"rx": 6, "ry": 4, "rxx": 1}, ) @@ -937,7 +937,7 @@ def test_weyl_specialize_mirror_ctrl(self, aaa=-0.456): self.check_two_qubit_weyl_specialization( k1 @ Ud(a + da, b + db, c + dc) @ k2, 0.999, - Specializations.MirrorControlledEquiv, + Specialization.MirrorControlledEquiv, {"rz": 6, "ry": 4, "rzz": 1, "swap": 1}, ) @@ -951,7 +951,7 @@ def test_weyl_specialize_general(self, aaa=0.456, bbb=0.345, ccc=0.123): self.check_two_qubit_weyl_specialization( k1 @ Ud(a + da, b + db, c + dc) @ k2, 0.999, - Specializations.General, + Specialization.General, {"rz": 8, "ry": 4, "rxx": 1, "ryy": 1, "rzz": 1}, ) From 1a8be8549588fdbeab259e6db5070e016439e086 Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Tue, 12 Mar 2024 05:39:14 +0900 Subject: [PATCH 27/46] Add back specialize() method and deprecate --- qiskit/synthesis/two_qubit/two_qubit_decompose.py | 10 ++++++++++ .../rust-two-qubit-weyl-ec551f3f9c812124.yaml | 15 ++++++++++++++- 2 files changed, 24 insertions(+), 1 deletion(-) diff --git a/qiskit/synthesis/two_qubit/two_qubit_decompose.py b/qiskit/synthesis/two_qubit/two_qubit_decompose.py index 22b4f1e4efa3..ebc1def81784 100644 --- a/qiskit/synthesis/two_qubit/two_qubit_decompose.py +++ b/qiskit/synthesis/two_qubit/two_qubit_decompose.py @@ -43,6 +43,7 @@ OneQubitEulerDecomposer, DEFAULT_ATOL, ) +from qiskit.utils.deprecation import deprecate_func from qiskit._accelerate import two_qubit_decompose logger = logging.getLogger(__name__) @@ -150,6 +151,15 @@ class TwoQubitWeylDecomposition(two_qubit_decompose.TwoQubitWeylDecomposition): requested_fidelity: Optional[float] # None means no automatic specialization calculated_fidelity: float # Fidelity after specialization + @deprecate_func(since="1.1.0", removal_timeline="in the 2.0.0 release") + def specialize(self): + """Make changes to the decomposition to comply with any specializations. + + This method will always raise a ``NotImplementedError`` because + there are no specializations to comply with in the current implementation. + """ + raise NotImplementedError + def circuit( self, *, euler_basis: str | None = None, simplify: bool = False, atol: float = DEFAULT_ATOL ) -> QuantumCircuit: diff --git a/releasenotes/notes/rust-two-qubit-weyl-ec551f3f9c812124.yaml b/releasenotes/notes/rust-two-qubit-weyl-ec551f3f9c812124.yaml index 5ea9d6012948..57379286f82f 100644 --- a/releasenotes/notes/rust-two-qubit-weyl-ec551f3f9c812124.yaml +++ b/releasenotes/notes/rust-two-qubit-weyl-ec551f3f9c812124.yaml @@ -5,10 +5,23 @@ features_synthesis: in Rust for better performance. upgrade_synthesis: - | - The :class:`.TwoQubitWeylDecomposition` no longer will self specialize into + The :class:`.TwoQubitWeylDecomposition` no longer will self-specialize into a subclass on creation. This was an internal detail of the :class:`.TwoQubitWeylDecomposition` previously, and was not a documented public behavior as all the subclasses behaved the same and were only used for internal dispatch. However, as it was discoverable behavior this release note is to document that this will no longer occur and all instances of :class:`.TwoQubitWeylDecomposition` will be of the same type. There is no change in behavior for public methods of the class. + +deprecations_synthesis: + - | + The :meth:`.TwoQubitWeylDecomposition.specialize` method is now deprecated + and will be removed in the Qiskit 2.0.0 release. This method never had + a public purpose and was unsafe for an end user to call as it would + mutate the calculated decomposition in the object and produce invalid + fields in the object. It was only used internally to construct a new + :class:`.TwoQubitWeylDecomposition` object. Despite this it was still a + documented part of the public API for the class and is now being + deprecated without any potential replacement. This release it always will + raise a ``NotImplementedError` when called because the specialization + subclassing has been removed as part of the Rust rewrite of the class. From a2fcc58a9edd2c1b1253d88a33d8054d0eca47f7 Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Tue, 12 Mar 2024 06:11:02 +0900 Subject: [PATCH 28/46] Reorganize Python/Rust split to wrap rust instead of subclass This commit reworks the handoff between python and rust space for the TwoQubitWeylDecomposition class. Previously the Python side TwoQubitWeylDecomposition was a subclass of the Rust struct/pyclass of the same name. This was originally done to deduplicate the the getters for the attributes. However, because of the overhead associated with some of the rust getters and needing to do some import normalization to a numpy array the deduplication wasn't worth the cost. --- .../two_qubit/two_qubit_decompose.py | 31 +++++++++++++++++-- test/python/synthesis/test_synthesis.py | 4 ++- 2 files changed, 31 insertions(+), 4 deletions(-) diff --git a/qiskit/synthesis/two_qubit/two_qubit_decompose.py b/qiskit/synthesis/two_qubit/two_qubit_decompose.py index ebc1def81784..5a1ca636e9bc 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(two_qubit_decompose.TwoQubitWeylDecomposition): +class TwoQubitWeylDecomposition: r"""Two-qubit Weyl decomposition. Decompose two-qubit unitary @@ -151,6 +151,29 @@ class TwoQubitWeylDecomposition(two_qubit_decompose.TwoQubitWeylDecomposition): requested_fidelity: Optional[float] # None means no automatic specialization calculated_fidelity: float # Fidelity after specialization + def __init__( + self, + unitary_matrix: np.ndarray, + fidelity: float | None = 1.0 - 1.0e-9, + *, + _specialization: two_qubit_decompose.Specialization | None = None, + ): + unitary_matrix = np.asarray(unitary_matrix, dtype=complex) + self._inner_decomposition = two_qubit_decompose.TwoQubitWeylDecomposition( + unitary_matrix, fidelity=fidelity, _specialization=_specialization + ) + self.a = self._inner_decomposition.a + self.b = self._inner_decomposition.b + self.c = self._inner_decomposition.c + self.global_phase = self._inner_decomposition.global_phase + self.K1l = self._inner_decomposition.K1l + self.K1r = self._inner_decomposition.K1r + self.K2l = self._inner_decomposition.K2l + self.K2r = self._inner_decomposition.K2r + self.unitary_matrix = unitary_matrix + self.requested_fidelity = fidelity + self.calculated_fidelity = self._inner_decomposition.calculated_fidelity + @deprecate_func(since="1.1.0", removal_timeline="in the 2.0.0 release") def specialize(self): """Make changes to the decomposition to comply with any specializations. @@ -164,7 +187,9 @@ def circuit( self, *, euler_basis: str | None = None, simplify: bool = False, atol: float = DEFAULT_ATOL ) -> QuantumCircuit: """Returns Weyl decomposition in circuit form.""" - circuit_sequence = super().circuit(euler_basis=euler_basis, simplify=simplify, atol=atol) + circuit_sequence = self._inner_decomposition.circuit( + euler_basis=euler_basis, simplify=simplify, atol=atol + ) circ = QuantumCircuit(2, global_phase=circuit_sequence.global_phase) for name, params, qubits in circuit_sequence: getattr(circ, name)(*params, *qubits) @@ -192,7 +217,7 @@ def __repr__(self): + b64ascii + [ f"requested_fidelity={self.requested_fidelity},", - f"_specialization={self.specialization}," + f"_specialization={self._inner_decomposition.specialization}," f"calculated_fidelity={self.calculated_fidelity},", f"actual_fidelity={self.actual_fidelity()},", f"abc={(self.a, self.b, self.c)})", diff --git a/test/python/synthesis/test_synthesis.py b/test/python/synthesis/test_synthesis.py index bea64f9449f2..0e891e4b0392 100644 --- a/test/python/synthesis/test_synthesis.py +++ b/test/python/synthesis/test_synthesis.py @@ -219,7 +219,9 @@ def check_two_qubit_weyl_specialization( "Incorrect saved unitary in the decomposition.", ) self.assertEqual( - decomp.specialization, expected_specialization, "Incorrect Weyl specialization." + decomp._inner_decomposition.specialization, + expected_specialization, + "Incorrect Weyl specialization.", ) circ = decomp.circuit(simplify=True) self.assertDictEqual( From e3c18ae23a32f2efa669fd90d0953454195f3a94 Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Tue, 12 Mar 2024 06:51:03 +0900 Subject: [PATCH 29/46] Remove unecessary allocations --- crates/accelerate/src/two_qubit_decompose.rs | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/crates/accelerate/src/two_qubit_decompose.rs b/crates/accelerate/src/two_qubit_decompose.rs index c134cab8ed6b..0ff25f583a2c 100644 --- a/crates/accelerate/src/two_qubit_decompose.rs +++ b/crates/accelerate/src/two_qubit_decompose.rs @@ -182,7 +182,7 @@ fn decompose_two_qubit_product_gate( let eye = aview2(&ONE_QUBIT_IDENTITY); let mut temp = kron(&eye, &r_t_conj); temp = unitary.dot(&temp); - let mut l = temp.slice_mut(s![..;2, ..;2]).to_owned(); + let mut l = temp.slice(s![..;2, ..;2]).to_owned(); let det_l = det_one_qubit(l.view()); assert!( det_l.abs() >= 0.9, @@ -585,15 +585,14 @@ impl TwoQubitWeylDecomposition { rand_a = state.sample(StandardNormal); rand_b = state.sample(StandardNormal); } - let m2_real = m2.mapv(|val| rand_a * val.re + rand_b * val.im).to_owned(); + let m2_real = m2.mapv(|val| rand_a * val.re + rand_b * val.im); let p_inner = m2_real .view() .into_faer() .selfadjoint_eigendecomposition(Lower) .u() .into_ndarray() - .mapv(Complex64::from) - .to_owned(); + .mapv(Complex64::from); let d_inner = p_inner.t().dot(&m2).dot(&p_inner).diag().to_owned(); let mut diag_d: Array2 = Array2::zeros((4, 4)); diag_d @@ -602,7 +601,7 @@ impl TwoQubitWeylDecomposition { .enumerate() .for_each(|(index, x)| *x = d_inner[index]); - let compare = p_inner.dot(&diag_d).dot(&p_inner.t()).to_owned(); + let compare = p_inner.dot(&diag_d).dot(&p_inner.t()); found = abs_diff_eq!(compare.view(), m2, epsilon = 1.0e-13); if found { p = p_inner; @@ -618,7 +617,7 @@ impl TwoQubitWeylDecomposition { } let mut d = -d.map(|x| x.arg() / 2.); d[3] = -d[0] - d[1] - d[2]; - let mut cs: Array1 = (0..3) + let mut cs: SmallVec<[f64; 3]> = (0..3) .map(|i| ((d[i] + d[3]) / 2.0).rem_euclid(TWO_PI)) .collect(); let cstemp: SmallVec<[f64; 3]> = cs From 376b2e7d9aa3746b6e028697b4b3e84a62f8bf08 Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Tue, 12 Mar 2024 06:55:06 +0900 Subject: [PATCH 30/46] Rename DEFAULT_ATOL to ANGLE_ZERO_EPSILON --- .../src/euler_one_qubit_decomposer.rs | 22 +++++++++---------- crates/accelerate/src/two_qubit_decompose.rs | 4 ++-- 2 files changed, 13 insertions(+), 13 deletions(-) diff --git a/crates/accelerate/src/euler_one_qubit_decomposer.rs b/crates/accelerate/src/euler_one_qubit_decomposer.rs index 460474cea113..ba1dfe2a711c 100644 --- a/crates/accelerate/src/euler_one_qubit_decomposer.rs +++ b/crates/accelerate/src/euler_one_qubit_decomposer.rs @@ -28,7 +28,7 @@ use numpy::PyReadonlyArray2; use crate::utils::SliceOrInt; -pub const DEFAULT_ATOL: f64 = 1e-12; +pub const ANGLE_ZERO_EPSILON: f64 = 1e-12; #[pyclass(module = "qiskit._accelerate.euler_one_qubit_decomposer")] pub struct OneQubitGateErrorMap { @@ -151,7 +151,7 @@ fn circuit_kak( let mut circuit: Vec<(String, SmallVec<[f64; 3]>)> = Vec::with_capacity(3); let mut atol = match atol { Some(atol) => atol, - None => DEFAULT_ATOL, + None => ANGLE_ZERO_EPSILON, }; if !simplify { atol = -1.0; @@ -210,7 +210,7 @@ fn circuit_u3( let mut circuit = Vec::new(); let atol = match atol { Some(atol) => atol, - None => DEFAULT_ATOL, + None => ANGLE_ZERO_EPSILON, }; let phi = mod_2pi(phi, atol); let lam = mod_2pi(lam, atol); @@ -234,7 +234,7 @@ fn circuit_u321( let mut circuit = Vec::new(); let mut atol = match atol { Some(atol) => atol, - None => DEFAULT_ATOL, + None => ANGLE_ZERO_EPSILON, }; if !simplify { atol = -1.0; @@ -272,7 +272,7 @@ fn circuit_u( let mut circuit = Vec::new(); let mut atol = match atol { Some(atol) => atol, - None => DEFAULT_ATOL, + None => ANGLE_ZERO_EPSILON, }; if !simplify { atol = -1.0; @@ -313,7 +313,7 @@ where }; let mut atol = match atol { Some(atol) => atol, - None => DEFAULT_ATOL, + None => ANGLE_ZERO_EPSILON, }; if !simplify { atol = -1.0; @@ -373,7 +373,7 @@ fn circuit_rr( let mut circuit = Vec::new(); let mut atol = match atol { Some(atol) => atol, - None => DEFAULT_ATOL, + None => ANGLE_ZERO_EPSILON, }; if !simplify { atol = -1.0; @@ -428,7 +428,7 @@ pub fn generate_circuit( "PSX" => { let mut inner_atol = match atol { Some(atol) => atol, - None => DEFAULT_ATOL, + None => ANGLE_ZERO_EPSILON, }; if !simplify { inner_atol = -1.0; @@ -458,7 +458,7 @@ pub fn generate_circuit( "ZSX" => { let mut inner_atol = match atol { Some(atol) => atol, - None => DEFAULT_ATOL, + None => ANGLE_ZERO_EPSILON, }; if !simplify { inner_atol = -1.0; @@ -488,7 +488,7 @@ pub fn generate_circuit( "U1X" => { let mut inner_atol = match atol { Some(atol) => atol, - None => DEFAULT_ATOL, + None => ANGLE_ZERO_EPSILON, }; if !simplify { inner_atol = -1.0; @@ -518,7 +518,7 @@ pub fn generate_circuit( "ZSXX" => { let mut inner_atol = match atol { Some(atol) => atol, - None => DEFAULT_ATOL, + None => ANGLE_ZERO_EPSILON, }; if !simplify { inner_atol = -1.0; diff --git a/crates/accelerate/src/two_qubit_decompose.rs b/crates/accelerate/src/two_qubit_decompose.rs index 0ff25f583a2c..0b7ca0e68ca1 100644 --- a/crates/accelerate/src/two_qubit_decompose.rs +++ b/crates/accelerate/src/two_qubit_decompose.rs @@ -41,7 +41,7 @@ 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, + angles_from_unitary, det_one_qubit, unitary_to_gate_sequence_inner, ANGLE_ZERO_EPSILON, }; use crate::utils; @@ -1034,7 +1034,7 @@ impl TwoQubitWeylDecomposition { self.weyl_gate( simplify, &mut gate_sequence, - atol.unwrap_or(DEFAULT_ATOL), + atol.unwrap_or(ANGLE_ZERO_EPSILON), &mut global_phase, ); let c1r = unitary_to_gate_sequence_inner( From 89070e445e4cf60c3e825634ddc139d9c97b5a8d Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Tue, 12 Mar 2024 06:58:53 +0900 Subject: [PATCH 31/46] Stop obsessing over -0 --- crates/accelerate/src/two_qubit_decompose.rs | 14 ++++---------- 1 file changed, 4 insertions(+), 10 deletions(-) diff --git a/crates/accelerate/src/two_qubit_decompose.rs b/crates/accelerate/src/two_qubit_decompose.rs index 0b7ca0e68ca1..439da22d888b 100644 --- a/crates/accelerate/src/two_qubit_decompose.rs +++ b/crates/accelerate/src/two_qubit_decompose.rs @@ -82,7 +82,7 @@ const B_NON_NORMALIZED: [[Complex64; 4]; 4] = [ ], [ Complex64::new(1., 0.), - Complex64::new(-0., -1.), + Complex64::new(0., -1.), Complex64::new(0., 0.), Complex64::new(0., 0.), ], @@ -99,7 +99,7 @@ const B_NON_NORMALIZED_DAGGER: [[Complex64; 4]; 4] = [ Complex64::new(0., -0.5), Complex64::new(0., 0.), Complex64::new(0., 0.), - Complex64::new(-0., 0.5), + Complex64::new(0., 0.5), ], [ Complex64::new(0., 0.), @@ -110,7 +110,7 @@ const B_NON_NORMALIZED_DAGGER: [[Complex64; 4]; 4] = [ [ Complex64::new(0., 0.), Complex64::new(0.5, 0.), - Complex64::new(-0.5, -0.), + Complex64::new(-0.5, 0.), Complex64::new(0., 0.), ], ]; @@ -539,18 +539,12 @@ impl TwoQubitWeylDecomposition { let mut u = unitary_matrix.as_array().to_owned(); let unitary_matrix = unitary_matrix.as_array().to_owned(); - // Faer sometimes returns negative 0s which will throw off the signs - // after the powf we do below, normalize to 0. instead by adding a - // zero complex. let det_u = - u.view().into_faer_complex().determinant().to_num_complex() + Complex64::new(0., 0.); + u.view().into_faer_complex().determinant().to_num_complex(); let det_pow = det_u.powf(-0.25); u.mapv_inplace(|x| x * det_pow); let mut global_phase = det_u.arg() / 4.; let u_p = 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"; From 5fac9acd72b937b19e5f08fce5305a83682a6890 Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Tue, 12 Mar 2024 07:16:35 +0900 Subject: [PATCH 32/46] Handle enum to int conversion as method --- crates/accelerate/src/two_qubit_decompose.rs | 63 +++++++++++--------- 1 file changed, 36 insertions(+), 27 deletions(-) diff --git a/crates/accelerate/src/two_qubit_decompose.rs b/crates/accelerate/src/two_qubit_decompose.rs index 439da22d888b..81080da75841 100644 --- a/crates/accelerate/src/two_qubit_decompose.rs +++ b/crates/accelerate/src/two_qubit_decompose.rs @@ -347,6 +347,39 @@ enum Specialization { SimabmbEquiv, } +impl Specialization { + fn to_u8(self) -> u8 { + match self { + Specialization::General => 0, + Specialization::IdEquiv => 1, + Specialization::SWAPEquiv => 2, + Specialization::PartialSWAPEquiv => 3, + Specialization::PartialSWAPFlipEquiv => 4, + Specialization::ControlledEquiv => 5, + Specialization::MirrorControlledEquiv => 6, + Specialization::SimaabEquiv => 7, + Specialization::SimabbEquiv => 8, + Specialization::SimabmbEquiv => 9, + } + } + + fn from_u8(val: u8) -> Self { + match val { + 0 => Specialization::General, + 1 => Specialization::IdEquiv, + 2 => Specialization::SWAPEquiv, + 3 => Specialization::PartialSWAPEquiv, + 4 => Specialization::PartialSWAPFlipEquiv, + 5 => Specialization::ControlledEquiv, + 6 => Specialization::MirrorControlledEquiv, + 7 => Specialization::SimaabEquiv, + 8 => Specialization::SimabbEquiv, + 9 => Specialization::SimabmbEquiv, + _ => unreachable!("Invalid specialization value"), + } + } +} + #[derive(Clone, Debug)] #[allow(non_snake_case)] #[pyclass(module = "qiskit._accelerate.two_qubit_decompose", subclass)] @@ -426,18 +459,7 @@ const IPX: [[Complex64; 2]; 2] = [ #[pymethods] impl TwoQubitWeylDecomposition { fn __getstate__(&self, py: Python) -> ([f64; 5], [PyObject; 5], u8, String, Option) { - let specialization = match self.specialization { - Specialization::General => 0, - Specialization::IdEquiv => 1, - Specialization::SWAPEquiv => 2, - Specialization::PartialSWAPEquiv => 3, - Specialization::PartialSWAPFlipEquiv => 4, - Specialization::ControlledEquiv => 5, - Specialization::MirrorControlledEquiv => 6, - Specialization::SimaabEquiv => 7, - Specialization::SimabbEquiv => 8, - Specialization::SimabmbEquiv => 9, - }; + let specialization = self.specialization.to_u8(); ( [ self.a, @@ -481,19 +503,7 @@ impl TwoQubitWeylDecomposition { self.unitary_matrix = state.1[4].as_array().to_owned(); self.default_euler_basis = state.3; self.requested_fidelity = state.4; - self.specialization = match state.2 { - 0 => Specialization::General, - 1 => Specialization::IdEquiv, - 2 => Specialization::SWAPEquiv, - 3 => Specialization::PartialSWAPEquiv, - 4 => Specialization::PartialSWAPFlipEquiv, - 5 => Specialization::ControlledEquiv, - 6 => Specialization::MirrorControlledEquiv, - 7 => Specialization::SimaabEquiv, - 8 => Specialization::SimabbEquiv, - 9 => Specialization::SimabmbEquiv, - _ => unreachable!("Invalid specialization value"), - }; + self.specialization = Specialization::from_u8(state.2); } fn __getnewargs__(&self, py: Python) -> (PyObject, Option, Option, bool) { @@ -539,8 +549,7 @@ impl TwoQubitWeylDecomposition { let mut u = unitary_matrix.as_array().to_owned(); let unitary_matrix = unitary_matrix.as_array().to_owned(); - let det_u = - u.view().into_faer_complex().determinant().to_num_complex(); + let det_u = u.view().into_faer_complex().determinant().to_num_complex(); let det_pow = det_u.powf(-0.25); u.mapv_inplace(|x| x * det_pow); let mut global_phase = det_u.arg() / 4.; From cbbdb53f6478aad603c47808177fd1bfbd6a2dcd Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Tue, 12 Mar 2024 07:21:49 +0900 Subject: [PATCH 33/46] Cleanup decompose_two_qubit_product_gate() --- crates/accelerate/src/two_qubit_decompose.rs | 36 +++++++++++--------- 1 file changed, 20 insertions(+), 16 deletions(-) diff --git a/crates/accelerate/src/two_qubit_decompose.rs b/crates/accelerate/src/two_qubit_decompose.rs index 81080da75841..51d6f561e989 100644 --- a/crates/accelerate/src/two_qubit_decompose.rs +++ b/crates/accelerate/src/two_qubit_decompose.rs @@ -165,32 +165,36 @@ 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(); + special_unitary: ArrayView2, +) -> PyResult<(Array2, Array2, f64)> { + let mut r: Array2 = special_unitary.slice(s![..2, ..2]).to_owned(); let mut det_r = det_one_qubit(r.view()); if det_r.abs() < 0.1 { - r = unitary.slice(s![2.., ..2]).to_owned(); + r = special_unitary.slice(s![2.., ..2]).to_owned(); det_r = det_one_qubit(r.view()); } - assert!( - det_r.abs() >= 0.1, - "decompose_two_qubit_product_gate: unable to decompose: detR < 0.1" - ); + if det_r.abs() < 0.1 { + import_exception!(qiskit, QiskitError); + return Err(QiskitError::new_err( + "decompose_two_qubit_product_gate: unable to decompose: detR < 0.1", + )); + } r.mapv_inplace(|x| x / det_r.sqrt()); let r_t_conj: Array2 = r.t().mapv(|x| x.conj()); let eye = aview2(&ONE_QUBIT_IDENTITY); let mut temp = kron(&eye, &r_t_conj); - temp = unitary.dot(&temp); + temp = special_unitary.dot(&temp); let mut l = temp.slice(s![..;2, ..;2]).to_owned(); let det_l = det_one_qubit(l.view()); - assert!( - det_l.abs() >= 0.9, - "decompose_two_qubit_product_gate: unable to decompose: detL < 0.9" - ); + if det_l.abs() < 0.9 { + import_exception!(qiskit, QiskitError); + return Err(QiskitError::new_err( + "decompose_two_qubit_product_gate: unable to decompose: detL < 0.9", + )); + } l.mapv_inplace(|x| x / det_l.sqrt()); let phase = det_l.arg() / 2.; - (l, r, phase) + Ok((l, r, phase)) } fn __weyl_coordinates(unitary: MatRef) -> [f64; 3] { @@ -650,9 +654,9 @@ impl TwoQubitWeylDecomposition { 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()); + 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()); + let (K2l, mut K2r, phase_r) = decompose_two_qubit_product_gate(k2.view())?; global_phase += phase_l + phase_r; // Flip into Weyl chamber From 9a914c2a08fba7597490222aea06698be0b8b299 Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Tue, 12 Mar 2024 07:28:18 +0900 Subject: [PATCH 34/46] Use a trait for trace_to_fid() --- crates/accelerate/src/two_qubit_decompose.rs | 39 +++++++++++--------- 1 file changed, 21 insertions(+), 18 deletions(-) diff --git a/crates/accelerate/src/two_qubit_decompose.rs b/crates/accelerate/src/two_qubit_decompose.rs index 51d6f561e989..11df3b5e6b6c 100644 --- a/crates/accelerate/src/two_qubit_decompose.rs +++ b/crates/accelerate/src/two_qubit_decompose.rs @@ -164,6 +164,24 @@ impl Arg for c64 { } } +pub trait TraceToFidelity { + /// Average gate fidelity is :math:`Fbar = (d + |Tr (Utarget \\cdot U^dag)|^2) / d(d+1)` + /// M. Horodecki, P. Horodecki and R. Horodecki, PRA 60, 1888 (1999) + fn trace_to_fid(self) -> f64; +} + +impl TraceToFidelity for Complex64 { + fn trace_to_fid(self) -> f64 { + (4.0 + self.abs().powi(2)) / 20.0 + } +} + +impl TraceToFidelity for c64 { + fn trace_to_fid(self) -> f64 { + (4.0 + self.faer_abs2()) / 20.0 + } +} + fn decompose_two_qubit_product_gate( special_unitary: ArrayView2, ) -> PyResult<(Array2, Array2, f64)> { @@ -279,27 +297,12 @@ fn __num_basis_gates(basis_b: f64, basis_fidelity: f64, unitary: MatRef) -> traces .into_iter() .enumerate() - .map(|(idx, trace)| { - ( - idx, - trace_to_fid_c64(&trace) * basis_fidelity.powi(idx as i32), - ) - }) + .map(|(idx, trace)| (idx, trace.trace_to_fid() * basis_fidelity.powi(idx as i32))) .min_by(|(_idx1, fid1), (_idx2, fid2)| fid2.partial_cmp(fid1).unwrap()) .unwrap() .0 } -/// Average gate fidelity is :math:`Fbar = (d + |Tr (Utarget \\cdot U^dag)|^2) / d(d+1)` -/// M. Horodecki, P. Horodecki and R. Horodecki, PRA 60, 1888 (1999) -fn trace_to_fid_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 { @@ -720,7 +723,7 @@ impl TwoQubitWeylDecomposition { da.sin() * db.sin() * dc.sin(), ); match fidelity { - Some(fid) => trace_to_fid(tr) >= fid, + Some(fid) => tr.trace_to_fid() >= fid, // Set to false here to default to general specialization in the absence of a // fidelity and provided specialization. None => false, @@ -954,7 +957,7 @@ impl TwoQubitWeylDecomposition { da.sin() * db.sin() * dc.sin(), ) }; - specialized.calculated_fidelity = trace_to_fid(tr); + specialized.calculated_fidelity = tr.trace_to_fid(); if let Some(fid) = specialized.requested_fidelity { if specialized.calculated_fidelity + 1.0e-13 < fid { import_exception!(qiskit, QiskitError); From 33d073d2e05ecaa326a9bf70f02c1d104d29d703 Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Tue, 12 Mar 2024 08:30:10 +0900 Subject: [PATCH 35/46] Use an enum instead of string for euler basis --- .../src/euler_one_qubit_decomposer.rs | 140 ++++++++++++------ crates/accelerate/src/two_qubit_decompose.rs | 48 +++--- 2 files changed, 122 insertions(+), 66 deletions(-) diff --git a/crates/accelerate/src/euler_one_qubit_decomposer.rs b/crates/accelerate/src/euler_one_qubit_decomposer.rs index ba1dfe2a711c..6af60400a910 100644 --- a/crates/accelerate/src/euler_one_qubit_decomposer.rs +++ b/crates/accelerate/src/euler_one_qubit_decomposer.rs @@ -11,6 +11,7 @@ // that they have been altered from the originals. #![allow(clippy::too_many_arguments)] +#![allow(clippy::upper_case_acronyms)] use hashbrown::HashMap; use num_complex::{Complex64, ComplexFloat}; @@ -409,7 +410,7 @@ fn circuit_rr( #[pyfunction] pub fn generate_circuit( - target_basis: &str, + target_basis: &EulerBasis, theta: f64, phi: f64, lam: f64, @@ -418,14 +419,14 @@ pub fn generate_circuit( atol: Option, ) -> PyResult { let res = match target_basis { - "ZYZ" => circuit_kak(theta, phi, lam, phase, "rz", "ry", simplify, atol), - "ZXZ" => circuit_kak(theta, phi, lam, phase, "rz", "rx", simplify, atol), - "XZX" => circuit_kak(theta, phi, lam, phase, "rx", "rz", simplify, atol), - "XYX" => circuit_kak(theta, phi, lam, phase, "rx", "ry", simplify, atol), - "U3" => circuit_u3(theta, phi, lam, phase, simplify, atol), - "U321" => circuit_u321(theta, phi, lam, phase, simplify, atol), - "U" => circuit_u(theta, phi, lam, phase, simplify, atol), - "PSX" => { + EulerBasis::ZYZ => circuit_kak(theta, phi, lam, phase, "rz", "ry", simplify, atol), + EulerBasis::ZXZ => circuit_kak(theta, phi, lam, phase, "rz", "rx", simplify, atol), + EulerBasis::XZX => circuit_kak(theta, phi, lam, phase, "rx", "rz", simplify, atol), + EulerBasis::XYX => circuit_kak(theta, phi, lam, phase, "rx", "ry", simplify, atol), + EulerBasis::U3 => circuit_u3(theta, phi, lam, phase, simplify, atol), + EulerBasis::U321 => circuit_u321(theta, phi, lam, phase, simplify, atol), + EulerBasis::U => circuit_u(theta, phi, lam, phase, simplify, atol), + EulerBasis::PSX => { let mut inner_atol = match atol { Some(atol) => atol, None => ANGLE_ZERO_EPSILON, @@ -455,7 +456,7 @@ pub fn generate_circuit( None::>, ) } - "ZSX" => { + EulerBasis::ZSX => { let mut inner_atol = match atol { Some(atol) => atol, None => ANGLE_ZERO_EPSILON, @@ -485,7 +486,7 @@ pub fn generate_circuit( None::>, ) } - "U1X" => { + EulerBasis::U1X => { let mut inner_atol = match atol { Some(atol) => atol, None => ANGLE_ZERO_EPSILON, @@ -515,7 +516,7 @@ pub fn generate_circuit( None::>, ) } - "ZSXX" => { + EulerBasis::ZSXX => { let mut inner_atol = match atol { Some(atol) => atol, None => ANGLE_ZERO_EPSILON, @@ -548,32 +549,88 @@ pub fn generate_circuit( Some(fnxpi), ) } - "RR" => circuit_rr(theta, phi, lam, phase, simplify, atol), - other => { - return Err(PyTypeError::new_err(format!( - "Invalid target basis: {other}" - ))) - } + EulerBasis::RR => circuit_rr(theta, phi, lam, phase, simplify, atol), }; Ok(res) } +#[derive(Clone, Debug, Copy)] +#[pyclass(module = "qiskit._accelerate.euler_one_qubit_decomposer")] +pub enum EulerBasis { + U321, + U3, + U, + PSX, + ZSX, + ZSXX, + U1X, + RR, + ZYZ, + ZXZ, + XYX, + XZX, +} + +#[pymethods] +impl EulerBasis { + #![allow(clippy::wrong_self_convention)] + pub fn to_str(&self) -> String { + match self { + Self::U321 => "U321".to_string(), + Self::U3 => "U3".to_string(), + Self::U => "U".to_string(), + Self::PSX => "PSX".to_string(), + Self::ZSX => "ZSX".to_string(), + Self::ZSXX => "ZSXX".to_string(), + Self::U1X => "U1X".to_string(), + Self::RR => "RR".to_string(), + Self::ZYZ => "ZYZ".to_string(), + Self::ZXZ => "ZXZ".to_string(), + Self::XYX => "XYX".to_string(), + Self::XZX => "XZX".to_string(), + } + } + + #[staticmethod] + pub fn from_string(input: &str) -> PyResult { + let res = match input { + "U321" => EulerBasis::U321, + "U3" => EulerBasis::U3, + "U" => EulerBasis::U, + "PSX" => EulerBasis::PSX, + "ZSX" => EulerBasis::ZSX, + "ZSXX" => EulerBasis::ZSXX, + "U1X" => EulerBasis::U1X, + "RR" => EulerBasis::RR, + "ZYZ" => EulerBasis::ZYZ, + "ZXZ" => EulerBasis::ZXZ, + "XYX" => EulerBasis::XYX, + "XZX" => EulerBasis::XZX, + basis => { + return Err(PyTypeError::new_err(format!( + "Invalid target basis {basis}" + ))); + } + }; + Ok(res) + } +} + #[inline] -pub fn angles_from_unitary(unitary: ArrayView2, target_basis: &str) -> [f64; 4] { +pub fn angles_from_unitary(unitary: ArrayView2, target_basis: EulerBasis) -> [f64; 4] { match target_basis { - "U321" => params_u3_inner(unitary), - "U3" => params_u3_inner(unitary), - "U" => params_u3_inner(unitary), - "PSX" => params_u1x_inner(unitary), - "ZSX" => params_u1x_inner(unitary), - "ZSXX" => params_u1x_inner(unitary), - "U1X" => params_u1x_inner(unitary), - "RR" => params_zyz_inner(unitary), - "ZYZ" => params_zyz_inner(unitary), - "ZXZ" => params_zxz_inner(unitary), - "XYX" => params_xyx_inner(unitary), - "XZX" => params_xzx_inner(unitary), - &_ => unreachable!(), + EulerBasis::U321 => params_u3_inner(unitary), + EulerBasis::U3 => params_u3_inner(unitary), + EulerBasis::U => params_u3_inner(unitary), + EulerBasis::PSX => params_u1x_inner(unitary), + EulerBasis::ZSX => params_u1x_inner(unitary), + EulerBasis::ZSXX => params_u1x_inner(unitary), + EulerBasis::U1X => params_u1x_inner(unitary), + EulerBasis::RR => params_zyz_inner(unitary), + EulerBasis::ZYZ => params_zyz_inner(unitary), + EulerBasis::ZXZ => params_zxz_inner(unitary), + EulerBasis::XYX => params_xyx_inner(unitary), + EulerBasis::XZX => params_xzx_inner(unitary), } } @@ -643,20 +700,15 @@ pub fn unitary_to_gate_sequence( simplify: bool, atol: Option, ) -> PyResult> { - const VALID_BASES: [&str; 12] = [ - "U321", "U3", "U", "PSX", "ZSX", "ZSXX", "U1X", "RR", "ZYZ", "ZXZ", "XYX", "XZX", - ]; - for basis in &target_basis_list { - if !VALID_BASES.contains(basis) { - return Err(PyTypeError::new_err(format!( - "Invalid target basis {basis}" - ))); - } + let mut target_basis_vec: Vec = Vec::with_capacity(target_basis_list.len()); + for basis in target_basis_list { + let basis_enum = EulerBasis::from_string(basis)?; + target_basis_vec.push(basis_enum) } let unitary_mat = unitary.as_array(); Ok(unitary_to_gate_sequence_inner( unitary_mat, - &target_basis_list, + &target_basis_vec, qubit, error_map, simplify, @@ -667,7 +719,7 @@ pub fn unitary_to_gate_sequence( #[inline] pub fn unitary_to_gate_sequence_inner( unitary_mat: ArrayView2, - target_basis_list: &[&str], + target_basis_list: &[EulerBasis], qubit: usize, error_map: Option<&OneQubitGateErrorMap>, simplify: bool, @@ -676,7 +728,7 @@ pub fn unitary_to_gate_sequence_inner( target_basis_list .iter() .map(|target_basis| { - let [theta, phi, lam, phase] = angles_from_unitary(unitary_mat, target_basis); + let [theta, phi, lam, phase] = angles_from_unitary(unitary_mat, *target_basis); generate_circuit(target_basis, theta, phi, lam, phase, simplify, atol).unwrap() }) .min_by(|a, b| { diff --git a/crates/accelerate/src/two_qubit_decompose.rs b/crates/accelerate/src/two_qubit_decompose.rs index 11df3b5e6b6c..8c00dfc73354 100644 --- a/crates/accelerate/src/two_qubit_decompose.rs +++ b/crates/accelerate/src/two_qubit_decompose.rs @@ -41,7 +41,8 @@ use numpy::PyReadonlyArray2; use numpy::ToPyArray; use crate::euler_one_qubit_decomposer::{ - angles_from_unitary, det_one_qubit, unitary_to_gate_sequence_inner, ANGLE_ZERO_EPSILON, + angles_from_unitary, det_one_qubit, unitary_to_gate_sequence_inner, EulerBasis, + ANGLE_ZERO_EPSILON, }; use crate::utils; @@ -405,7 +406,7 @@ pub struct TwoQubitWeylDecomposition { K2r: Array2, #[pyo3(get)] specialization: Specialization, - default_euler_basis: String, + default_euler_basis: EulerBasis, #[pyo3(get)] requested_fidelity: Option, #[pyo3(get)] @@ -483,7 +484,7 @@ impl TwoQubitWeylDecomposition { self.unitary_matrix.to_pyarray(py).into(), ], specialization, - self.default_euler_basis.clone(), + self.default_euler_basis.to_str(), self.requested_fidelity, ) } @@ -508,7 +509,7 @@ impl TwoQubitWeylDecomposition { self.K2l = state.1[2].as_array().to_owned(); self.K2r = state.1[3].as_array().to_owned(); self.unitary_matrix = state.1[4].as_array().to_owned(); - self.default_euler_basis = state.3; + self.default_euler_basis = EulerBasis::from_string(&state.3).unwrap(); self.requested_fidelity = state.4; self.specialization = Specialization::from_u8(state.2); } @@ -544,7 +545,7 @@ impl TwoQubitWeylDecomposition { K1r: Array2::zeros((0, 0)), K2r: Array2::zeros((0, 0)), specialization: Specialization::General, - default_euler_basis: "ZYZ".to_string(), + default_euler_basis: EulerBasis::ZYZ, requested_fidelity: fidelity, calculated_fidelity: 0., unitary_matrix: Array2::zeros((0, 0)), @@ -562,7 +563,7 @@ impl TwoQubitWeylDecomposition { let mut global_phase = det_u.arg() / 4.; let u_p = transform_from_magic_basis(u.view(), true); let m2 = u_p.t().dot(&u_p); - let default_euler_basis = "ZYZ"; + let default_euler_basis = EulerBasis::ZYZ; // M2 is a symmetric complex matrix. We need to decompose it as M2 = P D P^T where // P ∈ SO(4), D is diagonal with unit-magnitude elements. @@ -769,7 +770,7 @@ impl TwoQubitWeylDecomposition { K2l, K2r, specialization: Specialization::General, - default_euler_basis: default_euler_basis.to_owned(), + default_euler_basis, requested_fidelity: fidelity, calculated_fidelity: -1.0, unitary_matrix, @@ -848,11 +849,11 @@ impl TwoQubitWeylDecomposition { } } Specialization::ControlledEquiv => { - let euler_basis = "XYX".to_owned(); + let euler_basis = EulerBasis::XYX; let [k2ltheta, k2lphi, k2llambda, k2lphase] = - angles_from_unitary(general.K2l.view(), &euler_basis); + angles_from_unitary(general.K2l.view(), euler_basis); let [k2rtheta, k2rphi, k2rlambda, k2rphase] = - angles_from_unitary(general.K2r.view(), &euler_basis); + angles_from_unitary(general.K2r.view(), euler_basis); TwoQubitWeylDecomposition { specialization, a, @@ -869,9 +870,9 @@ impl TwoQubitWeylDecomposition { } Specialization::MirrorControlledEquiv => { let [k2ltheta, k2lphi, k2llambda, k2lphase] = - angles_from_unitary(general.K2l.view(), "ZYZ"); + angles_from_unitary(general.K2l.view(), EulerBasis::ZYZ); let [k2rtheta, k2rphi, k2rlambda, k2rphase] = - angles_from_unitary(general.K2r.view(), "ZYZ"); + angles_from_unitary(general.K2r.view(), EulerBasis::ZYZ); TwoQubitWeylDecomposition { specialization, a: PI4, @@ -887,7 +888,7 @@ impl TwoQubitWeylDecomposition { } Specialization::SimaabEquiv => { let [k2ltheta, k2lphi, k2llambda, k2lphase] = - angles_from_unitary(general.K2l.view(), "ZYZ"); + angles_from_unitary(general.K2l.view(), EulerBasis::ZYZ); TwoQubitWeylDecomposition { specialization, a: (a + b) / 2., @@ -902,9 +903,9 @@ impl TwoQubitWeylDecomposition { } } Specialization::SimabbEquiv => { - let euler_basis = "XYX".to_owned(); + let euler_basis = EulerBasis::XYX; let [k2ltheta, k2lphi, k2llambda, k2lphase] = - angles_from_unitary(general.K2l.view(), &euler_basis); + angles_from_unitary(general.K2l.view(), euler_basis); TwoQubitWeylDecomposition { specialization, a, @@ -920,9 +921,9 @@ impl TwoQubitWeylDecomposition { } } Specialization::SimabmbEquiv => { - let euler_basis = "XYX".to_owned(); + let euler_basis = EulerBasis::XYX; let [k2ltheta, k2lphi, k2llambda, k2lphase] = - angles_from_unitary(general.K2l.view(), &euler_basis); + angles_from_unitary(general.K2l.view(), euler_basis); TwoQubitWeylDecomposition { specialization, a, @@ -1008,9 +1009,12 @@ impl TwoQubitWeylDecomposition { euler_basis: Option<&str>, simplify: bool, atol: Option, - ) -> TwoQubitGateSequence { - let euler_basis: &str = euler_basis.unwrap_or(&self.default_euler_basis); - let target_1q_basis_list: Vec<&str> = vec![euler_basis]; + ) -> PyResult { + let euler_basis: EulerBasis = match euler_basis { + Some(basis) => EulerBasis::from_string(basis)?, + None => self.default_euler_basis, + }; + let target_1q_basis_list: Vec = vec![euler_basis]; let mut gate_sequence = Vec::new(); let mut global_phase: f64 = self.global_phase; @@ -1072,10 +1076,10 @@ impl TwoQubitWeylDecomposition { for gate in c1l.gates { gate_sequence.push((gate.0, gate.1, smallvec![1])) } - TwoQubitGateSequence { + Ok(TwoQubitGateSequence { gates: gate_sequence, global_phase, - } + }) } } From 160d70552fcc4a3a29f5d3801849b82b6bebdb20 Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Tue, 12 Mar 2024 08:32:04 +0900 Subject: [PATCH 36/46] Fix release note typo --- releasenotes/notes/rust-two-qubit-weyl-ec551f3f9c812124.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/releasenotes/notes/rust-two-qubit-weyl-ec551f3f9c812124.yaml b/releasenotes/notes/rust-two-qubit-weyl-ec551f3f9c812124.yaml index 57379286f82f..12d5c27d9af8 100644 --- a/releasenotes/notes/rust-two-qubit-weyl-ec551f3f9c812124.yaml +++ b/releasenotes/notes/rust-two-qubit-weyl-ec551f3f9c812124.yaml @@ -23,5 +23,5 @@ deprecations_synthesis: :class:`.TwoQubitWeylDecomposition` object. Despite this it was still a documented part of the public API for the class and is now being deprecated without any potential replacement. This release it always will - raise a ``NotImplementedError` when called because the specialization + raise a ``NotImplementedError`` when called because the specialization subclassing has been removed as part of the Rust rewrite of the class. From 9bf2f155888399d1686ef0b2cf94c9e229f01720 Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Tue, 12 Mar 2024 21:57:59 +0900 Subject: [PATCH 37/46] Improve magic basis transformation functions Co-authored-by: Jake Lishman --- crates/accelerate/src/two_qubit_decompose.rs | 41 ++++++++++---------- 1 file changed, 21 insertions(+), 20 deletions(-) diff --git a/crates/accelerate/src/two_qubit_decompose.rs b/crates/accelerate/src/two_qubit_decompose.rs index 8c00dfc73354..16ced4bf8a10 100644 --- a/crates/accelerate/src/two_qubit_decompose.rs +++ b/crates/accelerate/src/two_qubit_decompose.rs @@ -116,28 +116,29 @@ const B_NON_NORMALIZED_DAGGER: [[Complex64; 4]; 4] = [ ], ]; -fn transform_from_magic_basis(unitary: ArrayView2, reverse: bool) -> Array2 { +enum MagicBasisTransform { + Into, + OutOf, +} + +fn magic_basis_transform( + unitary: ArrayView2, + direction: MagicBasisTransform, +) -> Array2 { let _b_nonnormalized = aview2(&B_NON_NORMALIZED); let _b_nonnormalized_dagger = aview2(&B_NON_NORMALIZED_DAGGER); - if reverse { - _b_nonnormalized_dagger.dot(&unitary).dot(&_b_nonnormalized) - } else { - _b_nonnormalized.dot(&unitary).dot(&_b_nonnormalized_dagger) + match direction { + MagicBasisTransform::OutOf => _b_nonnormalized_dagger.dot(&unitary).dot(&_b_nonnormalized), + MagicBasisTransform::Into => _b_nonnormalized.dot(&unitary).dot(&_b_nonnormalized_dagger), } } -fn transform_from_magic_basis_faer(u: Mat, reverse: bool) -> Mat { +fn transform_from_magic_basis(u: Mat) -> 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() + magic_basis_transform(unitary, MagicBasisTransform::OutOf) + .view() + .into_faer_complex() + .to_owned() } // faer::c64 and num_complex::Complex are both structs @@ -218,7 +219,7 @@ fn decompose_two_qubit_product_gate( fn __weyl_coordinates(unitary: MatRef) -> [f64; 3] { let uscaled = scale(C1 / unitary.determinant().powf(0.25)) * unitary; - let uup = transform_from_magic_basis_faer(uscaled, true); + let uup = transform_from_magic_basis(uscaled); let mut darg: Vec<_> = (uup.transpose() * &uup) .complex_eigenvalues() .into_iter() @@ -561,7 +562,7 @@ impl TwoQubitWeylDecomposition { let det_pow = det_u.powf(-0.25); u.mapv_inplace(|x| x * det_pow); let mut global_phase = det_u.arg() / 4.; - let u_p = transform_from_magic_basis(u.view(), true); + let u_p = magic_basis_transform(u.view(), MagicBasisTransform::OutOf); let m2 = u_p.t().dot(&u_p); let default_euler_basis = EulerBasis::ZYZ; @@ -654,8 +655,8 @@ impl TwoQubitWeylDecomposition { .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); + let k1 = magic_basis_transform(u_p.dot(&p).dot(&temp).view(), MagicBasisTransform::Into); + let k2 = magic_basis_transform(p.t(), MagicBasisTransform::Into); #[allow(non_snake_case)] let (mut K1l, mut K1r, phase_l) = decompose_two_qubit_product_gate(k1.view())?; From 6712ade9d4b57a03efea7cdf00de48693a8c3052 Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Wed, 13 Mar 2024 06:13:34 +0900 Subject: [PATCH 38/46] Remove eigh() util function --- crates/accelerate/src/utils.rs | 21 +-------------------- 1 file changed, 1 insertion(+), 20 deletions(-) diff --git a/crates/accelerate/src/utils.rs b/crates/accelerate/src/utils.rs index 90d74a9bd516..bc8a66568f9b 100644 --- a/crates/accelerate/src/utils.rs +++ b/crates/accelerate/src/utils.rs @@ -15,9 +15,7 @@ use pyo3::types::PySlice; use faer::prelude::*; use faer::IntoFaerComplex; -use faer::IntoNdarray; -use faer::Side::Lower; -use num_complex::{Complex, Complex64}; +use num_complex::Complex; use numpy::{IntoPyArray, PyReadonlyArray2}; /// A private enumeration type used to extract arguments to pymethod @@ -55,25 +53,8 @@ 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(()) } From 27a7587a0c79b6038b07124dca14b7419a43da57 Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Wed, 13 Mar 2024 06:17:35 +0900 Subject: [PATCH 39/46] Revert unecessary changes to callers of TwoQubitWeylDecomposition --- qiskit/synthesis/two_qubit/two_qubit_decompose.py | 4 ++-- test/python/synthesis/test_synthesis.py | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/qiskit/synthesis/two_qubit/two_qubit_decompose.py b/qiskit/synthesis/two_qubit/two_qubit_decompose.py index 5a1ca636e9bc..41168547f5b4 100644 --- a/qiskit/synthesis/two_qubit/two_qubit_decompose.py +++ b/qiskit/synthesis/two_qubit/two_qubit_decompose.py @@ -272,7 +272,7 @@ 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).to_matrix()) + decomp = TwoQubitWeylDecomposition(rxx_equivalent_gate(test_angle)) circ = QuantumCircuit(2) circ.rxx(test_angle, 0, 1) @@ -317,7 +317,7 @@ def __call__(self, unitary, *, atol=DEFAULT_ATOL) -> QuantumCircuit: """ # pylint: disable=attribute-defined-outside-init - self.decomposer = TwoQubitWeylDecomposition(np.asarray(unitary, dtype=complex)) + self.decomposer = TwoQubitWeylDecomposition(unitary) oneq_decompose = OneQubitEulerDecomposer("ZYZ") c1l, c1r, c2l, c2r = ( diff --git a/test/python/synthesis/test_synthesis.py b/test/python/synthesis/test_synthesis.py index 0e891e4b0392..bba28f280600 100644 --- a/test/python/synthesis/test_synthesis.py +++ b/test/python/synthesis/test_synthesis.py @@ -620,13 +620,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.data, fidelity=0.99) + weyl1 = TwoQubitWeylDecomposition(target, 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.data, fidelity=0.99) + weyl1 = TwoQubitWeylDecomposition(target, fidelity=0.99) self.assertRoundTripPickle(weyl1) def test_two_qubit_weyl_decomposition_cnot(self): From 03c4af53523f4fabbb72d82fe4a0bea61d09ecc4 Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Wed, 13 Mar 2024 06:28:52 +0900 Subject: [PATCH 40/46] Restore debug logging and add test assertions for it --- .../two_qubit/two_qubit_decompose.py | 13 ++++++ test/python/synthesis/test_synthesis.py | 41 ++++++++++++++----- 2 files changed, 44 insertions(+), 10 deletions(-) diff --git a/qiskit/synthesis/two_qubit/two_qubit_decompose.py b/qiskit/synthesis/two_qubit/two_qubit_decompose.py index 41168547f5b4..41e2d18bb4ed 100644 --- a/qiskit/synthesis/two_qubit/two_qubit_decompose.py +++ b/qiskit/synthesis/two_qubit/two_qubit_decompose.py @@ -173,6 +173,19 @@ def __init__( self.unitary_matrix = unitary_matrix self.requested_fidelity = fidelity self.calculated_fidelity = self._inner_decomposition.calculated_fidelity + 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, + ) @deprecate_func(since="1.1.0", removal_timeline="in the 2.0.0 release") def specialize(self): diff --git a/test/python/synthesis/test_synthesis.py b/test/python/synthesis/test_synthesis.py index bba28f280600..864bbd334239 100644 --- a/test/python/synthesis/test_synthesis.py +++ b/test/python/synthesis/test_synthesis.py @@ -14,6 +14,8 @@ import pickle import unittest +import contextlib +import logging import math import numpy as np import scipy @@ -130,10 +132,24 @@ 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) - weyl2: TwoQubitWeylDecomposition = eval(repr1) # pylint: disable=eval-used + with self.assertDebugOnly(): + 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)) @@ -176,7 +192,8 @@ 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 - decomp = TwoQubitWeylDecomposition(target_unitary, fidelity=None) + with self.assertDebugOnly(): + 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 ( @@ -204,12 +221,14 @@ def check_two_qubit_weyl_specialization( # Loop to check both for implicit and explicity specialization for decomposer in (TwoQubitWeylDecomposition, expected_specialization): if isinstance(decomposer, TwoQubitWeylDecomposition): - decomp = decomposer(target_unitary, fidelity=fidelity) + with self.assertDebugOnly(): + decomp = decomposer(target_unitary, fidelity=fidelity) decomp_name = decomp.specialization else: - decomp = TwoQubitWeylDecomposition( - target_unitary, fidelity=None, _specialization=expected_specialization - ) + with self.assertDebugOnly(): + decomp = TwoQubitWeylDecomposition( + target_unitary, fidelity=None, _specialization=expected_specialization + ) decomp_name = expected_specialization self.assertRoundTrip(decomp) self.assertRoundTripPickle(decomp) @@ -233,9 +252,10 @@ def check_two_qubit_weyl_specialization( actual_unitary = Operator(circ).data trace = np.trace(actual_unitary.T.conj() @ target_unitary) self.assertAlmostEqual(trace.imag, 0, places=13, msg=f"Real trace for {decomp_name}") - decomp2 = TwoQubitWeylDecomposition( - target_unitary, fidelity=None, _specialization=expected_specialization - ) # Shouldn't raise + with self.assertDebugOnly(): + decomp2 = TwoQubitWeylDecomposition( + target_unitary, fidelity=None, _specialization=expected_specialization + ) # Shouldn't raise self.assertRoundTrip(decomp2) self.assertRoundTripPickle(decomp2) if expected_specialization != Specialization.General: @@ -1064,7 +1084,8 @@ def test_approx_supercontrolled_decompose_random(self, seed): tgt = random_unitary(4, seed=state).data tgt *= np.exp(1j * tgt_phase) - traces_pred = decomposer.traces(TwoQubitWeylDecomposition(tgt)) + with self.assertDebugOnly(): + traces_pred = decomposer.traces(TwoQubitWeylDecomposition(tgt)) for i in range(4): with self.subTest(i=i): From 84e9c6752cf1693d56e54e05ed03ceb01ead3fe8 Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Thu, 14 Mar 2024 11:01:23 +0900 Subject: [PATCH 41/46] Update qiskit/synthesis/two_qubit/two_qubit_decompose.py Co-authored-by: Lev Bishop <18673315+levbishop@users.noreply.github.com> --- qiskit/synthesis/two_qubit/two_qubit_decompose.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/qiskit/synthesis/two_qubit/two_qubit_decompose.py b/qiskit/synthesis/two_qubit/two_qubit_decompose.py index 41e2d18bb4ed..5067b27f9258 100644 --- a/qiskit/synthesis/two_qubit/two_qubit_decompose.py +++ b/qiskit/synthesis/two_qubit/two_qubit_decompose.py @@ -230,7 +230,7 @@ def __repr__(self): + b64ascii + [ f"requested_fidelity={self.requested_fidelity},", - f"_specialization={self._inner_decomposition.specialization}," + f"_specialization={self._inner_decomposition.specialization},", f"calculated_fidelity={self.calculated_fidelity},", f"actual_fidelity={self.actual_fidelity()},", f"abc={(self.a, self.b, self.c)})", From c52bffa7c071da01d1a9a21e0f15c73748617e84 Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Thu, 14 Mar 2024 10:59:32 +0900 Subject: [PATCH 42/46] Add specialization to __str__ --- qiskit/synthesis/two_qubit/two_qubit_decompose.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/qiskit/synthesis/two_qubit/two_qubit_decompose.py b/qiskit/synthesis/two_qubit/two_qubit_decompose.py index 5067b27f9258..ac7e37fe5406 100644 --- a/qiskit/synthesis/two_qubit/two_qubit_decompose.py +++ b/qiskit/synthesis/two_qubit/two_qubit_decompose.py @@ -256,7 +256,8 @@ def from_bytes( return cls(arr, fidelity=requested_fidelity, _specialization=_specialization) def __str__(self): - pre = f"{self.__class__.__name__}(\n\t" + specialization = str(self._inner_decomposition.specialization).split(".")[1] + pre = f"{self.__class__.__name__} [specialization={specialization}] (\n\t" circ_indent = "\n\t".join(self.circuit(simplify=True).draw("text").lines(-1)) return f"{pre}{circ_indent}\n)" From 76d2ccf42849e29d96d46e15dd1eb58c1802a7bf Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Thu, 14 Mar 2024 11:12:57 +0900 Subject: [PATCH 43/46] Add previous specialized class docstrings as inline rust code comments --- crates/accelerate/src/two_qubit_decompose.rs | 51 ++++++++++++++++++++ 1 file changed, 51 insertions(+) diff --git a/crates/accelerate/src/two_qubit_decompose.rs b/crates/accelerate/src/two_qubit_decompose.rs index 16ced4bf8a10..be3a32255293 100644 --- a/crates/accelerate/src/two_qubit_decompose.rs +++ b/crates/accelerate/src/two_qubit_decompose.rs @@ -777,6 +777,10 @@ impl TwoQubitWeylDecomposition { unitary_matrix, }; let mut specialized: TwoQubitWeylDecomposition = match specialization { + // :math:`U \sim U_d(0,0,0) \sim Id` + // + // This gate binds 0 parameters, we make it canonical by setting + // :math:`K2_l = Id` , :math:`K2_r = Id`. Specialization::IdEquiv => TwoQubitWeylDecomposition { specialization, a: 0., @@ -788,6 +792,10 @@ impl TwoQubitWeylDecomposition { K2r: Array2::eye(2), ..general }, + // :math:`U \sim U_d(\pi/4, \pi/4, \pi/4) \sim U(\pi/4, \pi/4, -\pi/4) \sim \text{SWAP}` + // + // This gate binds 0 parameters, we make it canonical by setting + // :math:`K2_l = Id` , :math:`K2_r = Id`. Specialization::SWAPEquiv => { if c > 0. { TwoQubitWeylDecomposition { @@ -817,6 +825,11 @@ impl TwoQubitWeylDecomposition { } } } + // :math:`U \sim U_d(\alpha\pi/4, \alpha\pi/4, \alpha\pi/4) \sim \text{SWAP}^\alpha` + // + // This gate binds 3 parameters, we make it canonical by setting: + // + // :math:`K2_l = Id`. Specialization::PartialSWAPEquiv => { let closest = closest_partial_swap(a, b, c); let mut k2l_dag = general.K2l.t().to_owned(); @@ -833,6 +846,14 @@ impl TwoQubitWeylDecomposition { ..general } } + // :math:`U \sim U_d(\alpha\pi/4, \alpha\pi/4, -\alpha\pi/4) \sim \text{SWAP}^\alpha` + // + // (a non-equivalent root of SWAP from the TwoQubitWeylPartialSWAPEquiv + // similar to how :math:`x = (\pm \sqrt(x))^2`) + // + // This gate binds 3 parameters, we make it canonical by setting: + // + // :math:`K2_l = Id` Specialization::PartialSWAPFlipEquiv => { let closest = closest_partial_swap(a, b, -c); let mut k2l_dag = general.K2l.t().to_owned(); @@ -849,6 +870,12 @@ impl TwoQubitWeylDecomposition { ..general } } + // :math:`U \sim U_d(\alpha, 0, 0) \sim \text{Ctrl-U}` + // + // This gate binds 4 parameters, we make it canonical by setting: + // + // :math:`K2_l = Ry(\theta_l) Rx(\lambda_l)` , + // :math:`K2_r = Ry(\theta_r) Rx(\lambda_r)` . Specialization::ControlledEquiv => { let euler_basis = EulerBasis::XYX; let [k2ltheta, k2lphi, k2llambda, k2lphase] = @@ -869,6 +896,11 @@ impl TwoQubitWeylDecomposition { ..general } } + // :math:`U \sim U_d(\pi/4, \pi/4, \alpha) \sim \text{SWAP} \cdot \text{Ctrl-U}` + // + // This gate binds 4 parameters, we make it canonical by setting: + // + // :math:`K2_l = Ry(\theta_l)\cdot Rz(\lambda_l)` , :math:`K2_r = Ry(\theta_r)\cdot Rz(\lambda_r)` Specialization::MirrorControlledEquiv => { let [k2ltheta, k2lphi, k2llambda, k2lphase] = angles_from_unitary(general.K2l.view(), EulerBasis::ZYZ); @@ -887,6 +919,11 @@ impl TwoQubitWeylDecomposition { ..general } } + // :math:`U \sim U_d(\alpha, \alpha, \beta), \alpha \geq |\beta|` + // + // This gate binds 5 parameters, we make it canonical by setting: + // + // :math:`K2_l = Ry(\theta_l)\cdot Rz(\lambda_l)`. Specialization::SimaabEquiv => { let [k2ltheta, k2lphi, k2llambda, k2lphase] = angles_from_unitary(general.K2l.view(), EulerBasis::ZYZ); @@ -903,6 +940,11 @@ impl TwoQubitWeylDecomposition { ..general } } + // :math:`U \sim U_d(\alpha, \beta, -\beta), \alpha \geq \beta \geq 0` + // + // This gate binds 5 parameters, we make it canonical by setting: + // + // :math:`K2_l = Ry(\theta_l)Rx(\lambda_l)` Specialization::SimabbEquiv => { let euler_basis = EulerBasis::XYX; let [k2ltheta, k2lphi, k2llambda, k2lphase] = @@ -921,6 +963,11 @@ impl TwoQubitWeylDecomposition { ..general } } + // :math:`U \sim U_d(\alpha, \beta, -\beta), \alpha \geq \beta \geq 0` + // + // This gate binds 5 parameters, we make it canonical by setting: + // + // :math:`K2_l = Ry(\theta_l)Rx(\lambda_l)` Specialization::SimabmbEquiv => { let euler_basis = EulerBasis::XYX; let [k2ltheta, k2lphi, k2llambda, k2lphase] = @@ -939,6 +986,10 @@ impl TwoQubitWeylDecomposition { ..general } } + // U has no special symmetry. + // + // This gate binds all 6 possible parameters, so there is no need to make the single-qubit + // pre-/post-gates canonical. Specialization::General => general, }; From ade9c8010bc358f6abe2bb8ee180790a24be2dd4 Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Thu, 14 Mar 2024 11:16:30 +0900 Subject: [PATCH 44/46] Rename fSim variants and suprress rustc warning about camel case --- crates/accelerate/src/two_qubit_decompose.rs | 33 +++++++++++--------- 1 file changed, 18 insertions(+), 15 deletions(-) diff --git a/crates/accelerate/src/two_qubit_decompose.rs b/crates/accelerate/src/two_qubit_decompose.rs index be3a32255293..cb3c1308829e 100644 --- a/crates/accelerate/src/two_qubit_decompose.rs +++ b/crates/accelerate/src/two_qubit_decompose.rs @@ -351,9 +351,12 @@ enum Specialization { MirrorControlledEquiv, // These next 3 gates use the definition of fSim from eq (1) in: // https://arxiv.org/pdf/2001.08343.pdf - SimaabEquiv, - SimabbEquiv, - SimabmbEquiv, + #[allow(non_camel_case_types)] + fSimaabEquiv, + #[allow(non_camel_case_types)] + fSimabbEquiv, + #[allow(non_camel_case_types)] + fSimabmbEquiv, } impl Specialization { @@ -366,9 +369,9 @@ impl Specialization { Specialization::PartialSWAPFlipEquiv => 4, Specialization::ControlledEquiv => 5, Specialization::MirrorControlledEquiv => 6, - Specialization::SimaabEquiv => 7, - Specialization::SimabbEquiv => 8, - Specialization::SimabmbEquiv => 9, + Specialization::fSimaabEquiv => 7, + Specialization::fSimabbEquiv => 8, + Specialization::fSimabmbEquiv => 9, } } @@ -381,9 +384,9 @@ impl Specialization { 4 => Specialization::PartialSWAPFlipEquiv, 5 => Specialization::ControlledEquiv, 6 => Specialization::MirrorControlledEquiv, - 7 => Specialization::SimaabEquiv, - 8 => Specialization::SimabbEquiv, - 9 => Specialization::SimabmbEquiv, + 7 => Specialization::fSimaabEquiv, + 8 => Specialization::fSimabbEquiv, + 9 => Specialization::fSimabmbEquiv, _ => unreachable!("Invalid specialization value"), } } @@ -751,11 +754,11 @@ impl TwoQubitWeylDecomposition { } else if is_close(PI4, PI4, c) { Specialization::MirrorControlledEquiv } else if is_close((a + b) / 2., (a + b) / 2., c) { - Specialization::SimaabEquiv + Specialization::fSimaabEquiv } else if is_close(a, (b + c) / 2., (b + c) / 2.) { - Specialization::SimabbEquiv + Specialization::fSimabbEquiv } else if is_close(a, (b - c) / 2., (c - b) / 2.) { - Specialization::SimabmbEquiv + Specialization::fSimabmbEquiv } else { Specialization::General } @@ -924,7 +927,7 @@ impl TwoQubitWeylDecomposition { // This gate binds 5 parameters, we make it canonical by setting: // // :math:`K2_l = Ry(\theta_l)\cdot Rz(\lambda_l)`. - Specialization::SimaabEquiv => { + Specialization::fSimaabEquiv => { let [k2ltheta, k2lphi, k2llambda, k2lphase] = angles_from_unitary(general.K2l.view(), EulerBasis::ZYZ); TwoQubitWeylDecomposition { @@ -945,7 +948,7 @@ impl TwoQubitWeylDecomposition { // This gate binds 5 parameters, we make it canonical by setting: // // :math:`K2_l = Ry(\theta_l)Rx(\lambda_l)` - Specialization::SimabbEquiv => { + Specialization::fSimabbEquiv => { let euler_basis = EulerBasis::XYX; let [k2ltheta, k2lphi, k2llambda, k2lphase] = angles_from_unitary(general.K2l.view(), euler_basis); @@ -968,7 +971,7 @@ impl TwoQubitWeylDecomposition { // This gate binds 5 parameters, we make it canonical by setting: // // :math:`K2_l = Ry(\theta_l)Rx(\lambda_l)` - Specialization::SimabmbEquiv => { + Specialization::fSimabmbEquiv => { let euler_basis = EulerBasis::XYX; let [k2ltheta, k2lphi, k2llambda, k2lphase] = angles_from_unitary(general.K2l.view(), euler_basis); From d9fb03d40cceee8eb374202e82b84c70491de642 Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Thu, 14 Mar 2024 13:02:57 +0900 Subject: [PATCH 45/46] Update tests for correct specialization enum name --- test/python/synthesis/test_synthesis.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/python/synthesis/test_synthesis.py b/test/python/synthesis/test_synthesis.py index 864bbd334239..d3f8560cb9b2 100644 --- a/test/python/synthesis/test_synthesis.py +++ b/test/python/synthesis/test_synthesis.py @@ -903,7 +903,7 @@ def test_weyl_specialize_fsim_aab(self, aaa=0.456, bbb=0.132): self.check_two_qubit_weyl_specialization( k1 @ Ud(a + da, b + db, c + dc) @ k2, 0.999, - Specialization.SimaabEquiv, + Specialization.fSimaabEquiv, {"rz": 7, "ry": 4, "rxx": 1, "ryy": 1, "rzz": 1}, ) @@ -917,7 +917,7 @@ def test_weyl_specialize_fsim_abb(self, aaa=0.456, bbb=0.132): self.check_two_qubit_weyl_specialization( k1 @ Ud(a + da, b + db, c + dc) @ k2, 0.999, - Specialization.SimabbEquiv, + Specialization.fSimabbEquiv, {"rx": 7, "ry": 4, "rxx": 1, "ryy": 1, "rzz": 1}, ) @@ -931,7 +931,7 @@ def test_weyl_specialize_fsim_abmb(self, aaa=0.456, bbb=0.132): self.check_two_qubit_weyl_specialization( k1 @ Ud(a + da, b + db, c + dc) @ k2, 0.999, - Specialization.SimabmbEquiv, + Specialization.fSimabmbEquiv, {"rx": 7, "ry": 4, "rxx": 1, "ryy": 1, "rzz": 1}, ) From 4dda4beb6682fac945dfee2d5a6f6a7621a5e076 Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Thu, 14 Mar 2024 21:31:56 +0900 Subject: [PATCH 46/46] Expose specialization enum via private class attr and use for __repr__ --- qiskit/synthesis/two_qubit/two_qubit_decompose.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/qiskit/synthesis/two_qubit/two_qubit_decompose.py b/qiskit/synthesis/two_qubit/two_qubit_decompose.py index ac7e37fe5406..d8d6b0d8033d 100644 --- a/qiskit/synthesis/two_qubit/two_qubit_decompose.py +++ b/qiskit/synthesis/two_qubit/two_qubit_decompose.py @@ -151,6 +151,8 @@ class TwoQubitWeylDecomposition: requested_fidelity: Optional[float] # None means no automatic specialization calculated_fidelity: float # Fidelity after specialization + _specializations = two_qubit_decompose.Specialization + def __init__( self, unitary_matrix: np.ndarray, @@ -224,13 +226,15 @@ def __repr__(self): b64ascii[-1] += "," pretty = [f"# {x.rstrip()}" for x in str(self).splitlines()] indent = "\n" + " " * 4 + specialization_variant = str(self._inner_decomposition.specialization).split(".")[1] + specialization_repr = f"{type(self).__qualname__}._specializations.{specialization_variant}" lines = ( [prefix] + pretty + b64ascii + [ f"requested_fidelity={self.requested_fidelity},", - f"_specialization={self._inner_decomposition.specialization},", + f"_specialization={specialization_repr},", f"calculated_fidelity={self.calculated_fidelity},", f"actual_fidelity={self.actual_fidelity()},", f"abc={(self.a, self.b, self.c)})",