Add infrastructure for gates, instruction, and operations in Rust

.github/workflows/tests.yml

@@ -36,6 +36,15 @@ jobs:
           python -m pip install -U -r requirements.txt -c constraints.txt
           python -m pip install -U -r requirements-dev.txt -c constraints.txt
           python -m pip install -c constraints.txt -e .
+        if: matrix.python-version == '3.10'
+        env:
+          QISKIT_NO_CACHE_GATES: 1
+      - name: 'Install dependencies'
+        run: |
+          python -m pip install -U -r requirements.txt -c constraints.txt
+          python -m pip install -U -r requirements-dev.txt -c constraints.txt
+          python -m pip install -c constraints.txt -e .
+        if: matrix.python-version == '3.12'
       - name: 'Install optionals'
         run: |
           python -m pip install -r requirements-optional.txt -c constraints.txt

CONTRIBUTING.md

@@ -135,6 +135,18 @@ Note that in order to run `python setup.py ...` commands you need have build
 dependency packages installed in your environment, which are listed in the
 `pyproject.toml` file under the `[build-system]` section.
 
+### Compile time options
+
+When building qiskit from source there are options available to control how
+Qiskit is built. Right now the only option is if you set the environment
+variable `QISKIT_NO_CACHE_GATES=1` this will disable runtime caching of
+Python gate objects when accessing them from a `QuantumCircuit` or `DAGCircuit`.
+This makes a tradeoff between runtime performance for Python access and memory
+overhead. Caching gates will result in better runtime for users of Python at
+the cost of increased memory consumption. If you're working with any custom
+transpiler passes written in python or are otherwise using a workflow that
+repeatedly accesses the `operation` attribute of a `CircuitInstruction` or `op`
+attribute of `DAGOpNode` enabling caching is recommended.
 
 ## Issues and pull requests
 

Cargo.toml

@@ -16,6 +16,11 @@ license = "Apache-2.0"
 [workspace.dependencies]
 indexmap.version = "2.2.6"
 hashbrown.version = "0.14.0"
+num-complex = "0.4"
+ndarray = "^0.15.6"
+numpy = "0.21.0"
+smallvec = "1.13"
+
 # Most of the crates don't need the feature `extension-module`, since only `qiskit-pyext` builds an
 # actual C extension (the feature disables linking in `libpython`, which is forbidden in Python
 # distributions).  We only activate that feature when building the C extension module; we still need

crates/accelerate/Cargo.toml

@@ -11,29 +11,29 @@ doctest = false
 
 [dependencies]
 rayon = "1.10"
-numpy = "0.21.0"
+numpy.workspace = true
 rand = "0.8"
 rand_pcg = "0.3"
 rand_distr = "0.4.3"
 ahash = "0.8.11"
 num-traits = "0.2"
-num-complex = "0.4"
+num-complex.workspace = true
 num-bigint = "0.4"
 rustworkx-core = "0.14"
 faer = "0.19.0"
 itertools = "0.13.0"
 qiskit-circuit.workspace = true
 
 [dependencies.smallvec]
-version = "1.13"
+workspace = true
 features = ["union"]
 
 [dependencies.pyo3]
 workspace = true
 features = ["hashbrown", "indexmap", "num-complex", "num-bigint", "smallvec"]
 
 [dependencies.ndarray]
-version = "^0.15.6"
+workspace = true
 features = ["rayon", "approx-0_5"]
 
 [dependencies.approx]

crates/accelerate/src/isometry.rs

@@ -23,7 +23,7 @@ use itertools::Itertools;
 use ndarray::prelude::*;
 use numpy::{IntoPyArray, PyReadonlyArray1, PyReadonlyArray2};
 
-use crate::two_qubit_decompose::ONE_QUBIT_IDENTITY;
+use qiskit_circuit::gate_matrix::ONE_QUBIT_IDENTITY;
 
 /// Find special unitary matrix that maps [c0,c1] to [r,0] or [0,r] if basis_state=0 or
 /// basis_state=1 respectively

crates/accelerate/src/two_qubit_decompose.rs

@@ -51,6 +51,7 @@ use rand::prelude::*;
 use rand_distr::StandardNormal;
 use rand_pcg::Pcg64Mcg;
 
