aphreco/optimizer/
genetic_algorithm.rs

1use super::base::{ConcreteOptimizer, OptOptions};
2use super::result::OptResult;
3
4use crate::model::OptModelTrait;
5use crate::objective::Objective;
6
7use ndarray::Array1;
8use rand::distributions::{Distribution, WeightedIndex};
9use rand::rngs::ThreadRng;
10use rand::{thread_rng, Rng};
11use std::thread;
12
13type Individual = (f64, Array1<f64>);
14type Population = Vec<Individual>;
15
16#[allow(dead_code)]
17pub struct GeneticAlgorithm {
18  max_gen: u64,
19  n_pop: usize,
20  n_elite: usize,
21  mutation_rate: f64,
22  len_x: usize,
23  verbose: bool,
24}
25
26impl ConcreteOptimizer for GeneticAlgorithm {
27  fn new(len_x: usize, options: &OptOptions) -> Self {
28    let (max_gen, n_pop, mutation_rate, verbose) = match options {
29      OptOptions::Default => (100, 10, 0.8, false),
30
31      OptOptions::GeneticAlgorithm {
32        max_gen,
33        n_pop,
34        mutation_rate,
35        verbose,
36      } => (*max_gen, *n_pop, *mutation_rate, *verbose),
37
38      _ => panic!("Invalid OptOptions variant."),
39    };
40
41    let n_elite = if n_pop / 10 == 0 { 1 } else { n_pop / 10 };
42
43    Self {
44      max_gen,
45      n_pop,
46      n_elite,
47      mutation_rate,
48      len_x,
49      verbose,
50    }
51  }
52
53  fn run<M, const LEN_Y: usize, const LEN_P: usize, const LEN_B: usize, const LEN_X: usize>(
54    &self,
55    objective: &mut Objective<M, LEN_Y, LEN_P, LEN_B, LEN_X>,
56  ) -> OptResult
57  where
58    M: OptModelTrait<LEN_Y, LEN_P, LEN_B, LEN_X>,
59  {
60    let mut fcall: u64 = 0;
61    let mut rng = thread_rng();
62
63    // bounds
64    let log10_bounds = self.get_bounds(objective);
65
66    // make initial population
67    let mut pop = self.make_initial_pop(&log10_bounds, &mut rng);
68    let mut next_pop = pop.clone();
69
70    for n_gen in 0..self.max_gen {
71      // vector for join-handles
72      let mut handles = Vec::new();
73
74      for ind in pop.iter() {
75        let thread_ind = ind.clone();
76        let mut thread_objective = objective.clone();
77
78        // ===== FORK =====
79        let handle = thread::spawn(move || {
80          // if individual is an elite, the fitness has already been
81          // evaluated in the previous generation.
82          if thread_ind.0 == f64::INFINITY {
83            let thread_f = thread_objective.obj(&thread_ind.1);
84            thread_f
85          } else {
86            thread_ind.0
87          }
88        });
89        // ================
90
91        fcall += 1;
92        handles.push(handle);
93      }
94
95      // ===== JOIN =====
96      for (ind, handle) in pop.iter_mut().zip(handles) {
97        ind.0 = handle.join().unwrap();
98      }
99      // ================
100
101      // sort individuals in descending order (fmax objective)
102      pop.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap());
103
104      // print
105      if self.verbose {
106        println!("{:5}:   f:{:.4e}   x:{:10.8}", n_gen, pop[0].0, pop[0].1);
107      }
108
109      // choose elite
110      for i in 0..self.n_elite {
111        next_pop[i] = pop[i].clone();
112      }
113
114      // crossover
115      for i in self.n_elite..self.n_pop {
116        // current WeightedIndex
117        let dist = self.weight_distribution(&pop);
118        let (index_p1, index_p2) = self.select(dist, &mut rng);
119        next_pop[i] = self.crossover(&pop[index_p1], &pop[index_p2], &mut rng);
120      }
121
122      // mutate
123      // TODO: in parallel
124      for i in 1..self.n_pop {
125        self.mutate(&mut next_pop[i], &log10_bounds, &mut rng);
126      }
127
128      // alternate
129      pop = next_pop.clone();
130    }
131
132    println!("Finished. fcall = {}", fcall);
133    OptResult::new(pop[0].1.clone(), objective.x_index.clone(), pop[0].0)
134  }
135}
136
137impl GeneticAlgorithm {
138  fn get_bounds<M, const LEN_Y: usize, const LEN_P: usize, const LEN_B: usize, const LEN_X: usize>(
139    &self,
140    objective: &Objective<M, LEN_Y, LEN_P, LEN_B, LEN_X>,
141  ) -> Vec<(f64, f64)>
142  where
143    M: OptModelTrait<LEN_Y, LEN_P, LEN_B, LEN_X>,
144  {
145    let mut log10_bounds: Vec<(f64, f64)> = Vec::new();
146    let x_bounds = objective
147      .x_bounds
148      .as_ref()
149      .expect("please define lower and upper bounds.");
150
151    for &(lb, ub) in x_bounds.iter() {
152      log10_bounds.push((f64::log10(lb), f64::log10(ub)));
153    }
154
155    log10_bounds
156  }
157
158  fn make_initial_pop(&self, log10_bounds: &Vec<(f64, f64)>, rng: &mut ThreadRng) -> Population {
159    let mut pop = Vec::new();
160    let mut ind: Individual = (f64::INFINITY, Array1::zeros(self.len_x));
161
162    for _ in 0..self.n_pop {
163      for (i, &(log10_lb, log10_ub)) in log10_bounds.iter().enumerate() {
164        ind.1[i] = 10f64.powf(rng.gen_range(log10_lb..log10_ub));
165      }
166      pop.push(ind.clone());
167    }
168    pop
169  }
170
171  fn weight_distribution(&self, pop: &Population) -> WeightedIndex<f64> {
172    // make a roulette (weights) for selection
173    let f_worst = pop[self.n_pop - 1].0;
174    let f_best = pop[0].0;
175
176    let weights: Vec<f64> = pop
177      .iter()
178      .map(|x| (x.0 - f_worst) / (f_best - f_worst))
179      .collect();
180
181    // distribution
182    WeightedIndex::new(&weights).unwrap()
183  }
184
185  fn crossover(
186    &self,
187    parent1: &Individual,
188    parent2: &Individual,
189    rng: &mut ThreadRng,
190  ) -> Individual {
191    let mut child: Individual = (f64::INFINITY, Array1::from(vec![0.0; self.len_x]));
192    let mut i_rand_0_100: isize;
193
194    for i in 0..self.len_x {
195      i_rand_0_100 = rng.gen_range(0..=100);
196      if i_rand_0_100 % 2 == 0 {
197        child.1[i] = parent1.1[i];
198      } else {
199        child.1[i] = parent2.1[i];
200      }
201    }
202    child
203  }
204
205  fn mutate(&self, ind: &mut Individual, log10_bounds: &Vec<(f64, f64)>, rng: &mut ThreadRng) {
206    let mut f_rand_0_1: f64;
207    let mut new_x: f64;
208
209    for i in 0..self.len_x {
210      f_rand_0_1 = rng.gen_range(0.0..1.0);
211
212      if f_rand_0_1 < self.mutation_rate {
213        new_x = ind.1[i] * rng.gen_range(0.8..1.25);
214        let lb = 10f64.powf(log10_bounds[i].0);
215        let ub = 10f64.powf(log10_bounds[i].1);
216
217        if new_x < lb {
218          ind.1[i] = lb;
219        } else if new_x > ub {
220          ind.1[i] = ub;
221        } else {
222          ind.1[i] = new_x;
223        }
224      }
225    }
226  }
227
228  fn select(&self, dist: WeightedIndex<f64>, rng: &mut ThreadRng) -> (usize, usize) {
229    let index_p1 = dist.sample(rng);
230    let index_p2 = dist.sample(rng);
231    (index_p1, index_p2)
232  }
233}
aphreco/optimizer/genetic_algorithm.rs

aphreco/optimizer/
genetic_algorithm.rs
