Fully port Optimize1qGatesDecomposition to Rust

This commit builds off of Qiskit#12550 and the other data model in Rust
infrastructure and migrates the Optimize1qGatesDecomposition pass to
operate fully in Rust. The full path of the transpiler pass now never
leaves rust until it has finished modifying the DAGCircuit. There is
still some python interaction necessary to handle parts of the data
model that are still in Python, mainly calibrations and parameter
expressions (for global phase). But otherwise the entirety of the pass
operates in rust now.

This is just a first pass at the migration here, it moves the pass to be
a single for loop in rust. The next steps here are to look at operating
the pass in parallel. There is no data dependency between the
optimizations being done by the pass so we should be able to the
throughput of the pass by leveraging multithreading to handle each run
in parallel. This commit does not attempt this though, because of the
Python dependency and also the data structures around gates and the
dag aren't really setup for multithreading yet and there likely will
need to be some work to support that (this pass is a good candidate to
work through the bugs on that).

Part of Qiskit#12208

crates/accelerate/src/euler_one_qubit_decomposer.rs

@@ -13,7 +13,8 @@
 #![allow(clippy::too_many_arguments)]
 #![allow(clippy::upper_case_acronyms)]
 
-use hashbrown::HashMap;
+use hashbrown::{HashMap, HashSet};
+use indexmap::IndexSet;
 use num_complex::{Complex64, ComplexFloat};
 use smallvec::{smallvec, SmallVec};
 use std::cmp::Ordering;
@@ -29,14 +30,19 @@ use pyo3::Python;
 use ndarray::prelude::*;
 use numpy::PyReadonlyArray2;
 use pyo3::pybacked::PyBackedStr;
+use rustworkx_core::petgraph::stable_graph::NodeIndex;
 
 use qiskit_circuit::circuit_data::CircuitData;
+use qiskit_circuit::dag_circuit::{DAGCircuit, NodeType};
 use qiskit_circuit::dag_node::DAGOpNode;
 use qiskit_circuit::operations::{Operation, Param, StandardGate};
 use qiskit_circuit::slice::{PySequenceIndex, SequenceIndex};
 use qiskit_circuit::util::c64;
 use qiskit_circuit::Qubit;
 
+use crate::nlayout::PhysicalQubit;
+use crate::target_transpiler::Target;
+
 pub const ANGLE_ZERO_EPSILON: f64 = 1e-12;
 
 #[pyclass(module = "qiskit._accelerate.euler_one_qubit_decomposer")]
@@ -69,6 +75,7 @@ impl OneQubitGateErrorMap {
     }
 }
 
