diff --git a/crates/accelerate/src/two_qubit_decompose.rs b/crates/accelerate/src/two_qubit_decompose.rs index 16ced4bf8a10..847a7069ab3f 100644 --- a/crates/accelerate/src/two_qubit_decompose.rs +++ b/crates/accelerate/src/two_qubit_decompose.rs @@ -18,15 +18,16 @@ // of real components and one of imaginary components. // In order to avoid copying we want to use `MatRef` or `MatMut`. -use approx::abs_diff_eq; +use approx::{abs_diff_eq, relative_eq}; use num_complex::{Complex, Complex64, ComplexFloat}; +use num_traits::Zero; use pyo3::exceptions::PyIndexError; use pyo3::import_exception; use pyo3::prelude::*; use pyo3::wrap_pyfunction; use pyo3::Python; use smallvec::{smallvec, SmallVec}; -use std::f64::consts::PI; +use std::f64::consts::{FRAC_1_SQRT_2, PI}; use faer::IntoFaerComplex; use faer::IntoNdarray; @@ -38,11 +39,11 @@ use ndarray::linalg::kron; use ndarray::prelude::*; use ndarray::Zip; use numpy::PyReadonlyArray2; -use numpy::ToPyArray; +use numpy::{IntoPyArray, ToPyArray}; use crate::euler_one_qubit_decomposer::{ angles_from_unitary, det_one_qubit, unitary_to_gate_sequence_inner, EulerBasis, - ANGLE_ZERO_EPSILON, + OneQubitGateSequence, ANGLE_ZERO_EPSILON, }; use crate::utils; @@ -166,6 +167,11 @@ impl Arg for c64 { } } +#[inline(always)] +fn transpose_conjugate(mat: ArrayView2) -> Array2 { + mat.t().mapv(|x| x.conj()) +} + 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) @@ -200,7 +206,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 r_t_conj: Array2 = transpose_conjugate(r.view()); let eye = aview2(&ONE_QUBIT_IDENTITY); let mut temp = kron(&eye, &r_t_conj); temp = special_unitary.dot(&temp); @@ -336,6 +342,17 @@ fn rz_matrix(theta: f64) -> Array2 { ] } +const HGATE: [[Complex64; 2]; 2] = [ + [ + Complex64::new(FRAC_1_SQRT_2, 0.), + Complex64::new(FRAC_1_SQRT_2, 0.), + ], + [ + Complex64::new(FRAC_1_SQRT_2, 0.), + Complex64::new(-FRAC_1_SQRT_2, 0.), + ], +]; + const DEFAULT_FIDELITY: f64 = 1.0 - 1.0e-9; const C1_IM: Complex64 = Complex64::new(0.0, 1.0); @@ -1159,6 +1176,686 @@ impl TwoQubitGateSequence { } } } +#[allow(non_snake_case)] +#[pyclass(module = "qiskit._accelerate.two_qubit_decompose", subclass)] +pub struct TwoQubitBasisDecomposer { + gate: String, + basis_fidelity: f64, + euler_basis: EulerBasis, + pulse_optimize: Option, + basis_decomposer: TwoQubitWeylDecomposition, + #[pyo3(get)] + super_controlled: bool, + u0l: Array2, + u0r: Array2, + u1l: Array2, + u1ra: Array2, + u1rb: Array2, + u2la: Array2, + u2lb: Array2, + u2ra: Array2, + u2rb: Array2, + u3l: Array2, + u3r: Array2, + q0l: Array2, + q0r: Array2, + q1la: Array2, + q1lb: Array2, + q1ra: Array2, + q1rb: Array2, + q2l: Array2, + q2r: Array2, +} +impl TwoQubitBasisDecomposer { + fn decomp1_inner( + &self, + target: &TwoQubitWeylDecomposition, + ) -> SmallVec<[Array2; 8]> { + // FIXME: fix for z!=0 and c!=0 using closest reflection (not always in the Weyl chamber) + smallvec![ + transpose_conjugate(self.basis_decomposer.K2r.view()).dot(&target.K2r), + transpose_conjugate(self.basis_decomposer.K2l.view()).dot(&target.K2l), + target + .K1r + .dot(&transpose_conjugate(self.basis_decomposer.K1r.view())), + target + .K1l + .dot(&transpose_conjugate(self.basis_decomposer.K1l.view())), + ] + } + + fn decomp2_supercontrolled_inner( + &self, + target: &TwoQubitWeylDecomposition, + ) -> SmallVec<[Array2; 8]> { + smallvec![ + self.q2r.dot(&target.K2r), + self.q2l.dot(&target.K2l), + self.q1ra.dot(&rz_matrix(2. * target.b)).dot(&self.q1rb), + self.q1la.dot(&rz_matrix(-2. * target.a)).dot(&self.q1lb), + target.K1r.dot(&self.q0r), + target.K1l.dot(&self.q0l), + ] + } + + fn decomp3_supercontrolled_inner( + &self, + target: &TwoQubitWeylDecomposition, + ) -> SmallVec<[Array2; 8]> { + smallvec![ + self.u3r.dot(&target.K2r), + self.u3l.dot(&target.K2l), + self.u2ra.dot(&rz_matrix(2. * target.b)).dot(&self.u2rb), + self.u2la.dot(&rz_matrix(-2. * target.a)).dot(&self.u2lb), + self.u1ra.dot(&rz_matrix(-2. * target.c)).dot(&self.u1rb), + self.u1l.clone(), + target.K1r.dot(&self.u0r), + target.K1l.dot(&self.u0l), + ] + } + + fn append_1q_sequence( + &self, + gates: &mut TwoQubitSequenceVec, + global_phase: &mut f64, + unitary: ArrayView2, + qubit: u8, + ) { + let target_1q_basis_list = vec![self.euler_basis]; + let sequence = unitary_to_gate_sequence_inner( + unitary, + &target_1q_basis_list, + qubit as usize, + None, + true, + None, + ); + if let Some(sequence) = sequence { + *global_phase += sequence.global_phase; + for gate in sequence.gates { + gates.push((gate.0, gate.1, smallvec![qubit])); + } + } + } + + fn pulse_optimal_chooser( + &self, + best_nbasis: u8, + decomposition: &SmallVec<[Array2; 8]>, + target_decomposed: &TwoQubitWeylDecomposition, + ) -> Option { + if self.pulse_optimize.is_some() + && (best_nbasis == 0 || best_nbasis == 1 || best_nbasis > 3) + { + return None; + } + match self.euler_basis { + EulerBasis::ZSX => (), + EulerBasis::ZSXX => (), + _ => return None, + } + if self.gate != "cx" { + return None; + } + let mut gates = Vec::new(); + let mut global_phase = target_decomposed.global_phase; + global_phase -= best_nbasis as f64 * self.basis_decomposer.global_phase; + if best_nbasis == 3 { + let atol = 1e-10; // absolute tolerance for floats + let euler_q0: Vec<[f64; 3]> = decomposition + .iter() + .step_by(2) + .map(|decomp| { + let euler_angles = angles_from_unitary(decomp.view(), EulerBasis::ZXZ); + global_phase += euler_angles[3]; + [euler_angles[1], euler_angles[2], euler_angles[0]] + }) + .collect(); + let euler_q1: Vec<[f64; 3]> = decomposition + .iter() + .skip(1) + .step_by(2) + .map(|decomp| { + let euler_angles = angles_from_unitary(decomp.view(), EulerBasis::XZX); + global_phase += euler_angles[3]; + [euler_angles[1], euler_angles[2], euler_angles[0]] + }) + .collect(); + let x12 = euler_q0[1][2] + euler_q0[2][0]; + let x12_is_non_zero = !abs_diff_eq!(x12, PI4, epsilon = atol); + let mut x12_is_old_mult = None; + let mut x12_phase = 0.; + let x12_is_pi_mult = abs_diff_eq!(x12.sin(), 0., epsilon = atol); + if x12_is_pi_mult { + x12_is_old_mult = Some(abs_diff_eq!(x12.cos(), -1., epsilon = atol)); + x12_phase = PI * x12.cos(); + } + let x02_add = x12 - euler_q0[1][0]; + let x12_is_half_pi = abs_diff_eq!(x12, PI2, epsilon = atol); + + let mut euler_matrix_q0 = rz_matrix(euler_q0[0][0]).dot(&rx_matrix(euler_q0[0][1])); + if x12_is_non_zero && x12_is_pi_mult { + euler_matrix_q0 = euler_matrix_q0.dot(&rz_matrix(euler_q0[0][2] - x02_add)); + } else { + euler_matrix_q0 = euler_matrix_q0.dot(&rz_matrix(euler_q0[0][2] + euler_q0[1][0])); + } + euler_matrix_q0 = euler_matrix_q0.dot(&aview2(&HGATE)); + self.append_1q_sequence(&mut gates, &mut global_phase, euler_matrix_q0.view(), 0); + let euler_matrix_q1 = rx_matrix(euler_q1[0][0]) + .dot(&rz_matrix(euler_q1[0][1])) + .dot(&rx_matrix(euler_q1[0][2] + euler_q1[1][0])); + self.append_1q_sequence(&mut gates, &mut global_phase, euler_matrix_q1.view(), 1); + gates.push(("cx".to_string(), smallvec![], smallvec![1, 0])); + if x12_is_pi_mult { + if x12_is_non_zero { + global_phase += x12_phase; + } + if x12_is_non_zero && x12_is_old_mult.unwrap() { + gates.push(("rz".to_string(), smallvec![-euler_q0[1][1]], smallvec![0])); + } else { + gates.push(("rz".to_string(), smallvec![euler_q0[1][1]], smallvec![0])); + global_phase += PI; + } + } + if x12_is_half_pi { + gates.push(("sx".to_string(), smallvec![], smallvec![0])); + global_phase -= PI4; + } else if x12_is_non_zero && !x12_is_pi_mult { + if self.pulse_optimize.is_none() { + self.append_1q_sequence( + &mut gates, + &mut global_phase, + rx_matrix(x12).view(), + 0, + ); + } else { + return None; + } + } + if abs_diff_eq!(euler_q1[1][1], PI2, epsilon = atol) { + gates.push(("sx".to_string(), smallvec![], smallvec![1])); + global_phase -= PI4 + } else if self.pulse_optimize.is_none() { + self.append_1q_sequence( + &mut gates, + &mut global_phase, + rx_matrix(euler_q1[1][1]).view(), + 1, + ); + } else { + return None; + } + gates.push(( + "rz".to_string(), + smallvec![euler_q1[1][2] + euler_q1[2][0]], + smallvec![1], + )); + gates.push(("cx".to_string(), smallvec![], smallvec![1, 0])); + gates.push(("rz".to_string(), smallvec![euler_q0[2][1]], smallvec![0])); + if abs_diff_eq!(euler_q1[2][1], PI2, epsilon = atol) { + gates.push(("sx".to_string(), smallvec![], smallvec![1])); + global_phase -= PI4; + } else if self.pulse_optimize.is_none() { + self.append_1q_sequence( + &mut gates, + &mut global_phase, + rx_matrix(euler_q1[2][1]).view(), + 1, + ); + } else { + return None; + } + gates.push(("cx".to_string(), smallvec![], smallvec![1, 0])); + let euler_matrix = aview2(&HGATE) + .dot(&rz_matrix(euler_q0[2][2] + euler_q0[3][0])) + .dot(&rx_matrix(euler_q0[3][1])) + .dot(&rz_matrix(euler_q0[3][2])); + self.append_1q_sequence(&mut gates, &mut global_phase, euler_matrix.view(), 0); + let euler_matrix = aview2(&HGATE) + .dot(&rx_matrix(euler_q1[2][2] + euler_q1[3][0])) + .dot(&rz_matrix(euler_q1[3][1])) + .dot(&rx_matrix(euler_q1[3][2])); + self.append_1q_sequence(&mut gates, &mut global_phase, euler_matrix.view(), 1); + // TODO: fix the sign problem to avoid correction here + // if cmath.isclose( + // target_decomposed.unitary_matrix[0, 0], -(Operator(qc).data[0,0]), abs_tol=atol + // ): + // qc.global_phase += math.pi + } else if best_nbasis == 2 { + let euler_q0: Vec<[f64; 3]> = decomposition + .iter() + .step_by(2) + .map(|decomp| { + let euler_angles = angles_from_unitary(decomp.view(), EulerBasis::ZXZ); + global_phase += euler_angles[3]; + [euler_angles[1], euler_angles[2], euler_angles[0]] + }) + .collect(); + let euler_q1: Vec<[f64; 3]> = decomposition + .iter() + .skip(1) + .step_by(2) + .map(|decomp| { + let euler_angles = angles_from_unitary(decomp.view(), EulerBasis::XZX); + global_phase += euler_angles[3]; + [euler_angles[1], euler_angles[2], euler_angles[0]] + }) + .collect(); + let euler_matrix_q0 = rz_matrix(euler_q0[0][0]) + .dot(&rx_matrix(euler_q0[0][1])) + .dot(&rz_matrix(euler_q0[0][2] + euler_q0[1][0] + PI2)); + self.append_1q_sequence(&mut gates, &mut global_phase, euler_matrix_q0.view(), 0); + let euler_matrix_q1 = rx_matrix(euler_q1[0][0]) + .dot(&rz_matrix(euler_q1[0][1])) + .dot(&rx_matrix(euler_q1[0][2] + euler_q1[1][0])); + self.append_1q_sequence(&mut gates, &mut global_phase, euler_matrix_q1.view(), 1); + gates.push(("cx".to_string(), smallvec![], smallvec![0, 1])); + gates.push(("sx".to_string(), smallvec![], smallvec![0])); + gates.push(( + "rz".to_string(), + smallvec![euler_q0[1][1] - PI], + smallvec![0], + )); + gates.push(("sx".to_string(), smallvec![], smallvec![0])); + gates.push(("rz".to_string(), smallvec![euler_q1[1][1]], smallvec![1])); + global_phase += PI2; + gates.push(("cx".to_string(), smallvec![], smallvec![0, 1])); + let euler_matrix_q0 = rz_matrix(euler_q0[1][2] + euler_q0[2][0] + PI2) + .dot(&rx_matrix(euler_q0[2][1])) + .dot(&rz_matrix(euler_q0[2][2])); + self.append_1q_sequence(&mut gates, &mut global_phase, euler_matrix_q0.view(), 0); + let euler_matrix_q1 = rx_matrix(euler_q1[1][2] + euler_q1[2][0]) + .dot(&rz_matrix(euler_q1[2][1])) + .dot(&rx_matrix(euler_q1[2][2])); + self.append_1q_sequence(&mut gates, &mut global_phase, euler_matrix_q1.view(), 1); + } else { + return None; + } + + Some(TwoQubitGateSequence { + gates, + global_phase, + }) + } +} + +fn decomp0_inner(target: &TwoQubitWeylDecomposition) -> SmallVec<[Array2; 8]> { + smallvec![target.K1r.dot(&target.K2r), target.K1l.dot(&target.K2l),] +} + +#[pymethods] +impl TwoQubitBasisDecomposer { + fn __getnewargs__(&self, py: Python) -> (String, PyObject, f64, String, Option) { + ( + self.gate.clone(), + self.basis_decomposer.unitary_matrix.to_pyarray(py).into(), + self.basis_fidelity, + self.euler_basis.to_str(), + self.pulse_optimize, + ) + } + + #[new] + #[pyo3(signature=(gate, gate_matrix, basis_fidelity=1.0, euler_basis="U", pulse_optimize=None))] + fn new( + gate: String, + gate_matrix: PyReadonlyArray2, + basis_fidelity: f64, + euler_basis: &str, + pulse_optimize: Option, + ) -> PyResult { + let ipz: ArrayView2 = aview2(&IPZ); + let basis_decomposer = + TwoQubitWeylDecomposition::new(gate_matrix, Some(DEFAULT_FIDELITY), None, false)?; + let super_controlled = relative_eq!(basis_decomposer.a, PI4, max_relative = 1e-09) + && relative_eq!(basis_decomposer.c, 0.0, max_relative = 1e-09); + + // Create some useful matrices U1, U2, U3 are equivalent to the basis, + // expand as Ui = Ki1.Ubasis.Ki2 + let b = basis_decomposer.b; + let temp = Complex64::new(0.5, -0.5); + let k11l = array![ + [ + temp * (Complex64::new(0., -1.) * Complex64::new(0., -b).exp()), + temp * Complex64::new(0., -b).exp() + ], + [ + temp * (Complex64::new(0., -1.) * Complex64::new(0., b).exp()), + temp * -(Complex64::new(0., b).exp()) + ], + ]; + let k11r = array![ + [ + FRAC_1_SQRT_2 * (Complex64::new(0., 1.) * Complex64::new(0., -b).exp()), + FRAC_1_SQRT_2 * -Complex64::new(0., -b).exp() + ], + [ + FRAC_1_SQRT_2 * Complex64::new(0., b).exp(), + FRAC_1_SQRT_2 * (Complex64::new(0., -1.) * Complex64::new(0., b).exp()) + ], + ]; + let k12l = array![ + [Complex64::new(0.5, 0.5), Complex64::new(0.5, 0.5)], + [Complex64::new(-0.5, 0.5), Complex64::new(0.5, -0.5)], + ]; + let k12r = array![ + [ + Complex64::new(0., FRAC_1_SQRT_2), + Complex64::new(FRAC_1_SQRT_2, 0.) + ], + [ + Complex64::new(-FRAC_1_SQRT_2, 0.), + Complex64::new(0., -FRAC_1_SQRT_2) + ], + ]; + let k32l_k21l = array![ + [ + FRAC_1_SQRT_2 * Complex64::new(1., (2. * b).cos()), + FRAC_1_SQRT_2 * (Complex64::new(0., 1.) * (2. * b).sin()) + ], + [ + FRAC_1_SQRT_2 * (Complex64::new(0., 1.) * (2. * b).sin()), + FRAC_1_SQRT_2 * Complex64::new(1., -(2. * b).cos()) + ], + ]; + let temp = Complex64::new(0.5, 0.5); + let k21r = array![ + [ + temp * (Complex64::new(0., -1.) * Complex64::new(0., -2. * b).exp()), + temp * Complex64::new(0., -2. * b).exp() + ], + [ + temp * (Complex64::new(0., 1.) * Complex64::new(0., 2. * b).exp()), + temp * Complex64::new(0., 2. * b).exp() + ], + ]; + let k22l = array![ + [ + Complex64::new(FRAC_1_SQRT_2, 0.), + Complex64::new(-FRAC_1_SQRT_2, 0.) + ], + [ + Complex64::new(FRAC_1_SQRT_2, 0.), + Complex64::new(FRAC_1_SQRT_2, 0.) + ], + ]; + let k22r = array![ + [Complex64::zero(), Complex64::new(1., 0.)], + [Complex64::new(-1., 0.), Complex64::zero()], + ]; + let k31l = array![ + [ + FRAC_1_SQRT_2 * Complex64::new(0., -b).exp(), + FRAC_1_SQRT_2 * Complex64::new(0., -b).exp() + ], + [ + FRAC_1_SQRT_2 * -Complex64::new(0., b).exp(), + FRAC_1_SQRT_2 * Complex64::new(0., b).exp() + ], + ]; + let temp = Complex64::new(0., 1.); + let k31r = array![ + [temp * Complex64::new(0., b).exp(), Complex64::zero()], + [Complex64::zero(), temp * -Complex64::new(0., -b).exp()], + ]; + let temp = Complex64::new(0.5, 0.5); + let k32r = array![ + [ + temp * Complex64::new(0., b).exp(), + temp * -Complex64::new(0., -b).exp() + ], + [ + temp * (Complex64::new(0., -1.) * Complex64::new(0., b).exp()), + temp * (Complex64::new(0., -1.) * Complex64::new(0., -b).exp()) + ], + ]; + let k1ld = transpose_conjugate(basis_decomposer.K1l.view()); + let k1rd = transpose_conjugate(basis_decomposer.K1r.view()); + let k2ld = transpose_conjugate(basis_decomposer.K2l.view()); + let k2rd = transpose_conjugate(basis_decomposer.K2r.view()); + // Pre-build the fixed parts of the matrices used in 3-part decomposition + let u0l = k31l.dot(&k1ld); + let u0r = k31r.dot(&k1rd); + let u1l = k2ld.dot(&k32l_k21l).dot(&k1ld); + let u1ra = k2rd.dot(&k32r); + let u1rb = k21r.dot(&k1rd); + let u2la = k2ld.dot(&k22l); + let u2lb = k11l.dot(&k1ld); + let u2ra = k2rd.dot(&k22r); + let u2rb = k11r.dot(&k1rd); + let u3l = k2ld.dot(&k12l); + let u3r = k2rd.dot(&k12r); + // Pre-build the fixed parts of the matrices used in the 2-part decomposition + let q0l = transpose_conjugate(k12l.view()).dot(&k1ld); + let q0r = transpose_conjugate(k12r.view()).dot(&ipz).dot(&k1rd); + let q1la = k2ld.dot(&transpose_conjugate(k11l.view())); + let q1lb = k11l.dot(&k1ld); + let q1ra = k2rd.dot(&ipz).dot(&transpose_conjugate(k11r.view())); + let q1rb = k11r.dot(&k1rd); + let q2l = k2ld.dot(&k12l); + let q2r = k2rd.dot(&k12r); + + Ok(TwoQubitBasisDecomposer { + gate, + basis_fidelity, + euler_basis: EulerBasis::from_string(euler_basis)?, + pulse_optimize, + basis_decomposer, + super_controlled, + u0l, + u0r, + u1l, + u1ra, + u1rb, + u2la, + u2lb, + u2ra, + u2rb, + u3l, + u3r, + q0l, + q0r, + q1la, + q1lb, + q1ra, + q1rb, + q2l, + q2r, + }) + } + + fn traces(&self, target: &TwoQubitWeylDecomposition) -> [Complex64; 4] { + [ + 4. * Complex64::new( + target.a.cos() * target.b.cos() * target.c.cos(), + target.a.sin() * target.b.sin() * target.c.sin(), + ), + 4. * Complex64::new( + (PI4 - target.a).cos() + * (self.basis_decomposer.b - target.b).cos() + * target.c.cos(), + (PI4 - target.a).sin() + * (self.basis_decomposer.b - target.b).sin() + * target.c.sin(), + ), + Complex64::new(4. * target.c.cos(), 0.), + Complex64::new(4., 0.), + ] + } + + /// Decompose target :math:`\sim U_d(x, y, z)` with :math:`0` uses of the basis gate. + /// Result :math:`U_r` has trace: + /// + /// .. math:: + /// + /// \Big\vert\text{Tr}(U_r\cdot U_\text{target}^{\dag})\Big\vert = + /// 4\Big\vert (\cos(x)\cos(y)\cos(z)+ j \sin(x)\sin(y)\sin(z)\Big\vert + /// + /// which is optimal for all targets and bases + #[staticmethod] + fn decomp0(py: Python, target: &TwoQubitWeylDecomposition) -> SmallVec<[PyObject; 2]> { + decomp0_inner(target) + .into_iter() + .map(|x| x.into_pyarray(py).into()) + .collect() + } + + /// Decompose target :math:`\sim U_d(x, y, z)` with :math:`1` use of the basis gate + /// math:`\sim U_d(a, b, c)`. + /// Result :math:`U_r` has trace: + /// + /// .. math:: + /// + /// \Big\vert\text{Tr}(U_r \cdot U_\text{target}^{\dag})\Big\vert = + /// 4\Big\vert \cos(x-a)\cos(y-b)\cos(z-c) + j \sin(x-a)\sin(y-b)\sin(z-c)\Big\vert + /// + /// which is optimal for all targets and bases with ``z==0`` or ``c==0``. + fn decomp1(&self, py: Python, target: &TwoQubitWeylDecomposition) -> SmallVec<[PyObject; 4]> { + self.decomp1_inner(target) + .into_iter() + .map(|x| x.into_pyarray(py).into()) + .collect() + } + + /// Decompose target :math:`\sim U_d(x, y, z)` with :math:`2` uses of the basis gate. + /// + /// For supercontrolled basis :math:`\sim U_d(\pi/4, b, 0)`, all b, result :math:`U_r` has trace + /// + /// .. math:: + /// + /// \Big\vert\text{Tr}(U_r \cdot U_\text{target}^\dag) \Big\vert = 4\cos(z) + /// + /// which is the optimal approximation for basis of CNOT-class :math:`\sim U_d(\pi/4, 0, 0)` + /// or DCNOT-class :math:`\sim U_d(\pi/4, \pi/4, 0)` and any target. It may + /// be sub-optimal for :math:`b \neq 0` (i.e. there exists an exact decomposition for any target + /// using :math:`B \sim U_d(\pi/4, \pi/8, 0)`, but it may not be this decomposition). + /// This is an exact decomposition for supercontrolled basis and target :math:`\sim U_d(x, y, 0)`. + /// No guarantees for non-supercontrolled basis. + fn decomp2_supercontrolled( + &self, + py: Python, + target: &TwoQubitWeylDecomposition, + ) -> SmallVec<[PyObject; 6]> { + self.decomp2_supercontrolled_inner(target) + .into_iter() + .map(|x| x.into_pyarray(py).into()) + .collect() + } + + /// Decompose target with :math:`3` uses of the basis. + /// + /// This is an exact decomposition for supercontrolled basis :math:`\sim U_d(\pi/4, b, 0)`, all b, + /// and any target. No guarantees for non-supercontrolled basis. + fn decomp3_supercontrolled( + &self, + py: Python, + target: &TwoQubitWeylDecomposition, + ) -> SmallVec<[PyObject; 8]> { + self.decomp3_supercontrolled_inner(target) + .into_iter() + .map(|x| x.into_pyarray(py).into()) + .collect() + } + + /// Decompose a two-qubit ``unitary`` over fixed basis and :math:`SU(2)` using the best + /// approximation given that each basis application has a finite ``basis_fidelity``. + #[pyo3(signature = (unitary, basis_fidelity=None, approximate=true, _num_basis_uses=None))] + fn __call__( + &self, + unitary: PyReadonlyArray2, + basis_fidelity: Option, + approximate: bool, + _num_basis_uses: Option, + ) -> PyResult { + let basis_fidelity = if !approximate { + 1.0 + } else { + basis_fidelity.unwrap_or(self.basis_fidelity) + }; + let target_decomposed = + TwoQubitWeylDecomposition::new(unitary, Some(DEFAULT_FIDELITY), None, false)?; + let traces = self.traces(&target_decomposed); + let best_nbasis = traces + .into_iter() + .enumerate() + .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; + let best_nbasis = _num_basis_uses.unwrap_or(best_nbasis as u8); + let decomposition = match best_nbasis { + 0 => decomp0_inner(&target_decomposed), + 1 => self.decomp1_inner(&target_decomposed), + 2 => self.decomp2_supercontrolled_inner(&target_decomposed), + 3 => self.decomp3_supercontrolled_inner(&target_decomposed), + _ => unreachable!("Invalid basis to use"), + }; + let pulse_optimize = self.pulse_optimize.unwrap_or(true); + let sequence = if pulse_optimize { + self.pulse_optimal_chooser(best_nbasis, &decomposition, &target_decomposed) + } else { + None + }; + if let Some(seq) = sequence { + return Ok(seq); + } + let target_1q_basis_list = vec![self.euler_basis]; + let euler_decompositions: SmallVec<[Option; 8]> = decomposition + .iter() + .map(|decomp| { + unitary_to_gate_sequence_inner( + decomp.view(), + &target_1q_basis_list, + 0, + None, + true, + None, + ) + }) + .collect(); + let mut gates = Vec::new(); + let mut global_phase = target_decomposed.global_phase; + global_phase -= best_nbasis as f64 * self.basis_decomposer.global_phase; + if best_nbasis == 2 { + global_phase += PI; + } + for i in 0..best_nbasis as usize { + if let Some(euler_decomp) = &euler_decompositions[2 * i] { + for gate in &euler_decomp.gates { + gates.push((gate.0.clone(), gate.1.clone(), smallvec![0])); + } + global_phase += euler_decomp.global_phase + } + if let Some(euler_decomp) = &euler_decompositions[2 * i + 1] { + for gate in &euler_decomp.gates { + gates.push((gate.0.clone(), gate.1.clone(), smallvec![1])); + } + global_phase += euler_decomp.global_phase + } + gates.push((self.gate.clone(), smallvec![], smallvec![0, 1])); + } + if let Some(euler_decomp) = &euler_decompositions[2 * best_nbasis as usize] { + for gate in &euler_decomp.gates { + gates.push((gate.0.clone(), gate.1.clone(), smallvec![0])); + } + global_phase += euler_decomp.global_phase + } + if let Some(euler_decomp) = &euler_decompositions[2 * best_nbasis as usize + 1] { + for gate in &euler_decomp.gates { + gates.push((gate.0.clone(), gate.1.clone(), smallvec![1])); + } + global_phase += euler_decomp.global_phase + } + Ok(TwoQubitGateSequence { + gates, + global_phase, + }) + } + + fn num_basis_gates(&self, unitary: PyReadonlyArray2) -> usize { + _num_basis_gates(self.basis_decomposer.b, self.basis_fidelity, unitary) + } +} #[pymodule] pub fn two_qubit_decompose(_py: Python, m: &PyModule) -> PyResult<()> { @@ -1166,5 +1863,6 @@ pub fn two_qubit_decompose(_py: Python, m: &PyModule) -> PyResult<()> { 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 41e2d18bb4ed..c848033b9cce 100644 --- a/qiskit/synthesis/two_qubit/two_qubit_decompose.py +++ b/qiskit/synthesis/two_qubit/two_qubit_decompose.py @@ -36,7 +36,7 @@ import numpy as np from qiskit.circuit import QuantumRegister, QuantumCircuit, Gate -from qiskit.circuit.library.standard_gates import CXGate, RXGate +from qiskit.circuit.library.standard_gates import CXGate, U3Gate, U2Gate, U1Gate from qiskit.exceptions import QiskitError from qiskit.quantum_info.operators import Operator from qiskit.synthesis.one_qubit.one_qubit_decompose import ( @@ -494,234 +494,25 @@ def __init__( self.basis_fidelity = basis_fidelity self.pulse_optimize = pulse_optimize - basis = self.basis = TwoQubitWeylDecomposition(Operator(gate).data) - self._decomposer1q = OneQubitEulerDecomposer(euler_basis) - - # FIXME: find good tolerances - self.is_supercontrolled = math.isclose(basis.a, np.pi / 4) and math.isclose(basis.c, 0.0) - - # Create some useful matrices U1, U2, U3 are equivalent to the basis, - # expand as Ui = Ki1.Ubasis.Ki2 - b = basis.b - K11l = ( - 1 - / (1 + 1j) - * np.array( - [ - [-1j * cmath.exp(-1j * b), cmath.exp(-1j * b)], - [-1j * cmath.exp(1j * b), -cmath.exp(1j * b)], - ], - dtype=complex, - ) - ) - K11r = ( - 1 - / math.sqrt(2) - * np.array( - [ - [1j * cmath.exp(-1j * b), -cmath.exp(-1j * b)], - [cmath.exp(1j * b), -1j * cmath.exp(1j * b)], - ], - dtype=complex, - ) - ) - K12l = 1 / (1 + 1j) * np.array([[1j, 1j], [-1, 1]], dtype=complex) - K12r = 1 / math.sqrt(2) * np.array([[1j, 1], [-1, -1j]], dtype=complex) - K32lK21l = ( - 1 - / math.sqrt(2) - * np.array( - [ - [1 + 1j * np.cos(2 * b), 1j * np.sin(2 * b)], - [1j * np.sin(2 * b), 1 - 1j * np.cos(2 * b)], - ], - dtype=complex, - ) - ) - K21r = ( - 1 - / (1 - 1j) - * np.array( - [ - [-1j * cmath.exp(-2j * b), cmath.exp(-2j * b)], - [1j * cmath.exp(2j * b), cmath.exp(2j * b)], - ], - dtype=complex, - ) - ) - K22l = 1 / math.sqrt(2) * np.array([[1, -1], [1, 1]], dtype=complex) - K22r = np.array([[0, 1], [-1, 0]], dtype=complex) - K31l = ( - 1 - / math.sqrt(2) - * np.array( - [[cmath.exp(-1j * b), cmath.exp(-1j * b)], [-cmath.exp(1j * b), cmath.exp(1j * b)]], - dtype=complex, - ) + self._inner_decomposer = two_qubit_decompose.TwoQubitBasisDecomposer( + "USER_GATE", + Operator(gate).data, + basis_fidelity=basis_fidelity, + euler_basis=euler_basis, + pulse_optimize=pulse_optimize, ) - K31r = 1j * np.array([[cmath.exp(1j * b), 0], [0, -cmath.exp(-1j * b)]], dtype=complex) - K32r = ( - 1 - / (1 - 1j) - * np.array( - [ - [cmath.exp(1j * b), -cmath.exp(-1j * b)], - [-1j * cmath.exp(1j * b), -1j * cmath.exp(-1j * b)], - ], - dtype=complex, - ) - ) - k1ld = basis.K1l.T.conj() - k1rd = basis.K1r.T.conj() - k2ld = basis.K2l.T.conj() - k2rd = basis.K2r.T.conj() - - # Pre-build the fixed parts of the matrices used in 3-part decomposition - self.u0l = K31l.dot(k1ld) - self.u0r = K31r.dot(k1rd) - self.u1l = k2ld.dot(K32lK21l).dot(k1ld) - self.u1ra = k2rd.dot(K32r) - self.u1rb = K21r.dot(k1rd) - self.u2la = k2ld.dot(K22l) - self.u2lb = K11l.dot(k1ld) - self.u2ra = k2rd.dot(K22r) - self.u2rb = K11r.dot(k1rd) - self.u3l = k2ld.dot(K12l) - self.u3r = k2rd.dot(K12r) - - # Pre-build the fixed parts of the matrices used in the 2-part decomposition - self.q0l = K12l.T.conj().dot(k1ld) - self.q0r = K12r.T.conj().dot(_ipz).dot(k1rd) - self.q1la = k2ld.dot(K11l.T.conj()) - self.q1lb = K11l.dot(k1ld) - self.q1ra = k2rd.dot(_ipz).dot(K11r.T.conj()) - self.q1rb = K11r.dot(k1rd) - self.q2l = k2ld.dot(K12l) - self.q2r = k2rd.dot(K12r) - - # Decomposition into different number of gates - # In the future could use different decomposition functions for different basis classes, etc + self.is_supercontrolled = self._inner_decomposer.super_controlled if not self.is_supercontrolled: warnings.warn( - "Only know how to decompose properly for supercontrolled basis gate. " - "This gate is ~Ud({}, {}, {})".format(basis.a, basis.b, basis.c), + "Only know how to decompose properly for a supercontrolled basis gate.", stacklevel=2, ) - self.decomposition_fns = [ - self.decomp0, - self.decomp1, - self.decomp2_supercontrolled, - self.decomp3_supercontrolled, - ] - self._rqc = None - - def traces(self, target): - r""" - Give the expected traces :math:`\Big\vert\text{Tr}(U \cdot U_\text{target}^{\dag})\Big\vert` - for a different number of basis gates. - """ - # Future gotcha: extending this to non-supercontrolled basis. - # Careful: closest distance between a1,b1,c1 and a2,b2,c2 may be between reflections. - # This doesn't come up if either c1==0 or c2==0 but otherwise be careful. - ta, tb, tc = target.a, target.b, target.c - bb = self.basis.b - return [ - 4 - * complex( - math.cos(ta) * math.cos(tb) * math.cos(tc), - math.sin(ta) * math.sin(tb) * math.sin(tc), - ), - 4 - * complex( - math.cos(math.pi / 4 - ta) * math.cos(bb - tb) * math.cos(tc), - math.sin(math.pi / 4 - ta) * math.sin(bb - tb) * math.sin(tc), - ), - 4 * math.cos(tc), - 4, - ] - - @staticmethod - def decomp0(target): - r""" - Decompose target :math:`\sim U_d(x, y, z)` with :math:`0` uses of the basis gate. - Result :math:`U_r` has trace: - - .. math:: - - \Big\vert\text{Tr}(U_r\cdot U_\text{target}^{\dag})\Big\vert = - 4\Big\vert (\cos(x)\cos(y)\cos(z)+ j \sin(x)\sin(y)\sin(z)\Big\vert - - which is optimal for all targets and bases - """ - - U0l = target.K1l.dot(target.K2l) - U0r = target.K1r.dot(target.K2r) - return U0r, U0l - - def decomp1(self, target): - r"""Decompose target :math:`\sim U_d(x, y, z)` with :math:`1` use of the basis gate - :math:`\sim U_d(a, b, c)`. - Result :math:`U_r` has trace: - - .. math:: - - \Big\vert\text{Tr}(U_r \cdot U_\text{target}^{\dag})\Big\vert = - 4\Big\vert \cos(x-a)\cos(y-b)\cos(z-c) + j \sin(x-a)\sin(y-b)\sin(z-c)\Big\vert - - which is optimal for all targets and bases with ``z==0`` or ``c==0``. - """ - # FIXME: fix for z!=0 and c!=0 using closest reflection (not always in the Weyl chamber) - U0l = target.K1l.dot(self.basis.K1l.T.conj()) - U0r = target.K1r.dot(self.basis.K1r.T.conj()) - U1l = self.basis.K2l.T.conj().dot(target.K2l) - U1r = self.basis.K2r.T.conj().dot(target.K2r) - - return U1r, U1l, U0r, U0l - - def decomp2_supercontrolled(self, target): - r""" - Decompose target :math:`\sim U_d(x, y, z)` with :math:`2` uses of the basis gate. - - For supercontrolled basis :math:`\sim U_d(\pi/4, b, 0)`, all b, result :math:`U_r` has trace - - .. math:: - - \Big\vert\text{Tr}(U_r \cdot U_\text{target}^\dag) \Big\vert = 4\cos(z) - - which is the optimal approximation for basis of CNOT-class :math:`\sim U_d(\pi/4, 0, 0)` - or DCNOT-class :math:`\sim U_d(\pi/4, \pi/4, 0)` and any target. It may - be sub-optimal for :math:`b \neq 0` (i.e. there exists an exact decomposition for any target - using :math:`B \sim U_d(\pi/4, \pi/8, 0)`, but it may not be this decomposition). - This is an exact decomposition for supercontrolled basis and target :math:`\sim U_d(x, y, 0)`. - No guarantees for non-supercontrolled basis. - """ - - U0l = target.K1l.dot(self.q0l) - U0r = target.K1r.dot(self.q0r) - U1l = self.q1la.dot(rz_array(-2 * target.a)).dot(self.q1lb) - U1r = self.q1ra.dot(rz_array(2 * target.b)).dot(self.q1rb) - U2l = self.q2l.dot(target.K2l) - U2r = self.q2r.dot(target.K2r) - - return U2r, U2l, U1r, U1l, U0r, U0l - - def decomp3_supercontrolled(self, target): - r""" - Decompose target with :math:`3` uses of the basis. - This is an exact decomposition for supercontrolled basis :math:`\sim U_d(\pi/4, b, 0)`, all b, - and any target. No guarantees for non-supercontrolled basis. - """ - - U0l = target.K1l.dot(self.u0l) - U0r = target.K1r.dot(self.u0r) - U1l = self.u1l - U1r = self.u1ra.dot(rz_array(-2 * target.c)).dot(self.u1rb) - U2l = self.u2la.dot(rz_array(-2 * target.a)).dot(self.u2lb) - U2r = self.u2ra.dot(rz_array(2 * target.b)).dot(self.u2rb) - U3l = self.u3l.dot(target.K2l) - U3r = self.u3r.dot(target.K2r) - - return U3r, U3l, U2r, U2l, U1r, U1l, U0r, U0l + self.num_basis_gates = self._inner_decomposer.num_basis_gates + self.traces = self._inner_decomposer.traces + self.decomp0 = staticmethod(self._inner_decomposer.decomp0) + self.decomp1 = self._inner_decomposer.decomp1 + self.decomp2_supercontrolled = self._inner_decomposer.decomp2_supercontrolled + self.decomp3_supercontrolled = self._inner_decomposer.decomp3_supercontrolled def __call__( self, @@ -747,313 +538,34 @@ def __call__( Raises: QiskitError: if ``pulse_optimize`` is True but we don't know how to do it. """ - basis_fidelity = basis_fidelity or self.basis_fidelity - if approximate is False: - basis_fidelity = 1.0 - unitary = np.asarray(unitary, dtype=complex) - - target_decomposed = TwoQubitWeylDecomposition(unitary) - traces = self.traces(target_decomposed) - expected_fidelities = [trace_to_fid(traces[i]) * basis_fidelity**i for i in range(4)] - - best_nbasis = int(np.argmax(expected_fidelities)) - if _num_basis_uses is not None: - best_nbasis = _num_basis_uses - decomposition = self.decomposition_fns[best_nbasis](target_decomposed) - - # attempt pulse optimal decomposition - try: - if self.pulse_optimize in {None, True}: - return_circuit = self._pulse_optimal_chooser( - best_nbasis, decomposition, target_decomposed - ) - if return_circuit: - return return_circuit - except QiskitError: - if self.pulse_optimize: - raise - # do default decomposition + sequence = self._inner_decomposer( + np.asarray(unitary, dtype=complex), + basis_fidelity, + approximate, + _num_basis_uses=_num_basis_uses, + ) q = QuantumRegister(2) - decomposition_euler = [self._decomposer1q._decompose(x) for x in decomposition] - return_circuit = QuantumCircuit(q) - return_circuit.global_phase = target_decomposed.global_phase - return_circuit.global_phase -= best_nbasis * self.basis.global_phase - if best_nbasis == 2: - return_circuit.global_phase += np.pi - for i in range(best_nbasis): - return_circuit.compose(decomposition_euler[2 * i], [q[0]], inplace=True) - return_circuit.compose(decomposition_euler[2 * i + 1], [q[1]], inplace=True) - return_circuit.append(self.gate, [q[0], q[1]]) - return_circuit.compose(decomposition_euler[2 * best_nbasis], [q[0]], inplace=True) - return_circuit.compose(decomposition_euler[2 * best_nbasis + 1], [q[1]], inplace=True) - return return_circuit - - def _pulse_optimal_chooser( - self, best_nbasis, decomposition, target_decomposed - ) -> QuantumCircuit: - """Determine method to find pulse optimal circuit. This method may be - removed once a more general approach is used. - - Returns: - QuantumCircuit: pulse optimal quantum circuit. - None: Probably ``nbasis==1`` and original circuit is fine. - - Raises: - QiskitError: Decomposition for selected basis not implemented. - """ - circuit = None - if self.pulse_optimize and best_nbasis in {0, 1}: - # already pulse optimal - return None - elif self.pulse_optimize and best_nbasis > 3: - raise QiskitError( - f"Unexpected number of entangling gates ({best_nbasis}) in decomposition." - ) - if self._decomposer1q.basis in {"ZSX", "ZSXX"}: - if isinstance(self.gate, CXGate): - if best_nbasis == 3: - circuit = self._get_sx_vz_3cx_efficient_euler(decomposition, target_decomposed) - elif best_nbasis == 2: - circuit = self._get_sx_vz_2cx_efficient_euler(decomposition, target_decomposed) - else: - raise QiskitError("pulse_optimizer currently only works with CNOT entangling gate") - else: - raise QiskitError( - '"pulse_optimize" currently only works with ZSX basis ' - f"({self._decomposer1q.basis} used)" - ) - return circuit - - def _get_sx_vz_2cx_efficient_euler(self, decomposition, target_decomposed): - """ - Decomposition of SU(4) gate for device with SX, virtual RZ, and CNOT gates assuming - two CNOT gates are needed. - - This first decomposes each unitary from the KAK decomposition into ZXZ on the source - qubit of the CNOTs and XZX on the targets in order to commute operators to beginning and - end of decomposition. The beginning and ending single qubit gates are then - collapsed and re-decomposed with the single qubit decomposer. This last step could be avoided - if performance is a concern. - """ - best_nbasis = 2 # by assumption - num_1q_uni = len(decomposition) - # list of euler angle decompositions on qubits 0 and 1 - euler_q0 = np.empty((num_1q_uni // 2, 3), dtype=float) - euler_q1 = np.empty((num_1q_uni // 2, 3), dtype=float) - global_phase = 0.0 - - # decompose source unitaries to zxz - zxz_decomposer = OneQubitEulerDecomposer("ZXZ") - for iqubit, decomp in enumerate(decomposition[0::2]): - euler_angles = zxz_decomposer.angles_and_phase(decomp) - euler_q0[iqubit, [1, 2, 0]] = euler_angles[:3] - global_phase += euler_angles[3] - # decompose target unitaries to xzx - xzx_decomposer = OneQubitEulerDecomposer("XZX") - for iqubit, decomp in enumerate(decomposition[1::2]): - euler_angles = xzx_decomposer.angles_and_phase(decomp) - euler_q1[iqubit, [1, 2, 0]] = euler_angles[:3] - global_phase += euler_angles[3] - qc = QuantumCircuit(2) - qc.global_phase = target_decomposed.global_phase - qc.global_phase -= best_nbasis * self.basis.global_phase - qc.global_phase += global_phase - - # TODO: make this more effecient to avoid double decomposition - # prepare beginning 0th qubit local unitary - circ = QuantumCircuit(1) - circ.rz(euler_q0[0][0], 0) - circ.rx(euler_q0[0][1], 0) - circ.rz(euler_q0[0][2] + euler_q0[1][0] + math.pi / 2, 0) - # re-decompose to basis of 1q decomposer - qceuler = self._decomposer1q(Operator(circ).data) - qc.compose(qceuler, [0], inplace=True) - - # prepare beginning 1st qubit local unitary - circ = QuantumCircuit(1) - circ.rx(euler_q1[0][0], 0) - circ.rz(euler_q1[0][1], 0) - circ.rx(euler_q1[0][2] + euler_q1[1][0], 0) - qceuler = self._decomposer1q(Operator(circ).data) - qc.compose(qceuler, [1], inplace=True) - - qc.cx(0, 1) - # the central decompositions are dependent on the specific form of the - # unitaries coming out of the two qubit decomposer which have some flexibility - # of choice. - qc.sx(0) - qc.rz(euler_q0[1][1] - math.pi, 0) - qc.sx(0) - qc.rz(euler_q1[1][1], 1) - qc.global_phase += math.pi / 2 - - qc.cx(0, 1) - - circ = QuantumCircuit(1) - circ.rz(euler_q0[1][2] + euler_q0[2][0] + math.pi / 2, 0) - circ.rx(euler_q0[2][1], 0) - circ.rz(euler_q0[2][2], 0) - qceuler = self._decomposer1q(Operator(circ).data) - qc.compose(qceuler, [0], inplace=True) - circ = QuantumCircuit(1) - circ.rx(euler_q1[1][2] + euler_q1[2][0], 0) - circ.rz(euler_q1[2][1], 0) - circ.rx(euler_q1[2][2], 0) - qceuler = self._decomposer1q(Operator(circ).data) - qc.compose(qceuler, [1], inplace=True) - - return qc - - def _get_sx_vz_3cx_efficient_euler(self, decomposition, target_decomposed): - """ - Decomposition of SU(4) gate for device with SX, virtual RZ, and CNOT gates assuming - three CNOT gates are needed. - - This first decomposes each unitary from the KAK decomposition into ZXZ on the source - qubit of the CNOTs and XZX on the targets in order commute operators to beginning and - end of decomposition. Inserting Hadamards reverses the direction of the CNOTs and transforms - a variable Rx -> variable virtual Rz. The beginning and ending single qubit gates are then - collapsed and re-decomposed with the single qubit decomposer. This last step could be avoided - if performance is a concern. - """ - best_nbasis = 3 # by assumption - num_1q_uni = len(decomposition) - # create structure to hold euler angles: 1st index represents unitary "group" wrt cx - # 2nd index represents index of euler triple. - euler_q0 = np.empty((num_1q_uni // 2, 3), dtype=float) - euler_q1 = np.empty((num_1q_uni // 2, 3), dtype=float) - global_phase = 0.0 - atol = 1e-10 # absolute tolerance for floats - - # decompose source unitaries to zxz - zxz_decomposer = OneQubitEulerDecomposer("ZXZ") - for iqubit, decomp in enumerate(decomposition[0::2]): - euler_angles = zxz_decomposer.angles_and_phase(decomp) - euler_q0[iqubit, [1, 2, 0]] = euler_angles[:3] - global_phase += euler_angles[3] - # decompose target unitaries to xzx - xzx_decomposer = OneQubitEulerDecomposer("XZX") - for iqubit, decomp in enumerate(decomposition[1::2]): - euler_angles = xzx_decomposer.angles_and_phase(decomp) - euler_q1[iqubit, [1, 2, 0]] = euler_angles[:3] - global_phase += euler_angles[3] - - qc = QuantumCircuit(2) - qc.global_phase = target_decomposed.global_phase - qc.global_phase -= best_nbasis * self.basis.global_phase - qc.global_phase += global_phase - - x12 = euler_q0[1][2] + euler_q0[2][0] - x12_isNonZero = not math.isclose(x12, 0, abs_tol=atol) - x12_isOddMult = None - x12_isPiMult = math.isclose(math.sin(x12), 0, abs_tol=atol) - if x12_isPiMult: - x12_isOddMult = math.isclose(math.cos(x12), -1, abs_tol=atol) - x12_phase = math.pi * math.cos(x12) - x02_add = x12 - euler_q0[1][0] - x12_isHalfPi = math.isclose(x12, math.pi / 2, abs_tol=atol) - - # TODO: make this more effecient to avoid double decomposition - circ = QuantumCircuit(1) - circ.rz(euler_q0[0][0], 0) - circ.rx(euler_q0[0][1], 0) - if x12_isNonZero and x12_isPiMult: - circ.rz(euler_q0[0][2] - x02_add, 0) - else: - circ.rz(euler_q0[0][2] + euler_q0[1][0], 0) - circ.h(0) - qceuler = self._decomposer1q(Operator(circ).data) - qc.compose(qceuler, [0], inplace=True) - - circ = QuantumCircuit(1) - circ.rx(euler_q1[0][0], 0) - circ.rz(euler_q1[0][1], 0) - circ.rx(euler_q1[0][2] + euler_q1[1][0], 0) - circ.h(0) - qceuler = self._decomposer1q(Operator(circ).data) - qc.compose(qceuler, [1], inplace=True) - - qc.cx(1, 0) - - if x12_isPiMult: - # even or odd multiple - if x12_isNonZero: - qc.global_phase += x12_phase - if x12_isNonZero and x12_isOddMult: - qc.rz(-euler_q0[1][1], 0) - else: - qc.rz(euler_q0[1][1], 0) - qc.global_phase += math.pi - if x12_isHalfPi: - qc.sx(0) - qc.global_phase -= math.pi / 4 - elif x12_isNonZero and not x12_isPiMult: - # this is non-optimal but doesn't seem to occur currently - if self.pulse_optimize is None: - qc.compose(self._decomposer1q(Operator(RXGate(x12)).data), [0], inplace=True) - else: - raise QiskitError("possible non-pulse-optimal decomposition encountered") - if math.isclose(euler_q1[1][1], math.pi / 2, abs_tol=atol): - qc.sx(1) - qc.global_phase -= math.pi / 4 - else: - # this is non-optimal but doesn't seem to occur currently - if self.pulse_optimize is None: - qc.compose( - self._decomposer1q(Operator(RXGate(euler_q1[1][1])).data), [1], inplace=True - ) - else: - raise QiskitError("possible non-pulse-optimal decomposition encountered") - qc.rz(euler_q1[1][2] + euler_q1[2][0], 1) - - qc.cx(1, 0) - - qc.rz(euler_q0[2][1], 0) - if math.isclose(euler_q1[2][1], math.pi / 2, abs_tol=atol): - qc.sx(1) - qc.global_phase -= math.pi / 4 - else: - # this is non-optimal but doesn't seem to occur currently - if self.pulse_optimize is None: - qc.compose( - self._decomposer1q(Operator(RXGate(euler_q1[2][1])).data), [1], inplace=True - ) - else: - raise QiskitError("possible non-pulse-optimal decomposition encountered") - - qc.cx(1, 0) - - circ = QuantumCircuit(1) - circ.h(0) - circ.rz(euler_q0[2][2] + euler_q0[3][0], 0) - circ.rx(euler_q0[3][1], 0) - circ.rz(euler_q0[3][2], 0) - qceuler = self._decomposer1q(Operator(circ).data) - qc.compose(qceuler, [0], inplace=True) - - circ = QuantumCircuit(1) - circ.h(0) - circ.rx(euler_q1[2][2] + euler_q1[3][0], 0) - circ.rz(euler_q1[3][1], 0) - circ.rx(euler_q1[3][2], 0) - qceuler = self._decomposer1q(Operator(circ).data) - qc.compose(qceuler, [1], inplace=True) - - # TODO: fix the sign problem to avoid correction here - if cmath.isclose( - target_decomposed.unitary_matrix[0, 0], -(Operator(qc).data[0, 0]), abs_tol=atol - ): - qc.global_phase += math.pi - return qc + circ = QuantumCircuit(q, global_phase=sequence.global_phase) + for name, params, qubits in sequence: + try: + getattr(circ, name)(*params, *qubits) + except AttributeError as exc: + if name == "USER_GATE": + circ.append(self.gate, qubits) + elif name == "u3": + gate = U3Gate(*params) + circ.append(gate, qubits) + elif name == "u2": + gate = U2Gate(*params) + circ.append(gate, qubits) + elif name == "u1": + gate = U1Gate(*params) + circ.append(gate, qubits) + else: + raise QiskitError(f"Unknown gate {name}") from exc - def num_basis_gates(self, unitary): - """Computes the number of basis gates needed in - a decomposition of input unitary - """ - return two_qubit_decompose._num_basis_gates( - self.basis.b, self.basis_fidelity, np.asarray(unitary, dtype=complex) - ) + return circ class TwoQubitDecomposeUpToDiagonal: diff --git a/qiskit/transpiler/passes/optimization/consolidate_blocks.py b/qiskit/transpiler/passes/optimization/consolidate_blocks.py index 41407ac2af28..8dca049aa9ac 100644 --- a/qiskit/transpiler/passes/optimization/consolidate_blocks.py +++ b/qiskit/transpiler/passes/optimization/consolidate_blocks.py @@ -131,7 +131,7 @@ def run(self, dag): if ( # pylint: disable=too-many-boolean-expressions self.force_consolidate or unitary.num_qubits > 2 - or self.decomposer.num_basis_gates(unitary) < basis_count + or self.decomposer.num_basis_gates(matrix) < basis_count or len(block) > max_2q_depth or ((self.basis_gates is not None) and outside_basis) or ((self.target is not None) and outside_basis) diff --git a/qiskit/transpiler/passes/synthesis/unitary_synthesis.py b/qiskit/transpiler/passes/synthesis/unitary_synthesis.py index 5d9ba1e6cde7..e4bad121dacb 100644 --- a/qiskit/transpiler/passes/synthesis/unitary_synthesis.py +++ b/qiskit/transpiler/passes/synthesis/unitary_synthesis.py @@ -811,7 +811,7 @@ def is_controlled(gate): pi2_decomposer = TwoQubitBasisDecomposer( pi2_basis, euler_basis=basis_1q, - basis_fidelity=basis_2q_fidelity, + # basis_fidelity=basis_2q_fidelity, pulse_optimize=True, ) embodiments.update({pi / 2: XXEmbodiments[pi2_decomposer.gate.base_class]}) diff --git a/test/python/transpiler/test_unitary_synthesis.py b/test/python/transpiler/test_unitary_synthesis.py index 0f4423426a8a..ab31da639fb9 100644 --- a/test/python/transpiler/test_unitary_synthesis.py +++ b/test/python/transpiler/test_unitary_synthesis.py @@ -210,6 +210,8 @@ def test_swap_synthesis_to_directional_cx(self): pm_nat = PassManager([triv_layout_pass, unisynth_pass_nat]) qc_out_nat = pm_nat.run(qc) + print(qc) + print(qc_out) self.assertEqual(Operator(qc), Operator(qc_out)) self.assertEqual(Operator(qc), Operator(qc_out_nat))