open-standards/wai/wai-quantum/src/quantum_phasor.rs

//! The phasor bridge — same interference, two substrates (`wai.quantum.phasor`).
//!
//! Schuld, Sweke & Meyer (2020) proved that a quantum model with data encoded by
//! Pauli-rotation gates is *exactly* a truncated Fourier series of phasors:
//!
//! ```text
//!   f(x) = ⟨0| U†(x) M U(x) |0⟩ = Σ_{ω ∈ Ω} c_ω · e^{iωx}
//! ```
//!
//! where the encoding fixes the frequency set `Ω` and the trainable gates fix the
//! complex coefficients `c_ω`. That is the whole bridge: a **quantum feature map is
//! a phasor bank**, and phasor binding (complex multiplication = phase addition) is
//! *interference*. So the identical function can be computed two ways — on a
//! **quantum state vector**, or by a **classical phasor sum** (the FHRR / HDC
//! algebra) — and where the spectrum is compact (the common case, per the
//! dequantization theorems) the classical phasor substrate reproduces it with none
//! of the qubits. The **quantized-phase** variant is qFHRR: the same phasors at a
//! few bits of phase, the cheapest deployable realization of what the circuit
//! computes.
//!
//! This module builds a single-qubit **data-reuploading** model, extracts its
//! Fourier spectrum, and shows the phasor sum reproduces it exactly — the quantum
//! end and the deployable end computing one interference function. Deterministic
//! f64.

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

#[derive(Clone, Copy)]
struct C {
    re: f64,
    im: f64,
}
impl C {
    #[inline]
    fn new(re: f64, im: f64) -> C {
        C { re, im }
    }
    #[inline]
    fn mul(self, o: C) -> C {
        C { re: self.re * o.re - self.im * o.im, im: self.re * o.im + self.im * o.re }
    }
    #[inline]
    fn add(self, o: C) -> C {
        C { re: self.re + o.re, im: self.im + o.im }
    }
}
#[inline]
fn expi(t: f64) -> C {
    C::new(t.cos(), t.sin())
}

// ---------------------------------------------------------------------------
// The quantum model: single-qubit data reuploading  ⟨Z⟩ = f(x)
// ---------------------------------------------------------------------------

/// Trainable block `W = RZ(φ)·RY(θ)` applied to a single-qubit amplitude pair.
fn apply_w(a: &mut [C; 2], theta: f64, phi: f64) {
    // RY(θ)
    let (c, s) = ((theta * 0.5).cos(), (theta * 0.5).sin());
    let (a0, a1) = (a[0], a[1]);
    let r0 = C::new(a0.re * c - a1.re * s, a0.im * c - a1.im * s);
    let r1 = C::new(a0.re * s + a1.re * c, a0.im * s + a1.im * c);
    // RZ(φ)
    a[0] = r0.mul(expi(-phi * 0.5));
    a[1] = r1.mul(expi(phi * 0.5));
}
/// Data-encoding gate `S(x) = RZ(x)` (Pauli-Z generator → integer frequencies).
fn apply_encode(a: &mut [C; 2], x: f64) {
    a[0] = a[0].mul(expi(-x * 0.5));
    a[1] = a[1].mul(expi(x * 0.5));
}

/// The single-qubit data-reuploading model with `r` encoding repetitions and
/// `params.len() == 2·(r+1)` trainable angles (θ,φ per block): `f(x) = ⟨Z⟩`, a real
/// trig polynomial of degree `r`.
pub fn model_value(x: f64, r: usize, params: &[f64]) -> f64 {
    let mut a = [C::new(1.0, 0.0), C::new(0.0, 0.0)];
    for i in 0..r {
        apply_w(&mut a, params[2 * i], params[2 * i + 1]);
        apply_encode(&mut a, x);
    }
    apply_w(&mut a, params[2 * r], params[2 * r + 1]);
    // ⟨Z⟩ = |a0|² − |a1|²
    (a[0].re * a[0].re + a[0].im * a[0].im) - (a[1].re * a[1].re + a[1].im * a[1].im)
}

// ---------------------------------------------------------------------------
// Spectrum extraction (the model IS a Fourier series) + phasor evaluation
// ---------------------------------------------------------------------------

/// Extract the Fourier coefficients `c_k` for `k = 0..=r` (the model has spectrum
/// `{−r..r}` with Hermitian symmetry `c_{−k} = c_k*`). Since `f` is a degree-`r`
/// trig polynomial, `2r+1` samples recover it exactly (a DFT).
pub fn spectrum(r: usize, params: &[f64]) -> Vec<(f64, f64)> {
    let n = 2 * r + 1;
    let samples: Vec<f64> = (0..n).map(|j| model_value(2.0 * PI * j as f64 / n as f64, r, params)).collect();
    (0..=r)
        .map(|k| {
            let mut acc = C::new(0.0, 0.0);
            for (j, &f) in samples.iter().enumerate() {
                let ph = -(k as f64) * 2.0 * PI * j as f64 / n as f64;
                acc = acc.add(C::new(f, 0.0).mul(expi(ph)));
            }
            (acc.re / n as f64, acc.im / n as f64)
        })
        .collect()
}

