Merge branch 'main' into schedule-follow-ip

.github/workflows/coverage.yml

@@ -71,7 +71,7 @@ jobs:
           lcov --add-tracefile python.info --add-tracefile rust.info --output-file coveralls.info
 
       - name: Coveralls
-        uses: coverallsapp/github-action@master
+        uses: coverallsapp/github-action@v2.3.0
         with:
           github-token: ${{ secrets.GITHUB_TOKEN }}
           path-to-lcov: coveralls.info

Cargo.toml

@@ -15,7 +15,7 @@ license = "Apache-2.0"
 # Each crate can add on specific features freely as it inherits.
 [workspace.dependencies]
 bytemuck = "1.18"
-indexmap.version = "2.5.0"
+indexmap.version = "2.6.0"
 hashbrown.version = "0.14.5"
 num-bigint = "0.4"
 num-complex = "0.4"

crates/accelerate/src/circuit_library/mod.rs

@@ -14,9 +14,11 @@ use pyo3::prelude::*;
 
 mod entanglement;
 mod pauli_feature_map;
+mod quantum_volume;
 
 pub fn circuit_library(m: &Bound<PyModule>) -> PyResult<()> {
     m.add_wrapped(wrap_pyfunction!(pauli_feature_map::pauli_feature_map))?;
     m.add_wrapped(wrap_pyfunction!(entanglement::get_entangler_map))?;
+    m.add_wrapped(wrap_pyfunction!(quantum_volume::quantum_volume))?;
     Ok(())
 }

crates/accelerate/src/circuit_library/pauli_feature_map.rs

@@ -138,7 +138,7 @@ fn pauli_evolution(
 ///     insert_barriers: Whether to insert barriers in between the Hadamard and evolution layers.
 ///     data_map_func: An accumulation function that takes as input a vector of parameters the
 ///         current gate acts on and returns a scalar.
-///     
+///
 /// Returns:
 ///     The ``CircuitData`` to construct the Pauli feature map.
 #[pyfunction]
@@ -207,7 +207,13 @@ pub fn pauli_feature_map(
         }
     }
 
