wai-quantum 0.3.18

//! Learning from quantum **data** — the surviving quantum advantage
//! (`wai.quantum.qdata`).
//!
//! Most "quantum ML" evaporates under dequantization: on *classical* input the
//! quantum kernel is a phasor layer ([[crate::quantum_kernel]]) and the generative
//! model is a phasor sum ([[crate::quantum_phasor]]). What does **not** evaporate is
//! an advantage in *access to quantum data* — and the sharpest provable instance is
//! **Pauli-channel learning** (Chen–Cotler–Huang–Li 2022; Huang et al., *Science*
//! 2022, 40-qubit demo):
//!
//! - To estimate an unknown `n`-qubit Pauli channel's `2^n` Pauli fidelities `{λ_x}`
//!   to accuracy `ε`, **any entanglement-free strategy** — single copies, even
//!   adaptive — needs **Ω(2^n / ε²)** channel uses. It must spread its shots across
//!   `2^n` probes, so each fidelity is seen `M / 2^n` times.
//! - **With a two-copy quantum memory** you Bell-measure the channel's Choi state.
//!   Each measurement returns an error pattern `y ~ p(y)` *directly*, and `M` such
//!   samples estimate **every** fidelity `λ_x = Σ_y p(y)·(−1)^{x·y}` at once, each to
//!   error `~1/√M` — **independent of `n`**. So `ε` costs only **O(1/ε²)** copies.
//!
//! This is an unconditional, exponential separation in *sample* complexity, and it
//! is what makes [[crate::quantum_noise]]'s Bell / Cycle-Benchmarking sampling
//! (`wai.quantum.cal.noise_learn`) exponentially data-efficient — the real moat,
//! surfaced here as a first-class quantum-ML capability. Both estimators are
//! simulated exactly; deterministic given `seed` (reproducible f64).

fn splitmix64(s: &mut u64) -> u64 {
    *s = s.wrapping_add(0x9E37_79B9_7F4A_7C15);
    let mut z = *s;
    z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
    z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
    z ^ (z >> 31)
}
#[inline]
fn u01(s: &mut u64) -> f64 {
    (splitmix64(s) >> 11) as f64 / (1u64 << 53) as f64
}
#[inline]
fn sign_bit(x: u32, y: u32) -> f64 {
    if (x & y).count_ones() & 1 == 0 { 1.0 } else { -1.0 }
}

/// A sparse Pauli channel given by its error distribution `p(y)` over `n`-bit
/// error patterns — the structure `noise_learn` assumes. `weights` sum to 1 and
/// include the identity pattern (`0`). The channel's Pauli fidelities are
/// `λ_x = Σ_y p(y)·(−1)^{x·y}` — its Walsh–Hadamard spectrum.
#[derive(Clone, Debug)]
pub struct SparseChannel {
    pub patterns: Vec<u32>,
    pub weights: Vec<f64>,
    pub n: u32,
}

impl SparseChannel {
    /// The Pauli fidelity `λ_x` (exact, from the model).
    pub fn fidelity(&self, x: u32) -> f64 {
        self.patterns
            .iter()
            .zip(&self.weights)
            .map(|(&y, &w)| w * sign_bit(x, y))
            .sum()
    }
    fn sample(&self, st: &mut u64) -> u32 {
        let u = u01(st);
        let mut c = 0.0;
        for (&y, &w) in self.patterns.iter().zip(&self.weights) {
            c += w;
            if u < c {
                return y;
            }
        }
        *self.patterns.last().unwrap_or(&0)
    }
}

/// A deterministic demonstration channel on `n` qubits: identity with probability
/// `~0.5`, plus a few random weight-≤3 error patterns — sparse and non-trivially
/// entangling in its fidelity spectrum. Weights normalized to sum to 1.
pub fn demo_channel(n: u32, seed: u64) -> SparseChannel {
    let n = n.clamp(1, 16);
    let mask = ((1u64 << n) - 1) as u32;
    let mut st = seed.wrapping_mul(0x2545_F491).wrapping_add(1);
    let mut patterns = vec![0u32]; // identity
    let mut weights = vec![0.0f64];
    let k = 3 + (splitmix64(&mut st) % 3) as usize; // 3..=5 error patterns
    for _ in 0..k {
        // a low-weight error pattern: pick 1..=3 random bit positions
        let bits = 1 + (splitmix64(&mut st) % 3);
        let mut y = 0u32;
        for _ in 0..bits {
            let pos = (splitmix64(&mut st) % n as u64) as u32;
            y |= 1 << pos;
        }
        y &= mask;
        patterns.push(y);
        weights.push(0.02 + 0.10 * u01(&mut st));
    }
    // identity takes the remaining mass (kept ≥ 0.4 so the channel is benign)
    let err_mass: f64 = weights.iter().skip(1).sum();
    weights[0] = (1.0 - err_mass).max(0.4);
    let z: f64 = weights.iter().sum();
    for w in &mut weights {
        *w /= z;
    }
    SparseChannel { patterns, weights, n }
}

