wai-quantum 0.3.40

//! Noisy circuits — `wai.quantum.noisy`.
//!
//! Every quantum-advantage experiment to date runs on noisy hardware, and its
//! headline is a fidelity. For random-circuit sampling, the estimate is the
//! linear cross-entropy benchmark,
//!
//! ```text
//! F = 2ⁿ · ⟨p_ideal(x)⟩ − 1
//! ```
//!
//! averaged over the measured bitstrings. Behind it sits a simple error model:
//! a circuit's fidelity is the product of its gates' and readouts' success
//! probabilities. Whether that model holds for a given circuit and noise is a
//! question a classical computer can answer. This module answers it three
//! ways:
//!
//! - **Trajectories.** [`trajectories`] runs the circuit on a pure state many
//!   times. After each gate it draws one branch of each noise channel on the
//!   qubits that gate touched.
//!   - Pauli channels (depolarizing noise, dephasing) pick a Pauli and stay
//!     unitary.
//!   - Amplitude damping jumps with the state's own probability `γ·P(1)` and
//!     is renormalised.
//!
//!   Averages over trajectories converge to the noisy expectation, with a
//!   standard error.
//! - **The exact density matrix.** [`DensityMatrix`] holds `ρ` as a state on
//!   `2n` qubits: rows are qubits `0..n`, columns `n..2n`. A gate `U` acts as
//!   `U ⊗ U*`, and a channel as its Kraus sum. It reuses the state vector's
//!   kernels, and it is the referee for small registers.
//! - **The digital error model.** [`digital_fidelity`] is the product of
//!   success probabilities. [`linear_xeb`] scores samples against the ideal
//!   distribution. [`exact_xeb`] gives the value an infinite sample would
//!   reach.
//!
//! Each trajectory is seeded on its own and results are reduced in trajectory
//! order, so a run is the same bits at any thread count and on every machine.
//!
//! # Checked
//!
//! - **Trajectories against the exact density matrix.** On 5-qubit random
//!   circuits with every noise channel, 4,000 trajectories land within four
//!   standard errors of the exact expectations.
//! - **Sampling.** Sampled bitstrings, with readout flips, follow the exact
//!   distribution, and the sampled XEB converges to [`exact_xeb`].
//!
//! # What the digital error model gets wrong
//!
//! The circuits here are 8-qubit brickwork circuits: random rotations, then CZ
//! on alternating pairs. Under depolarizing noise, the exact XEB (normalised
//! by the noiseless circuit's own) sits *above* the digital model:
//!
//! | layers | two-qubit error 0.005 | 0.01  |
//! |--------|-----------------------|-------|
//! | 4      | 1.035                 | 1.073 |
//! | 8      | 1.045                 | 1.095 |
//! | 16     | 1.034                 | 1.073 |
//! | 32     | 1.036                 | 1.077 |
//! | 48     | 1.041                 | 1.089 |
//!
//! The margin grows with the error rate but not with depth. A per-gate bias
//! would compound with every layer; this one does not. It behaves like a
//! fixed set of errors near the measurement that the model wrongly counts as
//! destroying the signal. A Z error just before a computational-basis
//! measurement, for one, changes no outcome at all.
//!
//! The model therefore understates the fidelity of a shallow-noise circuit by
//! a few percent, independent of size. An experiment that reports fidelity
//! from that model inherits the same offset.
//!
//! # Honest boundaries
//!
//! - **Markovian, gate-local noise.** Errors follow the gate that caused them,
//!   on that gate's qubits. Crosstalk, leakage, coherent over-rotation and
//!   correlated noise in time are not modelled.
//! - **The exact referee is small.** The density matrix is `4ⁿ` amplitudes,
//!   so it stops at 12 qubits.

use crate::gates::{complex, Gate, C};
use crate::quantum_sv::{StateVector, SvError};

