diff --git a/src/lib.rs b/src/lib.rs index 6c844bf..122046c 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,3 +1,5 @@ +// https://docs.rs/ndarray/latest/ndarray/doc/ndarray_for_numpy_users/index.html + pub mod c; pub mod diffusions; // add to c bindings diff --git a/src/main.rs b/src/main.rs index e889db5..802a6e3 100644 --- a/src/main.rs +++ b/src/main.rs @@ -5,31 +5,44 @@ use plotly::{Plot, Scatter}; use stochastic_rs::{prelude::*, processes::fbm::Fbm}; fn main() { - // let start = Instant::now(); - // let pb = ProgressBar::new(10000); - // for _ in 0..10000 { - // 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 _ = Fbm::new( + let fbm = Fbm::new( 0.7, - 5000, + 10000, None, Some(10000), Some(NoiseGenerationMethod::Fft( - FractionalNoiseGenerationMethod::DaviesHarte, + FractionalNoiseGenerationMethod::Kroese, )), - ) - .sample_par(); - println!("Time elapsed: {:?}", start.elapsed()); + ); + + let pb = ProgressBar::new(1); + // let mut plot = Plot::new(); + for _ in 0..1 { + let path = fbm.sample(); + // plot.add_trace(Scatter::new((0..5000).collect::>(), path)); + // plot.show(); + pb.inc(1); + } + pb.finish(); + println!("Time elapsed: {:?}", start.elapsed().as_secs_f64()); + + // let start = Instant::now(); + // let _ = par_fbm(10000, 0.7, 2500, None, None); + // println!("Time elapsed: {:?}", start.elapsed()); + + // let start = Instant::now(); + // let _ = Fbm::new( + // 0.7, + // 5000, + // None, + // Some(10000), + // Some(NoiseGenerationMethod::Fft( + // FractionalNoiseGenerationMethod::DaviesHarte, + // )), + // ) + // .sample_par(); + // println!("Time elapsed: {:?}", start.elapsed()); // CIR // let mut plot = Plot::new(); diff --git a/src/noises/fgn.rs b/src/noises/fgn.rs index aa71e38..27e4145 100644 --- a/src/noises/fgn.rs +++ b/src/noises/fgn.rs @@ -1,12 +1,15 @@ use crate::utils::{FractionalNoiseGenerationMethod, Generator}; use nalgebra::{ComplexField, DMatrix, DVector, Dim, Dyn, RowDVector}; use ndarray::{concatenate, prelude::*}; -use ndrustfft::{ndfft, ndifft, FftHandler}; +use ndrustfft::{ndfft, ndfft_par, ndifft, FftHandler}; use num_complex::Complex; use rand::{thread_rng, Rng}; use rand_distr::{Distribution, Normal, StandardNormal}; use rayon::prelude::*; -use std::cmp::Ordering::{Equal, Greater, Less}; +use std::{ + cmp::Ordering::{Equal, Greater, Less}, + sync::{Arc, Mutex}, +}; pub struct FgnFft { hurst: f64, @@ -31,39 +34,36 @@ impl FgnFft { match method.unwrap_or(FractionalNoiseGenerationMethod::Kroese) { FractionalNoiseGenerationMethod::Kroese => { - let r = Array::from_shape_fn((n + 1,), |i| { - if i == 0 { + let mut r = Array1::linspace(0.0, n as f64, n + 1); + r.par_mapv_inplace(|x| { + if x == 0.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)) + * ((x as f64 + 1.0).powf(2.0 * hurst) - 2.0 * (x as f64).powf(2.0 * hurst) + + (x 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 sqrt_eigenvalues = Array1::>::zeros(r.len()); - ndfft(&data, &mut sqrt_eigenvalues, &mut r_fft, 0); - let sqrt_eigenvalues = sqrt_eigenvalues - .mapv(|x| x.re / (2.0 * n as f64)) - .mapv(f64::sqrt); + ndfft_par(&data, &mut sqrt_eigenvalues, &mut r_fft, 0); + sqrt_eigenvalues.par_mapv_inplace(|x| Complex::new((x.re / (2.0 * n as f64)).sqrt(), x.im)); Self { hurst, n, t: t.unwrap_or(1.0), - sqrt_eigenvalues, + sqrt_eigenvalues: sqrt_eigenvalues.mapv(|x| x.re), m, method: method.unwrap_or(FractionalNoiseGenerationMethod::Kroese), } @@ -104,26 +104,25 @@ impl Generator for FgnFft { fn sample(&self) -> Vec { match self.method { FractionalNoiseGenerationMethod::Kroese => { - 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); - } - + rnd.par_mapv_inplace(|_| { + Complex::new( + rand::thread_rng().sample(StandardNormal), + rand::thread_rng().sample(StandardNormal), + ) + }); let mut fgn = Array1::>::zeros(2 * self.n); for (i, v) in rnd.iter().enumerate() { fgn[i] = self.sqrt_eigenvalues[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); + ndfft_par(&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() + let mut fgn: ArrayBase, ndarray::Dim<[usize; 1]>> = + 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)); + fgn.to_vec() } FractionalNoiseGenerationMethod::DaviesHarte => { // let m = 2usize.pow((self.n - 2).next_power_of_two() as u32) + 1;