Merge pull request #194 from mahf-opt/add-swarm-operators

Add more swarm-based metaheuristics

src/components/mapping/mod.rs

@@ -22,8 +22,8 @@
 //! using [`ValueOf<Progress<ValueOf<Iterations>>>`] as input lens and
 //! [`ValueOf<InertiaWeight>`] as output lens:
 //!
-//! [`InertiaWeight`]: crate::components::swarm::InertiaWeight
-//! [`ParticleVelocitiesUpdate`]: crate::components::swarm::ParticleVelocitiesUpdate
+//! [`InertiaWeight`]: crate::components::swarm::pso::InertiaWeight
+//! [`ParticleVelocitiesUpdate`]: crate::components::swarm::pso::ParticleVelocitiesUpdate
 //! [`ValueOf<Progress<ValueOf<Iterations>>>`]: crate::lens::ValueOf
 //! [`ValueOf<InertiaWeight>`]: crate::lens::ValueOf
 //!
@@ -40,7 +40,7 @@
 //!     0.4,
 //!     0.9,
 //!     ValueOf::<common::Progress<ValueOf<common::Iterations>>>::new(),
-//!     ValueOf::<swarm::InertiaWeight<swarm::ParticleVelocitiesUpdate>>::new(),
+//!     ValueOf::<swarm::pso::InertiaWeight<swarm::pso::ParticleVelocitiesUpdate>>::new(),
 //! )
 //! # }
 //! ```

src/components/replacement/bh.rs

@@ -0,0 +1,61 @@
+//! Replacement components for the Black Hole algorithm (BH).
+
+use rand::Rng;
+use serde::{Deserialize, Serialize};
+
+use crate::{
+    problems::LimitedVectorProblem, utils::squared_euclidean, Component, ExecResult,
+    SingleObjectiveProblem, State,
+};
+
+#[derive(Clone, Serialize, Deserialize)]
+pub struct EventHorizon;
+
+impl EventHorizon {
+    pub fn new<P>() -> Box<dyn Component<P>>
+    where
+        P: LimitedVectorProblem<Element = f64>,
+        P: SingleObjectiveProblem,
+    {
+        Box::new(Self)
+    }
+}
+
+impl<P> Component<P> for EventHorizon
+where
+    P: LimitedVectorProblem<Element = f64>,
+    P: SingleObjectiveProblem,
+{
+    fn execute(&self, problem: &P, state: &mut State<P>) -> ExecResult<()> {
+        let mut populations = state.populations_mut();
+        let offspring = populations.current_mut();
+
+        let f_bh = state.best_objective_value().unwrap().value();
+
+        let fitness_sum = offspring.iter().map(|x| x.objective().value()).sum::<f64>();
+        let radius = f_bh / fitness_sum;
+
+        let (index, best) = offspring
+            .iter()
+            .enumerate()
+            .min_by_key(|(_u, i)| i.objective())
+            .unwrap();
+        let distances = offspring
+            .iter()
+            .map(|o| squared_euclidean(o.solution(), best.solution()).sqrt())
+            .collect::<Vec<f64>>();
+
+        for (u, i) in offspring.iter_mut().enumerate() {
+            // do not replace best individual
+            if distances[u] < radius && u != index {
+                let rand: Vec<f64> = problem
+                    .domain()
+                    .iter()
+                    .map(|d| state.random_mut().gen_range(d.clone()))
+                    .collect();
+                *i.solution_mut() = rand;
+            }
+        }
+        Ok(())
+    }
+}

src/components/replacement/mod.rs

@@ -8,6 +8,7 @@ use crate::{
     Individual, Problem, State,
 };
 
+pub mod bh;
 pub mod common;
 pub mod sa;
 

src/components/replacement/sa.rs

@@ -75,7 +75,6 @@ impl<P: SingleObjectiveProblem> Component<P> for ExponentialAnnealingAcceptance
         } else {
             populations.pop();
         }
-
         Ok(())
     }
 }

src/components/swarm/bh.rs

