Remove qiskit_circuit from QI matrix composition and Clifford (#15907)

This is the beginning of a thread of work to separate out dependence on
`qiskit-circuit` from `qiskit-quantum-information`, to better allow the
base types in the QI crate to be embedded in first-part circuit IR
types.

In general, the definitions of the quantum-information objects will
remain in `qiskit-quantum-information`.  First-party core-IR attributes
of them that related to circuit objects will move to `qiskit-circuit`,
while general circuit-level synthesis will move to `qiskit-synthesis`.

Author: jakelishman

crates/README.md

@@ -23,7 +23,10 @@ The intention is that (much) longer term, we might be wanting to expose more of
   This depends on `qiskit-circuit`, but is otherwise pretty standalone, and it's unlikely that other things will need to interact with it.
 * `qiskit-qasm3` is the Qiskit-specific side of the OpenQASM 3 importer.
   The actual parser lives at https://github.com/Qiskit/openqasm3_parser, and is its own set of Rust-only crates.
-* `qiskit-quantum-info` is the crate for quantum information.
+* `qiskit-quantum-info` contains definitions and operations on quantum-information-related objects,
+  without reference to circuit-level objects.  This is a base crate.  Interactions between
+  circuit-level objects and quantum-info objects are either in `qiskit-circuit` (if they are part of
+  core IR definitions) or `qiskit_synthesis::quantum_info` (if not).
 * `qiskit-synthesis` is the crate for synthesis functionality.
 * `qiskit-circuit-library` is the crate for circuit library functions. It contains constructors or other
   circuit functionality that builds on the core circuit data model defined in `qiskit-circuit`.

crates/quantum_info/src/clifford.rs

@@ -13,7 +13,6 @@ use std::fmt;
 
 use fixedbitset::FixedBitSet;
 use ndarray::{Array2, ArrayView2};
-use qiskit_circuit::Qubit;
 
 /// Symplectic matrix.
 pub struct SymplecticMatrix {
@@ -291,7 +290,7 @@ impl Clifford {
 
     /// Evolving the single-qubit Pauli-Z with Z on qubit qbit.
     /// Returns the evolved Pauli in the a sparse ZX format: (sign, z, x, indices).
-    pub fn get_inverse_z(&self, qbit: usize) -> (bool, Vec<bool>, Vec<bool>, Vec<Qubit>) {
+    pub fn get_inverse_z(&self, qbit: usize) -> (bool, Vec<bool>, Vec<bool>, Vec<usize>) {
         let mut z = Vec::with_capacity(self.num_qubits);
         let mut x = Vec::with_capacity(self.num_qubits);
         let mut indices = Vec::with_capacity(self.num_qubits);
@@ -304,7 +303,7 @@ impl Clifford {
             if z_bit || x_bit {
                 z.push(z_bit);
                 x.push(x_bit);
-                indices.push(Qubit::new(i));
+                indices.push(i);
                 if x_bit {
                     pauli_indices.push(i);
                 }

crates/quantum_info/src/unitary_compose.rs

@@ -13,7 +13,6 @@
 use ndarray::{Array, Array2, ArrayView, ArrayView2, IxDyn};
 use ndarray_einsum::*;
 use num_complex::{Complex, Complex64};
-use qiskit_circuit::Qubit;
 
 static LOWERCASE: [u8; 26] = [
     b'a', b'b', b'c', b'd', b'e', b'f', b'g', b'h', b'i', b'j', b'k', b'l', b'm', b'n', b'o', b'p',
@@ -30,13 +29,13 @@ static _UPPERCASE: [u8; 26] = [
 pub fn compose(
     gate_unitary: &ArrayView2<Complex64>,
     overall_unitary: &ArrayView2<Complex64>,
-    qubits: &[Qubit],
+    indices: &[usize],
     front: bool,
 ) -> Result<Array2<Complex64>, &'static str> {
     let gate_qubits = gate_unitary.shape()[0].ilog2() as usize;
 
     // Full composition of operators
-    if qubits.is_empty() {
+    if indices.is_empty() {
         if front {
             return Ok(gate_unitary.dot(overall_unitary));
         } else {
@@ -54,9 +53,9 @@ pub fn compose(
     //product, which is the last tensor wire index.
     let tensor = per_qubit_shaped(gate_unitary);
     let mat = per_qubit_shaped(overall_unitary);
-    let indices = qubits
+    let indices = indices
         .iter()
-        .map(|q| num_indices - 1 - q.index())
+        .map(|i| num_indices - 1 - i)
         .collect::<Vec<usize>>();
     let num_rows = usize::pow(2, num_indices as u32);
 
@@ -173,11 +172,11 @@ fn _ind(i: usize, reversed: bool) -> usize {
 pub fn matmul_2q(
     left: &ArrayView2<Complex64>,
     right: &ArrayView2<Complex64>,
-    qargs: &[Qubit],
+    indices: &[usize],
 ) -> Array2<Complex64> {
     let mut out = Array2::zeros((4, 4));
 
-    let rev = qargs[0].0 == 1;
+    let rev = indices[0] == 1;
     for i in 0..4usize {
         for j in 0..4usize {
             for k in 0..4usize {

crates/synthesis/src/matrix/sim.rs

@@ -48,6 +48,7 @@ pub fn sim_unitary_circuit(circuit: &CircuitData) -> Result<Array2<Complex64>, S
     // Product matrix holding the result
     let mut product_mat: Array2<Complex64> =
         Array2::<Complex64>::eye(2_usize.pow(num_qubits as u32)) * global_phase_exp;
+    let mut qubits = Vec::new();
 
     for inst in circuit.data() {
         if !circuit.get_cargs(inst.clbits).is_empty() {
@@ -61,13 +62,14 @@ pub fn sim_unitary_circuit(circuit: &CircuitData) -> Result<Array2<Complex64>, S
             continue;
         }
 
-        let qubits = circuit.get_qargs(inst.qubits);
+        qubits.clear();
+        qubits.extend(circuit.get_qargs(inst.qubits).iter().map(|q| q.index()));
 
         let mat = inst
             .try_matrix()
             .ok_or_else(|| format!("Cannot extract matrix for operation {:?}.", inst.op.name()))?;
 
-        product_mat = unitary_compose::compose(&product_mat.view(), &mat.view(), qubits, false)?;
+        product_mat = unitary_compose::compose(&product_mat.view(), &mat.view(), &qubits, false)?;
     }
 
     Ok(product_mat)

crates/transpiler/src/commutation_checker.rs

@@ -679,10 +679,13 @@ impl CommutationChecker {
             }
         }
 
-        let first_qarg: Vec<Qubit> = Vec::from_iter((0..first_qargs.len() as u32).map(Qubit));
-        let second_qarg: Vec<Qubit> = second_qargs.iter().map(|q| qarg[q]).collect();
+        let first_indices = Vec::from_iter(0..first_qargs.len());
+        let second_indices = second_qargs
+            .iter()
+            .map(|q| qarg[q].index())
+            .collect::<Vec<_>>();
 
-        if first_qarg.len() > second_qarg.len() {
+        if first_indices.len() > second_indices.len() {
             return Err(CommutationError::FirstInstructionTooLarge);
         };
         let first_mat = match try_matrix_with_definition(first_op, first_params, None) {
@@ -701,7 +704,7 @@ impl CommutationChecker {
         //  2. This code here expands the first op to match the second -- hence we always
         //     match the operator sizes.
         // This whole extension logic could be avoided since we know the second one is larger.
-        let extra_qarg2 = num_qubits - first_qarg.len() as u32;
+        let extra_qarg2 = num_qubits - first_indices.len() as u32;
         let first_mat = if extra_qarg2 > 0 {
             let id_op = Array2::<Complex64>::eye(usize::pow(2, extra_qarg2));
             kron(&id_op, &first_mat)
@@ -715,7 +718,7 @@ impl CommutationChecker {
         let op12 = match unitary_compose::compose(
             &first_mat.view(),
             &second_mat.view(),
-            &second_qarg,
+            &second_indices,
             false,
         ) {
             Ok(matrix) => matrix,
@@ -724,13 +727,13 @@ impl CommutationChecker {
         let op21 = match unitary_compose::compose(
             &first_mat.view(),
             &second_mat.view(),
-            &second_qarg,
+            &second_indices,
             true,
         ) {
             Ok(matrix) => matrix,
             Err(e) => return Err(CommutationError::EinsumError(e)),
         };
-        let (fid, phase) = gate_metrics::gate_fidelity(&op12.view(), &op21.view(), None);
+        let (fid, phase) = gate_metrics::gate_fidelity(&op12.view(), &op21.view());
 
         // we consider the gates as commuting if the process fidelity of
         // AB (BA)^\dagger is approximately the identity and there is no global phase difference

crates/transpiler/src/gate_metrics.rs

@@ -12,7 +12,7 @@
 
 use ndarray::ArrayView2;
 use num_complex::{Complex64, ComplexFloat};
-use qiskit_circuit::{Qubit, operations::StandardGate};
+use qiskit_circuit::operations::StandardGate;
 
 use qiskit_quantum_info::unitary_compose;
 
@@ -46,19 +46,13 @@ pub fn rotation_trace_and_dim(rotation: StandardGate, angle: f64) -> Option<(Com
     Some((trace_over_dim, dim))
 }
 
-pub fn gate_fidelity(
-    left: &ArrayView2<Complex64>,
-    right: &ArrayView2<Complex64>,
-    qargs: Option<&[Qubit]>,
-) -> (f64, f64) {
+pub fn gate_fidelity(left: &ArrayView2<Complex64>, right: &ArrayView2<Complex64>) -> (f64, f64) {
     let dim = left.nrows();
 
     let left = left.t().mapv(|el| el.conj());
     let product = match dim {
         2 => unitary_compose::matmul_1q(&left.view(), right),
-        4 => {
-            unitary_compose::matmul_2q(&left.view(), right, qargs.unwrap_or(&[Qubit(0), Qubit(1)]))
-        }
+        4 => unitary_compose::matmul_2q(&left.view(), right, &[0, 1]),
         _ => left.dot(right),
     };
     let trace = product.diag().sum();

crates/transpiler/src/passes/litinski_transformation.rs

@@ -86,6 +86,8 @@ pub fn run_litinski_transformation(
     let num_qubits = dag.num_qubits();
     let mut clifford = Clifford::identity(num_qubits);
 
+    let mut qargs = Vec::new();
+
     // Keep track of the update to the global phase (produced when converting T/Tdg gates
     // to RZ-rotations).
     let mut global_phase_update = Param::Float(0.);
@@ -251,6 +253,8 @@ pub fn run_litinski_transformation(
                     // where signs `true` and `false` correspond to coefficients `-1` and `+1` respectively.
                     let (sign, z, x, indices) =
                         clifford.get_inverse_z(dag.get_qargs(inst.qubits)[0].index());
+                    qargs.clear();
+                    qargs.extend(indices.iter().map(|i| Qubit::new(*i)));
 
                     // In the legacy path, we add PauliEvolutionGate as rotation gates, otherwise
                     // we add PauliProductRotation. The new path should not call Python at any
@@ -279,7 +283,7 @@ pub fn run_litinski_transformation(
                                 .get_bound(py)
                                 .call1((obs, time.clone()))?;
                             Ok(PyOperationTypes::Gate(PyInstruction {
-                                qubits: indices.len() as u32,
+                                qubits: qargs.len() as u32,
                                 clbits: 0,
                                 params: 1,
                                 op_name: "PauliEvolution".to_string(),
@@ -291,7 +295,7 @@ pub fn run_litinski_transformation(
 
                     new_dag.apply_operation_back(
                         packed_op,
-                        &indices,
+                        &qargs,
                         &[],
                         Some(Parameters::Params(smallvec![param])),
                         None,
@@ -304,13 +308,16 @@ pub fn run_litinski_transformation(
                     // where signs `true` and `false` correspond to coefficients `-1` and `+1` respectively.
                     let (sign, z, x, indices) =
                         clifford.get_inverse_z(dag.get_qargs(inst.qubits)[0].index());
+                    qargs.clear();
+                    qargs.extend(indices.iter().map(|i| Qubit::new(*i)));
+
                     let ppm = PauliProductMeasurement { z, x, neg: sign };
 
                     let ppm_clbits = dag.get_cargs(inst.clbits);
 
                     new_dag.apply_operation_back(
                         PauliBased::PauliProductMeasurement(ppm).into(),
-                        &indices,
+                        &qargs,
                         ppm_clbits,
                         None,
                         None,