refactor: remove nalgebra and cholesky method

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "stochastic-rs"
-version = "0.4.0"
+version = "0.4.1"
 edition = "2021"
 license = "MIT"
 description = "A Rust library for stochastic processes"
@@ -12,7 +12,6 @@ keywords = ["stochastic", "process", "random", "simulation", "monte-carlo"]
 
 [dependencies]
 linreg = "0.2.0"
-nalgebra = { version = "0.32.2", features = ["rayon"] }
 ndarray = { version = "0.15.6", features = ["rayon", "matrixmultiply-threading"] }
 num-complex = { version = "0.4.4", features = ["rand"] }
 rand = "0.8.5"

src/jumps/vg.rs

@@ -1,7 +1,7 @@
 use crate::noises::gn;
 use ndarray::Array1;
+use ndarray_rand::rand_distr::Gamma;
 use ndarray_rand::RandomExt;
-use rand_distr::Gamma;
 
 pub fn vg(mu: f64, sigma: f64, nu: f64, n: usize, x0: Option<f64>, t: Option<f64>) -> Vec<f64> {
   let dt = t.unwrap_or(1.0) / n as f64;

src/noises/fgn.rs

@@ -1,13 +1,10 @@
 use crate::utils::Generator;
-use nalgebra::{DMatrix, DVector, Dim, Dyn, RowDVector};
 use ndarray::parallel::prelude::*;
 use ndarray::{concatenate, prelude::*};
+use ndarray_rand::rand_distr::StandardNormal;
 use ndarray_rand::RandomExt;
 use ndrustfft::{ndfft_par, FftHandler};
 use num_complex::{Complex, ComplexDistribution};
-use rand::{thread_rng, Rng};
-use rand_distr::StandardNormal;
-use std::cmp::Ordering::{Equal, Greater, Less};
 
 pub struct FgnFft {
   hurst: f64,
@@ -87,81 +84,3 @@ impl Generator for FgnFft {
       .collect()
   }
 }
-
-pub struct FgnCholesky {
-  hurst: f64,
-  n: usize,
-  t: Option<f64>,
-  afc_sqrt: DMatrix<f64>,
-  m: Option<usize>,
-}
-
-impl FgnCholesky {
-  pub fn new(hurst: f64, n: usize, t: Option<f64>, m: Option<usize>) -> Self {
-    if !(0.0..1.0).contains(&hurst) {
-      panic!("Hurst parameter must be between 0 and 1");
-    }
-
-    let afc_matrix_sqrt = |n: usize, hurst: f64| -> DMatrix<f64> {
-      let mut acf_v = RowDVector::<f64>::zeros(n);
-      acf_v[0] = 1.0;
-
-      for i in 1..n {
-        let idx = i as f64;
-
-        acf_v[i] = 0.5
-          * ((idx + 1.0).powf(2.0 * hurst) - 2.0 * idx.powf(2.0 * hurst)
-            + (idx - 1.0).powf(2.0 * hurst))
-      }
-      let mut m: nalgebra::Matrix<f64, Dyn, Dyn, nalgebra::VecStorage<f64, Dyn, Dyn>> =
-        DMatrix::<f64>::zeros_generic(Dyn::from_usize(n), Dyn::from_usize(n));
-
-      for i in 0..n {
-        for j in 0..n {
-          match i.cmp(&j) {
-            Equal => m[(i, j)] = acf_v[0],
-            Greater => m[(i, j)] = acf_v[i - j],
-            Less => continue,
-          }
-        }
-      }
-
-      m.cholesky().unwrap().l()
-    };
-
-    Self {
-      hurst,
-      n,
-      t,
-      afc_sqrt: afc_matrix_sqrt(n, hurst),
-      m,
-    }
-  }
-}
-
-impl Generator for FgnCholesky {
-  fn sample(&self) -> Vec<f64> {
-    let noise = thread_rng()
-      .sample_iter::<f64, StandardNormal>(StandardNormal)
-      .take(self.n)
-      .collect();
-    let noise = DVector::<f64>::from_vec(noise);
-
-    ((self.afc_sqrt.clone() * noise).transpose()
-      * (self.n as f64).powf(-self.hurst)
-      * self.t.unwrap_or(1.0).powf(self.hurst))
-    .data
-    .as_vec()
-    .clone()
-  }
-
-  fn sample_par(&self) -> Vec<Vec<f64>> {
-    if self.m.is_none() {
-      panic!("m must be specified for parallel sampling")
-    }
-    (0..self.m.unwrap())
-      .into_par_iter()
-      .map(|_| self.sample())
-      .collect()
-  }
-}

src/noises/gn.rs

@@ -1,11 +1,23 @@
-use rand::{thread_rng, Rng};
-use rand_distr::StandardNormal;
+use ndarray::Array1;
+use ndarray_rand::rand_distr::Normal;
+use ndarray_rand::RandomExt;
 
 pub fn gn(n: usize, t: Option<f64>) -> Vec<f64> {
   let sqrt_dt = (t.unwrap_or(1.0) / n as f64).sqrt();
-  thread_rng()
-    .sample_iter::<f64, StandardNormal>(StandardNormal)
-    .take(n)
-    .map(|z| z * sqrt_dt)
-    .collect()
+  let gn = Array1::random(n, Normal::new(0.0, sqrt_dt).unwrap());
+
+  gn.to_vec()
+}
+
+#[cfg(test)]
+mod tests {
+  use super::*;
+
+  #[test]
+  fn test_gn() {
+    let n = 10;
+    let t = 1.0;
+    let gn = gn(n, Some(t));
+    assert_eq!(gn.len(), n);
+  }
 }

src/processes/fbm.rs

@@ -1,7 +1,4 @@
-use crate::{
-  noises::fgn::{FgnCholesky, FgnFft},
-  utils::{Generator, NoiseGenerationMethod},
-};
+use crate::{noises::fgn::FgnFft, utils::Generator};
 use ndarray::{Array1, Axis};
 use rayon::prelude::*;
 
@@ -11,50 +8,27 @@ pub struct Fbm {
   #[allow(dead_code)]
   n: usize,
   m: Option<usize>,
-  method: NoiseGenerationMethod,
-  fgn_fft: Option<FgnFft>,
-  fgn_cholesky: Option<FgnCholesky>,
+  fgn: Option<FgnFft>,
 }
 
 impl Fbm {
-  pub fn new(
-    hurst: f64,
-    n: usize,
-    t: Option<f64>,
-    m: Option<usize>,
-    method: Option<NoiseGenerationMethod>,
-  ) -> Self {
+  pub fn new(hurst: f64, n: usize, t: Option<f64>, m: Option<usize>) -> Self {
     if !(0.0..1.0).contains(&hurst) {
       panic!("Hurst parameter must be in (0, 1)")
     }
 
-    match method.unwrap_or(NoiseGenerationMethod::Fft) {
-      NoiseGenerationMethod::Fft => Self {
-        hurst,
-        n,
-        m,
-        method: NoiseGenerationMethod::Fft,
-        fgn_fft: Some(FgnFft::new(hurst, n - 1, t, None)),
-        fgn_cholesky: None,
-      },
-      NoiseGenerationMethod::Cholesky => Self {
-        hurst,
-        n,
-        m,
-        method: NoiseGenerationMethod::Cholesky,
-        fgn_fft: None,
-        fgn_cholesky: Some(FgnCholesky::new(hurst, n - 1, t, None)),
-      },
+    Self {
+      hurst,
+      n,
+      m,
+      fgn: Some(FgnFft::new(hurst, n - 1, t, None)),
     }
   }
 }
 
 impl Generator for Fbm {
   fn sample(&self) -> Vec<f64> {
-    let fgn = match self.method {
-      NoiseGenerationMethod::Fft => self.fgn_fft.as_ref().unwrap().sample(),
-      NoiseGenerationMethod::Cholesky => self.fgn_cholesky.as_ref().unwrap().sample(),
-    };
+    let fgn = self.fgn.as_ref().unwrap().sample();
     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).collect()

src/processes/poisson.rs

@@ -1,8 +1,8 @@
 use ndarray::Array1;
+use ndarray_rand::rand_distr::{Distribution, Exp};
 use ndarray_rand::rand_distr::{Normal, Poisson};
 use ndarray_rand::RandomExt;
 use rand::thread_rng;
-use rand_distr::{Distribution, Exp};
 
 pub fn poisson(n: usize, lambda: usize, t_max: Option<f64>) -> Vec<f64> {
   if n == 0 || lambda == 0 {

src/utils.rs

@@ -1,12 +1,5 @@
 use std::error::Error;
 
-#[repr(C)]
-#[derive(Clone, Copy)]
-pub enum NoiseGenerationMethod {
-  Cholesky,
-  Fft,
-}
-
 pub trait Generator: Sync + Send {
   fn sample(&self) -> Vec<f64>;
   fn sample_par(&self) -> Vec<Vec<f64>>;