diff --git a/Cargo.lock b/Cargo.lock index 3de6bc619612..4d6f7284abc5 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -1227,6 +1227,7 @@ version = "1.3.0" dependencies = [ "ahash 0.8.11", "approx 0.5.1", + "bytemuck", "faer", "faer-ext", "hashbrown 0.14.5", diff --git a/crates/accelerate/Cargo.toml b/crates/accelerate/Cargo.toml index df8c2436a5c9..85f3b7b3bea4 100644 --- a/crates/accelerate/Cargo.toml +++ b/crates/accelerate/Cargo.toml @@ -26,6 +26,7 @@ qiskit-circuit.workspace = true thiserror.workspace = true ndarray_einsum_beta = "0.7" once_cell = "1.20.0" +bytemuck.workspace = true [dependencies.smallvec] workspace = true diff --git a/crates/accelerate/src/lib.rs b/crates/accelerate/src/lib.rs index 7c347c0a3d75..44a016f2c608 100644 --- a/crates/accelerate/src/lib.rs +++ b/crates/accelerate/src/lib.rs @@ -39,6 +39,7 @@ pub mod remove_diagonal_gates_before_measure; pub mod results; pub mod sabre; pub mod sampled_exp_val; +pub mod sparse_observable; pub mod sparse_pauli_op; pub mod split_2q_unitaries; pub mod star_prerouting; diff --git a/crates/accelerate/src/sparse_observable.rs b/crates/accelerate/src/sparse_observable.rs new file mode 100644 index 000000000000..5bc77e825510 --- /dev/null +++ b/crates/accelerate/src/sparse_observable.rs @@ -0,0 +1,1674 @@ +// This code is part of Qiskit. +// +// (C) Copyright IBM 2024 +// +// This code is licensed under the Apache License, Version 2.0. You may +// obtain a copy of this license in the LICENSE.txt file in the root directory +// of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +// +// Any modifications or derivative works of this code must retain this +// copyright notice, and modified files need to carry a notice indicating +// that they have been altered from the originals. + +use std::collections::btree_map; + +use num_complex::Complex64; +use thiserror::Error; + +use numpy::{ + PyArray1, PyArrayDescr, PyArrayDescrMethods, PyArrayLike1, PyReadonlyArray1, PyReadonlyArray2, + PyUntypedArrayMethods, +}; + +use pyo3::exceptions::{PyTypeError, PyValueError}; +use pyo3::intern; +use pyo3::prelude::*; +use pyo3::sync::GILOnceCell; +use pyo3::types::{IntoPyDict, PyList, PyType}; + +use qiskit_circuit::imports::{ImportOnceCell, NUMPY_COPY_ONLY_IF_NEEDED}; +use qiskit_circuit::slice::{PySequenceIndex, SequenceIndex}; + +static PAULI_TYPE: ImportOnceCell = ImportOnceCell::new("qiskit.quantum_info", "Pauli"); +static SPARSE_PAULI_OP_TYPE: ImportOnceCell = + ImportOnceCell::new("qiskit.quantum_info", "SparsePauliOp"); + +/// Named handle to the alphabet of single-qubit terms. +/// +/// This is just the Rust-space representation. We make a separate Python-space `enum.IntEnum` to +/// represent the same information, since we enforce strongly typed interactions in Rust, including +/// not allowing the stored values to be outside the valid `BitTerm`s, but doing so in Python would +/// make it very difficult to use the class efficiently with Numpy array views. We attach this +/// sister class of `BitTerm` to `SparseObservable` as a scoped class. +/// +/// # Representation +/// +/// The `u8` representation and the exact numerical values of these are part of the public API. The +/// low two bits are the symplectic Pauli representation of the required measurement basis with Z in +/// the Lsb0 and X in the Lsb1 (e.g. X and its eigenstate projectors all have their two low bits as +/// `0b10`). The high two bits are `00` for the operator, `10` for the projector to the positive +/// eigenstate, and `01` for the projector to the negative eigenstate. +/// +/// The `0b00_00` representation thus ends up being the natural representation of the `I` operator, +/// but this is never stored, and is not named in the enumeration. +/// +/// This operator does not store phase terms of $-i$. `BitTerm::Y` has `(1, 1)` as its `(z, x)` +/// representation, and represents exactly the Pauli Y operator; any additional phase is stored only +/// in a corresponding coefficient. +/// +/// # Dev notes +/// +/// This type is required to be `u8`, but it's a subtype of `u8` because not all `u8` are valid +/// `BitTerm`s. For interop with Python space, we accept Numpy arrays of `u8` to represent this, +/// which we transmute into slices of `BitTerm`, after checking that all the values are correct (or +/// skipping the check if Python space promises that it upheld the checks). +/// +/// We deliberately _don't_ impl `numpy::Element` for `BitTerm` (which would let us accept and +/// return `PyArray1` at Python-space boundaries) so that it's clear when we're doing +/// the transmute, and we have to be explicit about the safety of that. +#[repr(u8)] +#[derive(Clone, Copy, Debug, Hash, PartialEq, Eq)] +pub enum BitTerm { + /// Pauli X operator. + X = 0b00_10, + /// Projector to the positive eigenstate of Pauli X. + Plus = 0b10_10, + /// Projector to the negative eigenstate of Pauli X. + Minus = 0b01_10, + /// Pauli Y operator. + Y = 0b00_11, + /// Projector to the positive eigenstate of Pauli Y. + Right = 0b10_11, + /// Projector to the negative eigenstate of Pauli Y. + Left = 0b01_11, + /// Pauli Z operator. + Z = 0b00_01, + /// Projector to the positive eigenstate of Pauli Z. + Zero = 0b10_01, + /// Projector to the negative eigenstate of Pauli Z. + One = 0b01_01, +} +impl<'py> FromPyObject<'py> for BitTerm { + fn extract_bound(ob: &Bound<'py, PyAny>) -> PyResult { + ob.extract::()?.try_into().map_err(PyErr::from) + } +} +unsafe impl ::bytemuck::CheckedBitPattern for BitTerm { + type Bits = u8; + + #[inline(always)] + fn is_valid_bit_pattern(bits: &Self::Bits) -> bool { + *bits <= 0b11_11 && (*bits & 0b11_00) < 0b11_00 && (*bits & 0b00_11) != 0 + } +} +unsafe impl ::bytemuck::NoUninit for BitTerm {} + +impl BitTerm { + /// Get the label of this `BitTerm` used in Python-space applications. This is a single-letter + /// string. + #[inline] + fn py_label(&self) -> &'static str { + match self { + Self::X => "X", + Self::Plus => "+", + Self::Minus => "-", + Self::Y => "Y", + Self::Right => "r", + Self::Left => "l", + Self::Z => "Z", + Self::Zero => "0", + Self::One => "1", + } + } + + /// Get the name of this `BitTerm`, which is how Python space refers to the integer constant. + #[inline] + fn py_name(&self) -> &'static str { + match self { + Self::X => "X", + Self::Plus => "PLUS", + Self::Minus => "MINUS", + Self::Y => "Y", + Self::Right => "RIGHT", + Self::Left => "LEFT", + Self::Z => "Z", + Self::Zero => "ZERO", + Self::One => "ONE", + } + } + + /// Attempt to convert a `u8` into `BitTerm`. + /// + /// Unlike the implementation of `TryFrom`, this allows `b'I'` as an alphabet letter, + /// returning `Ok(None)` for it. All other letters outside the alphabet return the complete + /// error condition. + #[inline] + fn try_from_u8(value: u8) -> Result, BitTermFromU8Error> { + match value { + b'+' => Ok(Some(BitTerm::Plus)), + b'-' => Ok(Some(BitTerm::Minus)), + b'0' => Ok(Some(BitTerm::Zero)), + b'1' => Ok(Some(BitTerm::One)), + b'I' => Ok(None), + b'X' => Ok(Some(BitTerm::X)), + b'Y' => Ok(Some(BitTerm::Y)), + b'Z' => Ok(Some(BitTerm::Z)), + b'l' => Ok(Some(BitTerm::Left)), + b'r' => Ok(Some(BitTerm::Right)), + _ => Err(BitTermFromU8Error(value)), + } + } +} + +static BIT_TERM_PY_ENUM: GILOnceCell> = GILOnceCell::new(); +static BIT_TERM_INTO_PY: GILOnceCell<[Option>; 16]> = GILOnceCell::new(); + +/// Construct the Python-space `IntEnum` that represents the same values as the Rust-spce `BitTerm`. +/// +/// We don't make `BitTerm` a direct `pyclass` because we want the behaviour of `IntEnum`, which +/// specifically also makes its variants subclasses of the Python `int` type; we use a type-safe +/// enum in Rust, but from Python space we expect people to (carefully) deal with the raw ints in +/// Numpy arrays for efficiency. +/// +/// The resulting class is attached to `SparseObservable` as a class attribute, and its +/// `__qualname__` is set to reflect this. +fn make_py_bit_term(py: Python) -> PyResult> { + let terms = [ + BitTerm::X, + BitTerm::Plus, + BitTerm::Minus, + BitTerm::Y, + BitTerm::Right, + BitTerm::Left, + BitTerm::Z, + BitTerm::Zero, + BitTerm::One, + ] + .into_iter() + .flat_map(|term| { + let mut out = vec![(term.py_name(), term as u8)]; + if term.py_name() != term.py_label() { + // Also ensure that the labels are created as aliases. These can't be (easily) accessed + // by attribute-getter (dot) syntax, but will work with the item-getter (square-bracket) + // syntax, or programmatically with `getattr`. + out.push((term.py_label(), term as u8)); + } + out + }) + .collect::>(); + let obj = py.import_bound("enum")?.getattr("IntEnum")?.call( + ("BitTerm", terms), + Some( + &[ + ("module", "qiskit.quantum_info"), + ("qualname", "SparseObservable.BitTerm"), + ] + .into_py_dict_bound(py), + ), + )?; + Ok(obj.downcast_into::()?.unbind()) +} + +// Return the relevant value from the Python-space sister enumeration. These are Python-space +// singletons and subclasses of Python `int`. We only use this for interaction with "high level" +// Python space; the efficient Numpy-like array paths use `u8` directly so Numpy can act on it +// efficiently. +impl IntoPy> for BitTerm { + fn into_py(self, py: Python) -> Py { + let terms = BIT_TERM_INTO_PY.get_or_init(py, || { + let py_enum = BIT_TERM_PY_ENUM + .get_or_try_init(py, || make_py_bit_term(py)) + .expect("creating a simple Python enum class should be infallible") + .bind(py); + ::std::array::from_fn(|val| { + ::bytemuck::checked::try_cast(val as u8) + .ok() + .map(|term: BitTerm| { + py_enum + .getattr(term.py_name()) + .expect("the created `BitTerm` enum should have matching attribute names to the terms") + .unbind() + }) + }) + }); + terms[self as usize] + .as_ref() + .expect("the lookup table initialiser populated a 'Some' in all valid locations") + .clone_ref(py) + } +} +impl ToPyObject for BitTerm { + fn to_object(&self, py: Python) -> Py { + self.into_py(py) + } +} + +/// The error type for a failed conversion into `BitTerm`. +#[derive(Error, Debug)] +#[error("{0} is not a valid letter of the single-qubit alphabet")] +pub struct BitTermFromU8Error(u8); +impl From for PyErr { + fn from(value: BitTermFromU8Error) -> PyErr { + PyValueError::new_err(value.to_string()) + } +} + +// `BitTerm` allows safe `as` casting into `u8`. This is the reverse, which is fallible, because +// `BitTerm` is a value-wise subtype of `u8`. +impl ::std::convert::TryFrom for BitTerm { + type Error = BitTermFromU8Error; + + fn try_from(value: u8) -> Result { + ::bytemuck::checked::try_cast(value).map_err(|_| BitTermFromU8Error(value)) + } +} + +/// Error cases stemming from data coherence at the point of entry into `SparseObservable` from raw +/// arrays. +/// +/// These are generally associated with the Python-space `ValueError` because all of the +/// `TypeError`-related ones are statically forbidden (within Rust) by the language, and conversion +/// failures on entry to Rust from Python space will automatically raise `TypeError`. +#[derive(Error, Debug)] +pub enum CoherenceError { + #[error("`boundaries` ({boundaries}) must be one element longer than `coeffs` ({coeffs})")] + MismatchedTermCount { coeffs: usize, boundaries: usize }, + #[error("`bit_terms` ({bit_terms}) and `indices` ({indices}) must be the same length")] + MismatchedItemCount { bit_terms: usize, indices: usize }, + #[error("the first item of `boundaries` ({0}) must be 0")] + BadInitialBoundary(usize), + #[error("the last item of `boundaries` ({last}) must match the length of `bit_terms` and `indices` ({items})")] + BadFinalBoundary { last: usize, items: usize }, + #[error("all qubit indices must be less than the number of qubits")] + BitIndexTooHigh, + #[error("the values in `boundaries` include backwards slices")] + DecreasingBoundaries, + #[error("the values in `indices` are not term-wise increasing")] + UnsortedIndices, +} +impl From for PyErr { + fn from(value: CoherenceError) -> PyErr { + PyValueError::new_err(value.to_string()) + } +} + +/// An error related to processing of a string label (both dense and sparse). +#[derive(Error, Debug)] +pub enum LabelError { + #[error("label with length {label} cannot be added to a {num_qubits}-qubit operator")] + WrongLengthDense { num_qubits: u32, label: usize }, + #[error("label with length {label} does not match indices of length {indices}")] + WrongLengthIndices { label: usize, indices: usize }, + #[error("index {index} is out of range for a {num_qubits}-qubit operator")] + BadIndex { index: u32, num_qubits: u32 }, + #[error("index {index} is duplicated in a single specifier")] + DuplicateIndex { index: u32 }, + #[error("labels must only contain letters from the alphabet 'IXYZ+-rl01'")] + OutsideAlphabet, +} +impl From for PyErr { + fn from(value: LabelError) -> PyErr { + PyValueError::new_err(value.to_string()) + } +} + +/// An observable over Pauli bases that stores its data in a qubit-sparse format. +/// +/// Mathematics +/// =========== +/// +/// This observable represents a sum over strings of the Pauli operators and Pauli-eigenstate +/// projectors, with each term weighted by some complex number. That is, the full observable is +/// +/// .. math:: +/// +/// \text{\texttt{SparseObservable}} = \sum_i c_i \bigotimes_n A^{(n)}_i +/// +/// for complex numbers :math:`c_i` and single-qubit operators acting on qubit :math:`n` from a +/// restricted alphabet :math:`A^{(n)}_i`. The sum over :math:`i` is the sum of the individual +/// terms, and the tensor product produces the operator strings. +/// +/// The alphabet of allowed single-qubit operators that the :math:`A^{(n)}_i` are drawn from is the +/// Pauli operators and the Pauli-eigenstate projection operators. Explicitly, these are: +/// +/// .. _sparse-observable-alphabet: +/// .. table:: Alphabet of single-qubit terms used in :class:`SparseObservable` +/// +/// ======= ======================================= =============== =========================== +/// Label Operator Numeric value :class:`.BitTerm` attribute +/// ======= ======================================= =============== =========================== +/// ``"I"`` :math:`I` (identity) Not stored. Not stored. +/// +/// ``"X"`` :math:`X` (Pauli X) ``0b0010`` (2) :attr:`~.BitTerm.X` +/// +/// ``"Y"`` :math:`Y` (Pauli Y) ``0b0011`` (3) :attr:`~.BitTerm.Y` +/// +/// ``"Z"`` :math:`Z` (Pauli Z) ``0b0001`` (1) :attr:`~.BitTerm.Z` +/// +/// ``"+"`` :math:`\lvert+\rangle\langle+\rvert` ``0b1010`` (10) :attr:`~.BitTerm.PLUS` +/// (projector to positive eigenstate of X) +/// +/// ``"-"`` :math:`\lvert-\rangle\langle-\rvert` ``0b0110`` (6) :attr:`~.BitTerm.MINUS` +/// (projector to negative eigenstate of X) +/// +/// ``"r"`` :math:`\lvert r\rangle\langle r\rvert` ``0b1011`` (11) :attr:`~.BitTerm.RIGHT` +/// (projector to positive eigenstate of Y) +/// +/// ``"l"`` :math:`\lvert l\rangle\langle l\rvert` ``0b0111`` (7) :attr:`~.BitTerm.LEFT` +/// (projector to negative eigenstate of Y) +/// +/// ``"0"`` :math:`\lvert0\rangle\langle0\rvert` ``0b1001`` (9) :attr:`~.BitTerm.ZERO` +/// (projector to positive eigenstate of Z) +/// +/// ``"1"`` :math:`\lvert1\rangle\langle1\rvert` ``0b0101`` (5) :attr:`~.BitTerm.ONE` +/// (projector to negative eigenstate of Z) +/// ======= ======================================= =============== =========================== +/// +/// The allowed alphabet forms an overcomplete basis of the operator space. This means that there +/// is not a unique summation to represent a given observable. By comparison, +/// :class:`.SparsePauliOp` uses a precise basis of the operator space, so (after combining terms of +/// the same Pauli string, removing zeros, and sorting the terms to some canonical order) there is +/// only one representation of any operator. +/// +/// :class:`SparseObservable` uses its particular overcomplete basis with the aim of making +/// "efficiency of measurement" equivalent to "efficiency of representation". For example, the +/// observable :math:`{\lvert0\rangle\langle0\rvert}^{\otimes n}` can be efficiently measured on +/// hardware with simple :math:`Z` measurements, but can only be represented by +/// :class:`.SparsePauliOp` as :math:`{(I + Z)}^{\otimes n}/2^n`, which requires :math:`2^n` stored +/// terms. :class:`SparseObservable` requires only a single term to store this. +/// +/// The downside to this is that it is impractical to take an arbitrary matrix or +/// :class:`.SparsePauliOp` and find the *best* :class:`SparseObservable` representation. You +/// typically will want to construct a :class:`SparseObservable` directly, rather than trying to +/// decompose into one. +/// +/// +/// Representation +/// ============== +/// +/// The internal representation of a :class:`SparseObservable` stores only the non-identity qubit +/// operators. This makes it significantly more efficient to represent observables such as +/// :math:`\sum_{n\in \text{qubits}} Z^{(n)}`; :class:`SparseObservable` requires an amount of +/// memory linear in the total number of qubits, while :class:`.SparsePauliOp` scales quadratically. +/// +/// The terms are stored compressed, similar in spirit to the compressed sparse row format of sparse +/// matrices. In this analogy, the terms of the sum are the "rows", and the qubit terms are the +/// "columns", where an absent entry represents the identity rather than a zero. More explicitly, +/// the representation is made up of four contiguous arrays: +/// +/// .. _sparse-observable-arrays: +/// .. table:: Data arrays used to represent :class:`.SparseObservable` +/// +/// ================== =========== ============================================================= +/// Attribute Length Description +/// ================== =========== ============================================================= +/// :attr:`coeffs` :math:`t` The complex scalar multiplier for each term. +/// +/// :attr:`bit_terms` :math:`s` Each of the non-identity single-qubit terms for all of the +/// operators, in order. These correspond to the non-identity +/// :math:`A^{(n)}_i` in the sum description, where the entries +/// are stored in order of increasing :math:`i` first, and in +/// order of increasing :math:`n` within each term. +/// +/// :attr:`indices` :math:`s` The corresponding qubit (:math:`n`) for each of the operators +/// in :attr:`bit_terms`. :class:`SparseObservable` requires +/// that this list is term-wise sorted, and algorithms can rely +/// on this invariant being upheld. +/// +/// :attr:`boundaries` :math:`t+1` The indices that partition :attr:`bit_terms` and +/// :attr:`indices` into complete terms. For term number +/// :math:`i`, its complex coefficient is ``coeffs[i]``, and its +/// non-identity single-qubit operators and their corresponding +/// qubits are the slice ``boundaries[i] : boundaries[i+1]`` into +/// :attr:`bit_terms` and :attr:`indices` respectively. +/// :attr:`boundaries` always has an explicit 0 as its first +/// element. +/// ================== =========== ============================================================= +/// +/// The length parameter :math:`t` is the number of terms in the sum, and the parameter :math:`s` is +/// the total number of non-identity single-qubit terms. +/// +/// As illustrative examples: +/// +/// * in the case of a zero operator, :attr:`boundaries` is length 1 (a single 0) and all other +/// vectors are empty. +/// * in the case of a fully simplified identity operator, :attr:`boundaries` is ``[0, 0]``, +/// :attr:`coeffs` has a single entry, and :attr:`bit_terms` and :attr:`indices` are empty. +/// +/// These two cases are not special, they're fully consistent with the rules and should not need +/// special handling. +/// +/// The scalar item of the :attr:`bit_terms` array is stored as a numeric byte. The numeric values +/// are related to the symplectic Pauli representation that :class:`.SparsePauliOp` uses, and are +/// accessible with named access by an enumeration: +/// +/// .. +/// This is documented manually here because the Python-space `Enum` is generated +/// programmatically from Rust - it'd be _more_ confusing to try and write a docstring somewhere +/// else in this source file. The use of `autoattribute` is because it pulls in the numeric +/// value. +/// +/// .. py:class:: SparseObservable.BitTerm +/// +/// An :class:`~enum.IntEnum` that provides named access to the numerical values used to +/// represent each of the single-qubit alphabet terms enumerated in +/// :ref:`sparse-observable-alphabet`. +/// +/// This class is attached to :class:`.SparseObservable`. Access it as +/// :class:`.SparseObservable.BitTerm`. +/// +/// The numeric structure of these is that they are all four-bit values of which the low two +/// bits are the (phase-less) symplectic representation of the Pauli operator related to the +/// object, where the low bit denotes a contribution by :math:`Z` and the second lowest a +/// contribution by :math:`X`, while the upper two bits are ``00`` for a Pauli operator, ``01`` +/// for the negative-eigenstate projector, and ``10`` for teh positive-eigenstate projector. +/// +/// .. autoattribute:: qiskit.quantum_info::SparseObservable.BitTerm.X +/// +/// The Pauli :math:`X` operator. +/// +/// .. autoattribute:: qiskit.quantum_info::SparseObservable.BitTerm.PLUS +/// +/// The projector to the positive eigenstate of the :math:`X` operator: +/// :math:`\lvert+\rangle\langle+\rvert`. Uses the single-letter label ``"+"``. +/// +/// .. autoattribute:: qiskit.quantum_info::SparseObservable.BitTerm.MINUS +/// +/// The projector to the negative eigenstate of the :math:`X` operator: +/// :math:`\lvert-\rangle\langle-\rvert`. Uses the single-letter label ``"-"``. +/// +/// .. autoattribute:: qiskit.quantum_info::SparseObservable.BitTerm.Y +/// +/// The Pauli :math:`Y` operator. +/// +/// .. autoattribute:: qiskit.quantum_info::SparseObservable.BitTerm.RIGHT +/// +/// The projector to the positive eigenstate of the :math:`Y` operator: +/// :math:`\lvert r\rangle\langle r\rvert`. Uses the single-letter label ``"r"``. +/// +/// .. autoattribute:: qiskit.quantum_info::SparseObservable.BitTerm.LEFT +/// +/// The projector to the negative eigenstate of the :math:`Y` operator: +/// :math:`\lvert l\rangle\langle l\rvert`. Uses the single-letter label ``"l"``. +/// +/// .. autoattribute:: qiskit.quantum_info::SparseObservable.BitTerm.Z +/// +/// The Pauli :math:`Z` operator. +/// +/// .. autoattribute:: qiskit.quantum_info::SparseObservable.BitTerm.ZERO +/// +/// The projector to the positive eigenstate of the :math:`Z` operator: +/// :math:`\lvert0\rangle\langle0\rvert`. Uses the single-letter label ``"0"``. +/// +/// .. autoattribute:: qiskit.quantum_info::SparseObservable.BitTerm.ONE +/// +/// The projector to the negative eigenstate of the :math:`Z` operator: +/// :math:`\lvert1\rangle\langle1\rvert`. Uses the single-letter label ``"1"``. +/// +/// Each of the array-like attributes behaves like a Python sequence. You can index and slice these +/// with standard :class:`list`-like semantics. Slicing an attribute returns a Numpy +/// :class:`~numpy.ndarray` containing a copy of the relevant data with the natural ``dtype`` of the +/// field; this lets you easily do mathematics on the results, like bitwise operations on +/// :attr:`bit_terms`. You can assign to indices or slices of each of the attributes, but beware +/// that you must uphold :ref:`the data coherence rules ` while doing +/// this. For example:: +/// +/// >>> obs = SparseObservable.from_list([("XZY", 1.5j), ("+1r", -0.5)]) +/// >>> assert isinstance(obs.coeffs[:], np.ndarray) +/// >>> # Reduce all single-qubit terms to the relevant Pauli operator, if they are a projector. +/// >>> obs.bit_terms[:] = obs.bit_terms[:] & 0b00_11 +/// >>> assert obs == SparseObservable.from_list([("XZY", 1.5j), ("XZY", -0.5)]) +/// +/// +/// Construction +/// ============ +/// +/// :class:`SparseObservable` defines several constructors. The default constructor will attempt to +/// delegate to one of the more specific constructors, based on the type of the input. You can +/// always use the specific constructors to have more control over the construction. +/// +/// .. _sparse-observable-convert-constructors: +/// .. table:: Construction from other objects +/// +/// ============================ ================================================================ +/// Method Summary +/// ============================ ================================================================ +/// :meth:`from_label` Convert a dense string label into a single-term +/// :class:`.SparseObservable`. +/// +/// :meth:`from_list` Sum a list of tuples of dense string labels and the associated +/// coefficients into an observable. +/// +/// :meth:`from_sparse_list` Sum a list of tuples of sparse string labels, the qubits they +/// apply to, and their coefficients into an observable. +/// +/// :meth:`from_pauli` Raise a single :class:`.Pauli` into a single-term +/// :class:`.SparseObservable`. +/// +/// :meth:`from_sparse_pauli_op` Raise a :class:`.SparsePauliOp` into a :class:`SparseObservable`. +/// +/// :meth:`from_raw_parts` Build the observable from :ref:`the raw data arrays +/// `. +/// ============================ ================================================================ +/// +/// .. py:function:: SparseObservable.__new__(data, /, num_qubits=None) +/// +/// The default constructor of :class:`SparseObservable`. +/// +/// This delegates to one of :ref:`the explicit conversion-constructor methods +/// `, based on the type of the ``data`` argument. If +/// ``num_qubits`` is supplied and constructor implied by the type of ``data`` does not accept a +/// number, the given integer must match the input. +/// +/// :param data: The data type of the input. This can be another :class:`SparseObservable`, in +/// which case the input is copied, a :class:`.Pauli` or :class:`.SparsePauliOp`, in which +/// case :meth:`from_pauli` or :meth:`from_sparse_pauli_op` are called as appropriate, or it +/// can be a list in a valid format for either :meth:`from_list` or +/// :meth:`from_sparse_list`. +/// :param int|None num_qubits: Optional number of qubits for the operator. For most data +/// inputs, this can be inferred and need not be passed. It is only necessary for empty +/// lists or the sparse-list format. If given unnecessarily, it must match the data input. +/// +/// In addition to the conversion-based constructors, there are also helper methods that construct +/// special forms of observables. +/// +/// .. table:: Construction of special observables +/// +/// ============================ ================================================================ +/// Method Summary +/// ============================ ================================================================ +/// :meth:`zero` The zero operator on a given number of qubits. +/// +/// :meth:`identity` The identity operator on a given number of qubits. +/// ============================ ================================================================ +#[pyclass(module = "qiskit.quantum_info")] +#[derive(Clone, Debug, PartialEq)] +pub struct SparseObservable { + /// The number of qubits the operator acts on. This is not inferable from any other shape or + /// values, since identities are not stored explicitly. + num_qubits: u32, + /// The coefficients of each abstract term in in the sum. This has as many elements as terms in + /// the sum. + coeffs: Vec, + /// A flat list of single-qubit terms. This is more naturally a list of lists, but is stored + /// flat for memory usage and locality reasons, with the sublists denoted by `boundaries.` + bit_terms: Vec, + /// A flat list of the qubit indices that the corresponding entries in `bit_terms` act on. This + /// list must always be term-wise sorted, where a term is a sublist as denoted by `boundaries`. + indices: Vec, + /// Indices that partition `bit_terms` and `indices` into sublists for each individual term in + /// the sum. `boundaries[0]..boundaries[1]` is the range of indices into `bit_terms` and + /// `indices` that correspond to the first term of the sum. All unspecified qubit indices are + /// implicitly the identity. This is one item longer than `coeffs`, since `boundaries[0]` is + /// always an explicit zero (for algorithmic ease). + boundaries: Vec, +} + +impl SparseObservable { + /// Create a new observable from the raw components that make it up. + /// + /// This checks the input values for data coherence on entry. If you are certain you have the + /// correct values, you can call `new_unchecked` instead. + pub fn new( + num_qubits: u32, + coeffs: Vec, + bit_terms: Vec, + indices: Vec, + boundaries: Vec, + ) -> Result { + if coeffs.len() + 1 != boundaries.len() { + return Err(CoherenceError::MismatchedTermCount { + coeffs: coeffs.len(), + boundaries: boundaries.len(), + }); + } + if bit_terms.len() != indices.len() { + return Err(CoherenceError::MismatchedItemCount { + bit_terms: bit_terms.len(), + indices: indices.len(), + }); + } + // We already checked that `boundaries` is at least length 1. + if *boundaries.first().unwrap() != 0 { + return Err(CoherenceError::BadInitialBoundary(boundaries[0])); + } + if *boundaries.last().unwrap() != indices.len() { + return Err(CoherenceError::BadFinalBoundary { + last: *boundaries.last().unwrap(), + items: indices.len(), + }); + } + for (&left, &right) in boundaries[..].iter().zip(&boundaries[1..]) { + if right < left { + return Err(CoherenceError::DecreasingBoundaries); + } + let indices = &indices[left..right]; + if !indices.is_empty() { + for (index_left, index_right) in indices[..].iter().zip(&indices[1..]) { + if index_left >= index_right { + return Err(CoherenceError::UnsortedIndices); + } + } + } + if indices.last().map(|&ix| ix >= num_qubits).unwrap_or(false) { + return Err(CoherenceError::BitIndexTooHigh); + } + } + // SAFETY: we've just done the coherence checks. + Ok(unsafe { Self::new_unchecked(num_qubits, coeffs, bit_terms, indices, boundaries) }) + } + + /// Create a new observable from the raw components without checking data coherence. + /// + /// # Safety + /// + /// It is up to the caller to ensure that the data-coherence requirements, as enumerated in the + /// struct-level documentation, have been upheld. + #[inline(always)] + pub unsafe fn new_unchecked( + num_qubits: u32, + coeffs: Vec, + bit_terms: Vec, + indices: Vec, + boundaries: Vec, + ) -> Self { + Self { + num_qubits, + coeffs, + bit_terms, + indices, + boundaries, + } + } + + /// Get an iterator over the individual terms of the operator. + pub fn iter(&'_ self) -> impl ExactSizeIterator> + '_ { + self.coeffs.iter().enumerate().map(|(i, coeff)| { + let start = self.boundaries[i]; + let end = self.boundaries[i + 1]; + SparseTermView { + coeff: *coeff, + bit_terms: &self.bit_terms[start..end], + indices: &self.indices[start..end], + } + }) + } + + /// Add the term implied by a dense string label onto this observable. + pub fn add_dense_label>( + &mut self, + label: L, + coeff: Complex64, + ) -> Result<(), LabelError> { + let label = label.as_ref(); + if label.len() != self.num_qubits() as usize { + return Err(LabelError::WrongLengthDense { + num_qubits: self.num_qubits(), + label: label.len(), + }); + } + // The only valid characeters in the alphabet are ASCII, so if we see something other than + // ASCII, we're already in the failure path. + for (i, letter) in label.iter().rev().enumerate() { + match BitTerm::try_from_u8(*letter) { + Ok(Some(term)) => { + self.bit_terms.push(term); + self.indices.push(i as u32); + } + Ok(None) => (), + Err(_) => { + // Undo any modifications to ourselves so we stay in a consistent state. + let num_single_terms = self.boundaries[self.boundaries.len() - 1]; + self.bit_terms.truncate(num_single_terms); + self.indices.truncate(num_single_terms); + return Err(LabelError::OutsideAlphabet); + } + } + } + self.coeffs.push(coeff); + self.boundaries.push(self.bit_terms.len()); + Ok(()) + } +} + +#[pymethods] +impl SparseObservable { + #[pyo3(signature = (data, /, num_qubits=None))] + #[new] + fn py_new(data: Bound, num_qubits: Option) -> PyResult { + let py = data.py(); + let check_num_qubits = |data: &Bound| -> PyResult<()> { + let Some(num_qubits) = num_qubits else { + return Ok(()); + }; + let other_qubits = data.getattr(intern!(py, "num_qubits"))?.extract::()?; + if num_qubits == other_qubits { + return Ok(()); + } + Err(PyValueError::new_err(format!( + "explicitly given 'num_qubits' ({num_qubits}) does not match operator ({other_qubits})" + ))) + }; + + if data.is_instance(PAULI_TYPE.get_bound(py))? { + check_num_qubits(&data)?; + return Self::py_from_pauli(data); + } + if data.is_instance(SPARSE_PAULI_OP_TYPE.get_bound(py))? { + check_num_qubits(&data)?; + return Self::py_from_sparse_pauli_op(data); + } + if let Ok(label) = data.extract::() { + let num_qubits = num_qubits.unwrap_or(label.len() as u32); + if num_qubits as usize != label.len() { + return Err(PyValueError::new_err(format!( + "explicitly given 'num_qubits' ({}) does not match label ({})", + num_qubits, + label.len(), + ))); + } + return Self::py_from_label(&label).map_err(PyErr::from); + } + if let Ok(observable) = data.downcast::() { + check_num_qubits(&data)?; + return Ok(observable.borrow().clone()); + } + if let Ok(vec) = data.extract() { + return Self::py_from_list(vec, num_qubits); + } + if let Ok(vec) = data.extract() { + let Some(num_qubits) = num_qubits else { + return Err(PyValueError::new_err( + "if using the sparse-list form, 'num_qubits' must be provided", + )); + }; + return Self::py_from_sparse_list(vec, num_qubits).map_err(PyErr::from); + } + Err(PyTypeError::new_err(format!( + "unknown input format for 'SparseObservable': {}", + data.get_type().repr()?, + ))) + } + + /// The number of qubits the operator acts on. + /// + /// This is not inferable from any other shape or values, since identities are not stored + /// explicitly. + #[getter] + #[inline] + pub fn num_qubits(&self) -> u32 { + self.num_qubits + } + + /// The number of terms in the sum this operator is tracking. + #[getter] + #[inline] + pub fn num_ops(&self) -> usize { + self.coeffs.len() + } + + /// The coefficients of each abstract term in in the sum. This has as many elements as terms in + /// the sum. + #[getter] + fn get_coeffs(slf_: Py) -> ArrayView { + ArrayView { + base: slf_, + slot: ArraySlot::Coeffs, + } + } + + /// A flat list of single-qubit terms. This is more naturally a list of lists, but is stored + /// flat for memory usage and locality reasons, with the sublists denoted by `boundaries.` + #[getter] + fn get_bit_terms(slf_: Py) -> ArrayView { + ArrayView { + base: slf_, + slot: ArraySlot::BitTerms, + } + } + + /// A flat list of the qubit indices that the corresponding entries in :attr:`bit_terms` act on. + /// This list must always be term-wise sorted, where a term is a sublist as denoted by + /// :attr:`boundaries`. + /// + /// .. warning:: + /// + /// If writing to this attribute from Python space, you *must* ensure that you only write in + /// indices that are term-wise sorted. + #[getter] + fn get_indices(slf_: Py) -> ArrayView { + ArrayView { + base: slf_, + slot: ArraySlot::Indices, + } + } + + /// Indices that partition :attr:`bit_terms` and :attr:`indices` into sublists for each + /// individual term in the sum. ``boundaries[0] : boundaries[1]`` is the range of indices into + /// :attr:`bit_terms` and :attr:`indices` that correspond to the first term of the sum. All + /// unspecified qubit indices are implicitly the identity. This is one item longer than + /// :attr:`coeffs`, since ``boundaries[0]`` is always an explicit zero (for algorithmic ease). + #[getter] + fn get_boundaries(slf_: Py) -> ArrayView { + ArrayView { + base: slf_, + slot: ArraySlot::Boundaries, + } + } + + // The documentation for this is inlined into the class-level documentation of + // `SparseObservable`. + #[allow(non_snake_case)] + #[classattr] + fn BitTerm(py: Python) -> PyResult> { + BIT_TERM_PY_ENUM + .get_or_try_init(py, || make_py_bit_term(py)) + .map(|obj| obj.clone_ref(py)) + } + + /// Get the zero operator over the given number of qubits. + /// + /// Examples: + /// + /// Get the zero operator for 100 qubits:: + /// + /// >>> SparseObservable.zero(100) + /// + #[pyo3(signature = (/, num_qubits))] + #[staticmethod] + pub fn zero(num_qubits: u32) -> Self { + Self { + num_qubits, + coeffs: vec![], + bit_terms: vec![], + indices: vec![], + boundaries: vec![0], + } + } + + /// Get the identity operator over the given number of qubits. + /// + /// Examples: + /// + /// Get the identity operator for 100 qubits:: + /// + /// >>> SparseObservable.identity(100) + /// + #[pyo3(signature = (/, num_qubits))] + #[staticmethod] + pub fn identity(num_qubits: u32) -> Self { + Self { + num_qubits, + coeffs: vec![Complex64::new(1.0, 0.0)], + bit_terms: vec![], + indices: vec![], + boundaries: vec![0, 0], + } + } + + fn __repr__(&self) -> String { + let num_terms = format!( + "{} term{}", + self.num_ops(), + if self.num_ops() == 1 { "" } else { "s" } + ); + let qubits = format!( + "{} qubit{}", + self.num_qubits(), + if self.num_qubits() == 1 { "" } else { "s" } + ); + let terms = if self.num_ops() == 0 { + "0.0".to_owned() + } else { + self.iter() + .map(|term| { + let coeff = format!("{}", term.coeff).replace("i", "j"); + let paulis = term + .indices + .iter() + .zip(term.bit_terms) + .rev() + .map(|(i, op)| format!("{}_{}", op.py_label(), i)) + .collect::>() + .join(" "); + format!("({})({})", coeff, paulis) + }) + .collect::>() + .join(" + ") + }; + format!( + "", + num_terms, qubits, terms + ) + } + + fn __reduce__(&self, py: Python) -> PyResult> { + let bit_terms: &[u8] = ::bytemuck::cast_slice(&self.bit_terms); + Ok(( + py.get_type_bound::().getattr("from_raw_parts")?, + ( + self.num_qubits, + PyArray1::from_slice_bound(py, &self.coeffs), + PyArray1::from_slice_bound(py, bit_terms), + PyArray1::from_slice_bound(py, &self.indices), + PyArray1::from_slice_bound(py, &self.boundaries), + false, + ), + ) + .into_py(py)) + } + + fn __eq__(slf: Bound, other: Bound) -> bool { + if slf.is(&other) { + return true; + } + let Ok(other) = other.downcast_into::() else { + return false; + }; + slf.borrow().eq(&other.borrow()) + } + + // This doesn't actually have any interaction with Python space, but uses the `py_` prefix on + // its name to make it clear it's different to the Rust concept of `Copy`. + /// Get a copy of this observable. + /// + /// Examples: + /// + /// .. code-block:: python + /// + /// >>> obs = SparseObservable.from_list([("IXZ+lr01", 2.5), ("ZXI-rl10", 0.5j)]) + /// >>> assert obs == obs.copy() + /// >>> assert obs is not obs.copy() + #[pyo3(name = "copy")] + fn py_copy(&self) -> Self { + self.clone() + } + + /// Construct a single-term observable from a dense string label. + /// + /// The resulting operator will have a coefficient of 1. The label must be a sequence of the + /// alphabet ``'IXYZ+-rl01'``. The label is interpreted analagously to a bitstring. In other + /// words, the right-most letter is associated with qubit 0, and so on. This is the same as the + /// labels for :class:`.Pauli` and :class:`.SparsePauliOp`. + /// + /// Args: + /// label (str): the dense label. + /// + /// Examples: + /// + /// .. code-block:: python + /// + /// >>> SparseObservable.from_label("IIII+ZI") + /// + /// >>> label = "IYXZI" + /// >>> pauli = Pauli(label) + /// >>> assert SparseObservable.from_label(label) == SparseObservable.from_pauli(pauli) + /// + /// See also: + /// :meth:`from_list` + /// A generalization of this method that constructs a sum operator from multiple labels + /// and their corresponding coefficients. + #[staticmethod] + #[pyo3(name = "from_label", signature = (label, /))] + fn py_from_label(label: &str) -> Result { + let mut out = Self::zero(label.len() as u32); + out.add_dense_label(label, Complex64::new(1.0, 0.0))?; + Ok(out) + } + + /// Construct an observable from a list of dense labels and coefficients. + /// + /// This is analagous to :meth:`.SparsePauliOp.from_list`, except it uses + /// :ref:`the extended alphabet ` of :class:`.SparseObservable`. In + /// this dense form, you must supply all identities explicitly in each label. + /// + /// The label must be a sequence of the alphabet ``'IXYZ+-rl01'``. The label is interpreted + /// analagously to a bitstring. In other words, the right-most letter is associated with qubit + /// 0, and so on. This is the same as the labels for :class:`.Pauli` and + /// :class:`.SparsePauliOp`. + /// + /// Args: + /// iter (list[tuple[str, complex]]): Pairs of labels and their associated coefficients to + /// sum. The labels are interpreted the same way as in :meth:`from_label`. + /// num_qubits (int | None): It is not necessary to specify this if you are sure that + /// ``iter`` is not an empty sequence, since it can be inferred from the label lengths. + /// If ``iter`` may be empty, you must specify this argument to disambiguate how many + /// qubits the observable is for. If this is given and ``iter`` is not empty, the value + /// must match the label lengths. + /// + /// Examples: + /// + /// Construct an observable from a list of labels of the same length:: + /// + /// >>> SparseObservable.from_list([ + /// ... ("III++", 1.0), + /// ... ("II--I", 1.0j), + /// ... ("I++II", -0.5), + /// ... ("--III", -0.25j), + /// ... ]) + /// + /// + /// Use ``num_qubits`` to disambiguate potentially empty inputs:: + /// + /// >>> SparseObservable.from_list([], num_qubits=10) + /// + /// + /// This method is equivalent to calls to :meth:`from_sparse_list` with the explicit + /// qubit-arguments field set to decreasing integers:: + /// + /// >>> labels = ["XY+Z", "rl01", "-lXZ"] + /// >>> coeffs = [1.5j, 2.0, -0.5] + /// >>> from_list = SparseObservable.from_list(list(zip(labels, coeffs))) + /// >>> from_sparse_list = SparseObservable.from_sparse_list([ + /// ... (label, (3, 2, 1, 0), coeff) + /// ... for label, coeff in zip(labels, coeffs) + /// ... ]) + /// >>> assert from_list == from_sparse_list + /// + /// See also: + /// :meth:`from_label` + /// A similar constructor, but takes only a single label and always has its coefficient + /// set to ``1.0``. + /// + /// :meth:`from_sparse_list` + /// Construct the observable from a list of labels without explicit identities, but with + /// the qubits each single-qubit term applies to listed explicitly. + #[staticmethod] + #[pyo3(name = "from_list", signature = (iter, /, *, num_qubits=None))] + fn py_from_list(iter: Vec<(String, Complex64)>, num_qubits: Option) -> PyResult { + if iter.is_empty() && num_qubits.is_none() { + return Err(PyValueError::new_err( + "cannot construct an observable from an empty list without knowing `num_qubits`", + )); + } + let num_qubits = match num_qubits { + Some(num_qubits) => num_qubits, + None => iter[0].0.len() as u32, + }; + let mut out = Self { + num_qubits, + coeffs: Vec::with_capacity(iter.len()), + bit_terms: Vec::new(), + indices: Vec::new(), + boundaries: Vec::with_capacity(iter.len() + 1), + }; + out.boundaries.push(0); + for (label, coeff) in iter { + out.add_dense_label(&label, coeff)?; + } + Ok(out) + } + + /// Construct an observable from a list of labels, the qubits each item applies to, and the + /// coefficient of the whole term. + /// + /// This is analagous to :meth:`.SparsePauliOp.from_sparse_list`, except it uses + /// :ref:`the extended alphabet ` of :class:`.SparseObservable`. + /// + /// The "labels" and "indices" fields of the triples are associated by zipping them together. + /// For example, this means that a call to :meth:`from_list` can be converted to the form used + /// by this method by setting the "indices" field of each triple to ``(num_qubits-1, ..., 1, + /// 0)``. + /// + /// Args: + /// iter (list[tuple[str, Sequence[int], complex]]): triples of labels, the qubits + /// each single-qubit term applies to, and the coefficient of the entire term. + /// + /// num_qubits (int): the number of qubits in the operator. + /// + /// Examples: + /// + /// Construct a simple operator:: + /// + /// >>> SparseObservable.from_sparse_list( + /// ... [("ZX", (1, 4), 1.0), ("YY", (0, 3), 2j)], + /// ... num_qubits=5, + /// ... ) + /// + /// + /// Construct the identity observable (though really, just use :meth:`identity`):: + /// + /// >>> SparseObservable.from_sparse_list([("", (), 1.0)], num_qubits=100) + /// + /// + /// This method can replicate the behavior of :meth:`from_list`, if the qubit-arguments + /// field of the triple is set to decreasing integers:: + /// + /// >>> labels = ["XY+Z", "rl01", "-lXZ"] + /// >>> coeffs = [1.5j, 2.0, -0.5] + /// >>> from_list = SparseObservable.from_list(list(zip(labels, coeffs))) + /// >>> from_sparse_list = SparseObservable.from_sparse_list([ + /// ... (label, (3, 2, 1, 0), coeff) + /// ... for label, coeff in zip(labels, coeffs) + /// ... ]) + /// >>> assert from_list == from_sparse_list + #[staticmethod] + #[pyo3(name = "from_sparse_list", signature = (iter, /, num_qubits))] + fn py_from_sparse_list( + iter: Vec<(String, Vec, Complex64)>, + num_qubits: u32, + ) -> Result { + let coeffs = iter.iter().map(|(_, _, coeff)| *coeff).collect(); + let mut boundaries = Vec::with_capacity(iter.len() + 1); + boundaries.push(0); + let mut indices = Vec::new(); + let mut bit_terms = Vec::new(); + let mut sorted = btree_map::BTreeMap::new(); + for (label, qubits, _) in iter { + sorted.clear(); + let label: &[u8] = label.as_ref(); + if label.len() != qubits.len() { + return Err(LabelError::WrongLengthIndices { + label: label.len(), + indices: indices.len(), + }); + } + for (letter, index) in label.iter().zip(qubits) { + if index >= num_qubits { + return Err(LabelError::BadIndex { index, num_qubits }); + } + let btree_map::Entry::Vacant(entry) = sorted.entry(index) else { + return Err(LabelError::DuplicateIndex { index }); + }; + entry.insert( + BitTerm::try_from_u8(*letter).map_err(|_| LabelError::OutsideAlphabet)?, + ); + } + for (index, term) in sorted.iter() { + let Some(term) = term else { + continue; + }; + indices.push(*index); + bit_terms.push(*term); + } + boundaries.push(bit_terms.len()); + } + Ok(Self { + num_qubits, + coeffs, + indices, + bit_terms, + boundaries, + }) + } + + /// Construct a :class:`.SparseObservable` from a single :class:`.Pauli` instance. + /// + /// The output observable will have a single term, with a unitary coefficient dependent on the + /// phase. + /// + /// Args: + /// pauli (:class:`.Pauli`): the single Pauli to convert. + /// + /// Examples: + /// + /// .. code-block:: python + /// + /// >>> label = "IYXZI" + /// >>> pauli = Pauli(label) + /// >>> SparseObservable.from_pauli(pauli) + /// + /// >>> assert SparseObservable.from_label(label) == SparseObservable.from_pauli(pauli) + #[staticmethod] + #[pyo3(name = "from_pauli", signature = (pauli, /))] + fn py_from_pauli(pauli: Bound) -> PyResult { + let py = pauli.py(); + let num_qubits = pauli.getattr(intern!(py, "num_qubits"))?.extract::()?; + let z = pauli + .getattr(intern!(py, "z"))? + .extract::>()?; + let x = pauli + .getattr(intern!(py, "x"))? + .extract::>()?; + let mut bit_terms = Vec::new(); + let mut indices = Vec::new(); + let mut num_ys = 0; + for (i, (x, z)) in x.as_array().iter().zip(z.as_array().iter()).enumerate() { + // The only failure case possible here is the identity, because of how we're + // constructing the value to convert. + let Ok(term) = ::bytemuck::checked::try_cast((*x as u8) << 1 | (*z as u8)) else { + continue; + }; + num_ys += (term == BitTerm::Y) as isize; + indices.push(i as u32); + bit_terms.push(term); + } + let boundaries = vec![0, indices.len()]; + // The "empty" state of a `Pauli` represents the identity, which isn't our empty state + // (that's zero), so we're always going to have a coefficient. + let group_phase = pauli + // `Pauli`'s `_phase` is a Numpy array ... + .getattr(intern!(py, "_phase"))? + // ... that should have exactly 1 element ... + .call_method0(intern!(py, "item"))? + // ... which is some integral type. + .extract::()?; + let phase = match (group_phase - num_ys).rem_euclid(4) { + 0 => Complex64::new(1.0, 0.0), + 1 => Complex64::new(0.0, -1.0), + 2 => Complex64::new(-1.0, 0.0), + 3 => Complex64::new(0.0, 1.0), + _ => unreachable!("`x % 4` has only four values"), + }; + let coeffs = vec![phase]; + Ok(Self { + num_qubits, + coeffs, + bit_terms, + indices, + boundaries, + }) + } + + /// Construct a :class:`.SparseObservable` from a :class:`.SparsePauliOp` instance. + /// + /// This will be a largely direct translation of the :class:`.SparsePauliOp`; in particular, + /// there is no on-the-fly summing of like terms, nor any attempt to refactorize sums of Pauli + /// terms into equivalent projection operators. + /// + /// Args: + /// op (:class:`.SparsePauliOp`): the operator to convert. + /// + /// Examples: + /// + /// .. code-block:: python + /// + /// >>> spo = SparsePauliOp.from_list([("III", 1.0), ("IIZ", 0.5), ("IZI", 0.5)]) + /// >>> SparseObservable.from_sparse_pauli_op(spo) + /// + #[staticmethod] + #[pyo3(name = "from_sparse_pauli_op", signature = (op, /))] + fn py_from_sparse_pauli_op(op: Bound) -> PyResult { + let py = op.py(); + let pauli_list_ob = op.getattr(intern!(py, "paulis"))?; + let coeffs = op + .getattr(intern!(py, "coeffs"))? + .extract::>() + .map_err(|_| PyTypeError::new_err("only 'SparsePauliOp' with complex-typed coefficients can be converted to 'SparseObservable"))? + .as_array() + .to_vec(); + let op_z = pauli_list_ob + .getattr(intern!(py, "z"))? + .extract::>()?; + let op_x = pauli_list_ob + .getattr(intern!(py, "x"))? + .extract::>()?; + // We don't extract the `phase`, because that's supposed to be 0 for all `SparsePauliOp` + // instances - they use the symplectic convention in the representation with any phase term + // absorbed into the coefficients (like us). + let [num_ops, num_qubits] = *op_z.shape() else { + unreachable!("shape is statically known to be 2D") + }; + if op_x.shape() != [num_ops, num_qubits] { + return Err(PyValueError::new_err(format!( + "'x' and 'z' have different shapes ({:?} and {:?})", + op_x.shape(), + op_z.shape() + ))); + } + if num_ops != coeffs.len() { + return Err(PyValueError::new_err(format!( + "'x' and 'z' have a different number of operators to 'coeffs' ({} and {})", + num_ops, + coeffs.len(), + ))); + } + + let mut bit_terms = Vec::new(); + let mut indices = Vec::new(); + let mut boundaries = Vec::with_capacity(num_ops + 1); + boundaries.push(0); + for (term_x, term_z) in op_x + .as_array() + .rows() + .into_iter() + .zip(op_z.as_array().rows()) + { + for (i, (x, z)) in term_x.iter().zip(term_z.iter()).enumerate() { + // The only failure case possible here is the identity, because of how we're + // constructing the value to convert. + let Ok(term) = ::bytemuck::checked::try_cast((*x as u8) << 1 | (*z as u8)) else { + continue; + }; + indices.push(i as u32); + bit_terms.push(term); + } + boundaries.push(indices.len()); + } + + Ok(Self { + num_qubits: num_qubits as u32, + coeffs, + bit_terms, + indices, + boundaries, + }) + } + + // SAFETY: this cannot invoke undefined behaviour if `check = true`, but if `check = false` then + // the `bit_terms` must all be valid `BitTerm` representations. + /// Construct a :class:`.SparseObservable` from raw Numpy arrays that match :ref:`the required + /// data representation described in the class-level documentation `. + /// + /// The data from each array is copied into fresh, growable Rust-space allocations. + /// + /// Args: + /// num_qubits: number of qubits in the observable. + /// coeffs: complex coefficients of each term of the observable. This should be a Numpy + /// array with dtype :attr:`~numpy.complex128`. + /// bit_terms: flattened list of the single-qubit terms comprising all complete terms. This + /// should be a Numpy array with dtype :attr:`~numpy.uint8` (which is compatible with + /// :class:`.BitTerm`). + /// indices: flattened term-wise sorted list of the qubits each single-qubit term corresponds + /// to. This should be a Numpy array with dtype :attr:`~numpy.uint32`. + /// boundaries: the indices that partition ``bit_terms`` and ``indices`` into terms. This + /// should be a Numpy array with dtype :attr:`~numpy.uintp`. + /// check: if ``True`` (the default), validate that the data satisfies all coherence + /// guarantees. If ``False``, no checks are done. + /// + /// .. warning:: + /// + /// If ``check=False``, the ``bit_terms`` absolutely *must* be all be valid values + /// of :class:`.SparseObservable.BitTerm`. If they are not, Rust-space undefined + /// behavior may occur, entirely invalidating the program execution. + /// + /// Examples: + /// + /// Construct a sum of :math:`Z` on each individual qubit:: + /// + /// >>> num_qubits = 100 + /// >>> terms = np.full((num_qubits,), SparseObservable.BitTerm.Z, dtype=np.uint8) + /// >>> indices = np.arange(num_qubits, dtype=np.uint32) + /// >>> coeffs = np.ones((num_qubits,), dtype=complex) + /// >>> boundaries = np.arange(num_qubits + 1, dtype=np.uintp) + /// >>> SparseObservable.from_raw_parts(num_qubits, coeffs, terms, indices, boundaries) + /// + #[deny(unsafe_op_in_unsafe_fn)] + #[staticmethod] + #[pyo3( + signature = (/, num_qubits, coeffs, bit_terms, indices, boundaries, check=true), + name = "from_raw_parts", + )] + unsafe fn py_from_raw_parts<'py>( + num_qubits: u32, + coeffs: PyArrayLike1<'py, Complex64>, + bit_terms: PyArrayLike1<'py, u8>, + indices: PyArrayLike1<'py, u32>, + boundaries: PyArrayLike1<'py, usize>, + check: bool, + ) -> PyResult { + let coeffs = coeffs.as_array().to_vec(); + let bit_terms = if check { + bit_terms + .as_array() + .into_iter() + .copied() + .map(BitTerm::try_from) + .collect::>()? + } else { + let bit_terms_as_u8 = bit_terms.as_array().to_vec(); + // SAFETY: the caller enforced that each `u8` is a valid `BitTerm`, and `BitTerm` is be + // represented by a `u8`. We can't use `bytemuck` because we're casting a `Vec`. + unsafe { ::std::mem::transmute::, Vec>(bit_terms_as_u8) } + }; + let indices = indices.as_array().to_vec(); + let boundaries = boundaries.as_array().to_vec(); + + if check { + Self::new(num_qubits, coeffs, bit_terms, indices, boundaries).map_err(PyErr::from) + } else { + // SAFETY: the caller promised they have upheld the coherence guarantees. + Ok(unsafe { Self::new_unchecked(num_qubits, coeffs, bit_terms, indices, boundaries) }) + } + } +} + +/// A view object onto a single term of a `SparseObservable`. +/// +/// The lengths of `bit_terms` and `indices` are guaranteed to be created equal, but might be zero +/// (in the case that the term is proportional to the identity). +#[derive(Clone, Copy, Debug)] +pub struct SparseTermView<'a> { + pub coeff: Complex64, + pub bit_terms: &'a [BitTerm], + pub indices: &'a [u32], +} + +/// Helper class of `ArrayView` that denotes the slot of the `SparseObservable` we're looking at. +#[derive(Clone, Copy, PartialEq, Eq)] +enum ArraySlot { + Coeffs, + BitTerms, + Indices, + Boundaries, +} + +/// Custom wrapper sequence class to get safe views onto the Rust-space data. We can't directly +/// expose Python-managed wrapped pointers without introducing some form of runtime exclusion on the +/// ability of `SparseObservable` to re-allocate in place; we can't leave dangling pointers for +/// Python space. +#[pyclass(frozen, sequence)] +struct ArrayView { + base: Py, + slot: ArraySlot, +} +#[pymethods] +impl ArrayView { + fn __repr__(&self, py: Python) -> PyResult { + let obs = self.base.borrow(py); + let data = match self.slot { + // Simple integers look the same in Rust-space debug as Python. + ArraySlot::Indices => format!("{:?}", obs.indices), + ArraySlot::Boundaries => format!("{:?}", obs.boundaries), + // Complexes don't have a nice repr in Rust, so just delegate the whole load to Python + // and convert back. + ArraySlot::Coeffs => PyList::new_bound(py, &obs.coeffs).repr()?.to_string(), + ArraySlot::BitTerms => format!( + "[{}]", + obs.bit_terms + .iter() + .map(BitTerm::py_label) + .collect::>() + .join(", ") + ), + }; + Ok(format!( + "", + match self.slot { + ArraySlot::Coeffs => "coeffs", + ArraySlot::BitTerms => "bit_terms", + ArraySlot::Indices => "indices", + ArraySlot::Boundaries => "boundaries", + }, + data, + )) + } + + fn __getitem__(&self, py: Python, index: PySequenceIndex) -> PyResult> { + fn get_from_slice( + py: Python, + slice: &[T], + index: PySequenceIndex, + ) -> PyResult> { + match index.with_len(slice.len())? { + SequenceIndex::Int(index) => Ok(slice[index].to_object(py)), + indices => Ok(PyArray1::from_iter_bound( + py, + indices.iter().map(|index| slice[index]), + ) + .into_any() + .unbind()), + } + } + + let obs = self.base.borrow(py); + match self.slot { + ArraySlot::Coeffs => get_from_slice(py, &obs.coeffs, index), + ArraySlot::BitTerms => get_from_slice( + py, + ::bytemuck::cast_slice::(&obs.bit_terms), + index, + ), + ArraySlot::Indices => get_from_slice(py, &obs.indices, index), + ArraySlot::Boundaries => get_from_slice(py, &obs.boundaries, index), + } + } + + fn __setitem__(&self, index: PySequenceIndex, values: &Bound) -> PyResult<()> { + /// Set values of a slice according to the indexer, using `extract` to retrieve the + /// Rust-space object from the collection of Python-space values. + /// + /// This indirects the Python extraction through an intermediate type to marginally improve + /// the error messages for things like `BitTerm`, where Python-space extraction might fail + /// because the user supplied an invalid alphabet letter. + /// + /// This allows broadcasting a single item into many locations in a slice (like Numpy), but + /// otherwise requires that the index and values are the same length (unlike Python's + /// `list`) because that would change the length. + fn set_in_slice<'py, T, S>( + slice: &mut [T], + index: PySequenceIndex<'py>, + values: &Bound<'py, PyAny>, + ) -> PyResult<()> + where + T: Copy + TryFrom, + S: FromPyObject<'py>, + PyErr: From<>::Error>, + { + match index.with_len(slice.len())? { + SequenceIndex::Int(index) => { + slice[index] = values.extract::()?.try_into()?; + Ok(()) + } + indices => { + if let Ok(value) = values.extract::() { + let value = value.try_into()?; + for index in indices { + slice[index] = value; + } + } else { + let values = values + .iter()? + .map(|value| value?.extract::()?.try_into().map_err(PyErr::from)) + .collect::>>()?; + if indices.len() != values.len() { + return Err(PyValueError::new_err(format!( + "tried to set a slice of length {} with a sequence of length {}", + indices.len(), + values.len(), + ))); + } + for (index, value) in indices.into_iter().zip(values) { + slice[index] = value; + } + } + Ok(()) + } + } + } + + let mut obs = self.base.borrow_mut(values.py()); + match self.slot { + ArraySlot::Coeffs => set_in_slice::<_, Complex64>(&mut obs.coeffs, index, values), + ArraySlot::BitTerms => set_in_slice::(&mut obs.bit_terms, index, values), + ArraySlot::Indices => set_in_slice::<_, u32>(&mut obs.indices, index, values), + ArraySlot::Boundaries => set_in_slice::<_, usize>(&mut obs.boundaries, index, values), + } + } + + fn __len__(&self, py: Python) -> usize { + let obs = self.base.borrow(py); + match self.slot { + ArraySlot::Coeffs => obs.coeffs.len(), + ArraySlot::BitTerms => obs.bit_terms.len(), + ArraySlot::Indices => obs.indices.len(), + ArraySlot::Boundaries => obs.boundaries.len(), + } + } + + #[pyo3(signature = (/, dtype=None, copy=None))] + fn __array__<'py>( + &self, + py: Python<'py>, + dtype: Option<&Bound<'py, PyAny>>, + copy: Option, + ) -> PyResult> { + // This method always copies, so we don't leave dangling pointers lying around in Numpy + // arrays; it's not enough just to set the `base` of the Numpy array to the + // `SparseObservable`, since the `Vec` we're referring to might re-allocate and invalidate + // the pointer the Numpy array is wrapping. + if !copy.unwrap_or(true) { + return Err(PyValueError::new_err( + "cannot produce a safe view onto movable memory", + )); + } + let obs = self.base.borrow(py); + match self.slot { + ArraySlot::Coeffs => { + cast_array_type(py, PyArray1::from_slice_bound(py, &obs.coeffs), dtype) + } + ArraySlot::Indices => { + cast_array_type(py, PyArray1::from_slice_bound(py, &obs.indices), dtype) + } + ArraySlot::Boundaries => { + cast_array_type(py, PyArray1::from_slice_bound(py, &obs.boundaries), dtype) + } + ArraySlot::BitTerms => { + let bit_terms: &[u8] = ::bytemuck::cast_slice(&obs.bit_terms); + cast_array_type(py, PyArray1::from_slice_bound(py, bit_terms), dtype) + } + } + } +} + +/// Use the Numpy Python API to convert a `PyArray` into a dynamically chosen `dtype`, copying only +/// if required. +fn cast_array_type<'py, T>( + py: Python<'py>, + array: Bound<'py, PyArray1>, + dtype: Option<&Bound<'py, PyAny>>, +) -> PyResult> { + let base_dtype = array.dtype(); + let dtype = dtype + .map(|dtype| PyArrayDescr::new_bound(py, dtype)) + .unwrap_or_else(|| Ok(base_dtype.clone()))?; + if dtype.is_equiv_to(&base_dtype) { + return Ok(array.into_any()); + } + PyModule::import_bound(py, intern!(py, "numpy"))? + .getattr(intern!(py, "array"))? + .call( + (array,), + Some( + &[ + (intern!(py, "copy"), NUMPY_COPY_ONLY_IF_NEEDED.get_bound(py)), + (intern!(py, "dtype"), dtype.as_any()), + ] + .into_py_dict_bound(py), + ), + ) + .map(|obj| obj.into_any()) +} + +pub fn sparse_observable(m: &Bound) -> PyResult<()> { + m.add_class::()?; + Ok(()) +} + +#[cfg(test)] +mod test { + use super::*; + + #[test] + fn test_error_safe_add_dense_label() { + let base = SparseObservable::identity(5); + let mut modified = base.clone(); + assert!(matches!( + modified.add_dense_label("IIZ$X", Complex64::new(1.0, 0.0)), + Err(LabelError::OutsideAlphabet) + )); + // `modified` should have been left in a valid state. + assert_eq!(base, modified); + } +} diff --git a/crates/circuit/src/imports.rs b/crates/circuit/src/imports.rs index 6e87bcb4f101..201aaa4e0ad7 100644 --- a/crates/circuit/src/imports.rs +++ b/crates/circuit/src/imports.rs @@ -116,6 +116,8 @@ pub static UNITARY_GATE: ImportOnceCell = ImportOnceCell::new( "qiskit.circuit.library.generalized_gates.unitary", "UnitaryGate", ); +pub static NUMPY_COPY_ONLY_IF_NEEDED: ImportOnceCell = + ImportOnceCell::new("qiskit._numpy_compat", "COPY_ONLY_IF_NEEDED"); /// A mapping from the enum variant in crate::operations::StandardGate to the python /// module path and class name to import it. This is used to populate the conversion table diff --git a/crates/pyext/src/lib.rs b/crates/pyext/src/lib.rs index 2e0929543f1b..7448ca5de1f1 100644 --- a/crates/pyext/src/lib.rs +++ b/crates/pyext/src/lib.rs @@ -52,6 +52,7 @@ fn _accelerate(m: &Bound) -> PyResult<()> { add_submodule(m, ::qiskit_accelerate::results::results, "results")?; add_submodule(m, ::qiskit_accelerate::sabre::sabre, "sabre")?; add_submodule(m, ::qiskit_accelerate::sampled_exp_val::sampled_exp_val, "sampled_exp_val")?; + add_submodule(m, ::qiskit_accelerate::sparse_observable::sparse_observable, "sparse_observable")?; add_submodule(m, ::qiskit_accelerate::sparse_pauli_op::sparse_pauli_op, "sparse_pauli_op")?; add_submodule(m, ::qiskit_accelerate::split_2q_unitaries::split_2q_unitaries_mod, "split_2q_unitaries")?; add_submodule(m, ::qiskit_accelerate::star_prerouting::star_prerouting, "star_prerouting")?; diff --git a/qiskit/__init__.py b/qiskit/__init__.py index 1566db0be2c2..5e954d65feb6 100644 --- a/qiskit/__init__.py +++ b/qiskit/__init__.py @@ -77,6 +77,7 @@ sys.modules["qiskit._accelerate.results"] = _accelerate.results sys.modules["qiskit._accelerate.sabre"] = _accelerate.sabre sys.modules["qiskit._accelerate.sampled_exp_val"] = _accelerate.sampled_exp_val +sys.modules["qiskit._accelerate.sparse_observable"] = _accelerate.sparse_observable sys.modules["qiskit._accelerate.sparse_pauli_op"] = _accelerate.sparse_pauli_op sys.modules["qiskit._accelerate.star_prerouting"] = _accelerate.star_prerouting sys.modules["qiskit._accelerate.stochastic_swap"] = _accelerate.stochastic_swap diff --git a/qiskit/quantum_info/__init__.py b/qiskit/quantum_info/__init__.py index ba3299c9e100..2f5d707c8802 100644 --- a/qiskit/quantum_info/__init__.py +++ b/qiskit/quantum_info/__init__.py @@ -28,6 +28,7 @@ Pauli Clifford ScalarOp + SparseObservable SparsePauliOp CNOTDihedral PauliList @@ -113,6 +114,9 @@ """ from __future__ import annotations + +from qiskit._accelerate.sparse_observable import SparseObservable + from .analysis import hellinger_distance, hellinger_fidelity, Z2Symmetries from .operators import ( Clifford, diff --git a/releasenotes/notes/sparse-observable-7de70dcdf6962a64.yaml b/releasenotes/notes/sparse-observable-7de70dcdf6962a64.yaml new file mode 100644 index 000000000000..48617e6d1968 --- /dev/null +++ b/releasenotes/notes/sparse-observable-7de70dcdf6962a64.yaml @@ -0,0 +1,32 @@ +--- +features_quantum_info: + - | + A new observable class has been added. :class:`.SparseObservable` represents observables as a + sum of terms, similar to :class:`.SparsePauliOp`, but with two core differences: + + 1. Each complete term is stored as (effectively) a series of ``(qubit, bit_term)`` pairs, + without storing qubits that undergo the identity for that term. This significantly improves + the memory usage of observables such as the weighted sum of Paulis :math:`\sum_i c_i Z_i`. + + 2. The single-qubit term alphabet is overcomplete for the operator space; it can represent Pauli + operators (like :class:`.SparsePauliOp`), but also projectors onto the eigenstates of the + Pauli operators, like :math:`\lvert 0\rangle\langle 0\rangle`. Such projectors can be + measured on hardware equally as efficiently as their corresponding Pauli operator, but + :class:`.SparsePauliOp` would require an exponential number of terms to represent + :math:`{\lvert0\rangle\langle0\rvert}^{\otimes n}` over :math:`n` qubits, while + :class:`.SparseObservable` needs only a single term. + + You can construct and manipulate :class:`.SparseObservable` using an interface familiar to users + of :class:`.SparsePauliOp`:: + + from qiskit.quantum_info import SparseObservable + + obs = SparseObservable.from_sparse_list([ + ("XZY", (2, 1, 0), 1.5j), + ("+-", (100, 99), 0.5j), + ("01", (50, 49), 0.5), + ]) + + :class:`.SparseObservable` is not currently supported as an input format to the primitives + (:mod:`qiskit.primitives`), but we expect to expand these interfaces to include them in the + future. diff --git a/test/python/quantum_info/test_sparse_observable.py b/test/python/quantum_info/test_sparse_observable.py new file mode 100644 index 000000000000..b8510282b943 --- /dev/null +++ b/test/python/quantum_info/test_sparse_observable.py @@ -0,0 +1,924 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2024 +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +# pylint: disable=missing-module-docstring,missing-class-docstring,missing-function-docstring + +import copy +import pickle +import unittest + +import ddt +import numpy as np + +from qiskit.circuit import Parameter +from qiskit.quantum_info import SparseObservable, SparsePauliOp, Pauli + +from test import QiskitTestCase # pylint: disable=wrong-import-order + + +@ddt.ddt +class TestSparseObservable(QiskitTestCase): + def test_default_constructor_pauli(self): + data = Pauli("IXYIZ") + self.assertEqual(SparseObservable(data), SparseObservable.from_pauli(data)) + self.assertEqual( + SparseObservable(data, num_qubits=data.num_qubits), SparseObservable.from_pauli(data) + ) + with self.assertRaisesRegex(ValueError, "explicitly given 'num_qubits'"): + SparseObservable(data, num_qubits=data.num_qubits + 1) + + with_phase = Pauli("-jIYYXY") + self.assertEqual(SparseObservable(with_phase), SparseObservable.from_pauli(with_phase)) + self.assertEqual( + SparseObservable(with_phase, num_qubits=data.num_qubits), + SparseObservable.from_pauli(with_phase), + ) + + self.assertEqual(SparseObservable(Pauli("")), SparseObservable.from_pauli(Pauli(""))) + + def test_default_constructor_sparse_pauli_op(self): + data = SparsePauliOp.from_list([("IIXIY", 1.0), ("XYYZI", -0.25), ("XYIYY", -0.25 + 0.75j)]) + self.assertEqual(SparseObservable(data), SparseObservable.from_sparse_pauli_op(data)) + self.assertEqual( + SparseObservable(data, num_qubits=data.num_qubits), + SparseObservable.from_sparse_pauli_op(data), + ) + with self.assertRaisesRegex(ValueError, "explicitly given 'num_qubits'"): + SparseObservable(data, num_qubits=data.num_qubits + 1) + with self.assertRaisesRegex(TypeError, "complex-typed coefficients"): + SparseObservable(SparsePauliOp(["XX"], [Parameter("x")])) + + def test_default_constructor_label(self): + data = "IIXI+-I01rlIYZ" + self.assertEqual(SparseObservable(data), SparseObservable.from_label(data)) + self.assertEqual( + SparseObservable(data, num_qubits=len(data)), SparseObservable.from_label(data) + ) + with self.assertRaisesRegex(ValueError, "explicitly given 'num_qubits'"): + SparseObservable(data, num_qubits=len(data) + 1) + + def test_default_constructor_list(self): + data = [("IXIIZ", 0.5), ("+I-II", 1.0 - 0.25j), ("IIrlI", -0.75)] + self.assertEqual(SparseObservable(data), SparseObservable.from_list(data)) + self.assertEqual(SparseObservable(data, num_qubits=5), SparseObservable.from_list(data)) + with self.assertRaisesRegex(ValueError, "label with length 5 cannot be added"): + SparseObservable(data, num_qubits=4) + with self.assertRaisesRegex(ValueError, "label with length 5 cannot be added"): + SparseObservable(data, num_qubits=6) + self.assertEqual( + SparseObservable([], num_qubits=5), SparseObservable.from_list([], num_qubits=5) + ) + + def test_default_constructor_sparse_list(self): + data = [("ZX", (0, 3), 0.5), ("-+", (2, 4), 1.0 - 0.25j), ("rl", (2, 1), -0.75)] + self.assertEqual( + SparseObservable(data, num_qubits=5), + SparseObservable.from_sparse_list(data, num_qubits=5), + ) + self.assertEqual( + SparseObservable(data, num_qubits=10), + SparseObservable.from_sparse_list(data, num_qubits=10), + ) + with self.assertRaisesRegex(ValueError, "'num_qubits' must be provided"): + SparseObservable(data) + self.assertEqual( + SparseObservable([], num_qubits=5), SparseObservable.from_sparse_list([], num_qubits=5) + ) + + def test_default_constructor_copy(self): + base = SparseObservable.from_list([("IXIZIY", 1.0), ("+-rl01", -1.0j)]) + copied = SparseObservable(base) + self.assertEqual(base, copied) + self.assertIsNot(base, copied) + + # Modifications to `copied` don't propagate back. + copied.coeffs[1] = -0.5j + self.assertNotEqual(base, copied) + + with self.assertRaisesRegex(ValueError, "explicitly given 'num_qubits'"): + SparseObservable(base, num_qubits=base.num_qubits + 1) + + def test_default_constructor_failed_inference(self): + with self.assertRaises(TypeError): + # Mixed dense/sparse list. + SparseObservable([("IIXIZ", 1.0), ("+-", (2, 3), -1.0)], num_qubits=5) + + def test_num_qubits(self): + self.assertEqual(SparseObservable.zero(0).num_qubits, 0) + self.assertEqual(SparseObservable.zero(10).num_qubits, 10) + + self.assertEqual(SparseObservable.identity(0).num_qubits, 0) + self.assertEqual(SparseObservable.identity(1_000_000).num_qubits, 1_000_000) + + def test_num_ops(self): + self.assertEqual(SparseObservable.zero(0).num_ops, 0) + self.assertEqual(SparseObservable.zero(10).num_ops, 0) + self.assertEqual(SparseObservable.identity(0).num_ops, 1) + self.assertEqual(SparseObservable.identity(1_000_000).num_ops, 1) + self.assertEqual(SparseObservable.from_list([("IIIXIZ", 1.0), ("YY+-II", 0.5j)]).num_ops, 2) + + def test_bit_term_enum(self): + # These are very explicit tests that effectively just duplicate magic numbers, but the point + # is that those magic numbers are required to be constant as their values are part of the + # public interface. + + self.assertEqual( + set(SparseObservable.BitTerm), + { + SparseObservable.BitTerm.X, + SparseObservable.BitTerm.Y, + SparseObservable.BitTerm.Z, + SparseObservable.BitTerm.PLUS, + SparseObservable.BitTerm.MINUS, + SparseObservable.BitTerm.RIGHT, + SparseObservable.BitTerm.LEFT, + SparseObservable.BitTerm.ZERO, + SparseObservable.BitTerm.ONE, + }, + ) + # All the enumeration items should also be integers. + self.assertIsInstance(SparseObservable.BitTerm.X, int) + values = { + "X": 0b00_10, + "Y": 0b00_11, + "Z": 0b00_01, + "PLUS": 0b10_10, + "MINUS": 0b01_10, + "RIGHT": 0b10_11, + "LEFT": 0b01_11, + "ZERO": 0b10_01, + "ONE": 0b01_01, + } + self.assertEqual({name: getattr(SparseObservable.BitTerm, name) for name in values}, values) + + # The single-character label aliases can be accessed with index notation. + labels = { + "X": SparseObservable.BitTerm.X, + "Y": SparseObservable.BitTerm.Y, + "Z": SparseObservable.BitTerm.Z, + "+": SparseObservable.BitTerm.PLUS, + "-": SparseObservable.BitTerm.MINUS, + "r": SparseObservable.BitTerm.RIGHT, + "l": SparseObservable.BitTerm.LEFT, + "0": SparseObservable.BitTerm.ZERO, + "1": SparseObservable.BitTerm.ONE, + } + self.assertEqual({label: SparseObservable.BitTerm[label] for label in labels}, labels) + + @ddt.data( + SparseObservable.zero(0), + SparseObservable.zero(10), + SparseObservable.identity(0), + SparseObservable.identity(1_000), + SparseObservable.from_label("IIXIZI"), + SparseObservable.from_list([("YIXZII", -0.25), ("01rl+-", 0.25 + 0.5j)]), + ) + def test_pickle(self, observable): + self.assertEqual(observable, copy.copy(observable)) + self.assertIsNot(observable, copy.copy(observable)) + self.assertEqual(observable, copy.deepcopy(observable)) + self.assertEqual(observable, pickle.loads(pickle.dumps(observable))) + + @ddt.data( + # This is every combination of (0, 1, many) for (terms, qubits, non-identites per term). + SparseObservable.zero(0), + SparseObservable.zero(1), + SparseObservable.zero(10), + SparseObservable.identity(0), + SparseObservable.identity(1), + SparseObservable.identity(1_000), + SparseObservable.from_label("IIXIZI"), + SparseObservable.from_label("X"), + SparseObservable.from_list([("YIXZII", -0.25), ("01rl+-", 0.25 + 0.5j)]), + ) + def test_repr(self, data): + # The purpose of this is just to test that the `repr` doesn't crash, rather than asserting + # that it has any particular form. + self.assertIsInstance(repr(data), str) + self.assertIn("SparseObservable", repr(data)) + + def test_from_raw_parts(self): + # Happiest path: exactly typed inputs. + num_qubits = 100 + terms = np.full((num_qubits,), SparseObservable.BitTerm.Z, dtype=np.uint8) + indices = np.arange(num_qubits, dtype=np.uint32) + coeffs = np.ones((num_qubits,), dtype=complex) + boundaries = np.arange(num_qubits + 1, dtype=np.uintp) + observable = SparseObservable.from_raw_parts(num_qubits, coeffs, terms, indices, boundaries) + self.assertEqual(observable.num_qubits, num_qubits) + np.testing.assert_equal(observable.bit_terms, terms) + np.testing.assert_equal(observable.indices, indices) + np.testing.assert_equal(observable.coeffs, coeffs) + np.testing.assert_equal(observable.boundaries, boundaries) + + self.assertEqual( + observable, + SparseObservable.from_raw_parts( + num_qubits, coeffs, terms, indices, boundaries, check=False + ), + ) + + # At least the initial implementation of `SparseObservable` requires `from_raw_parts` to be + # a copy constructor in order to allow it to be resized by Rust space. This is checking for + # that, but if the implementation changes, it could potentially be relaxed. + self.assertFalse(np.may_share_memory(observable.coeffs, coeffs)) + + # Conversion from array-likes, including mis-typed but compatible arrays. + observable = SparseObservable.from_raw_parts( + num_qubits, list(coeffs), tuple(terms), observable.indices, boundaries.astype(np.int16) + ) + self.assertEqual(observable.num_qubits, num_qubits) + np.testing.assert_equal(observable.bit_terms, terms) + np.testing.assert_equal(observable.indices, indices) + np.testing.assert_equal(observable.coeffs, coeffs) + np.testing.assert_equal(observable.boundaries, boundaries) + + # Construction of zero operator. + self.assertEqual( + SparseObservable.from_raw_parts(10, [], [], [], [0]), SparseObservable.zero(10) + ) + + # Construction of an operator with an intermediate identity term. For the initial + # constructor tests, it's hard to check anything more than the construction succeeded. + self.assertEqual( + SparseObservable.from_raw_parts( + 10, [1.0j, 0.5, 2.0], [1, 3, 2], [0, 1, 2], [0, 1, 1, 3] + ).num_ops, + # The three are [(1.0j)*(Z_1), 0.5, 2.0*(X_2 Y_1)] + 3, + ) + + def test_from_raw_parts_checks_coherence(self): + with self.assertRaisesRegex(ValueError, "not a valid letter"): + SparseObservable.from_raw_parts(2, [1.0j], [ord("$")], [0], [0, 1]) + with self.assertRaisesRegex(ValueError, r"boundaries.*must be one element longer"): + SparseObservable.from_raw_parts(2, [1.0j], [SparseObservable.BitTerm.Z], [0], [0]) + with self.assertRaisesRegex(ValueError, r"`bit_terms` \(1\) and `indices` \(0\)"): + SparseObservable.from_raw_parts(2, [1.0j], [SparseObservable.BitTerm.Z], [], [0, 1]) + with self.assertRaisesRegex(ValueError, r"`bit_terms` \(0\) and `indices` \(1\)"): + SparseObservable.from_raw_parts(2, [1.0j], [], [1], [0, 1]) + with self.assertRaisesRegex(ValueError, r"the first item of `boundaries` \(1\) must be 0"): + SparseObservable.from_raw_parts(2, [1.0j], [SparseObservable.BitTerm.Z], [0], [1, 1]) + with self.assertRaisesRegex(ValueError, r"the last item of `boundaries` \(2\)"): + SparseObservable.from_raw_parts(2, [1.0j], [1], [0], [0, 2]) + with self.assertRaisesRegex(ValueError, r"the last item of `boundaries` \(1\)"): + SparseObservable.from_raw_parts(2, [1.0j], [1, 2], [0, 1], [0, 1]) + with self.assertRaisesRegex(ValueError, r"all qubit indices must be less than the number"): + SparseObservable.from_raw_parts(4, [1.0j], [1, 2], [0, 4], [0, 2]) + with self.assertRaisesRegex(ValueError, r"all qubit indices must be less than the number"): + SparseObservable.from_raw_parts(4, [1.0j, -0.5j], [1, 2], [0, 4], [0, 1, 2]) + with self.assertRaisesRegex(ValueError, "the values in `boundaries` include backwards"): + SparseObservable.from_raw_parts( + 5, [1.0j, -0.5j, 2.0], [1, 2, 3, 2], [0, 1, 2, 3], [0, 2, 1, 4] + ) + with self.assertRaisesRegex( + ValueError, "the values in `indices` are not term-wise increasing" + ): + SparseObservable.from_raw_parts(4, [1.0j], [1, 2], [1, 0], [0, 2]) + + # There's no test of attempting to pass incoherent data and `check=False` because that + # permits undefined behaviour in Rust (it's unsafe), so all bets would be off. + + def test_from_label(self): + # The label is interpreted like a bitstring, with the right-most item associated with qubit + # 0, and increasing as we move to the left (like `Pauli`, and other bitstring conventions). + self.assertEqual( + # Ruler for counting terms: dcba9876543210 + SparseObservable.from_label("I+-II01IrlIXYZ"), + SparseObservable.from_raw_parts( + 14, + [1.0], + [ + SparseObservable.BitTerm.Z, + SparseObservable.BitTerm.Y, + SparseObservable.BitTerm.X, + SparseObservable.BitTerm.LEFT, + SparseObservable.BitTerm.RIGHT, + SparseObservable.BitTerm.ONE, + SparseObservable.BitTerm.ZERO, + SparseObservable.BitTerm.MINUS, + SparseObservable.BitTerm.PLUS, + ], + [0, 1, 2, 4, 5, 7, 8, 11, 12], + [0, 9], + ), + ) + + self.assertEqual(SparseObservable.from_label("I" * 10), SparseObservable.identity(10)) + + # The empty label case is a 0-qubit identity, since `from_label` always sets a coefficient + # of 1. + self.assertEqual(SparseObservable.from_label(""), SparseObservable.identity(0)) + + def test_from_label_failures(self): + with self.assertRaisesRegex(ValueError, "labels must only contain letters from"): + # Bad letters that are still ASCII. + SparseObservable.from_label("I+-$%I") + with self.assertRaisesRegex(ValueError, "labels must only contain letters from"): + # Unicode shenangigans. + SparseObservable.from_label("🐍") + + def test_from_list(self): + label = "IXYI+-0lr1ZZY" + self.assertEqual( + SparseObservable.from_list([(label, 1.0)]), SparseObservable.from_label(label) + ) + self.assertEqual( + SparseObservable.from_list([(label, 1.0)], num_qubits=len(label)), + SparseObservable.from_label(label), + ) + + self.assertEqual( + SparseObservable.from_list([("IIIXZI", 1.0j), ("+-IIII", -0.5)]), + SparseObservable.from_raw_parts( + 6, + [1.0j, -0.5], + [ + SparseObservable.BitTerm.Z, + SparseObservable.BitTerm.X, + SparseObservable.BitTerm.MINUS, + SparseObservable.BitTerm.PLUS, + ], + [1, 2, 4, 5], + [0, 2, 4], + ), + ) + + self.assertEqual(SparseObservable.from_list([], num_qubits=5), SparseObservable.zero(5)) + self.assertEqual(SparseObservable.from_list([], num_qubits=0), SparseObservable.zero(0)) + + def test_from_list_failures(self): + with self.assertRaisesRegex(ValueError, "labels must only contain letters from"): + # Bad letters that are still ASCII. + SparseObservable.from_list([("XZIIZY", 0.5), ("I+-$%I", 1.0j)]) + with self.assertRaisesRegex(ValueError, "labels must only contain letters from"): + # Unicode shenangigans. + SparseObservable.from_list([("🐍", 0.5)]) + with self.assertRaisesRegex(ValueError, "label with length 4 cannot be added"): + SparseObservable.from_list([("IIZ", 0.5), ("IIXI", 1.0j)]) + with self.assertRaisesRegex(ValueError, "label with length 2 cannot be added"): + SparseObservable.from_list([("IIZ", 0.5), ("II", 1.0j)]) + with self.assertRaisesRegex(ValueError, "label with length 3 cannot be added"): + SparseObservable.from_list([("IIZ", 0.5), ("IXI", 1.0j)], num_qubits=2) + with self.assertRaisesRegex(ValueError, "label with length 3 cannot be added"): + SparseObservable.from_list([("IIZ", 0.5), ("IXI", 1.0j)], num_qubits=4) + with self.assertRaisesRegex(ValueError, "cannot construct.*without knowing `num_qubits`"): + SparseObservable.from_list([]) + + def test_from_sparse_list(self): + self.assertEqual( + SparseObservable.from_sparse_list( + [ + ("XY", (0, 1), 0.5), + ("+-", (1, 3), -0.25j), + ("rl0", (0, 2, 4), 1.0j), + ], + num_qubits=5, + ), + SparseObservable.from_list([("IIIYX", 0.5), ("I-I+I", -0.25j), ("0IlIr", 1.0j)]), + ) + + # The indices should be allowed to be given in unsorted order, but they should be term-wise + # sorted in the output. + from_unsorted = SparseObservable.from_sparse_list( + [ + ("XYZ", (2, 1, 0), 1.5j), + ("+rl", (2, 0, 1), -0.5), + ], + num_qubits=3, + ) + self.assertEqual(from_unsorted, SparseObservable.from_list([("XYZ", 1.5j), ("+lr", -0.5)])) + np.testing.assert_equal( + from_unsorted.indices, np.array([0, 1, 2, 0, 1, 2], dtype=np.uint32) + ) + + # Explicit identities should still work, just be skipped over. + explicit_identity = SparseObservable.from_sparse_list( + [ + ("ZXI", (0, 1, 2), 1.0j), + ("XYIII", (0, 1, 2, 3, 8), -0.5j), + ], + num_qubits=10, + ) + self.assertEqual( + explicit_identity, + SparseObservable.from_sparse_list( + [("XZ", (1, 0), 1.0j), ("YX", (1, 0), -0.5j)], num_qubits=10 + ), + ) + np.testing.assert_equal(explicit_identity.indices, np.array([0, 1, 0, 1], dtype=np.uint32)) + + self.assertEqual( + SparseObservable.from_sparse_list([("", (), 1.0)], num_qubits=5), + SparseObservable.identity(5), + ) + self.assertEqual( + SparseObservable.from_sparse_list([("", (), 1.0)], num_qubits=0), + SparseObservable.identity(0), + ) + + self.assertEqual( + SparseObservable.from_sparse_list([], num_qubits=1_000_000), + SparseObservable.zero(1_000_000), + ) + self.assertEqual( + SparseObservable.from_sparse_list([], num_qubits=0), + SparseObservable.zero(0), + ) + + def test_from_sparse_list_failures(self): + with self.assertRaisesRegex(ValueError, "labels must only contain letters from"): + # Bad letters that are still ASCII. + SparseObservable.from_sparse_list( + [("XZZY", (5, 3, 1, 0), 0.5), ("+$", (2, 1), 1.0j)], num_qubits=8 + ) + # Unicode shenangigans. These two should fail with a `ValueError`, but the exact message + # isn't important. "\xff" is "ÿ", which is two bytes in UTF-8 (so has a length of 2 in + # Rust), but has a length of 1 in Python, so try with both a length-1 and length-2 index + # sequence, and both should still raise `ValueError`. + with self.assertRaises(ValueError): + SparseObservable.from_sparse_list([("\xff", (1,), 0.5)], num_qubits=5) + with self.assertRaises(ValueError): + SparseObservable.from_sparse_list([("\xff", (1, 2), 0.5)], num_qubits=5) + + with self.assertRaisesRegex(ValueError, "label with length 2 does not match indices"): + SparseObservable.from_sparse_list([("XZ", (0,), 1.0)], num_qubits=5) + with self.assertRaisesRegex(ValueError, "label with length 2 does not match indices"): + SparseObservable.from_sparse_list([("XZ", (0, 1, 2), 1.0)], num_qubits=5) + + with self.assertRaisesRegex(ValueError, "index 3 is out of range for a 3-qubit operator"): + SparseObservable.from_sparse_list([("XZY", (0, 1, 3), 1.0)], num_qubits=3) + with self.assertRaisesRegex(ValueError, "index 4 is out of range for a 3-qubit operator"): + SparseObservable.from_sparse_list([("XZY", (0, 1, 4), 1.0)], num_qubits=3) + with self.assertRaisesRegex(ValueError, "index 3 is out of range for a 3-qubit operator"): + # ... even if it's for an explicit identity. + SparseObservable.from_sparse_list([("+-I", (0, 1, 3), 1.0)], num_qubits=3) + + with self.assertRaisesRegex(ValueError, "index 3 is duplicated"): + SparseObservable.from_sparse_list([("XZ", (3, 3), 1.0)], num_qubits=5) + with self.assertRaisesRegex(ValueError, "index 3 is duplicated"): + SparseObservable.from_sparse_list([("XYZXZ", (3, 0, 1, 2, 3), 1.0)], num_qubits=5) + + def test_from_pauli(self): + # This function should be infallible provided `Pauli` doesn't change its interface and the + # user doesn't violate the typing. + + # Simple check that the labels are interpreted in the same order. + self.assertEqual( + SparseObservable.from_pauli(Pauli("IIXZI")), SparseObservable.from_label("IIXZI") + ) + + # `Pauli` accepts a phase in its label, which we can't (because of clashes with the other + # alphabet letters), and we should get that right. + self.assertEqual( + SparseObservable.from_pauli(Pauli("iIXZIX")), + SparseObservable.from_list([("IXZIX", 1.0j)]), + ) + self.assertEqual( + SparseObservable.from_pauli(Pauli("-iIXZIX")), + SparseObservable.from_list([("IXZIX", -1.0j)]), + ) + self.assertEqual( + SparseObservable.from_pauli(Pauli("-IXZIX")), + SparseObservable.from_list([("IXZIX", -1.0)]), + ) + + # `Pauli` has its internal phase convention for how it stores `Y`; we should get this right + # regardless of how many Ys are in the label, or if there's a phase. + paulis = {"IXYZ" * n: Pauli("IXYZ" * n) for n in range(1, 5)} + from_paulis, from_labels = zip( + *( + (SparseObservable.from_pauli(pauli), SparseObservable.from_label(label)) + for label, pauli in paulis.items() + ) + ) + self.assertEqual(from_paulis, from_labels) + + phased_paulis = {"IXYZ" * n: Pauli("j" + "IXYZ" * n) for n in range(1, 5)} + from_paulis, from_lists = zip( + *( + (SparseObservable.from_pauli(pauli), SparseObservable.from_list([(label, 1.0j)])) + for label, pauli in phased_paulis.items() + ) + ) + self.assertEqual(from_paulis, from_lists) + + self.assertEqual(SparseObservable.from_pauli(Pauli("III")), SparseObservable.identity(3)) + self.assertEqual(SparseObservable.from_pauli(Pauli("")), SparseObservable.identity(0)) + + def test_from_sparse_pauli_op(self): + self.assertEqual( + SparseObservable.from_sparse_pauli_op(SparsePauliOp.from_list([("IIIII", 1.0)])), + SparseObservable.identity(5), + ) + + data = [("ZXZXZ", 0.25), ("IYXZI", 1.0j), ("IYYZX", 0.5), ("YYYXI", -0.5), ("IYYYY", 2j)] + self.assertEqual( + SparseObservable.from_sparse_pauli_op(SparsePauliOp.from_list(data)), + SparseObservable.from_list(data), + ) + + # These two _should_ produce the same structure as `SparseObservable.zero(num_qubits)`, but + # because `SparsePauliOp` doesn't represent the zero operator "natively" - with an empty sum + # - they actually come out looking like `0.0` times the identity, which is less efficient + # but acceptable. + self.assertEqual( + SparseObservable.from_sparse_pauli_op(SparsePauliOp.from_list([], num_qubits=1)), + SparseObservable.from_list([("I", 0.0)]), + ) + self.assertEqual( + SparseObservable.from_sparse_pauli_op(SparsePauliOp.from_list([], num_qubits=0)), + SparseObservable.from_list([("", 0.0)]), + ) + + def test_from_sparse_pauli_op_failures(self): + parametric = SparsePauliOp.from_list([("IIXZ", Parameter("x"))], dtype=object) + with self.assertRaisesRegex(TypeError, "complex-typed coefficients"): + SparseObservable.from_sparse_pauli_op(parametric) + + def test_zero(self): + zero_5 = SparseObservable.zero(5) + self.assertEqual(zero_5.num_qubits, 5) + np.testing.assert_equal(zero_5.coeffs, np.array([], dtype=complex)) + np.testing.assert_equal(zero_5.bit_terms, np.array([], dtype=np.uint8)) + np.testing.assert_equal(zero_5.indices, np.array([], dtype=np.uint32)) + np.testing.assert_equal(zero_5.boundaries, np.array([0], dtype=np.uintp)) + + zero_0 = SparseObservable.zero(0) + self.assertEqual(zero_0.num_qubits, 0) + np.testing.assert_equal(zero_0.coeffs, np.array([], dtype=complex)) + np.testing.assert_equal(zero_0.bit_terms, np.array([], dtype=np.uint8)) + np.testing.assert_equal(zero_0.indices, np.array([], dtype=np.uint32)) + np.testing.assert_equal(zero_0.boundaries, np.array([0], dtype=np.uintp)) + + def test_identity(self): + id_5 = SparseObservable.identity(5) + self.assertEqual(id_5.num_qubits, 5) + np.testing.assert_equal(id_5.coeffs, np.array([1], dtype=complex)) + np.testing.assert_equal(id_5.bit_terms, np.array([], dtype=np.uint8)) + np.testing.assert_equal(id_5.indices, np.array([], dtype=np.uint32)) + np.testing.assert_equal(id_5.boundaries, np.array([0, 0], dtype=np.uintp)) + + id_0 = SparseObservable.identity(0) + self.assertEqual(id_0.num_qubits, 0) + np.testing.assert_equal(id_0.coeffs, np.array([1], dtype=complex)) + np.testing.assert_equal(id_0.bit_terms, np.array([], dtype=np.uint8)) + np.testing.assert_equal(id_0.indices, np.array([], dtype=np.uint32)) + np.testing.assert_equal(id_0.boundaries, np.array([0, 0], dtype=np.uintp)) + + @ddt.data( + SparseObservable.zero(0), + SparseObservable.zero(5), + SparseObservable.identity(0), + SparseObservable.identity(1_000_000), + SparseObservable.from_list([("+-rl01", -0.5), ("IXZYZI", 1.0j)]), + ) + def test_copy(self, obs): + self.assertEqual(obs, obs.copy()) + self.assertIsNot(obs, obs.copy()) + + def test_equality(self): + sparse_data = [("XZ", (1, 0), 0.5j), ("+lr", (3, 1, 0), -0.25)] + op = SparseObservable.from_sparse_list(sparse_data, num_qubits=5) + self.assertEqual(op, op.copy()) + # Take care that Rust space allows multiple views onto the same object. + self.assertEqual(op, op) + + # No costly automatic simplification (mathematically, these operators _are_ the same). + self.assertNotEqual( + SparseObservable.from_list([("+", 1.0), ("-", 1.0)]), SparseObservable.from_label("X") + ) + + # Difference in qubit count. + self.assertNotEqual( + op, SparseObservable.from_sparse_list(sparse_data, num_qubits=op.num_qubits + 1) + ) + self.assertNotEqual(SparseObservable.zero(2), SparseObservable.zero(3)) + self.assertNotEqual(SparseObservable.identity(2), SparseObservable.identity(3)) + + # Difference in coeffs. + self.assertNotEqual( + SparseObservable.from_list([("IIXZI", 1.0), ("+-rl0", -0.5j)]), + SparseObservable.from_list([("IIXZI", 1.0), ("+-rl0", 0.5j)]), + ) + self.assertNotEqual( + SparseObservable.from_list([("IIXZI", 1.0), ("+-rl0", -0.5j)]), + SparseObservable.from_list([("IIXZI", 1.0j), ("+-rl0", -0.5j)]), + ) + + # Difference in bit terms. + self.assertNotEqual( + SparseObservable.from_list([("IIXZI", 1.0), ("+-rl0", -0.5j)]), + SparseObservable.from_list([("IIYZI", 1.0), ("+-rl0", -0.5j)]), + ) + self.assertNotEqual( + SparseObservable.from_list([("IIXZI", 1.0), ("+-rl0", -0.5j)]), + SparseObservable.from_list([("IIXZI", 1.0), ("+-rl1", -0.5j)]), + ) + + # Difference in indices. + self.assertNotEqual( + SparseObservable.from_list([("IIXZI", 1.0), ("+Irl0", -0.5j)]), + SparseObservable.from_list([("IXIZI", 1.0), ("+Irl0", -0.5j)]), + ) + self.assertNotEqual( + SparseObservable.from_list([("IIXZI", 1.0), ("+Irl0", -0.5j)]), + SparseObservable.from_list([("IIXZI", 1.0), ("I+rl0", -0.5j)]), + ) + + # Difference in boundaries. + self.assertNotEqual( + SparseObservable.from_sparse_list( + [("XZ", (0, 1), 1.5), ("+-", (2, 3), -0.5j)], num_qubits=5 + ), + SparseObservable.from_sparse_list( + [("XZ+", (0, 1, 2), 1.5), ("-", (3,), -0.5j)], num_qubits=5 + ), + ) + + def test_write_into_attributes_scalar(self): + coeffs = SparseObservable.from_sparse_list( + [("XZ", (1, 0), 1.5j), ("+-", (3, 2), -1.5j)], num_qubits=8 + ) + coeffs.coeffs[0] = -2.0 + self.assertEqual( + coeffs, + SparseObservable.from_sparse_list( + [("XZ", (1, 0), -2.0), ("+-", (3, 2), -1.5j)], num_qubits=8 + ), + ) + coeffs.coeffs[1] = 1.5 + 0.25j + self.assertEqual( + coeffs, + SparseObservable.from_sparse_list( + [("XZ", (1, 0), -2.0), ("+-", (3, 2), 1.5 + 0.25j)], num_qubits=8 + ), + ) + + bit_terms = SparseObservable.from_sparse_list( + [("XZ", (0, 1), 1.5j), ("+-", (2, 3), -1.5j)], num_qubits=8 + ) + bit_terms.bit_terms[0] = SparseObservable.BitTerm.Y + bit_terms.bit_terms[3] = SparseObservable.BitTerm.LEFT + self.assertEqual( + bit_terms, + SparseObservable.from_sparse_list( + [("YZ", (0, 1), 1.5j), ("+l", (2, 3), -1.5j)], num_qubits=8 + ), + ) + + indices = SparseObservable.from_sparse_list( + [("XZ", (0, 1), 1.5j), ("+-", (2, 3), -1.5j)], num_qubits=8 + ) + # These two sets keep the observable in term-wise increasing order. We don't test what + # happens if somebody violates the Rust-space requirement to be term-wise increasing. + indices.indices[1] = 4 + indices.indices[3] = 7 + self.assertEqual( + indices, + SparseObservable.from_sparse_list( + [("XZ", (0, 4), 1.5j), ("+-", (2, 7), -1.5j)], num_qubits=8 + ), + ) + + boundaries = SparseObservable.from_sparse_list( + [("XZ", (0, 1), 1.5j), ("+-", (2, 3), -1.5j)], num_qubits=8 + ) + # Move a single-qubit term from the second summand into the first (the particular indices + # ensure we remain term-wise sorted). + boundaries.boundaries[1] += 1 + self.assertEqual( + boundaries, + SparseObservable.from_sparse_list( + [("XZ+", (0, 1, 2), 1.5j), ("-", (3,), -1.5j)], num_qubits=8 + ), + ) + + def test_write_into_attributes_broadcast(self): + coeffs = SparseObservable.from_list([("XIIZI", 1.5j), ("IIIl0", -0.25), ("1IIIY", 0.5)]) + coeffs.coeffs[:] = 1.5j + np.testing.assert_array_equal(coeffs.coeffs, [1.5j, 1.5j, 1.5j]) + coeffs.coeffs[1:] = 1.0j + np.testing.assert_array_equal(coeffs.coeffs, [1.5j, 1.0j, 1.0j]) + coeffs.coeffs[:2] = -0.5 + np.testing.assert_array_equal(coeffs.coeffs, [-0.5, -0.5, 1.0j]) + coeffs.coeffs[::2] = 1.5j + np.testing.assert_array_equal(coeffs.coeffs, [1.5j, -0.5, 1.5j]) + coeffs.coeffs[::-1] = -0.5j + np.testing.assert_array_equal(coeffs.coeffs, [-0.5j, -0.5j, -0.5j]) + + # It's hard to broadcast into `indices` without breaking data coherence; the broadcasting is + # more meant for fast modifications to `coeffs` and `bit_terms`. + indices = SparseObservable.from_list([("XIIZI", 1.5j), ("IIlI0", -0.25), ("1IIIY", 0.5)]) + indices.indices[::2] = 1 + self.assertEqual( + indices, SparseObservable.from_list([("XIIZI", 1.5j), ("IIl0I", -0.25), ("1IIYI", 0.5)]) + ) + + bit_terms = SparseObservable.from_list([("XIIZI", 1.5j), ("IIlI0", -0.25), ("1IIIY", 0.5)]) + bit_terms.bit_terms[::2] = SparseObservable.BitTerm.Z + self.assertEqual( + bit_terms, + SparseObservable.from_list([("XIIZI", 1.5j), ("IIlIZ", -0.25), ("1IIIZ", 0.5)]), + ) + bit_terms.bit_terms[3:1:-1] = SparseObservable.BitTerm.PLUS + self.assertEqual( + bit_terms, + SparseObservable.from_list([("XIIZI", 1.5j), ("II+I+", -0.25), ("1IIIZ", 0.5)]), + ) + bit_terms.bit_terms[bit_terms.boundaries[2] : bit_terms.boundaries[3]] = ( + SparseObservable.BitTerm.MINUS + ) + self.assertEqual( + bit_terms, + SparseObservable.from_list([("XIIZI", 1.5j), ("II+I+", -0.25), ("-III-", 0.5)]), + ) + + boundaries = SparseObservable.from_list([("IIIIZX", 1j), ("II+-II", -0.5), ("rlIIII", 0.5)]) + boundaries.boundaries[1:3] = 1 + self.assertEqual( + boundaries, + SparseObservable.from_list([("IIIIIX", 1j), ("IIIIII", -0.5), ("rl+-ZI", 0.5)]), + ) + + def test_write_into_attributes_slice(self): + coeffs = SparseObservable.from_list([("XIIZI", 1.5j), ("IIIl0", -0.25), ("1IIIY", 0.5)]) + coeffs.coeffs[:] = [2.0, 0.5, -0.25] + self.assertEqual( + coeffs, SparseObservable.from_list([("XIIZI", 2.0), ("IIIl0", 0.5), ("1IIIY", -0.25)]) + ) + # This should assign the coefficients in reverse order - we more usually spell it + # `coeffs[:] = coeffs{::-1]`, but the idea is to check the set-item slicing order. + coeffs.coeffs[::-1] = coeffs.coeffs[:] + self.assertEqual( + coeffs, SparseObservable.from_list([("XIIZI", -0.25), ("IIIl0", 0.5), ("1IIIY", 2.0)]) + ) + + indices = SparseObservable.from_list([("IIIIZX", 0.25), ("II+-II", 1j), ("rlIIII", 0.5)]) + indices.indices[:4] = [4, 5, 1, 2] + self.assertEqual( + indices, SparseObservable.from_list([("ZXIIII", 0.25), ("III+-I", 1j), ("rlIIII", 0.5)]) + ) + + bit_terms = SparseObservable.from_list([("IIIIZX", 0.25), ("II+-II", 1j), ("rlIIII", 0.5)]) + bit_terms.bit_terms[::2] = [ + SparseObservable.BitTerm.Y, + SparseObservable.BitTerm.RIGHT, + SparseObservable.BitTerm.ZERO, + ] + self.assertEqual( + bit_terms, + SparseObservable.from_list([("IIIIZY", 0.25), ("II+rII", 1j), ("r0IIII", 0.5)]), + ) + + operators = SparseObservable.from_list([("XZY", 1.5j), ("+1r", -0.5)]) + # Reduce all single-qubit terms to the relevant Pauli operator, if they are a projector. + operators.bit_terms[:] = operators.bit_terms[:] & 0b00_11 + self.assertEqual(operators, SparseObservable.from_list([("XZY", 1.5j), ("XZY", -0.5)])) + + boundaries = SparseObservable.from_list([("IIIIZX", 0.25), ("II+-II", 1j), ("rlIIII", 0.5)]) + boundaries.boundaries[1:-1] = [1, 5] + self.assertEqual( + boundaries, + SparseObservable.from_list([("IIIIIX", 0.25), ("Il+-ZI", 1j), ("rIIIII", 0.5)]), + ) + + def test_attributes_reject_bad_writes(self): + obs = SparseObservable.from_list([("XZY", 1.5j), ("+-r", -0.5)]) + with self.assertRaises(TypeError): + obs.coeffs[0] = [0.25j, 0.5j] + with self.assertRaises(TypeError): + obs.bit_terms[0] = [SparseObservable.BitTerm.PLUS] * 4 + with self.assertRaises(TypeError): + obs.indices[0] = [0, 1] + with self.assertRaises(TypeError): + obs.boundaries[0] = (0, 1) + with self.assertRaisesRegex(ValueError, "not a valid letter"): + obs.bit_terms[0] = 0 + with self.assertRaisesRegex(ValueError, "not a valid letter"): + obs.bit_terms[:] = 0 + with self.assertRaisesRegex( + ValueError, "tried to set a slice of length 2 with a sequence of length 1" + ): + obs.coeffs[:] = [1.0j] + with self.assertRaisesRegex( + ValueError, "tried to set a slice of length 6 with a sequence of length 8" + ): + obs.bit_terms[:] = [SparseObservable.BitTerm.Z] * 8 + + def test_attributes_sequence(self): + """Test attributes of the `Sequence` protocol.""" + # Length + obs = SparseObservable.from_list([("XZY", 1.5j), ("+-r", -0.5)]) + self.assertEqual(len(obs.coeffs), 2) + self.assertEqual(len(obs.indices), 6) + self.assertEqual(len(obs.bit_terms), 6) + self.assertEqual(len(obs.boundaries), 3) + + # Iteration + self.assertEqual(list(obs.coeffs), [1.5j, -0.5]) + self.assertEqual(tuple(obs.indices), (0, 1, 2, 0, 1, 2)) + self.assertEqual(next(iter(obs.boundaries)), 0) + # multiple iteration through same object + bit_terms = obs.bit_terms + self.assertEqual(set(bit_terms), {SparseObservable.BitTerm[x] for x in "XYZ+-r"}) + self.assertEqual(set(bit_terms), {SparseObservable.BitTerm[x] for x in "XYZ+-r"}) + + # Implicit iteration methods. + self.assertIn(SparseObservable.BitTerm.PLUS, obs.bit_terms) + self.assertNotIn(4, obs.indices) + self.assertEqual(list(reversed(obs.coeffs)), [-0.5, 1.5j]) + + # Index by scalar + self.assertEqual(obs.coeffs[1], -0.5) + self.assertEqual(obs.indices[-1], 2) + self.assertEqual(obs.bit_terms[0], SparseObservable.BitTerm.Y) + self.assertEqual(obs.boundaries[-2], 3) + with self.assertRaises(IndexError): + _ = obs.coeffs[10] + with self.assertRaises(IndexError): + _ = obs.boundaries[-4] + + # Index by slice. This is API guaranteed to be a Numpy array to make it easier to + # manipulate subslices with mathematic operations. + self.assertIsInstance(obs.coeffs[:], np.ndarray) + np.testing.assert_array_equal( + obs.coeffs[:], np.array([1.5j, -0.5], dtype=np.complex128), strict=True + ) + self.assertIsInstance(obs.indices[::-1], np.ndarray) + np.testing.assert_array_equal( + obs.indices[::-1], np.array([2, 1, 0, 2, 1, 0], dtype=np.uint32), strict=True + ) + self.assertIsInstance(obs.bit_terms[2:4], np.ndarray) + np.testing.assert_array_equal( + obs.bit_terms[2:4], + np.array([SparseObservable.BitTerm.X, SparseObservable.BitTerm.RIGHT], dtype=np.uint8), + strict=True, + ) + self.assertIsInstance(obs.boundaries[-2:-3:-1], np.ndarray) + np.testing.assert_array_equal( + obs.boundaries[-2:-3:-1], np.array([3], dtype=np.uintp), strict=True + ) + + def test_attributes_to_array(self): + obs = SparseObservable.from_list([("XZY", 1.5j), ("+-r", -0.5)]) + + # Natural dtypes. + np.testing.assert_array_equal( + obs.coeffs, np.array([1.5j, -0.5], dtype=np.complex128), strict=True + ) + np.testing.assert_array_equal( + obs.indices, np.array([0, 1, 2, 0, 1, 2], dtype=np.uint32), strict=True + ) + np.testing.assert_array_equal( + obs.bit_terms, + np.array([SparseObservable.BitTerm[x] for x in "YZXr-+"], dtype=np.uint8), + strict=True, + ) + np.testing.assert_array_equal( + obs.boundaries, np.array([0, 3, 6], dtype=np.uintp), strict=True + ) + + # Cast dtypes. + np.testing.assert_array_equal( + np.array(obs.indices, dtype=np.uint8), + np.array([0, 1, 2, 0, 1, 2], dtype=np.uint8), + strict=True, + ) + np.testing.assert_array_equal( + np.array(obs.boundaries, dtype=np.int64), + np.array([0, 3, 6], dtype=np.int64), + strict=True, + ) + + @unittest.skipIf( + int(np.__version__.split(".", maxsplit=1)[0]) < 2, + "Numpy 1.x did not have a 'copy' keyword parameter to 'numpy.asarray'", + ) + def test_attributes_reject_no_copy_array(self): + obs = SparseObservable.from_list([("XZY", 1.5j), ("+-r", -0.5)]) + with self.assertRaisesRegex(ValueError, "cannot produce a safe view"): + np.asarray(obs.coeffs, copy=False) + with self.assertRaisesRegex(ValueError, "cannot produce a safe view"): + np.asarray(obs.indices, copy=False) + with self.assertRaisesRegex(ValueError, "cannot produce a safe view"): + np.asarray(obs.bit_terms, copy=False) + with self.assertRaisesRegex(ValueError, "cannot produce a safe view"): + np.asarray(obs.boundaries, copy=False) + + def test_attributes_repr(self): + # We're not testing much about the outputs here, just that they don't crash. + obs = SparseObservable.from_list([("XZY", 1.5j), ("+-r", -0.5)]) + self.assertIn("coeffs", repr(obs.coeffs)) + self.assertIn("bit_terms", repr(obs.bit_terms)) + self.assertIn("indices", repr(obs.indices)) + self.assertIn("boundaries", repr(obs.boundaries))