/// Per-gate noise, applied after each gate on the qubits it touched, plus
/// readout errors.
#[derive(Clone, Copy, Debug, PartialEq, Default)]
pub struct NoiseModel {
    /// Depolarizing probability after a one-qubit gate (X, Y, Z each `p/3`);
    /// also applied to each qubit of a gate on three or more qubits.
    pub depol1: f64,
    /// Depolarizing probability after a two-qubit gate (each of the 15
    /// non-identity two-qubit Paulis `p/15`).
    pub depol2: f64,
    /// Amplitude-damping strength `γ` per touched qubit per gate.
    pub damping: f64,
    /// Phase-flip probability per touched qubit per gate.
    pub dephasing: f64,
    /// Bit-flip probability of each measured bit.
    pub readout: f64,
}

/// The fidelity the digital error model predicts: the product of `1 − e`
/// over every error location. The locations are:
/// - each gate's Pauli error;
/// - each touched qubit's dephasing, and its damping at the average jump
///   probability `γ/2`;
/// - each readout.
///
/// Each error is taken to destroy the signal entirely; see the module notes
/// for what that misses.
pub fn digital_fidelity(gates: &[Gate], n: u32, noise: &NoiseModel) -> f64 {
    let mut f = 1.0;
    for g in gates {
        match g.qubits.len() {
            1 => f *= 1.0 - noise.depol1,
            2 => f *= 1.0 - noise.depol2,
            k => {
                for _ in 0..k {
                    f *= 1.0 - noise.depol1;
                }
            }
        }
        for _ in &g.qubits {
            f *= (1.0 - noise.dephasing) * (1.0 - noise.damping / 2.0);
        }
    }
    for _ in 0..n {
        f *= 1.0 - noise.readout;
    }
    f
}

// ---------------------------------------------------------------------------
// Randomness
// ---------------------------------------------------------------------------

#[derive(Clone)]
struct Rng(u64);

impl Rng {
    fn next(&mut self) -> u64 {
        self.0 = self.0.wrapping_add(0x9e37_79b9_7f4a_7c15);
        let mut z = self.0;
        z = (z ^ (z >> 30)).wrapping_mul(0xbf58_476d_1ce4_e5b9);
        z = (z ^ (z >> 27)).wrapping_mul(0x94d0_49bb_1331_11eb);
        z ^ (z >> 31)
    }
    /// Uniform on `[0, 1)`.
    fn unit(&mut self) -> f64 {
        (self.next() >> 11) as f64 * (1.0 / 9_007_199_254_740_992.0)
    }
}

fn trajectory_rng(seed: u64, t: usize) -> Rng {
    let mut r = Rng(seed ^ (t as u64).wrapping_mul(0xd1b5_4a32_d192_ed03));
    r.next();
    r
}

const PAULI_I: [C; 4] = [C { re: 1.0, im: 0.0 }, C { re: 0.0, im: 0.0 }, C { re: 0.0, im: 0.0 }, C { re: 1.0, im: 0.0 }];

fn pauli(q: u32, p: u8) -> Gate {
    match p {
        1 => Gate::x(q),
        2 => Gate::y(q),
        3 => Gate::z(q),
        _ => Gate { qubits: vec![q], matrix: PAULI_I.to_vec() },
    }
}

/// The two Kraus operators of amplitude damping.
fn damping_kraus(q: u32, gamma: f64) -> [Gate; 2] {
    let (z, o) = (complex(0.0, 0.0), complex(1.0, 0.0));
    [
        Gate { qubits: vec![q], matrix: vec![o, z, z, complex((1.0 - gamma).sqrt(), 0.0)] },
        Gate { qubits: vec![q], matrix: vec![z, complex(gamma.sqrt(), 0.0), z, z] },
    ]
}

// ---------------------------------------------------------------------------
// Trajectories
// ---------------------------------------------------------------------------

/// Probability that qubit `q` reads 1.
fn prob_one(sv: &StateVector, q: u32) -> f64 {
    let bit = 1usize << (sv.n - 1 - q);
    let mut p = 0.0;
    for (x, (r, i)) in sv.re.iter().zip(&sv.im).enumerate() {
        if x & bit != 0 {
            p += r * r + i * i;
        }
    }
    p
}

fn scale(sv: &mut StateVector, s: f64) {
    sv.re.iter_mut().for_each(|v| *v *= s);
    sv.im.iter_mut().for_each(|v| *v *= s);
}

