From 220e0d58575b26b6bbc5a41c59b1669794a11abc Mon Sep 17 00:00:00 2001 From: Jon Date: Sun, 7 Jan 2024 11:57:40 +0100 Subject: [PATCH] allow user to set ParticleSwarm's random number generator (#383) * allow user to set ParticleSwarm random number generator * use R for Random and make default nondeterministic use T where R is not possible * update comment on default rng generator Co-authored-by: Stefan Kroboth --------- Co-authored-by: Stefan Kroboth --- argmin-math/src/lib.rs | 3 +- argmin-math/src/nalgebra_m/random.rs | 17 +++-- argmin-math/src/ndarray_m/random.rs | 17 +++-- argmin-math/src/primitives/random.rs | 8 ++- argmin-math/src/vec/random.rs | 15 +++-- argmin/src/solver/particleswarm/mod.rs | 93 ++++++++++++++++++++------ 6 files changed, 106 insertions(+), 47 deletions(-) diff --git a/argmin-math/src/lib.rs b/argmin-math/src/lib.rs index d54e98ee6..f686fd658 100644 --- a/argmin-math/src/lib.rs +++ b/argmin-math/src/lib.rs @@ -234,6 +234,7 @@ mod vec; pub use crate::vec::*; use anyhow::Error; +use rand::Rng; /// Dot/scalar product of `T` and `self` pub trait ArgminDot { @@ -340,7 +341,7 @@ pub trait ArgminInv { /// Create a random number pub trait ArgminRandom { /// Get a random element between min and max, - fn rand_from_range(min: &Self, max: &Self) -> Self; + fn rand_from_range(min: &Self, max: &Self, rng: &mut R) -> Self; } /// Minimum and Maximum of type `T` diff --git a/argmin-math/src/nalgebra_m/random.rs b/argmin-math/src/nalgebra_m/random.rs index 37d89097f..46e039661 100644 --- a/argmin-math/src/nalgebra_m/random.rs +++ b/argmin-math/src/nalgebra_m/random.rs @@ -22,12 +22,10 @@ where DefaultAllocator: Allocator, { #[inline] - fn rand_from_range(min: &Self, max: &Self) -> OMatrix { + fn rand_from_range(min: &Self, max: &Self, rng: &mut T) -> OMatrix { assert!(!min.is_empty()); assert_eq!(min.shape(), max.shape()); - let mut rng = rand::thread_rng(); - Self::from_iterator_generic( R::from_usize(min.nrows()), C::from_usize(min.ncols()), @@ -53,6 +51,7 @@ mod tests { use super::*; use nalgebra::{Matrix2x3, Vector3}; use paste::item; + use rand::SeedableRng; macro_rules! make_test { ($t:ty) => { @@ -61,7 +60,8 @@ mod tests { fn []() { let a = Vector3::new(1 as $t, 2 as $t, 3 as $t); let b = Vector3::new(2 as $t, 3 as $t, 4 as $t); - let random = Vector3::<$t>::rand_from_range(&a, &b); + let mut rng = rand::rngs::StdRng::seed_from_u64(42); + let random = Vector3::<$t>::rand_from_range(&a, &b, &mut rng); for i in 0..3 { assert!(random[i] >= a[i]); assert!(random[i] <= b[i]); @@ -74,7 +74,8 @@ mod tests { fn []() { let a = Vector3::new(1 as $t, 2 as $t, 3 as $t); let b = Vector3::new(1 as $t, 2 as $t, 3 as $t); - let random = Vector3::<$t>::rand_from_range(&a, &b); + let mut rng = rand::rngs::StdRng::seed_from_u64(42); + let random = Vector3::<$t>::rand_from_range(&a, &b, &mut rng); for i in 0..3 { assert!((random[i] as f64 - a[i] as f64).abs() < std::f64::EPSILON); assert!((random[i] as f64 - b[i] as f64).abs() < std::f64::EPSILON); @@ -87,7 +88,8 @@ mod tests { fn []() { let b = Vector3::new(1 as $t, 2 as $t, 3 as $t); let a = Vector3::new(2 as $t, 3 as $t, 4 as $t); - let random = Vector3::<$t>::rand_from_range(&a, &b); + let mut rng = rand::rngs::StdRng::seed_from_u64(42); + let random = Vector3::<$t>::rand_from_range(&a, &b, &mut rng); for i in 0..3 { assert!(random[i] >= b[i]); assert!(random[i] <= a[i]); @@ -106,7 +108,8 @@ mod tests { 2 as $t, 4 as $t, 6 as $t, 3 as $t, 5 as $t, 7 as $t ); - let random = Matrix2x3::<$t>::rand_from_range(&a, &b); + let mut rng = rand::rngs::StdRng::seed_from_u64(42); + let random = Matrix2x3::<$t>::rand_from_range(&a, &b, &mut rng); for i in 0..3 { for j in 0..2 { assert!(random[(j, i)] >= a[(j, i)]); diff --git a/argmin-math/src/ndarray_m/random.rs b/argmin-math/src/ndarray_m/random.rs index a13b45c00..5f44b4b16 100644 --- a/argmin-math/src/ndarray_m/random.rs +++ b/argmin-math/src/ndarray_m/random.rs @@ -5,19 +5,17 @@ // http://opensource.org/licenses/MIT>, at your option. This file may not be // copied, modified, or distributed except according to those terms. -use rand::Rng; +use rand::{Rng, SeedableRng}; use crate::ArgminRandom; macro_rules! make_random { ($t:ty) => { impl ArgminRandom for ndarray::Array1<$t> { - fn rand_from_range(min: &Self, max: &Self) -> ndarray::Array1<$t> { + fn rand_from_range(min: &Self, max: &Self, rng: &mut R) -> ndarray::Array1<$t> { assert!(!min.is_empty()); assert_eq!(min.len(), max.len()); - let mut rng = rand::thread_rng(); - ndarray::Array1::from_iter(min.iter().zip(max.iter()).map(|(a, b)| { // Do not require a < b: @@ -35,12 +33,10 @@ macro_rules! make_random { } impl ArgminRandom for ndarray::Array2<$t> { - fn rand_from_range(min: &Self, max: &Self) -> ndarray::Array2<$t> { + fn rand_from_range(min: &Self, max: &Self, rng: &mut R) -> ndarray::Array2<$t> { assert!(!min.is_empty()); assert_eq!(min.raw_dim(), max.raw_dim()); - let mut rng = rand::thread_rng(); - ndarray::Array2::from_shape_fn(min.raw_dim(), |(i, j)| { let a = min.get((i, j)).unwrap(); let b = max.get((i, j)).unwrap(); @@ -78,6 +74,7 @@ mod tests { use super::*; use ndarray::{array, Array1, Array2}; use paste::item; + use rand::SeedableRng; macro_rules! make_test { ($t:ty) => { @@ -86,7 +83,8 @@ mod tests { fn []() { let a = array![1 as $t, 2 as $t, 4 as $t]; let b = array![2 as $t, 3 as $t, 5 as $t]; - let random = Array1::<$t>::rand_from_range(&a, &b); + let mut rng = rand::rngs::StdRng::seed_from_u64(42); + let random = Array1::<$t>::rand_from_range(&a, &b, &mut rng); for i in 0..3usize { assert!(random[i] >= a[i]); assert!(random[i] <= b[i]); @@ -105,7 +103,8 @@ mod tests { [2 as $t, 3 as $t, 5 as $t], [3 as $t, 4 as $t, 6 as $t] ]; - let random = Array2::<$t>::rand_from_range(&a, &b); + let mut rng = rand::rngs::StdRng::seed_from_u64(42); + let random = Array2::<$t>::rand_from_range(&a, &b, &mut rng); for i in 0..3 { for j in 0..2 { assert!(random[(j, i)] >= a[(j, i)]); diff --git a/argmin-math/src/primitives/random.rs b/argmin-math/src/primitives/random.rs index eb052da0e..af370cee5 100644 --- a/argmin-math/src/primitives/random.rs +++ b/argmin-math/src/primitives/random.rs @@ -12,8 +12,8 @@ macro_rules! make_random { ($t:ty) => { impl ArgminRandom for $t { #[inline] - fn rand_from_range(min: &Self, max: &Self) -> $t { - rand::thread_rng().gen_range(*min..*max) + fn rand_from_range(min: &Self, max: &Self, rng: &mut R) -> $t { + rng.gen_range(*min..*max) } } }; @@ -36,6 +36,7 @@ make_random!(usize); mod tests { use super::*; use paste::item; + use rand::SeedableRng; macro_rules! make_test { ($t:ty) => { @@ -44,7 +45,8 @@ mod tests { fn []() { let a = 1 as $t; let b = 2 as $t; - let random = $t::rand_from_range(&a, &b); + let mut rng = rand::rngs::StdRng::seed_from_u64(42); + let random = $t::rand_from_range(&a, &b, &mut rng); assert!(random >= a); assert!(random <= b); } diff --git a/argmin-math/src/vec/random.rs b/argmin-math/src/vec/random.rs index 044d664b0..60e34225a 100644 --- a/argmin-math/src/vec/random.rs +++ b/argmin-math/src/vec/random.rs @@ -11,12 +11,10 @@ use rand::Rng; macro_rules! make_random { ($t:ty) => { impl ArgminRandom for Vec<$t> { - fn rand_from_range(min: &Self, max: &Self) -> Vec<$t> { + fn rand_from_range(min: &Self, max: &Self, rng: &mut R) -> Vec<$t> { assert!(!min.is_empty()); assert_eq!(min.len(), max.len()); - let mut rng = rand::thread_rng(); - min.iter() .zip(max.iter()) .map(|(a, b)| { @@ -37,12 +35,12 @@ macro_rules! make_random { } impl ArgminRandom for Vec> { - fn rand_from_range(min: &Self, max: &Self) -> Vec> { + fn rand_from_range(min: &Self, max: &Self, rng: &mut R) -> Vec> { assert!(!min.is_empty()); assert_eq!(min.len(), max.len()); min.iter() .zip(max.iter()) - .map(|(a, b)| Vec::<$t>::rand_from_range(a, b)) + .map(|(a, b)| Vec::<$t>::rand_from_range(a, b, rng)) .collect() } } @@ -66,6 +64,7 @@ make_random!(usize); mod tests { use super::*; use paste::item; + use rand::SeedableRng; macro_rules! make_test { ($t:ty) => { @@ -74,7 +73,8 @@ mod tests { fn []() { let a = vec![1 as $t, 2 as $t, 4 as $t]; let b = vec![2 as $t, 3 as $t, 5 as $t]; - let random = Vec::<$t>::rand_from_range(&a, &b); + let mut rng = rand::rngs::StdRng::seed_from_u64(42); + let random = Vec::<$t>::rand_from_range(&a, &b, &mut rng); for i in 0..3usize { assert!(random[i] >= a[i]); assert!(random[i] <= b[i]); @@ -93,7 +93,8 @@ mod tests { vec![2 as $t, 3 as $t, 5 as $t], vec![3 as $t, 4 as $t, 6 as $t] ]; - let random = Vec::>::rand_from_range(&a, &b); + let mut rng = rand::rngs::StdRng::seed_from_u64(42); + let random = Vec::>::rand_from_range(&a, &b, &mut rng); for i in 0..3 { for j in 0..2 { assert!(random[j][i] >= a[j][i]); diff --git a/argmin/src/solver/particleswarm/mod.rs b/argmin/src/solver/particleswarm/mod.rs index 42e1021fb..83a3ef932 100644 --- a/argmin/src/solver/particleswarm/mod.rs +++ b/argmin/src/solver/particleswarm/mod.rs @@ -25,6 +25,7 @@ use crate::core::{ KV, }; use argmin_math::{ArgminAdd, ArgminMinMax, ArgminMul, ArgminRandom, ArgminSub, ArgminZeroLike}; +use rand::{Rng, SeedableRng}; #[cfg(feature = "serde1")] use serde::{Deserialize, Serialize}; @@ -50,7 +51,7 @@ use serde::{Deserialize, Serialize}; /// \[1\] #[derive(Clone)] #[cfg_attr(feature = "serde1", derive(Serialize, Deserialize))] -pub struct ParticleSwarm { +pub struct ParticleSwarm { /// Inertia weight weight_inertia: F, /// Cognitive acceleration coefficient @@ -61,9 +62,11 @@ pub struct ParticleSwarm { bounds: (P, P), /// Number of particles num_particles: usize, + /// Random number generator + rng_generator: R, } -impl ParticleSwarm +impl ParticleSwarm where P: Clone + SyncAlias + ArgminSub + ArgminMul + ArgminRandom + ArgminZeroLike, F: ArgminFloat, @@ -91,7 +94,7 @@ where /// # use argmin::solver::particleswarm::ParticleSwarm; /// # let lower_bound: Vec = vec![-1.0, -1.0]; /// # let upper_bound: Vec = vec![1.0, 1.0]; - /// let pso: ParticleSwarm<_, f64> = ParticleSwarm::new((lower_bound, upper_bound), 40); + /// let pso: ParticleSwarm<_, f64, _> = ParticleSwarm::new((lower_bound, upper_bound), 40); /// ``` pub fn new(bounds: (P, P), num_particles: usize) -> Self { ParticleSwarm { @@ -100,9 +103,52 @@ where weight_social: float!(0.5 + 2.0f64.ln()), bounds, num_particles, + rng_generator: rand::rngs::StdRng::from_entropy(), } } +} +impl ParticleSwarm +where + P: Clone + SyncAlias + ArgminSub + ArgminMul + ArgminRandom + ArgminZeroLike, + F: ArgminFloat, + R0: Rng, +{ + /// Set the random number generator + /// + /// Defaults to `rand::rngs::StdRng::from_entropy()` + /// + /// # Example + /// ``` + /// # use argmin::solver::particleswarm::ParticleSwarm; + /// # use argmin::core::Error; + /// # use rand::SeedableRng; + /// # fn main() -> Result<(), Error> { + /// # let lower_bound: Vec = vec![-1.0, -1.0]; + /// # let upper_bound: Vec = vec![1.0, 1.0]; + /// let pso: ParticleSwarm<_, f64, _> = + /// ParticleSwarm::new((lower_bound, upper_bound), 40) + /// .with_rng_generator(rand_xoshiro::Xoroshiro128Plus::seed_from_u64(1729)); + /// # Ok(()) + /// # } + /// ``` + pub fn with_rng_generator(self, generator: R1) -> ParticleSwarm { + ParticleSwarm { + weight_inertia: self.weight_inertia, + weight_cognitive: self.weight_cognitive, + weight_social: self.weight_social, + bounds: self.bounds, + num_particles: self.num_particles, + rng_generator: generator, + } + } +} +impl ParticleSwarm +where + P: Clone + SyncAlias + ArgminSub + ArgminMul + ArgminRandom + ArgminZeroLike, + F: ArgminFloat, + R: Rng, +{ /// Set inertia factor on particle velocity /// /// Defaults to `1/(2 * ln(2))`. @@ -115,7 +161,7 @@ where /// # fn main() -> Result<(), Error> { /// # let lower_bound: Vec = vec![-1.0, -1.0]; /// # let upper_bound: Vec = vec![1.0, 1.0]; - /// let pso: ParticleSwarm<_, f64> = + /// let pso: ParticleSwarm<_, f64, _> = /// ParticleSwarm::new((lower_bound, upper_bound), 40).with_inertia_factor(0.5)?; /// # Ok(()) /// # } @@ -143,7 +189,7 @@ where /// # fn main() -> Result<(), Error> { /// # let lower_bound: Vec = vec![-1.0, -1.0]; /// # let upper_bound: Vec = vec![1.0, 1.0]; - /// let pso: ParticleSwarm<_, f64> = + /// let pso: ParticleSwarm<_, f64, _> = /// ParticleSwarm::new((lower_bound, upper_bound), 40).with_cognitive_factor(1.1)?; /// # Ok(()) /// # } @@ -171,7 +217,7 @@ where /// # fn main() -> Result<(), Error> { /// # let lower_bound: Vec = vec![-1.0, -1.0]; /// # let upper_bound: Vec = vec![1.0, 1.0]; - /// let pso: ParticleSwarm<_, f64> = + /// let pso: ParticleSwarm<_, f64, _> = /// ParticleSwarm::new((lower_bound, upper_bound), 40).with_social_factor(1.1)?; /// # Ok(()) /// # } @@ -214,23 +260,23 @@ where } /// Initializes positions and velocities for all particles - fn initialize_positions_and_velocities(&self) -> (Vec

