diffusionx 0.9.1

A multi-threaded crate for random number generation and stochastic process simulation, with optional CUDA GPU acceleration.
diffusionx-0.9.1 has been yanked.
Visit the last successful build: diffusionx-0.11.2

Implemented

Random Number Generation

  • Normal distribution
  • Uniform distribution
  • Exponential distribution
  • Poisson distribution
  • $\alpha$-stable distribution

GPU Acceleration (CUDA)

  • Brownian motion moment calculations
  • $\alpha$-stable Lévy process moment calculations
  • Ornstein-Uhlenbeck process moment calculations
  • $\alpha$-stable random number generation

[!NOTE] DiffusionX uses the high-quality Xoshiro256++ random number generator as the common entropy source across all distributions.

Stochastic Processes Simulation

  • Brownian motion
  • $\alpha$-stable Lévy process
  • Cauchy process
  • $\alpha$-stable subordinator
  • Inverse $\alpha$-stable subordinator
  • Poisson process
  • Fractional Brownian motion
  • Continuous-time random walk
  • Ornstein-Uhlenbeck process
  • Langevin equation
  • Generalized Langevin equation
  • Subordinated Langevin equation
  • Lévy walk
  • Birth-death process
  • Random walk
  • Brownian excursion
  • Brownian meander
  • Gamma process
  • Geometric Brownian motion
  • Brownian yet non-Gaussian process

Usage

Random Number Generation

use diffusionx::random::{normal, uniform, stable};
fn main() -> Result<(), Box<dyn std::error::Error>> {
    // Generate a normal random number with mean 0.0 and std 1.0
    let normal_sample = normal::rand(0.0, 1.0)?;
    // Generate 1000 standard normal random numbers
    let std_normal_samples = normal::standard_rands::<f64>(1000);

    // Generate a uniform random number in range [0, 10)
    let uniform_sample = uniform::range_rand(0..10)?;
    // Generate 1000 uniform random numbers in range [0, 1)
    let std_uniform_samples = uniform::standard_rands(1000);

    // Generate 1000 standard stable random numbers
    let stable_samples = stable::standard_rands(1.5, 0.5, 1000)?;

    Ok(())
}

Stochastic Process Simulation

use diffusionx::simulation::{prelude::*, continuous::Bm};
fn main() -> Result<(), Box<dyn std::error::Error>> {
    // Create standard Brownian motion object
    let bm = Bm::default();
    // Create trajectory with duration 1.0
    let traj = bm.duration(1.0)?;
    // Simulate Brownian motion trajectory with time step 0.01
    let (times, positions) = traj.simulate(0.01)?;
    println!("times: {:?}", times);
    println!("positions: {:?}", positions);

    // Calculate first-order raw moment with 1000 particles and time step 0.01
    let mean = traj.raw_moment(1, 1000, 0.01)?;
    println!("mean: {mean}");
    // Calculate second-order central moment with 1000 particles and time step 0.01
    let msd = traj.central_moment(2, 1000, 0.01)?;
    println!("MSD: {msd}");
    // Calculate EATAMSD with duration 100.0, delta 1.0, 10000 particles, time step 0.1,
    // and Gauss-Legendre quadrature order 10
    let eatamsd = bm.eatamsd(100.0, 1.0, 10000, 0.1, 10)?;
    println!("EATAMSD: {eatamsd}");
    // Calculate first passage time of Brownian motion with boundaries at -1.0 and 1.0
    let fpt = bm.fpt((-1.0, 1.0), 1000, 0.01)?;
    println!("fpt: {fpt}");
    Ok(())
}

Visualization

[!NOTE] The visualization requires the visualize feature to be enabled.

# In your Cargo.toml

[dependencies]

diffusionx = { version = "*", features = ["visualize"] }