/// Run one trajectory of `gates` on `n` qubits under `noise`.
fn run_trajectory(gates: &[Gate], n: u32, noise: &NoiseModel, rng: &mut Rng) -> Result<StateVector, SvError> {
    let mut sv = StateVector::zero(n)?;
    for g in gates {
        sv.apply(g)?;
        let k = g.qubits.len();
        if k == 2 && noise.depol2 > 0.0 && rng.unit() < noise.depol2 {
            let which = 1 + (rng.next() % 15) as u8;
            sv.apply(&pauli(g.qubits[0], which >> 2))?;
            sv.apply(&pauli(g.qubits[1], which & 3))?;
        }
        for &q in &g.qubits {
            if k != 2 && noise.depol1 > 0.0 && rng.unit() < noise.depol1 {
                sv.apply(&pauli(q, 1 + (rng.next() % 3) as u8))?;
            }
            if noise.dephasing > 0.0 && rng.unit() < noise.dephasing {
                sv.apply(&Gate::z(q))?;
            }
            if noise.damping > 0.0 {
                let p1 = prob_one(&sv, q);
                let [k0, k1] = damping_kraus(q, noise.damping);
                let jump = noise.damping * p1;
                if rng.unit() < jump {
                    sv.apply(&k1)?;
                    scale(&mut sv, 1.0 / jump.sqrt());
                } else {
                    sv.apply(&k0)?;
                    scale(&mut sv, 1.0 / (1.0 - jump).sqrt());
                }
            }
        }
    }
    Ok(sv)
}

/// A Pauli observable: `Σ c·P`, each string a list of `(qubit, 'X'|'Y'|'Z')`.
pub type Observable = Vec<(f64, Vec<(u32, char)>)>;

/// One trajectory's observable values and measured bitstring.
type Outcome = Result<(Vec<f64>, u64), SvError>;

/// A mean over trajectories and its standard error.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct Estimate {
    pub mean: f64,
    pub stderr: f64,
}

fn estimate(values: &[f64]) -> Estimate {
    let n = values.len() as f64;
    let mean = values.iter().sum::<f64>() / n;
    let var = if values.len() > 1 { values.iter().map(|v| (v - mean) * (v - mean)).sum::<f64>() / (n - 1.0) } else { 0.0 };
    Estimate { mean, stderr: (var / n).sqrt() }
}

/// The outcome of a trajectory run.
#[derive(Clone, Debug, PartialEq)]
pub struct TrajectoryRun {
    /// Each observable's mean over trajectories.
    pub expectations: Vec<Estimate>,
    /// One measured bitstring per trajectory, readout noise included.
    pub samples: Vec<u64>,
}

/// Run `count` noisy trajectories of `gates` on `n` qubits. Each yields the
/// observables' values and one measured bitstring (with readout flips).
/// Trajectory `t` is seeded from `(seed, t)`; results come back in trajectory
/// order whatever `threads` is (wasm32 runs one).
pub fn trajectories(gates: &[Gate], n: u32, noise: &NoiseModel, observables: &[Observable], count: usize, seed: u64, threads: usize) -> Result<TrajectoryRun, SvError> {
    let run = |t: usize| -> Outcome {
        let mut rng = trajectory_rng(seed, t);
        let sv = run_trajectory(gates, n, noise, &mut rng)?;
        let values: Vec<f64> = observables.iter().map(|o| sv.expectation(o)).collect();
        let mut x = sv.sample(1, rng.next())[0];
        if noise.readout > 0.0 {
            for q in 0..n {
                if rng.unit() < noise.readout {
                    x ^= 1 << (n - 1 - q);
                }
            }
        }
        Ok((values, x))
    };
    let threads = if cfg!(target_arch = "wasm32") { 1 } else { threads.max(1).min(count.max(1)) };
    let results: Vec<Outcome> = if threads == 1 {
        (0..count).map(run).collect()
    } else {
        let next = std::sync::atomic::AtomicUsize::new(0);
        let slots: Vec<std::sync::Mutex<Option<Outcome>>> = (0..count).map(|_| std::sync::Mutex::new(None)).collect();
        std::thread::scope(|scope| {
            for _ in 0..threads {
                scope.spawn(|| loop {
                    let t = next.fetch_add(1, std::sync::atomic::Ordering::Relaxed);
                    if t >= count {
                        break;
                    }
                    *slots[t].lock().unwrap() = Some(run(t));
                });
            }
        });
        slots.into_iter().map(|m| m.into_inner().unwrap().expect("every trajectory ran")).collect()
    };
    let mut per_obs: Vec<Vec<f64>> = vec![Vec::with_capacity(count); observables.len()];
    let mut samples = Vec::with_capacity(count);
    for r in results {
        let (values, x) = r?;
        for (k, v) in values.into_iter().enumerate() {
            per_obs[k].push(v);
        }
        samples.push(x);
    }
    Ok(TrajectoryRun { expectations: per_obs.iter().map(|v| estimate(v)).collect(), samples })
}

