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<T, U> {
@@ -340,7 +341,7 @@ pub trait ArgminInv<T> {
 /// 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<R: Rng>(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<N, R, C>,
 {
     #[inline]
-    fn rand_from_range(min: &Self, max: &Self) -> OMatrix<N, R, C> {
+    fn rand_from_range<T: Rng>(min: &Self, max: &Self, rng: &mut T) -> OMatrix<N, R, C> {
         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 [<test_random_vec_ $t>]() {
                     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 [<test_random_vec_equal $t>]() {
                     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 [<test_random_vec_reverse_ $t>]() {
                     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<R: Rng>(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<R: Rng>(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 [<test_random_vec_ $t>]() {
                     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<R: Rng>(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 [<test_random_vec_ $t>]() {
                     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<R: Rng>(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<Vec<$t>> {
-            fn rand_from_range(min: &Self, max: &Self) -> Vec<Vec<$t>> {
+            fn rand_from_range<R: Rng>(min: &Self, max: &Self, rng: &mut R) -> Vec<Vec<$t>> {
                 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 [<test_random_vec_ $t>]() {
                     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::<Vec<$t>>::rand_from_range(&a, &b);
+                    let mut rng = rand::rngs::StdRng::seed_from_u64(42);
+                    let random = Vec::<Vec<$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/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\] <https://en.wikipedia.org/wiki/Particle_swarm_optimization>
 #[derive(Clone)]
 #[cfg_attr(feature = "serde1", derive(Serialize, Deserialize))]
-pub struct ParticleSwarm<P, F> {
+pub struct ParticleSwarm<P, F, R> {
     /// Inertia weight
     weight_inertia: F,
     /// Cognitive acceleration coefficient
@@ -61,9 +62,11 @@ pub struct ParticleSwarm<P, F> {
     bounds: (P, P),
     /// Number of particles
     num_particles: usize,
+    /// Random number generator
+    rng_generator: R,
 }
 
-impl<P, F> ParticleSwarm<P, F>
+impl<P, F> ParticleSwarm<P, F, rand::rngs::StdRng>
 where
     P: Clone + SyncAlias + ArgminSub<P, P> + ArgminMul<F, P> + ArgminRandom + ArgminZeroLike,
     F: ArgminFloat,
@@ -91,7 +94,7 @@ where
     /// # use argmin::solver::particleswarm::ParticleSwarm;
     /// # let lower_bound: Vec<f64> = vec![-1.0, -1.0];
     /// # let upper_bound: Vec<f64> = 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<P, F, R0> ParticleSwarm<P, F, R0>
+where
+    P: Clone + SyncAlias + ArgminSub<P, P> + ArgminMul<F, P> + 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<f64> = vec![-1.0, -1.0];
+    /// # let upper_bound: Vec<f64> = 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<R1: Rng>(self, generator: R1) -> ParticleSwarm<P, F, R1> {
+        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<P, F, R> ParticleSwarm<P, F, R>
+where
+    P: Clone + SyncAlias + ArgminSub<P, P> + ArgminMul<F, P> + 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<f64> = vec![-1.0, -1.0];
     /// # let upper_bound: Vec<f64> = 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<f64> = vec![-1.0, -1.0];
     /// # let upper_bound: Vec<f64> = 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<f64> = vec![-1.0, -1.0];
     /// # let upper_bound: Vec<f64> = 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<P>, Vec<P>) {
+    fn initialize_positions_and_velocities(&mut self) -> (Vec<P>, Vec<P>) {
         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<O, P, F> Solver<O, PopulationState<Particle<P, F>, F>> for ParticleSwarm<P, F>
+impl<O, P, F, R> Solver<O, PopulationState<Particle<P, F>, F>> for ParticleSwarm<P, F, R>
 where
     O: CostFunction<Param = P, Output = F> + 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<Vec<f64>, f64>);
+    test_trait_impl!(
+        particleswarm,
+        ParticleSwarm<Vec<f64>, f64, rand::rngs::StdRng>
+    );
 
     #[test]
     fn test_new() {
         let lower_bound: Vec<f64> = vec![-1.0, -1.0];
         let upper_bound: Vec<f64> = 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<f64> = vec![-1.0, -1.0];
         let upper_bound: Vec<f64> = 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<f64> = vec![-1.0, -1.0];
         let upper_bound: Vec<f64> = 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<f64> = vec![-1.0, -1.0];
         let upper_bound: Vec<f64> = 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<Particle<Vec<f64>, 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<f64> = 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<Particle<Vec<f64>, 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<f64> = vec![-1.0, -1.0];
         let upper_bound: Vec<f64> = 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<Particle<Vec<f64>, 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<f64> = vec![-1.0, -1.0];
         let upper_bound: Vec<f64> = 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<Particle<Vec<f64>, f64>, f64> = PopulationState::new();
 
         // init