, Vec

) { + fn initialize_positions_and_velocities(&mut self) -> (Vec

, Vec

) { let (min, max) = &self.bounds; let delta = max.sub(min); let delta_neg = delta.mul(&float!(-1.0)); ( (0..self.num_particles) - .map(|_| P::rand_from_range(min, max)) + .map(|_| P::rand_from_range(min, max, &mut self.rng_generator)) .collect(), (0..self.num_particles) - .map(|_| P::rand_from_range(&delta_neg, &delta)) + .map(|_| P::rand_from_range(&delta_neg, &delta, &mut self.rng_generator)) .collect(), ) } } -impl Solver, F>> for ParticleSwarm +impl Solver, F>> for ParticleSwarm where O: CostFunction + SyncAlias, P: SerializeAlias @@ -243,6 +289,7 @@ where + ArgminRandom + ArgminMinMax, F: ArgminFloat, + R: Rng, { const NAME: &'static str = "Particle Swarm Optimization"; @@ -315,13 +362,15 @@ where // ad 2) let to_optimum = p.best_position.sub(&p.position); - let pull_to_optimum = P::rand_from_range(&zero, &to_optimum); + let pull_to_optimum = + P::rand_from_range(&zero, &to_optimum, &mut self.rng_generator); let pull_to_optimum = pull_to_optimum.mul(&self.weight_cognitive); // ad 3) let to_global_optimum = best_particle.position.sub(&p.position); let pull_to_global_optimum = - P::rand_from_range(&zero, &to_global_optimum).mul(&self.weight_social); + P::rand_from_range(&zero, &to_global_optimum, &mut self.rng_generator) + .mul(&self.weight_social); p.velocity = momentum.add(&pull_to_optimum).add(&pull_to_global_optimum); let new_position = p.position.add(&p.velocity); @@ -408,13 +457,16 @@ mod tests { use crate::test_trait_impl; use approx::assert_relative_eq; - test_trait_impl!(particleswarm, ParticleSwarm, f64>); + test_trait_impl!( + particleswarm, + ParticleSwarm, f64, rand::rngs::StdRng> + ); #[test] fn test_new() { let lower_bound: Vec = vec![-1.0, -1.0]; let upper_bound: Vec = vec![1.0, 1.0]; - let pso: ParticleSwarm<_, f64> = + let pso: ParticleSwarm<_, f64, rand::rngs::StdRng> = ParticleSwarm::new((lower_bound.clone(), upper_bound.clone()), 40); let ParticleSwarm { weight_inertia, @@ -422,6 +474,7 @@ mod tests { weight_social, bounds, num_particles, + .. } = pso; assert_relative_eq!( @@ -538,7 +591,7 @@ mod tests { let lower_bound: Vec = vec![-1.0, -1.0]; let upper_bound: Vec = vec![1.0, 1.0]; let num_particles = 100; - let pso: ParticleSwarm<_, f64> = + let mut pso: ParticleSwarm<_, f64, _> = ParticleSwarm::new((lower_bound, upper_bound), num_particles); let (positions, velocities) = pso.initialize_positions_and_velocities(); @@ -565,7 +618,7 @@ mod tests { let lower_bound: Vec = vec![-1.0, -1.0]; let upper_bound: Vec = vec![1.0, 1.0]; let num_particles = 10; - let mut pso: ParticleSwarm<_, f64> = + let mut pso: ParticleSwarm<_, f64, _> = ParticleSwarm::new((lower_bound, upper_bound), num_particles); struct PsoProblem { @@ -630,7 +683,7 @@ mod tests { fn test_init_provided_population_wrong_size() { let lower_bound: Vec = vec![-1.0, -1.0]; let upper_bound: Vec = vec![1.0, 1.0]; - let mut pso: ParticleSwarm<_, f64> = ParticleSwarm::new((lower_bound, upper_bound), 40); + let mut pso: ParticleSwarm<_, f64, _> = ParticleSwarm::new((lower_bound, upper_bound), 40); let state: PopulationState, f64>, f64> = PopulationState::new() .population(vec![Particle::new(vec![1.0, 2.0], 12.0, vec![0.1, 0.3])]); let res = pso.init(&mut Problem::new(TestProblem::new()), state); @@ -650,7 +703,7 @@ mod tests { let upper_bound: Vec = vec![1.0, 1.0]; let particle_a = Particle::new(vec![1.0, 2.0], 12.0, vec![0.1, 0.3]); let particle_b = Particle::new(vec![2.0, 3.0], 10.0, vec![0.2, 0.4]); - let mut pso: ParticleSwarm<_, f64> = ParticleSwarm::new((lower_bound, upper_bound), 2); + let mut pso: ParticleSwarm<_, f64, _> = ParticleSwarm::new((lower_bound, upper_bound), 2); let state: PopulationState, f64>, f64> = PopulationState::new().population(vec![particle_a.clone(), particle_b.clone()]); let res = pso.init(&mut Problem::new(TestProblem::new()), state); @@ -668,7 +721,7 @@ mod tests { fn test_init_random_population() { let lower_bound: Vec = vec![-1.0, -1.0]; let upper_bound: Vec = vec![1.0, 1.0]; - let mut pso: ParticleSwarm<_, f64> = ParticleSwarm::new((lower_bound, upper_bound), 40); + let mut pso: ParticleSwarm<_, f64, _> = ParticleSwarm::new((lower_bound, upper_bound), 40); let state: PopulationState, f64>, f64> = PopulationState::new(); let res = pso.init(&mut Problem::new(TestProblem::new()), state); assert!(res.is_ok()); @@ -707,7 +760,7 @@ mod tests { // setup let lower_bound: Vec = vec![-1.0, -1.0]; let upper_bound: Vec = vec![1.0, 1.0]; - let mut pso: ParticleSwarm<_, f64> = ParticleSwarm::new((lower_bound, upper_bound), 100); + let mut pso: ParticleSwarm<_, f64, _> = ParticleSwarm::new((lower_bound, upper_bound), 100); let state: PopulationState, f64>, f64> = PopulationState::new(); // init