Add Rust matrix support for PauliProductRotation

crates/circuit/Cargo.toml

@@ -64,4 +64,4 @@ features = ["nalgebra"]
 cache_pygates = []
 
 [dev-dependencies]
-pyo3 = { workspace = true, features = ["auto-initialize"] }
+pyo3 = { workspace = true, features = ["auto-initialize"] }

crates/circuit/src/circuit_drawer.rs

@@ -1216,6 +1216,8 @@ mod tests {
     };
     use crate::parameter::parameter_expression::ParameterExpression;
     use crate::parameter::symbol_expr::Symbol;
+    #[cfg(feature = "cache_pygates")]
+    use std::sync::OnceLock;
 
     fn basic_circuit() -> CircuitData {
         let qreg = QuantumRegister::new_owning("q", 2);
@@ -1407,6 +1409,8 @@ q_1: ┤ H ├┤1            ├┤1               ├┤ my_ch ├
             clbits: circuit.cargs_interner().get_default(),
             params: None,
             label: None,
+            #[cfg(feature = "cache_pygates")]
+            py_op: OnceLock::new(),
         };
         circuit.push(inst).unwrap();
         let inst = PackedInstruction {
@@ -1418,6 +1422,8 @@ q_1: ┤ H ├┤1            ├┤1               ├┤ my_ch ├
             clbits: circuit.cargs_interner().get_default(),
             params: None,
             label: Some(Box::new("my little identity".to_string())),
+            #[cfg(feature = "cache_pygates")]
+            py_op: OnceLock::new(),
         };
         circuit.push(inst).unwrap();
         let inst = PackedInstruction {
@@ -1429,6 +1435,8 @@ q_1: ┤ H ├┤1            ├┤1               ├┤ my_ch ├
             clbits: circuit.cargs_interner().get_default(),
             params: None,
             label: None,
+            #[cfg(feature = "cache_pygates")]
+            py_op: OnceLock::new(),
         };
         circuit.push(inst).unwrap();
         let inst = PackedInstruction {
@@ -1440,6 +1448,8 @@ q_1: ┤ H ├┤1            ├┤1               ├┤ my_ch ├
             clbits: circuit.cargs_interner().get_default(),
             params: None,
             label: Some(Box::new("my small identity".to_string())),
+            #[cfg(feature = "cache_pygates")]
+            py_op: OnceLock::new(),
         };
         circuit.push(inst).unwrap();
         let inst = PackedInstruction {
@@ -1451,6 +1461,8 @@ q_1: ┤ H ├┤1            ├┤1               ├┤ my_ch ├
             clbits: circuit.cargs_interner().get_default(),
             params: None,
             label: None,
+            #[cfg(feature = "cache_pygates")]
+            py_op: OnceLock::new(),
         };
         circuit.push(inst).unwrap();
         let inst = PackedInstruction {
@@ -1462,6 +1474,8 @@ q_1: ┤ H ├┤1            ├┤1               ├┤ my_ch ├
             clbits: circuit.cargs_interner().get_default(),
             params: None,
             label: Some(Box::new("my medium identity".to_string())),
+            #[cfg(feature = "cache_pygates")]
+            py_op: OnceLock::new(),
         };
         circuit.push(inst).unwrap();
         let inst = PackedInstruction {
@@ -1473,6 +1487,8 @@ q_1: ┤ H ├┤1            ├┤1               ├┤ my_ch ├
             clbits: circuit.cargs_interner().get_default(),
             params: None,
             label: None,
+            #[cfg(feature = "cache_pygates")]
+            py_op: OnceLock::new(),
         };
         circuit.push(inst).unwrap();
         let inst = PackedInstruction {
@@ -1484,6 +1500,8 @@ q_1: ┤ H ├┤1            ├┤1               ├┤ my_ch ├
             clbits: circuit.cargs_interner().get_default(),
             params: None,
             label: Some(Box::new("my bigger identity".to_string())),
+            #[cfg(feature = "cache_pygates")]
+            py_op: OnceLock::new(),
         };
         circuit.push(inst).unwrap();
         let result = draw_circuit(&circuit, false, false, Some(80)).unwrap();
@@ -1535,7 +1553,7 @@ q_3: ──────────────────────┤1
                 .map(|x| Qubit(x + (i as u32 % (num_wires - num_qubits))))
                 .collect::<Vec<_>>();
             let params = (0..num_params)
-                .map(|_x| 3.141.into())
+                .map(|_x| std::f64::consts::PI.into())
                 .collect::<Vec<Param>>();
             circuit.push_standard_gate(op, &params, &qubits).unwrap();
         }
@@ -1624,7 +1642,9 @@ q_4: ─────────────────────────
     #[test]
     fn test_global_phase() {
         let mut circuit = basic_circuit();
-        circuit.set_global_phase_param(3.14.into()).unwrap();
+        circuit
+            .set_global_phase_param(std::f64::consts::PI.into())
+            .unwrap();
         let result = draw_circuit(&circuit, true, false, None).unwrap();
 
         let expected = "
@@ -1729,6 +1749,8 @@ q_1: ┤1         ├┤1            ├┤1         ├
                 clbits: circuit.cargs_interner().get_default(),
                 params: None,
                 label: None,
+                #[cfg(feature = "cache_pygates")]
+                py_op: OnceLock::new(),
             };
             circuit.push(inst).unwrap();
         }
