diff --git a/.gitignore b/.gitignore index 4fffb2f..fd3233f 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1,3 @@ /target /Cargo.lock +flamegraph.svg \ No newline at end of file diff --git a/Cargo.toml b/Cargo.toml index e296c2d..9d02490 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "stochastic-rs" -version = "0.6.5" +version = "0.7.0" edition = "2021" license = "MIT" description = "A Rust library for stochastic processes" @@ -20,16 +20,27 @@ num-complex = { version = "0.4.6", features = ["rand"] } rand = "0.8.5" rand_distr = "0.4.3" rayon = "1.10.0" -indicatif = "0.17.7" plotly = "0.9.0" ndarray-rand = "0.15.0" ndrustfft = "0.5.0" -derive_builder = "0.20.0" +derive_builder = "0.20.1" +tikv-jemallocator = { version = "0.6.0", optional = true } +mimalloc = { version = "0.1.43", optional = true } [dev-dependencies] +[features] +mimalloc = ["dep:mimalloc"] +jemalloc = ["dep:tikv-jemallocator"] +default = ["jemalloc"] + [lib] name = "stochastic_rs" crate-type = ["cdylib", "lib"] path = "src/lib.rs" +doctest = false + +[profile.release] +debug = false +codegen-units = 1 diff --git a/flamegraph.svg b/flamegraph.svg deleted file mode 100644 index 8923c83..0000000 --- a/flamegraph.svg +++ /dev/null @@ -1,491 +0,0 @@ -Flame Graph Reset ZoomSearch libdyld.dylib`dyld4::LibSystemHelpers::getenv (1 samples, 0.02%)libsystem_kernel.dylib`__exit (1 samples, 0.02%)libsystem_malloc.dylib`_malloc_zone_malloc (1 samples, 0.02%)libsystem_platform.dylib`_platform_memmove (390 samples, 6.32%)libsyste..stochastic-rs`<alloc::vec::Vec<T> as alloc::vec::spec_from_iter::SpecFromIter<T,I>>::from_iter (158 samples, 2.56%)st..libsystem_malloc.dylib`free_medium (2 samples, 0.03%)libsystem_kernel.dylib`madvise (1 samples, 0.02%)libsystem_malloc.dylib`free_medium (65 samples, 1.05%)libsystem_kernel.dylib`madvise (65 samples, 1.05%)libsystem_platform.dylib`__bzero (198 samples, 3.21%)lib..libsystem_platform.dylib`_platform_memmove (166 samples, 2.69%)li..stochastic-rs`<rand_distr::normal::StandardNormal as rand::distributions::distribution::Distribution<f64>>::sample (997 samples, 16.16%)stochastic-rs`<rand_distr..libsystem_malloc.dylib`szone_malloc_should_clear (1 samples, 0.02%)libsystem_malloc.dylib`medium_malloc_should_clear (1 samples, 0.02%)libsystem_malloc.dylib`medium_free_list_remove_ptr (1 samples, 0.02%)stochastic-rs`<rand_distr::normal::StandardNormal as rand::distributions::distribution::Distribution<f64>>::sample (1,530 samples, 24.81%)stochastic-rs`<rand_distr::normal::Stand..libsystem_m.dylib`exp (13 samples, 0.21%)stochastic-rs`<ndarray::ArrayBase<S,D> as ndarray_rand::RandomExt<S,A,D>>::random (2,687 samples, 43.56%)stochastic-rs`<ndarray::ArrayBase<S,D> as ndarray_rand::RandomExt<S,A,D..stochastic-rs`ndarray::iterators::to_vec_mapped (1,690 samples, 27.40%)stochastic-rs`ndarray::iterators::to_vec_map..stochastic-rs`DYLD-STUB$$exp (2 samples, 0.03%)stochastic-rs`ndarray::impl_methods::_<impl ndarray::ArrayBase<S,D>>::map (25 samples, 0.41%)stochastic-rs`ndarray::impl_constructors::_<impl ndarray::ArrayBase<S,D>>::build_uninit (1 samples, 0.02%)libsystem_malloc.dylib`szone_malloc_should_clear (1 samples, 0.02%)libsystem_malloc.dylib`medium_malloc_should_clear (1 samples, 0.02%)libsystem_malloc.dylib`medium_malloc_from_free_list (1 samples, 0.02%)stochastic-rs`ndarray::impl_ops::arithmetic_ops::_<impl core::ops::arith::Mul<&ndarray::ArrayBase<S2,E>> for &ndarray::ArrayBase<S,D>>::mul (349 samples, 5.66%)stochas..stochastic-rs`ndarray::zip::Zip<(P1,P2,PLast),D>::collect_with_partial (348 samples, 5.64%)stochas..libdyld.dylib`tlv_get_addr (1 samples, 0.02%)stochastic-rs`rayon::iter::plumbing::bridge_unindexed (2 samples, 0.03%)stochastic-rs`rayon_core::current_num_threads (1 samples, 0.02%)libsystem_malloc.dylib`medium_malloc_should_clear (1 samples, 0.02%)stochastic-rs`rustfft::array_utils::iter_chunks (159 samples, 2.58%)st..stochastic-rs`rustfft::algorithm::radix4::bitreversed_transpose (234 samples, 3.79%)stoc..stochastic-rs`<stochastic_rs::noises::fgn::FgnFft as stochastic_rs::utils::Generator>::sample (5,264 samples, 85.34%)stochastic-rs`<stochastic_rs::noises::fgn::FgnFft as stochastic_rs::utils::Generator>::samplestochastic-rs`rustfft::Fft::process (1,771 samples, 28.71%)stochastic-rs`rustfft::Fft::processstochastic-rs`rustfft::neon::neon_radix4::Neon64Radix4<T>::perform_fft_out_of_place (1,611 samples, 26.12%)stochastic-rs`rustfft::neon::neon_radix4::..stochastic-rs`rustfft::algorithm::radix4::reverse_bits (99 samples, 1.61%)stochastic-rs`DYLD-STUB$$free (1 samples, 0.02%)dyld`start (6,151 samples, 99.72%)dyld`startstochastic-rs`main (6,150 samples, 99.71%)stochastic-rs`mainstochastic-rs`std::rt::lang_start_internal (6,150 samples, 99.71%)stochastic-rs`std::rt::lang_start_internalstochastic-rs`std::rt::lang_start::_{{closure}} (6,150 samples, 99.71%)stochastic-rs`std::rt::lang_start::_{{closure}}stochastic-rs`std::sys::backtrace::__rust_begin_short_backtrace (6,150 samples, 99.71%)stochastic-rs`std::sys::backtrace::__rust_begin_short_backtracestochastic-rs`stochastic_rs::main (6,150 samples, 99.71%)stochastic-rs`stochastic_rs::mainstochastic-rs`<stochastic_rs::processes::fbm::Fbm as stochastic_rs::utils::Generator>::sample (6,150 samples, 99.71%)stochastic-rs`<stochastic_rs::processes::fbm::Fbm as stochastic_rs::utils::Generator>::samplestochastic-rs`ndarray::impl_methods::_<impl ndarray::ArrayBase<S,D>>::accumulate_axis_inplace (336 samples, 5.45%)stochas..libsystem_kernel.dylib`swtch_pri (7 samples, 0.11%)stochastic-rs`<rayon_core::latch::LatchRef<L> as rayon_core::latch::Latch>::set (1 samples, 0.02%)libsystem_kernel.dylib`__psynch_cvbroad (1 samples, 0.02%)stochastic-rs`rayon::iter::plumbing::bridge_unindexed_producer_consumer (1 samples, 0.02%)stochastic-rs`rayon_core::join::join_context::_{{closure}} (1 samples, 0.02%)stochastic-rs`rayon_core::join::join_context::_{{closure}} (1 samples, 0.02%)stochastic-rs`rayon_core::join::join_context::_{{closure}} (1 samples, 0.02%)stochastic-rs`rayon_core::join::join_context::_{{closure}} (1 samples, 0.02%)stochastic-rs`rayon_core::registry::WorkerThread::wait_until_cold (1 samples, 0.02%)stochastic-rs`<rayon_core::job::StackJob<L,F,R> as rayon_core::job::Job>::execute (1 samples, 0.02%)stochastic-rs`rayon_core::join::join_context::_{{closure}} (1 samples, 0.02%)stochastic-rs`rayon_core::registry::WorkerThread::wait_until_cold (1 samples, 0.02%)stochastic-rs`<rayon_core::job::StackJob<L,F,R> as rayon_core::job::Job>::execute (1 samples, 0.02%)stochastic-rs`rayon_core::join::join_context::_{{closure}} (1 samples, 0.02%)stochastic-rs`rayon_core::registry::WorkerThread::wait_until_cold (1 samples, 0.02%)stochastic-rs`rayon_core::registry::WorkerThread::find_work (1 samples, 0.02%)stochastic-rs`<core::iter::adapters::chain::Chain<A,B> as core::iter::traits::iterator::Iterator>::try_fold (1 samples, 0.02%)stochastic-rs`crossbeam_deque::deque::Stealer<T>::steal (1 samples, 0.02%)stochastic-rs`rayon_core::job::StackJob<L,F,R>::run_inline (1 samples, 0.02%)stochastic-rs`rayon::iter::plumbing::bridge_unindexed_producer_consumer (1 samples, 0.02%)stochastic-rs`rayon_core::join::join_context::_{{closure}} (1 samples, 0.02%)stochastic-rs`rayon_core::registry::WorkerThread::wait_until_cold (1 samples, 0.02%)stochastic-rs`<rayon_core::job::StackJob<L,F,R> as rayon_core::job::Job>::execute (1 samples, 0.02%)stochastic-rs`rayon_core::join::join_context::_{{closure}} (1 samples, 0.02%)stochastic-rs`rayon_core::registry::WorkerThread::wait_until_cold (1 samples, 0.02%)stochastic-rs`<rayon_core::job::StackJob<L,F,R> as rayon_core::job::Job>::execute (1 samples, 0.02%)stochastic-rs`rayon_core::join::join_context::_{{closure}} (1 samples, 0.02%)stochastic-rs`rayon_core::registry::WorkerThread::wait_until_cold (1 samples, 0.02%)stochastic-rs`rayon_core::registry::WorkerThread::find_work (1 samples, 0.02%)stochastic-rs`<core::iter::adapters::chain::Chain<A,B> as core::iter::traits::iterator::Iterator>::try_fold (1 samples, 0.02%)stochastic-rs`crossbeam_deque::deque::Stealer<T>::steal (1 samples, 0.02%)stochastic-rs`rayon_core::job::StackJob<L,F,R>::run_inline (1 samples, 0.02%)stochastic-rs`rayon::iter::plumbing::bridge_unindexed_producer_consumer (1 samples, 0.02%)stochastic-rs`rayon_core::join::join_context::_{{closure}} (1 samples, 0.02%)stochastic-rs`rayon_core::job::StackJob<L,F,R>::run_inline (1 samples, 0.02%)stochastic-rs`rayon::iter::plumbing::bridge_unindexed_producer_consumer (1 samples, 0.02%)stochastic-rs`rayon_core::join::join_context::_{{closure}} (1 samples, 0.02%)stochastic-rs`rayon_core::registry::WorkerThread::wait_until_cold (1 samples, 0.02%)stochastic-rs`<rayon_core::job::StackJob<L,F,R> as rayon_core::job::Job>::execute (1 samples, 0.02%)stochastic-rs`rayon_core::join::join_context::_{{closure}} (1 samples, 0.02%)stochastic-rs`rayon_core::registry::WorkerThread::wait_until_cold (1 samples, 0.02%)stochastic-rs`<rayon_core::job::StackJob<L,F,R> as rayon_core::job::Job>::execute (1 samples, 0.02%)stochastic-rs`rayon_core::join::join_context::_{{closure}} (1 samples, 0.02%)stochastic-rs`rayon_core::join::join_context::_{{closure}} (1 samples, 0.02%)stochastic-rs`rayon_core::registry::WorkerThread::wait_until_cold (1 samples, 0.02%)stochastic-rs`rayon_core::registry::WorkerThread::find_work (1 samples, 0.02%)stochastic-rs`<core::iter::adapters::chain::Chain<A,B> as core::iter::traits::iterator::Iterator>::try_fold (1 samples, 0.02%)stochastic-rs`crossbeam_deque::deque::Stealer<T>::steal (1 samples, 0.02%)stochastic-rs`rayon_core::job::StackJob<L,F,R>::run_inline (2 samples, 0.03%)stochastic-rs`rayon::iter::plumbing::bridge_unindexed_producer_consumer (2 samples, 0.03%)stochastic-rs`rayon_core::join::join_context::_{{closure}} (2 samples, 0.03%)stochastic-rs`rayon_core::registry::WorkerThread::wait_until_cold (1 samples, 0.02%)stochastic-rs`rayon_core::join::join_context::_{{closure}} (3 samples, 0.05%)stochastic-rs`rayon_core::join::join_context::_{{closure}} (1 samples, 0.02%)stochastic-rs`rayon_core::registry::WorkerThread::wait_until_cold (1 samples, 0.02%)stochastic-rs`<rayon_core::job::StackJob<L,F,R> as rayon_core::job::Job>::execute (1 samples, 0.02%)stochastic-rs`rayon_core::join::join_context::_{{closure}} (1 samples, 0.02%)stochastic-rs`rayon_core::join::join_context::_{{closure}} (1 samples, 0.02%)stochastic-rs`rayon_core::registry::WorkerThread::wait_until_cold (1 samples, 0.02%)stochastic-rs`<rayon_core::job::StackJob<L,F,R> as rayon_core::job::Job>::execute (1 samples, 0.02%)stochastic-rs`rayon_core::join::join_context::_{{closure}} (1 samples, 0.02%)stochastic-rs`<rayon_core::job::StackJob<L,F,R> as rayon_core::job::Job>::execute (1 samples, 0.02%)stochastic-rs`rayon_core::join::join_context::_{{closure}} (1 samples, 0.02%)stochastic-rs`rayon_core::registry::WorkerThread::wait_until_cold (1 samples, 0.02%)libsystem_kernel.dylib`swtch_pri (1 samples, 0.02%)stochastic-rs`rayon_core::join::join_context::_{{closure}} (2 samples, 0.03%)stochastic-rs`rayon_core::registry::WorkerThread::wait_until_cold (2 samples, 0.03%)stochastic-rs`rayon_core::registry::WorkerThread::find_work (1 samples, 0.02%)stochastic-rs`<core::iter::adapters::chain::Chain<A,B> as core::iter::traits::iterator::Iterator>::try_fold (1 samples, 0.02%)stochastic-rs`crossbeam_deque::deque::Stealer<T>::steal (1 samples, 0.02%)stochastic-rs`crossbeam_epoch::default::with_handle (1 samples, 0.02%)stochastic-rs`crossbeam_epoch::internal::Global::collect (1 samples, 0.02%)stochastic-rs`crossbeam_epoch::internal::Global::try_advance (1 samples, 0.02%)stochastic-rs`<rayon_core::job::StackJob<L,F,R> as rayon_core::job::Job>::execute (9 samples, 0.15%)stochastic-rs`rayon_core::join::join_context::_{{closure}} (8 samples, 0.13%)stochastic-rs`rayon_core::registry::WorkerThread::wait_until_cold (3 samples, 0.05%)stochastic-rs`<rayon_core::job::StackJob<L,F,R> as rayon_core::job::Job>::execute (3 samples, 0.05%)stochastic-rs`rayon_core::join::join_context::_{{closure}} (3 samples, 0.05%)stochastic-rs`rayon_core::registry::WorkerThread::wait_until_cold (1 samples, 0.02%)stochastic-rs`<rayon_core::job::StackJob<L,F,R> as rayon_core::job::Job>::execute (1 samples, 0.02%)stochastic-rs`rayon_core::join::join_context::_{{closure}} (1 samples, 0.02%)stochastic-rs`rayon_core::join::join_context::_{{closure}} (1 samples, 0.02%)stochastic-rs`rayon_core::job::StackJob<L,F,R>::run_inline (1 samples, 0.02%)stochastic-rs`rayon::iter::plumbing::bridge_unindexed_producer_consumer (1 samples, 0.02%)stochastic-rs`rayon_core::join::join_context::_{{closure}} (1 samples, 0.02%)stochastic-rs`rayon_core::registry::WorkerThread::wait_until_cold (1 samples, 0.02%)stochastic-rs`<rayon_core::job::StackJob<L,F,R> as rayon_core::job::Job>::execute (1 samples, 0.02%)stochastic-rs`rayon_core::join::join_context::_{{closure}} (1 samples, 0.02%)stochastic-rs`rayon_core::join::join_context::_{{closure}} (1 samples, 0.02%)stochastic-rs`rayon_core::registry::WorkerThread::wait_until_cold (1 samples, 0.02%)stochastic-rs`<rayon_core::job::StackJob<L,F,R> as rayon_core::job::Job>::execute (1 samples, 0.02%)stochastic-rs`rayon_core::join::join_context::_{{closure}} (1 samples, 0.02%)stochastic-rs`rayon_core::registry::WorkerThread::wait_until_cold (1 samples, 0.02%)libsystem_kernel.dylib`swtch_pri (1 samples, 0.02%)all (6,168 samples, 100%)libsystem_pthread.dylib`thread_start (17 samples, 0.28%)libsystem_pthread.dylib`_pthread_start (17 samples, 0.28%)stochastic-rs`std::sys::pal::unix::thread::Thread::new::thread_start (17 samples, 0.28%)stochastic-rs`core::ops::function::FnOnce::call_once{{vtable.shim}} (17 samples, 0.28%)stochastic-rs`std::sys::backtrace::__rust_begin_short_backtrace (17 samples, 0.28%)stochastic-rs`rayon_core::registry::ThreadBuilder::run (17 samples, 0.28%)stochastic-rs`rayon_core::registry::WorkerThread::wait_until_cold (17 samples, 0.28%)stochastic-rs`rayon_core::sleep::Sleep::sleep (1 samples, 0.02%)stochastic-rs`std::sync::condvar::Condvar::wait (1 samples, 0.02%)libsystem_kernel.dylib`__psynch_cvwait (1 samples, 0.02%) \ No newline at end of file diff --git a/src/lib.rs b/src/lib.rs index 3bdfd1c..441e27f 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -53,6 +53,16 @@ //! - Developed by [dancixx](https://github.com/dancixx). //! - Contributions and feedback are welcome! +#![feature(portable_simd)] + +#[cfg(feature = "mimalloc")] +#[global_allocator] +static GLOBAL: mimalloc::MiMalloc = mimalloc::MiMalloc; + +#[cfg(feature = "jemalloc")] +#[global_allocator] +static GLOBAL: tikv_jemallocator::Jemalloc = tikv_jemallocator::Jemalloc; + pub mod diffusions; pub mod interest; pub mod jumps; diff --git a/src/main.rs b/src/main.rs index 5249aab..16dbb49 100644 --- a/src/main.rs +++ b/src/main.rs @@ -1,12 +1,24 @@ -use std::time::Instant; - use stochastic_rs::{processes::fbm::Fbm, utils::Generator}; fn main() { - let start = Instant::now(); - let fbm = Fbm::new(0.95, 10000, Some(1.0), None); - for _ in 0..10000 { - fbm.sample(); + let fbm = Fbm::new(0.75, 10000, Some(1.0), None); + let mut runs = Vec::new(); + + for _ in 0..20 { + let start = std::time::Instant::now(); + for _ in 0..1000 { + let _ = fbm.sample(); + } + + let duration = start.elapsed(); + println!( + "Time elapsed in expensive_function() is: {:?}", + duration.as_secs_f32() + ); + runs.push(duration.as_secs_f32()); } - println!("{}", start.elapsed().as_secs_f64()); + + let sum: f32 = runs.iter().sum(); + let average = sum / runs.len() as f32; + println!("Average time: {}", average); } diff --git a/src/noises/fgn.rs b/src/noises/fgn.rs index d816cde..306cbd3 100644 --- a/src/noises/fgn.rs +++ b/src/noises/fgn.rs @@ -3,12 +3,14 @@ //! The `FgnFft` struct provides methods to generate fractional Gaussian noise (FGN) //! using the Fast Fourier Transform (FFT) approach. +use std::sync::Arc; + use crate::utils::Generator; use ndarray::parallel::prelude::*; use ndarray::{concatenate, prelude::*}; use ndarray_rand::rand_distr::StandardNormal; use ndarray_rand::RandomExt; -use ndrustfft::{ndfft_par, FftHandler}; +use ndrustfft::{ndfft, FftHandler}; use num_complex::{Complex, ComplexDistribution}; /// Struct for generating Fractional Gaussian Noise (FGN) using FFT. @@ -17,10 +19,9 @@ pub struct FgnFft { n: usize, offset: usize, t: f64, - sqrt_eigenvalues: Array1>, + sqrt_eigenvalues: Arc>>, m: Option, - fft_handler: FftHandler, - fft_fgn: Array1>, + fft_handler: Arc>, } impl FgnFft { @@ -54,7 +55,7 @@ impl FgnFft { let offset = n_ - n; let n = n_; let mut r = Array1::linspace(0.0, n as f64, n + 1); - r.par_mapv_inplace(|x| { + r.mapv_inplace(|x| { if x == 0.0 { 1.0 } else { @@ -71,18 +72,17 @@ impl FgnFft { let data = r.mapv(|v| Complex::new(v, 0.0)); let r_fft = FftHandler::new(r.len()); let mut sqrt_eigenvalues = Array1::>::zeros(r.len()); - ndfft_par(&data, &mut sqrt_eigenvalues, &r_fft, 0); - sqrt_eigenvalues.par_mapv_inplace(|x| Complex::new((x.re / (2.0 * n as f64)).sqrt(), x.im)); + ndfft(&data, &mut sqrt_eigenvalues, &r_fft, 0); + sqrt_eigenvalues.mapv_inplace(|x| Complex::new((x.re / (2.0 * n as f64)).sqrt(), x.im)); Self { hurst, n, offset, t: t.unwrap_or(1.0), - sqrt_eigenvalues, + sqrt_eigenvalues: Arc::new(sqrt_eigenvalues), m, - fft_handler: FftHandler::new(2 * n), - fft_fgn: Array1::>::zeros(2 * n), + fft_handler: Arc::new(FftHandler::new(2 * n)), } } } @@ -105,14 +105,13 @@ impl Generator for FgnFft { 2 * self.n, ComplexDistribution::new(StandardNormal, StandardNormal), ); - let fgn = &self.sqrt_eigenvalues * &rnd; - let fft_handler = self.fft_handler.clone(); - let mut fgn_fft = self.fft_fgn.clone(); - ndfft_par(&fgn, &mut fgn_fft, &fft_handler, 0); - let scale_factor = self.t.powf(self.hurst) * (self.n as f64).powf(-self.hurst); + let fgn = &*self.sqrt_eigenvalues * &rnd; + let mut fgn_fft = Array1::>::zeros(2 * self.n); + ndfft(&fgn, &mut fgn_fft, &*self.fft_handler, 0); + let scale = (self.n as f64).powf(-self.hurst) * self.t.powf(self.hurst); let fgn = fgn_fft .slice(s![1..self.n - self.offset + 1]) - .mapv(|x: Complex| x.re * scale_factor); + .mapv(|x: Complex| x.re * scale); fgn } @@ -137,12 +136,38 @@ impl Generator for FgnFft { panic!("m must be specified for parallel sampling"); } - let mut xs = Array2::zeros((self.m.unwrap(), self.n)); + let mut xs = Array2::zeros((self.m.unwrap(), self.n - self.offset)); xs.axis_iter_mut(Axis(0)).into_par_iter().for_each(|mut x| { - x.assign(&self.sample()); + x.assign(&self.sample().slice(s![..self.n - self.offset])); }); xs } } + +#[cfg(test)] +mod tests { + use plotly::{common::Line, Plot, Scatter}; + + use super::*; + + #[test] + fn plot() { + let fgn = FgnFft::new(0.9, 1000, Some(1.0), Some(1)); + let mut plot = Plot::new(); + let d = fgn.sample_par(); + for data in d.axis_iter(Axis(0)) { + let trace = Scatter::new((0..data.len()).collect::>(), data.to_vec()) + .mode(plotly::common::Mode::Lines) + .line( + Line::new() + .color("orange") + .shape(plotly::common::LineShape::Linear), + ) + .name("Fgn"); + plot.add_trace(trace); + } + plot.show(); + } +} diff --git a/src/processes/fbm.rs b/src/processes/fbm.rs index b470d24..690a23c 100644 --- a/src/processes/fbm.rs +++ b/src/processes/fbm.rs @@ -1,7 +1,7 @@ //! Fractional Brownian Motion (fBM) generator. use crate::{noises::fgn::FgnFft, utils::Generator}; -use ndarray::{Array1, Array2, Axis}; +use ndarray::{s, Array1, Array2, Axis}; use rayon::prelude::*; /// Struct for generating Fractional Brownian Motion (fBM). @@ -67,9 +67,15 @@ impl Generator for Fbm { /// ``` fn sample(&self) -> Array1 { let fgn = self.fgn.as_ref().unwrap().sample(); - let mut fbm = Array1::::from(fgn); - fbm.accumulate_axis_inplace(Axis(0), |&x, y| *y += x); - vec![0.0].into_iter().chain(fbm).collect() + + let mut fbm = Array1::::zeros(self.n); + fbm.slice_mut(s![1..]).assign(&fgn); + + for i in 1..self.n { + fbm[i] += fbm[i - 1]; + } + + fbm } /// Generates parallel samples of fractional Brownian motion (fBM). @@ -102,3 +108,29 @@ impl Generator for Fbm { xs } } + +#[cfg(test)] +mod tests { + use plotly::{common::Line, Plot, Scatter}; + + use super::*; + + #[test] + fn plot() { + let fbm = Fbm::new(0.45, 1000, Some(1.0), Some(10)); + let mut plot = Plot::new(); + let d = fbm.sample_par(); + for data in d.axis_iter(Axis(0)) { + let trace = Scatter::new((0..data.len()).collect::>(), data.to_vec()) + .mode(plotly::common::Mode::Lines) + .line( + Line::new() + .color("orange") + .shape(plotly::common::LineShape::Linear), + ) + .name("Fbm"); + plot.add_trace(trace); + } + plot.show(); + } +}