@@ -0,0 +1,85 @@
+use rand::distributions::{Distribution, Uniform};
+use serde::Serialize;
+
+use crate::{
+    component::ExecResult,
+    components::Component,
+    identifier::{Global, Identifier, PhantomId},
+    population::{AsSolutionsMut, BestIndividual},
+    problems::LimitedVectorProblem,
+    SingleObjectiveProblem, State,
+};
+
+/// Updates the positions in the black hole algorithm.
+///
+/// Originally proposed for, and used as operator in [`bh`].
+/// The operator is similar to the [`ParticleVelocitiesUpdate`] in [`pso`].
+/// Specifically, they are identical when for [`pso`]:
+/// inertia weight = 0,
+/// c_1 = 0,
+/// c_2 = 1,
+/// v_max = 1
+///
+/// [`bh`]: crate::heuristics::bh
+/// [`ParticleVelocitiesUpdate`]: crate::components::swarm::pso::ParticleVelocitiesUpdate`
+/// [`pso`]: crate::heuristics::pso
+#[derive(Clone, Serialize)]
+pub struct BlackHoleParticlesUpdate<I: Identifier = Global> {
+    id: PhantomId<I>,
+}
+
+impl<I: Identifier> BlackHoleParticlesUpdate<I> {
+    pub fn from_params() -> Self {
+        Self {
+            id: PhantomId::default(),
+        }
+    }
+
+    pub fn new_with_id<P>() -> Box<dyn Component<P>>
+    where
+        P: SingleObjectiveProblem + LimitedVectorProblem<Element = f64>,
+    {
+        Box::new(Self::from_params())
+    }
+}
+
+impl BlackHoleParticlesUpdate<Global> {
+    pub fn new<P>() -> Box<dyn Component<P>>
+    where
+        P: SingleObjectiveProblem + LimitedVectorProblem<Element = f64>,
+    {
+        Self::new_with_id()
+    }
+}
+
+impl<P, I> Component<P> for BlackHoleParticlesUpdate<I>
+where
+    P: SingleObjectiveProblem + LimitedVectorProblem<Element = f64>,
+    I: Identifier,
+{
+    fn init(&self, _problem: &P, _state: &mut State<P>) -> ExecResult<()> {
+        Ok(())
+    }
+
+    fn execute(&self, _problem: &P, state: &mut State<P>) -> ExecResult<()> {
+        let distribution = Uniform::new(0.0, 1.0);
+
+        // Get necessary state like global best `xg`
+        let best = state.populations().current().best_individual().cloned();
+        let best_ind = best.unwrap();
+        let xg = best_ind.solution();
+        let mut population = state.populations_mut();
+        let xs = population.current_mut().as_solutions_mut();
+
+        // Perform the update step.
+        for x in xs {
+            for i in 0..x.len() {
+                // Calculate change in position
+                let pos = distribution.sample(&mut *state.random_mut()) * (xg[i] - x[i]);
+                // Add value to particle position
+                x[i] += pos;
+            }
+        }
+        Ok(())
+    }
+}

src/components/swarm/fa.rs

