paulistrings/channel/
ptm.rs

1//! [`PtmChannel`], every built-in one- and two-qubit gate and noise channel as a compact Pauli-transfer matrix (PTM).
2
3mod clifford;
4mod noise;
5mod unitary;
6
7use super::prepared::{PrepareToken, Prepared, LOCAL_DIM};
8use super::{qubit_loc, read_pauli, support_mask, write_pauli, Channel, OutputBuffer};
9use crate::pauli_sum::hash::Gf2Hash;
10use num_complex::Complex64;
11
12const ZERO: Complex64 = Complex64::new(0.0, 0.0);
13const ONE: Complex64 = Complex64::new(1.0, 0.0);
14
15/// A one- or two-qubit gate or noise channel, stored as its Pauli-transfer matrix.
16///
17/// Built by the named constructors: the Cliffords [`hadamard`](Self::hadamard), [`phase_gate`](Self::phase_gate), [`pauli_x`](Self::pauli_x), [`pauli_y`](Self::pauli_y), [`pauli_z`](Self::pauli_z), [`cnot`](Self::cnot), [`cz`](Self::cz), [`swap`](Self::swap);
18/// the general unitaries [`unitary_1q`](Self::unitary_1q), [`unitary_2q`](Self::unitary_2q);
19/// and the noise channels [`depolarizing`](Self::depolarizing), [`dephasing`](Self::dephasing), [`pauli_channel`](Self::pauli_channel), [`depolarizing_2q`](Self::depolarizing_2q), [`amplitude_damping`](Self::amplitude_damping).
20///
21/// On the support, Paulis are indexed `x | (z << 1)` per qubit (`I = 0, X = 1, Z = 2, Y = 3`), the constructor's first qubit in the low two bits.
22/// [`Channel::apply`] maps input Pauli `s` to `Σ_t ptm[s][t] · P_t`, and [`Channel::apply_adjoint`] reads the transpose.
23/// [`Channel::debug_name`] is the constructor's name (`"hadamard"`, `"cnot"`, …), with `_adjoint` appended after [`adjoint`](Self::adjoint).
24#[derive(Clone, Debug)]
25pub struct PtmChannel {
26    /// Support qubits in constructor order; only the first `num_support_qubits` are meaningful.
27    qubits: [u32; 2],
28    num_support_qubits: u8,
29    kind: PtmKind,
30    /// Whether this is the [`adjoint`](Self::adjoint) of the constructed channel; [`Channel::debug_name`] appends `_adjoint`.
31    adjoint: bool,
32    storage: PtmStorage,
33}
34
35/// The constructor a [`PtmChannel`] came from, which [`Channel::debug_name`] reports.
36#[derive(Clone, Copy, Debug, PartialEq, Eq)]
37enum PtmKind {
38    Hadamard,
39    PhaseGate,
40    PauliX,
41    PauliY,
42    PauliZ,
43    Cnot,
44    Cz,
45    Swap,
46    Unitary1q,
47    Unitary2q,
48    Depolarizing,
49    Dephasing,
50    PauliChannel,
51    Depolarizing2q,
52    AmplitudeDamping,
53}
54
55impl PtmKind {
56    /// The constructor's name, with `_adjoint` appended for the adjoint.
57    fn name(self, adjoint: bool) -> &'static str {
58        let [plain, adjoint_name] = match self {
59            Self::Hadamard => ["hadamard", "hadamard_adjoint"],
60            Self::PhaseGate => ["phase_gate", "phase_gate_adjoint"],
61            Self::PauliX => ["pauli_x", "pauli_x_adjoint"],
62            Self::PauliY => ["pauli_y", "pauli_y_adjoint"],
63            Self::PauliZ => ["pauli_z", "pauli_z_adjoint"],
64            Self::Cnot => ["cnot", "cnot_adjoint"],
65            Self::Cz => ["cz", "cz_adjoint"],
66            Self::Swap => ["swap", "swap_adjoint"],
67            Self::Unitary1q => ["unitary_1q", "unitary_1q_adjoint"],
68            Self::Unitary2q => ["unitary_2q", "unitary_2q_adjoint"],
69            Self::Depolarizing => ["depolarizing", "depolarizing_adjoint"],
70            Self::Dephasing => ["dephasing", "dephasing_adjoint"],
71            Self::PauliChannel => ["pauli_channel", "pauli_channel_adjoint"],
72            Self::Depolarizing2q => ["depolarizing_2q", "depolarizing_2q_adjoint"],
73            Self::AmplitudeDamping => ["amplitude_damping", "amplitude_damping_adjoint"],
74        };
75        if adjoint {
76            adjoint_name
77        } else {
78            plain
79        }
80    }
81}
82
83/// The nonzero PTM entries, in the cheapest form the constructor knows to be exact.
84#[derive(Clone, Debug)]
85enum PtmStorage {
86    /// A Clifford: input `s` goes to `image[s]` with coefficient `-1` where bit `s` of `negated` is set, `+1` otherwise.
87    SignedPermutation { image: [u8; LOCAL_DIM], negated: u16 },
88    /// A Pauli-diagonal channel: input `s` keeps its key and is scaled by `scale[s]`.
89    Diagonal { scale: Box<[f64]> },
90    /// Everything else, ascending by `(input, output)`.
91    Sparse { entries: Box<[SparseEntry]> },
92}
93
94/// One nonzero PTM entry: `amplitude` takes local Pauli `input` to `output`.
95#[derive(Clone, Copy, Debug)]
96struct SparseEntry {
97    input: u8,
98    output: u8,
99    amplitude: Complex64,
100}
101
102impl PtmChannel {
103    /// A channel on `qubits` (one or two, distinct) from its storage.
104    fn new(qubits: &[usize], storage: PtmStorage, kind: PtmKind) -> Self {
105        debug_assert!(qubits.len() == 1 || qubits.len() == 2);
106        let narrow = |q: usize| u32::try_from(q).expect("PtmChannel: qubit index exceeds u32");
107        let mut padded = [narrow(qubits[0]); 2];
108        for (slot, &q) in padded.iter_mut().zip(qubits) {
109            *slot = narrow(q);
110        }
111        Self {
112            qubits: padded,
113            num_support_qubits: qubits.len() as u8,
114            kind,
115            adjoint: false,
116            storage,
117        }
118    }
119
120    /// The support qubits, in constructor order.
121    pub fn qubits(&self) -> impl ExactSizeIterator<Item = usize> + '_ {
122        self.qubits[..self.num_support_qubits as usize].iter().map(|&q| q as usize)
123    }
124
125    /// The adjoint channel, whose PTM is this one's transpose.
126    pub fn adjoint(&self) -> Self {
127        let storage = match &self.storage {
128            PtmStorage::SignedPermutation { image, negated } => {
129                let mut inverse = [0u8; LOCAL_DIM];
130                let mut inverse_negated = 0u16;
131                for (s, &t) in image[..self.dimension()].iter().enumerate() {
132                    let t = t as usize;
133                    inverse[t] = s as u8;
134                    inverse_negated |= ((negated >> s) & 1) << t;
135                }
136                PtmStorage::SignedPermutation { image: inverse, negated: inverse_negated }
137            }
138            PtmStorage::Diagonal { scale } => PtmStorage::Diagonal { scale: scale.clone() },
139            PtmStorage::Sparse { entries } => {
140                let mut transposed: Vec<((u8, u8), Complex64)> = entries
141                    .iter()
142                    .map(|entry| ((entry.output, entry.input), entry.amplitude))
143                    .collect();
144                transposed.sort_by_key(|&(position, _)| position);
145                sparse_storage(transposed)
146            }
147        };
148        Self { storage, adjoint: !self.adjoint, ..*self }
149    }
150
151    /// `4^k` for `k` support qubits.
152    fn dimension(&self) -> usize {
153        1 << (2 * self.num_support_qubits)
154    }
155
156    /// The local Pauli index of the input's restriction to the support.
157    fn local_index<const W: usize>(&self, x: &[u64; W], z: &[u64; W]) -> usize {
158        self.qubits()
159            .enumerate()
160            .map(|(j, q)| {
161                let (word, bit, _) = qubit_loc(q);
162                read_pauli(x, z, word, bit) << (2 * j)
163            })
164            .fold(0, |index, part| index | part)
165    }
166
167    /// The input with its support overwritten by local Pauli index `t`.
168    fn with_local_index<const W: usize>(
169        &self,
170        x: &[u64; W],
171        z: &[u64; W],
172        t: usize,
173    ) -> ([u64; W], [u64; W]) {
174        let mut new_x = *x;
175        let mut new_z = *z;
176        for (j, q) in self.qubits().enumerate() {
177            let (word, bit, mask) = qubit_loc(q);
178            write_pauli(&mut new_x, &mut new_z, word, bit, mask, (t >> (2 * j)) & 3);
179        }
180        (new_x, new_z)
181    }
182
183    /// Body of `apply` (`transpose == false`) and `apply_adjoint`.
184    fn apply_rows<const W: usize>(
185        &self,
186        transpose: bool,
187        x: &[u64; W],
188        z: &[u64; W],
189        coeff: Complex64,
190        out: &mut OutputBuffer<W>,
191    ) {
192        let s = self.local_index(x, z);
193        match &self.storage {
194            PtmStorage::SignedPermutation { image, negated } => {
195                let (source, t) = if transpose {
196                    let source = image[..self.dimension()]
197                        .iter()
198                        .position(|&image| image as usize == s)
199                        .expect("a signed permutation is a bijection");
200                    (source, source)
201                } else {
202                    (s, image[s] as usize)
203                };
204                let (new_x, new_z) = self.with_local_index(x, z, t);
205                let signed = if (negated >> source) & 1 == 1 { -coeff } else { coeff };
206                out.push(new_x, new_z, signed);
207            }
208            PtmStorage::Diagonal { scale } => out.push(*x, *z, coeff * scale[s]),
209            PtmStorage::Sparse { entries } => {
210                for entry in entries.iter() {
211                    let (from, to) = if transpose {
212                        (entry.output, entry.input)
213                    } else {
214                        (entry.input, entry.output)
215                    };
216                    if from as usize == s {
217                        let (new_x, new_z) = self.with_local_index(x, z, to as usize);
218                        out.push(new_x, new_z, coeff * entry.amplitude);
219                    }
220                }
221            }
222        }
223    }
224
225    /// Calls `visit(s, t, amplitude)` for every stored entry.
226    fn for_each_entry(&self, mut visit: impl FnMut(usize, usize, Complex64)) {
227        match &self.storage {
228            PtmStorage::SignedPermutation { image, negated } => {
229                for (s, &t) in image[..self.dimension()].iter().enumerate() {
230                    let amplitude = if (negated >> s) & 1 == 1 { -ONE } else { ONE };
231                    visit(s, t as usize, amplitude);
232                }
233            }
234            PtmStorage::Diagonal { scale } => {
235                for (s, &value) in scale.iter().enumerate() {
236                    visit(s, s, Complex64::new(value, 0.0));
237                }
238            }
239            PtmStorage::Sparse { entries } => {
240                for entry in entries.iter() {
241                    visit(entry.input as usize, entry.output as usize, entry.amplitude);
242                }
243            }
244        }
245    }
246
247    /// The dense `table[s][t]` (transposed when `adjoint`) in ascending-qubit local indices, the order [`Prepared`] uses.
248    fn transfer_matrix(&self, adjoint: bool) -> [[Complex64; LOCAL_DIM]; LOCAL_DIM] {
249        let reversed = self.num_support_qubits == 2 && self.qubits[0] > self.qubits[1];
250        let reorder = |index: usize| {
251            if reversed {
252                ((index & 3) << 2) | (index >> 2)
253            } else {
254                index
255            }
256        };
257        let mut table = [[ZERO; LOCAL_DIM]; LOCAL_DIM];
258        self.for_each_entry(|s, t, amplitude| {
259            let (row, column) = if adjoint { (t, s) } else { (s, t) };
260            table[reorder(row)][reorder(column)] += amplitude;
261        });
262        table
263    }
264}
265
266/// Sparse storage from `(position, amplitude)` pairs ascending by position, exact zeros dropped.
267fn sparse_storage(entries: Vec<((u8, u8), Complex64)>) -> PtmStorage {
268    PtmStorage::Sparse {
269        entries: entries
270            .into_iter()
271            .filter(|&(_, amplitude)| amplitude != ZERO)
272            .map(|((input, output), amplitude)| SparseEntry { input, output, amplitude })
273            .collect(),
274    }
275}
276
277/// Panics unless the two support qubits differ.
278fn assert_distinct(constructor: &str, q0: usize, q1: usize) {
279    assert!(q0 != q1, "PtmChannel::{constructor}: the two qubits must differ (both are {q0})");
280}
281
282/// Panics unless `value` lies in `[0, 1]`.
283fn assert_probability(constructor: &str, name: &str, value: f64) {
284    assert!(
285        (0.0..=1.0).contains(&value),
286        "PtmChannel::{constructor}: {name} must lie in [0, 1] (got {value})"
287    );
288}
289
290impl<const W: usize> Channel<W> for PtmChannel {
291    fn support(&self) -> [u64; W] {
292        let mut qubits = [0usize; 2];
293        for (slot, q) in qubits.iter_mut().zip(self.qubits()) {
294            *slot = q;
295        }
296        support_mask(&qubits[..self.num_support_qubits as usize])
297    }
298
299    fn debug_name(&self) -> &'static str {
300        self.kind.name(self.adjoint)
301    }
302
303    fn apply(
304        &self,
305        input_x: &[u64; W],
306        input_z: &[u64; W],
307        coeff: Complex64,
308        out: &mut OutputBuffer<W>,
309    ) {
310        self.apply_rows(false, input_x, input_z, coeff, out);
311    }
312
313    fn apply_adjoint(
314        &self,
315        input_x: &[u64; W],
316        input_z: &[u64; W],
317        coeff: Complex64,
318        out: &mut OutputBuffer<W>,
319    ) {
320        self.apply_rows(true, input_x, input_z, coeff, out);
321    }
322
323    /// Tabulated straight from the stored PTM, without probing `apply`; a test holds it equal to the probe's table.
324    fn prepare_for_layer(
325        &self,
326        hash: &Gf2Hash<W>,
327        adjoint: bool,
328        _: PrepareToken,
329    ) -> Option<Prepared<W>> {
330        Some(Prepared::from_transfer_matrix(&self.support(), &self.transfer_matrix(adjoint), hash))
331    }
332}
333
334#[cfg(test)]
335mod tests;
paulistrings/channel/ptm.rs

paulistrings/channel/
ptm.rs
