diff --git a/Cargo.toml b/Cargo.toml index 2f651cf..aa8ce6b 100644 --- a/Cargo.toml +++ b/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" diff --git a/src/jumps/vg.rs b/src/jumps/vg.rs index cd19c4b..7f27db8 100644 --- a/src/jumps/vg.rs +++ b/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, t: Option) -> Vec { let dt = t.unwrap_or(1.0) / n as f64; diff --git a/src/noises/fgn.rs b/src/noises/fgn.rs index 41335fd..41558ee 100644 --- a/src/noises/fgn.rs +++ b/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, - 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_matrix_sqrt = |n: usize, hurst: f64| -> DMatrix { - let mut acf_v = RowDVector::::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> = - DMatrix::::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 { - 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() - } -} diff --git a/src/noises/gn.rs b/src/noises/gn.rs index 156090d..c16d26f 100644 --- a/src/noises/gn.rs +++ b/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) -> Vec { let sqrt_dt = (t.unwrap_or(1.0) / n as f64).sqrt(); - thread_rng() - .sample_iter::(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); + } } diff --git a/src/processes/fbm.rs b/src/processes/fbm.rs index 5eab2ae..79218ae 100644 --- a/src/processes/fbm.rs +++ b/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, - method: NoiseGenerationMethod, - fgn_fft: Option, - fgn_cholesky: Option, + fgn: Option, } impl Fbm { - pub fn new( - hurst: f64, - n: usize, - t: Option, - m: Option, - method: Option, - ) -> Self { + pub fn new(hurst: f64, n: usize, t: Option, m: 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 { - 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 { - 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::::from_vec(fgn); fbm.accumulate_axis_inplace(Axis(0), |&x, y| *y += x); vec![0.0].into_iter().chain(fbm).collect() diff --git a/src/processes/poisson.rs b/src/processes/poisson.rs index 16080d0..959d550 100644 --- a/src/processes/poisson.rs +++ b/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) -> Vec { if n == 0 || lambda == 0 { diff --git a/src/utils.rs b/src/utils.rs index 0ff3610..f1dc705 100644 --- a/src/utils.rs +++ b/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; fn sample_par(&self) -> Vec>;