/// Linear cross-entropy fidelity of `samples` against the ideal state:
/// `2ⁿ · mean(p_ideal(x)) − 1`.
pub fn linear_xeb(ideal: &StateVector, samples: &[u64]) -> Estimate {
    let dim = (1u64 << ideal.n) as f64;
    let v: Vec<f64> = samples
        .iter()
        .map(|&x| {
            let a = ideal.amplitude(x as usize);
            dim * (a.re * a.re + a.im * a.im) - 1.0
        })
        .collect();
    estimate(&v)
}

// ---------------------------------------------------------------------------
// The exact density matrix
// ---------------------------------------------------------------------------

/// Largest register the exact density matrix accepts.
pub const MAX_DENSITY_QUBITS: u32 = 12;

/// `ρ` on `n` qubits, stored as a state on `2n` qubits: entry `(r, c)` at
/// index `r·2ⁿ + c`.
#[derive(Clone, Debug, PartialEq)]
pub struct DensityMatrix {
    pub n: u32,
    vec: StateVector,
}

fn conj(g: &Gate, shift: u32) -> Gate {
    Gate { qubits: g.qubits.iter().map(|&q| q + shift).collect(), matrix: g.matrix.iter().map(|c| complex(c.re, -c.im)).collect() }
}

impl DensityMatrix {
    /// `|0…0⟩⟨0…0|`.
    pub fn zero(n: u32) -> Result<DensityMatrix, SvError> {
        if n > MAX_DENSITY_QUBITS {
            return Err(SvError::TooLarge);
        }
        Ok(DensityMatrix { n, vec: StateVector::zero(2 * n)? })
    }

    /// `ρ → K ρ K†`.
    fn sandwich(&mut self, k: &Gate) -> Result<(), SvError> {
        self.vec.apply(k)?;
        self.vec.apply(&conj(k, self.n))
    }

    /// `ρ → Σ_i w_i K_i ρ K_i†`.
    fn channel(&mut self, terms: &[(f64, Gate)]) -> Result<(), SvError> {
        let mut acc = StateVector { n: self.vec.n, re: vec![0.0; self.vec.re.len()], im: vec![0.0; self.vec.im.len()] };
        for (w, k) in terms {
            let mut t = DensityMatrix { n: self.n, vec: self.vec.clone() };
            t.sandwich(k)?;
            for (a, b) in acc.re.iter_mut().zip(&t.vec.re) {
                *a += w * b;
            }
            for (a, b) in acc.im.iter_mut().zip(&t.vec.im) {
                *a += w * b;
            }
        }
        self.vec = acc;
        Ok(())
    }

    fn depolarize1(&mut self, q: u32, p: f64) -> Result<(), SvError> {
        let mut terms = vec![(1.0 - p, pauli(q, 0))];
        for k in 1..4 {
            terms.push((p / 3.0, pauli(q, k)));
        }
        self.channel(&terms)
    }

    fn depolarize2(&mut self, a: u32, b: u32, p: f64) -> Result<(), SvError> {
        let mut terms = Vec::with_capacity(16);
        for k in 0..16u8 {
            let w = if k == 0 { 1.0 - p } else { p / 15.0 };
            let (pa, pb) = (pauli(a, k >> 2), pauli(b, k & 3));
            terms.push((w, two_qubit(&pa, &pb)));
        }
        self.channel(&terms)
    }

