From 3420dc6bb62dfc057f87b7115631f09c793a3197 Mon Sep 17 00:00:00 2001
From: Daniel Boros <dancixx@gmail.com>
Date: Sat, 28 Oct 2023 10:33:08 +0200
Subject: [PATCH] feat: improve performance

---
 Cargo.toml           |  3 ++-
 src/main.rs          | 23 ++++++++++++++---------
 src/noises/fgn.rs    | 34 +++++++++++++++++-----------------
 src/processes/bm.rs  | 12 ++++--------
 src/processes/fbm.rs | 14 +++++---------
 5 files changed, 42 insertions(+), 44 deletions(-)

diff --git a/Cargo.toml b/Cargo.toml
index fe8de17..4c1e176 100644
--- a/Cargo.toml
+++ b/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]
 
diff --git a/src/main.rs b/src/main.rs
index d088f94..cb8ee27 100644
--- a/src/main.rs
+++ b/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();
diff --git a/src/noises/fgn.rs b/src/noises/fgn.rs
index 1cc8e2d..d01c8df 100644
--- a/src/noises/fgn.rs
+++ b/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,7 +51,7 @@ impl FgnFft {
       hurst,
       n,
       t: t.unwrap_or(1.0),
-      sqrt_eigenvalues: sqrt_eigenvalues.mapv(|x| x.re),
+      sqrt_eigenvalues,
       m,
     }
   }
@@ -58,22 +59,21 @@ impl FgnFft {
 
 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()
   }
 
diff --git a/src/processes/bm.rs b/src/processes/bm.rs
index 044928f..b53b48a 100644
--- a/src/processes/bm.rs
+++ b/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()
 }
diff --git a/src/processes/fbm.rs b/src/processes/fbm.rs
index 8ad04f2..7b2f4e1 100644
--- a/src/processes/fbm.rs
+++ b/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>> {