@@ -1747,18 +1769,22 @@ q_1: ┤1         ├┤1            ├┤1         ├
                     clbits: circuit.cargs_interner().get_default(),
                     params: Some(Box::new(Parameters::Params(smallvec![param]))),
                     label: None,
+                    #[cfg(feature = "cache_pygates")]
+                    py_op: OnceLock::new(),
                 };
                 circuit.push(inst).unwrap();
             }
         }
         for i in [1, 2, 3, 4] {
-            let qubits = (0..i).map(|x| Qubit::new(x)).collect::<Vec<_>>();
+            let qubits = (0..i).map(Qubit::new).collect::<Vec<_>>();
             let inst = PackedInstruction {
                 op: StandardInstruction::Barrier(i as u32).into(),
                 qubits: circuit.add_qargs(&qubits),
                 clbits: circuit.cargs_interner().get_default(),
                 params: None,
                 label: None,
+                #[cfg(feature = "cache_pygates")]
+                py_op: OnceLock::new(),
             };
             circuit.push(inst).unwrap();
         }
@@ -1769,6 +1795,8 @@ q_1: ┤1         ├┤1            ├┤1         ├
                 clbits: circuit.add_cargs(&[Clbit::new(i)]),
                 params: None,
                 label: None,
+                #[cfg(feature = "cache_pygates")]
+                py_op: OnceLock::new(),
             };
             circuit.push(inst).unwrap();
         }

crates/circuit/src/gate_matrix.rs

@@ -10,6 +10,7 @@
 // copyright notice, and modified files need to carry a notice indicating
 // that they have been altered from the originals.
 
+use ndarray::{Array2, array};
 use num_complex::Complex64;
 use std::f64::consts::FRAC_1_SQRT_2;
 
@@ -516,3 +517,56 @@ pub fn xx_plus_yy_gate(theta: f64, beta: f64) -> GateArray2Q {
         [C_ZERO, C_ZERO, C_ZERO, C_ONE],
     ]
 }
+
+pub fn pauli_zx_to_matrix(z: &[bool], x: &[bool]) -> Array2<Complex64> {
+    debug_assert_eq!(z.len(), x.len());
+
+    let i = array![
+        [Complex64::new(1.0, 0.0), Complex64::new(0.0, 0.0)],
+        [Complex64::new(0.0, 0.0), Complex64::new(1.0, 0.0)],
+    ];
+    let x_mat = array![
+        [Complex64::new(0.0, 0.0), Complex64::new(1.0, 0.0)],
+        [Complex64::new(1.0, 0.0), Complex64::new(0.0, 0.0)],
+    ];
+    let z_mat = array![
+        [Complex64::new(1.0, 0.0), Complex64::new(0.0, 0.0)],
+        [Complex64::new(0.0, 0.0), Complex64::new(-1.0, 0.0)],
+    ];
+    let y_mat = array![
+        [Complex64::new(0.0, 0.0), Complex64::new(0.0, -1.0)],
+        [Complex64::new(0.0, 1.0), Complex64::new(0.0, 0.0)],
+    ];
+
+    fn kron(a: &Array2<Complex64>, b: &Array2<Complex64>) -> Array2<Complex64> {
+        let (ar, ac) = a.dim();
+        let (br, bc) = b.dim();
+        let mut out = Array2::<Complex64>::zeros((ar * br, ac * bc));
+
+        for i in 0..ar {
+            for j in 0..ac {
+                let coeff = a[(i, j)];
+                for k in 0..br {
+                    for l in 0..bc {
+                        out[(i * br + k, j * bc + l)] = coeff * b[(k, l)];
+                    }
+                }
+            }
+        }
+        out
+    }
+
+    let mut result = array![[Complex64::new(1.0, 0.0)]];
+
+    for (&z_bit, &x_bit) in z.iter().zip(x.iter()) {
+        let factor = match (z_bit, x_bit) {
+            (false, false) => &i,
+            (false, true) => &x_mat,
+            (true, false) => &z_mat,
+            (true, true) => &y_mat,
+        };
+        result = kron(&result, factor);
+    }
+
+    result
+}

crates/circuit/src/operations.rs

@@ -9,7 +9,6 @@
 // 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 approx::relative_eq;
 use std::f64::consts::PI;
 use std::fmt::Debug;
@@ -3262,24 +3261,6 @@ impl PyInstruction {
         })
     }
 