use diffusionx::{
    simulation::{continuous::Bm, prelude::*},
};
fn main() -> Result<(), Box<dyn std::error::Error>> {
    // Create Brownian motion trajectory
    let bm = Bm::default();
    let traj = bm.duration(10.0)?;

    // Configure and create visualization
    let config = PlotConfigBuilder::default()
    .time_step(0.01)
    .output_path("brownian_motion.png")
    .caption("Brownian Motion Trajectory")
    .x_label("t")
    .y_label("B")
    .legend("bm")
    .size((800, 600))
    .backend(PlotterBackend::BitMap)
    .build()?;

    // Generate plot
    traj.plot(&config)?;

    Ok(())
}

Architecture and Extensibility

DiffusionX is designed with a trait-based system for high extensibility and performance:

Core Traits

  • ContinuousProcess: Base trait for continuous stochastic processes
  • PointProcess: Base trait for point processes
  • DiscreteProcess: Base trait for discrete stochastic processes
  • Moment: Trait for statistical moments calculation, including raw and central moments
  • Visualize: Trait for plotting process trajectories

Extending with Custom Processes

  1. Adding a New Continuous Process:

    #[derive(Debug, Clone)]
struct MyProcess {
    // Your parameters
    // Should be `Send + Sync` for parallel computation
    // and `Clone`
}

impl ContinuousProcess for MyProcess {
    fn start(&self) -> f64 {
        0.0  // or any default value
    }
    fn simulate(
         &self,
         duration: f64,
         time_step: f64
     ) -> XResult<(Vec<f64>, Vec<f64>)> {
        // Implement your simulation logic
        todo!()
    }
}

  2. Implementing ContinuousProcess trait automatically provides

    • mean mean
    • msd msd
    • raw moment raw_moment
    • central moment central_moment
    • first passage time fpt
    • occupation time occupation_time
    • TAMSD tamsd
    • visualization plot

Example:

Run the following Cargo command in your project directory:

cargo add diffusionx --features io,visualize

or add the following line to your Cargo.toml:

[dependencies]

diffusionx = { version = "*", features = ["io", "visualize"] }


#[cfg(feature = "io")]
use diffusionx::utils::write_csv;
use diffusionx::{
    XError, XResult, check_duration_time_step,
    random::normal,
    simulation::prelude::*,
    utils::{diff, linspace},
};

/// CIR
#[allow(clippy::upper_case_acronyms)]
#[derive(Clone)]
struct CIR {
    speed: f64,
    mean: f64,
    volatility: f64,
    start_position: f64,
}

impl CIR {
    fn new(
        speed: impl Into<f64>,
        mean: impl Into<f64>,
        volatility: impl Into<f64>,
        start_position: impl Into<f64>,
    ) -> XResult<Self> {
        let speed: f64 = speed.into();
        if speed <= 0.0 {
            return Err(XError::InvalidParameters(format!(
                "speed must be greater than 0, but got {speed}"
            )));
        }
        Ok(Self {
            speed,
            mean: mean.into(),
            volatility: volatility.into(),
            start_position: start_position.into(),
        })
    }
}

impl ContinuousProcess for CIR {
    fn start(&self) -> f64 {
        self.start_position
    }

    fn simulate(&self, duration: f64, time_step: f64) -> XResult<Pair> {
        check_duration_time_step(duration, time_step)?;

        let t = linspace(0.0, duration, time_step);
        let num_steps = t.len() - 1;
        let initial_x = self.start_position.max(0.0);
        let noises = normal::standard_rands::<f64>(num_steps);
        let delta = diff(&t);

        let x = std::iter::once(initial_x)
            .chain(
                noises
                    .iter()
                    .zip(delta)
                    .scan(initial_x, |state, (&xi, delta_t)| {
                        let current_x = *state;
                        let drift = self.speed * (self.mean - current_x);
                        let diffusion = self.volatility * current_x.sqrt().max(0.0);

                        let next_x = current_x + drift * delta_t + diffusion * xi * delta_t.sqrt();
                        *state = next_x.max(0.0);

                        Some(*state)
                    }),
            )
            .collect();

        Ok((t, x))
    }
}