/// Evaluate the model as a **phasor sum** — the FHRR / HDC form — from its
/// coefficients: `f(x) = c_0 + 2·Σ_{k≥1} Re(c_k · e^{ikx})`. Pure classical
/// interference (complex phasor addition), no qubits.
pub fn phasor_value(x: f64, coeffs: &[(f64, f64)]) -> f64 {
    let mut acc = coeffs[0].0; // c_0 (real)
    for (k, &(re, im)) in coeffs.iter().enumerate().skip(1) {
        let e = expi(k as f64 * x);
        acc += 2.0 * (re * e.re - im * e.im); // 2·Re(c_k e^{ikx})
    }
    acc
}

/// The **qFHRR** form: the same phasor sum with each phase quantized to `bits`
/// bits (`2^bits` phase levels) — integer-native, the cheapest deployable
/// realization. Approximates [`phasor_value`] with an error bounded by the phase
/// resolution.
pub fn phasor_value_quantized(x: f64, coeffs: &[(f64, f64)], bits: u32) -> f64 {
    let levels = (1u64 << bits) as f64;
    let quant = |t: f64| -> f64 {
        let frac = t / (2.0 * PI);
        (frac * levels).round() / levels * 2.0 * PI
    };
    let mut acc = coeffs[0].0;
    for (k, &(re, im)) in coeffs.iter().enumerate().skip(1) {
        let e = expi(quant(k as f64 * x));
        acc += 2.0 * (re * e.re - im * e.im);
    }
    acc
}

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)
}
/// Deterministic random trainable params for a demo model of `r` reuploads.
pub fn random_params(r: usize, seed: u64) -> Vec<f64> {
    let mut st = seed.wrapping_mul(0x2545_F491).wrapping_add(1);
    (0..2 * (r + 1))
        .map(|_| {
            let u = (splitmix64(&mut st) >> 11) as f64 / (1u64 << 53) as f64;
            (u * 2.0 - 1.0) * PI
        })
        .collect()
}

/// Sample the quantum model and the (full-precision) phasor model on a grid, plus
/// the max discrepancy — the proof they are one function.
pub fn compare(r: usize, params: &[f64], grid: usize) -> (Vec<f64>, Vec<f64>, Vec<f64>, f64) {
    let coeffs = spectrum(r, params);
    let mut xs = Vec::with_capacity(grid);
    let mut q = Vec::with_capacity(grid);
    let mut ph = Vec::with_capacity(grid);
    let mut maxerr = 0.0f64;
    for j in 0..grid {
        let x = 2.0 * PI * j as f64 / grid as f64;
        let qv = model_value(x, r, params);
        let pv = phasor_value(x, &coeffs);
        xs.push(x);
        q.push(qv);
        ph.push(pv);
        maxerr = maxerr.max((qv - pv).abs());
    }
    (xs, q, ph, maxerr)
}

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

    #[test]
    fn phasor_sum_reproduces_the_quantum_model_exactly() {
        // the core claim: the quantum feature map IS its Fourier/phasor series.
        for (r, seed) in [(2usize, 1u64), (3, 7), (4, 11), (5, 3)] {
            let params = random_params(r, seed);
            let (_, _, _, maxerr) = compare(r, &params, 200);
            assert!(maxerr < 1e-9, "r={r}: phasor sum must equal the quantum model, err={maxerr}");
        }
    }

    #[test]
    fn spectrum_is_compact_and_hermitian_real_output() {
        // r reuploads → exactly r+1 stored coeffs (c_0..c_r); the output is real.
        let r = 4;
        let params = random_params(r, 5);
        let coeffs = spectrum(r, &params);
        assert_eq!(coeffs.len(), r + 1);
        // c_0 is real (imag ~ 0) since f is real
        assert!(coeffs[0].1.abs() < 1e-9, "c_0 must be real");
        // model output is in [-1, 1] (it's ⟨Z⟩)
        for j in 0..50 {
            let v = model_value(j as f64 * 0.1, r, &params);
            assert!(v <= 1.0 + 1e-9 && v >= -1.0 - 1e-9);
        }
    }

    #[test]
    fn qfhrr_quantized_phasors_approximate_with_bounded_error() {
        let r = 3;
        let params = random_params(r, 9);
        let coeffs = spectrum(r, &params);
        // more bits → smaller error, monotonically
        let err = |bits: u32| {
            let mut m = 0.0f64;
            for j in 0..400 {
                let x = 2.0 * PI * j as f64 / 400.0;
                m = m.max((phasor_value(x, &coeffs) - phasor_value_quantized(x, &coeffs, bits)).abs());
            }
            m
        };
        let e3 = err(3);
        let e6 = err(6);
        assert!(e6 < e3, "more phase bits ⇒ smaller error: {e6} < {e3}");
        assert!(e3 < 1.0, "even 3-bit qFHRR phases track the model: {e3}");
        assert!(err(10) < 0.02, "10-bit qFHRR is essentially exact");
    }

    #[test]
    fn deterministic() {
        let a = compare(3, &random_params(3, 42), 64).3;
        let b = compare(3, &random_params(3, 42), 64).3;
        assert_eq!(a.to_bits(), b.to_bits());
    }
}