multivariate_optimization/optimize.rs
1//! Multivariate optimization (through estimation of distribution).
2//!
3//! # Example
4//!
5//! ```
6//! use multivariate_optimization::optimize::*;
7//! use multivariate_optimization::testfuncs::rastrigin;
8//! use rand::{Rng, thread_rng};
9//! let mut rng = thread_rng();
10//!
11//! const DIM: usize = 3;
12//! let search_space: Vec<SearchRange> = (0..DIM).map(|_| {
13//! SearchRange::Finite {
14//! low: rng.gen_range(-7.68..=-2.56),
15//! high: rng.gen_range(2.56..=7.68),
16//! }
17//! }).collect();
18//!
19//! const POPULATION: usize = 1000;
20//! const MAX_ITERATIONS: usize = 1000;
21//! let mut solver = Solver::new(search_space, |params| {
22//! let cost = rastrigin(¶ms);
23//! BasicSpecimen { params, cost }
24//! });
25//! solver.set_speed_factor(0.5);
26//!
27//! let initial_specimens = solver.random_specimens(POPULATION);
28//! solver.extend_specimens(initial_specimens);
29//! for iter in 0..MAX_ITERATIONS {
30//! let specimens = solver.specimens();
31//! println!(
32//! "{} {}",
33//! specimens[0].cost,
34//! specimens[specimens.len()-1].cost,
35//! );
36//! if solver.converged() {
37//! break;
38//! }
39//! let new_specimens = solver.recombined_specimens(POPULATION, 0.0);
40//! solver.replace_worst_specimens(new_specimens);
41//! }
42//! let specimen = solver.into_specimen();
43//! assert_eq!(specimen.cost, 0.0);
44//! ```
45//!
46//! See also [`Solver`].
47
48use crate::distributions::{MultivarNormDist, NormDist};
49use crate::splitter::Splitter;
50use crate::triangular::Triangular;
51
52use futures::stream::{FuturesOrdered, StreamExt};
53use rand::distributions::{Distribution, Standard, Uniform};
54use rand::seq::SliceRandom;
55use rand::Rng;
56use rayon::prelude::*;
57
58use std::cmp::Ordering;
59use std::future::Future;
60use std::ops::RangeInclusive;
61
62/// Search range for a single dimension.
63///
64/// `SearchRange` describes a search range for a single dimension and
65/// `Vec<SearchRange>` describes a multidimensional search space (which can be
66/// passed to [`Solver::new`]).
67#[derive(Copy, Clone, Debug)]
68pub enum SearchRange {
69 /// Finite search range.
70 Finite {
71 /// Lower bound.
72 low: f64,
73 /// Upper bound.
74 high: f64,
75 },
76 /// Infinite search range.
77 Infinite {
78 /// Initial average value for starting search.
79 average: f64,
80 /// Initial standard deviation for starting search.
81 stddev: f64,
82 },
83}
84
85impl From<RangeInclusive<f64>> for SearchRange {
86 fn from(range: RangeInclusive<f64>) -> Self {
87 SearchRange::Finite {
88 low: *range.start(),
89 high: *range.end(),
90 }
91 }
92}
93
94/// Enum with one-dimensional distribution used for initial search.
95#[derive(Clone, Debug)]
96enum SearchDist {
97 Finite(Uniform<f64>),
98 Infinite(NormDist<f64>),
99}
100
101impl From<SearchRange> for SearchDist {
102 /// Create one-dimensional distribution for given SearchRange.
103 fn from(search_range: SearchRange) -> SearchDist {
104 match search_range {
105 SearchRange::Finite { low, high } => {
106 SearchDist::Finite(Uniform::new_inclusive(low, high))
107 }
108 SearchRange::Infinite { average, stddev } => {
109 SearchDist::Infinite(NormDist::new(average, stddev))
110 }
111 }
112 }
113}
114
115impl Distribution<f64> for SearchDist {
116 fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> f64 {
117 match self {
118 SearchDist::Finite(dist) => dist.sample(rng),
119 SearchDist::Infinite(dist) => dist.sample(rng),
120 }
121 }
122}
123
124/// Specimen in evolutionary process.
125///
126/// Two methods are required to be implemented:
127/// * The [`params`](Specimen::params) method needs to return the parameters
128/// originally used to create the specimen (see `constructor` argument to
129/// function [`Solver::new`]).
130/// * The [`cmp_cost`](Specimen::cmp_cost) method compares specimen by their
131/// fitness (lower cost values are better).
132///
133/// For most cases, the module provides a [`BasicSpecimen`] which is an
134/// implementation only storing the `params` and `cost` values.
135pub trait Specimen {
136 /// Parameters used for creation.
137 fn params(&self) -> &[f64];
138 /// Compare specimen's fitness ([`Less`] means better).
139 ///
140 /// [`Less`]: Ordering::Less
141 fn cmp_cost(&self, other: &Self) -> Ordering;
142 /// Euclidean distance between two specimens' parameters
143 fn params_dist(&self, other: &Self) -> f64 {
144 self.params()
145 .iter()
146 .copied()
147 .zip(other.params().iter().copied())
148 .map(|(a, b)| (a - b).powf(2.0))
149 .sum()
150 }
151}
152
153/// Most simple implementation of a [`Specimen`].
154#[derive(Clone, Debug)]
155pub struct BasicSpecimen {
156 /// Coefficients used to create specimen
157 pub params: Vec<f64>,
158 /// Fitness (smaller values are better)
159 pub cost: f64,
160}
161
162impl Specimen for BasicSpecimen {
163 fn params(&self) -> &[f64] {
164 &self.params
165 }
166 fn cmp_cost(&self, other: &Self) -> Ordering {
167 let a = self.cost;
168 let b = other.cost;
169 a.partial_cmp(&b)
170 .unwrap_or_else(|| a.is_nan().cmp(&b.is_nan()).then(Ordering::Equal))
171 }
172}
173
174/// Parallel solver for multidimensional problems.
175///
176/// Usual workflow:
177///
178/// * [`Solver::new`]
179/// * [`Solver::set_speed_factor`]
180/// * Pass result of [`Solver::random_specimens`] to [`Solver::extend_specimens`]
181/// * In a loop:
182/// * Inspect first element (or more elements) of [`Solver::specimens`],
183/// e.g. for a break condition
184/// * Pass result of [`Solver::recombined_specimens`] to [`Solver::replace_worst_specimens`]
185/// * Extract best specimen, e.g. using [`Solver::into_specimen`]
186///
187/// See [module level documentation] for a code example.
188///
189/// [module level documentation]: self
190#[derive(Clone, Debug)]
191pub struct Solver<S, C> {
192 search_space: Vec<SearchRange>,
193 search_dists: Vec<SearchDist>,
194 constructor: C,
195 division_count: usize,
196 min_population: usize,
197 is_sorted: bool,
198 specimens: Vec<S>,
199}
200
201impl<S, C> Solver<S, C> {
202 /// Create `Solver` for search space and [`Specimen`] `constructor` closure.
203 ///
204 /// The closure takes a [`Vec<f64>`] as argument, which contains the
205 /// coefficients/parameters, and it returns an [`S: Specimen`].
206 /// See [module level documentation] for a code example.
207 ///
208 /// Alternatively, [`Future`]s or [`Result`]s may be returned by the
209 /// closure. Note that when the closure returns `Future`s, then the methods
210 /// [`extend_specimens_async`] and [`replace_worst_specimens_async`] must
211 /// be used instead of their synchronous equivalents.
212 ///
213 /// [`S: Specimen`]: Specimen
214 /// [module level documentation]: self
215 /// [`extend_specimens_async`]: Self::extend_specimens_async
216 /// [`replace_worst_specimens_async`]: Self::replace_worst_specimens_async
217 pub fn new<T>(search_space: Vec<SearchRange>, constructor: C) -> Self
218 where
219 C: Fn(Vec<f64>) -> T + Sync,
220 {
221 let search_dists: Vec<SearchDist> = search_space
222 .iter()
223 .copied()
224 .map(|search_range| SearchDist::from(search_range))
225 .collect();
226 let mut solver = Solver {
227 search_space,
228 search_dists,
229 constructor,
230 division_count: Default::default(),
231 min_population: Default::default(),
232 is_sorted: true,
233 specimens: vec![],
234 };
235 solver.set_division_count(1);
236 solver
237 }
238 /// Dimensionality of search space.
239 pub fn dim(&self) -> usize {
240 self.search_space.len()
241 }
242 /// Simplify calculation by dividing dimensions into a given number of
243 /// groups.
244 ///
245 /// Divides dimensions into (almost) equally sized groups when calculating
246 /// covariances to reduce computation time. The number of groups is given
247 /// as an integer.
248 ///
249 /// See [`Solver::set_speed_factor`] for a high-level interface.
250 pub fn set_division_count(&mut self, mut division_count: usize) {
251 let dim = self.dim();
252 if dim == 0 {
253 self.division_count = 1;
254 self.min_population = 0;
255 } else {
256 if division_count < 1 {
257 division_count = 1;
258 } else if division_count > dim {
259 division_count = dim;
260 }
261 self.division_count = division_count;
262 self.min_population = (dim - 1) / division_count + 1;
263 }
264 }
265 /// Simplify calculation by dividing dimensions according to speed factor.
266 ///
267 /// This method is a high-level interface for
268 /// [`Solver::set_division_count`].
269 ///
270 /// Uses a `speed_factor` to divide dimensions into (almost) equally sized
271 /// groups when calculating covariances to reduce computation time. Higher
272 /// speed factors generally result in faster (but possibly inaccurate)
273 /// calculation.
274 ///
275 /// Speed factors range from `0.0` to `1.0`, however, a factor greater than
276 /// `0.75` is not recommended due to the introduced overhead.
277 pub fn set_speed_factor(&mut self, speed_factor: f64) {
278 assert!(
279 speed_factor >= 0.0 && speed_factor <= 1.0,
280 "speed_factor must be between 0.0 and 1.0"
281 );
282 self.set_division_count((self.dim() as f64).powf(speed_factor).round() as usize);
283 }
284 /// Simply calculation by dividing dimensions into groups of specified
285 /// maximum size.
286 ///
287 /// This method is an alternative interface for
288 /// [`Solver::set_division_count`] where the maximum size of each division
289 /// is specified.
290 /// See [`Solver::set_speed_factor`] for a high-level interface.
291 pub fn set_max_division_size(&mut self, max_division_size: usize) {
292 self.set_division_count((self.dim() as f64 / max_division_size as f64).ceil() as usize);
293 }
294 /// Number of groups into which dimensions are split when calculating
295 /// covariances.
296 pub fn division_count(&self) -> usize {
297 self.division_count
298 }
299 /// Minimum required population.
300 ///
301 /// The number depends on the number of dimensions and the number of groups
302 /// into which the dimensions are split when calculating covariances.
303 pub fn min_population(&self) -> usize {
304 self.min_population
305 }
306}
307
308impl<S, C> Solver<S, C>
309where
310 S: Specimen + Send,
311{
312 /// Ensures that speciments are sorted based on cost (best first).
313 fn sort(&mut self) {
314 if !self.is_sorted {
315 self.specimens.par_sort_by(S::cmp_cost);
316 self.is_sorted = true;
317 }
318 }
319 /// Add specimens to population.
320 pub fn extend_specimens<I: IntoIterator<Item = S>>(&mut self, iter: I) {
321 self.is_sorted = false;
322 self.specimens.extend(iter);
323 }
324 /// Replace worst specimens in population.
325 pub fn replace_worst_specimens<I: IntoIterator<Item = S>>(&mut self, iter: I) {
326 let count = self.specimens.len();
327 self.extend_specimens(iter);
328 self.truncate_specimens(count);
329 }
330 /// Add specimens to population asynchronously.
331 pub async fn extend_specimens_async<F, I>(&mut self, iter: I)
332 where
333 F: Future<Output = S> + Send,
334 I: IntoIterator<Item = F>,
335 {
336 let new_specimens = FuturesOrdered::from_iter(iter);
337 self.specimens.reserve(new_specimens.len());
338 self.is_sorted = false;
339 new_specimens
340 .for_each(|specimen| {
341 self.specimens.push(specimen);
342 async { () }
343 })
344 .await;
345 }
346 /// Replace worst specimens in population asynchronously.
347 pub async fn replace_worst_specimens_async<F, I>(&mut self, iter: I)
348 where
349 F: Future<Output = S> + Send,
350 I: IntoIterator<Item = F>,
351 {
352 let count = self.specimens.len();
353 self.extend_specimens_async(iter).await;
354 self.truncate_specimens(count);
355 }
356 /// Truncate population of specimens to given `count`
357 /// (drops worst fitting specimens).
358 pub fn truncate_specimens(&mut self, count: usize) {
359 self.sort();
360 self.specimens.truncate(count);
361 }
362 /// Return true if specimens have converged.
363 ///
364 /// Note that this method only returns true if the [cost] of all specimens
365 /// is exactly equal. Additional, more sensitive termination conditions may
366 /// be required in practice in order to avoid endless optimization.
367 ///
368 /// [cost]: Specimen::cmp_cost
369 pub fn converged(&mut self) -> bool {
370 let len = self.specimens.len();
371 if len == 0 {
372 true
373 } else {
374 self.sort();
375 self.specimens[0].cmp_cost(&self.specimens[len - 1]) == Ordering::Equal
376 }
377 }
378 /// Sorted population of [`Specimen`]s as shared slice (best first).
379 pub fn specimens(&mut self) -> &[S] {
380 self.sort();
381 &self.specimens
382 }
383 /// Sorted population of [`Specimen`]s as mutable [`Vec`] (best first).
384 ///
385 /// The `Vec` may be modified and doesn't need to be (re-)sorted by the
386 /// caller after modifying.
387 pub fn specimens_mut(&mut self) -> &mut Vec<S> {
388 self.sort();
389 self.is_sorted = false;
390 &mut self.specimens
391 }
392 /// Consume [`Solver`] and return all [`Specimen`]s, ordered by fitness
393 /// (best first).
394 pub fn into_specimens(mut self) -> Vec<S> {
395 self.sort();
396 self.specimens
397 }
398 /// Consume [`Solver`] and return best [`Specimen`].
399 pub fn into_specimen(self) -> S {
400 self.into_specimens()
401 .into_iter()
402 .next()
403 .expect("solver contains no specimen")
404 }
405}
406
407impl<S, C, T> Solver<S, C>
408where
409 S: Specimen + Send + Sync,
410 C: Fn(Vec<f64>) -> T + Sync,
411 T: Send,
412{
413 /// Create random specimen (optionally async if `T` is a [`Future`]).
414 fn random_specimen<R>(&self, rng: &mut R) -> T
415 where
416 R: Rng + ?Sized,
417 {
418 (self.constructor)(
419 self.search_dists
420 .iter()
421 .map(|dist| dist.sample(rng))
422 .collect::<Vec<f64>>(),
423 )
424 }
425 /// Create random specimens
426 ///
427 /// Note that depending on the return type of the closure passed to
428 /// [`Solver::new`], a [`Vec`] of [`Future`]s or [`Result`]s may be
429 /// returned.
430 pub fn random_specimens(&self, count: usize) -> Vec<T> {
431 (0..count)
432 .into_par_iter()
433 .map_init(|| rand::thread_rng(), |rng, _| self.random_specimen(rng))
434 .collect()
435 }
436 /// Create recombined specimens.
437 ///
438 /// Note that depending on the return type of the closure passed to
439 /// [`Solver::new`], a [`Vec`] of [`Future`]s or [`Result`]s may be
440 /// returned.
441 ///
442 /// Setting the `local_factor` to a value greater than `0.0` (but smaller
443 /// than `1.0`) selects a particular specimen with a correspondingly
444 /// proportional chance to be modified. This allows performing more
445 /// localized searches. A reasonable value seems to be
446 /// `0.01 / self.dim() as f64`.
447 pub fn recombined_specimens(&mut self, children_count: usize, local_factor: f64) -> Vec<T> {
448 self.sort();
449 let total_count = self.specimens.len();
450 let total_weight = total_count as f64;
451 let splitter = Splitter::new(&mut rand::thread_rng(), self.dim(), self.division_count());
452 let sub_averages = splitter
453 .groups()
454 .par_iter()
455 .map(|group| {
456 (0..group.len())
457 .into_par_iter()
458 .map(|i| {
459 let i_orig = group[i];
460 self.specimens
461 .par_iter()
462 .map(|specimen| specimen.params()[i_orig])
463 .sum::<f64>()
464 / total_weight
465 })
466 .collect::<Vec<_>>()
467 })
468 .collect::<Vec<_>>(); // TODO: use boxed slice when supported by rayon
469 let sub_dists = splitter
470 .groups()
471 .par_iter()
472 .zip(sub_averages.into_par_iter())
473 .map(|(group, averages)| {
474 let covariances = Triangular::<f64>::par_new(group.len(), |(i, j)| {
475 let i_orig = group[i];
476 let j_orig = group[j];
477 self.specimens
478 .iter()
479 .map(|specimen| {
480 let a = specimen.params()[i_orig] - averages[i];
481 let b = specimen.params()[j_orig] - averages[j];
482 a * b
483 })
484 .sum::<f64>()
485 / total_weight
486 });
487 MultivarNormDist::new(averages, covariances)
488 })
489 .collect::<Vec<_>>(); // TODO: use boxed slice when supported by rayon
490 let local_exp = if local_factor > 0.0 {
491 1.0 / local_factor
492 } else {
493 f64::INFINITY
494 };
495 (0..children_count)
496 .into_par_iter()
497 .map_init(
498 || rand::thread_rng(),
499 |rng, _| {
500 let param_groups_iters: Box<[_]> = sub_dists
501 .iter()
502 .map(|dist| dist.sample(rng).into_iter())
503 .collect();
504 let mut params: Vec<_> = splitter.merge(param_groups_iters).collect();
505 let specimen = self.specimens.choose(rng).unwrap();
506 let parent_params = specimen.params();
507 let factor1: f64 = Standard.sample(rng);
508 let factor1 = 2.0 * factor1.powf(local_exp);
509 let factor2: f64 = 1.0 - factor1;
510 for i in 0..params.len() {
511 params[i] = factor1 * parent_params[i] + factor2 * params[i];
512 if let SearchRange::Finite { low, high } = self.search_space[i] {
513 if !(low..=high).contains(¶ms[i]) {
514 return self.random_specimen(rng);
515 }
516 }
517 }
518 (self.constructor)(params)
519 },
520 )
521 .collect()
522 }
523}
524
525#[cfg(test)]
526mod tests {
527 use super::{BasicSpecimen, SearchRange, Solver, Specimen as _};
528 use rand::{thread_rng, Rng};
529 #[test]
530 fn test_solver() {
531 let mut rng = thread_rng();
532 const PARAMCNT: usize = 3;
533 let ranges = vec![-1.0..=1.0; PARAMCNT];
534 let search_space: Vec<SearchRange> = ranges.iter().cloned().map(Into::into).collect();
535 let goals: Vec<f64> = ranges
536 .iter()
537 .cloned()
538 .map(|range| rng.gen_range(range) * 0.75)
539 .collect();
540 let mut solver = Solver::new(search_space, |params: Vec<f64>| {
541 let mut cost: f64 = 0.0;
542 for (param, goal) in params.iter().zip(goals.iter()) {
543 cost += (param - goal) * (param - goal);
544 }
545 BasicSpecimen { params, cost }
546 });
547 let initial_specimens = solver.random_specimens(200);
548 solver.extend_specimens(initial_specimens);
549 for _ in 0..1000 {
550 let new_specimens = solver.recombined_specimens(10, 0.0);
551 solver.replace_worst_specimens(new_specimens);
552 }
553 for (param, goal) in solver
554 .specimens
555 .first()
556 .unwrap()
557 .params()
558 .iter()
559 .zip(goals.iter())
560 {
561 assert!((param - goal).abs() < 0.01);
562 }
563 }
564}