Add Rustiq-based synthesis for PauliEvolutionGate (#13301)

* py version for expand

* starting to write some code

* implementing

* cleanup

* cleanup

* expand fully & simplify lie trotter

* use examples that actually do not commute

* add plugin structure

* fixing global phase for all-I rotations

* fixes

* fixing plugin names

* minor

* removing a random print statement

* additional improvements

* improving rustiq plugin

* merge with #13239

* Adding pauli evolution plugins to docstrings

* adding documentation on rustiq plugin

* fixes after refactoring

* typo

* more merges with #13295; adding more Rustiq tests

* more efficient append_sx and append_sxdg gates for cliffords

* review comments

* moving the pauli network synthesis logic into a separate file

* some code review suggestions

* simplifying the code by merging the oredered and unorderd version of rotation injection

* more review comments

* adding python tests

* more code review suggestions

* more review comments

* more review comments

* test for preserve_order

* lint

* upgrading rustiq-core to 0.0.10

* clippy: removing mutable ref

* Improving PauliEvolution synthesis tests.

Making sure that the number of rotation gates in the synthesized
circuit equals the number of non-trivial Pauli rotation gates.

* documentation fixes after the merge

---------

Co-authored-by: Julien Gacon <[email protected]>

crates/accelerate/Cargo.toml

@@ -26,6 +26,7 @@ qiskit-circuit.workspace = true
 thiserror.workspace = true
 ndarray_einsum_beta = "0.7"
 once_cell = "1.20.2"
+rustiq-core = "0.0.10"
 bytemuck.workspace = true
 
 [dependencies.smallvec]

crates/accelerate/src/synthesis/clifford/greedy_synthesis.rs

@@ -15,8 +15,8 @@ use indexmap::IndexSet;
 use ndarray::{s, ArrayView2};
 use smallvec::smallvec;
 
-use crate::synthesis::clifford::utils::CliffordGatesVec;
 use crate::synthesis::clifford::utils::{adjust_final_pauli_gates, SymplecticMatrix};
+use crate::synthesis::clifford::utils::{Clifford, CliffordGatesVec};
 use qiskit_circuit::operations::StandardGate;
 use qiskit_circuit::Qubit;
 
@@ -437,3 +437,14 @@ impl GreedyCliffordSynthesis<'_> {
         Ok((self.num_qubits, clifford_gates))
     }
 }
+
+/// Resynthesizes a clifford circuit using the greedy Clifford synthesis algorithm.
+pub fn resynthesize_clifford_circuit(
+    num_qubits: usize,
+    gates: &CliffordGatesVec,
+) -> Result<CliffordGatesVec, String> {
+    let sim_clifford = Clifford::from_gate_sequence(gates, num_qubits)?;
+    let mut synthesis = GreedyCliffordSynthesis::new(sim_clifford.tableau.view())?;
+    let (_, new_gates) = synthesis.run()?;
+    Ok(new_gates)
+}

crates/accelerate/src/synthesis/clifford/mod.rs

@@ -11,9 +11,9 @@
 // that they have been altered from the originals.
 
 mod bm_synthesis;
-mod greedy_synthesis;
+pub(crate) mod greedy_synthesis;
 mod random_clifford;
-mod utils;
+pub(crate) mod utils;
 
 use crate::synthesis::clifford::bm_synthesis::synth_clifford_bm_inner;
 use crate::synthesis::clifford::greedy_synthesis::GreedyCliffordSynthesis;

crates/accelerate/src/synthesis/clifford/utils.rs

@@ -149,6 +149,30 @@ impl Clifford {
         azip!((z in &mut z, &x in &x) *z ^= x);
     }
 