    /// Run `gates` under `noise` from `|0…0⟩⟨0…0|`, exactly.
    pub fn run(gates: &[Gate], n: u32, noise: &NoiseModel) -> Result<DensityMatrix, SvError> {
        let mut rho = DensityMatrix::zero(n)?;
        for g in gates {
            rho.sandwich(g)?;
            let k = g.qubits.len();
            if k == 2 && noise.depol2 > 0.0 {
                rho.depolarize2(g.qubits[0], g.qubits[1], noise.depol2)?;
            }
            for &q in &g.qubits {
                if k != 2 && noise.depol1 > 0.0 {
                    rho.depolarize1(q, noise.depol1)?;
                }
                if noise.dephasing > 0.0 {
                    rho.channel(&[(1.0 - noise.dephasing, pauli(q, 0)), (noise.dephasing, pauli(q, 3))])?;
                }
                if noise.damping > 0.0 {
                    let [k0, k1] = damping_kraus(q, noise.damping);
                    rho.channel(&[(1.0, k0), (1.0, k1)])?;
                }
            }
        }
        Ok(rho)
    }

    /// Entry `ρ_{r,c}`.
    pub fn entry(&self, r: usize, c: usize) -> C {
        self.vec.amplitude((r << self.n) | c)
    }

    /// `Tr ρ`.
    pub fn trace(&self) -> f64 {
        (0..1usize << self.n).map(|x| self.entry(x, x).re).sum()
    }

    /// `Tr(ρ Σ c·P)`.
    pub fn expectation(&self, observable: &[(f64, Vec<(u32, char)>)]) -> f64 {
        let n = self.n;
        let mut total = 0.0;
        for (coef, string) in observable {
            let (mut flip, mut zmask, mut ys) = (0usize, 0usize, 0u32);
            for &(q, p) in string {
                let bit = 1usize << (n - 1 - q);
                match p {
                    'X' => flip |= bit,
                    'Y' => {
                        flip |= bit;
                        zmask |= bit;
                        ys += 1;
                    }
                    'Z' => zmask |= bit,
                    _ => {}
                }
            }
            // Tr(Pρ) = Σ_y ⟨y⊕f| P |y⟩ ρ_{y, y⊕f} = Σ_y i^{ys} (−1)^{|y∧z|} ρ_{y, y⊕f}.
            let mut v = 0.0;
            for y in 0..1usize << n {
                let e = self.entry(y, y ^ flip);
                let sign = if (y & zmask).count_ones().is_multiple_of(2) { 1.0 } else { -1.0 };
                v += sign
                    * match ys % 4 {
                        0 => e.re,
                        1 => -e.im,
                        2 => -e.re,
                        _ => e.im,
                    };
            }
            total += coef * v;
        }
        total
    }

    /// The measured distribution, readout flips included.
    pub fn probabilities(&self, readout: f64) -> Vec<f64> {
        let n = self.n;
        let mut p: Vec<f64> = (0..1usize << n).map(|x| self.entry(x, x).re).collect();
        if readout > 0.0 {
            for q in 0..n {
                let bit = 1usize << (n - 1 - q);
                for x in 0..p.len() {
                    if x & bit == 0 {
                        let (a, b) = (p[x], p[x | bit]);
                        p[x] = (1.0 - readout) * a + readout * b;
                        p[x | bit] = readout * a + (1.0 - readout) * b;
                    }
                }
            }
        }
        p
    }
}

/// The tensor product `a ⊗ b` of two one-qubit gates on different qubits.
fn two_qubit(a: &Gate, b: &Gate) -> Gate {
    let mut m = vec![complex(0.0, 0.0); 16];
    for r in 0..4 {
        for c in 0..4 {
            let (x, y) = (a.matrix[(r >> 1) * 2 + (c >> 1)], b.matrix[(r & 1) * 2 + (c & 1)]);
            m[r * 4 + c] = complex(x.re * y.re - x.im * y.im, x.re * y.im + x.im * y.re);
        }
    }
    Gate { qubits: vec![a.qubits[0], b.qubits[0]], matrix: m }
}

