Source Code for Module park.optim.diffev

  1  """Differential Evolution Optimization 
  2   
  3  :Author: Robert Kern 
  4   
  5  Copyright 2005 by Robert Kern. 
  6  """ 
  7   
  8  import numpy as np 
  9   
 10   
 11  # Notes: for future modifications: 
 12  # Ali, M. M., and A. Toern. Topographical differential evolution using 
 13  # pre-calculated differentials. _Stochastic and Global Optimization_. 1--17. 
 14  # 
 15  #  A good scale value: 
 16  #    F = max(l_min, 1-min(abs(f_min/f_max), abs(f_max/f_min))) 
 17  #      ~ 0.3 <= l_min <= 0.4 
 18  #      ~ f_min and f_max are the minimum and maximum values in the initial 
 19  #        population. 
 20  # 
 21  #  Pre-calculated differentials: 
 22  #    Keep a set of differentials A. 
 23  #    For each x_i of the population S: 
 24  #      Every M steps, randomly select 3 points x_r1, x_r2, x_r3 from S (not x_i). 
 25  #        Compute new x_i using x_r1 + F*(x_r2-x_r3). 
 26  #        Store differential in A. 
 27  #      Each other step: 
 28  #        Randomly select x_r1 from S and a differential vector from A. 
 29  #      Crossover. 
 30  # 
 31  #  Convergence criterion: 
 32  #    f_max - f_min < eps 
 33  # 
 34  #  Topographical DEPD: 
 35  #    Two populations S and Sa (auxiliary). 
 36  #    Phase counter t = 0 and array shift[:] = False. 
 37  #    Stopping condition: e.g. t >= 4. 
 38  #    g << N, number of nearest neighbors to search for graph minima. 
 39  #    Ng << N, number of points for graph. 
 40  #    For each x_i in S, do DEPD as described above to get y_i. 
 41  #    if f(y_i) < f(x_i): 
 42  #      if shift[i] is False: 
 43  #        shift[i] = True 
 44  #        S[i] = y_i 
 45  #      else: 
 46  #        Sa[i] = y_i 
 47  #      if alltrue(shift,axis=0): 
 48  #        Find graph minima of f(x) using the Ng best points in S. 
 49  #        Do local search from each minimum. 
 50  #        Replace worst Ng points in S with best Ng points in Sa. 
 51  #        If best from this phase is better than previous best, t=0. 
 52  #        Else: t += 1. 
 53  #        shift[:] = False 
 54  #    Next generation. 
 55   
 56 -class DiffEvolver(object): 
 57      """Minimize a function using differential evolution. 
 58   
 59      Constructors 
 60      ------------ 
 61      DiffEvolver(func, pop0, args=(), crossover_rate=0.5, scale=None, 
 62          strategy=('rand', 2, 'bin'), eps=1e-6) 
 63        func -- function to minimize 
 64        pop0 -- sequence of initial vectors 
 65        args -- additional arguments to apply to func 
 66        crossover_rate -- crossover probability [0..1] usually 0.5 or so 
 67        scale -- scaling factor to apply to differences [0..1] usually > 0.5 
 68          if None, then calculated from pop0 using a heuristic 
 69        strategy -- tuple specifying the differencing/crossover strategy 
 70          The first element is one of 'rand', 'best', 'rand-to-best' to specify 
 71          how to obtain an initial trial vector. 
 72          The second element is either 1 or 2 (or only 1 for 'rand-to-best') to 
 73          specify the number of difference vectors to add to the initial trial. 
 74          The third element is (currently) 'bin' to specify binomial crossover. 
 75        eps -- if the maximum and minimum function values of a given generation are 
 76          with eps of each other, convergence has been achieved. 
 77        prng -- a RandomState instance. By default, this is the global 
 78          numpy.random instance. 
 79   
 80      DiffEvolver.frombounds(func, lbound, ubound, npop, crossover_rate=0.5, 
 81          scale=None, strategy=('rand', 2, 'bin'), eps=1e-6) 
 82        Randomly initialize the population within given rectangular bounds. 
 83        lbound -- lower bound vector 
 84        ubound -- upper bound vector 
 85        npop -- size of population 
 86   
 87      Public Methods 
 88      -------------- 
 89      solve(newgens=100) 
 90        Run the minimizer for newgens more generations. Return the best parameter 
 91        vector from the whole run. 
 92   
 93      Public Members 
 94      -------------- 
 95      best_value -- lowest function value in the history 
 96      best_vector -- minimizing vector 
 97      best_val_history -- list of best_value's for each generation 
 98      best_vec_history -- list of best_vector's for each generation 
 99      population -- current population 
100      pop_values -- respective function values for each of the current population 
101      generations -- number of generations already computed 
102      func, args, crossover_rate, scale, strategy, eps -- from constructor 
103      """ 
104 -    def __init__(self, func, pop0, args=(), crossover_rate=0.5, scale=None, 
105              strategy=('rand', 2, 'bin'), eps=1e-6, prng=np.random): 
106          self.func = func 
107          self.population = np.array(pop0) 
108          self.npop, self.ndim = self.population.shape 
109          self.args = args 
110          self.crossover_rate = crossover_rate 
111          self.strategy = strategy 
112          self.eps = eps 
113   
114          self.pop_values = [self.func(m, *args) for m in self.population] 
115          bestidx = np.argmin(self.pop_values) 
116          self.best_vector = self.population[bestidx] 
117          self.best_value = self.pop_values[bestidx] 
118   
119          if scale is None: 
120              self.scale = self.calculate_scale() 
121          else: 
122              self.scale = scale 
123   
124          self.generations = 0 
125          self.best_val_history = [] 
126          self.best_vec_history = [] 
127   
128          self.jump_table = { 
129              ('rand', 1, 'bin'): (self.choose_rand, self.diff1, self.bin_crossover), 
130              ('rand', 2, 'bin'): (self.choose_rand, self.diff2, self.bin_crossover), 
131              ('best', 1, 'bin'): (self.choose_best, self.diff1, self.bin_crossover), 
132              ('best', 2, 'bin'): (self.choose_best, self.diff2, self.bin_crossover), 
133              ('rand-to-best', 1, 'bin'): 
134                  (self.choose_rand_to_best, self.diff1, self.bin_crossover), 
135              } 
136   
137 -    def clear(self): 
138          self.best_val_history = [] 
139          self.best_vec_history = [] 
140          self.generations = 0 
141          self.pop_values = [self.func(m, *self.args) for m in self.population] 
142   
143 -    def frombounds(cls, func, lbound, ubound, npop, crossover_rate=0.5, 
144              scale=None, strategy=('rand', 2, 'bin'), eps=1e-6, prng=np.random): 
145          lbound = np.asarray(lbound) 
146          ubound = np.asarray(ubound) 
147          pop0 = prng.uniform(lbound, ubound, size=(npop, len(lbound))) 
148          return cls(func, pop0, crossover_rate=crossover_rate, scale=scale, 
149              strategy=strategy, eps=eps) 
150      frombounds = classmethod(frombounds) 
151   
152 -    def calculate_scale(self): 
153          rat = abs(max(self.pop_values)/self.best_value) 
154          rat = min(rat, 1./rat) 
155          return max(0.3, 1.-rat) 
156   
157 -    def bin_crossover(self, oldgene, newgene): 
158          mask = self.prng.rand(self.ndim) < self.crossover_rate 
159          return np.where(mask, newgene, oldgene) 
160   
161 -    def select_samples(self, candidate, nsamples): 
162          possibilities = range(self.npop) 
163          possibilities.remove(candidate) 
164          return self.prng.permutation(possibilities)[:nsamples] 
165   
166 -    def diff1(self, candidate): 
167          i1, i2 = self.select_samples(candidate, 2) 
168          return self.scale * (self.population[i1] - self.population[i2]) 
169   
170 -    def diff2(self, candidate): 
171          i1, i2, i3, i4 = self.select_samples(candidate, 4) 
172          return self.scale * (self.population[i1] - self.population[i2] + 
173                               self.population[i3] - self.population[i4]) 
174   
175 -    def choose_best(self, candidate): 
176          return self.best_vector 
177   
178 -    def choose_rand(self, candidate): 
179          i = self.select_samples(candidate, 1)[0] 
180          return self.population[i] 
181   
182 -    def choose_rand_to_best(self, candidate): 
183          return ((1-self.scale) * self.population[candidate] + 
184                  self.scale * self.best_vector) 
185   
186 -    def get_trial(self, candidate): 
187          chooser, differ, crosser = self.jump_table[self.strategy] 
188          trial = crosser(self.population[candidate], 
189              chooser(candidate) + differ(candidate)) 
190          return trial 
191   
192 -    def converged(self): 
193          return max(self.pop_values) - min(self.pop_values) <= self.eps 
194   
195 -    def solve(self, newgens=100): 
196          """ 
197          Run for newgens more generations. 
198   
199          Return best parameter vector from the entire run. 
200          """ 
201          for gen in xrange(self.generations+1, self.generations+newgens+1): 
202              # Set this up as a mapper for simple parallelism 
203              # Better strategies exist. 
204              trials = [self.get_trial(c) for c in range(self.npop)] 
205              values = map(self.func,trials) 
206              for candidate, trial_value in enumerate(values): 
207                  if value < self.pop_values[candidate]: 
208                      self.population[candidate] = trials[candidate] 
209                      self.pop_values[candidate] = trial_value 
210                      if trial_value < self.best_value: 
211                          self.best_vector = trial 
212                          self.best_value = trial_value 
213              self.best_val_history.append(self.best_value) 
214              self.best_vec_history.append(self.best_vector) 
215              if self.converged(): 
216                  break 
217          self.generations = gen 
218          return self.best_vector 
219