Port synth_cnot_phase_aam to Rust

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

@@ -14,7 +14,7 @@ use crate::QiskitError;
 use numpy::{IntoPyArray, PyArray2, PyReadonlyArray2, PyReadwriteArray2};
 use pyo3::prelude::*;
 
-mod pmh;
+pub mod pmh;
 pub mod utils;
 
 #[pyfunction]

crates/accelerate/src/synthesis/linear/pmh.rs

@@ -13,7 +13,7 @@
 use hashbrown::HashMap;
 use ndarray::{s, Array1, Array2, ArrayViewMut2, Axis};
 use numpy::PyReadonlyArray2;
-use smallvec::smallvec;
+use smallvec::{smallvec, SmallVec};
 use std::cmp;
 
 use qiskit_circuit::circuit_data::CircuitData;
@@ -159,7 +159,16 @@ pub fn synth_cnot_count_full_pmh(
     section_size: Option<i64>,
 ) -> PyResult<CircuitData> {
     let arrayview = matrix.as_array();
-    let mut mat: Array2<bool> = arrayview.to_owned();
+    let mat: Array2<bool> = arrayview.to_owned();
+    let num_qubits = mat.nrows();
+
+    let instructions = _synth_cnot_count_full_pmh(mat, section_size);
+    CircuitData::from_standard_gates(py, num_qubits as u32, instructions, Param::Float(0.0))
+}
+
+type Instructions = (StandardGate, SmallVec<[Param; 3]>, SmallVec<[Qubit; 2]>);
+pub fn _synth_cnot_count_full_pmh(mat: Array2<bool>, sec_size: Option<i64>) -> Vec<Instructions> {
+    let mut mat = mat;
     let num_qubits = mat.nrows(); // is a quadratic matrix
 
     // If given, use the user-specified input size. If None, we default to
@@ -168,7 +177,7 @@ pub fn synth_cnot_count_full_pmh(
     // In addition, we at least set a block size of 2, which, in practice, minimizes the bound
     // until ~100 qubits.
     let alpha = 0.56;
-    let blocksize = match section_size {
+    let blocksize = match sec_size {
         Some(section_size) => section_size as usize,
         None => std::cmp::max(2, (alpha * (num_qubits as f64).log2()).floor() as usize),
     };
@@ -190,6 +199,5 @@ pub fn synth_cnot_count_full_pmh(
                 smallvec![Qubit(ctrl as u32), Qubit(target as u32)],
             )
         });
-
-    CircuitData::from_standard_gates(py, num_qubits as u32, instructions, Param::Float(0.0))
+    instructions.collect()
 }

crates/accelerate/src/synthesis/linear_phase/cnot_phase_synth.rs