+#[derive(Debug)]
 #[pyclass(sequence)]
 pub struct OneQubitGateSequence {
     pub gates: Vec<(StandardGate, SmallVec<[f64; 3]>)>,
@@ -571,7 +578,7 @@ pub fn generate_circuit(
     Ok(res)
 }
 
-#[derive(Clone, Debug, Copy)]
+#[derive(Clone, Debug, Copy, Eq, Hash, PartialEq)]
 #[pyclass(module = "qiskit._accelerate.euler_one_qubit_decomposer")]
 pub enum EulerBasis {
     U321,
@@ -965,72 +972,221 @@ pub fn params_zxz(unitary: PyReadonlyArray2<Complex64>) -> [f64; 4] {
     params_zxz_inner(mat)
 }
 
-type OptimizeDecompositionReturn = Option<((f64, usize), (f64, usize), OneQubitGateSequence)>;
+fn compute_error_term_from_target(gate: &str, target: &Target, qubit: PhysicalQubit) -> f64 {
+    1. - target.get_error(gate, &[qubit]).unwrap_or(0.)
+}
+
+fn compute_error_from_target_one_qubit_sequence(
+    circuit: &OneQubitGateSequence,
+    qubit: PhysicalQubit,
+    target: Option<&Target>,
+) -> (f64, usize) {
+    match target {
+        Some(target) => {
+            let num_gates = circuit.gates.len();
+            let gate_fidelities: f64 = circuit
+                .gates
+                .iter()
+                .map(|gate| compute_error_term_from_target(gate.0.name(), target, qubit))
+                .product();
+            (1. - gate_fidelities, num_gates)
+        }
+        None => (circuit.gates.len() as f64, circuit.gates.len()),
+    }
+}
 
 #[pyfunction]
-pub fn optimize_1q_gates_decomposition(
-    runs: Vec<Vec<PyRef<DAGOpNode>>>,
-    qubits: Vec<usize>,
-    bases: Vec<Vec<PyBackedStr>>,
-    simplify: bool,
-    error_map: Option<&OneQubitGateErrorMap>,
-    atol: Option<f64>,
-) -> Vec<OptimizeDecompositionReturn> {
-    runs.iter()
-        .enumerate()
-        .map(|(index, raw_run)| -> OptimizeDecompositionReturn {
-            let mut error = match error_map {
-                Some(_) => 1.,
-                None => raw_run.len() as f64,
-            };
-            let qubit = qubits[index];
-            let operator = &raw_run
+#[pyo3(signature = (dag, *, target=None, basis_gates=None, global_decomposers=None))]
+pub(crate) fn optimize_1q_gates_decomposition(
+    py: Python,
+    dag: &mut DAGCircuit,
+    target: Option<&Target>,
+    basis_gates: Option<HashSet<String>>,
+    global_decomposers: Option<Vec<String>>,
+) -> PyResult<()> {
+    let runs: Vec<Vec<NodeIndex>> = dag.collect_1q_runs().unwrap().collect();
+    for raw_run in runs {
+        let mut error = match target {
+            Some(_) => 1.,
+            None => raw_run.len() as f64,
+        };
+        let qubit = if let NodeType::Operation(inst) = &dag.dag[raw_run[0]] {
+            dag.get_qubits(inst.qubits)[0].into()
+        } else {
+            unreachable!("nodes in runs will always be op nodes")
+        };
+        if !dag.calibrations_empty() {
+            let mut has_calibration = false;
+            for node in &raw_run {
+                if dag.has_calibration_for_index(py, *node)? {
+                    has_calibration = true;
+                    break;
+                }
+            }
+            if has_calibration {
+                continue;
+            }
+        }
+        let basis_gates = match target {
+            Some(target) => Some(
+                target
+                    .operation_names_for_qargs(Some(&smallvec![qubit]))
+                    .unwrap(),
+            ),
+            None => {
+                let basis = basis_gates.as_ref();
+                basis.map(|basis| basis.iter().map(|x| x.as_str()).collect())
+            }
+        };
+
+        let mut target_basis_set: IndexSet<EulerBasis> = match target {
+            Some(_target) => EULER_BASIS_MAP
                 .iter()
-                .map(|node| {
-                    if let Some(err_map) = error_map {
-                        error *=
-                            compute_error_term(node.instruction.operation.name(), err_map, qubit)
+                .enumerate()
+                .filter_map(|(idx, gates)| {
+                    if !gates
+                        .iter()
+                        .all(|gate| basis_gates.as_ref().unwrap().contains(gate))
+                    {
+                        return None;
                     }
-                    node.instruction
-                        .operation
-                        .matrix(&node.instruction.params)
-                        .expect("No matrix defined for operation")
+                    let basis = EULER_BASIS_NAMES[idx];
+                    Some(basis)
                 })
-                .fold(
-                    [
-                        [Complex64::new(1., 0.), Complex64::new(0., 0.)],
-                        [Complex64::new(0., 0.), Complex64::new(1., 0.)],
-                    ],
-                    |mut operator, node| {
-                        matmul_1q(&mut operator, node);
-                        operator
-                    },
-                );
-            let old_error = if error_map.is_some() {
-                (1. - error, raw_run.len())
-            } else {
-                (error, raw_run.len())
-            };
-            let target_basis_vec: Vec<EulerBasis> = bases[index]
-                .iter()
-                .map(|basis| EulerBasis::__new__(basis).unwrap())
-                .collect();
-            unitary_to_gate_sequence_inner(
-                aview2(operator),
-                &target_basis_vec,
-                qubit,
-                error_map,
-                simplify,
-                atol,
-            )
-            .map(|out_seq| {
-                let new_error = compute_error_one_qubit_sequence(&out_seq, qubit, error_map);
-                (old_error, new_error, out_seq)
+                .collect(),
+            None => match &global_decomposers {
+                Some(bases) => bases
+                    .iter()
+                    .map(|basis| EulerBasis::__new__(basis).unwrap())
+                    .collect(),
+                None => match basis_gates {
+                    Some(ref gates) => EULER_BASIS_MAP
+                        .iter()
+                        .enumerate()
+                        .filter_map(|(idx, basis_gates)| {
+                            if !gates.iter().all(|gate| basis_gates.as_ref().contains(gate)) {
+                                return None;
+                            }
+                            let basis = EULER_BASIS_NAMES[idx];
+                            Some(basis)
+                        })
+                        .collect(),
+                    None => EULER_BASIS_NAMES.iter().copied().collect(),
+                },
+            },
+        };
+        if target_basis_set.contains(&EulerBasis::U3)
+            && target_basis_set.contains(&EulerBasis::U321)
+        {
+            target_basis_set.swap_remove(&EulerBasis::U3);
+        }
+        if target_basis_set.contains(&EulerBasis::ZSX)
+            && target_basis_set.contains(&EulerBasis::ZSXX)
+        {
+            target_basis_set.swap_remove(&EulerBasis::ZSX);
+        }
+        let target_basis_vec: Vec<EulerBasis> = target_basis_set.into_iter().collect();
+        let operator = raw_run
+            .iter()
+            .map(|node_index| {
+                let node = &dag.dag[*node_index];
+                if let NodeType::Operation(inst) = node {
+                    if let Some(target) = target {
+                        error *= compute_error_term_from_target(inst.op.name(), target, qubit);
+                    }
+                    inst.op.matrix(inst.params_view()).unwrap()
+                } else {
+                    unreachable!("Can only have op nodes here")
+                }
             })
-        })
-        .collect()
+            .fold(
+                [
+                    [Complex64::new(1., 0.), Complex64::new(0., 0.)],
+                    [Complex64::new(0., 0.), Complex64::new(1., 0.)],
+                ],
+                |mut operator, node| {
+                    matmul_1q(&mut operator, node);
+                    operator
+                },
+            );
+
+        let old_error = if target.is_some() {
+            (1. - error, raw_run.len())
+        } else {
+            (error, raw_run.len())
+        };
+        let sequence = unitary_to_gate_sequence_inner(
+            aview2(&operator),
+            &target_basis_vec,
+            qubit.0 as usize,
+            None,
+            true,
+            None,
+        );
+        let sequence = match sequence {
+            Some(seq) => seq,
+            None => continue,
+        };
+        let new_error = compute_error_from_target_one_qubit_sequence(&sequence, qubit, target);
+
+        let mut outside_basis = false;
+        if let Some(basis) = basis_gates {
+            for node in &raw_run {
+                if let NodeType::Operation(inst) = &dag.dag[*node] {
+                    if !basis.contains(inst.op.name()) {
+                        outside_basis = true;
+                        break;
+                    }
+                }
+            }
+        } else {
+            outside_basis = false;
+        }
+        if outside_basis
+            || new_error < old_error
+            || new_error.0.abs() < 1e-9 && old_error.0.abs() >= 1e-9
+        {
+            for gate in sequence.gates {
+                dag.insert_1q_on_incoming_qubit((gate.0, &gate.1), raw_run[0]);
+            }
+            dag.add_global_phase(py, &Param::Float(sequence.global_phase))?;
+            for node in raw_run {
+                dag.remove_op_node(node);
+            }
+        }
+    }
+    Ok(())
 }
 
+static EULER_BASIS_MAP: [&[&str]; 12] = [
+    &["u3"],
+    &["u3", "u2", "u1"],
+    &["u"],
+    &["p", "sx"],
+    &["u1", "rx"],
+    &["r"],
+    &["rz", "ry"],
+    &["rz", "rx"],
+    &["rz", "rx"],
+    &["rx", "ry"],
+    &["rz", "sx", "x"],
+    &["rz", "sx"],
+];
+static EULER_BASIS_NAMES: [EulerBasis; 12] = [
+    EulerBasis::U3,
+    EulerBasis::U321,
+    EulerBasis::U,
+    EulerBasis::PSX,
+    EulerBasis::U1X,
+    EulerBasis::RR,
+    EulerBasis::ZYZ,
+    EulerBasis::ZXZ,
+    EulerBasis::XZX,
+    EulerBasis::XYX,
+    EulerBasis::ZSXX,
+    EulerBasis::ZSX,
+];
+
 fn matmul_1q(operator: &mut [[Complex64; 2]; 2], other: Array2<Complex64>) {
     *operator = [
         [

crates/accelerate/src/nlayout.rs

@@ -14,6 +14,7 @@ use pyo3::prelude::*;
 use pyo3::types::PyList;
 
 use hashbrown::HashMap;
+use qiskit_circuit::Qubit;
 
 /// A newtype for the different categories of qubits used within layouts.  This is to enforce
 /// significantly more type safety when dealing with mixtures of physical and virtual qubits, as we
@@ -24,7 +25,7 @@ use hashbrown::HashMap;
 macro_rules! qubit_newtype {
     ($id: ident) => {
         #[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)]
-        pub struct $id(u32);
+        pub struct $id(pub u32);
 
         impl $id {
             #[inline]
@@ -72,6 +73,25 @@ impl PhysicalQubit {
         layout.phys_to_virt[self.index()]
     }
 }
+
+/// This is only safe in the context of a physical circuit during transpilation
+/// after the qubit indices of the circuit/dag circuit refer to the physical
+/// qubits (once we've run a layout pass and applied it).
+impl From<Qubit> for PhysicalQubit {
+    fn from(s: Qubit) -> PhysicalQubit {
+        PhysicalQubit::new(s.0)
+    }
+}
+
+/// This is only safe in the context of a physical circuit during transpilation
+/// after the qubit indices of the circuit/dag circuit refer to the physical
+/// qubits (once we've run a layout pass and applied it).
+impl From<&Qubit> for PhysicalQubit {
+    fn from(s: &Qubit) -> PhysicalQubit {
+        PhysicalQubit::new(s.0)
+    }
+}
+
 qubit_newtype!(VirtualQubit);
 impl VirtualQubit {
     /// Get the physical qubit that currently corresponds to this index of virtual qubit in the

crates/accelerate/src/target_transpiler/mod.rs

@@ -940,6 +940,17 @@ impl Target {
         });
     }
 
+    /// Get the error rate of a given instruction in the target
+    pub fn get_error(&self, name: &str, qargs: &[PhysicalQubit]) -> Option<f64> {
+        self.gate_map.get(name).and_then(|gate_props| {
+            let qargs_key: Qargs = qargs.iter().cloned().collect();
+            match gate_props.get(Some(&qargs_key)) {
+                Some(props) => props.as_ref().and_then(|inst_props| inst_props.error),
+                None => None,
+            }
+        })
+    }
+
     /// Get an iterator over the indices of all physical qubits of the target
     pub fn physical_qubits(&self) -> impl ExactSizeIterator<Item = usize> {
         0..self.num_qubits.unwrap_or_default()