Oxidize synth_permutation_acg

crates/accelerate/src/permutation.rs

@@ -14,8 +14,13 @@ use ndarray::{Array1, ArrayView1};
 use numpy::PyArrayLike1;
 use pyo3::exceptions::PyValueError;
 use pyo3::prelude::*;
+use smallvec::smallvec;
 use std::vec::Vec;
 
+use qiskit_circuit::circuit_data::CircuitData;
+use qiskit_circuit::operations::{Param, StandardGate};
+use qiskit_circuit::Qubit;
+
 fn validate_permutation(pattern: &ArrayView1<i64>) -> PyResult<()> {
     let n = pattern.len();
     let mut seen: Vec<bool> = vec![false; n];
@@ -77,6 +82,86 @@ fn get_ordered_swap(pattern: &ArrayView1<i64>) -> Vec<(i64, i64)> {
     swaps
 }
 
+fn pattern_to_cycles(pattern: &ArrayView1<i64>, invert_order: &bool) -> Vec<Vec<usize>> {
+    // vector keeping track of which elements in the permutation pattern have been visited
+    let mut explored: Vec<bool> = vec![false; pattern.len()];
+
+    // vector to store the cycles
+    let mut cycles: Vec<Vec<usize>> = Vec::new();
+
+    // cast the input pattern in terms of integers to usize, such that it can be used as index
+    // also invert the bit ordering if ``invert_order`` is true
+    let permutation: Array1<usize> = if *invert_order {
+        invert(pattern) // implies cast to usize
+    } else {
+        pattern.mapv(|x| x as usize)
+    };
+
+    for mut i in permutation.clone() {
+        let mut cycle: Vec<usize> = Vec::new();
+
+        // follow the cycle until we reached an entry we saw before
+        while !explored[i] {
+            cycle.push(i);
+            explored[i] = true;
+            i = permutation[i];
+        }
+        // cycles must have more than 1 element
+        if cycle.len() > 1 {
+            cycles.push(cycle);
+        }
+    }
+
+    cycles
+}
+
+/// Given a disjoint cycle decomposition, decomposes every cycle into a SWAP
+fn decompose_cycles(cycles: &Vec<Vec<usize>>) -> Vec<(usize, usize)> {
+    let mut swaps: Vec<(usize, usize)> = Vec::new();
+
+    for cycle in cycles {
+        let length = cycle.len();
+
+        if length > 2 {
+            // handle first element separately, which accesses the last element
+            swaps.push((cycle[length - 1], cycle[length - 3]));
+            for i in 1..(length - 1) / 2 {
+                swaps.push((cycle[i - 1], cycle[length - 3 - i]));
+            }
+        }
+
+        // no size check needed, cycles always have at least 2 elements
+        swaps.push((cycle[length - 1], cycle[length - 2]));
+        for i in 1..length / 2 {
+            swaps.push((cycle[i - 1], cycle[length - 2 - i]));
+        }
+    }
+
+    swaps
+}
+
+#[pyfunction]
+#[pyo3(signature = (pattern))]
+fn _synth_permutation_acg(py: Python, pattern: PyArrayLike1<i64>) -> PyResult<CircuitData> {
+    let view = pattern.as_array();
+    let num_qubits = view.len();
+    let cycles = pattern_to_cycles(&view, &true);
+    let swaps = decompose_cycles(&cycles);
+
+    CircuitData::from_standard_gates(
+        py,
+        num_qubits as u32,
+        swaps.iter().map(|(i, j)| {
+            (
+                StandardGate::SwapGate,
+                smallvec![],
+                smallvec![Qubit(*i as u32), Qubit(*j as u32)],
+            )
+        }),
+        Param::Float(0.0),
+    )
+}
+
 /// Checks whether an array of size N is a permutation of 0, 1, ..., N - 1.
 #[pyfunction]
 #[pyo3(signature = (pattern))]
@@ -111,10 +196,31 @@ fn _get_ordered_swap(py: Python, permutation_in: PyArrayLike1<i64>) -> PyResult<
     Ok(get_ordered_swap(&view).to_object(py))
 }
 
