Rewrite default UnitarySynthesis to cache decomposers (Qiskit#15492)

* Rewrite default `UnitarySynthesis` to cache decomposers

This is a complete rewrite of the default `UnitarySynthesis` plugin to
vastly improve the efficiency of 2q decomposition.

This is good for approximately a 2x runtime improvement in the
`UnitarySynthesis` pass, since before this commit, we were
reconstructing each relevant 2q decomposer each time we had a new matrix
to decompose.  However, constructing a 2q KAK decomposer is about as
expensive as using a constructed KAK decomposer on a single matrix.
This commit caches the available decomposers for each encountered pair
of qubits, and allows the cache to be persisted between calls to
`run_unitary_synthesis`.

Before the move of the default unitary synthesis plugin to Rust, we
effectively had decomposer caching for the "loose constraint" (basis
gates + coupling map) hardware description, since we just chose a single
decomposer on initialisation and used it throughout.  The `Target` form
in Python space cached at the `qargs` level, which meant that multiple
unitaries on the same qargs pair would use the same set of decomposers,
but each qargs pair would be calculated separately on first access,
still (generally) leading to multiple constructions of the same
decomposer.  Both of these types of caching were lost in the move to
Rust, but the effects were largely masked by the total runtime still
being drastically better than the Python-space versions.

This new form reinstates all the previous caching, and additionally
caches at the level of individual decomposer construction as well (by
caching the arguments used to construct a decomposer), so that (mostly)
homogeneous `Target`s will re-use the same decomposer whenever it is
valid on more than one 2q link.

* Give synthesis/state structs explicit names

There's lots of types of "synthesis" in Qiskit...

* Correct copyright years

* Avoid magic numbers in `ApproximationDegree::is_approximate`

Co-authored-by: Alexander Ivrii <alexi@il.ibm.com>

---------

Co-authored-by: Alexander Ivrii <alexi@il.ibm.com>

Author: jakelishman

crates/cext/src/transpiler/passes/unitary_synthesis.rs

@@ -10,14 +10,15 @@
 // copyright notice, and modified files need to carry a notice indicating
 // that they have been altered from the originals.
 
-use hashbrown::HashSet;
-
 use crate::pointers::{const_ptr_as_ref, mut_ptr_as_ref};
 
+use qiskit_circuit::PhysicalQubit;
 use qiskit_circuit::circuit_data::CircuitData;
 use qiskit_circuit::converters::dag_to_circuit;
 use qiskit_circuit::dag_circuit::DAGCircuit;
-use qiskit_transpiler::passes::run_unitary_synthesis;
+use qiskit_transpiler::passes::{
+    UnitarySynthesisConfig, UnitarySynthesisState, run_unitary_synthesis,
+};
 use qiskit_transpiler::target::Target;
 
 /// @ingroup QkTranspilerPasses
@@ -75,7 +76,7 @@ pub unsafe extern "C" fn qk_transpiler_pass_standalone_unitary_synthesis(
     // SAFETY: Per documentation, the pointer is non-null and aligned.
     let circuit = unsafe { mut_ptr_as_ref(circuit) };
     let target = unsafe { const_ptr_as_ref(target) };
-    let mut dag = match DAGCircuit::from_circuit_data(circuit, false, None, None, None, None) {
+    let dag = match DAGCircuit::from_circuit_data(circuit, false, None, None, None, None) {
         Ok(dag) => dag,
         Err(e) => panic!("{}", e),
     };
@@ -84,19 +85,21 @@ pub unsafe extern "C" fn qk_transpiler_pass_standalone_unitary_synthesis(
     } else {
         Some(approximation_degree)
     };
-    let qubit_indices = (0..dag.num_qubits()).collect();
+    let physical_qubits = (0..dag.num_qubits() as u32)
+        .map(PhysicalQubit::new)
+        .collect::<Vec<_>>();
+    let mut synthesis_state = UnitarySynthesisState::new(UnitarySynthesisConfig {
+        approximation_degree,
+        run_python_decomposers: false,
+        ..Default::default()
+    });
     let out_dag = match run_unitary_synthesis(
-        &mut dag,
-        qubit_indices,
+        &dag,
+        &["unitary".to_string()].into_iter().collect(),
         min_qubits,
-        Some(target),
-        HashSet::new(),
-        ["unitary".to_string()].into_iter().collect(),
-        HashSet::new(),
-        approximation_degree,
-        None,
-        None,
-        false,
+        &physical_qubits,
+        &mut synthesis_state,
+        target.into(),
     ) {
         Ok(dag) => dag,
         Err(e) => panic!("{}", e),

crates/cext/src/transpiler/transpile_function.rs

@@ -15,6 +15,7 @@ use std::ffi::c_char;
 use qiskit_circuit::circuit_data::CircuitData;
 use qiskit_circuit::dag_circuit::DAGCircuit;
 use qiskit_transpiler::commutation_checker::get_standard_commutation_checker;
+use qiskit_transpiler::passes::{UnitarySynthesisConfig, UnitarySynthesisState};
 use qiskit_transpiler::standard_equivalence_library::generate_standard_equivalence_library;
 use qiskit_transpiler::target::Target;
 use qiskit_transpiler::transpile;
@@ -164,13 +165,19 @@ pub unsafe extern "C" fn qk_transpile_stage_init(
         }
         Some(options.approximation_degree)
     };
+    let mut synthesis_state = UnitarySynthesisState::new(UnitarySynthesisConfig {
+        approximation_degree,
+        run_python_decomposers: false,
+        ..Default::default()
+    });
     let mut commutation_checker = get_standard_commutation_checker();
 
     match init_stage(
         dag,
         target,
         options.optimization_level.into(),
         approximation_degree,
+        &mut synthesis_state,
         &mut out_layout,
         &mut commutation_checker,
     ) {
@@ -399,6 +406,11 @@ pub unsafe extern "C" fn qk_transpile_stage_optimization(
         }
         Some(options.approximation_degree)
     };
+    let mut synthesis_state = UnitarySynthesisState::new(UnitarySynthesisConfig {
+        approximation_degree,
+        run_python_decomposers: false,
+        ..Default::default()
+    });
     let mut equiv_lib = generate_standard_equivalence_library();
     let mut commutation_checker = get_standard_commutation_checker();
 
@@ -407,6 +419,7 @@ pub unsafe extern "C" fn qk_transpile_stage_optimization(
         target,
         options.optimization_level.into(),
         approximation_degree,
+        &mut synthesis_state,
         &mut commutation_checker,
         &mut equiv_lib,
     ) {
@@ -500,9 +513,14 @@ pub unsafe extern "C" fn qk_transpile_stage_translation(
         }
         Some(options.approximation_degree)
     };
+    let mut synthesis_state = UnitarySynthesisState::new(UnitarySynthesisConfig {
+        approximation_degree,
+        run_python_decomposers: false,
+        ..Default::default()
+    });
     let mut equiv_lib = generate_standard_equivalence_library();
 
-    match translation_stage(dag, target, approximation_degree, &mut equiv_lib) {
+    match translation_stage(dag, target, &mut synthesis_state, &mut equiv_lib) {
         Ok(_) => ExitCode::Success,
         Err(e) => {
             if !error.is_null() {

crates/circuit/src/dag_circuit.rs

@@ -8634,11 +8634,35 @@ impl DAGCircuitBuilder {
     pub fn insert_qargs(&mut self, qargs: &[Qubit]) -> Interned<[Qubit]> {
         self.dag.qargs_interner.insert(qargs)
     }
+    /// Merge the `qargs` in a different [Interner] into this DAG, remapping the qubits.
+    ///
+    /// This is useful for simplifying the direct mapping of [PackedInstruction]s from one DAG to
+    /// another, like in substitution methods, or rebuilding a new DAG out of a lot of smaller ones.
+    /// See [Interner::merge_map_slice] for more information on the mapping function.
+    pub fn merge_qargs(
+        &mut self,
+        other: &Interner<[Qubit]>,
+        map_fn: impl FnMut(&Qubit) -> Option<Qubit>,
+    ) -> InternedMap<[Qubit]> {
+        self.dag.merge_qargs(other, map_fn)
+    }
 
     /// Packs qargs into the circuit.
     pub fn insert_cargs(&mut self, cargs: &[Clbit]) -> Interned<[Clbit]> {
         self.dag.cargs_interner.insert(cargs)
     }
+    /// Merge the `cargs` in a different [Interner] into this DAG, remapping the clbits.
+    ///
+    /// This is useful for simplifying the direct mapping of [PackedInstruction]s from one DAG to
+    /// another, like in substitution methods, or rebuilding a new DAG out of a lot of smaller ones.
+    /// See [Interner::merge_map_slice] for more information on the mapping function.
+    pub fn merge_cargs(
+        &mut self,
+        other: &Interner<[Clbit]>,
+        map_fn: impl FnMut(&Clbit) -> Option<Clbit>,
+    ) -> InternedMap<[Clbit]> {
+        self.dag.merge_cargs(other, map_fn)
+    }
 
     /// Adds a new value to the global phase of the inner [DAGCircuit].
     pub fn add_global_phase(&mut self, param: &Param) -> PyResult<()> {

crates/circuit/src/operations.rs

@@ -93,6 +93,13 @@ impl<'a, 'py> FromPyObject<'a, 'py> for Param {
 }
 
 impl Param {
+    /// Get the float value, if one is stored.
+    pub fn try_float(&self) -> Option<f64> {
+        match self {
+            Self::Float(f) => Some(*f),
+            _ => None,
+        }
+    }
     pub fn eq(&self, other: &Param) -> PyResult<bool> {
         match [self, other] {
             [Self::Float(a), Self::Float(b)] => Ok(a == b),

crates/circuit/src/packed_instruction.rs

@@ -26,7 +26,7 @@ use crate::operations::{
 use crate::{Block, Clbit, Qubit};
 use hashbrown::HashMap;
 use nalgebra::Matrix2;
-use ndarray::Array2;
+use ndarray::{Array2, CowArray, Ix2};
 use num_complex::Complex64;
 use pyo3::prelude::*;
 use pyo3::types::{PyDict, PyType};
@@ -428,6 +428,17 @@ impl PackedOperation {
         }
     }
 
+    /// Does this [PackedOperation] represent an explicit gate?
+    ///
+    /// This can be either a [StandardGate] or a [PyGate].
+    #[inline]
+    pub fn is_gate(&self) -> bool {
+        matches!(
+            self.discriminant(),
+            PackedOperationType::StandardGate | PackedOperationType::PyGate
+        )
+    }
+
     /// Create a `PackedOperation` from a `StandardGate`.
     #[inline]
     pub fn from_standard_gate(standard: StandardGate) -> Self {
@@ -824,6 +835,10 @@ impl PackedInstruction {
         Ok(())
     }
 
+    /// Extract an owned `ndarray` matrix from this instruction, if available.
+    ///
+    /// The returned value is always owned.  If you may be able to handle a read-only reference, see
+    /// [try_cow_array] instead.
     pub fn try_matrix(&self) -> Option<Array2<Complex64>> {
         match self.op.view() {
             OperationRef::StandardGate(g) => g.matrix(self.params_view()),
@@ -833,6 +848,19 @@ impl PackedInstruction {
         }
     }
 
+    /// Extract an `ndarray` matrix from this instruciton, if available.
+    ///
+    /// The returned value will preferentially be a view, if the matrix already exists (e.g. for
+    /// `Unitary`).
+    pub fn try_cow_array(&self) -> Option<CowArray<Complex64, Ix2>> {
+        match self.op.view() {
+            OperationRef::StandardGate(g) => g.matrix(self.params_view()).map(CowArray::from),
+            OperationRef::Gate(g) => g.matrix().map(CowArray::from),
+            OperationRef::Unitary(u) => Some(CowArray::from(u.matrix_view())),
+            _ => None,
+        }
+    }
+
     /// Returns a static matrix for 1-qubit gates. Will return `None` when the gate is not 1-qubit.
     #[inline]
     pub fn try_matrix_as_static_1q(&self) -> Option<[[Complex64; 2]; 2]> {

crates/synthesis/src/discrete_basis/basic_approximations.rs

@@ -355,7 +355,7 @@ impl Point for BasicPoint {
 ///
 /// This struct allows to construct a tree of basic approximations and to query the closest
 /// sequence given an target sequence (or SO(3) matrix).
-#[derive(Debug)]
+#[derive(Clone, Debug)]
 pub struct BasicApproximations {
     /// All points as flattened SO(3) matrix stored in a R* tree. This does not include the
     /// sequence of gates, see ``approximations``.

crates/synthesis/src/discrete_basis/solovay_kitaev.rs

@@ -33,6 +33,7 @@ use super::math::{self, group_commutator_decomposition};
 /// This generates the basic approximation set once as R-tree and re-uses it for
 /// each queried decomposition.
 #[pyclass]
+#[derive(Clone, Debug)]
 pub struct SolovayKitaevSynthesis {
     /// The set of basic approximations.
     basic_approximations: BasicApproximations,

crates/synthesis/src/euler_one_qubit_decomposer.rs

@@ -605,7 +605,7 @@ pub static EULER_BASIS_NAMES: [EulerBasis; EULER_BASIS_SIZE] = [
 ];
 
 /// A structure containing a set of supported `EulerBasis` for running 1q synthesis
-#[derive(Debug, Clone)]
+#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
 pub struct EulerBasisSet {
     basis: [bool; EULER_BASIS_SIZE],
     initialized: bool,
@@ -620,6 +620,24 @@ impl EulerBasisSet {
         }
     }
 
+    /// Get the set of supported bases given a function that marks whether each basis gate (by name)
+    /// is supported.
+    ///
+    /// The `is_supported` function may be called more than once for the same gate name.
+    pub fn from_support(mut is_supported: impl FnMut(&str) -> bool) -> Self {
+        let mut out = Self {
+            basis: EULER_BASES.map(|basis| basis.iter().all(|gate| is_supported(gate))),
+            initialized: true,
+        };
+        if out.basis_supported(EulerBasis::U321) && out.basis_supported(EulerBasis::U3) {
+            out.remove(EulerBasis::U3);
+        }
+        if out.basis_supported(EulerBasis::ZSXX) && out.basis_supported(EulerBasis::ZSX) {
+            out.remove(EulerBasis::ZSX);
+        }
+        out
+    }
+
     /// Return true if this has been initialized any basis is supported
     pub fn initialized(&self) -> bool {
         self.initialized

crates/synthesis/src/two_qubit_decompose.rs

@@ -1279,6 +1279,13 @@ impl TwoQubitGateSequence {
             global_phase: 0.,
         }
     }
+    /// Create this sequence from the consituent parts.
+    pub fn from_sequence(gates: TwoQubitSequenceVec, global_phase: f64) -> Self {
+        Self {
+            gates,
+            global_phase,
+        }
+    }
 }
 
 impl Default for TwoQubitGateSequence {

crates/transpiler/src/passes/mod.rs

@@ -92,7 +92,9 @@ pub use remove_diagonal_gates_before_measure::{
 pub use remove_identity_equiv::{remove_identity_equiv_mod, run_remove_identity_equiv};
 pub use split_2q_unitaries::{run_split_2q_unitaries, split_2q_unitaries_mod};
 pub use substitute_pi4_rotations::{py_run_substitute_pi4_rotations, substitute_pi4_rotations_mod};
-pub use unitary_synthesis::{run_unitary_synthesis, unitary_synthesis_mod};
+pub use unitary_synthesis::{
+    UnitarySynthesisConfig, UnitarySynthesisState, run_unitary_synthesis, unitary_synthesis_mod,
+};
 pub use unroll_3q_or_more::{run_unroll_3q_or_more, unroll_3q_or_more_mod};
 pub use vf2::{error_map_mod, vf2_layout_mod, vf2_layout_pass_average, vf2_layout_pass_exact};
 pub use wrap_angles::{run_wrap_angles, wrap_angles_mod};