diff --git a/crates/circuit/Cargo.toml b/crates/circuit/Cargo.toml index 204f1eede5d0..b0e439888892 100644 --- a/crates/circuit/Cargo.toml +++ b/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"] } \ No newline at end of file diff --git a/crates/circuit/src/circuit_drawer.rs b/crates/circuit/src/circuit_drawer.rs index 45ba2306e05f..2fb5e8380db3 100644 --- a/crates/circuit/src/circuit_drawer.rs +++ b/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::>(); let params = (0..num_params) - .map(|_x| 3.141.into()) + .map(|_x| std::f64::consts::PI.into()) .collect::>(); circuit.push_standard_gate(op, ¶ms, &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::>(); + let qubits = (0..i).map(Qubit::new).collect::>(); 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(); } diff --git a/crates/circuit/src/gate_matrix.rs b/crates/circuit/src/gate_matrix.rs index f361b204eab2..b7a4608de09e 100644 --- a/crates/circuit/src/gate_matrix.rs +++ b/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 { + 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, b: &Array2) -> Array2 { + let (ar, ac) = a.dim(); + let (br, bc) = b.dim(); + let mut out = Array2::::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 +} diff --git a/crates/circuit/src/operations.rs b/crates/circuit/src/operations.rs index 082797e9a8a3..5fb650f81234 100644 --- a/crates/circuit/src/operations.rs +++ b/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> { - Python::attach(|py| -> Option> { - match self.instruction.getattr(py, intern!(py, "to_matrix")) { - Ok(to_matrix) => { - let res: Option> = to_matrix.call0(py).ok()?.extract(py).ok(); - match res { - Some(x) => { - let array: PyReadonlyArray2 = x.extract(py).ok()?; - Some(array.as_array().to_owned()) - } - None => None, - } - } - Err(_) => None, - } - }) - } - pub fn definition(&self) -> Option { Python::attach(|py| -> Option { match self.instruction.getattr(py, intern!(py, "definition")) { @@ -3311,6 +3292,24 @@ impl PyInstruction { }) } + pub fn matrix(&self) -> Option> { + Python::attach(|py| -> Option> { + match self.instruction.getattr(py, intern!(py, "to_matrix")) { + Ok(to_matrix) => { + let res: Option> = to_matrix.call0(py).ok()?.extract(py).ok(); + match res { + Some(x) => { + let array: PyReadonlyArray2 = 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> { let z = self.z.to_pyarray(py); @@ -3588,6 +3569,26 @@ impl PauliProductRotation { Ok(gate.unbind()) } + pub fn matrix(&self) -> Option> { + 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 diff --git a/crates/circuit/src/packed_instruction.rs b/crates/circuit/src/packed_instruction.rs index 0577a0fd49fa..809fcaca5861 100644 --- a/crates/circuit/src/packed_instruction.rs +++ b/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}; diff --git a/crates/qpy/src/consts.rs b/crates/qpy/src/consts.rs index b7442847965e..947421168ca2 100644 --- a/crates/qpy/src/consts.rs +++ b/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(), diff --git a/crates/quantum_info/src/sparse_pauli_op.rs b/crates/quantum_info/src/sparse_pauli_op.rs index f26140961356..d7379f5e1229 100644 --- a/crates/quantum_info/src/sparse_pauli_op.rs +++ b/crates/quantum_info/src/sparse_pauli_op.rs @@ -992,6 +992,56 @@ fn to_matrix_dense_inner(paulis: &MatrixCompressedPaulis, parallel: bool) -> Vec out } +pub fn pauli_zx_to_matrix(z: &[bool], x: &[bool]) -> ndarray::Array2 { + use ndarray::Array2; + use num_complex::Complex64; + + assert_eq!(z.len(), x.len()); + + let n = z.len(); + let dim = 1usize << n; + let mut out = Array2::::zeros((dim, dim)); + + for col in 0..dim { + let mut row = col; + let mut phase = Complex64::new(1.0, 0.0); + + for i in 0..n { + let bit_index = n - 1 - i; + let bit = (col >> bit_index) & 1; + + match (z[i], x[i]) { + (false, false) => { + // I + } + (false, true) => { + // X + row ^= 1 << bit_index; + } + (true, false) => { + // Z + if bit == 1 { + phase = -phase; + } + } + (true, true) => { + // Y = iXZ + row ^= 1 << bit_index; + phase *= if bit == 0 { + Complex64::new(0.0, 1.0) + } else { + Complex64::new(0.0, -1.0) + }; + } + } + } + + out[(row, col)] = phase; + } + + out +} + type CSRData = (Vec, Vec, Vec); type ToCSRData = fn(&MatrixCompressedPaulis) -> CSRData; diff --git a/crates/transpiler/src/passes/split_2q_unitaries.rs b/crates/transpiler/src/passes/split_2q_unitaries.rs index a65d53ac0908..10939afe57b8 100644 --- a/crates/transpiler/src/passes/split_2q_unitaries.rs +++ b/crates/transpiler/src/passes/split_2q_unitaries.rs @@ -175,7 +175,7 @@ pub fn run_split_2q_unitaries( inst.params.as_deref().cloned(), inst.label.as_ref().map(|x| x.to_string()), #[cfg(feature = "cache_pygates")] - inst.py_op.get().map(|x| x.clone()), + inst.py_op.get().cloned(), )?; } Ok(Some((new_dag.build(), mapping))) diff --git a/releasenotes/notes/ppr-rust-matrix-support-15869.rst b/releasenotes/notes/ppr-rust-matrix-support-15869.rst new file mode 100644 index 000000000000..3c7ed529d820 --- /dev/null +++ b/releasenotes/notes/ppr-rust-matrix-support-15869.rst @@ -0,0 +1,6 @@ +--- +features: + - | + Added Rust-side matrix support for ``PauliProductRotationGate`` by reusing + the existing Pauli matrix construction logic from the Rust + ``SparsePauliOp`` implementation. \ No newline at end of file diff --git a/test/python/circuit/library/test_pauli_product_rotation.py b/test/python/circuit/library/test_pauli_product_rotation.py index 3545052a6ffc..932c81784f69 100644 --- a/test/python/circuit/library/test_pauli_product_rotation.py +++ b/test/python/circuit/library/test_pauli_product_rotation.py @@ -24,7 +24,7 @@ PauliEvolutionGate, PauliProductMeasurement, ) -from qiskit.quantum_info import Pauli, Operator, SparsePauliOp +from qiskit.quantum_info import Pauli, Operator, SparsePauliOp, Statevector from qiskit.qpy import dump, load from qiskit.compiler import transpile from test import QiskitTestCase # pylint: disable=wrong-import-order @@ -223,6 +223,119 @@ def test_matrix_conventions(self): evo = PauliEvolutionGate(pauli, time=angle / 2) self.assertTrue(np.allclose(evo.to_matrix(), ppr.to_matrix())) + def test_single_qubit_x_rotation_matrix(self): + """Test matrix support for a single-qubit X rotation.""" + theta = np.pi / 3 + gate = PauliProductRotationGate(Pauli("X"), theta) + mat = np.asarray(Operator(gate).data) + + expected = np.array( + [ + [np.cos(theta / 2), -1j * np.sin(theta / 2)], + [-1j * np.sin(theta / 2), np.cos(theta / 2)], + ], + dtype=complex, + ) + + self.assertEqual(mat.shape, (2, 2)) + self.assertTrue(np.allclose(mat, expected)) + + def test_single_qubit_z_rotation_matrix(self): + """Test matrix support for a single-qubit Z rotation.""" + theta = np.pi / 4 + gate = PauliProductRotationGate(Pauli("Z"), theta) + mat = np.asarray(Operator(gate).data) + + expected = np.array( + [ + [np.exp(-1j * theta / 2), 0], + [0, np.exp(1j * theta / 2)], + ], + dtype=complex, + ) + + self.assertEqual(mat.shape, (2, 2)) + self.assertTrue(np.allclose(mat, expected)) + + def test_two_qubit_zz_matrix(self): + """Test matrix support for a two-qubit ZZ rotation.""" + theta = np.pi / 5 + gate = PauliProductRotationGate(Pauli("ZZ"), theta) + mat = np.asarray(Operator(gate).data) + + self.assertEqual(mat.shape, (4, 4)) + self.assertTrue(np.allclose(mat.conj().T @ mat, np.eye(4))) + + def test_three_qubit_xyz_matrix(self): + """Test matrix support for a three-qubit XYZ rotation.""" + theta = 0.7 + gate = PauliProductRotationGate(Pauli("XYZ"), theta) + mat = np.asarray(Operator(gate).data) + + self.assertEqual(mat.shape, (8, 8)) + self.assertTrue(np.allclose(mat.conj().T @ mat, np.eye(8))) + + def test_identity_pauli_is_global_phase(self): + """Test that an all-identity Pauli only contributes a global phase.""" + theta = np.pi / 6 + gate = PauliProductRotationGate(Pauli("II"), theta) + mat = np.asarray(Operator(gate).data) + expected = np.exp(-1j * theta / 2) * np.eye(4) + + self.assertTrue(np.allclose(mat, expected)) + + def test_theta_zero_is_identity(self): + """Test that zero angle gives the identity.""" + gate = PauliProductRotationGate(Pauli("XZ"), 0.0) + mat = np.asarray(Operator(gate).data) + + self.assertTrue(np.allclose(mat, np.eye(4))) + + def test_theta_2pi_is_negative_identity(self): + """Test that a 2π rotation gives -I up to global phase.""" + gate = PauliProductRotationGate(Pauli("Z"), 2 * np.pi) + mat = np.asarray(Operator(gate).data) + + self.assertTrue(np.allclose(mat, -np.eye(2))) + + def test_matrix_consistent_with_simulation(self): + """Test matrix action is consistent with statevector simulation.""" + theta = np.pi / 3 + gate = PauliProductRotationGate(Pauli("ZZ"), theta) + + qc = QuantumCircuit(2) + qc.append(gate, [0, 1]) + + psi0 = Statevector.from_label("00") + evolved = psi0.evolve(qc) + mat = np.asarray(Operator(gate).data) + expected = Statevector(mat @ psi0.data) + + self.assertTrue(np.allclose(evolved.data, expected.data)) + + @data( + ("X", 1), + ("Y", 1), + ("Z", 1), + ("XX", 2), + ("YY", 2), + ("ZZ", 2), + ("XY", 2), + ("YZ", 2), + ("XXX", 3), + ("ZZZ", 3), + ) + def test_matrix_unitary_for_common_paulis(self, case): + """Test matrix support for common Pauli strings.""" + pauli, n_qubits = case + theta = 1.23 + + gate = PauliProductRotationGate(Pauli(pauli), theta) + mat = np.asarray(Operator(gate).data) + + self.assertEqual(mat.shape, (2**n_qubits, 2**n_qubits)) + self.assertTrue(np.allclose(mat.conj().T @ mat, np.eye(2**n_qubits))) + @data(0, 1, 2) def test_control(self, num_ctrl_qubits): """Check calling the control method."""