@@ -0,0 +1,292 @@
+// 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 crate::synthesis::linear::pmh::_synth_cnot_count_full_pmh;
+use ndarray::Array2;
+use numpy::PyReadonlyArray2;
+use pyo3::{prelude::*, types::PyList};
+use qiskit_circuit::circuit_data::CircuitData;
+use qiskit_circuit::operations::{Param, StandardGate};
+use qiskit_circuit::Qubit;
+use smallvec::smallvec;
+use std::f64::consts::PI;
+
+/// This function implements a Gray-code inspired algorithm of synthesizing a circuit
+/// over CNOT and phase-gates with minimal-CNOT for a given phase-polynomial.
+/// The algorithm is described as "Gray-Synth" algorithm in Algorithm-1, page 12
+/// of paper "https://arxiv.org/abs/1712.01859".
+#[pyfunction]
+#[pyo3(signature = (cnots, angles, section_size=2))]
+pub fn synth_cnot_phase_aam(
+    py: Python,
+    cnots: PyReadonlyArray2<u8>,
+    angles: &Bound<PyList>,
+    section_size: Option<i64>,
+) -> PyResult<CircuitData> {
+    let s = cnots.as_array().to_owned();
+    let num_qubits = s.nrows();
+    let mut s_cpy = s.clone();
+    let mut instructions = vec![];
+
+    let mut rust_angles: Vec<String> = angles
+        .iter()
+        .filter_map(|data| data.extract::<String>().ok())
+        .collect();
+    let mut state = Array2::<u8>::eye(num_qubits);
+
+    for qubit_idx in 0..num_qubits {
+        let mut index = 0_usize;
+        let mut swtch: bool = true;
+
+        while index < s_cpy.ncols() {
+            let icnot = s_cpy.column(index).to_vec();
+            if icnot == state.row(qubit_idx).to_vec() {
+                match rust_angles.remove(index) {
+                    gate if gate == "t" => instructions.push((
+                        StandardGate::TGate,
+                        smallvec![],
+                        smallvec![Qubit(qubit_idx as u32)],
+                    )),
+                    gate if gate == "tdg" => instructions.push((
+                        StandardGate::TdgGate,
+                        smallvec![],
+                        smallvec![Qubit(qubit_idx as u32)],
+                    )),
+                    gate if gate == "s" => instructions.push((
+                        StandardGate::SGate,
+                        smallvec![],
+                        smallvec![Qubit(qubit_idx as u32)],
+                    )),
+                    gate if gate == "sdg" => instructions.push((
+                        StandardGate::SdgGate,
+                        smallvec![],
+                        smallvec![Qubit(qubit_idx as u32)],
+                    )),
+                    gate if gate == "z" => instructions.push((
+                        StandardGate::ZGate,
+                        smallvec![],
+                        smallvec![Qubit(qubit_idx as u32)],
+                    )),
+                    angles_in_pi => instructions.push((
+                        StandardGate::PhaseGate,
+                        smallvec![Param::Float((angles_in_pi.parse::<f64>()?) % PI)],
+                        smallvec![Qubit(qubit_idx as u32)],
+                    )),
+                };
+                s_cpy.remove_index(numpy::ndarray::Axis(1), index);
+                if index == s_cpy.ncols() {
+                    break;
+                }
+                if index == 0 {
+                    swtch = false;
+                } else {
+                    index -= 1;
+                }
+            }
+            if swtch {
+                index += 1;
+            } else {
+                swtch = true;
+            }
+        }
+    }
+
+    let epsilion: usize = num_qubits;
+    let mut q = vec![(s, (0..num_qubits).collect::<Vec<usize>>(), epsilion)];
+
+    while !q.is_empty() {
+        let (mut _s, mut _i, mut _ep) = q.pop().unwrap();
+
+        if _s.is_empty() {
+            continue;
+        }
+
+        if _ep < num_qubits {
+            let mut condition = true;
+            while condition {
+                condition = false;
+
+                for _j in 0..num_qubits {
+                    if (_j != _ep) && (_s.row(_j).sum() as usize == _s.row(_j).len()) {
+                        condition = true;
+                        instructions.push((
+                            StandardGate::CXGate,
+                            smallvec![],
+                            smallvec![Qubit(_j as u32), Qubit(_ep as u32)],
+                        ));
+
+                        for _k in 0..state.ncols() {
+                            state[(_ep, _k)] ^= state[(_j, _k)];
+                        }
+
+                        let mut index = 0_usize;
+                        let mut swtch: bool = true;
+                        while index < s_cpy.ncols() {
+                            let icnot = s_cpy.column(index).to_vec();
+                            if icnot == state.row(_ep).to_vec() {
+                                match rust_angles.remove(index) {
+                                    gate if gate == "t" => instructions.push((
+                                        StandardGate::TGate,
+                                        smallvec![],
+                                        smallvec![Qubit(_ep as u32)],
+                                    )),
+                                    gate if gate == "tdg" => instructions.push((
+                                        StandardGate::TdgGate,
+                                        smallvec![],
+                                        smallvec![Qubit(_ep as u32)],
+                                    )),
+                                    gate if gate == "s" => instructions.push((
+                                        StandardGate::SGate,
+                                        smallvec![],
+                                        smallvec![Qubit(_ep as u32)],
+                                    )),
+                                    gate if gate == "sdg" => instructions.push((
+                                        StandardGate::SdgGate,
+                                        smallvec![],
+                                        smallvec![Qubit(_ep as u32)],
+                                    )),
+                                    gate if gate == "z" => instructions.push((
+                                        StandardGate::ZGate,
+                                        smallvec![],
+                                        smallvec![Qubit(_ep as u32)],
+                                    )),
+                                    angles_in_pi => instructions.push((
+                                        StandardGate::PhaseGate,
+                                        smallvec![Param::Float(
+                                            (angles_in_pi.parse::<f64>()?) % PI
+                                        )],
+                                        smallvec![Qubit(_ep as u32)],
+                                    )),
+                                };
+                                s_cpy.remove_index(numpy::ndarray::Axis(1), index);
+                                if index == s_cpy.ncols() {
+                                    break;
+                                }
+                                if index == 0 {
+                                    swtch = false;
+                                } else {
+                                    index -= 1;
+                                }
+                            }
+                            if swtch {
+                                index += 1;
+                            } else {
+                                swtch = true;
+                            }
+                        }
+
+                        q.push((_s, _i, _ep));
+                        let mut unique_q = vec![];
+                        for data in q.into_iter() {
+                            if !unique_q.contains(&data) {
+                                unique_q.push(data);
+                            }
+                        }
+
+                        q = unique_q;
+
+                        for data in &mut q {
+                            let (ref mut _temp_s, _, _) = data;
+
+                            if _temp_s.is_empty() {
+                                continue;
+                            }
+
+                            for idx in 0.._temp_s.row(_j).len() {
+                                _temp_s[(_j, idx)] ^= _temp_s[(_ep, idx)];
+                            }
+                        }
+
+                        (_s, _i, _ep) = q.pop().unwrap();
+                    }
+                }
+            }
+        }
+
+        if _i.is_empty() {
+            continue;
+        }
+
+        let maxes: Vec<usize> = _s
+            .axis_iter(numpy::ndarray::Axis(0))
+            .map(|row| {
+                std::cmp::max(
+                    row.iter().filter(|&&x| x == 0).count(),
+                    row.iter().filter(|&&x| x == 1).count(),
+                )
+            })
+            .collect();
+
+        let maxes2: Vec<usize> = _i.iter().map(|&_i_idx| maxes[_i_idx]).collect();
+
+        let _temp_argmax = maxes2
+            .iter()
+            .enumerate()
+            .max_by(|(_, x), (_, y)| x.cmp(y))
+            .map(|(idx, _)| idx)
+            .unwrap();
+
+        let _j = _i[_temp_argmax];
+
+        let mut cnots0_t = vec![];
+        let mut cnots1_t = vec![];
+
+        let mut cnots0_t_shape = (0_usize, _s.column(0).len());
+        let mut cnots1_t_shape = (0_usize, 0_usize);
+        cnots1_t_shape.1 = cnots0_t_shape.1;
+        for cols in _s.columns() {
+            if cols[_j] == 0 {
+                cnots0_t_shape.0 += 1;
+                cnots0_t.append(&mut cols.to_vec());
+            } else if cols[_j] == 1 {
+                cnots1_t_shape.0 += 1;
+                cnots1_t.append(&mut cols.to_vec());
+            }
+        }
+
+        let cnots0 =
+            Array2::from_shape_vec((cnots0_t_shape.0, cnots0_t_shape.1), cnots0_t).unwrap();
+        let cnots1 =
+            Array2::from_shape_vec((cnots1_t_shape.0, cnots1_t_shape.1), cnots1_t).unwrap();
+
+        let cnots0 = cnots0.reversed_axes().to_owned();
+        let cnots1 = cnots1.reversed_axes().to_owned();
+
+        if _ep == num_qubits {
+            q.push((
+                cnots1,
+                _i.clone().into_iter().filter(|&x| x != _j).collect(),
+                _j,
+            ));
+        } else {
+            q.push((
+                cnots1,
+                _i.clone().into_iter().filter(|&x| x != _j).collect(),
+                _ep,
+            ));
+        }
+
+        q.push((
+            cnots0,
+            _i.clone().into_iter().filter(|&x| x != _j).collect(),
+            _ep,
+        ));
+    }
+
+    let state_bool = state.mapv(|x| x != 0);
+    let mut instrs = _synth_cnot_count_full_pmh(state_bool, section_size)
+        .into_iter()
+        .rev()
+        .collect();
+    instructions.append(&mut instrs);
+    CircuitData::from_standard_gates(py, num_qubits as u32, instructions, Param::Float(0.0))
+}

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

@@ -10,6 +10,7 @@
 // copyright notice, and modified files need to carry a notice indicating
 // that they have been altered from the originals.
 
+mod cnot_phase_synth;
 use numpy::PyReadonlyArray2;
 use pyo3::{
     prelude::*,
@@ -42,5 +43,6 @@ fn synth_cz_depth_line_mr(py: Python, mat: PyReadonlyArray2<bool>) -> PyResult<C
 
 pub fn linear_phase(m: &Bound<PyModule>) -> PyResult<()> {
     m.add_wrapped(wrap_pyfunction!(synth_cz_depth_line_mr))?;
+    m.add_wrapped(wrap_pyfunction!(cnot_phase_synth::synth_cnot_phase_aam))?;
     Ok(())
 }