fn main() -> XResult<()> {
    let duration = 10.0;
    let particles = 10_000;
    let time_step = 0.01;
    let cir = CIR::new(1, 1, 1, 0.5)?;

    #[allow(unused)]
    let (t, x) = cir.simulate(duration, time_step)?;
    #[cfg(feature = "io")]
    write_csv("tmp/CIR.csv", &t, &x)?;
    // mean
    let mean = cir.mean(duration, particles, time_step)?; // or let mean = traj.raw_moment(1, particles, time_step)?;
    println!("mean: {mean}");
    // msd
    let msd = cir.msd(duration, particles, time_step)?; // or let msd = traj.central_moment(2, particles, time_step)?;
    println!("MSD: {msd}");
    // FPT
    let max_duration = 1000.0;
    let fpt = cir
        .fpt((-1.0, 1.0), max_duration, time_step)?
        .unwrap_or(-1.0);
    println!("FPT: {fpt}");
    // occupation time
    let occupation_time = cir.occupation_time((-1.0, 1.0), duration, time_step)?;
    println!("Occupation Time: {occupation_time}");
    // TAMSD
    let slag = 1.0;
    let quad_order = 10;
    let tamsd = TAMSD::new(&cir, duration, slag)?;
    let eatamsd = tamsd.mean(particles, time_step, quad_order)?;
    println!("EATAMSD: {eatamsd}");

    #[cfg(feature = "visualize")]
    {
        let traj = cir.duration(duration)?;
        // Visualization
        let config = PlotConfigBuilder::default()
            .time_step(time_step)
            .output_path("tmp/CIR.svg")
            .caption("CIR")
            .show_grid(false)
            .x_label("t")
            .y_label("r")
            .legend("CIR")
            .backend(PlotterBackend::SVG)
            .build()
            .unwrap();
        traj.plot(&config)?;
    }
    Ok(())
}

Result:

mean: 0.9957644815350275
MSD: 0.7441251895881059
FPT: 0.38
Occupation Time: 4.719999999999995
EATAMSD: 0.6085042089895467

GPU Acceleration

[!NOTE] GPU acceleration requires the cuda feature and a CUDA-capable GPU.

# In your Cargo.toml

[dependencies]

diffusionx = { version = "*", features = ["cuda"] }


use diffusionx::{
    simulation::continuous::Bm,
    gpu::GPUMoment,
};
fn main() -> Result<(), Box<dyn std::error::Error>> {
    let bm = Bm::<f32>::default();

    // GPU-accelerated moment calculations
    let mean = bm.mean_gpu(1.0, 100_000, 0.01)?;
    let msd = bm.msd_gpu(1.0, 100_000, 0.01)?;
    let raw_moment = bm.raw_moment_gpu(1.0, 2, 100_000, 0.01)?;
    let central_moment = bm.central_moment_gpu(1.0, 2, 100_000, 0.01)?;

    // Fractional moments are also supported
    let frac_raw = bm.frac_raw_moment_gpu(1.0, 1.5, 100_000, 0.01)?;
    let frac_central = bm.frac_central_moment_gpu(1.0, 1.5, 100_000, 0.01)?;

    println!("Mean: {mean}, MSD: {msd}");
    Ok(())
}

GPU-accelerated stable random number generation:

use diffusionx::gpu::stable::standard_stable_rands;
fn main() -> Result<(), Box<dyn std::error::Error>> {
    // Generate 1 million stable random numbers on GPU
    let samples = standard_stable_rands(1.5, 0.5, 1_000_000)?;
    Ok(())
}

Benchmark

Performance benchmark tests compare the Rust, C++, Julia, and Python implementations, which can be found here.

License

Licensed under either of:

at your option.

Contribution

Unless you explicitly state otherwise, any contribution intentionally submitted for inclusion in the work by you, as defined in the Apache-2.0 license, shall be dual licensed as above, without any additional terms or conditions.

Dedicated to my brief yet unforgettable years in Lanzhou and to XX.