Merge branch 'main' into goodbye-extensions

.github/CODEOWNERS

@@ -30,7 +30,7 @@ pulse/                @Qiskit/terra-core @eggerdj @wshanks
 synthesis/            @Qiskit/terra-core @alexanderivrii @ShellyGarion
 scheduler/            @Qiskit/terra-core @eggerdj @wshanks
 visualization/        @Qiskit/terra-core @nonhermitian
-primitives/           @Qiskit/terra-core @ikkoham @t-imamichi
+primitives/           @Qiskit/terra-core @Qiskit/qiskit-primitives
 # Override the release notes directories to have _no_ code owners, so any review
 # from somebody with write access is acceptable.
 /releasenotes/notes

crates/accelerate/Cargo.toml

@@ -21,12 +21,13 @@ rand = "0.8"
 rand_pcg = "0.3"
 rand_distr = "0.4.3"
 ahash = "0.8.6"
+num-traits = "0.2"
 num-complex = "0.4"
 num-bigint = "0.4"
 rustworkx-core = "0.13"
 
 [dependencies.smallvec]
-version = "1.11"
+version = "1.12"
 features = ["union"]
 
 [dependencies.pyo3]

crates/accelerate/src/quantum_circuit/circuit_data.rs

@@ -213,32 +213,33 @@ impl CircuitData {
         Ok((ty, args, None::<()>, self_.iter()?).into_py(py))
     }
 
-    /// Returns the current sequence of registered :class:`.Qubit`
-    /// instances as a list.
+    /// Returns the current sequence of registered :class:`.Qubit` instances as a list.
     ///
-    /// .. note::
+    /// .. warning::
     ///
-    ///     This list is not kept in sync with the container.
+    ///     Do not modify this list yourself.  It will invalidate the :class:`CircuitData` data
+    ///     structures.
     ///
     /// Returns:
     ///     list(:class:`.Qubit`): The current sequence of registered qubits.
     #[getter]
