//! Phasor-feature classifier — the quantum kernel, deployed (`wai.quantum.kernel`).
//!
//! A quantum kernel `K(x,x') = |⟨ψ(x)|ψ(x')⟩|²` is, by the Schuld–Sweke–Meyer
//! identity, a shift-invariant kernel whose Fourier spectrum is set by the data
//! encoding — and **random Fourier features** are its classical realization: draw
//! random frequencies `ω`, map `x ↦ φ(x) = √(2/D)·cos(ω·x + b)` (a bank of phasors,
//! the FHRR feature map), and `⟨φ(x),φ(x')⟩ ≈ K(x,x')` (Rahimi–Recht). The
//! dequantization theorems say this matches the quantum kernel whenever the
//! spectrum is compact. So the "quantum kernel" ships as an ordinary phasor layer:
//! train a linear model on the phasor features and it separates data a linear model
//! cannot — the supervised counterpart to the generative Born machine, deployable
//! on a laptop, and at **qFHRR** phase resolution for the edge.
//!
//! Every phasor evaluation is a counted operation (`feature_ops`) — the hook for
//! pricing interference computing in joules. Deterministic f64.

// Indexed loops are deliberate in the ridge-fit / linear-solve matrix code.
#![allow(clippy::needless_range_loop)]

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

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
}
/// Standard normal via Box–Muller (deterministic f64).
fn gauss(s: &mut u64) -> f64 {
    let u1 = u01(s).max(1e-300);
    let u2 = u01(s);
    (-2.0 * crate::repro::ln(u1)).sqrt() * crate::repro::sin_cos(2.0 * PI * u2).1
}

/// A random-Fourier-feature (phasor) map approximating a Gaussian kernel
/// `exp(-γ·|x−x'|²)`.
pub struct Rff {
    /// `omegas[i]` is a `dim`-vector frequency; `phases[i]` a bias in `[0,2π)`.
    omegas: Vec<Vec<f64>>,
    phases: Vec<f64>,
    scale: f64,
}

impl Rff {
    pub fn dim_features(&self) -> usize {
        self.phases.len()
    }
}

/// Sample a phasor feature map: `d_features` frequencies `ω ~ N(0, 2γ·I)` and
/// uniform phase biases. Deterministic in `seed`.
pub fn sample_rff(dim: usize, d_features: usize, gamma: f64, seed: u64) -> Rff {
    let mut st = seed.wrapping_mul(0xA24B_AED4).wrapping_add(1);
    let sd = (2.0 * gamma).sqrt();
    let omegas = (0..d_features)
        .map(|_| (0..dim).map(|_| gauss(&mut st) * sd).collect())
        .collect();
    let phases = (0..d_features).map(|_| u01(&mut st) * 2.0 * PI).collect();
    Rff { omegas, phases, scale: (2.0 / d_features as f64).sqrt() }
}

/// The phasor feature vector `φ(x)`.
pub fn features(rff: &Rff, x: &[f64]) -> Vec<f64> {
    rff.omegas
        .iter()
        .zip(&rff.phases)
        .map(|(w, &b)| {
            let dot: f64 = w.iter().zip(x).map(|(wi, xi)| wi * xi).sum();
            rff.scale * crate::repro::sin_cos(dot + b).1
        })
        .collect()
}

/// `φ(x)` with each phasor phase snapped to `bits`-bit qFHRR resolution.
pub fn features_quantized(rff: &Rff, x: &[f64], bits: u32) -> Vec<f64> {
    let levels = (1u64 << bits) as f64;
    rff.omegas
        .iter()
        .zip(&rff.phases)
        .map(|(w, &b)| {
            let dot: f64 = w.iter().zip(x).map(|(wi, xi)| wi * xi).sum();
            let t = dot + b;
            let tq = (t / (2.0 * PI) * levels).round() / levels * 2.0 * PI;
            rff.scale * crate::repro::sin_cos(tq).1
        })
        .collect()
}