+/// Decompose a SWAP pattern into a series of SWAP gate indices to implement them.
+/// For example, let the pattern be [1, 2, 3, 4, 0, 6, 5], which contains the two cycles
+/// [1, 2, 3, 4, 0] and [6, 5]. These can then be implemented with the SWAPs
+/// [(0, 3), (1, 2), (0, 4), (1, 3)], respectively [(6, 5)].
+/// If ``invert`` is True, reverse the indices before computing the SWAPs.
+#[pyfunction]
+#[pyo3(signature = (pattern, invert_order=true))]
+fn _pattern_to_swaps(
+    py: Python,
+    pattern: PyArrayLike1<i64>,
+    invert_order: bool,
+) -> PyResult<PyObject> {
+    let view = pattern.as_array();
+    let cycles = pattern_to_cycles(&view, &invert_order);
+    let swaps = decompose_cycles(&cycles);
+    let swaps_i64: Vec<(i64, i64)> = swaps.iter().map(|&x| (x.0 as i64, x.1 as i64)).collect();
+    Ok(swaps_i64.to_object(py))
+}
+
 #[pymodule]
 pub fn permutation(m: &Bound<PyModule>) -> PyResult<()> {
     m.add_function(wrap_pyfunction!(_validate_permutation, m)?)?;
     m.add_function(wrap_pyfunction!(_inverse_pattern, m)?)?;
     m.add_function(wrap_pyfunction!(_get_ordered_swap, m)?)?;
+    m.add_function(wrap_pyfunction!(_pattern_to_swaps, m)?)?;
+    m.add_function(wrap_pyfunction!(_synth_permutation_acg, m)?)?;
     Ok(())
 }

qiskit/synthesis/permutation/permutation_full.py

@@ -16,12 +16,8 @@
 
 import numpy as np
 from qiskit.circuit.quantumcircuit import QuantumCircuit
-from .permutation_utils import (
-    _get_ordered_swap,
-    _inverse_pattern,
-    _pattern_to_cycles,
-    _decompose_cycles,
-)
+from qiskit._accelerate.permutation import _synth_permutation_acg
+from .permutation_utils import _get_ordered_swap
 
 
 def synth_permutation_basic(pattern: list[int] | np.ndarray[int]) -> QuantumCircuit:
@@ -87,16 +83,4 @@ def synth_permutation_acg(pattern: list[int] | np.ndarray[int]) -> QuantumCircui
            *Routing Permutations on Graphs Via Matchings.*,
            `(Full paper) <https://www.cs.tau.ac.il/~nogaa/PDFS/r.pdf>`_
     """
-
-    num_qubits = len(pattern)
-    qc = QuantumCircuit(num_qubits)
-
-    # invert pattern (Qiskit notation is opposite)
-    cur_pattern = _inverse_pattern(pattern)
-    cycles = _pattern_to_cycles(cur_pattern)
-    swaps = _decompose_cycles(cycles)
-
-    for swap in swaps:
-        qc.swap(swap[0], swap[1])
-
-    return qc
+    return QuantumCircuit._from_circuit_data(_synth_permutation_acg(pattern))

qiskit/synthesis/permutation/permutation_utils.py

@@ -17,33 +17,5 @@
     _inverse_pattern,
     _get_ordered_swap,
     _validate_permutation,
+    _pattern_to_swaps,
 )
-
-
-def _pattern_to_cycles(pattern):
-    """Given a permutation pattern, creates its disjoint cycle decomposition."""
-    nq = len(pattern)
-    explored = [False] * nq
-    cycles = []
-    for i in pattern:
-        cycle = []
-        while not explored[i]:
-            cycle.append(i)
-            explored[i] = True
-            i = pattern[i]
-        if len(cycle) >= 2:
-            cycles.append(cycle)
-    return cycles
-
-
-def _decompose_cycles(cycles):
-    """Given a disjoint cycle decomposition, decomposes every cycle into a SWAP
-    circuit of depth 2."""
-    swap_list = []
-    for cycle in cycles:
-        m = len(cycle)
-        for i in range((m - 1) // 2):
-            swap_list.append((cycle[i - 1], cycle[m - 3 - i]))
-        for i in range(m // 2):
-            swap_list.append((cycle[i - 1], cycle[m - 2 - i]))
-    return swap_list

releasenotes/notes/oxidize-acg-0294a87c0d5974fa.yaml

@@ -0,0 +1,5 @@
+---
+upgrade_synthesis:
+  - |
+    Port :func:`.synth_permutation_acg`, used to synthesize qubit permutations, to Rust.
+    This produces an approximate 3x performance improvement on 1000 qubit circuits.