/// The noiseless circuit's own cross-entropy, `2ⁿ Σ_x p_ideal(x)² − 1`: 1 for
/// a fully scrambling circuit, larger for one that is not. Noisy values are
/// compared with the digital model after dividing by it.
pub fn ideal_xeb(ideal: &StateVector) -> f64 {
    let dim = (1u64 << ideal.n) as f64;
    let mut s = 0.0;
    for (r, i) in ideal.re.iter().zip(&ideal.im) {
        let q = r * r + i * i;
        s += q * q;
    }
    dim * s - 1.0
}

/// The cross-entropy fidelity an infinite sample would give:
/// `2ⁿ Σ_x p_noisy(x) p_ideal(x) − 1`.
pub fn exact_xeb(ideal: &StateVector, noisy: &DensityMatrix, readout: f64) -> f64 {
    let dim = (1u64 << ideal.n) as f64;
    let p = noisy.probabilities(readout);
    let mut s = 0.0;
    for (x, px) in p.iter().enumerate() {
        let a = ideal.amplitude(x);
        s += px * (a.re * a.re + a.im * a.im);
    }
    dim * s - 1.0
}

#[cfg(test)]
mod tests {
    use super::*;

    /// A brickwork random circuit: random one-qubit rotations, then CZ on
    /// alternating pairs, `layers` times.
    fn random_circuit(n: u32, layers: usize, seed: u64) -> Vec<Gate> {
        let mut rng = Rng(seed);
        let mut g = Vec::new();
        for l in 0..layers {
            for q in 0..n {
                let a = core::f64::consts::TAU * rng.unit();
                g.push(match rng.next() % 3 {
                    0 => Gate::rx(q, a),
                    1 => Gate::ry(q, a),
                    _ => Gate::sqrt_w(q),
                });
            }
            let mut q = (l % 2) as u32;
            while q + 1 < n {
                g.push(Gate::cz(q, q + 1));
                q += 2;
            }
        }
        g
    }

    fn observables(n: u32) -> Vec<Observable> {
        vec![vec![(1.0, vec![(0, 'Z')])], vec![(1.0, vec![(1, 'X'), (2, 'X')])], vec![(1.0, vec![(n - 1, 'Y'), (0, 'Z')])], vec![(0.5, vec![(2, 'Z')]), (0.5, vec![(3, 'Y')])]]
    }

    #[test]
    fn noiseless_runs_match_the_pure_state() {
        let g = random_circuit(5, 4, 1);
        let mut sv = StateVector::zero(5).unwrap();
        sv.run(&g, 0).unwrap();
        let rho = DensityMatrix::run(&g, 5, &NoiseModel::default()).unwrap();
        let tr = trajectories(&g, 5, &NoiseModel::default(), &observables(5), 3, 2, 2).unwrap();
        for (k, o) in observables(5).iter().enumerate() {
            let pure = sv.expectation(o);
            assert!((rho.expectation(o) - pure).abs() < 1e-12);
            assert!((tr.expectations[k].mean - pure).abs() < 1e-12);
            assert!(tr.expectations[k].stderr < 1e-12);
        }
        assert!((rho.trace() - 1.0).abs() < 1e-12);
    }

    #[test]
    fn trajectories_converge_to_the_density_matrix() {
        let n = 5;
        let g = random_circuit(n, 5, 7);
        let noise = NoiseModel { depol1: 0.01, depol2: 0.03, damping: 0.02, dephasing: 0.01, readout: 0.0 };
        let rho = DensityMatrix::run(&g, n, &noise).unwrap();
        assert!((rho.trace() - 1.0).abs() < 1e-12);
        let tr = trajectories(&g, n, &noise, &observables(n), 4000, 11, 4).unwrap();
        for (k, o) in observables(n).iter().enumerate() {
            let exact = rho.expectation(o);
            let e = tr.expectations[k];
            assert!((e.mean - exact).abs() < 4.0 * e.stderr + 1e-9, "observable {k}: {} ± {} vs {exact}", e.mean, e.stderr);
        }
    }