/// Solve `A w = b` for a symmetric positive-definite `A` (Gaussian elimination
/// with partial pivoting). `a` is consumed.
fn solve(mut a: Vec<Vec<f64>>, mut b: Vec<f64>) -> Vec<f64> {
    let n = b.len();
    for col in 0..n {
        // partial pivot
        let mut piv = col;
        for r in (col + 1)..n {
            if a[r][col].abs() > a[piv][col].abs() {
                piv = r;
            }
        }
        a.swap(col, piv);
        b.swap(col, piv);
        let d = a[col][col];
        if d.abs() < 1e-15 {
            continue;
        }
        for r in (col + 1)..n {
            let f = a[r][col] / d;
            if f != 0.0 {
                for c in col..n {
                    a[r][c] -= f * a[col][c];
                }
                b[r] -= f * b[col];
            }
        }
    }
    let mut w = vec![0.0; n];
    for i in (0..n).rev() {
        let mut s = b[i];
        for c in (i + 1)..n {
            s -= a[i][c] * w[c];
        }
        w[i] = if a[i][i].abs() < 1e-15 { 0.0 } else { s / a[i][i] };
    }
    w
}

/// Fit a ridge classifier on phasor features: `w = (ΦᵀΦ + λI)⁻¹ Φᵀ y`. Labels `y`
/// in `{−1,+1}`; predict `sign(w·φ(x))`.
pub fn fit_ridge(rff: &Rff, x: &[Vec<f64>], y: &[f64], lambda: f64) -> Vec<f64> {
    let d = rff.dim_features();
    let mut a = vec![vec![0.0; d]; d];
    let mut b = vec![0.0; d];
    for (xi, &yi) in x.iter().zip(y) {
        let phi = features(rff, xi);
        for i in 0..d {
            b[i] += phi[i] * yi;
            for j in i..d {
                a[i][j] += phi[i] * phi[j];
            }
        }
    }
    for i in 0..d {
        for j in 0..i {
            a[i][j] = a[j][i];
        }
        a[i][i] += lambda;
    }
    solve(a, b)
}

/// The raw decision score `w·φ(x)`.
pub fn score(rff: &Rff, w: &[f64], x: &[f64]) -> f64 {
    features(rff, x).iter().zip(w).map(|(f, wi)| f * wi).sum()
}

/// Classification accuracy on `(x, y)`.
pub fn accuracy(rff: &Rff, w: &[f64], x: &[Vec<f64>], y: &[f64]) -> f64 {
    let mut correct = 0usize;
    for (xi, yi) in x.iter().zip(y) {
        if score(rff, w, xi).signum() == yi.signum() {
            correct += 1;
        }
    }
    correct as f64 / x.len() as f64
}

/// Accuracy using `bits`-bit qFHRR features (the edge-deployable path).
pub fn accuracy_quantized(rff: &Rff, w: &[f64], x: &[Vec<f64>], y: &[f64], bits: u32) -> f64 {
    let mut correct = 0usize;
    for (xi, yi) in x.iter().zip(y) {
        let s: f64 = features_quantized(rff, xi, bits).iter().zip(w).map(|(f, wi)| f * wi).sum();
        if s.signum() == yi.signum() {
            correct += 1;
        }
    }
    correct as f64 / x.len() as f64
}

/// Phasor operations to featurize `n_samples`: `D · (dim + 1)` per sample — the
/// count a joule meter prices.
pub fn feature_ops(rff: &Rff, n_samples: usize) -> u64 {
    let dim = rff.omegas.first().map(|w| w.len()).unwrap_or(0);
    (rff.dim_features() as u64) * (dim as u64 + 1) * n_samples as u64
}

// ---------------------------------------------------------------------------
// Deterministic 2-D benchmark datasets (nonlinearly separable)
// ---------------------------------------------------------------------------

