added SIMD support

Author: hexqnt

.cargo/config.toml

@@ -1,2 +1,7 @@
 [target.wasm32-unknown-unknown]
-rustflags = ["--cfg", "getrandom_backend=\"wasm_js\""]
+rustflags = [
+    "--cfg",
+    "getrandom_backend=\"wasm_js\"",
+    "-C",
+    "target-feature=+simd128",
+]

.github/workflows/ci.yml

@@ -22,7 +22,7 @@ jobs:
       - name: Setup Rust
         uses: dtolnay/rust-toolchain@v1
         with:
-          toolchain: stable
+          toolchain: nightly
 
       - name: Cache cargo
         uses: Swatinem/rust-cache@v2

.github/workflows/release.yml

@@ -29,7 +29,7 @@ jobs:
       - name: Setup Rust
         uses: dtolnay/rust-toolchain@v1
         with:
-          toolchain: stable
+          toolchain: nightly
 
       - name: Cache cargo
         uses: Swatinem/rust-cache@v2

Cargo.lock

@@ -11,6 +11,10 @@ readme = "README.md"
 keywords = ["curve-fitting", "spline", "optimization", "egui", "wasm"]
 categories = ["science", "mathematics", "visualization", "gui"]
 
+[features]
+default = ["portable-simd"]
+portable-simd = []
+
 [dependencies]
 argmin = { version = "0.11.0", default-features = false }
 argmin-math = { version = "0.5.1", default-features = false, features = [

Cargo.toml

@@ -0,0 +1,2 @@
+[toolchain]
+channel = "nightly"

rust-toolchain.toml

@@ -13,10 +13,11 @@ use argmin::solver::quasinewton::LBFGS;
 
 use crate::domain::{
     CurveFamily, CurveParams, FitResult, InputError, LbfgsConfig, NelderMeadConfig,
-    OptimizerConfig, Point, Points, SteepestDescentConfig,
+    OptimizerConfig, Points, SteepestDescentConfig,
 };
 
 mod curve;
+mod simd;
 mod spline;
 
 #[cfg(not(target_arch = "wasm32"))]
@@ -34,6 +35,94 @@ pub use spline::{
     fit_natural_cubic_spline_with_config, fit_natural_cubic_spline_with_optimizer_config,
 };
 
+#[cfg(feature = "portable-simd")]
+#[doc(hidden)]
+pub mod simd_bench {
+    use super::{OptimizationLossMetric, simd};
+
+    pub fn polynomial_cost_scalar(
+        param: &[f64],
+        x_values: &[f64],
+        y_values: &[f64],
+        loss_metric: OptimizationLossMetric,
+    ) -> f64 {
+        simd::polynomial_cost_scalar(param, x_values, y_values, loss_metric)
+    }
+
+    pub fn polynomial_cost_simd(
+        param: &[f64],
+        x_values: &[f64],
+        y_values: &[f64],
+        loss_metric: OptimizationLossMetric,
+    ) -> f64 {
+        simd::polynomial_cost_simd(param, x_values, y_values, loss_metric)
+    }
+
+    pub fn inverse_cost_scalar(
+        param: &[f64],
+        x_values: &[f64],
+        y_values: &[f64],
+        loss_metric: OptimizationLossMetric,
+    ) -> f64 {
+        simd::inverse_cost_scalar(param, x_values, y_values, loss_metric)
+    }
+
+    pub fn inverse_cost_simd(
+        param: &[f64],
+        x_values: &[f64],
+        y_values: &[f64],
+        loss_metric: OptimizationLossMetric,
+    ) -> f64 {
+        simd::inverse_cost_simd(param, x_values, y_values, loss_metric)
+    }
+
+    pub fn polynomial_gradient_scalar(
+        x_values: &[f64],
+        y_values: &[f64],
+        param: &[f64],
+        loss_metric: OptimizationLossMetric,
+        gradient: &mut [f64],
+    ) {
+        simd::accumulate_polynomial_gradient_scalar(
+            x_values,
+            y_values,
+            param,
+            loss_metric,
+            gradient,
+        );
+    }
+
+    pub fn polynomial_gradient_simd(
+        x_values: &[f64],
+        y_values: &[f64],
+        param: &[f64],
+        loss_metric: OptimizationLossMetric,
+        gradient: &mut [f64],
+    ) {
+        simd::accumulate_polynomial_gradient_simd(x_values, y_values, param, loss_metric, gradient);
+    }
+
+    pub fn inverse_gradient_scalar(
+        x_values: &[f64],
+        y_values: &[f64],
+        param: &[f64],
+        loss_metric: OptimizationLossMetric,
+        gradient: &mut [f64],
+    ) {
+        simd::accumulate_inverse_gradient_scalar(x_values, y_values, param, loss_metric, gradient);
+    }
+
+    pub fn inverse_gradient_simd(
+        x_values: &[f64],
+        y_values: &[f64],
+        param: &[f64],
+        loss_metric: OptimizationLossMetric,
+        gradient: &mut [f64],
+    ) {
+        simd::accumulate_inverse_gradient_simd(x_values, y_values, param, loss_metric, gradient);
+    }
+}
+
 const PARAM_EPS: f64 = 1e-9;
 const LARGE_COST: f64 = 1e24;
 const LN_2: f64 = std::f64::consts::LN_2;
@@ -174,42 +263,29 @@ pub struct IterationMetricSnapshot {
 struct CurveProblem {
     family: CurveFamily,
     points: Points,
+    point_x: Box<[f64]>,
+    point_y: Box<[f64]>,
     loss_metric: OptimizationLossMetric,
 }
 
 impl CurveProblem {
     fn new(family: CurveFamily, points: &Points, loss_metric: OptimizationLossMetric) -> Self {
+        let mut point_x = Vec::with_capacity(points.len());
+        let mut point_y = Vec::with_capacity(points.len());
+        for point in points.as_slice() {
+            point_x.push(point.x());
+            point_y.push(point.y());
+        }
         Self {
             family,
             points: points.clone(),
+            point_x: point_x.into_boxed_slice(),
+            point_y: point_y.into_boxed_slice(),
             loss_metric,
         }
     }
 }
 
-fn accumulate_polynomial_gradient(
-    points: &[Point],
-    param: &[f64],
-    loss_metric: OptimizationLossMetric,
-    gradient: &mut [f64],
-) {
-    debug_assert_eq!(gradient.len(), param.len());
-    for point in points {
-        let x = point.x();
-        let model = param
-            .iter()
-            .copied()
-            .fold(0.0, |acc, coefficient| acc * x + coefficient);
-        let residual = loss_metric.residual_derivative(model - point.y());
-
-        let mut basis = 1.0;
-        for gradient_value in gradient.iter_mut().rev() {
-            *gradient_value += residual * basis;
-            basis *= x;
-        }
-    }
-}
-
 #[derive(Debug, Clone, PartialEq)]
 /// Подробный результат подгонки сплайна.
 pub struct SplineResult {
@@ -259,6 +335,21 @@ fn soft_l1_from_residuals(residuals: &[[f64; 2]]) -> f64 {
         / residuals.len() as f64
 }
 
+fn is_polynomial_family(family: CurveFamily) -> bool {
+    matches!(
+        family,
+        CurveFamily::Linear
+            | CurveFamily::Quadratic
+            | CurveFamily::Cubic
+            | CurveFamily::Quartic
+            | CurveFamily::Quintic
+            | CurveFamily::Sextic
+            | CurveFamily::Septic
+            | CurveFamily::Octic
+            | CurveFamily::Nonic
+    )
+}
+
 /// Число узлов сплайна по умолчанию.
 pub const DEFAULT_SPLINE_KNOTS: usize = 8;
 /// Число сэмплов кривой для визуализации по умолчанию.
@@ -1069,6 +1160,24 @@ impl CostFunction for CurveProblem {
     type Output = f64;
 
     fn cost(&self, param: &Self::Param) -> Result<Self::Output, argmin::core::Error> {
+        if is_polynomial_family(self.family) {
+            return Ok(simd::polynomial_cost(
+                param,
+                self.point_x.as_ref(),
+                self.point_y.as_ref(),
+                self.loss_metric,
+            ));
+        }
+
+        if self.family == CurveFamily::Inverse {
+            return Ok(simd::inverse_cost(
+                param,
+                self.point_x.as_ref(),
+                self.point_y.as_ref(),
+                self.loss_metric,
+            ));
+        }
+
         let sample_count = self.points.len() as f64;
         let mut sum = 0.0;
 
@@ -1106,9 +1215,13 @@ impl Gradient for CurveProblem {
             | CurveFamily::Sextic
             | CurveFamily::Septic
             | CurveFamily::Octic
-            | CurveFamily::Nonic => {
-                accumulate_polynomial_gradient(points, param, self.loss_metric, &mut gradient);
-            }
+            | CurveFamily::Nonic => simd::accumulate_polynomial_gradient(
+                self.point_x.as_ref(),
+                self.point_y.as_ref(),
+                param,
+                self.loss_metric,
+                &mut gradient,
+            ),
             CurveFamily::Arrhenius => {
                 for point in points {
                     let x = positive_x(point.x());
@@ -1119,16 +1232,13 @@ impl Gradient for CurveProblem {
                     gradient[1] += residual * (param[0] * exp_term / x);
                 }
             }
-            CurveFamily::Inverse => {
-                for point in points {
-                    let x = positive_x(point.x());
-                    let residual = self
-                        .loss_metric
-                        .residual_derivative(param[0] + param[1] / x - point.y());
-                    gradient[0] += residual;
-                    gradient[1] += residual / x;
-                }
-            }
+            CurveFamily::Inverse => simd::accumulate_inverse_gradient(
+                self.point_x.as_ref(),
+                self.point_y.as_ref(),
+                param,
+                self.loss_metric,
+                &mut gradient,
+            ),
             CurveFamily::Logistic => {
                 for point in points {
                     let z = param[1] * (point.x() - param[2]);