Clean two_qubit_decompose code

crates/accelerate/src/two_qubit_decompose.rs

@@ -32,7 +32,7 @@ use ndarray::linalg::kron;
 use ndarray::prelude::*;
 use ndarray::Zip;
 use numpy::{IntoPyArray, ToPyArray};
-use numpy::{PyReadonlyArray1, PyReadonlyArray2};
+use numpy::{PyArray2, PyArrayLike2, PyReadonlyArray1, PyReadonlyArray2};
 
 use pyo3::exceptions::PyValueError;
 use pyo3::prelude::*;
@@ -154,6 +154,15 @@ impl TraceToFidelity for c64 {
     }
 }
 
+#[pyfunction]
+#[pyo3(name = "trace_to_fid")]
+/// 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 py_trace_to_fid(trace: Complex64) -> PyResult<f64> {
+    let fid = trace.trace_to_fid();
+    Ok(fid)
+}
+
 fn decompose_two_qubit_product_gate(
     special_unitary: ArrayView2<Complex64>,
 ) -> PyResult<(Array2<Complex64>, Array2<Complex64>, f64)> {
@@ -182,9 +191,31 @@ fn decompose_two_qubit_product_gate(
     }
     l.mapv_inplace(|x| x / det_l.sqrt());
     let phase = det_l.arg() / 2.;
+
     Ok((l, r, phase))
 }
 