+use qiskit_circuit::gate_matrix::{CX_GATE, H_GATE, ONE_QUBIT_IDENTITY, SX_GATE, X_GATE};
 use qiskit_circuit::SliceOrInt;
 
 const PI2: f64 = PI / 2.0;
@@ -60,11 +61,6 @@ const TWO_PI: f64 = 2.0 * PI;
 
 const C1: c64 = c64 { re: 1.0, im: 0.0 };
 
-pub static ONE_QUBIT_IDENTITY: [[Complex64; 2]; 2] = [
-    [Complex64::new(1., 0.), Complex64::new(0., 0.)],
-    [Complex64::new(0., 0.), Complex64::new(1., 0.)],
-];
-
 static B_NON_NORMALIZED: [[Complex64; 4]; 4] = [
     [
         Complex64::new(1.0, 0.),
@@ -342,54 +338,6 @@ fn rz_matrix(theta: f64) -> Array2<Complex64> {
     ]
 }
 
-static HGATE: [[Complex64; 2]; 2] = [
-    [
-        Complex64::new(FRAC_1_SQRT_2, 0.),
-        Complex64::new(FRAC_1_SQRT_2, 0.),
-    ],
-    [
-        Complex64::new(FRAC_1_SQRT_2, 0.),
-        Complex64::new(-FRAC_1_SQRT_2, 0.),
-    ],
-];
-
-static CXGATE: [[Complex64; 4]; 4] = [
-    [
-        Complex64::new(1., 0.),
-        Complex64::new(0., 0.),
-        Complex64::new(0., 0.),
-        Complex64::new(0., 0.),
-    ],
-    [
-        Complex64::new(0., 0.),
-        Complex64::new(0., 0.),
-        Complex64::new(0., 0.),
-        Complex64::new(1., 0.),
-    ],
-    [
-        Complex64::new(0., 0.),
-        Complex64::new(0., 0.),
-        Complex64::new(1., 0.),
-        Complex64::new(0., 0.),
-    ],
-    [
-        Complex64::new(0., 0.),
-        Complex64::new(1., 0.),
-        Complex64::new(0., 0.),
-        Complex64::new(0., 0.),
-    ],
-];
-
-static SXGATE: [[Complex64; 2]; 2] = [
-    [Complex64::new(0.5, 0.5), Complex64::new(0.5, -0.5)],
-    [Complex64::new(0.5, -0.5), Complex64::new(0.5, 0.5)],
-];
-
-static XGATE: [[Complex64; 2]; 2] = [
-    [Complex64::new(0., 0.), Complex64::new(1., 0.)],
-    [Complex64::new(1., 0.), Complex64::new(0., 0.)],
-];
-
 fn compute_unitary(sequence: &TwoQubitSequenceVec, global_phase: f64) -> Array2<Complex64> {
     let identity = aview2(&ONE_QUBIT_IDENTITY);
     let phase = Complex64::new(0., global_phase).exp();
@@ -402,10 +350,10 @@ fn compute_unitary(sequence: &TwoQubitSequenceVec, global_phase: f64) -> Array2<
             // sequence. If we get a different gate this is getting called
             // by something else and is invalid.
             let gate_matrix = match inst.0.as_ref() {
-                "sx" => aview2(&SXGATE).to_owned(),
+                "sx" => aview2(&SX_GATE).to_owned(),
                 "rz" => rz_matrix(inst.1[0]),
-                "cx" => aview2(&CXGATE).to_owned(),
-                "x" => aview2(&XGATE).to_owned(),
+                "cx" => aview2(&CX_GATE).to_owned(),
+                "x" => aview2(&X_GATE).to_owned(),
                 _ => unreachable!("Undefined gate"),
             };
             (gate_matrix, &inst.2)
@@ -1481,15 +1429,15 @@ impl TwoQubitBasisDecomposer {
         } else {
             euler_matrix_q0 = rz_matrix(euler_q0[0][2] + euler_q0[1][0]).dot(&euler_matrix_q0);
         }
-        euler_matrix_q0 = aview2(&HGATE).dot(&euler_matrix_q0);
+        euler_matrix_q0 = aview2(&H_GATE).dot(&euler_matrix_q0);
         self.append_1q_sequence(&mut gates, &mut global_phase, euler_matrix_q0.view(), 0);
 
         let rx_0 = rx_matrix(euler_q1[0][0]);
         let rz = rz_matrix(euler_q1[0][1]);
         let rx_1 = rx_matrix(euler_q1[0][2] + euler_q1[1][0]);
         let mut euler_matrix_q1 = rz.dot(&rx_0);
         euler_matrix_q1 = rx_1.dot(&euler_matrix_q1);
-        euler_matrix_q1 = aview2(&HGATE).dot(&euler_matrix_q1);
+        euler_matrix_q1 = aview2(&H_GATE).dot(&euler_matrix_q1);
         self.append_1q_sequence(&mut gates, &mut global_phase, euler_matrix_q1.view(), 1);
 
         gates.push(("cx".to_string(), smallvec![], smallvec![1, 0]));
@@ -1550,12 +1498,12 @@ impl TwoQubitBasisDecomposer {
             return None;
         }
         gates.push(("cx".to_string(), smallvec![], smallvec![1, 0]));
-        let mut euler_matrix = rz_matrix(euler_q0[2][2] + euler_q0[3][0]).dot(&aview2(&HGATE));
+        let mut euler_matrix = rz_matrix(euler_q0[2][2] + euler_q0[3][0]).dot(&aview2(&H_GATE));
         euler_matrix = rx_matrix(euler_q0[3][1]).dot(&euler_matrix);
         euler_matrix = rz_matrix(euler_q0[3][2]).dot(&euler_matrix);
         self.append_1q_sequence(&mut gates, &mut global_phase, euler_matrix.view(), 0);
 
-        let mut euler_matrix = rx_matrix(euler_q1[2][2] + euler_q1[3][0]).dot(&aview2(&HGATE));
+        let mut euler_matrix = rx_matrix(euler_q1[2][2] + euler_q1[3][0]).dot(&aview2(&H_GATE));
         euler_matrix = rz_matrix(euler_q1[3][1]).dot(&euler_matrix);
         euler_matrix = rx_matrix(euler_q1[3][2]).dot(&euler_matrix);
         self.append_1q_sequence(&mut gates, &mut global_phase, euler_matrix.view(), 1);

crates/circuit/Cargo.toml

@@ -11,4 +11,17 @@ doctest = false
 
 [dependencies]
 hashbrown.workspace = true
-pyo3.workspace = true
+num-complex.workspace = true
+ndarray.workspace = true
+numpy.workspace = true
+
+[dependencies.pyo3]
+workspace = true
+features = ["hashbrown", "indexmap", "num-complex", "num-bigint", "smallvec"]
+
+[dependencies.smallvec]
+workspace = true
+features = ["union"]
+
+[features]
+cache_pygates = []

crates/circuit/README.md

@@ -4,3 +4,66 @@ The Rust-based data structures for circuits.
 This currently defines the core data collections for `QuantumCircuit`, but may expand in the future to back `DAGCircuit` as well.
 
 This crate is a very low part of the Rust stack, if not the very lowest.
+
+The data model exposed by this crate is as follows.
+
+## CircuitData
+
+The core representation of a quantum circuit in Rust is the `CircuitData` struct. This containts the list
+of instructions that are comprising the circuit. Each element in this list is modeled by a
+`CircuitInstruction` struct. The `CircuitInstruction` contains the operation object and it's operands.
+This includes the parameters and bits. It also contains the potential mutable state of the Operation representation from the legacy Python data model; namely `duration`, `unit`, `condition`, and `label`.
+In the future we'll be able to remove all of that except for label.
+
+At rest a `CircuitInstruction` is compacted into a `PackedInstruction` which caches reused qargs
+in the instructions to reduce the memory overhead of `CircuitData`. The `PackedInstruction` objects
+get unpacked back to `CircuitInstruction` when accessed for a more convienent working form.
+
+Additionally the `CircuitData` contains a `param_table` field which is used to track parameterized
+instructions that are using python defined `ParameterExpression` objects for any parameters and also
+a global phase field which is used to track the global phase of the circuit.
+
+## Operation Model
+
+In the circuit crate all the operations used in a `CircuitInstruction` are part of the `OperationType`
+enum. The `OperationType` enum has four variants which are used to define the different types of
+operation objects that can be on a circuit:
+
+ - `StandardGate`: a rust native representation of a member of the Qiskit standard gate library. This is
+    an `enum` that enumerates all the gates in the library and statically defines all the gate properties
+    except for gates that take parameters,
+ - `PyGate`: A struct that wraps a gate outside the standard library defined in Python. This struct wraps
+    a `Gate` instance (or subclass) as a `PyObject`. The static properties of this object (such as name,
+    number of qubits, etc) are stored in Rust for performance but the dynamic properties such as
+    the matrix or definition are accessed by calling back into Python to get them from the stored
+    `PyObject`
+ - `PyInstruction`: A struct that wraps an instruction defined in Python. This struct wraps an
+    `Instruction` instance (or subclass) as a `PyObject`. The static properties of this object (such as
+    name, number of qubits, etc) are stored in Rust for performance but the dynamic properties such as
+    the definition are accessed by calling back into Python to get them from the stored `PyObject`. As
+    the primary difference between `Gate` and `Instruction` in the python data model are that `Gate` is a
+    specialized `Instruction` subclass that represents unitary operations the primary difference between
+    this and `PyGate` are that `PyInstruction` will always return `None` when it's matrix is accessed.
+ - `PyOperation`: A struct that wraps an operation defined in Python. This struct wraps an `Operation`
+    instance (or subclass) as a `PyObject`. The static properties of this object (such as name, number
+    of qubits, etc) are stored in Rust for performance. As `Operation` is the base abstract interface
+    definition of what can be put on a circuit this is mostly just a container for custom Python objects.
+    Anything that's operating on a bare operation will likely need to access it via the `PyObject`
+    manually because the interface doesn't define many standard properties outside of what's cached in
+    the struct.
+
+There is also an `Operation` trait defined which defines the common access pattern interface to these
+4 types along with the `OperationType` parent. This trait defines methods to access the standard data
+model attributes of operations in Qiskit. This includes things like the name, number of qubits, the matrix, the definition, etc.
+
+## ParamTable
+
+The `ParamTable` struct is used to track which circuit instructions are using `ParameterExpression`
+objects for any of their parameters. The Python space `ParameterExpression` is comprised of a symengine
+symbolic expression that defines operations using `Parameter` objects. Each `Parameter` is modeled by
+a uuid and a name to uniquely identify it. The parameter table maps the `Parameter` objects to the
+`CircuitInstruction` in the `CircuitData` that are using them. The `Parameter` comprised of 3 `HashMaps` internally that map the uuid (as `u128`, which is accesible in Python by using `uuid.int`) to the `ParamEntry`, the `name` to the uuid, and the uuid to the PyObject for the actual `Parameter`.
+
+The `ParamEntry` is just a `HashSet` of 2-tuples with usize elements. The two usizes represent the instruction index in the `CircuitData` and the index of the `CircuitInstruction.params` field of
+a give instruction where the given `Parameter` is used in the circuit. If the instruction index is
+`GLOBAL_PHASE_MAX`, that points to the global phase property of the circuit instead of a `CircuitInstruction`.

crates/circuit/src/bit_data.rs

@@ -12,7 +12,7 @@
 
 use crate::BitType;
 use hashbrown::HashMap;
-use pyo3::exceptions::{PyRuntimeError, PyValueError};
+use pyo3::exceptions::{PyKeyError, PyRuntimeError, PyValueError};
 use pyo3::prelude::*;
 use pyo3::types::PyList;
 use std::fmt::Debug;
@@ -83,6 +83,15 @@ pub(crate) struct BitData<T> {
 
 pub(crate) struct BitNotFoundError<'py>(pub(crate) Bound<'py, PyAny>);
 
+impl<'py> From<BitNotFoundError<'py>> for PyErr {
+    fn from(error: BitNotFoundError) -> Self {
+        PyKeyError::new_err(format!(
+            "Bit {:?} has not been added to this circuit.",
+            error.0
+        ))
+    }
+}
+
 impl<T> BitData<T>
 where
     T: From<BitType> + Copy,
@@ -142,7 +151,7 @@ where
     /// Map the provided native indices to the corresponding Python
     /// bit instances.
     /// Panics if any of the indices are out of range.
-    pub fn map_indices(&self, bits: &[T]) -> impl Iterator<Item = &Py<PyAny>> + ExactSizeIterator {
+    pub fn map_indices(&self, bits: &[T]) -> impl ExactSizeIterator<Item = &Py<PyAny>> {
         let v: Vec<_> = bits.iter().map(|i| self.get(*i).unwrap()).collect();
         v.into_iter()
     }