feat: add improve fractional base process generation

src/jumps/vg.rs

@@ -26,18 +26,3 @@ pub fn vg(mu: f64, sigma: f64, nu: f64, n: usize, x0: Option<f64>, t: Option<f64
 
     vg.to_vec()
 }
-
-#[cfg(test)]
-mod tests {
-    use super::*;
-    use plotly::{Plot, Scatter};
-
-    #[test]
-    fn plot_vg() {
-        let mut plot = Plot::new();
-        let path = || vg(0.0, 2.0, 1.0, 50, None, None);
-
-        plot.add_trace(Scatter::new((0..50).collect::<Vec<usize>>(), path()));
-        plot.show();
-    }
-}

src/main.rs

@@ -1,20 +1,27 @@
 use std::time::Instant;
 
 use indicatif::ProgressBar;
-use stochastic_rs::prelude::*;
+use stochastic_rs::{
+    prelude::*,
+    processes::fbm::{fbm, Fbm},
+};
 
 fn main() {
     // let start = Instant::now();
     // let pb = ProgressBar::new(10000);
     // for _ in 0..10000 {
-    //     let _ = par_fbm(rayon::max_num_threads(), 0.7, 10000, None, None);
+    //     let _ = fbm(0.7, 2500, None, None);
     //     pb.inc(1);
     // }
     // pb.finish();
     // println!("Time elapsed: {:?}", start.elapsed());
 
+    // let start = Instant::now();
+    // let _ = par_fbm(10000, 0.7, 2500, None, None);
+    // println!("Time elapsed: {:?}", start.elapsed());
+
     let start = Instant::now();
-    let _ = par_fbm(100000, 0.7, 2500, None, None);
+    let _ = Fbm::new(0.7, 5000, None, Some(100000), None).sample_par();
     println!("Time elapsed: {:?}", start.elapsed());
 
     // CIR

src/noises/fgn_cholesky.rs

@@ -1,8 +1,62 @@
+use crate::utils::Generator;
 use nalgebra::{DMatrix, DVector, Dim, Dyn, RowDVector};
 use rand::{thread_rng, Rng};
 use rand_distr::StandardNormal;
+use rayon::prelude::*;
 use std::cmp::Ordering::{Equal, Greater, Less};
 
+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_sqrt = afc_matrix_sqrt(n, hurst);
+
+        Self {
+            hurst,
+            n,
+            t,
+            afc_sqrt,
+            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()
+    }
+}
+
 pub fn fgn(hurst: f64, n: usize, t: Option<f64>) -> Vec<f64> {
     if !(0.0..1.0).contains(&hurst) {
         panic!("Hurst parameter must be in (0, 1)")

src/noises/fgn_fft.rs

@@ -1,8 +1,94 @@
+use crate::utils::Generator;
 use ndarray::{concatenate, prelude::*};
 use ndrustfft::{ndfft, FftHandler};
 use num_complex::Complex;
 use rand::Rng;
 use rand_distr::StandardNormal;
+use rayon::prelude::*;
+
+pub struct FgnFft {
+    hurst: f64,
+    n: usize,
+    t: f64,
+    lambda: Array1<f64>,
+    m: Option<usize>,
+}
+
+impl FgnFft {
+    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 r = Array::from_shape_fn((n + 1,), |i| {
+            if i == 0 {
+                1.0
+            } else {
+                0.5 * ((i as f64 + 1.0).powf(2.0 * hurst) - 2.0 * (i as f64).powf(2.0 * hurst)
+                    + (i as f64 - 1.0).powf(2.0 * hurst))
+            }
+        });
+
+        let r = concatenate(
+            Axis(0),
+            #[allow(clippy::reversed_empty_ranges)]
+            &[r.view(), r.slice(s![..;-1]).slice(s![1..-1]).view()],
+        )
+        .unwrap();
+
+        let mut data = Array1::<Complex<f64>>::zeros(r.len());
+        for (i, v) in r.iter().enumerate() {
+            data[i] = Complex::new(*v, 0.0);
+        }
+        let mut r_fft = FftHandler::new(r.len());
+        let mut lambda = Array1::<Complex<f64>>::zeros(r.len());
+        ndfft(&data, &mut lambda, &mut r_fft, 0);
+        let lambda = lambda.mapv(|x| x.re / (2.0 * n as f64)).mapv(f64::sqrt);
+
+        Self {
+            hurst,
+            n,
+            t: t.unwrap_or(1.0),
+            lambda,
+            m,
+        }
+    }
+}
+
+impl Generator for FgnFft {
+    fn sample(&self) -> Vec<f64> {
+        let mut rng = rand::thread_rng();
+        let mut rnd = Array1::<Complex<f64>>::zeros(2 * self.n);
+        for (_, v) in rnd.iter_mut().enumerate() {
+            let real: f64 = rng.sample(StandardNormal);
+            let imag: f64 = rng.sample(StandardNormal);
+            *v = Complex::new(real, imag);
+        }
+
+        let mut fgn = Array1::<Complex<f64>>::zeros(2 * self.n);
+        for (i, v) in rnd.iter().enumerate() {
+            fgn[i] = self.lambda[i] * v;
+        }
+        let mut fgn_fft_handler = FftHandler::new(2 * self.n);
+        let mut fgn_fft = Array1::<Complex<f64>>::zeros(2 * self.n);
+        ndfft(&fgn, &mut fgn_fft, &mut fgn_fft_handler, 0);
+
+        let fgn = fgn_fft.slice(s![1..self.n + 1]).mapv(|x| x.re);
+        fgn.mapv(|x| (x * (self.n as f64).powf(-self.hurst)) * self.t.powf(self.hurst))
+            .to_vec()
+    }
+
+    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()
+    }
+}
 
 // TODO: improve performance
 pub fn fgn(hurst: f64, n: usize, t: Option<f64>) -> Vec<f64> {

src/processes/fbm.rs

@@ -1,8 +1,79 @@
 use crate::{
-    noises::{fgn_cholesky, fgn_fft},
-    utils::NoiseGenerationMethod,
+    noises::{
+        fgn_cholesky::{self, FgnCholesky},
+        fgn_fft::{self, FgnFft},
+    },
+    utils::{Generator, NoiseGenerationMethod},
 };
 use ndarray::Array1;
+use rayon::prelude::*;
+
+pub struct Fbm {
+    n: usize,
+    m: Option<usize>,
+    method: NoiseGenerationMethod,
+    fgn_fft: Option<FgnFft>,
+    fgn_cholesky: Option<FgnCholesky>,
+}
+
+impl Fbm {
+    pub fn new(
+        hurst: f64,
+        n: usize,
+        t: Option<f64>,
+        m: Option<usize>,
+        method: Option<NoiseGenerationMethod>,
+    ) -> 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 {
+                n,
+                m,
+                method: NoiseGenerationMethod::Fft,
+                fgn_fft: Some(FgnFft::new(hurst, n - 1, t, None)),
+                fgn_cholesky: None,
+            },
+            NoiseGenerationMethod::Cholesky => Self {
+                n,
+                m,
+                method: NoiseGenerationMethod::Cholesky,
+                fgn_fft: None,
+                fgn_cholesky: Some(FgnCholesky::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 mut fbm = Array1::<f64>::zeros(self.n);
+
+        for i in 1..self.n {
+            fbm[i] = fbm[i - 1] + fgn[i - 1];
+        }
+
+        fbm.to_vec()
+    }
+
+    fn sample_par(&self) -> Vec<Vec<f64>> {
+        if self.m.is_none() {
+            panic!("Number of paths must be specified")
+        }
+
+        (0..self.m.unwrap())
+            .into_par_iter()
+            .map(|_| self.sample())
+            .collect()
+    }
+}
 
 pub fn fbm(
     hurst: f64,
@@ -26,29 +97,3 @@ pub fn fbm(
 
     fbm.to_vec()
 }
-
-#[cfg(test)]
-mod tests {
-    use super::*;
-    use crate::statistics::fractal_dim::higuchi_fd;
-    use crate::utils::NoiseGenerationMethod;
-
-    #[test]
-    fn test_fbm() {
-        let (hurst, n) = (0.7, 10000);
-
-        let path = fbm(hurst, n, None, Some(NoiseGenerationMethod::Fft));
-        assert_eq!(path.len(), n);
-        assert_eq!(path[0], 0.0);
-
-        let h = higuchi_fd(&path, 10);
-        println!("Higuchi FD: {}", h);
-        assert!((2.0 - h) < 10e-1);
-
-        let path = fbm(hurst, n, None, Some(NoiseGenerationMethod::Cholesky));
-        assert_eq!(path.len(), n);
-
-        let h = higuchi_fd(&path, 10);
-        assert!((2.0 - h) < 10e-1);
-    }
-}

src/utils.rs

@@ -4,3 +4,8 @@ pub enum NoiseGenerationMethod {
     Cholesky,
     Fft,
 }
+
+pub trait Generator: Sync + Send {
+    fn sample(&self) -> Vec<f64>;
+    fn sample_par(&self) -> Vec<Vec<f64>>;
+}