feat: add simd uniform

Author: dancixx

.gitignore

@@ -5,4 +5,5 @@ flamegraph.svg
 venv
 CI.yml
 /test
-.idea/
+.idea/
+.zed

Cargo.toml

@@ -81,6 +81,7 @@ default = []
 jemalloc = ["dep:tikv-jemallocator"]
 malliavin = []
 mimalloc = ["dep:mimalloc"]
+simd = []
 yahoo = ["dep:time", "dep:yahoo_finance_api", "dep:polars"]
 
 [lib]

src/lib.rs

@@ -3,7 +3,8 @@
 #![allow(clippy::type_complexity)]
 #![allow(clippy::too_many_arguments)]
 //#![warn(missing_docs)]
-
+//
+#![feature(portable_simd)]
 // TODO: this is just temporary
 #![allow(dead_code)]
 

src/stats/distr.rs

@@ -35,7 +35,8 @@ mod tests {
     beta::SimdBeta, binomial::SimdBinomial, cauchy::SimdCauchy, gamma::SimdGamma,
     geometric::SimdGeometric, hypergeometric::SimdHypergeometric, inverse_gauss::SimdInverseGauss,
     lognormal::SimdLogNormal, normal::SimdNormal, normal_inverse_gauss::SimdNormalInverseGauss,
-    pareto::SimdPareto, poisson::SimdPoisson, studentt::SimdStudentT, weibull::SimdWeibull,
+    pareto::SimdPareto, poisson::SimdPoisson, studentt::SimdStudentT, uniform::SimdUniform,
+    weibull::SimdWeibull,
   };
 
   use plotly::{
@@ -386,6 +387,22 @@ mod tests {
       plot.add_trace(trace);
     }
 
+    // 15) Uniform (0,1)
+    {
+      let (xa, ya) = subplot_axes(4, 3);
+      let mut rng = thread_rng();
+      let dist = SimdUniform::new(0.0, 1.0);
+      let samples: Vec<f32> = (0..sample_size).map(|_| dist.sample(&mut rng)).collect();
+      let (xs, bins) = make_histogram(&samples, 100, -4.0, 4.0);
+      let trace = Scatter::new(xs, bins)
+        .name("Uniform(0,1)")
+        .mode(Mode::Lines)
+        .line(Line::new().shape(LineShape::Linear))
+        .x_axis(&xa)
+        .y_axis(&ya);
+      plot.add_trace(trace);
+    }
+
     plot.show();
   }
 }

src/stochastic.rs

@@ -82,6 +82,12 @@ pub trait SamplingExt<T: Clone + Send + Sync + Zero>: Send + Sync {
   /// Sample the process
   fn sample(&self) -> Array1<T>;
 
+  /// Sample the process with simd acceleration
+  #[cfg(feature = "simd")]
+  fn sample_simd(&self) -> Array1<T> {
+    unimplemented!()
+  }
+
   /// Sample the process with CUDA support
   #[cfg(feature = "cuda")]
   fn sample_cuda(&self) -> Result<Either<Array1<T>, Array2<T>>> {

src/stochastic/diffusion/ou.rs

@@ -32,6 +32,25 @@ impl SamplingExt<f64> for OU {
     ou
   }
 
+  #[cfg(feature = "simd")]
+  // TODO: experimental
+  fn sample_simd(&self) -> Array1<f64> {
+    use crate::stats::distr::normal::SimdNormal;
+
+    let dt = self.t.unwrap_or(1.0) as f32 / (self.n - 1) as f32;
+    let gn = Array1::random(self.n, SimdNormal::new(0.0, dt.sqrt()));
+
+    let mut ou = Array1::<f64>::zeros(self.n);
+    ou[0] = self.x0.unwrap_or(0.0);
+
+    for i in 1..self.n {
+      ou[i] =
+        ou[i - 1] + self.theta * (self.mu - ou[i - 1]) * dt as f64 + self.sigma * gn[i - 1] as f64
+    }
+
+    ou
+  }
+
   /// Number of time steps
   fn n(&self) -> usize {
     self.n
@@ -73,6 +92,32 @@ mod tests {
     plot_1d!(ou.sample(), "Fractional Ornstein-Uhlenbeck (FOU) Process");
   }
 
+  #[cfg(feature = "simd")]
+  #[test]
+  fn sample_simd() {
+    use std::time::Instant;
+
+    let start = Instant::now();
+    let ou = OU::new(2.0, 1.0, 0.8, N, Some(X0), Some(1.0), None);
+
+    for _ in 0..100_000 {
+      ou.sample_simd();
+    }
+
+    let elapsed = start.elapsed();
+    println!("Elapsed time for sample_simd: {:?}", elapsed);
+
+    let start = Instant::now();
+    let ou = OU::new(2.0, 1.0, 0.8, N, Some(X0), Some(1.0), None);
+
+    for _ in 0..100_000 {
+      ou.sample();
+    }
+
+    let elapsed = start.elapsed();
+    println!("Elapsed time for sample: {:?}", elapsed);
+  }
+
   #[test]
   #[ignore = "Not implemented"]
   #[cfg(feature = "malliavin")]