Source code for park.diffev

from random import randint, random
import numpy

[docs]class DE: """ Differential evolution test implementation Implements the Scheme_DE_rand_1 variant F: float weighting factor which determines the magnitude of perturbation occurring in each generation. crossover_rate: float In general, 0 < crossover_rate < 1. Usually considerably lower than 1. Convergence slows as the value increases, but higher crossover rates may be necessary when func: w = f(p) The function to be minimized, of the form w = f(p), where p is a vector of length func_dim and w is a scalar func_dim: int The dimension of the objective function pop: array The initial population. This should be an iterable composed of Rank-1 numpy arrays. The population size must be at least 4, preferably much greater. l_bound, u_bound: vector arrays of size func_dim representing the upper and lower bounds for the parameters in the solution vectors tol: float if (max(pop_values) - min(pop_values) <= conv), the population has converged and the evolution will stop """ def __init__(self, F, crossover_rate, func, func_dim, pop, l_bound, u_bound, tol=1e-7): self.pop_size = len(pop) self.dimension = func_dim self.F_orig = F self.F = F self.crossover_rate = crossover_rate self.func = func self.tol = tol self.l_bound = l_bound self.u_bound = u_bound self.population = pop self.generations = 0 self.pop_values = [self.func(p) for p in self.population] self.best_gene,self.best_value = self.get_best_gene() self.best_gene_history = [self.best_gene] self.best_val_history = [self.best_value] self.ncalls = 0 #//////////////////////////////////////////////////
[docs] def evolve(self, numgens=1000, monitor=None): '''Evolve the population for numgens generations, or until it converges. Returns the best vector from the run''' try: import psyco psyco.full() except ImportError: pass start_gen = self.generations for gen in xrange(self.generations + 1, self.generations + numgens + 1): candidate_pop = [] for i in xrange(self.pop_size): (r1, r2, r3) = self.select_participants(i) #perturbation/mutation candidate = self.population[r1] + self.F*(self.population[r2] - self.population[r3]) #crossover candidate = self.crossover(candidate, i) #check bound constraints if not self.is_within_bounds(candidate): candidate = self.get_random_gene() #test fitness to determine membership in next gen candidate_val = self.func(candidate) if candidate_val < self.pop_values[i]: candidate_pop.append(candidate) self.pop_values[i] = float(candidate_val) else: candidate_pop.append(self.population[i]) self.ncalls += self.pop_size self.population = candidate_pop x,fx = self.get_best_gene() if fx < self.best_value: self.best_gene, self.best_value = x,fx if monitor is not None: monitor.improvement(x,fx,self.ncalls) self.best_val_history.append(self.best_value) self.best_gene_history.append(self.best_gene) self.generations = gen if monitor is not None: monitor.progress(gen-start_gen,numgens) if self.is_converged(): break return self.best_gene #////////////////////////////////////////////////
[docs] def get_random_gene(self): '''Generate a random gene within the bounds constraints. Used to replace out-of-bounds genes''' result = [numpy.random.uniform(low, high) for low, high in zip(self.l_bound, self.u_bound)] return numpy.array(result) #////////////////////////////////////////////////
[docs] def is_within_bounds(self, gene): '''Determine whether the gene meets the bounds constraints''' return numpy.all( (self.l_bound < gene) & (self.u_bound > gene) ) #//////////////////////////////////////////////// #This appears to be a total failure. I'll leave the code in, #but it's not called from the main evolution loop anymore.
[docs] def adjust_F(self): '''Adjust F to account for stagnation of the population ''' #has the best vector improved this generation? idx = len(self.best_val_history) if self.best_val_history[idx - 1] == self.best_val_history[idx - 2]: self.stagnant_gens += 1 else: self.stagnant_gens = 0 #adapt F to account for stagnation self.F = min( (self.F_orig + (self.stagnant_gens * 0.01)), 2.0) #////////////////////////////////////////////////
[docs] def get_best_gene(self): '''find the most fit gene in the current population''' #print "pop", self.pop_values best_index = numpy.argmin(self.pop_values) return self.population[best_index],self.pop_values[best_index], #////////////////////////////////////////////////
[docs] def select_participants(self, i): '''generate r1, r2, and r3 randomly from the range [0, NP-1] such that they are distinct values not equal to i''' choices = [i] for choice in xrange(3): while True: j = numpy.random.randint(0, self.pop_size-1) if j not in choices: break choices.append(j) return choices[1:] #////////////////////////////////////////////////
[docs] def crossover(self, candidate, i): '''Perform a crossover between the candidate and the ith member of the previous generation.''' result = [] #generate lower bound of crossover (this bound is inclusive) low = randint(0, self.dimension-1) #determine the size of the crossover L = 0 while random() < self.crossover_rate and L < (self.dimension - low): L += 1 #calculate the upper bound of crossover (this bound is exclusive) high = low + L high = numpy.repeat(high, self.dimension) low = numpy.repeat(low, self.dimension) seq = numpy.arange(0, self.dimension) result = numpy.choose( ((seq >= low)&(seq < high)), [candidate, self.population[i]] ) return result
[docs] def is_converged(self): '''check for convergence''' return max(self.pop_values) - min(self.pop_values) <= self.tol
[docs]def diffev(fn, bounds, x0=None, pop_scale=4, crossover_rate=0.8, Fscale=1, tolerance=1e-5, maxiter=1000, monitor=None): """ Run differential evolution, returning x,fx,ncalls """ lo, hi = numpy.asarray(bounds[0]), numpy.asarray(bounds[1]) dim = len(lo) pop = gen_pop(dim*pop_scale, lo, hi, dim) if x0 is not None: pop[0] = numpy.asarray(x0) evolver = DE(Fscale, crossover_rate, fn, dim, pop, lo, hi, tolerance) evolver.evolve(maxiter, monitor=monitor) return evolver.best_gene, evolver.best_value, evolver.ncalls ######################################################################## #end DE def ########################################################################
[docs]def gen_pop(size, l_bound, u_bound, dimension): '''generate a random population of vectors within the given bounds. dimension indicates the length of the vectors. l_bound and u_bound should be vectors of length dimension (any iterable container should work)''' result = [] for i in xrange(size): entry = [] for j in xrange(dimension): entry.append( ((u_bound[j] - l_bound[j]) * random()) + l_bound[j] ) result.append(numpy.array(entry)) return result