/// Concentric circles: inner disc (−1) vs outer ring (+1). Not linearly separable.
pub fn circles(n: usize, seed: u64) -> (Vec<Vec<f64>>, Vec<f64>) {
    let mut st = seed.wrapping_mul(0x2545_F491).wrapping_add(1);
    let mut x = Vec::with_capacity(n);
    let mut y = Vec::with_capacity(n);
    for i in 0..n {
        let inner = i % 2 == 0;
        let ang = u01(&mut st) * 2.0 * PI;
        let r = if inner { 0.35 * u01(&mut st) } else { 0.9 + 0.35 * u01(&mut st) };
        let (sin, cos) = crate::repro::sin_cos(ang);
        x.push(vec![r * cos + 0.03 * gauss(&mut st), r * sin + 0.03 * gauss(&mut st)]);
        y.push(if inner { -1.0 } else { 1.0 });
    }
    (x, y)
}

/// XOR clusters at `(±1,±1)`, label `sign(x1·x2)`. Not linearly separable.
pub fn xor_data(n: usize, seed: u64) -> (Vec<Vec<f64>>, Vec<f64>) {
    let mut st = seed.wrapping_mul(0x9E37_79B9).wrapping_add(3);
    let mut x = Vec::with_capacity(n);
    let mut y = Vec::with_capacity(n);
    for i in 0..n {
        let sx = if (i & 1) == 0 { 1.0 } else { -1.0 };
        let sy = if (i & 2) == 0 { 1.0 } else { -1.0 };
        x.push(vec![sx + 0.28 * gauss(&mut st), sy + 0.28 * gauss(&mut st)]);
        y.push(sx * sy);
    }
    (x, y)
}

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

    #[test]
    fn phasor_kernel_separates_circles() {
        let (xtr, ytr) = circles(300, 1);
        let (xte, yte) = circles(300, 2);
        let rff = sample_rff(2, 200, 4.0, 7);
        let w = fit_ridge(&rff, &xtr, &ytr, 1e-3);
        let acc = accuracy(&rff, &w, &xte, &yte);
        assert!(acc > 0.9, "phasor kernel must separate circles: {acc}");
    }

    #[test]
    fn phasor_kernel_separates_xor() {
        let (xtr, ytr) = xor_data(320, 5);
        let (xte, yte) = xor_data(320, 9);
        let rff = sample_rff(2, 160, 1.5, 3);
        let w = fit_ridge(&rff, &xtr, &ytr, 1e-3);
        assert!(accuracy(&rff, &w, &xte, &yte) > 0.9, "phasor kernel must solve XOR");
    }

    #[test]
    fn a_linear_model_fails_on_circles() {
        // a degenerate RFF with a single low-frequency feature ≈ linear — should fail
        // on a problem the full phasor bank solves, showing the features are the point.
        let (xtr, ytr) = circles(300, 1);
        let (xte, yte) = circles(300, 2);
        // near-linear: 1 feature, tiny gamma
        let rff = sample_rff(2, 2, 1e-4, 7);
        let w = fit_ridge(&rff, &xtr, &ytr, 1e-3);
        assert!(accuracy(&rff, &w, &xte, &yte) < 0.7, "a linear model can't separate circles");
    }

    #[test]
    fn qfhrr_quantized_features_still_classify() {
        let (xtr, ytr) = circles(300, 1);
        let (xte, yte) = circles(300, 2);
        let rff = sample_rff(2, 200, 4.0, 7);
        let w = fit_ridge(&rff, &xtr, &ytr, 1e-3);
        let full = accuracy(&rff, &w, &xte, &yte);
        let q4 = accuracy_quantized(&rff, &w, &xte, &yte, 4);
        assert!(q4 > 0.85, "4-bit qFHRR features still classify: {q4} (full {full})");
        assert!(accuracy_quantized(&rff, &w, &xte, &yte, 8) >= q4 - 0.02, "more bits don't hurt");
    }

    #[test]
    fn deterministic() {
        let (x, y) = circles(50, 1);
        let rff = sample_rff(2, 32, 3.0, 7);
        let a = fit_ridge(&rff, &x, &y, 1e-3);
        let b = fit_ridge(&rff, &x, &y, 1e-3);
        assert!(a.iter().zip(&b).all(|(p, q)| p.to_bits() == q.to_bits()));
    }
}