-    pub fn qubits(&self, py: Python<'_>) -> PyObject {
-        PyList::new(py, self.qubits.as_ref(py)).into_py(py)
+    pub fn qubits(&self, py: Python<'_>) -> Py<PyList> {
+        self.qubits.clone_ref(py)
     }
 
     /// Returns the current sequence of registered :class:`.Clbit`
     /// instances as a list.
     ///
-    /// .. note::
+    /// .. warning::
     ///
-    ///     This list is not kept in sync with the container.
+    ///     Do not modify this list yourself.  It will invalidate the :class:`CircuitData` data
+    ///     structures.
     ///
     /// Returns:
     ///     list(:class:`.Clbit`): The current sequence of registered clbits.
     #[getter]
-    pub fn clbits(&self, py: Python<'_>) -> PyObject {
-        PyList::new(py, self.clbits.as_ref(py)).into_py(py)
+    pub fn clbits(&self, py: Python<'_>) -> Py<PyList> {
+        self.clbits.clone_ref(py)
     }
 
     /// Registers a :class:`.Qubit` instance.
@@ -251,7 +252,13 @@ impl CircuitData {
     ///     ValueError: The specified ``bit`` is already present and flag ``strict``
     ///         was provided.
     #[pyo3(signature = (bit, *, strict=true))]
-    pub fn add_qubit(&mut self, py: Python<'_>, bit: &PyAny, strict: bool) -> PyResult<()> {
+    pub fn add_qubit(&mut self, py: Python, bit: &PyAny, strict: bool) -> PyResult<()> {
+        if self.qubits_native.len() != self.qubits.as_ref(bit.py()).len() {
+            return Err(PyRuntimeError::new_err(concat!(
+                "This circuit's 'qubits' list has become out of sync with the circuit data.",
+                " Did something modify it?"
+            )));
+        }
         let idx: BitType = self.qubits_native.len().try_into().map_err(|_| {
             PyRuntimeError::new_err(
                 "The number of qubits in the circuit has exceeded the maximum capacity",
@@ -263,7 +270,7 @@ impl CircuitData {
             .is_ok()
         {
             self.qubits_native.push(bit.into_py(py));
-            self.qubits = PyList::new(py, &self.qubits_native).into_py(py);
+            self.qubits.as_ref(py).append(bit)?;
         } else if strict {
             return Err(PyValueError::new_err(format!(
                 "Existing bit {:?} cannot be re-added in strict mode.",
@@ -283,7 +290,13 @@ impl CircuitData {
     ///     ValueError: The specified ``bit`` is already present and flag ``strict``
     ///         was provided.
     #[pyo3(signature = (bit, *, strict=true))]
-    pub fn add_clbit(&mut self, py: Python<'_>, bit: &PyAny, strict: bool) -> PyResult<()> {
+    pub fn add_clbit(&mut self, py: Python, bit: &PyAny, strict: bool) -> PyResult<()> {
+        if self.clbits_native.len() != self.clbits.as_ref(bit.py()).len() {
+            return Err(PyRuntimeError::new_err(concat!(
+                "This circuit's 'clbits' list has become out of sync with the circuit data.",
+                " Did something modify it?"
+            )));
+        }
         let idx: BitType = self.clbits_native.len().try_into().map_err(|_| {
             PyRuntimeError::new_err(
                 "The number of clbits in the circuit has exceeded the maximum capacity",
@@ -295,7 +308,7 @@ impl CircuitData {
             .is_ok()
         {
             self.clbits_native.push(bit.into_py(py));
-            self.clbits = PyList::new(py, &self.clbits_native).into_py(py);
+            self.clbits.as_ref(py).append(bit)?;
         } else if strict {
             return Err(PyValueError::new_err(format!(
                 "Existing bit {:?} cannot be re-added in strict mode.",

crates/accelerate/src/sparse_pauli_op.rs

@@ -10,13 +10,16 @@
 // copyright notice, and modified files need to carry a notice indicating
 // that they have been altered from the originals.
 
+use pyo3::exceptions::PyValueError;
 use pyo3::prelude::*;
 use pyo3::wrap_pyfunction;
 use pyo3::Python;
 
 use hashbrown::HashMap;
-use ndarray::{ArrayView1, Axis};
-use numpy::{IntoPyArray, PyReadonlyArray2};
+use ndarray::{s, Array1, Array2, ArrayView1, ArrayView2, Axis};
+use num_complex::Complex64;
+use num_traits::Zero;
+use numpy::{IntoPyArray, PyArray1, PyArray2, PyReadonlyArray2};
 
 /// Find the unique elements of an array.
 ///
@@ -60,8 +63,204 @@ pub fn unordered_unique(py: Python, array: PyReadonlyArray2<u16>) -> (PyObject,
     )
 }
 
+#[derive(Clone, Copy)]
+enum Pauli {
+    I,
+    X,
+    Y,
+    Z,
+}
+
+/// A complete ZX-convention representation of a Pauli decomposition.  This is all the components
+/// necessary to construct a Qiskit-space :class:`.SparsePauliOp`, where :attr:`phases` is in the
+/// ZX convention.
+#[pyclass(module = "qiskit._accelerate.sparse_pauli_op")]
+pub struct ZXPaulis {
+    #[pyo3(get)]
+    pub z: Py<PyArray2<bool>>,
+    #[pyo3(get)]
+    pub x: Py<PyArray2<bool>>,
+    #[pyo3(get)]
+    pub phases: Py<PyArray1<u8>>,
+    #[pyo3(get)]
+    pub coeffs: Py<PyArray1<Complex64>>,
+}
+
+/// Decompose a dense complex operator into the symplectic Pauli representation in the
+/// ZX-convention.
+///
+/// This is an implementation of the "tensorized Pauli decomposition" presented in
+/// `Hantzko, Binkowski and Gupta (2023) <https://arxiv.org/abs/2310.13421>`__.
+#[pyfunction]
+pub fn decompose_dense(
+    py: Python,
+    operator: PyReadonlyArray2<Complex64>,
+    tolerance: f64,
+) -> PyResult<ZXPaulis> {
+    let num_qubits = operator.shape()[0].ilog2() as usize;
+    let size = 1 << num_qubits;
+    if operator.shape() != [size, size] {
+        return Err(PyValueError::new_err(format!(
+            "input with shape {:?} cannot be interpreted as a multiqubit operator",
+            operator.shape()
+        )));
+    }
+    let mut paulis = vec![];
+    let mut coeffs = vec![];
+    if num_qubits > 0 {
+        decompose_dense_inner(
+            Complex64::new(1.0, 0.0),
+            num_qubits,
+            &[],
+            operator.as_array(),
+            &mut paulis,
+            &mut coeffs,
+            tolerance * tolerance,
+        );
+    }
+    if coeffs.is_empty() {
+        Ok(ZXPaulis {
+            z: PyArray2::zeros(py, [0, num_qubits], false).into(),
+            x: PyArray2::zeros(py, [0, num_qubits], false).into(),
+            phases: PyArray1::zeros(py, [0], false).into(),
+            coeffs: PyArray1::zeros(py, [0], false).into(),
+        })
+    } else {
+        // Constructing several arrays of different shapes at once is rather awkward in iterator
+        // logic, so we just loop manually.
+        let mut z = Array2::<bool>::uninit([paulis.len(), num_qubits]);
+        let mut x = Array2::<bool>::uninit([paulis.len(), num_qubits]);
+        let mut phases = Array1::<u8>::uninit(paulis.len());
+        for (i, paulis) in paulis.drain(..).enumerate() {
+            let mut phase = 0u8;
+            for (j, pauli) in paulis.into_iter().rev().enumerate() {
+                match pauli {
+                    Pauli::I => {
+                        z[[i, j]].write(false);
+                        x[[i, j]].write(false);
+                    }
+                    Pauli::X => {
+                        z[[i, j]].write(false);
+                        x[[i, j]].write(true);
+                    }
+                    Pauli::Y => {
+                        z[[i, j]].write(true);
+                        x[[i, j]].write(true);
+                        phase = phase.wrapping_add(1);
+                    }
+                    Pauli::Z => {
+                        z[[i, j]].write(true);
+                        x[[i, j]].write(false);
+                    }
+                }
+            }
+            phases[i].write(phase % 4);
+        }
+        // These are safe because the above loops write into every element.  It's guaranteed that
+        // each of the elements of the `paulis` vec will have `num_qubits` because they're all
+        // reading from the same base array.
+        let z = unsafe { z.assume_init() };
+        let x = unsafe { x.assume_init() };
+        let phases = unsafe { phases.assume_init() };
+        Ok(ZXPaulis {
+            z: z.into_pyarray(py).into(),
+            x: x.into_pyarray(py).into(),
+            phases: phases.into_pyarray(py).into(),
+            coeffs: PyArray1::from_vec(py, coeffs).into(),
+        })
+    }
+}
+
+/// Recurse worker routine of `decompose_dense`.  Should be called with at least one qubit.
+fn decompose_dense_inner(
+    factor: Complex64,
+    num_qubits: usize,
+    paulis: &[Pauli],
+    block: ArrayView2<Complex64>,
+    out_paulis: &mut Vec<Vec<Pauli>>,
+    out_coeffs: &mut Vec<Complex64>,
+    square_tolerance: f64,
+) {
+    if num_qubits == 0 {
+        // It would be safe to `return` here, but if it's unreachable then LLVM is allowed to
+        // optimise out this branch entirely in release mode, which is good for a ~2% speedup.
+        unreachable!("should not call this with an empty operator")
+    }
+    // Base recursion case.
+    if num_qubits == 1 {
+        let mut push_if_nonzero = |extra: Pauli, value: Complex64| {
+            if value.norm_sqr() <= square_tolerance {
+                return;
+            }
+            let paulis = {
+                let mut vec = Vec::with_capacity(paulis.len() + 1);
+                vec.extend_from_slice(paulis);
+                vec.push(extra);
+                vec
+            };
+            out_paulis.push(paulis);
+            out_coeffs.push(value);
+        };
+        push_if_nonzero(Pauli::I, 0.5 * factor * (block[[0, 0]] + block[[1, 1]]));
+        push_if_nonzero(Pauli::X, 0.5 * factor * (block[[0, 1]] + block[[1, 0]]));
+        push_if_nonzero(
+            Pauli::Y,
+            0.5 * Complex64::i() * factor * (block[[0, 1]] - block[[1, 0]]),
+        );
+        push_if_nonzero(Pauli::Z, 0.5 * factor * (block[[0, 0]] - block[[1, 1]]));
+        return;
+    }
+    let mut recurse_if_nonzero = |extra: Pauli, factor: Complex64, values: Array2<Complex64>| {
+        let mut is_zero = true;
+        for value in values.iter() {
+            if !value.is_zero() {
+                is_zero = false;
+                break;
+            }
+        }
+        if is_zero {
+            return;
+        }
+        let mut new_paulis = Vec::with_capacity(paulis.len() + 1);
+        new_paulis.extend_from_slice(paulis);
+        new_paulis.push(extra);
+        decompose_dense_inner(
+            factor,
+            num_qubits - 1,
+            &new_paulis,
+            values.view(),
+            out_paulis,
+            out_coeffs,
+            square_tolerance,
+        );
+    };
+    let mid = 1usize << (num_qubits - 1);
+    recurse_if_nonzero(
+        Pauli::I,
+        0.5 * factor,
+        &block.slice(s![..mid, ..mid]) + &block.slice(s![mid.., mid..]),
+    );
+    recurse_if_nonzero(
+        Pauli::X,
+        0.5 * factor,
+        &block.slice(s![..mid, mid..]) + &block.slice(s![mid.., ..mid]),
+    );
+    recurse_if_nonzero(
+        Pauli::Y,
+        0.5 * Complex64::i() * factor,
+        &block.slice(s![..mid, mid..]) - &block.slice(s![mid.., ..mid]),
+    );
+    recurse_if_nonzero(
+        Pauli::Z,
+        0.5 * factor,
+        &block.slice(s![..mid, ..mid]) - &block.slice(s![mid.., mid..]),
+    );
+}
+
 #[pymodule]
 pub fn sparse_pauli_op(_py: Python, m: &PyModule) -> PyResult<()> {
     m.add_wrapped(wrap_pyfunction!(unordered_unique))?;
+    m.add_wrapped(wrap_pyfunction!(decompose_dense))?;
+    m.add_class::<ZXPaulis>()?;
     Ok(())
 }

docs/apidoc/index.rst

@@ -33,8 +33,6 @@ API Reference
    qpy
    quantum_info
    result
-   tools
-   tools_jupyter
    transpiler
    transpiler_passes
    transpiler_preset

docs/conf.py

@@ -42,6 +42,8 @@
     "sphinx.ext.autosummary",
     "sphinx.ext.intersphinx",
     "sphinx.ext.doctest",
+    # This is used by qiskit/documentation to generate links to github.com.
+    "sphinx.ext.viewcode",
     "matplotlib.sphinxext.plot_directive",
     "reno.sphinxext",
     "sphinxcontrib.katex",

pyproject.toml

@@ -45,7 +45,6 @@ qasm3-import = [
 ]
 visualization = [
     "matplotlib >= 3.3",
-    "ipywidgets >= 7.3.0",
     "pydot",
     "Pillow >= 4.2.1",
     "pylatexenc >= 1.4",