← Back to branch summary

~ubuntu-branches/ubuntu/saucy/python-scipy/saucy

« back to all changes in this revision

Viewing changes to Lib/sandbox/rkern/diffev.py

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
http://bugs.debian.org/489149
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

Show diffs side-by-side

added added

removed removed

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