gled-model 2.30.0

//! Evaluating one quadratic Bézier segment of a curve at a beat position.
//!
//! The exact value needs a search along the curve. On native targets a
//! polynomial is fitted to each segment once, on a background thread, and used
//! from then on; until it is ready, and always in the browser, the segment is
//! sampled with a binary search.

use super::pos::{CurvePos, quadratic_bezier};
use std::hash::{Hash, Hasher};

#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub(crate) struct BezierCurve {
    pub(crate) start: CurvePos,
    pub(crate) bezier: CurvePos,
    pub(crate) end: CurvePos,
}

impl Hash for BezierCurve {
    fn hash<H: Hasher>(&self, state: &mut H) {
        self.start.x.to_be_bytes().hash(state);
        self.start.y.to_be_bytes().hash(state);
        self.bezier.x.to_be_bytes().hash(state);
        self.bezier.y.to_be_bytes().hash(state);
        self.end.x.to_be_bytes().hash(state);
        self.end.y.to_be_bytes().hash(state);
    }
}

impl BezierCurve {
    pub(crate) fn new(start: CurvePos, bezier: CurvePos, end: CurvePos) -> Self {
        Self { start, bezier, end }
    }

    /// The curve's y at `x`.
    pub(crate) fn value(&self, x: f32) -> f32 {
        #[cfg(not(target_arch = "wasm32"))]
        if let Some(value) = fitted::value(self, x) {
            return value;
        }
        self.sample(x)
    }

    /// The curve's y at `x`, found by a binary search along the curve.
    pub(crate) fn sample(&self, x: f32) -> f32 {
        let points = [self.start, self.bezier, self.end];
        let mut n = 2;
        let mut t = 0.5;
        loop {
            let sample = quadratic_bezier(points, t);
            let e = sample.x - x;
            if e.abs() < 0.005 {
                break sample.y;
            }
            let step = 1.0 / 2.0f32.powi(n);
            if e < 0.0 {
                t += step;
            } else {
                t -= step;
            }
            n += 1;
        }
    }
}

#[cfg(not(target_arch = "wasm32"))]
pub use fitted::{Polynomial, polynomials_fitting};

#[cfg(not(target_arch = "wasm32"))]
mod fitted {
    use super::BezierCurve;
    use crate::asset::curve::pos::quadratic_bezier;
    use kanal::{Sender, unbounded};
    use ndarray::arr1;
    use once_cell::sync::Lazy;
    use polyfit_residuals::{PolyFit, poly::NewtonPolynomial, try_fit_poly_with_residual};
    use rayon::iter::{IntoParallelIterator, ParallelIterator};
    use std::{
        collections::HashMap,
        sync::{
            Mutex, MutexGuard, PoisonError,
            atomic::{AtomicUsize, Ordering::Relaxed},
        },
        thread::spawn,
    };

    pub type Polynomial = NewtonPolynomial<f64, f64, Vec<f64>, Vec<f64>>;

    static POLYNOMIALS_FITTING: AtomicUsize = AtomicUsize::new(0);
    static POLYNOMIALS: Lazy<Mutex<HashMap<BezierCurve, Option<Polynomial>>>> =
        Lazy::new(|| Mutex::new(HashMap::new()));
    static POLYNOMIAL_QUEUE: Lazy<Sender<BezierCurve>> = Lazy::new(|| {
        let (sender, receiver) = unbounded::<BezierCurve>();
        spawn(move || {
            for bezier_curve in receiver {
                // Without a fit the entry stays `None`, so the segment is sampled.
                if let Some(polynomial) = fit(&bezier_curve) {
                    polynomials().insert(bezier_curve, Some(polynomial));
                }
                POLYNOMIALS_FITTING.fetch_sub(1, Relaxed);
            }
        });
        sender
    });

    /// The fitted polynomials. A poisoned lock only means a fit panicked; the
    /// map itself is still consistent.
    fn polynomials() -> MutexGuard<'static, HashMap<BezierCurve, Option<Polynomial>>> {
        POLYNOMIALS.lock().unwrap_or_else(PoisonError::into_inner)
    }

    /// The value from the fitted polynomial, or `None` while it is being
    /// fitted. The first call for a segment starts the fit.
    pub(super) fn value(curve: &BezierCurve, x: f32) -> Option<f32> {
        let mut polynomials = polynomials();
        match polynomials.get(curve) {
            Some(Some(polynomial)) => return Some(polynomial.eval(x as f64) as f32),
            None => {
                POLYNOMIALS_FITTING.fetch_add(1, Relaxed);
                POLYNOMIAL_QUEUE.send(*curve).ok();
                polynomials.insert(*curve, None);
            }
            Some(None) => (),
        }
        None
    }

    /// Number of segments whose polynomial is still being fitted.
    pub fn polynomials_fitting() -> usize {
        POLYNOMIALS_FITTING.load(Relaxed)
    }

    /// Fits polynomials of degree 1 to 5 to the segment and keeps the best.
    /// `None` when no fit succeeded; the segment is then sampled for ever.
    fn fit(curve: &BezierCurve) -> Option<Polynomial> {
        tracing::debug!("Fitting polynomial to curve: {curve:?}");

        let points = [curve.start, curve.bezier, curve.end];
        let (x_values, y_values): (Vec<_>, Vec<_>) = (0..1_000_000)
            .map(|t| t as f32 / 1_000_000.0)
            .map(|t| {
                let pos = quadratic_bezier(points, t);
                (pos.x as f64, pos.y as f64)
            })
            .unzip();

        let x_values = arr1(&x_values);
        let y_values = arr1(&y_values);

        let best = (1..6)
            .into_par_iter()
            .filter_map(|degree| {
                try_fit_poly_with_residual(x_values.view(), y_values.view(), degree)
                    .map_err(|err| {
                        tracing::warn!("Could not fit a degree {degree} polynomial: {err:?}")
                    })
                    .ok()
            })
            .collect::<Vec<_>>()
            .into_iter()
            .min_by(|a: &PolyFit<f64, f64>, b| a.residual.total_cmp(&b.residual))
            .map(|fit| fit.polynomial);

        tracing::debug!("Done fitting polynomial to curve: {curve:?}");
        best
    }
}