paulistrings/channel/
noise.rs

1//! Noise channels: Depolarizing, Dephasing, AmplitudeDamping. See §6.
2
3#![allow(unused)]
4
5use super::{Channel, OutputBuffer};
6use num_complex::Complex64;
7
8/// Single-qubit depolarizing noise with error probability `p`.
9///
10/// In the Heisenberg picture this is just a coefficient rescaling: the
11/// identity on the support qubit is preserved unchanged; every non-identity
12/// Pauli on the support is multiplied by `1 - 4p/3`. Self-adjoint, so the
13/// default `apply_adjoint` from the trait is correct.
14///
15/// # Examples
16///
17/// ```
18/// use paulistrings::channel::Depolarizing;
19/// let ch = Depolarizing { support: [3], p: 0.05 };
20/// # let _ = ch;
21/// ```
22pub struct Depolarizing {
23    /// The single qubit this channel acts on.
24    pub support: [u32; 1],
25    /// Error probability `p ∈ [0, 1]`.
26    pub p: f64,
27}
28
29impl<const W: usize> Channel<W> for Depolarizing {
30    #[inline]
31    fn max_fanout(&self) -> usize {
32        1
33    }
34
35    #[inline]
36    fn support(&self) -> &[u32] {
37        &self.support
38    }
39
40    #[inline]
41    fn apply(
42        &self,
43        input_x: &[u64; W],
44        input_z: &[u64; W],
45        coeff: Complex64,
46        out: &mut OutputBuffer<'_, W>,
47    ) {
48        let q = self.support[0] as usize;
49        debug_assert!(q < 64 * W);
50        let word = q / 64;
51        let bit = q % 64;
52        let x_bit = (input_x[word] >> bit) & 1;
53        let z_bit = (input_z[word] >> bit) & 1;
54        let scale = if (x_bit | z_bit) == 0 {
55            1.0
56        } else {
57            1.0 - 4.0 * self.p / 3.0
58        };
59        out.push(*input_x, *input_z, coeff * scale);
60    }
61}
62
63/// Single-qubit dephasing noise with error probability `p`.
64///
65/// Heisenberg dual: `E*(P) = (1-p) P + p Z P Z`. So I and Z are preserved;
66/// X and Y are scaled by `1 - 2p` (`Z·X·Z = -X`, `Z·Y·Z = -Y`). Equivalently,
67/// the scale fires iff the support qubit's `x_bit` is set. Self-adjoint.
68pub struct Dephasing {
69    /// The single qubit this channel acts on.
70    pub support: [u32; 1],
71    /// Error probability `p ∈ [0, 1]`.
72    pub p: f64,
73}
74
75impl<const W: usize> Channel<W> for Dephasing {
76    #[inline]
77    fn max_fanout(&self) -> usize {
78        1
79    }
80
81    #[inline]
82    fn support(&self) -> &[u32] {
83        &self.support
84    }
85
86    #[inline]
87    fn apply(
88        &self,
89        input_x: &[u64; W],
90        input_z: &[u64; W],
91        coeff: Complex64,
92        out: &mut OutputBuffer<'_, W>,
93    ) {
94        let q = self.support[0] as usize;
95        debug_assert!(q < 64 * W);
96        let word = q / 64;
97        let bit = q % 64;
98        let x_bit = (input_x[word] >> bit) & 1;
99        let scale = if x_bit == 1 { 1.0 - 2.0 * self.p } else { 1.0 };
100        out.push(*input_x, *input_z, coeff * scale);
101    }
102}
103
104/// Single-qubit amplitude damping with parameter `gamma`.
105///
106/// The only noise in the built-in set with genuine fan-out > 1. Heisenberg
107/// dual via Kraus operators `E_0 = |0⟩⟨0| + √(1-γ)|1⟩⟨1|` and
108/// `E_1 = √γ |0⟩⟨1|`:
109///
110/// - `I → I`
111/// - `X → √(1-γ) X`
112/// - `Y → √(1-γ) Y`
113/// - `Z → (1-γ) Z + γ I`   (the only fanout-2 case)
114///
115/// Not self-adjoint in general — the Heisenberg adjoint of amplitude damping
116/// is not amplitude damping itself. v0.1 uses the default `apply_adjoint =
117/// apply` placeholder; users propagating in `Direction::Heisenberg` mode
118/// with this channel should be aware of that approximation. (The
119/// non-Heisenberg use-case is the common one and is exact.)
120pub struct AmplitudeDamping {
121    /// The single qubit this channel acts on.
122    pub support: [u32; 1],
123    /// Damping parameter `γ ∈ [0, 1]`. The amplitude that escapes to the
124    /// `|0⟩` state per application.
125    pub gamma: f64,
126}
127
128impl<const W: usize> Channel<W> for AmplitudeDamping {
129    #[inline]
130    fn max_fanout(&self) -> usize {
131        2
132    }
133
134    #[inline]
135    fn support(&self) -> &[u32] {
136        &self.support
137    }
138
139    fn apply(
140        &self,
141        input_x: &[u64; W],
142        input_z: &[u64; W],
143        coeff: Complex64,
144        out: &mut OutputBuffer<'_, W>,
145    ) {
146        let q = self.support[0] as usize;
147        debug_assert!(q < 64 * W);
148        let word = q / 64;
149        let bit = q % 64;
150        let mask = 1u64 << bit;
151        let x_bit = (input_x[word] >> bit) & 1;
152        let z_bit = (input_z[word] >> bit) & 1;
153        match (x_bit, z_bit) {
154            (0, 0) => {
155                // I → I.
156                out.push(*input_x, *input_z, coeff);
157            }
158            (1, _) => {
159                // X or Y → √(1-γ) · same.
160                let scale = (1.0 - self.gamma).sqrt();
161                out.push(*input_x, *input_z, coeff * scale);
162            }
163            (0, 1) => {
164                // Z → (1-γ) Z + γ I. Emit Z first (matches the order in the
165                // doc-comment), then I (with the support's z-bit cleared).
166                out.push(*input_x, *input_z, coeff * (1.0 - self.gamma));
167                let mut nz = *input_z;
168                nz[word] &= !mask;
169                out.push(*input_x, nz, coeff * self.gamma);
170            }
171            _ => unreachable!(),
172        }
173    }
174}
175
176#[cfg(test)]
177mod tests {
178    use super::*;
179    use crate::pauli_string::PauliString;
180
181    const TOL: f64 = 1e-12;
182
183    #[allow(clippy::type_complexity)]
184    fn make_buf<const W: usize>(
185        cap: usize,
186    ) -> (Vec<[u64; W]>, Vec<[u64; W]>, Vec<Complex64>, usize) {
187        (
188            vec![[0u64; W]; cap],
189            vec![[0u64; W]; cap],
190            vec![Complex64::new(0.0, 0.0); cap],
191            0,
192        )
193    }
194
195    fn approx_eq(a: Complex64, b: Complex64, tol: f64) -> bool {
196        (a - b).norm() <= tol
197    }
198
199    /// Identity on the support qubit is preserved exactly — coefficient
200    /// rescaling does not touch the I sector.
201    #[test]
202    fn depolarizing_passes_identity_through() {
203        let ch = Depolarizing {
204            support: [0],
205            p: 0.1,
206        };
207        let p = PauliString::<1>::identity();
208        let (mut bx, mut bz, mut bc, mut len) = make_buf::<1>(1);
209        let mut buf = OutputBuffer::<1> {
210            x: &mut bx,
211            z: &mut bz,
212            coeff: &mut bc,
213            len: &mut len,
214        };
215        <Depolarizing as Channel<1>>::apply(&ch, &p.x, &p.z, Complex64::new(2.0, 0.0), &mut buf);
216        assert_eq!(len, 1);
217        assert_eq!(bx[0], p.x);
218        assert_eq!(bz[0], p.z);
219        assert!(approx_eq(bc[0], Complex64::new(2.0, 0.0), TOL));
220    }
221
222    /// X, Y, Z on the support qubit each get scaled by `1 - 4p/3`.
223    #[test]
224    fn depolarizing_scales_xyz() {
225        let p = 0.15;
226        let ch = Depolarizing { support: [0], p };
227        let scale = 1.0 - 4.0 * p / 3.0;
228        for pauli in [
229            PauliString::<1>::x(0),
230            PauliString::<1>::y(0),
231            PauliString::<1>::z(0),
232        ] {
233            let (mut bx, mut bz, mut bc, mut len) = make_buf::<1>(1);
234            let mut buf = OutputBuffer::<1> {
235                x: &mut bx,
236                z: &mut bz,
237                coeff: &mut bc,
238                len: &mut len,
239            };
240            <Depolarizing as Channel<1>>::apply(
241                &ch,
242                &pauli.x,
243                &pauli.z,
244                Complex64::new(1.0, 0.0),
245                &mut buf,
246            );
247            assert_eq!(len, 1);
248            assert_eq!(bx[0], pauli.x);
249            assert_eq!(bz[0], pauli.z);
250            assert!(approx_eq(bc[0], Complex64::new(scale, 0.0), TOL));
251        }
252    }
253
254    /// Off-support qubits are ignored: a Z on qubit 1 with the channel
255    /// supported on qubit 0 leaves the coefficient untouched (the support
256    /// qubit is in I-state).
257    #[test]
258    fn depolarizing_off_support_is_no_op() {
259        let ch = Depolarizing {
260            support: [0],
261            p: 0.2,
262        };
263        let p = PauliString::<1>::z(1);
264        let (mut bx, mut bz, mut bc, mut len) = make_buf::<1>(1);
265        let mut buf = OutputBuffer::<1> {
266            x: &mut bx,
267            z: &mut bz,
268            coeff: &mut bc,
269            len: &mut len,
270        };
271        <Depolarizing as Channel<1>>::apply(&ch, &p.x, &p.z, Complex64::new(3.0, 0.0), &mut buf);
272        assert_eq!(len, 1);
273        assert!(approx_eq(bc[0], Complex64::new(3.0, 0.0), TOL));
274    }
275
276    /// W=2: support qubit lives in word 1 (qubit 64+). The scale must
277    /// trigger when bits at qubit 64 are non-identity.
278    #[test]
279    fn depolarizing_w2_word_boundary() {
280        let p = 0.25;
281        let ch = Depolarizing { support: [64], p };
282        let scale = 1.0 - 4.0 * p / 3.0;
283        // X on qubit 64 → word-1 x-bit set.
284        let pauli = PauliString::<2>::x(64);
285        let (mut bx, mut bz, mut bc, mut len) = make_buf::<2>(1);
286        let mut buf = OutputBuffer::<2> {
287            x: &mut bx,
288            z: &mut bz,
289            coeff: &mut bc,
290            len: &mut len,
291        };
292        <Depolarizing as Channel<2>>::apply(
293            &ch,
294            &pauli.x,
295            &pauli.z,
296            Complex64::new(1.0, 0.0),
297            &mut buf,
298        );
299        assert_eq!(len, 1);
300        assert_eq!(bx[0], pauli.x);
301        assert_eq!(bz[0], pauli.z);
302        assert!(approx_eq(bc[0], Complex64::new(scale, 0.0), TOL));
303    }
304
305    /// I and Z commute with Z, so dephasing leaves their coefficients alone.
306    #[test]
307    fn dephasing_preserves_i_and_z() {
308        let ch = Dephasing {
309            support: [0],
310            p: 0.3,
311        };
312        for pauli in [PauliString::<1>::identity(), PauliString::<1>::z(0)] {
313            let (mut bx, mut bz, mut bc, mut len) = make_buf::<1>(1);
314            let mut buf = OutputBuffer::<1> {
315                x: &mut bx,
316                z: &mut bz,
317                coeff: &mut bc,
318                len: &mut len,
319            };
320            <Dephasing as Channel<1>>::apply(
321                &ch,
322                &pauli.x,
323                &pauli.z,
324                Complex64::new(2.5, 0.0),
325                &mut buf,
326            );
327            assert_eq!(len, 1);
328            assert_eq!(bx[0], pauli.x);
329            assert_eq!(bz[0], pauli.z);
330            assert!(approx_eq(bc[0], Complex64::new(2.5, 0.0), TOL));
331        }
332    }
333
334    /// X and Y both anticommute with Z, so dephasing scales them by 1 - 2p.
335    #[test]
336    fn dephasing_scales_x_and_y() {
337        let p = 0.2;
338        let ch = Dephasing { support: [0], p };
339        let scale = 1.0 - 2.0 * p;
340        for pauli in [PauliString::<1>::x(0), PauliString::<1>::y(0)] {
341            let (mut bx, mut bz, mut bc, mut len) = make_buf::<1>(1);
342            let mut buf = OutputBuffer::<1> {
343                x: &mut bx,
344                z: &mut bz,
345                coeff: &mut bc,
346                len: &mut len,
347            };
348            <Dephasing as Channel<1>>::apply(
349                &ch,
350                &pauli.x,
351                &pauli.z,
352                Complex64::new(1.0, 0.0),
353                &mut buf,
354            );
355            assert_eq!(len, 1);
356            assert_eq!(bx[0], pauli.x);
357            assert_eq!(bz[0], pauli.z);
358            assert!(approx_eq(bc[0], Complex64::new(scale, 0.0), TOL));
359        }
360    }
361
362    /// W=2: dephasing on qubit 64 scales an X@64 by 1 - 2p.
363    #[test]
364    fn dephasing_w2_word_boundary() {
365        let p = 0.4;
366        let ch = Dephasing { support: [64], p };
367        let scale = 1.0 - 2.0 * p;
368        let pauli = PauliString::<2>::x(64);
369        let (mut bx, mut bz, mut bc, mut len) = make_buf::<2>(1);
370        let mut buf = OutputBuffer::<2> {
371            x: &mut bx,
372            z: &mut bz,
373            coeff: &mut bc,
374            len: &mut len,
375        };
376        <Dephasing as Channel<2>>::apply(
377            &ch,
378            &pauli.x,
379            &pauli.z,
380            Complex64::new(1.0, 0.0),
381            &mut buf,
382        );
383        assert_eq!(len, 1);
384        assert!(approx_eq(bc[0], Complex64::new(scale, 0.0), TOL));
385    }
386
387    /// I on the support is fixed by amplitude damping (fanout 1, coeff
388    /// preserved).
389    #[test]
390    fn amplitude_damping_passes_identity_through() {
391        let ch = AmplitudeDamping {
392            support: [0],
393            gamma: 0.3,
394        };
395        let p = PauliString::<1>::identity();
396        let (mut bx, mut bz, mut bc, mut len) = make_buf::<1>(2);
397        let mut buf = OutputBuffer::<1> {
398            x: &mut bx,
399            z: &mut bz,
400            coeff: &mut bc,
401            len: &mut len,
402        };
403        <AmplitudeDamping as Channel<1>>::apply(
404            &ch,
405            &p.x,
406            &p.z,
407            Complex64::new(2.0, 0.0),
408            &mut buf,
409        );
410        assert_eq!(len, 1);
411        assert_eq!(bx[0], p.x);
412        assert_eq!(bz[0], p.z);
413        assert!(approx_eq(bc[0], Complex64::new(2.0, 0.0), TOL));
414    }
415
416    /// X and Y on the support each get scaled by √(1-γ), fanout 1.
417    #[test]
418    fn amplitude_damping_scales_x_and_y_by_sqrt() {
419        let gamma = 0.2;
420        let ch = AmplitudeDamping {
421            support: [0],
422            gamma,
423        };
424        let scale = (1.0 - gamma).sqrt();
425        for pauli in [PauliString::<1>::x(0), PauliString::<1>::y(0)] {
426            let (mut bx, mut bz, mut bc, mut len) = make_buf::<1>(2);
427            let mut buf = OutputBuffer::<1> {
428                x: &mut bx,
429                z: &mut bz,
430                coeff: &mut bc,
431                len: &mut len,
432            };
433            <AmplitudeDamping as Channel<1>>::apply(
434                &ch,
435                &pauli.x,
436                &pauli.z,
437                Complex64::new(1.0, 0.0),
438                &mut buf,
439            );
440            assert_eq!(len, 1);
441            assert_eq!(bx[0], pauli.x);
442            assert_eq!(bz[0], pauli.z);
443            assert!(approx_eq(bc[0], Complex64::new(scale, 0.0), TOL));
444        }
445    }
446
447    /// Z on the support fans out to `(1-γ)·Z + γ·I`. The first emit is the
448    /// Z term, the second is the I term (z-bit cleared on the support
449    /// qubit).
450    #[test]
451    fn amplitude_damping_z_fans_out_to_z_plus_i() {
452        let gamma = 0.25;
453        let ch = AmplitudeDamping {
454            support: [0],
455            gamma,
456        };
457        let p = PauliString::<1>::z(0);
458        let (mut bx, mut bz, mut bc, mut len) = make_buf::<1>(2);
459        let mut buf = OutputBuffer::<1> {
460            x: &mut bx,
461            z: &mut bz,
462            coeff: &mut bc,
463            len: &mut len,
464        };
465        <AmplitudeDamping as Channel<1>>::apply(
466            &ch,
467            &p.x,
468            &p.z,
469            Complex64::new(1.0, 0.0),
470            &mut buf,
471        );
472        assert_eq!(len, 2);
473        // First: (1-γ)·Z (z-bit kept).
474        assert_eq!(bx[0], p.x);
475        assert_eq!(bz[0], p.z);
476        assert!(approx_eq(bc[0], Complex64::new(1.0 - gamma, 0.0), TOL));
477        // Second: γ·I (z-bit cleared).
478        let id = PauliString::<1>::identity();
479        assert_eq!(bx[1], id.x);
480        assert_eq!(bz[1], id.z);
481        assert!(approx_eq(bc[1], Complex64::new(gamma, 0.0), TOL));
482    }
483
484    /// W=2: Z on qubit 64 fans out to (1-γ)·Z@64 + γ·I, with the I term's
485    /// z-bit cleared in word 1 only.
486    #[test]
487    fn amplitude_damping_w2_word_boundary() {
488        let gamma = 0.4;
489        let ch = AmplitudeDamping {
490            support: [64],
491            gamma,
492        };
493        let p = PauliString::<2>::z(64);
494        let (mut bx, mut bz, mut bc, mut len) = make_buf::<2>(2);
495        let mut buf = OutputBuffer::<2> {
496            x: &mut bx,
497            z: &mut bz,
498            coeff: &mut bc,
499            len: &mut len,
500        };
501        <AmplitudeDamping as Channel<2>>::apply(
502            &ch,
503            &p.x,
504            &p.z,
505            Complex64::new(1.0, 0.0),
506            &mut buf,
507        );
508        assert_eq!(len, 2);
509        assert_eq!(bx[0], p.x);
510        assert_eq!(bz[0], p.z);
511        assert!(approx_eq(bc[0], Complex64::new(1.0 - gamma, 0.0), TOL));
512        // I term: word 1 z-bit cleared.
513        assert_eq!(bx[1], [0u64, 0u64]);
514        assert_eq!(bz[1], [0u64, 0u64]);
515        assert!(approx_eq(bc[1], Complex64::new(gamma, 0.0), TOL));
516    }
517}
paulistrings/channel/noise.rs

paulistrings/channel/
noise.rs