/// The **quantum-memory (Bell)** estimator: draw `m` error patterns `y ~ p` (each a
/// two-copy Bell measurement of the Choi state), then estimate *every* fidelity
/// `λ̂_x = (1/m) Σ_i (−1)^{x·y_i}` and return the worst-case error
/// `max_x |λ̂_x − λ_x|`. The error scales `~1/√m`, **independent of `n`**.
pub fn bell_max_error(ch: &SparseChannel, m: u64, seed: u64) -> f64 {
    if m == 0 {
        return 1.0;
    }
    // tally the empirical error distribution (support is small)
    let mut st = seed.wrapping_mul(0x9E37_79B9).wrapping_add(7);
    let mut counts = vec![0u64; ch.patterns.len()];
    let idx: std::collections::HashMap<u32, usize> =
        ch.patterns.iter().enumerate().map(|(i, &y)| (y, i)).collect();
    for _ in 0..m {
        let y = ch.sample(&mut st);
        if let Some(&i) = idx.get(&y) {
            counts[i] += 1;
        }
    }
    let dim = 1u32 << ch.n;
    let mut worst = 0.0f64;
    for x in 0..dim {
        let est: f64 = ch
            .patterns
            .iter()
            .zip(&counts)
            .map(|(&y, &c)| (c as f64 / m as f64) * sign_bit(x, y))
            .sum();
        let e = (est - ch.fidelity(x)).abs();
        if e > worst {
            worst = e;
        }
    }
    worst
}

/// The **entanglement-free (conventional)** estimator: the `m` channel uses must be
/// split across all `2^n` fidelities, so each is measured `m / 2^n` times as `±1`
/// shots of bias `(1+λ_x)/2`. Returns `max_x |λ̂_x − λ_x|`, which scales
/// `~√(2^n / m)` — exponentially worse in `n`. If `m < 2^n` some fidelity is never
/// sampled and the error is maximal (`1.0`), the honest face of the Ω(2^n) bound.
pub fn conventional_max_error(ch: &SparseChannel, m: u64, seed: u64) -> f64 {
    let dim = 1u64 << ch.n;
    let per = m / dim;
    if per == 0 {
        return 1.0;
    }
    let mut worst = 0.0f64;
    for x in 0..dim as u32 {
        let lam = ch.fidelity(x);
        let p_plus = 0.5 * (1.0 + lam);
        let mut st = seed
            .wrapping_mul(0x100_0001)
            .wrapping_add(x as u64 + 1);
        let mut hits = 0u64;
        for _ in 0..per {
            if u01(&mut st) < p_plus {
                hits += 1;
            }
        }
        let est = 2.0 * (hits as f64 / per as f64) - 1.0;
        let e = (est - lam).abs();
        if e > worst {
            worst = e;
        }
    }
    worst
}

/// One point of the sample-complexity curve: at budget `m` channel uses, the
/// worst-case fidelity error of each strategy.
#[derive(Clone, Copy, Debug)]
pub struct SepPoint {
    pub copies: u64,
    pub bell_err: f64,
    pub conventional_err: f64,
}

/// The separation curve over a list of copy budgets `ms` for an `n`-qubit demo
/// channel — the Bell error (n-independent, `~1/√m`) against the entanglement-free
/// error (`~√(2^n/m)`), the exponential gap made visible.
pub fn separation_curve(n: u32, seed: u64, ms: &[u64]) -> Vec<SepPoint> {
    let ch = demo_channel(n, seed);
    ms.iter()
        .map(|&m| SepPoint {
            copies: m,
            bell_err: bell_max_error(&ch, m, seed ^ (m.wrapping_mul(0x51))),
            conventional_err: conventional_max_error(&ch, m, seed ^ (m.wrapping_mul(0x71))),
        })
        .collect()
}