+    /// Modifies the tableau in-place by appending SX-gate
+    pub fn append_sx(&mut self, qubit: usize) {
+        let (mut x, z, mut p) = self.tableau.multi_slice_mut((
+            s![.., qubit],
+            s![.., self.num_qubits + qubit],
+            s![.., 2 * self.num_qubits],
+        ));
+
+        azip!((p in &mut p, &x in &x, &z in &z)  *p ^= !x & z);
+        azip!((&z in &z, x in &mut x) *x ^= z);
+    }
+
+    /// Modifies the tableau in-place by appending SXDG-gate
+    pub fn append_sxdg(&mut self, qubit: usize) {
+        let (mut x, z, mut p) = self.tableau.multi_slice_mut((
+            s![.., qubit],
+            s![.., self.num_qubits + qubit],
+            s![.., 2 * self.num_qubits],
+        ));
+
+        azip!((p in &mut p, &x in &x, &z in &z)  *p ^= x & z);
+        azip!((&z in &z, x in &mut x) *x ^= z);
+    }
+
     /// Modifies the tableau in-place by appending H-gate
     pub fn append_h(&mut self, qubit: usize) {
         let (mut x, mut z, mut p) = self.tableau.multi_slice_mut((
@@ -227,6 +251,18 @@ impl Clifford {
                     clifford.append_s(qubits[0].index());
                     Ok(())
                 }
+                StandardGate::SdgGate => {
+                    clifford.append_sdg(qubits[0].0 as usize);
+                    Ok(())
+                }
+                StandardGate::SXGate => {
+                    clifford.append_sx(qubits[0].0 as usize);
+                    Ok(())
+                }
+                StandardGate::SXdgGate => {
+                    clifford.append_sxdg(qubits[0].0 as usize);
+                    Ok(())
+                }
                 StandardGate::HGate => {
                     clifford.append_h(qubits[0].index());
                     Ok(())
@@ -283,21 +319,18 @@ pub fn adjust_final_pauli_gates(
     // add pauli gates
     for qubit in 0..num_qubits {
         if delta_phase_pre[qubit] && delta_phase_pre[qubit + num_qubits] {
-            // println!("=> Adding Y-gate on {}", qubit);
             gate_seq.push((
                 StandardGate::YGate,
                 smallvec![],
                 smallvec![Qubit::new(qubit)],
             ));
         } else if delta_phase_pre[qubit] {
-            // println!("=> Adding Z-gate on {}", qubit);
             gate_seq.push((
                 StandardGate::ZGate,
                 smallvec![],
                 smallvec![Qubit::new(qubit)],
             ));
         } else if delta_phase_pre[qubit + num_qubits] {
-            // println!("=> Adding X-gate on {}", qubit);
             gate_seq.push((
                 StandardGate::XGate,
                 smallvec![],

crates/accelerate/src/synthesis/evolution/mod.rs

@@ -0,0 +1,79 @@
+// 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.
+
+mod pauli_network;
+
+use pyo3::prelude::*;
+use pyo3::types::PyList;
+
+use qiskit_circuit::circuit_data::CircuitData;
+
+use crate::synthesis::evolution::pauli_network::pauli_network_synthesis_inner;
+
+/// Calls Rustiq's pauli network synthesis algorithm and returns the
+/// Qiskit circuit data with Clifford gates and rotations.
+///
+/// # Arguments
+///
+/// * py: a GIL handle, needed to add and negate rotation parameters in Python space.
+/// * num_qubits: total number of qubits.
+/// * pauli_network: pauli network represented in sparse format. It's a list
+///     of triples such as `[("XX", [0, 3], theta), ("ZZ", [0, 1], 0.1)]`.
+/// * optimize_count: if `true`, Rustiq's synthesis algorithms aims to optimize
+///     the 2-qubit gate count; and if `false`, then the 2-qubit depth.
+/// * preserve_order: whether the order of paulis should be preserved, up to
+///     commutativity. If the order is not preserved, the returned circuit will
+///     generally not be equivalent to the given pauli network.
+/// * upto_clifford: if `true`, the final Clifford operator is not synthesized
+///     and the returned circuit will generally not be equivalent to the given
+///     pauli network. In addition, the argument `upto_phase` would be ignored.
+/// * upto_phase: if `true`, the global phase of the returned circuit may differ
+///     from the global phase of the given pauli network. The argument is considered
+///     to be `true` when `upto_clifford` is `true`.
+/// * resynth_clifford_method: describes the strategy to synthesize the final
+///     Clifford operator. If `0` a naive approach is used, which doubles the number
+///     of gates but preserves the global phase of the circuit. If `1`, the Clifford is
+///     resynthesized using Qiskit's greedy Clifford synthesis algorithm. If `2`, it
+///     is resynthesized by Rustiq itself. If `upto_phase` is `false`, the naive
+///     approach is used, as neither synthesis method preserves the global phase.
+///
+/// If `preserve_order` is `true` and both `upto_clifford` and `upto_phase` are `false`,
+/// the returned circuit is equivalent to the given pauli network.
+#[pyfunction]
+#[pyo3(signature = (num_qubits, pauli_network, optimize_count=true, preserve_order=true, upto_clifford=false, upto_phase=false, resynth_clifford_method=1))]
+#[allow(clippy::too_many_arguments)]
+pub fn pauli_network_synthesis(
+    py: Python,
+    num_qubits: usize,
+    pauli_network: &Bound<PyList>,
+    optimize_count: bool,
+    preserve_order: bool,
+    upto_clifford: bool,
+    upto_phase: bool,
+    resynth_clifford_method: usize,
+) -> PyResult<CircuitData> {
+    pauli_network_synthesis_inner(
+        py,
+        num_qubits,
+        pauli_network,
+        optimize_count,
+        preserve_order,
+        upto_clifford,
+        upto_phase,
+        resynth_clifford_method,
+    )
+}
+
+pub fn evolution(m: &Bound<PyModule>) -> PyResult<()> {
+    m.add_function(wrap_pyfunction!(pauli_network_synthesis, m)?)?;
+    Ok(())
+}