feat: improve performance

Cargo.toml

@@ -15,13 +15,14 @@ linreg = "0.2.0"
 nalgebra = { version = "0.32.2", features = ["rayon"] }
 ndarray = "0.15.6"
 ndrustfft = "0.4.1"
-num-complex = "0.4.4"
+num-complex = { version = "0.4.4", features = ["rand"] }
 rand = "0.8.5"
 rand_distr = "0.4.3"
 rayon = "1.8.0"
 rustfft = "6.1.0"
 indicatif = "0.17.7"
 plotly = "0.8.4"
+ndarray-rand = "0.14.0"
 
 [dev-dependencies]
 

src/main.rs

@@ -2,23 +2,28 @@ use std::time::Instant;
 
 use indicatif::ProgressBar;
 use plotly::{Plot, Scatter};
-use stochastic_rs::{prelude::*, processes::fbm::Fbm, statistics::fractal_dim::higuchi_fd};
+use stochastic_rs::{
+  prelude::{fgn::FgnFft, *},
+  processes::fbm::Fbm,
+  //processes::fbm::Fbm,
+  statistics::fractal_dim::higuchi_fd,
+};
 
 fn main() {
   let start = Instant::now();
   let fbm = Fbm::new(0.7, 10000, None, Some(10000), None);
-  let m = 10;
-  let pb = ProgressBar::new(m);
+  let m = 1;
+  // let pb = ProgressBar::new(m);
   // let mut plot = Plot::new();
   for _ in 0..m {
     let path = fbm.sample();
-    // let h = higuchi_fd(&path, 10);
-    // println!("Higuchi FD: {}", 2.0 - h);
-    // plot.add_trace(Scatter::new((0..5000).collect::<Vec<usize>>(), path));
-    // plot.show();
-    pb.inc(1);
+    // // let h = higuchi_fd(&path, 10);
+    // // println!("Higuchi FD: {}", 2.0 - h);
+    // // plot.add_trace(Scatter::new((0..10000).collect::<Vec<usize>>(), path));
+    // // plot.show();
+    //   pb.inc(1);
   }
-  pb.finish();
+  // pb.finish();
   println!("Time elapsed: {:?}", start.elapsed().as_secs_f64());
 
   // let start = Instant::now();

src/noises/fgn.rs

@@ -1,8 +1,9 @@
 use crate::utils::Generator;
 use nalgebra::{DMatrix, DVector, Dim, Dyn, RowDVector};
 use ndarray::{concatenate, prelude::*};
+use ndarray_rand::RandomExt;
 use ndrustfft::{ndfft_par, FftHandler};
-use num_complex::Complex;
+use num_complex::{Complex, ComplexDistribution};
 use rand::{thread_rng, Rng};
 use rand_distr::StandardNormal;
 use rayon::prelude::*;
@@ -12,7 +13,7 @@ pub struct FgnFft {
   hurst: f64,
   n: usize,
   t: f64,
-  sqrt_eigenvalues: Array1<f64>,
+  sqrt_eigenvalues: Array1<Complex<f64>>,
   m: Option<usize>,
 }
 
@@ -50,30 +51,29 @@ impl FgnFft {
       hurst,
       n,
       t: t.unwrap_or(1.0),
-      sqrt_eigenvalues: sqrt_eigenvalues.mapv(|x| x.re),
+      sqrt_eigenvalues,
       m,
     }
   }
 }
 
 impl Generator for FgnFft {
   fn sample(&self) -> Vec<f64> {
-    let mut rnd = Array1::<Complex<f64>>::zeros(2 * self.n);
-    rnd.par_mapv_inplace(|_| {
-      Complex::new(
-        rand::thread_rng().sample(StandardNormal),
-        rand::thread_rng().sample(StandardNormal),
-      )
-    });
-    let mut fgn = Array1::<Complex<f64>>::zeros(2 * self.n);
-    for (i, v) in rnd.iter().enumerate() {
-      fgn[i] = self.sqrt_eigenvalues[i] * v;
-    }
+    let rnd = Array1::<Complex<f64>>::random(
+      2 * self.n,
+      ComplexDistribution::new(StandardNormal, StandardNormal),
+    );
     let mut fgn_fft_handler = FftHandler::new(2 * self.n);
     let mut fgn_fft = Array1::<Complex<f64>>::zeros(2 * self.n);
-    ndfft_par(&fgn, &mut fgn_fft, &mut fgn_fft_handler, 0);
-    let mut fgn = fgn_fft.slice(s![1..self.n + 1]).mapv(|x| x.re);
-    fgn.par_mapv_inplace(|x| (x * (self.n as f64).powf(-self.hurst)) * self.t.powf(self.hurst));
+    ndfft_par(
+      &(&self.sqrt_eigenvalues * &rnd),
+      &mut fgn_fft,
+      &mut fgn_fft_handler,
+      0,
+    );
+    let fgn = fgn_fft
+      .slice(s![1..self.n + 1])
+      .mapv(|x| (x.re * (self.n as f64).powf(-self.hurst)) * self.t.powf(self.hurst));
     fgn.to_vec()
   }
 

src/processes/bm.rs

@@ -1,14 +1,10 @@
-use ndarray::Array1;
+use ndarray::{Array1, Axis};
 
 use crate::noises::gn;
 
 pub fn bm(n: usize, t: Option<f64>) -> Vec<f64> {
   let gn = gn::gn(n - 1, Some(t.unwrap_or(1.0)));
-  let mut bm = Array1::<f64>::zeros(n);
-
-  for i in 1..n {
-    bm[i] = bm[i - 1] + gn[i - 1];
-  }
-
-  bm.to_vec()
+  let mut bm = Array1::<f64>::from(gn);
+  bm.accumulate_axis_inplace(Axis(0), |&x, y| *y += x);
+  vec![0.0].into_iter().chain(bm.into_iter()).collect()
 }

src/processes/fbm.rs

@@ -2,12 +2,13 @@ use crate::{
   noises::fgn::{FgnCholesky, FgnFft},
   utils::{Generator, NoiseGenerationMethod},
 };
-use ndarray::Array1;
+use ndarray::{Array1, Axis};
 use rayon::prelude::*;
 
 pub struct Fbm {
   #[allow(dead_code)]
   hurst: f64,
+  #[allow(dead_code)]
   n: usize,
   m: Option<usize>,
   method: NoiseGenerationMethod,
@@ -54,14 +55,9 @@ impl Generator for Fbm {
       NoiseGenerationMethod::Fft => self.fgn_fft.as_ref().unwrap().sample(),
       NoiseGenerationMethod::Cholesky => self.fgn_cholesky.as_ref().unwrap().sample(),
     };
-
-    let mut fbm = Array1::<f64>::zeros(self.n);
-
-    for i in 1..self.n {
-      fbm[i] = fbm[i - 1] + fgn[i - 1];
-    }
-
-    fbm.to_vec()
+    let mut fbm = Array1::<f64>::from_vec(fgn);
+    fbm.accumulate_axis_inplace(Axis(0), |&x, y| *y += x);
+    vec![0.0].into_iter().chain(fbm.into_iter()).collect()
   }
 
   fn sample_par(&self) -> Vec<Vec<f64>> {