/// Empirically the smallest budget `m` on a geometric grid at which each strategy
/// reaches worst-case error `≤ eps`. Returns `(bell_copies, conventional_copies)`;
/// `conventional / bell ≈ 2^n` — the exponential price of no quantum memory. `0`
/// means the grid's ceiling did not suffice.
pub fn copies_to_epsilon(n: u32, seed: u64, eps: f64, ceil: u64) -> (u64, u64) {
    let ch = demo_channel(n, seed);
    let mut bell = 0u64;
    let mut conv = 0u64;
    let mut m = 16u64;
    while m <= ceil {
        if bell == 0 && bell_max_error(&ch, m, seed ^ (m.wrapping_mul(0x51))) <= eps {
            bell = m;
        }
        if conv == 0 && conventional_max_error(&ch, m, seed ^ (m.wrapping_mul(0x71))) <= eps {
            conv = m;
        }
        if bell != 0 && conv != 0 {
            break;
        }
        m = (m as f64 * 1.5).ceil() as u64;
    }
    (bell, conv)
}

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

    #[test]
    fn identity_channel_has_unit_fidelities() {
        let ch = SparseChannel { patterns: vec![0], weights: vec![1.0], n: 5 };
        for x in 0..(1u32 << 5) {
            assert!((ch.fidelity(x) - 1.0).abs() < 1e-12);
        }
    }

    #[test]
    fn fidelities_are_a_walsh_spectrum() {
        // a single error pattern y flips exactly the fidelities x that anticommute.
        let ch = SparseChannel { patterns: vec![0, 0b101], weights: vec![0.7, 0.3], n: 3 };
        // λ_x = 0.7 + 0.3·(−1)^{x·101}
        assert!((ch.fidelity(0b000) - 1.0).abs() < 1e-12); // 0.7+0.3
        assert!((ch.fidelity(0b100) - 0.4).abs() < 1e-12); // 0.7−0.3
        assert!((ch.fidelity(0b010) - 1.0).abs() < 1e-12); // even overlap
    }

    #[test]
    fn bell_error_shrinks_with_copies() {
        let ch = demo_channel(6, 42);
        let coarse = bell_max_error(&ch, 200, 1);
        let fine = bell_max_error(&ch, 20_000, 1);
        assert!(fine < coarse, "more Bell copies ⇒ smaller error ({fine} !< {coarse})");
        assert!(fine < 0.05, "20k Bell copies resolve all fidelities: {fine}");
    }

    #[test]
    fn the_exponential_separation() {
        // Bell error is ~n-independent; conventional error explodes with n at a
        // fixed copy budget — the whole point.
        let m = 4096u64;
        let bell4 = bell_max_error(&demo_channel(4, 7), m, 3);
        let bell8 = bell_max_error(&demo_channel(8, 7), m, 3);
        let conv4 = conventional_max_error(&demo_channel(4, 7), m, 3);
        let conv8 = conventional_max_error(&demo_channel(8, 7), m, 3);
        // Bell barely moves from 4→8 qubits…
        assert!(bell8 < 3.0 * bell4 + 0.02, "Bell ~n-independent: {bell4} → {bell8}");
        assert!(bell8 < 0.15, "Bell still resolves 8-qubit fidelities at {m} copies: {bell8}");
        // …while conventional collapses (m=4096 < 2^12 spread ⇒ many probes starved).
        assert!(conv8 > conv4, "conventional worsens with n: {conv4} → {conv8}");
        assert!(conv8 > 0.3, "8-qubit conventional is starved at {m} copies: {conv8}");
        assert!(conv8 > 5.0 * bell8, "Bell beats conventional at 8 qubits: {conv8} vs {bell8}");
    }

    #[test]
    fn copies_ratio_is_exponential() {
        let (bell, conv) = copies_to_epsilon(6, 11, 0.08, 5_000_000);
        assert!(bell > 0, "Bell reaches ε");
        assert!(conv > 0, "conventional reaches ε within the ceiling");
        // conventional needs ~2^6 = 64× more copies (loose bounds for sampling noise)
        assert!(conv > 8 * bell, "conventional needs exponentially more: {conv} vs {bell}");
    }

    #[test]
    fn deterministic() {
        let a = separation_curve(6, 99, &[128, 1024, 8192]);
        let b = separation_curve(6, 99, &[128, 1024, 8192]);
        for (p, q) in a.iter().zip(&b) {
            assert_eq!(p.bell_err.to_bits(), q.bell_err.to_bits());
            assert_eq!(p.conventional_err.to_bits(), q.conventional_err.to_bits());
        }
    }
}