wai-quantum 0.3.19

//! The full Born machine on the qFHRR phasor substrate (`wai.quantum.qfhrr`).
//!
//! The phasor bridge ([`crate::quantum_phasor`]) showed a *scalar* quantum
//! feature-map function equals a classical phasor sum. This runs the whole
//! **generative distribution**. A Born machine's amplitudes are complex numbers
//! `a = m·e^{iθ}` — and the phase `θ` is the interference-carrying quantity (`|a|²`
//! is what you sample; the phases are what make wrong answers cancel when
//! amplitudes add). So evolve the state with the phases held at **qFHRR
//! resolution** — quantized to `b` bits, `2^b` phase levels, integer-native — and
//! read `P(x) = |a_x|²`.
//!
//! The result: a *trained* Quantum Circuit Born Machine deploys end-to-end on the
//! phasor substrate. As the phase precision grows the qFHRR distribution converges
//! to the exact state vector (total-variation distance → 0), and even a few bits
//! reproduce the generative model well — the deployable form of interference-based
//! generation. Honest boundary: the phases are quantized (the qFHRR contract);
//! amplitude magnitudes are kept in `f64` here, so a fully-integer substrate would
//! quantize those too (a further step). Deterministic f64.

use std::f64::consts::PI;

#[derive(Clone, Copy)]
struct C {
    re: f64,
    im: f64,
}
impl C {
    #[inline]
    fn mag2(self) -> f64 {
        self.re * self.re + self.im * self.im
    }
}

/// Project an amplitude onto the qFHRR grid: keep magnitude, snap the phase to one
/// of `2^bits` levels. `bits == 0` means no quantization (exact).
#[inline]
fn qfhrr(a: C, bits: u32) -> C {
    if bits == 0 {
        return a;
    }
    let m = a.mag2().sqrt();
    if m <= 0.0 {
        return C { re: 0.0, im: 0.0 };
    }
    let theta = a.im.atan2(a.re);
    let levels = (1u64 << bits) as f64;
    let tq = (theta / (2.0 * PI) * levels).round() / levels * 2.0 * PI;
    C { re: m * tq.cos(), im: m * tq.sin() }
}

fn apply_ry(amps: &mut [C], q: u8, theta: f64, bits: u32) {
    let (c, s) = ((theta * 0.5).cos(), (theta * 0.5).sin());
    let bit = 1usize << q;
    let mut i = 0;
    while i < amps.len() {
        if i & bit == 0 {
            let j = i | bit;
            let (a, b) = (amps[i], amps[j]);
            amps[i] = qfhrr(C { re: a.re * c - b.re * s, im: a.im * c - b.im * s }, bits);
            amps[j] = qfhrr(C { re: a.re * s + b.re * c, im: a.im * s + b.im * c }, bits);
        }
        i += 1;
    }
}
fn apply_rz(amps: &mut [C], q: u8, theta: f64, bits: u32) {
    let (cm, sm) = ((theta * 0.5).cos(), (theta * 0.5).sin());
    let bit = 1usize << q;
    for (i, a) in amps.iter_mut().enumerate() {
        let (er, ei) = if i & bit == 0 { (cm, -sm) } else { (cm, sm) };
        *a = qfhrr(C { re: a.re * er - a.im * ei, im: a.re * ei + a.im * er }, bits);
    }
}
fn apply_cx(amps: &mut [C], c: u8, t: u8) {
    let (cb, tb) = (1usize << c, 1usize << t);
    for i in 0..amps.len() {
        if (i & cb) != 0 && (i & tb) == 0 {
            amps.swap(i, i | tb);
        }
    }
}