@@ -0,0 +1,137 @@
+use std::array::from_mut;
+
+use better_any::{Tid, TidAble};
+use derive_more::{Deref, DerefMut};
+use itertools::izip;
+use rand::Rng;
+use serde::Serialize;
+
+use crate::{
+    component::ExecResult,
+    components::Component,
+    identifier::{Global, Identifier, PhantomId},
+    problems::LimitedVectorProblem,
+    state::{common, common::Evaluator},
+    CustomState, State,
+};
+
+/// Updates the and firefly positions.
+///
+/// Originally proposed for, and used as operator in [`fa`].
+///
+/// Uses the [`RandomizationParameter`].
+///
+/// [`fa`]: crate::heuristics::fa
+#[derive(Clone, Serialize)]
+pub struct FireflyPositionsUpdate<I: Identifier = Global> {
+    pub alpha: f64,
+    pub beta: f64,
+    pub gamma: f64,
+    id: PhantomId<I>,
+}
+
+impl<I: Identifier> FireflyPositionsUpdate<I> {
+    pub fn from_params(alpha: f64, beta: f64, gamma: f64) -> Self {
+        Self {
+            alpha,
+            beta,
+            gamma,
+            id: PhantomId::default(),
+        }
+    }
+
+    pub fn new_with_id<P>(alpha: f64, beta: f64, gamma: f64) -> Box<dyn Component<P>>
+    where
+        P: LimitedVectorProblem<Element = f64>,
+    {
+        Box::new(Self::from_params(alpha, beta, gamma))
+    }
+}
+
+impl FireflyPositionsUpdate<Global> {
+    pub fn new<P>(alpha: f64, beta: f64, gamma: f64) -> Box<dyn Component<P>>
+    where
+        P: LimitedVectorProblem<Element = f64>,
+    {
+        Self::new_with_id(alpha, beta, gamma)
+    }
+}
+
+impl<P, I> Component<P> for FireflyPositionsUpdate<I>
+where
+    P: LimitedVectorProblem<Element = f64>,
+    I: Identifier,
+{
+    fn init(&self, _problem: &P, state: &mut State<P>) -> ExecResult<()> {
+        state.insert(RandomizationParameter(self.alpha));
+        Ok(())
+    }
+
+    fn execute(&self, problem: &P, state: &mut State<P>) -> ExecResult<()> {
+        // Prepare parameters
+        let &Self { beta, gamma, .. } = self;
+        let a = state.get_value::<RandomizationParameter>();
+
+        // Get population from state
+        let mut individuals = state.populations_mut().pop();
+
+        // scale for adapting to problem domain
+        let scales = problem
+            .domain()
+            .iter()
+            .map(|p| (p.end - p.start).abs())
+            .collect::<Vec<f64>>();
+
+        // Perform the update step.
+        for i in 0..individuals.len() {
+            for j in 0..individuals.len() {
+                // if individual j is "more attractive" (i.e. has lower fitness), move towards j
+                if individuals[i].objective() > individuals[j].objective() {
+                    // draw random values from uniform distribution between 0 and 1
+                    // according to paper: also possible to use normal distribution, depending on problem
+                    let rands: Vec<f64> = (0..problem.dimension())
+                        .map(|_| state.random_mut().gen_range(0.0..1.0))
+                        .collect();
+                    let mut current = individuals[i].clone();
+                    izip!(
+                        current.solution_mut(),
+                        individuals[j].solution(),
+                        &scales,
+                        rands
+                    )
+                    .map(|(xi, xj, scale, rand)| {
+                        let pos = beta * (-gamma * (*xi - xj).powf(2.0)).exp() * (xj - *xi)
+                            + a * (rand - 0.5) * scale;
+                        (xi, pos)
+                    })
+                    .for_each(|(xi, pos)| *xi += pos);
+                    individuals[i] = current;
+
+                    state.holding::<Evaluator<P, I>>(
+                        |evaluator: &mut Evaluator<P, I>, state| {
+                            evaluator.as_inner_mut().evaluate(
+                                problem,
+                                state,
+                                from_mut(&mut individuals[i]),
+                            );
+                            Ok(())
+                        },
+                    )?;
+                    *state.borrow_value_mut::<common::Evaluations>() += 1;
+                }
+            }
+        }
+        state.populations_mut().push(individuals);
+        Ok(())
+    }
+}
+
+/// The randomization parameter used to update the firefly positions.
+#[derive(Deref, DerefMut, Tid)]
+pub struct RandomizationParameter(
+    #[deref]
+    #[deref_mut]
+    pub f64,
+);
+
+impl CustomState<'_> for RandomizationParameter {}

src/components/swarm/mod.rs

@@ -0,0 +1,3 @@
+pub mod bh;
+pub mod fa;
+pub mod pso;

src/components/swarm.rs → src/components/swarm/pso.rs

@@ -125,9 +125,9 @@ pub struct ParticleVelocitiesUpdate<I: Identifier = Global> {
 
 impl<I: Identifier> ParticleVelocitiesUpdate<I> {
     pub fn from_params(weight: f64, c_1: f64, c_2: f64, v_max: f64) -> ExecResult<Self> {
-        ensure!(weight > 0., "`weight` must be > 0, but was {}", weight);
-        ensure!(c_1 > 0., "`c_1` must be > 0, but was {}", c_1);
-        ensure!(c_2 > 0., "`c_2` must be > 0, but was {}", c_2);
+        ensure!(weight >= 0., "`weight` must be >= 0, but was {}", weight);
+        ensure!(c_1 >= 0., "`c_1` must be >= 0, but was {}", c_1);
+        ensure!(c_2 >= 0., "`c_2` must be >= 0, but was {}", c_2);
         ensure!(v_max > 0., "`v_max` must be > 0, but was {}", v_max);
         Ok(Self {
             weight,