+#[pyfunction]
+#[pyo3(name = "decompose_two_qubit_product_gate")]
+/// Decompose :math:`U = U_l \otimes U_r` where :math:`U \in SU(4)`,
+/// and :math:`U_l,~U_r \in SU(2)`.
+/// Args:
+///     special_unitary_matrix: special unitary matrix to decompose
+/// Raises:
+///     QiskitError: if decomposition isn't possible.
+fn py_decompose_two_qubit_product_gate(
+    py: Python,
+    special_unitary: PyArrayLike2<Complex64>,
+) -> PyResult<(PyObject, PyObject, f64)> {
+    let view = special_unitary.as_array();
+    let (l, r, phase) = decompose_two_qubit_product_gate(view)?;
+    Ok((
+        l.into_pyarray_bound(py).unbind().into(),
+        r.into_pyarray_bound(py).unbind().into(),
+        phase,
+    ))
+}
+
 fn __weyl_coordinates(unitary: MatRef<c64>) -> [f64; 3] {
     let uscaled = scale(C1 / unitary.determinant().powf(0.25)) * unitary;
     let uup = transform_from_magic_basis(uscaled);
@@ -301,6 +332,43 @@ fn rz_matrix(theta: f64) -> Array2<Complex64> {
     array![[(-ilam2).exp(), C_ZERO], [C_ZERO, ilam2.exp()]]
 }
 
+/// Generates the array :math:`e^{(i a XX + i b YY + i c ZZ)}`
+fn ud(a: f64, b: f64, c: f64) -> Array2<Complex64> {
+    array![
+        [
+            (IM * c).exp() * (a - b).cos(),
+            C_ZERO,
+            C_ZERO,
+            IM * (IM * c).exp() * (a - b).sin()
+        ],
+        [
+            C_ZERO,
+            (M_IM * c).exp() * (a + b).cos(),
+            IM * (M_IM * c).exp() * (a + b).sin(),
+            C_ZERO
+        ],
+        [
+            C_ZERO,
+            IM * (M_IM * c).exp() * (a + b).sin(),
+            (M_IM * c).exp() * (a + b).cos(),
+            C_ZERO
+        ],
+        [
+            IM * (IM * c).exp() * (a - b).sin(),
+            C_ZERO,
+            C_ZERO,
+            (IM * c).exp() * (a - b).cos()
+        ]
+    ]
+}
+
+#[pyfunction]
+#[pyo3(name = "Ud")]
+fn py_ud(py: Python, a: f64, b: f64, c: f64) -> PyResult<Py<PyArray2<Complex64>>> {
+    let ud_mat = ud(a, b, c);
+    Ok(ud_mat.into_pyarray_bound(py).unbind())
+}
+
 fn compute_unitary(sequence: &TwoQubitSequenceVec, global_phase: f64) -> Array2<Complex64> {
     let identity = aview2(&ONE_QUBIT_IDENTITY);
     let phase = c64(0., global_phase).exp();
@@ -2278,9 +2346,12 @@ pub fn local_equivalence(weyl: PyReadonlyArray1<f64>) -> PyResult<[f64; 3]> {
 
 pub fn two_qubit_decompose(m: &Bound<PyModule>) -> PyResult<()> {
     m.add_wrapped(wrap_pyfunction!(_num_basis_gates))?;
+    m.add_wrapped(wrap_pyfunction!(py_decompose_two_qubit_product_gate))?;
     m.add_wrapped(wrap_pyfunction!(two_qubit_decompose_up_to_diagonal))?;
     m.add_wrapped(wrap_pyfunction!(two_qubit_local_invariants))?;
     m.add_wrapped(wrap_pyfunction!(local_equivalence))?;
+    m.add_wrapped(wrap_pyfunction!(py_trace_to_fid))?;
+    m.add_wrapped(wrap_pyfunction!(py_ud))?;
     m.add_class::<TwoQubitGateSequence>()?;
     m.add_class::<TwoQubitWeylDecomposition>()?;
     m.add_class::<Specialization>()?;

qiskit/synthesis/two_qubit/two_qubit_decompose.py

@@ -24,8 +24,6 @@
 arXiv:1811.12926 [quant-ph] (2018).
 """
 from __future__ import annotations
-import cmath
-import math
 import io
 import base64
 import warnings
@@ -91,41 +89,18 @@ def decompose_two_qubit_product_gate(special_unitary_matrix: np.ndarray):
         QiskitError: if decomposition isn't possible.
     """
     special_unitary_matrix = np.asarray(special_unitary_matrix, dtype=complex)
-    # extract the right component
-    R = special_unitary_matrix[:2, :2].copy()
-    detR = R[0, 0] * R[1, 1] - R[0, 1] * R[1, 0]
-    if abs(detR) < 0.1:
-        R = special_unitary_matrix[2:, :2].copy()
-        detR = R[0, 0] * R[1, 1] - R[0, 1] * R[1, 0]
-    if abs(detR) < 0.1:
-        raise QiskitError("decompose_two_qubit_product_gate: unable to decompose: detR < 0.1")
-    R /= np.sqrt(detR)
-
-    # extract the left component
-    temp = np.kron(np.eye(2), R.T.conj())
-    temp = special_unitary_matrix.dot(temp)
-    L = temp[::2, ::2]
-    detL = L[0, 0] * L[1, 1] - L[0, 1] * L[1, 0]
-    if abs(detL) < 0.9:
-        raise QiskitError("decompose_two_qubit_product_gate: unable to decompose: detL < 0.9")
-    L /= np.sqrt(detL)
-    phase = cmath.phase(detL) / 2
+    (L, R, phase) = two_qubit_decompose.decompose_two_qubit_product_gate(special_unitary_matrix)
 
     temp = np.kron(L, R)
     deviation = abs(abs(temp.conj().T.dot(special_unitary_matrix).trace()) - 4)
+
     if deviation > 1.0e-13:
         raise QiskitError(
             "decompose_two_qubit_product_gate: decomposition failed: "
             f"deviation too large: {deviation}"
         )
 
-    return L, R, phase
-
-
-_ipx = np.array([[0, 1j], [1j, 0]], dtype=complex)
-_ipy = np.array([[0, 1], [-1, 0]], dtype=complex)
-_ipz = np.array([[1j, 0], [0, -1j]], dtype=complex)
-_id = np.array([[1, 0], [0, 1]], dtype=complex)
+    return (L, R, phase)
 
 
 class TwoQubitWeylDecomposition:
@@ -239,7 +214,7 @@ 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)
+        return two_qubit_decompose.trace_to_fid(trace)
 
     def __repr__(self):
         """Represent with enough precision to allow copy-paste debugging of all corner cases"""
@@ -460,40 +435,6 @@ def _weyl_gate(self, circ: QuantumCircuit, atol=1.0e-13):
         return circ
 
 
-def Ud(a, b, c):
-    r"""Generates the array :math:`e^{(i a XX + i b YY + i c ZZ)}`"""
-    return np.array(
-        [
-            [cmath.exp(1j * c) * math.cos(a - b), 0, 0, 1j * cmath.exp(1j * c) * math.sin(a - b)],
-            [0, cmath.exp(-1j * c) * math.cos(a + b), 1j * cmath.exp(-1j * c) * math.sin(a + b), 0],
-            [0, 1j * cmath.exp(-1j * c) * math.sin(a + b), cmath.exp(-1j * c) * math.cos(a + b), 0],
-            [1j * cmath.exp(1j * c) * math.sin(a - b), 0, 0, cmath.exp(1j * c) * math.cos(a - b)],
-        ],
-        dtype=complex,
-    )
-
-
-def trace_to_fid(trace):
-    r"""Average gate fidelity is
-
-    .. math::
-
-        \bar{F} = \frac{d + |\mathrm{Tr} (U_\text{target} \cdot U^{\dag})|^2}{d(d+1)}
-
-    M. Horodecki, P. Horodecki and R. Horodecki, PRA 60, 1888 (1999)"""
-    return (4 + abs(trace) ** 2) / 20
-
-
-def rz_array(theta):
-    """Return numpy array for Rz(theta).
-
-    Rz(theta) = diag(exp(-i*theta/2),exp(i*theta/2))
-    """
-    return np.array(
-        [[cmath.exp(-1j * theta / 2.0), 0], [0, cmath.exp(1j * theta / 2.0)]], dtype=complex
-    )
-
-
 class TwoQubitBasisDecomposer:
     """A class for decomposing 2-qubit unitaries into minimal number of uses of a 2-qubit
     basis gate.

test/python/synthesis/test_synthesis.py

@@ -59,11 +59,11 @@
     two_qubit_cnot_decompose,
     TwoQubitBasisDecomposer,
     TwoQubitControlledUDecomposer,
-    Ud,
     decompose_two_qubit_product_gate,
 )
 from qiskit._accelerate.two_qubit_decompose import two_qubit_decompose_up_to_diagonal
 from qiskit._accelerate.two_qubit_decompose import Specialization
+from qiskit._accelerate.two_qubit_decompose import Ud
 from qiskit.synthesis.unitary import qsd
 from test import combine  # pylint: disable=wrong-import-order
 from test import QiskitTestCase  # pylint: disable=wrong-import-order

test/randomized/test_synthesis.py

@@ -23,8 +23,8 @@
 from qiskit.synthesis.two_qubit.two_qubit_decompose import (
     two_qubit_cnot_decompose,
     TwoQubitBasisDecomposer,
-    Ud,
 )
+from qiskit._accelerate.two_qubit_decompose import Ud
 
 
 class TestSynthesis(CheckDecompositions):