    #[test]
    fn readout_noise_and_sampling_match_the_exact_distribution() {
        let n = 4;
        let g = random_circuit(n, 3, 3);
        let noise = NoiseModel { depol2: 0.02, readout: 0.05, ..NoiseModel::default() };
        let p = DensityMatrix::run(&g, n, &noise).unwrap().probabilities(noise.readout);
        assert!((p.iter().sum::<f64>() - 1.0).abs() < 1e-12);
        let shots = 20_000;
        let tr = trajectories(&g, n, &noise, &[], shots, 5, 4).unwrap();
        let mut counts = vec![0usize; 1 << n];
        for &x in &tr.samples {
            counts[x as usize] += 1;
        }
        for (x, &c) in counts.iter().enumerate() {
            let expected = p[x] * shots as f64;
            let sd = (expected * (1.0 - p[x])).sqrt().max(1.0);
            assert!((c as f64 - expected).abs() < 5.0 * sd, "x={x}: {c} vs {expected:.1}");
        }
    }

    #[test]
    fn runs_do_not_depend_on_threads() {
        let g = random_circuit(6, 4, 9);
        let noise = NoiseModel { depol1: 0.02, depol2: 0.05, damping: 0.01, dephasing: 0.01, readout: 0.02 };
        let a = trajectories(&g, 6, &noise, &observables(6), 50, 3, 1).unwrap();
        let b = trajectories(&g, 6, &noise, &observables(6), 50, 3, 5).unwrap();
        assert_eq!(a, b);
    }

    fn self_xeb(ideal: &StateVector) -> f64 {
        ideal_xeb(ideal)
    }

    /// The exact XEB, normalised by the noiseless circuit's own, over the
    /// digital model's prediction.
    fn model_ratio(n: u32, layers: usize, p2: f64) -> f64 {
        let g = random_circuit(n, layers, 21);
        let mut ideal = StateVector::zero(n).unwrap();
        ideal.run(&g, 0).unwrap();
        let noise = NoiseModel { depol1: p2 / 10.0, depol2: p2, ..NoiseModel::default() };
        let exact = exact_xeb(&ideal, &DensityMatrix::run(&g, n, &noise).unwrap(), 0.0);
        exact / self_xeb(&ideal) / digital_fidelity(&g, n, &noise)
    }

    #[test]
    fn the_digital_model_is_low_by_a_depth_independent_margin() {
        // Under Pauli noise the exact (normalised) XEB sits above the digital
        // model by a margin that grows with the error rate but not with depth:
        // a fixed number of errors near the end, not a per-gate bias.
        for p2 in [0.005, 0.01] {
            let (shallow, deep) = (model_ratio(7, 8, p2), model_ratio(7, 24, p2));
            eprintln!("p2={p2}: {shallow:.4} {deep:.4}");
            assert!(shallow > 1.0 && deep > 1.0, "p2={p2}: {shallow} {deep}");
            assert!((shallow - deep).abs() < 0.03, "p2={p2}: {shallow} {deep}");
            assert!(deep < 1.0 + 12.0 * p2, "p2={p2}: {deep}");
        }
        // Noiseless: the exact XEB is the ideal circuit's own, 2ⁿ Σp² − 1.
        let g = random_circuit(7, 8, 21);
        let mut ideal = StateVector::zero(7).unwrap();
        ideal.run(&g, 0).unwrap();
        let noiseless = DensityMatrix::run(&g, 7, &NoiseModel::default()).unwrap();
        assert!((exact_xeb(&ideal, &noiseless, 0.0) - self_xeb(&ideal)).abs() < 1e-9);
    }

    #[test]
    fn sampled_xeb_converges_to_the_exact_value() {
        let n = 7;
        let g = random_circuit(n, 7, 4);
        let mut ideal = StateVector::zero(n).unwrap();
        ideal.run(&g, 0).unwrap();
        let noise = NoiseModel { depol1: 0.002, depol2: 0.015, readout: 0.01, ..NoiseModel::default() };
        let exact = exact_xeb(&ideal, &DensityMatrix::run(&g, n, &noise).unwrap(), noise.readout);
        let tr = trajectories(&g, n, &noise, &[], 20_000, 8, 4).unwrap();
        let est = linear_xeb(&ideal, &tr.samples);
        assert!((est.mean - exact).abs() < 4.0 * est.stderr, "{} ± {} vs {exact}", est.mean, est.stderr);
    }
}