-    pub fn matrix(&self) -> Option<Array2<Complex64>> {
-        Python::attach(|py| -> Option<Array2<Complex64>> {
-            match self.instruction.getattr(py, intern!(py, "to_matrix")) {
-                Ok(to_matrix) => {
-                    let res: Option<Py<PyAny>> = to_matrix.call0(py).ok()?.extract(py).ok();
-                    match res {
-                        Some(x) => {
-                            let array: PyReadonlyArray2<Complex64> = x.extract(py).ok()?;
-                            Some(array.as_array().to_owned())
-                        }
-                        None => None,
-                    }
-                }
-                Err(_) => None,
-            }
-        })
-    }
-
     pub fn definition(&self) -> Option<CircuitData> {
         Python::attach(|py| -> Option<CircuitData> {
             match self.instruction.getattr(py, intern!(py, "definition")) {
@@ -3311,6 +3292,24 @@ impl PyInstruction {
         })
     }
 
+    pub fn matrix(&self) -> Option<Array2<Complex64>> {
+        Python::attach(|py| -> Option<Array2<Complex64>> {
+            match self.instruction.getattr(py, intern!(py, "to_matrix")) {
+                Ok(to_matrix) => {
+                    let res: Option<Py<PyAny>> = to_matrix.call0(py).ok()?.extract(py).ok();
+                    match res {
+                        Some(x) => {
+                            let array: PyReadonlyArray2<Complex64> = x.extract(py).ok()?;
+                            Some(array.as_array().to_owned())
+                        }
+                        None => None,
+                    }
+                }
+                Err(_) => None,
+            }
+        })
+    }
+
     pub fn matrix_as_static_2q(&self) -> Option<[[Complex64; 4]; 4]> {
         if self.num_qubits() != 2 {
             return None;
@@ -3550,24 +3549,6 @@ pub struct PauliProductRotation {
     pub angle: Param,
 }
 
-impl Operation for PauliProductRotation {
-    fn name(&self) -> &str {
-        "pauli_product_rotation"
-    }
-    fn num_qubits(&self) -> u32 {
-        self.z.len() as u32
-    }
-    fn num_clbits(&self) -> u32 {
-        0
-    }
-    fn num_params(&self) -> u32 {
-        1
-    }
-    fn directive(&self) -> bool {
-        false
-    }
-}
-
 impl PauliProductRotation {
     pub fn create_py_op(&self, py: Python, label: Option<&str>) -> PyResult<Py<PyAny>> {
         let z = self.z.to_pyarray(py);
@@ -3588,6 +3569,26 @@ impl PauliProductRotation {
         Ok(gate.unbind())
     }
 
+    pub fn matrix(&self) -> Option<Array2<Complex64>> {
+        let angle = match self.angle {
+            Param::Float(f) => f,
+            _ => return None,
+        };
+
+        let pauli_mat = gate_matrix::pauli_zx_to_matrix(&self.z, &self.x);
+        let cos_val = (angle / 2.0).cos();
+        let sin_val = (angle / 2.0).sin();
+
+        let dim = pauli_mat.shape()[0];
+        let mut result = pauli_mat.mapv(|v| Complex64::new(0.0, -sin_val) * v);
+
+        for i in 0..dim {
+            result[[i, i]] += Complex64::new(cos_val, 0.0);
+        }
+
+        Some(result)
+    }
+
     /// Attempts to merge `self` and `other`.
     /// If successful, returns the merged [PauliProductRotation].
     /// If not successful, returns `None`.
@@ -3618,7 +3619,6 @@ impl PauliProductRotation {
             .count();
         let dim = 2u32.pow(num_qubits as u32);
         let tr_over_dim = if num_qubits == 0 {
-            // This is an identity Pauli rotation.
             (Complex64::new(0.0, -angle / 2.)).exp()
         } else {
             Complex64::new((angle / 2.).cos(), 0.)
@@ -3628,6 +3628,28 @@ impl PauliProductRotation {
     }
 }
 
+impl Operation for PauliProductRotation {
+    fn name(&self) -> &str {
+        "pauli_product_rotation"
+    }
+
+    fn num_qubits(&self) -> u32 {
+        self.z.len() as u32
+    }
+
+    fn num_clbits(&self) -> u32 {
+        0
+    }
+
+    fn num_params(&self) -> u32 {
+        1
+    }
+
+    fn directive(&self) -> bool {
+        false
+    }
+}
+
 impl PartialEq for PauliProductRotation {
     fn eq(&self, other: &Self) -> bool {
         self.x == other.x

crates/circuit/src/packed_instruction.rs

@@ -23,6 +23,7 @@ use crate::operations::{
     PyOperationTypes, PythonOperation, StandardGate, StandardInstruction, UnitaryGate,
 };
 use crate::{Block, Clbit, Qubit};
+
 use hashbrown::HashMap;
 use nalgebra::{Matrix2, Matrix4};
 use ndarray::{Array2, CowArray, Ix2};

crates/qpy/src/consts.rs

@@ -77,8 +77,8 @@ mod tests {
     #[test]
     fn test_name_coverage() {
         for gate in 1..STANDARD_GATE_SIZE as u8 {
-            let gate: StandardGate =
-                ::bytemuck::checked::try_cast::<_, StandardGate>(gate).unwrap();
+            let gate: StandardGate = ::bytemuck::checked::try_cast::<_, StandardGate>(gate)
+                .unwrap_or_else(|_| panic!("invalid StandardGate in test"));
             let gate_name = get_std_gate_class_name(&gate).clone();
             assert!(
                 standard_gate_from_gate_class_name(gate_name.as_str()).is_some(),