argmin/solver/particleswarm/
mod.rs

1// Copyright 2018-2020 argmin developers
2//
3// Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
4// http://apache.org/licenses/LICENSE-2.0> or the MIT license <LICENSE-MIT or
5// http://opensource.org/licenses/MIT>, at your option. This file may not be
6// copied, modified, or distributed except according to those terms.
7
8//! # References:
9//!
10//! TODO
11
12use crate::prelude::*;
13use serde::{de::DeserializeOwned, Deserialize, Serialize};
14use std;
15use std::default::Default;
16
17/// Particle Swarm Optimization (PSO)
18///
19/// [Example](https://github.com/argmin-rs/argmin/blob/master/examples/particleswarm.rs)
20///
21/// # References:
22///
23/// TODO
24#[derive(Serialize, Deserialize)]
25pub struct ParticleSwarm<P, F> {
26    particles: Vec<Particle<P, F>>,
27    best_position: P,
28    best_cost: F,
29
30    // Weights for particle updates
31    weight_momentum: F,
32    weight_particle: F,
33    weight_swarm: F,
34
35    search_region: (P, P),
36    num_particles: usize,
37}
38
39impl<P, F> ParticleSwarm<P, F>
40where
41    P: Position<F> + DeserializeOwned + Serialize,
42    F: ArgminFloat,
43{
44    /// Constructor
45    ///
46    /// Parameters:
47    ///
48    /// * `cost_function`: cost function
49    /// * `init_temp`: initial temperature
50    pub fn new(
51        search_region: (P, P),
52        num_particles: usize,
53        weight_momentum: F,
54        weight_particle: F,
55        weight_swarm: F,
56    ) -> Result<Self, Error> {
57        let particle_swarm = ParticleSwarm {
58            particles: vec![],
59            best_position: P::rand_from_range(
60                // FIXME: random smart?
61                &search_region.0,
62                &search_region.1,
63            ),
64            best_cost: F::infinity(),
65            weight_momentum,
66            weight_particle,
67            weight_swarm,
68            search_region,
69            num_particles,
70        };
71
72        Ok(particle_swarm)
73    }
74
75    fn initialize_particles<O: ArgminOp<Param = P, Output = F, Float = F>>(
76        &mut self,
77        op: &mut OpWrapper<O>,
78    ) {
79        self.particles = (0..self.num_particles)
80            .map(|_| self.initialize_particle(op))
81            .collect();
82
83        self.best_position = self.get_best_position();
84        self.best_cost = op.apply(&self.best_position).unwrap();
85        // TODO unwrap evil
86    }
87
88    fn initialize_particle<O: ArgminOp<Param = P, Output = F, Float = F>>(
89        &mut self,
90        op: &mut OpWrapper<O>,
91    ) -> Particle<P, F> {
92        let (min, max) = &self.search_region;
93        let delta = max.sub(min);
94        let delta_neg = delta.mul(&F::from_f64(-1.0).unwrap());
95
96        let initial_position = O::Param::rand_from_range(min, max);
97        let initial_cost = op.apply(&initial_position).unwrap(); // FIXME do not unwrap
98
99        Particle {
100            position: initial_position.clone(),
101            velocity: O::Param::rand_from_range(&delta_neg, &delta),
102            cost: initial_cost,
103            best_position: initial_position,
104            best_cost: initial_cost,
105        }
106    }
107
108    fn get_best_position(&self) -> P {
109        let mut best: Option<(&P, F)> = None;
110
111        for p in &self.particles {
112            match best {
113                Some(best_sofar) => {
114                    if p.cost < best_sofar.1 {
115                        best = Some((&p.position, p.cost))
116                    }
117                }
118                None => best = Some((&p.position, p.cost)),
119            }
120        }
121
122        match best {
123            Some(best_sofar) => best_sofar.0.clone(),
124            None => panic!("Particles not initialized"),
125        }
126    }
127}
128
129impl<O, P, F> Solver<O> for ParticleSwarm<P, F>
130where
131    O: ArgminOp<Output = F, Param = P, Float = F>,
132    O::Param: Position<F> + DeserializeOwned + Serialize,
133    O::Hessian: Clone + Default,
134    F: ArgminFloat,
135{
136    const NAME: &'static str = "Particle Swarm Optimization";
137
138    fn init(
139        &mut self,
140        _op: &mut OpWrapper<O>,
141        _state: &IterState<O>,
142    ) -> Result<Option<ArgminIterData<O>>, Error> {
143        self.initialize_particles(_op);
144
145        Ok(None)
146    }
147
148    /// Perform one iteration of algorithm
149    fn next_iter(
150        &mut self,
151        _op: &mut OpWrapper<O>,
152        _state: &IterState<O>,
153    ) -> Result<ArgminIterData<O>, Error> {
154        let zero = O::Param::zero_like(&self.best_position);
155
156        for p in self.particles.iter_mut() {
157            // New velocity is composed of
158            // 1) previous velocity (momentum),
159            // 2) motion toward particle optimum and
160            // 3) motion toward global optimum.
161
162            // ad 1)
163            let momentum = p.velocity.mul(&self.weight_momentum);
164
165            // ad 2)
166            let to_optimum = p.best_position.sub(&p.position);
167            let pull_to_optimum = O::Param::rand_from_range(&zero, &to_optimum);
168            let pull_to_optimum = pull_to_optimum.mul(&self.weight_particle);
169
170            // ad 3)
171            let to_global_optimum = self.best_position.sub(&p.position);
172            let pull_to_global_optimum =
173                O::Param::rand_from_range(&zero, &to_global_optimum).mul(&self.weight_swarm);
174
175            p.velocity = momentum.add(&pull_to_optimum).add(&pull_to_global_optimum);
176            let new_position = p.position.add(&p.velocity);
177
178            // Limit to search window:
179            p.position = O::Param::min(
180                &O::Param::max(&new_position, &self.search_region.0),
181                &self.search_region.1,
182            );
183
184            p.cost = _op.apply(&p.position)?;
185            if p.cost < p.best_cost {
186                p.best_position = p.position.clone();
187                p.best_cost = p.cost;
188
189                if p.cost < self.best_cost {
190                    self.best_position = p.position.clone();
191                    self.best_cost = p.cost;
192                }
193            }
194        }
195
196        // Store particles as population
197        let population = self
198            .particles
199            .iter()
200            .map(|particle| (particle.position.clone(), particle.cost))
201            .collect();
202
203        let out = ArgminIterData::new()
204            .param(self.best_position.clone())
205            .cost(self.best_cost)
206            .population(population)
207            .kv(make_kv!(
208                "particles" => &self.particles;
209            ));
210
211        Ok(out)
212    }
213}
214
215/// Position
216pub trait Position<F: ArgminFloat>:
217    Clone
218    + Default
219    + ArgminAdd<Self, Self>
220    + ArgminSub<Self, Self>
221    + ArgminMul<F, Self>
222    + ArgminZeroLike
223    + ArgminRandom
224    + ArgminMinMax
225    + std::fmt::Debug
226{
227}
228impl<T, F: ArgminFloat> Position<F> for T
229where
230    T: Clone
231        + Default
232        + ArgminAdd<Self, Self>
233        + ArgminSub<Self, Self>
234        + ArgminMul<F, Self>
235        + ArgminZeroLike
236        + ArgminRandom
237        + ArgminMinMax
238        + std::fmt::Debug,
239    F: ArgminFloat,
240{
241}
242
243// trait_bound!(Position<F>
244// ; Clone
245// , Default
246// , ArgminAdd<Self, Self>
247// , ArgminSub<Self, Self>
248// , ArgminMul<F, Self>
249// , ArgminZeroLike
250// , ArgminRandom
251// , ArgminMinMax
252// , std::fmt::Debug
253// );
254
255/// A single particle
256#[derive(Clone, Serialize, Deserialize, Debug)]
257pub struct Particle<T, F> {
258    /// Position of particle
259    pub position: T,
260    /// Velocity of particle
261    velocity: T,
262    /// Cost of particle
263    pub cost: F,
264    /// Best position of particle so far
265    best_position: T,
266    /// Best cost of particle so far
267    best_cost: F,
268}
argmin/solver/particleswarm/mod.rs

argmin/solver/particleswarm/
mod.rs
