Merge upstream/main and finalize PPR matrix fixes

Author: Adithyaphani

Cargo.lock

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

crates/circuit/Cargo.toml

@@ -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/circuit_drawer.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/gate_matrix.rs

@@ -3531,7 +3531,7 @@ impl UnitaryGate {
     }
 }
 
-/// A Pauli-based gate model, consisting of PauliProductRotation and PauliProductMeasurement ops.
+/// A Pauli-based gate model, consisting of [PauliProductRotation] and [PauliProductMeasurement] ops.
 #[derive(Clone, Debug, PartialEq, Eq)]
 pub enum PauliBased {
     PauliProductRotation(PauliProductRotation),
@@ -3568,20 +3568,63 @@ 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_dense_matrix(&self.z, &self.x);
+
+        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 identity = Array2::<Complex64>::eye(dim);
-        Some(
-            identity.mapv(|v| Complex64::new(v.re * cos_val, 0.))
-                + pauli_mat.mapv(|v| Complex64::new(0., -sin_val) * v),
-        )
+        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`.
+    pub fn merge_with(&self, other: &Self) -> Option<Self> {
+        if self.x == other.x && self.z == other.z {
+            Some(PauliProductRotation {
+                z: self.z.clone(),
+                x: self.x.clone(),
+                angle: radd_param(self.angle.clone(), other.angle.clone()),
+            })
+        } else {
+            None
+        }
+    }
+
+    /// For a [PauliProductRotation] gate with a floating-point angle return a tuple `(Tr(gate) / dim, dim)`.
+    /// Return `None` if the angle is parameterized.
+    pub fn rotation_trace_and_dim(&self) -> Option<(Complex64, f64)> {
+        let Param::Float(angle) = self.angle else {
+            return None;
+        };
+
+        let num_qubits = self
+            .z
+            .iter()
+            .zip(self.x.iter())
+            .filter(|(z, x)| **z || **x)
+            .count();
+        let dim = 2u32.pow(num_qubits as u32);
+        let tr_over_dim = if num_qubits == 0 {
+            (Complex64::new(0.0, -angle / 2.)).exp()
+        } else {
+            Complex64::new((angle / 2.).cos(), 0.)
+        };
+
+        Some((tr_over_dim, dim as f64))
     }
 }
 

crates/circuit/src/operations.rs

@@ -47,3 +47,12 @@ pub const C_ONE: Complex64 = c64(1., 0.);
 pub const C_M_ONE: Complex64 = c64(-1., 0.);
 pub const IM: Complex64 = c64(0., 1.);
 pub const M_IM: Complex64 = c64(0., -1.);
+
+use std::f64::consts::{FRAC_PI_2, FRAC_PI_4, FRAC_PI_8, SQRT_2};
+pub const C_FRAC_PI_2: Complex64 = c64(FRAC_PI_2, 0.0);
+pub const C_FRAC_PI_4: Complex64 = c64(FRAC_PI_4, 0.0);
+pub const C_FRAC_PI_8: Complex64 = c64(FRAC_PI_8, 0.0);
+pub const C_FRAC_PI_2_SQRT_2: Complex64 = c64(FRAC_PI_2 / SQRT_2, 0.0);
+pub const C_M_FRAC_PI_4: Complex64 = c64(-FRAC_PI_4, 0.0);
+pub const C_M_FRAC_PI_8: Complex64 = c64(-FRAC_PI_8, 0.0);
+pub const C_M_FRAC_PI_2_SQRT_2: Complex64 = c64(-FRAC_PI_2 / SQRT_2, 0.0);

crates/circuit/src/util.rs

@@ -72,6 +72,7 @@ fn _accelerate(m: &Bound<PyModule>) -> PyResult<()> {
     add_submodule(m, ::qiskit_transpiler::passes::sabre::sabre, "sabre")?;
     add_submodule(m, ::qiskit_accelerate::sampled_exp_val::sampled_exp_val, "sampled_exp_val")?;
     add_submodule(m, ::qiskit_quantum_info::sparse_observable::sparse_observable, "sparse_observable")?;
+    add_submodule(m, ::qiskit_quantum_info::sparse_observable::standard_generators::standard_generators, "standard_generators")?;
     add_submodule(m, ::qiskit_quantum_info::sparse_pauli_op::sparse_pauli_op, "sparse_pauli_op")?;
     add_submodule(m, ::qiskit_transpiler::passes::scheduling_mod, "scheduling")?;
     add_submodule(m, ::qiskit_quantum_info::unitary_sim::unitary_sim, "unitary_sim")?;

crates/pyext/src/lib.rs

@@ -21,6 +21,7 @@ qiskit-circuit.workspace = true
 qiskit-quantum-info.workspace = true
 num-bigint.workspace = true
 num-complex.workspace = true
+num-traits.workspace = true
 bytemuck.workspace = true
 binrw.workspace = true
 uuid = {version = "1", features = ["v4"]}