/// Run the hardware-efficient ansatz (`layers` blocks of RY/RZ + a CX ring) with
/// amplitude phases held at `bits`-bit qFHRR resolution, returning the (normalized)
/// Born distribution `P(x)`. `params.len() == layers·2·n`. `bits == 0` is exact.
pub fn born_probs(n: u8, layers: u8, params: &[f64], bits: u32) -> Vec<f64> {
    let dim = 1usize << n;
    let mut amps = vec![C { re: 0.0, im: 0.0 }; dim];
    amps[0] = C { re: 1.0, im: 0.0 };
    let mut p = 0;
    for _ in 0..layers {
        for q in 0..n {
            apply_ry(&mut amps, q, params[p], bits);
            p += 1;
            apply_rz(&mut amps, q, params[p], bits);
            p += 1;
        }
        for q in 0..n {
            apply_cx(&mut amps, q, (q + 1) % n);
        }
    }
    let mut probs: Vec<f64> = amps.iter().map(|a| a.mag2()).collect();
    let total: f64 = probs.iter().sum();
    if total > 0.0 {
        for x in probs.iter_mut() {
            *x /= total;
        }
    }
    probs
}

/// Total-variation distance `½·Σ|p−q|` between two distributions.
pub fn tv_distance(p: &[f64], q: &[f64]) -> f64 {
    0.5 * p.iter().zip(q).map(|(a, b)| (a - b).abs()).sum::<f64>()
}

/// TV distance of the `bits`-bit qFHRR distribution from the exact one — the
/// generative model's error on the phasor substrate.
pub fn tv_vs_exact(n: u8, layers: u8, params: &[f64], bits: u32) -> f64 {
    tv_distance(&born_probs(n, layers, params, 0), &born_probs(n, layers, params, bits))
}

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

    fn params(n: u8, layers: u8, seed: u64) -> Vec<f64> {
        let mut st = seed.wrapping_mul(0x2545_F491).wrapping_add(1);
        let mut sm = || {
            st = st.wrapping_add(0x9E37_79B9_7F4A_7C15);
            let mut z = st;
            z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
            z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
            (z ^ (z >> 31)) >> 11
        };
        (0..quantum_vml::num_params(n, layers))
            .map(|_| ((sm() as f64 / (1u64 << 53) as f64) * 2.0 - 1.0) * PI)
            .collect()
    }

    #[test]
    fn exact_mode_matches_the_vml_born_machine() {
        // bits=0 must reproduce the (unquantized) vml Born distribution exactly —
        // proving this is a faithful re-run of the same generative model.
        let (n, l) = (4u8, 3u8);
        let pr = params(n, l, 7);
        let a = born_probs(n, l, &pr, 0);
        let b = quantum_vml::born_probs(n, l, &pr);
        let maxerr = a.iter().zip(&b).map(|(x, y)| (x - y).abs()).fold(0.0, f64::max);
        assert!(maxerr < 1e-12, "exact qFHRR run must equal the vml Born machine: {maxerr}");
    }

    #[test]
    fn full_distribution_normalized() {
        let (n, l) = (4u8, 3u8);
        let pr = params(n, l, 3);
        for bits in [3, 5, 8, 0] {
            let p = born_probs(n, l, &pr, bits);
            assert!((p.iter().sum::<f64>() - 1.0).abs() < 1e-9, "bits={bits} must normalize");
        }
    }

    #[test]
    fn converges_to_exact_as_bits_grow() {
        // the deployable claim: the full generative distribution on the phasor
        // substrate approaches the exact state vector as phase precision grows.
        let (n, l) = (4u8, 4u8);
        let pr = params(n, l, 11);
        let tv3 = tv_vs_exact(n, l, &pr, 3);
        let tv6 = tv_vs_exact(n, l, &pr, 6);
        let tv10 = tv_vs_exact(n, l, &pr, 10);
        assert!(tv6 < tv3, "more phase bits ⇒ closer: {tv6} < {tv3}");
        assert!(tv10 < tv6, "…and closer still: {tv10} < {tv6}");
        assert!(tv10 < 0.01, "10-bit qFHRR runs the model essentially exactly: {tv10}");
        assert!(tv3 < 0.4, "even 3-bit phases give a usable distribution: {tv3}");
    }

    #[test]
    fn deterministic() {
        let pr = params(3, 3, 42);
        let a = born_probs(3, 3, &pr, 4);
        let b = born_probs(3, 3, &pr, 4);
        assert!(a.iter().zip(&b).all(|(x, y)| x.to_bits() == y.to_bits()));
    }
}