paulistrings/channel/ptm/
unitary.rs

1//! The general-unitary constructors of [`PtmChannel`], whose PTM is derived from the matrix.
2
3use super::{assert_distinct, sparse_storage, PtmChannel, PtmKind, ONE, ZERO};
4use num_complex::Complex64;
5
6/// Tolerance below which a derived PTM entry is treated as an exact zero, since every nonzero entry becomes an output term.
7const PTM_EPS: f64 = 1e-12;
8
9/// The largest `|(U†U − I)_ij|` the general-unitary constructors accept; the Python bindings check the same bound first.
10const UNITARITY_TOLERANCE: f64 = 1e-10;
11
12impl PtmChannel {
13    /// A general single-qubit unitary `U`: `ptm[s][t] = tr(P_t · U P_s U†) / 2`.
14    ///
15    /// # Panics
16    ///
17    /// Unless every entry of `U†U − I` is at most `1e-10` in absolute value.
18    pub fn unitary_1q(qubit: usize, u: [[Complex64; 2]; 2]) -> Self {
19        assert_unitary("unitary_1q", &u);
20        let u_dagger = dagger(&u);
21        let mut entries = Vec::new();
22        for s in 0..4 {
23            let conjugated = matmul(&matmul(&u, &pauli_matrix(s)), &u_dagger);
24            for t in 0..4 {
25                let value =
26                    trace(&matmul(&pauli_matrix(t), &conjugated)) / Complex64::new(2.0, 0.0);
27                entries.push(((s as u8, t as u8), clean(value)));
28            }
29        }
30        Self::new(&[qubit], sparse_storage(entries), PtmKind::Unitary1q)
31    }
32
33    /// A general two-qubit unitary `U` on `|q0 q1⟩`, `q0` the more significant tensor factor: `ptm[s][t] = tr(P_t · U P_s U†) / 4`.
34    ///
35    /// # Panics
36    ///
37    /// If `q0 == q1`, or unless every entry of `U†U − I` is at most `1e-10` in absolute value.
38    pub fn unitary_2q(q0: usize, q1: usize, u: [[Complex64; 4]; 4]) -> Self {
39        assert_distinct("unitary_2q", q0, q1);
40        assert_unitary("unitary_2q", &u);
41        let u_dagger = dagger(&u);
42        let two_qubit_pauli = |s: usize| kron(&pauli_matrix(s & 3), &pauli_matrix(s >> 2));
43        let mut entries = Vec::new();
44        for s in 0..16 {
45            let conjugated = matmul(&matmul(&u, &two_qubit_pauli(s)), &u_dagger);
46            for t in 0..16 {
47                let value =
48                    trace(&matmul(&two_qubit_pauli(t), &conjugated)) / Complex64::new(4.0, 0.0);
49                entries.push(((s as u8, t as u8), clean(value)));
50            }
51        }
52        Self::new(&[q0, q1], sparse_storage(entries), PtmKind::Unitary2q)
53    }
54}
55
56/// Panics unless `u` is unitary to [`UNITARITY_TOLERANCE`]; a NaN entry fails too.
57fn assert_unitary<const N: usize>(constructor: &str, u: &[[Complex64; N]; N]) {
58    let product = matmul(&dagger(u), u);
59    let deviation = (0..N)
60        .flat_map(|i| (0..N).map(move |j| (i, j)))
61        .map(|(i, j)| {
62            let identity = if i == j { ONE } else { ZERO };
63            (product[i][j] - identity).norm()
64        })
65        .fold(
66            0.0,
67            |largest: f64, entry| {
68                if entry > largest || entry.is_nan() {
69                    entry
70                } else {
71                    largest
72                }
73            },
74        );
75    assert!(
76        deviation <= UNITARITY_TOLERANCE,
77        "PtmChannel::{constructor}: the matrix is not unitary (largest entry of U†U − I is {deviation:e}, tolerance {UNITARITY_TOLERANCE:e})"
78    );
79}
80
81/// The four single-qubit Pauli matrices, in `I, X, Z, Y` index order.
82fn pauli_matrix(index: usize) -> [[Complex64; 2]; 2] {
83    let i = Complex64::new(0.0, 1.0);
84    match index {
85        0 => [[ONE, ZERO], [ZERO, ONE]],
86        1 => [[ZERO, ONE], [ONE, ZERO]],
87        2 => [[ONE, ZERO], [ZERO, -ONE]],
88        _ => [[ZERO, -i], [i, ZERO]],
89    }
90}
91
92fn matmul<const N: usize>(a: &[[Complex64; N]; N], b: &[[Complex64; N]; N]) -> [[Complex64; N]; N] {
93    let mut out = [[ZERO; N]; N];
94    for i in 0..N {
95        for k in 0..N {
96            let aik = a[i][k];
97            if aik == ZERO {
98                continue;
99            }
100            for j in 0..N {
101                out[i][j] += aik * b[k][j];
102            }
103        }
104    }
105    out
106}
107
108fn dagger<const N: usize>(a: &[[Complex64; N]; N]) -> [[Complex64; N]; N] {
109    core::array::from_fn(|i| core::array::from_fn(|j| a[j][i].conj()))
110}
111
112fn trace<const N: usize>(a: &[[Complex64; N]; N]) -> Complex64 {
113    a.iter().enumerate().map(|(i, row)| row[i]).sum()
114}
115
116/// `kron(a, b)` with `a` acting on the more significant qubit.
117fn kron(a: &[[Complex64; 2]; 2], b: &[[Complex64; 2]; 2]) -> [[Complex64; 4]; 4] {
118    core::array::from_fn(|row| {
119        core::array::from_fn(|column| a[row / 2][column / 2] * b[row % 2][column % 2])
120    })
121}
122
123/// Each part of `v` within [`PTM_EPS`] of zero rounded to exactly zero.
124fn clean(v: Complex64) -> Complex64 {
125    let snap = |part: f64| if part.abs() < PTM_EPS { 0.0 } else { part };
126    Complex64::new(snap(v.re), snap(v.im))
127}
128
129#[cfg(test)]
130mod tests;
paulistrings/channel/ptm/unitary.rs

paulistrings/channel/ptm/
unitary.rs