-    CircuitData::from_packed_operations(py, feature_dimension, 0, packed_insts, Param::Float(0.0))
+    CircuitData::from_packed_operations(
+        py,
+        feature_dimension,
+        0,
+        packed_insts.into_iter().map(Ok),
+        Param::Float(0.0),
+    )
 }
 
 fn _get_h_layer(feature_dimension: u32) -> impl Iterator<Item = Instruction> {

crates/accelerate/src/circuit_library/quantum_volume.rs

@@ -0,0 +1,176 @@
+// 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 pyo3::prelude::*;
+
+use crate::getenv_use_multiple_threads;
+use faer_ext::{IntoFaerComplex, IntoNdarrayComplex};
+use ndarray::prelude::*;
+use num_complex::Complex64;
+use numpy::IntoPyArray;
+use rand::prelude::*;
+use rand_distr::StandardNormal;
+use rand_pcg::Pcg64Mcg;
+use rayon::prelude::*;
+
+use qiskit_circuit::circuit_data::CircuitData;
+use qiskit_circuit::imports::UNITARY_GATE;
+use qiskit_circuit::operations::Param;
+use qiskit_circuit::operations::PyInstruction;
+use qiskit_circuit::packed_instruction::PackedOperation;
+use qiskit_circuit::{Clbit, Qubit};
+use smallvec::{smallvec, SmallVec};
+
+type Instruction = (
+    PackedOperation,
+    SmallVec<[Param; 3]>,
+    Vec<Qubit>,
+    Vec<Clbit>,
+);
+
+#[inline(always)]
+fn random_complex(rng: &mut Pcg64Mcg) -> Complex64 {
+    Complex64::new(rng.sample(StandardNormal), rng.sample(StandardNormal))
+        * std::f64::consts::FRAC_1_SQRT_2
+}
+
+// This function's implementation was modeled off of the algorithm used in the
+// `scipy.stats.unitary_group.rvs()` function defined here:
+//
+// https://github.com/scipy/scipy/blob/v1.14.1/scipy/stats/_multivariate.py#L4224-L4256
+#[inline]
+fn random_unitaries(seed: u64, size: usize) -> impl Iterator<Item = Array2<Complex64>> {
+    let mut rng = Pcg64Mcg::seed_from_u64(seed);
+
+    (0..size).map(move |_| {
+        let raw_numbers: [[Complex64; 4]; 4] = [
+            [
+                random_complex(&mut rng),
+                random_complex(&mut rng),
+                random_complex(&mut rng),
+                random_complex(&mut rng),
+            ],
+            [
+                random_complex(&mut rng),
+                random_complex(&mut rng),
+                random_complex(&mut rng),
+                random_complex(&mut rng),
+            ],
+            [
+                random_complex(&mut rng),
+                random_complex(&mut rng),
+                random_complex(&mut rng),
+                random_complex(&mut rng),
+            ],
+            [
+                random_complex(&mut rng),
+                random_complex(&mut rng),
+                random_complex(&mut rng),
+                random_complex(&mut rng),
+            ],
+        ];
+
+        let qr = aview2(&raw_numbers).into_faer_complex().qr();
+        let r = qr.compute_r();
+        let diag: [Complex64; 4] = [
+            r[(0, 0)].to_num_complex() / r[(0, 0)].abs(),
+            r[(1, 1)].to_num_complex() / r[(1, 1)].abs(),
+            r[(2, 2)].to_num_complex() / r[(2, 2)].abs(),
+            r[(3, 3)].to_num_complex() / r[(3, 3)].abs(),
+        ];
+        let mut q = qr.compute_q().as_ref().into_ndarray_complex().to_owned();
+        q.axis_iter_mut(Axis(0)).for_each(|mut row| {
+            row.iter_mut()
+                .enumerate()
+                .for_each(|(index, val)| *val *= diag[index])
+        });
+        q
+    })
+}
+
+const UNITARY_PER_SEED: usize = 50;
+
+#[pyfunction]
+pub fn quantum_volume(
+    py: Python,
+    num_qubits: u32,
+    depth: usize,
+    seed: Option<u64>,
+) -> PyResult<CircuitData> {
+    let width = num_qubits as usize / 2;
+    let num_unitaries = width * depth;
+    let mut permutation: Vec<Qubit> = (0..num_qubits).map(Qubit).collect();
+
+    let mut build_instruction = |(unitary_index, unitary_array): (usize, Array2<Complex64>),
+                                 rng: &mut Pcg64Mcg|
+     -> PyResult<Instruction> {
+        let layer_index = unitary_index % width;
+        if layer_index == 0 {
+            permutation.shuffle(rng);
+        }
+        let unitary = unitary_array.into_pyarray_bound(py);
+        let unitary_gate = UNITARY_GATE
+            .get_bound(py)
+            .call1((unitary.clone(), py.None(), false))?;
+        let instruction = PyInstruction {
+            qubits: 2,
+            clbits: 0,
+            params: 1,
+            op_name: "unitary".to_string(),
+            control_flow: false,
+            instruction: unitary_gate.unbind(),
+        };
+        let qubit = layer_index * 2;
+        Ok((
+            PackedOperation::from_instruction(Box::new(instruction)),
+            smallvec![Param::Obj(unitary.unbind().into())],
+            vec![permutation[qubit], permutation[qubit + 1]],
+            vec![],
+        ))
+    };
+
+    let mut per_thread = num_unitaries / UNITARY_PER_SEED;
+    if per_thread == 0 {
+        per_thread = 10;
+    }
+    let mut outer_rng = match seed {
+        Some(seed) => Pcg64Mcg::seed_from_u64(seed),
+        None => Pcg64Mcg::from_entropy(),
+    };
+    let seed_vec: Vec<u64> = rand::distributions::Standard
+        .sample_iter(&mut outer_rng)
+        .take(num_unitaries)
+        .collect();
+
+    let unitaries: Vec<Array2<Complex64>> = if getenv_use_multiple_threads() && num_unitaries > 200
+    {
+        seed_vec
+            .par_chunks(per_thread)
+            .flat_map_iter(|seeds| random_unitaries(seeds[0], seeds.len()))
+            .collect()
+    } else {
+        seed_vec
+            .chunks(per_thread)
+            .flat_map(|seeds| random_unitaries(seeds[0], seeds.len()))
+            .collect()
+    };
+    CircuitData::from_packed_operations(
+        py,
+        num_qubits,
+        0,
+        unitaries
+            .into_iter()
+            .enumerate()
+            .map(|x| build_instruction(x, &mut outer_rng)),
+        Param::Float(0.),
+    )
+}

crates/accelerate/src/elide_permutations.rs

@@ -0,0 +1,109 @@
+// 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 numpy::PyReadonlyArray1;
+use pyo3::prelude::*;
+
+use qiskit_circuit::dag_circuit::{DAGCircuit, NodeType};
+use qiskit_circuit::operations::{Operation, Param};
+use qiskit_circuit::Qubit;
+
+/// Run the ElidePermutations pass on `dag`.
+/// Args:
+///     dag (DAGCircuit): the DAG to be optimized.
+/// Returns:
+///     An `Option`: the value of `None` indicates that no optimization was
+///     performed and the original `dag` should be used, otherwise it's a
+///     tuple consisting of the optimized DAG and the induced qubit permutation.
+#[pyfunction]
+fn run(py: Python, dag: &mut DAGCircuit) -> PyResult<Option<(DAGCircuit, Vec<usize>)>> {
+    let permutation_gate_names = ["swap".to_string(), "permutation".to_string()];
+    let op_counts = dag.count_ops(py, false)?;
+    if !permutation_gate_names
+        .iter()
+        .any(|name| op_counts.contains_key(name))
+    {
+        return Ok(None);
+    }
+    let mut mapping: Vec<usize> = (0..dag.num_qubits()).collect();
+
+    // note that DAGCircuit::copy_empty_like clones the interners
+    let mut new_dag = dag.copy_empty_like(py, "alike")?;
+    for node_index in dag.topological_op_nodes()? {
+        if let NodeType::Operation(inst) = &dag.dag()[node_index] {
+            match (inst.op.name(), inst.condition()) {
+                ("swap", None) => {
+                    let qargs = dag.get_qargs(inst.qubits);
+                    let index0 = qargs[0].0 as usize;
+                    let index1 = qargs[1].0 as usize;
+                    mapping.swap(index0, index1);
+                }
+                ("permutation", None) => {
+                    if let Param::Obj(ref pyobj) = inst.params.as_ref().unwrap()[0] {
+                        let pyarray: PyReadonlyArray1<i32> = pyobj.extract(py)?;
+                        let pattern = pyarray.as_array();
+
+                        let qindices: Vec<usize> = dag
+                            .get_qargs(inst.qubits)
+                            .iter()
+                            .map(|q| q.0 as usize)
+                            .collect();
+
+                        let remapped_qindices: Vec<usize> = (0..qindices.len())
+                            .map(|i| pattern[i])
+                            .map(|i| qindices[i as usize])
+                            .collect();
+
+                        qindices
+                            .iter()
+                            .zip(remapped_qindices.iter())
+                            .for_each(|(old, new)| {
+                                mapping[*old] = *new;
+                            });
+                    } else {
+                        unreachable!();
+                    }
+                }
+                _ => {
+                    // General instruction
+                    let qargs = dag.get_qargs(inst.qubits);
+                    let cargs = dag.get_cargs(inst.clbits);
+                    let mapped_qargs: Vec<Qubit> = qargs
+                        .iter()
+                        .map(|q| q.0 as usize)
+                        .map(|q| mapping[q])
+                        .map(|q| Qubit(q.try_into().unwrap()))
+                        .collect();
+
+                    new_dag.apply_operation_back(
+                        py,
+                        inst.op.clone(),
+                        &mapped_qargs,
+                        cargs,
+                        inst.params.as_deref().cloned(),
+                        inst.extra_attrs.clone(),
+                        #[cfg(feature = "cache_pygates")]
+                        None,
+                    )?;
+                }
+            }
+        } else {
+            unreachable!();
+        }
+    }
+    Ok(Some((new_dag, mapping)))
+}
+
+pub fn elide_permutations(m: &Bound<PyModule>) -> PyResult<()> {
+    m.add_wrapped(wrap_pyfunction!(run))?;
+    Ok(())
+}

crates/accelerate/src/lib.rs

@@ -24,6 +24,7 @@ pub mod commutation_checker;
 pub mod convert_2q_block_matrix;
 pub mod dense_layout;
 pub mod edge_collections;
+pub mod elide_permutations;
 pub mod equivalence;
 pub mod error_map;
 pub mod euler_one_qubit_decomposer;

crates/accelerate/src/synthesis/mod.rs

@@ -13,6 +13,7 @@
 mod clifford;
 pub mod linear;
 pub mod linear_phase;
+mod multi_controlled;
 mod permutation;
 
 use pyo3::prelude::*;
@@ -34,5 +35,9 @@ pub fn synthesis(m: &Bound<PyModule>) -> PyResult<()> {
     clifford::clifford(&clifford_mod)?;
     m.add_submodule(&clifford_mod)?;
 
+    let mc_mod = PyModule::new_bound(m.py(), "multi_controlled")?;
+    multi_controlled::multi_controlled(&mc_mod)?;
+    m.add_submodule(&mc_mod)?;
+
     Ok(())
 }