argmin/solver/simulatedannealing/
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//! * [Simulated Annealing](struct.SimulatedAnnealing.html)
9//!
10//! # References
11//!
12//! [0] [Wikipedia](https://en.wikipedia.org/wiki/Simulated_annealing)
13//!
14//! [1] S Kirkpatrick, CD Gelatt Jr, MP Vecchi. (1983). "Optimization by Simulated Annealing".
15//! Science 13 May 1983, Vol. 220, Issue 4598, pp. 671-680
16//! DOI: 10.1126/science.220.4598.671
17
18use crate::prelude::*;
19use rand::prelude::*;
20use rand_xorshift::XorShiftRng;
21use serde::{Deserialize, Serialize};
22
23/// Temperature functions for Simulated Annealing.
24///
25/// Given the initial temperature `t_init` and the iteration number `i`, the current temperature
26/// `t_i` is given as follows:
27///
28/// * `SATempFunc::TemperatureFast`: `t_i = t_init / i`
29/// * `SATempFunc::Boltzmann`: `t_i = t_init / ln(i)`
30/// * `SATempFunc::Exponential`: `t_i = t_init * 0.95^i`
31#[derive(Clone, Copy, Serialize, Deserialize, Debug)]
32pub enum SATempFunc<F> {
33    /// `t_i = t_init / i`
34    TemperatureFast,
35    /// `t_i = t_init / ln(i)`
36    Boltzmann,
37    /// `t_i = t_init * x^i`
38    Exponential(F),
39    // /// User-provided temperature function. The first parameter must be the current temperature and
40    // /// the second parameter must be the iteration number.
41    // Custom(Box<Fn(f64, u64) -> f64>),
42}
43
44impl<F> std::default::Default for SATempFunc<F> {
45    fn default() -> Self {
46        SATempFunc::Boltzmann
47    }
48}
49
50/// Simulated Annealing
51///
52/// [Example](https://github.com/argmin-rs/argmin/blob/master/examples/simulatedannealing.rs)
53///
54/// # References
55///
56/// [0] [Wikipedia](https://en.wikipedia.org/wiki/Simulated_annealing)
57///
58/// [1] S Kirkpatrick, CD Gelatt Jr, MP Vecchi. (1983). "Optimization by Simulated Annealing".
59/// Science 13 May 1983, Vol. 220, Issue 4598, pp. 671-680
60/// DOI: 10.1126/science.220.4598.671  
61#[derive(Clone, Serialize, Deserialize)]
62pub struct SimulatedAnnealing<F> {
63    /// Initial temperature
64    init_temp: F,
65    /// which temperature function?
66    temp_func: SATempFunc<F>,
67    /// Number of iterations used for the caluclation of temperature. This is needed for
68    /// reannealing!
69    temp_iter: u64,
70    /// Iterations since the last accepted solution
71    stall_iter_accepted: u64,
72    /// Stop if stall_iter_accepted exceedes this number
73    stall_iter_accepted_limit: u64,
74    /// Iterations since the last best solution was found
75    stall_iter_best: u64,
76    /// Stop if stall_iter_best exceedes this number
77    stall_iter_best_limit: u64,
78    /// Reanneal after this number of iterations is reached
79    reanneal_fixed: u64,
80    /// Similar to `iter`, but will be reset to 0 when reannealing is performed
81    reanneal_iter_fixed: u64,
82    /// Reanneal after no accepted solution has been found for `reanneal_accepted` iterations
83    reanneal_accepted: u64,
84    /// Similar to `stall_iter_accepted`, but will be reset to 0 when reannealing  is performed
85    reanneal_iter_accepted: u64,
86    /// Reanneal after no new best solution has been found for `reanneal_best` iterations
87    reanneal_best: u64,
88    /// Similar to `stall_iter_best`, but will be reset to 0 when reannealing is performed
89    reanneal_iter_best: u64,
90    /// current temperature
91    cur_temp: F,
92    /// random number generator
93    rng: XorShiftRng,
94}
95
96impl<F> SimulatedAnnealing<F>
97where
98    F: ArgminFloat,
99{
100    /// Constructor
101    ///
102    /// Parameter:
103    ///
104    /// * `init_temp`: initial temperature
105    pub fn new(init_temp: F) -> Result<Self, Error> {
106        if init_temp <= F::from_f64(0.0).unwrap() {
107            Err(ArgminError::InvalidParameter {
108                text: "Initial temperature must be > 0.".to_string(),
109            }
110            .into())
111        } else {
112            Ok(SimulatedAnnealing {
113                init_temp,
114                temp_func: SATempFunc::TemperatureFast,
115                temp_iter: 0,
116                stall_iter_accepted: 0,
117                stall_iter_accepted_limit: std::u64::MAX,
118                stall_iter_best: 0,
119                stall_iter_best_limit: std::u64::MAX,
120                reanneal_fixed: std::u64::MAX,
121                reanneal_iter_fixed: 0,
122                reanneal_accepted: std::u64::MAX,
123                reanneal_iter_accepted: 0,
124                reanneal_best: std::u64::MAX,
125                reanneal_iter_best: 0,
126                cur_temp: init_temp,
127                rng: XorShiftRng::from_entropy(),
128            })
129        }
130    }
131
132    /// Set temperature function to one of the options in `SATempFunc`.
133    pub fn temp_func(mut self, temperature_func: SATempFunc<F>) -> Self {
134        self.temp_func = temperature_func;
135        self
136    }
137
138    /// The optimization stops after there has been no accepted solution after `iter` iterations
139    pub fn stall_accepted(mut self, iter: u64) -> Self {
140        self.stall_iter_accepted_limit = iter;
141        self
142    }
143
144    /// The optimization stops after there has been no new best solution after `iter` iterations
145    pub fn stall_best(mut self, iter: u64) -> Self {
146        self.stall_iter_best_limit = iter;
147        self
148    }
149
150    /// Start reannealing after `iter` iterations
151    pub fn reannealing_fixed(mut self, iter: u64) -> Self {
152        self.reanneal_fixed = iter;
153        self
154    }
155
156    /// Start reannealing after no accepted solution has been found for `iter` iterations
157    pub fn reannealing_accepted(mut self, iter: u64) -> Self {
158        self.reanneal_accepted = iter;
159        self
160    }
161
162    /// Start reannealing after no new best solution has been found for `iter` iterations
163    pub fn reannealing_best(mut self, iter: u64) -> Self {
164        self.reanneal_best = iter;
165        self
166    }
167
168    /// Update the temperature based on the current iteration number.
169    ///
170    /// Updates are performed based on specific update functions. See `SATempFunc` for details.
171    fn update_temperature(&mut self) {
172        self.cur_temp = match self.temp_func {
173            SATempFunc::TemperatureFast => {
174                self.init_temp / F::from_u64(self.temp_iter + 1).unwrap()
175            }
176            SATempFunc::Boltzmann => self.init_temp / F::from_u64(self.temp_iter + 1).unwrap().ln(),
177            SATempFunc::Exponential(x) => {
178                self.init_temp * x.powf(F::from_u64(self.temp_iter + 1).unwrap())
179            }
180        };
181    }
182
183    /// Perform reannealing
184    fn reanneal(&mut self) -> (bool, bool, bool) {
185        let out = (
186            self.reanneal_iter_fixed >= self.reanneal_fixed,
187            self.reanneal_iter_accepted >= self.reanneal_accepted,
188            self.reanneal_iter_best >= self.reanneal_best,
189        );
190        if out.0 || out.1 || out.2 {
191            self.reanneal_iter_fixed = 0;
192            self.reanneal_iter_accepted = 0;
193            self.reanneal_iter_best = 0;
194            self.cur_temp = self.init_temp;
195            self.temp_iter = 0;
196        }
197        out
198    }
199
200    /// Update the stall iter variables
201    fn update_stall_and_reanneal_iter(&mut self, accepted: bool, new_best: bool) {
202        self.stall_iter_accepted = if accepted {
203            0
204        } else {
205            self.stall_iter_accepted + 1
206        };
207
208        self.reanneal_iter_accepted = if accepted {
209            0
210        } else {
211            self.reanneal_iter_accepted + 1
212        };
213
214        self.stall_iter_best = if new_best {
215            0
216        } else {
217            self.stall_iter_best + 1
218        };
219
220        self.reanneal_iter_best = if new_best {
221            0
222        } else {
223            self.reanneal_iter_best + 1
224        };
225    }
226}
227
228impl<O, F> Solver<O> for SimulatedAnnealing<F>
229where
230    O: ArgminOp<Output = F, Float = F>,
231    F: ArgminFloat,
232{
233    const NAME: &'static str = "Simulated Annealing";
234    fn init(
235        &mut self,
236        _op: &mut OpWrapper<O>,
237        _state: &IterState<O>,
238    ) -> Result<Option<ArgminIterData<O>>, Error> {
239        Ok(Some(ArgminIterData::new().kv(make_kv!(
240            "initial_temperature" => self.init_temp;
241            "stall_iter_accepted_limit" => self.stall_iter_accepted_limit;
242            "stall_iter_best_limit" => self.stall_iter_best_limit;
243            "reanneal_fixed" => self.reanneal_fixed;
244            "reanneal_accepted" => self.reanneal_accepted;
245            "reanneal_best" => self.reanneal_best;
246        ))))
247    }
248
249    /// Perform one iteration of SA algorithm
250    fn next_iter(
251        &mut self,
252        op: &mut OpWrapper<O>,
253        state: &IterState<O>,
254    ) -> Result<ArgminIterData<O>, Error> {
255        // Careful: The order in here is *very* important, even if it may not seem so. Everything
256        // is linked to the iteration number, and getting things mixed up will lead to strange
257        // behaviour.
258
259        let prev_param = state.get_param();
260        let prev_cost = state.get_cost();
261
262        // Make a move
263        let new_param = op.modify(&prev_param, self.cur_temp)?;
264        // let new_param = op.modify(&prev_param, self.cur_temp)?;
265
266        // Evaluate cost function with new parameter vector
267        let new_cost = op.apply(&new_param)?;
268
269        // Acceptance function
270        //
271        // Decide whether new parameter vector should be accepted.
272        // If no, move on with old parameter vector.
273        //
274        // Any solution which satisfies `next_cost < prev_cost` will be accepted. Solutions worse
275        // than the previous one are accepted with a probability given as:
276        //
277        // `1 / (1 + exp((next_cost - prev_cost) / current_temperature))`,
278        //
279        // which will always be between 0 and 0.5.
280        let prob: f64 = self.rng.gen();
281        let prob = F::from_f64(prob).unwrap();
282        let accepted = (new_cost < state.get_prev_cost())
283            || (F::from_f64(1.0).unwrap()
284                / (F::from_f64(1.0).unwrap()
285                    + ((new_cost - state.get_prev_cost()) / self.cur_temp).exp())
286                > prob);
287
288        // Update stall iter variables
289        self.update_stall_and_reanneal_iter(accepted, new_cost <= state.get_best_cost());
290
291        let (r_fixed, r_accepted, r_best) = self.reanneal();
292
293        // Update temperature for next iteration.
294        self.temp_iter += 1;
295        // Todo: this variable may not be necessary (temp_iter does the same?)
296        self.reanneal_iter_fixed += 1;
297
298        self.update_temperature();
299
300        Ok(if accepted {
301            ArgminIterData::new().param(new_param).cost(new_cost)
302        } else {
303            ArgminIterData::new().param(prev_param).cost(prev_cost)
304        }
305        .kv(make_kv!(
306            "t" => self.cur_temp;
307            "new_be" => new_cost <= state.get_best_cost();
308            "acc" => accepted;
309            "st_i_be" => self.stall_iter_best;
310            "st_i_ac" => self.stall_iter_accepted;
311            "ra_i_fi" => self.reanneal_iter_fixed;
312            "ra_i_be" => self.reanneal_iter_best;
313            "ra_i_ac" => self.reanneal_iter_accepted;
314            "ra_fi" => r_fixed;
315            "ra_be" => r_best;
316            "ra_ac" => r_accepted;
317        )))
318    }
319
320    fn terminate(&mut self, _state: &IterState<O>) -> TerminationReason {
321        if self.stall_iter_accepted > self.stall_iter_accepted_limit {
322            return TerminationReason::AcceptedStallIterExceeded;
323        }
324        if self.stall_iter_best > self.stall_iter_best_limit {
325            return TerminationReason::BestStallIterExceeded;
326        }
327        TerminationReason::NotTerminated
328    }
329}
330
331#[cfg(test)]
332mod tests {
333    use super::*;
334    use crate::test_trait_impl;
335
336    test_trait_impl!(sa, SimulatedAnnealing<f64>);
337}
argmin/solver/simulatedannealing/mod.rs

argmin/solver/simulatedannealing/
mod.rs
