quantrs2_core/characterization/
gate_characterizer.rs

1//! Gate eigenstructure and characterization tools
2
3use crate::{
4    eigensolve::eigen_decompose_unitary,
5    error::{QuantRS2Error, QuantRS2Result},
6    gate::{single::*, GateOp},
7    qubit::QubitId,
8};
9use scirs2_core::ndarray::{Array1, Array2};
10use scirs2_core::Complex64 as Complex;
11
12/// Represents the eigenstructure of a quantum gate
13#[derive(Debug, Clone)]
14pub struct GateEigenstructure {
15    /// Eigenvalues of the gate
16    pub eigenvalues: Vec<Complex>,
17    /// Eigenvectors as columns of a matrix
18    pub eigenvectors: Array2<Complex>,
19    /// The original gate matrix
20    pub matrix: Array2<Complex>,
21}
22
23impl GateEigenstructure {
24    /// Check if the gate is diagonal in some basis
25    pub fn is_diagonal(&self, tolerance: f64) -> bool {
26        let n = self.matrix.nrows();
27        for i in 0..n {
28            for j in 0..n {
29                if i != j && self.matrix[(i, j)].norm() > tolerance {
30                    return false;
31                }
32            }
33        }
34        true
35    }
36
37    /// Get the phases of eigenvalues (assuming unitary gate)
38    pub fn eigenphases(&self) -> Vec<f64> {
39        self.eigenvalues
40            .iter()
41            .map(|&lambda| lambda.arg())
42            .collect()
43    }
44
45    /// Check if this represents a phase gate
46    pub fn is_phase_gate(&self, tolerance: f64) -> bool {
47        let magnitude = self.eigenvalues[0].norm();
48        self.eigenvalues
49            .iter()
50            .all(|&lambda| (lambda.norm() - magnitude).abs() < tolerance)
51    }
52
53    /// Get the rotation angle for single-qubit gates
54    pub fn rotation_angle(&self) -> Option<f64> {
55        if self.eigenvalues.len() != 2 {
56            return None;
57        }
58        let phase_diff = (self.eigenvalues[0] / self.eigenvalues[1]).arg();
59        Some(phase_diff.abs())
60    }
61
62    /// Get the rotation axis for single-qubit gates
63    pub fn rotation_axis(&self, tolerance: f64) -> Option<[f64; 3]> {
64        if self.eigenvalues.len() != 2 {
65            return None;
66        }
67
68        let v0 = self.eigenvectors.column(0);
69        let v1 = self.eigenvectors.column(1);
70
71        let bloch0 = eigenvector_to_bloch(&v0.to_owned());
72        let bloch1 = eigenvector_to_bloch(&v1.to_owned());
73
74        let axis = [
75            f64::midpoint(bloch0[0], bloch1[0]),
76            f64::midpoint(bloch0[1], bloch1[1]),
77            f64::midpoint(bloch0[2], bloch1[2]),
78        ];
79
80        let norm = axis[2]
81            .mul_add(axis[2], axis[0].mul_add(axis[0], axis[1] * axis[1]))
82            .sqrt();
83        if norm < tolerance {
84            None
85        } else {
86            Some([axis[0] / norm, axis[1] / norm, axis[2] / norm])
87        }
88    }
89}
90
91/// Convert an eigenvector to Bloch sphere coordinates
92fn eigenvector_to_bloch(v: &Array1<Complex>) -> [f64; 3] {
93    if v.len() != 2 {
94        return [0.0, 0.0, 0.0];
95    }
96
97    let v0 = v[0];
98    let v1 = v[1];
99    let rho00 = (v0 * v0.conj()).re;
100    let rho11 = (v1 * v1.conj()).re;
101    let rho01 = v0 * v1.conj();
102
103    [2.0 * rho01.re, -2.0 * rho01.im, rho00 - rho11]
104}
105
106/// Types of gates identified by characterization
107#[derive(Debug, Clone, PartialEq)]
108pub enum GateType {
109    /// Identity gate
110    Identity,
111    /// Pauli X gate
112    PauliX,
113    /// Pauli Y gate
114    PauliY,
115    /// Pauli Z gate
116    PauliZ,
117    /// Hadamard gate
118    Hadamard,
119    /// Phase gate
120    Phase { angle: f64 },
121    /// Rotation gate
122    Rotation { angle: f64, axis: [f64; 3] },
123    /// CNOT gate
124    CNOT,
125    /// Controlled phase gate
126    ControlledPhase { phase: f64 },
127    /// SWAP gate
128    SWAP,
129    /// General n-qubit gate
130    General { qubits: usize },
131}
132
133/// Gate characterization tools
134pub struct GateCharacterizer {
135    tolerance: f64,
136}
137
138impl GateCharacterizer {
139    /// Create a new gate characterizer
140    pub const fn new(tolerance: f64) -> Self {
141        Self { tolerance }
142    }
143
144    /// Compute the eigenstructure of a gate
145    pub fn eigenstructure(&self, gate: &dyn GateOp) -> QuantRS2Result<GateEigenstructure> {
146        let matrix_vec = gate.matrix()?;
147        let n = (matrix_vec.len() as f64).sqrt() as usize;
148
149        let mut matrix = Array2::zeros((n, n));
150        for i in 0..n {
151            for j in 0..n {
152                matrix[(i, j)] = matrix_vec[i * n + j];
153            }
154        }
155
156        let eigen = eigen_decompose_unitary(&matrix, self.tolerance)?;
157
158        Ok(GateEigenstructure {
159            eigenvalues: eigen.eigenvalues.to_vec(),
160            eigenvectors: eigen.eigenvectors,
161            matrix,
162        })
163    }
164
165    /// Identify the type of gate based on its eigenstructure
166    pub fn identify_gate_type(&self, gate: &dyn GateOp) -> QuantRS2Result<GateType> {
167        let eigen = self.eigenstructure(gate)?;
168        let n = eigen.eigenvalues.len();
169
170        match n {
171            2 => self.identify_single_qubit_gate(&eigen),
172            4 => self.identify_two_qubit_gate(&eigen),
173            _ => Ok(GateType::General {
174                qubits: (n as f64).log2() as usize,
175            }),
176        }
177    }
178
179    /// Identify single-qubit gate type
180    fn identify_single_qubit_gate(&self, eigen: &GateEigenstructure) -> QuantRS2Result<GateType> {
181        if self.is_pauli_gate(eigen) {
182            return Ok(self.identify_pauli_type(eigen));
183        }
184
185        if let Some(angle) = eigen.rotation_angle() {
186            if let Some(axis) = eigen.rotation_axis(self.tolerance) {
187                return Ok(GateType::Rotation { angle, axis });
188            }
189        }
190
191        if self.is_hadamard(eigen) {
192            return Ok(GateType::Hadamard);
193        }
194
195        Ok(GateType::General { qubits: 1 })
196    }
197
198    /// Identify two-qubit gate type
199    fn identify_two_qubit_gate(&self, eigen: &GateEigenstructure) -> QuantRS2Result<GateType> {
200        if self.is_cnot(eigen) {
201            return Ok(GateType::CNOT);
202        }
203
204        if self.is_controlled_phase(eigen) {
205            if let Some(phase) = self.extract_controlled_phase(eigen) {
206                return Ok(GateType::ControlledPhase { phase });
207            }
208        }
209
210        if self.is_swap_variant(eigen) {
211            return Ok(self.identify_swap_type(eigen));
212        }
213
214        Ok(GateType::General { qubits: 2 })
215    }
216
217    /// Check if gate is a Pauli gate
218    fn is_pauli_gate(&self, eigen: &GateEigenstructure) -> bool {
219        eigen.eigenvalues.iter().all(|&lambda| {
220            let is_plus_minus_one = (lambda - Complex::new(1.0, 0.0)).norm() < self.tolerance
221                || (lambda + Complex::new(1.0, 0.0)).norm() < self.tolerance;
222            let is_plus_minus_i = (lambda - Complex::new(0.0, 1.0)).norm() < self.tolerance
223                || (lambda + Complex::new(0.0, 1.0)).norm() < self.tolerance;
224            is_plus_minus_one || is_plus_minus_i
225        })
226    }
227
228    /// Identify which Pauli gate
229    fn identify_pauli_type(&self, eigen: &GateEigenstructure) -> GateType {
230        let matrix = &eigen.matrix;
231        let tolerance = self.tolerance;
232
233        if (matrix[(0, 1)] - Complex::new(1.0, 0.0)).norm() < tolerance
234            && (matrix[(1, 0)] - Complex::new(1.0, 0.0)).norm() < tolerance
235            && matrix[(0, 0)].norm() < tolerance
236            && matrix[(1, 1)].norm() < tolerance
237        {
238            GateType::PauliX
239        } else if (matrix[(0, 1)] - Complex::new(0.0, -1.0)).norm() < tolerance
240            && (matrix[(1, 0)] - Complex::new(0.0, 1.0)).norm() < tolerance
241            && matrix[(0, 0)].norm() < tolerance
242            && matrix[(1, 1)].norm() < tolerance
243        {
244            GateType::PauliY
245        } else if (matrix[(0, 0)] - Complex::new(1.0, 0.0)).norm() < tolerance
246            && (matrix[(1, 1)] - Complex::new(-1.0, 0.0)).norm() < tolerance
247            && matrix[(0, 1)].norm() < tolerance
248            && matrix[(1, 0)].norm() < tolerance
249        {
250            GateType::PauliZ
251        } else {
252            GateType::General { qubits: 1 }
253        }
254    }
255
256    /// Check if gate is Hadamard
257    fn is_hadamard(&self, eigen: &GateEigenstructure) -> bool {
258        eigen.eigenvalues.iter().all(|&lambda| {
259            (lambda - Complex::new(1.0, 0.0)).norm() < self.tolerance
260                || (lambda + Complex::new(1.0, 0.0)).norm() < self.tolerance
261        })
262    }
263
264    /// Check if gate is CNOT
265    fn is_cnot(&self, eigen: &GateEigenstructure) -> bool {
266        eigen.eigenvalues.len() == 4
267            && eigen.eigenvalues.iter().all(|&lambda| {
268                (lambda - Complex::new(1.0, 0.0)).norm() < self.tolerance
269                    || (lambda + Complex::new(1.0, 0.0)).norm() < self.tolerance
270            })
271    }
272
273    /// Check if gate is a controlled phase gate
274    fn is_controlled_phase(&self, eigen: &GateEigenstructure) -> bool {
275        let ones = eigen
276            .eigenvalues
277            .iter()
278            .filter(|&&lambda| (lambda - Complex::new(1.0, 0.0)).norm() < self.tolerance)
279            .count();
280        ones == 3
281    }
282
283    /// Extract phase from controlled phase gate
284    fn extract_controlled_phase(&self, eigen: &GateEigenstructure) -> Option<f64> {
285        eigen
286            .eigenvalues
287            .iter()
288            .find(|&&lambda| (lambda - Complex::new(1.0, 0.0)).norm() > self.tolerance)
289            .map(|&lambda| lambda.arg())
290    }
291
292    /// Check if gate is a SWAP variant
293    fn is_swap_variant(&self, eigen: &GateEigenstructure) -> bool {
294        let ones = eigen
295            .eigenvalues
296            .iter()
297            .filter(|&&lambda| (lambda - Complex::new(1.0, 0.0)).norm() < self.tolerance)
298            .count();
299        ones >= 2
300    }
301
302    /// Identify SWAP variant type
303    fn identify_swap_type(&self, eigen: &GateEigenstructure) -> GateType {
304        let matrix = &eigen.matrix;
305
306        if (matrix[(0, 0)] - Complex::new(1.0, 0.0)).norm() < self.tolerance
307            && (matrix[(1, 2)] - Complex::new(1.0, 0.0)).norm() < self.tolerance
308            && (matrix[(2, 1)] - Complex::new(1.0, 0.0)).norm() < self.tolerance
309            && (matrix[(3, 3)] - Complex::new(1.0, 0.0)).norm() < self.tolerance
310            && matrix[(0, 1)].norm() < self.tolerance
311            && matrix[(0, 2)].norm() < self.tolerance
312            && matrix[(0, 3)].norm() < self.tolerance
313            && matrix[(1, 0)].norm() < self.tolerance
314            && matrix[(1, 1)].norm() < self.tolerance
315            && matrix[(1, 3)].norm() < self.tolerance
316            && matrix[(2, 0)].norm() < self.tolerance
317            && matrix[(2, 2)].norm() < self.tolerance
318            && matrix[(2, 3)].norm() < self.tolerance
319            && matrix[(3, 0)].norm() < self.tolerance
320            && matrix[(3, 1)].norm() < self.tolerance
321            && matrix[(3, 2)].norm() < self.tolerance
322        {
323            GateType::SWAP
324        } else {
325            GateType::General { qubits: 2 }
326        }
327    }
328
329    /// Compare two matrices for equality
330    #[allow(dead_code)]
331    fn matrix_equals(a: &Array2<Complex>, b: &Array2<Complex>, tolerance: f64) -> bool {
332        a.shape() == b.shape()
333            && a.iter()
334                .zip(b.iter())
335                .all(|(a_ij, b_ij)| (a_ij - b_ij).norm() < tolerance)
336    }
337
338    /// Decompose a gate into rotation gates based on eigenstructure
339    pub fn decompose_to_rotations(
340        &self,
341        gate: &dyn GateOp,
342    ) -> QuantRS2Result<Vec<Box<dyn GateOp>>> {
343        let eigen = self.eigenstructure(gate)?;
344
345        match eigen.eigenvalues.len() {
346            2 => Self::decompose_single_qubit(&eigen),
347            _ => Err(QuantRS2Error::UnsupportedOperation(
348                "Rotation decomposition only supported for single-qubit gates".to_string(),
349            )),
350        }
351    }
352
353    /// Decompose single-qubit gate
354    fn decompose_single_qubit(eigen: &GateEigenstructure) -> QuantRS2Result<Vec<Box<dyn GateOp>>> {
355        let matrix = &eigen.matrix;
356
357        let alpha = matrix[(1, 1)].arg() - matrix[(1, 0)].arg();
358        let gamma = matrix[(1, 1)].arg() + matrix[(1, 0)].arg();
359        let beta = 2.0 * matrix[(1, 0)].norm().acos();
360
361        Ok(vec![
362            Box::new(RotationZ {
363                target: QubitId(0),
364                theta: alpha,
365            }),
366            Box::new(RotationY {
367                target: QubitId(0),
368                theta: beta,
369            }),
370            Box::new(RotationZ {
371                target: QubitId(0),
372                theta: gamma,
373            }),
374        ])
375    }
376
377    /// Find the closest Clifford gate to a given gate
378    pub fn find_closest_clifford(&self, gate: &dyn GateOp) -> QuantRS2Result<Box<dyn GateOp>> {
379        let eigen = self.eigenstructure(gate)?;
380
381        if eigen.eigenvalues.len() != 2 {
382            return Err(QuantRS2Error::UnsupportedOperation(
383                "Clifford approximation only supported for single-qubit gates".to_string(),
384            ));
385        }
386
387        let target = QubitId(0);
388        let clifford_gates: Vec<Box<dyn GateOp>> = vec![
389            Box::new(PauliX { target }),
390            Box::new(PauliY { target }),
391            Box::new(PauliZ { target }),
392            Box::new(Hadamard { target }),
393            Box::new(Phase { target }),
394        ];
395
396        let mut min_distance = f64::INFINITY;
397        let mut closest_gate = None;
398
399        for clifford in &clifford_gates {
400            let distance = self.gate_distance(gate, clifford.as_ref())?;
401            if distance < min_distance {
402                min_distance = distance;
403                closest_gate = Some(clifford.clone());
404            }
405        }
406
407        closest_gate.ok_or_else(|| {
408            QuantRS2Error::ComputationError("Failed to find closest Clifford gate".to_string())
409        })
410    }
411
412    /// Compute the distance between two gates (Frobenius norm)
413    pub fn gate_distance(&self, gate1: &dyn GateOp, gate2: &dyn GateOp) -> QuantRS2Result<f64> {
414        let m1_vec = gate1.matrix()?;
415        let m2_vec = gate2.matrix()?;
416
417        if m1_vec.len() != m2_vec.len() {
418            return Err(QuantRS2Error::InvalidInput(
419                "Gates must have the same dimensions".to_string(),
420            ));
421        }
422
423        let diff_sqr: f64 = m1_vec
424            .iter()
425            .zip(m2_vec.iter())
426            .map(|(a, b)| (a - b).norm_sqr())
427            .sum();
428        Ok(diff_sqr.sqrt())
429    }
430
431    /// Check if a gate is approximately equal to identity
432    pub fn is_identity(&self, gate: &dyn GateOp, tolerance: f64) -> bool {
433        let matrix_vec = match gate.matrix() {
434            Ok(m) => m,
435            Err(_) => return false,
436        };
437        let n = (matrix_vec.len() as f64).sqrt() as usize;
438
439        for i in 0..n {
440            for j in 0..n {
441                let idx = i * n + j;
442                let expected = if i == j {
443                    Complex::new(1.0, 0.0)
444                } else {
445                    Complex::new(0.0, 0.0)
446                };
447                if (matrix_vec[idx] - expected).norm() > tolerance {
448                    return false;
449                }
450            }
451        }
452        true
453    }
454
455    /// Extract the global phase of a gate
456    pub fn global_phase(&self, gate: &dyn GateOp) -> QuantRS2Result<f64> {
457        let eigen = self.eigenstructure(gate)?;
458
459        let det = eigen
460            .eigenvalues
461            .iter()
462            .fold(Complex::new(1.0, 0.0), |acc, &lambda| acc * lambda);
463        let n = eigen.eigenvalues.len() as f64;
464        Ok(det.arg() / n)
465    }
466}
quantrs2_core/characterization/gate_characterizer.rs

quantrs2_core/characterization/
gate_characterizer.rs
