feat!: use fastrand crate for random number generation

Cargo.toml

@@ -16,10 +16,10 @@ testing = []
 
 [dependencies]
 nalgebra = "0.31"
-rand = { version = "0.8", features = ["small_rng"] }
-rand_distr = "0.4"
 approx = "0.5"
 num-traits = "0.2"
 thiserror = "1"
 log = "0.4"
 getset = "0.1"
+fastrand = "2.0.1"
+fastrand-contrib = "0.1.0"

src/core/base.rs

@@ -1,3 +1,6 @@
+use fastrand::Rng;
+use fastrand_contrib::RngExt;
+
 use super::domain::Domain;
 
 /// Trait implemented by real numbers.
@@ -13,6 +16,16 @@ pub trait RealField: nalgebra::RealField {
     const EPSILON_CBRT: Self;
 }
 
+impl RealField for f32 {
+    const EPSILON_SQRT: Self = 0.00034526698;
+    const EPSILON_CBRT: Self = 0.0049215667;
+}
+
+impl RealField for f64 {
+    const EPSILON_SQRT: Self = 0.000000014901161193847656;
+    const EPSILON_CBRT: Self = 0.0000060554544523933395;
+}
+
 /// The base trait for [`System`](super::system::System) and
 /// [`Function`](super::function::Function).
 pub trait Problem {
@@ -24,12 +37,42 @@ pub trait Problem {
     fn domain(&self) -> Domain<Self::Field>;
 }
 
-impl RealField for f32 {
-    const EPSILON_SQRT: Self = 0.00034526698;
-    const EPSILON_CBRT: Self = 0.0049215667;
+/// Trait for types that can be sampled.
+pub trait Sample {
+    /// Sample value from the whole range of the type.
+    fn sample_any(rng: &mut Rng) -> Self;
+
+    /// Sample from uniform distribution (inclusive on both sides).
+    fn sample_uniform(lower: Self, upper: Self, rng: &mut Rng) -> Self;
+
+    /// Sample from normal distribution.
+    fn sample_normal(mu: Self, sigma: Self, rng: &mut Rng) -> Self;
 }
 
-impl RealField for f64 {
-    const EPSILON_SQRT: Self = 0.000000014901161193847656;
-    const EPSILON_CBRT: Self = 0.0000060554544523933395;
+impl Sample for f32 {
+    fn sample_any(rng: &mut Rng) -> Self {
+        rng.f32()
+    }
+
+    fn sample_uniform(lower: Self, upper: Self, rng: &mut Rng) -> Self {
+        rng.f32_range(lower..=upper)
+    }
+
+    fn sample_normal(mu: Self, sigma: Self, rng: &mut Rng) -> Self {
+        rng.f32_normal_approx(mu, sigma)
+    }
+}
+
+impl Sample for f64 {
+    fn sample_any(rng: &mut Rng) -> Self {
+        rng.f64()
+    }
+
+    fn sample_uniform(lower: Self, upper: Self, rng: &mut Rng) -> Self {
+        rng.f64_range(lower..=upper)
+    }
+
+    fn sample_normal(mu: Self, sigma: Self, rng: &mut Rng) -> Self {
+        rng.f64_normal_approx(mu, sigma)
+    }
 }

src/core/domain.rs

@@ -2,12 +2,12 @@
 
 use std::iter::FromIterator;
 
+use fastrand::Rng;
 use na::{Dim, DimName};
 use nalgebra as na;
 use nalgebra::{storage::StorageMut, OVector, RealField, Vector};
-use rand::Rng;
-use rand_distr::uniform::SampleUniform;
-use rand_distr::{Distribution, Standard, Uniform};
+
+use crate::core::Sample;
 
 fn estimate_magnitude<T: RealField + Copy>(lower: T, upper: T) -> T {
     let ten = T::from_subset(&10.0);
@@ -154,18 +154,17 @@ impl<T: RealField + Copy> Domain<T> {
     }
 
     /// Samples a point in the domain.
-    pub fn sample<D, Sx, R: Rng>(&self, x: &mut Vector<T, D, Sx>, rng: &mut R)
+    pub fn sample<D, Sx>(&self, x: &mut Vector<T, D, Sx>, rng: &mut Rng)
     where
         D: Dim,
         Sx: StorageMut<T, D> + na::IsContiguous,
-        T: SampleUniform,
-        Standard: Distribution<T>,
+        T: Sample,
     {
         x.iter_mut()
             .zip(self.lower.iter().copied().zip(self.upper.iter().copied()))
             .for_each(|(xi, (li, ui))| {
                 *xi = if !li.is_finite() || !ui.is_finite() {
-                    let random: T = rng.gen();
+                    let random = T::sample_any(rng);
 
                     if li.is_finite() || ui.is_finite() {
                         let clamped = random.max(li).min(ui);
@@ -175,7 +174,7 @@ impl<T: RealField + Copy> Domain<T> {
                         random
                     }
                 } else {
-                    Uniform::new_inclusive(li, ui).sample(rng)
+                    T::sample_uniform(li, ui, rng)
                 };
             });
     }

src/core/optimizer.rs

@@ -23,23 +23,21 @@ use super::{domain::Domain, function::Function};
 /// use gomez::nalgebra as na;
 /// use gomez::*;
 /// use na::{storage::StorageMut, Dynamic, IsContiguous, Vector};
-/// use rand::Rng;
-/// use rand_distr::{uniform::SampleUniform, Distribution, Standard};
+/// use fastrand::Rng;
 ///
-/// struct Random<R> {
-///     rng: R,
+/// struct Random {
+///     rng: Rng,
 /// }
 ///
-/// impl<R> Random<R> {
-///     fn new(rng: R) -> Self {
+/// impl Random {
+///     fn new(rng: Rng) -> Self {
 ///         Self { rng }
 ///     }
 /// }
 ///
-/// impl<F: Function, R: Rng> Optimizer<F> for Random<R>
+/// impl<F: Function> Optimizer<F> for Random
 /// where
-///     F::Field: SampleUniform,
-///     Standard: Distribution<F::Field>,
+///     F::Field: Sample,
 /// {
 ///     const NAME: &'static str = "Random";
 ///     type Error = std::convert::Infallible;

src/core/solver.rs

@@ -23,23 +23,21 @@ use super::{domain::Domain, system::System};
 /// use gomez::nalgebra as na;
 /// use gomez::*;
 /// use na::{storage::StorageMut, Dynamic, IsContiguous, Vector};
-/// use rand::Rng;
-/// use rand_distr::{uniform::SampleUniform, Distribution, Standard};
+/// use fastrand::Rng;
 ///
-/// struct Random<R> {
-///     rng: R,
+/// struct Random {
+///     rng: Rng,
 /// }
 ///
-/// impl<R> Random<R> {
-///     fn new(rng: R) -> Self {
+/// impl Random {
+///     fn new(rng: Rng) -> Self {
 ///         Self { rng }
 ///     }
 /// }
 ///
-/// impl<F: System, R: Rng> Solver<F> for Random<R>
+/// impl<F: System> Solver<F> for Random
 /// where
-///     F::Field: SampleUniform,
-///     Standard: Distribution<F::Field>,
+///     F::Field: Sample,
 /// {
 ///     const NAME: &'static str = "Random";
 ///     type Error = std::convert::Infallible;

src/solver/lipo.rs

@@ -9,18 +9,17 @@
 //! \[2\] [A Global Optimization Algorithm Worth
 //! Using](http://blog.dlib.net/2017/12/a-global-optimization-algorithm-worth.html)
 
+use fastrand::Rng;
 use getset::{CopyGetters, Setters};
 use log::{debug, trace};
 use nalgebra::{
     convert, try_convert, ComplexField, DimName, Dynamic, IsContiguous, OVector, StorageMut,
     Vector, U1,
 };
 use num_traits::{One, Zero};
-use rand::Rng;
-use rand_distr::{uniform::SampleUniform, Bernoulli, Distribution, Standard};
 use thiserror::Error;
 
-use crate::core::{Domain, Function, Optimizer, Problem, Solver, System};
+use crate::core::{Domain, Function, Optimizer, Problem, Sample, Solver, System};
 
 use super::NelderMead;
 
@@ -75,30 +74,30 @@ impl<F: Problem> Default for LipoOptions<F> {
 }
 
 /// LIPO solver. See [module](self) documentation for more details.
-pub struct Lipo<F: Problem, R> {
+pub struct Lipo<F: Problem> {
     options: LipoOptions<F>,
     alpha: F::Field,
     xs: Vec<OVector<F::Field, Dynamic>>,
     ys: Vec<F::Field>,
     best: usize,
     k: F::Field,
     k_inf: F::Field,
-    rng: R,
-    bernoulli: Bernoulli,
+    rng: Rng,
+    p_explore: f64,
     tmp: OVector<F::Field, Dynamic>,
     x_tmp: OVector<F::Field, Dynamic>,
     local_optimizer: NelderMead<F>,
     iter: usize,
 }
 
-impl<F: Problem, R: Rng> Lipo<F, R> {
+impl<F: Problem> Lipo<F> {
     /// Initializes LIPO solver with default options.
-    pub fn new(f: &F, dom: &Domain<F::Field>, rng: R) -> Self {
+    pub fn new(f: &F, dom: &Domain<F::Field>, rng: Rng) -> Self {
         Self::with_options(f, dom, LipoOptions::default(), rng)
     }
 
     /// Initializes LIPO solver with given options.
-    pub fn with_options(f: &F, dom: &Domain<F::Field>, options: LipoOptions<F>, rng: R) -> Self {
+    pub fn with_options(f: &F, dom: &Domain<F::Field>, options: LipoOptions<F>, rng: Rng) -> Self {
         let dim = Dynamic::new(dom.dim());
 
         let p_explore = options.p_explore.clamp(0.0, 1.0);
@@ -117,7 +116,7 @@ impl<F: Problem, R: Rng> Lipo<F, R> {
             k: F::Field::zero(),
             k_inf: F::Field::zero(),
             rng,
-            bernoulli: Bernoulli::new(p_explore).unwrap(),
+            p_explore,
             tmp: OVector::zeros_generic(dim, U1::name()),
             x_tmp: OVector::zeros_generic(dim, U1::name()),
             local_optimizer: NelderMead::new(f, dom),
@@ -209,10 +208,9 @@ pub enum LipoError {
     InfiniteLipschitzConstant,
 }
 
-impl<F: Function, R: Rng> Lipo<F, R>
+impl<F: Function> Lipo<F>
 where
-    F::Field: SampleUniform,
-    Standard: Distribution<F::Field>,
+    F::Field: Sample,
 {
     fn next_inner<Sx>(
         &mut self,
@@ -236,7 +234,7 @@ where
             best,
             k,
             rng,
-            bernoulli,
+            p_explore,
             tmp,
             x_tmp,
             local_optimizer,
@@ -268,7 +266,7 @@ where
             // Exploitation mode is allowed only when there is enough points
             // evaluated and the Lipschitz constant is estimated. Then there is
             // randomness involved in choosing whether we explore or exploit.
-            if !initialization && *k != F::Field::zero() && !bernoulli.sample(rng) {
+            if !initialization && *k != F::Field::zero() && rng.f64() >= *p_explore {
                 debug!("exploitation mode");
 
                 let mut tmp_best = ys[*best];
@@ -378,10 +376,9 @@ where
     }
 }
 
-impl<F: Function, R: Rng> Optimizer<F> for Lipo<F, R>
+impl<F: Function> Optimizer<F> for Lipo<F>
 where
-    F::Field: SampleUniform,
-    Standard: Distribution<F::Field>,
+    F::Field: Sample,
 {
     const NAME: &'static str = "LIPO";
 
@@ -400,10 +397,9 @@ where
     }
 }
 
-impl<F: System, R: Rng> Solver<F> for Lipo<F, R>
+impl<F: System> Solver<F> for Lipo<F>
 where
-    F::Field: SampleUniform,
-    Standard: Distribution<F::Field>,
+    F::Field: Sample,
 {
     const NAME: &'static str = "LIPO";
 
@@ -428,8 +424,6 @@ where
 
 #[cfg(test)]
 mod tests {
-    use rand::{rngs::SmallRng, SeedableRng};
-
     use super::*;
 
     use crate::testing::*;
@@ -441,7 +435,7 @@ mod tests {
         let f = Sphere::new(n);
         let dom = f.domain();
         let eps = convert(1e-3);
-        let rng = SmallRng::from_seed([3; 32]);
+        let rng = Rng::with_seed(3);
         let mut options = LipoOptions::default();
         options
             .set_local_optimization_iters(0)