diff --git a/src/jumps/vg.rs b/src/jumps/vg.rs index e52c294..d851278 100644 --- a/src/jumps/vg.rs +++ b/src/jumps/vg.rs @@ -26,18 +26,3 @@ pub fn vg(mu: f64, sigma: f64, nu: f64, n: usize, x0: Option, t: Option>(), path())); - plot.show(); - } -} diff --git a/src/main.rs b/src/main.rs index 043fd49..b204398 100644 --- a/src/main.rs +++ b/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 diff --git a/src/noises/fgn_cholesky.rs b/src/noises/fgn_cholesky.rs index 6c72711..9ff1867 100644 --- a/src/noises/fgn_cholesky.rs +++ b/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, + afc_sqrt: DMatrix, + m: Option, +} + +impl FgnCholesky { + pub fn new(hurst: f64, n: usize, t: Option, m: Option) -> 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 { + let noise = thread_rng() + .sample_iter::(StandardNormal) + .take(self.n) + .collect(); + let noise = DVector::::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> { + 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) -> Vec { if !(0.0..1.0).contains(&hurst) { panic!("Hurst parameter must be in (0, 1)") diff --git a/src/noises/fgn_fft.rs b/src/noises/fgn_fft.rs index e28007c..3e75820 100644 --- a/src/noises/fgn_fft.rs +++ b/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, + m: Option, +} + +impl FgnFft { + pub fn new(hurst: f64, n: usize, t: Option, m: Option) -> 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::>::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::>::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 { + let mut rng = rand::thread_rng(); + let mut rnd = Array1::>::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::>::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::>::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> { + 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) -> Vec { diff --git a/src/processes/fbm.rs b/src/processes/fbm.rs index c9d2726..8820dc7 100644 --- a/src/processes/fbm.rs +++ b/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, + method: NoiseGenerationMethod, + fgn_fft: Option, + fgn_cholesky: Option, +} + +impl Fbm { + pub fn new( + hurst: f64, + n: usize, + t: Option, + m: Option, + method: Option, + ) -> 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 { + 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::::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> { + 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); - } -} diff --git a/src/utils.rs b/src/utils.rs index 664227d..c043152 100644 --- a/src/utils.rs +++ b/src/utils.rs @@ -4,3 +4,8 @@ pub enum NoiseGenerationMethod { Cholesky, Fft, } + +pub trait Generator: Sync + Send { + fn sample(&self) -> Vec; + fn sample_par(&self) -> Vec>; +}