1
# Author: Travis Oliphant 2002
3
from __future__ import nested_scopes
6
from scipy_base import asarray, tan, exp, ones, squeeze, sign, \
26
def init(self, **options):
27
self.__dict__.update(options)
28
if self.lower == scipy_base.NINF:
29
self.lower = -scipy_base.limits.double_max
30
if self.upper == scipy_base.PINF:
31
self.upper = scipy_base.limits.double_max
37
def getstart_temp(self, best_state):
38
assert(not self.dims is None)
39
x0 = ones(self.dims,'d')
40
lrange = x0*self.lower
41
urange = x0*self.upper
44
for n in range(self.Ninit):
45
for k in range(self.dims):
46
x0[k] = random.uniform(lrange[k],urange[k])
47
fval = self.func(x0,*self.args)
53
best_state.cost = fval
54
best_state.x = asarray(x0).copy()
55
self.T0 = (fmax-fmin)*1.5
58
def accept_test(self, dE):
64
p = exp(-dE*1.0/self.boltzmann/T)
65
if (p > random.uniform(0.0,1.0)):
70
def update_guess(self, x0):
73
def update_temp(self, x0):
77
# A schedule due to Lester Ingber
78
class fast_sa(base_schedule):
79
def init(self, **options):
80
self.__dict__.update(options)
85
self.c = self.m * exp(-self.n * self.quench / self.dims)
87
def update_guess(self, x0):
89
u = squeeze(asarray([random.uniform(0.0,1.0) for x in x0]))
91
y = sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)
92
xc = y*(self.upper - self.lower)
96
def update_temp(self):
97
self.T = self.T0*exp(-self.c * self.k**(self.quench/self.dims))
101
class cauchy_sa(base_schedule):
102
def update_guess(self, x0):
104
numbers = squeeze(asarray([random.uniform(-pi/2,pi/2) for x in x0]))
105
xc = self.learn_rate * self.T * tan(numbers)
109
def update_temp(self):
110
self.T = self.T0/(1+self.k)
114
class boltzmann_sa(base_schedule):
115
def update_guess(self, x0):
116
std = min(sqrt(self.T), (self.upper-self.lower)/3.0/self.learn_rate)
118
xc = squeeze(asarray([random.normalvariate(0,std*self.learn_rate) \
123
def update_temp(self):
125
self.T = self.T0 / log(self.k+1.0)
134
# allow for general annealing temperature profile
135
# in that case use update given by alpha and omega and
136
# variation of all previous updates and temperature?
138
# Simulated annealing
140
def anneal(func, x0, args=(), schedule='fast', full_output=0,
141
T0=None, Tf=1e-12, maxeval=None, maxaccept=None, maxiter=400,
142
boltzmann=1.0, learn_rate=0.5, feps=1e-6, quench=1.0, m=1.0, n=1.0,
143
lower=-100, upper=100, dwell=50):
144
"""Minimize a function using simulated annealing.
146
Schedule is a schedule class implementing the annealing schedule.
147
Available ones are 'fast', 'cauchy', 'boltzmann'
151
func -- Function to be optimized
152
x0 -- Parameters to be optimized over
153
args -- Extra parameters to function
154
schedule -- Annealing schedule to use (a class)
155
full_output -- Return optional outputs
156
T0 -- Initial Temperature (estimated as 1.2 times the largest
157
cost-function deviation over random points in the range)
158
Tf -- Final goal temperature
159
maxeval -- Maximum function evaluations
160
maxaccept -- Maximum changes to accept
161
maxiter -- Maximum cooling iterations
162
learn_rate -- scale constant for adjusting guesses
163
boltzmann -- Boltzmann constant in acceptance test
164
(increase for less stringent test at each temperature).
165
feps -- Stopping relative error tolerance for the function value in
167
quench, m, n -- Parameters to alter fast_sa schedule
168
lower, upper -- lower and upper bounds on x0 (scalar or array).
169
dwell -- The number of times to search the space at each temperature.
171
Outputs: (xmin, {Jmin, T, feval, iter, accept,} retval)
173
xmin -- Point giving smallest value found
174
retval -- Flag indicating stopping condition:
175
0 : Cooled to global optimum
176
1 : Cooled to final temperature
177
2 : Maximum function evaluations
178
3 : Maximum cooling iterations reached
179
4 : Maximum accepted query locations reached
181
Jmin -- Minimum value of function found
182
T -- final temperature
183
feval -- Number of function evaluations
184
iter -- Number of cooling iterations
185
accept -- Number of tests accepted.
189
schedule = eval(schedule+'_sa()')
190
# initialize the schedule
191
schedule.init(dims=len(x0),func=func,args=args,boltzmann=boltzmann,T0=T0,
192
learn_rate=learn_rate, lower=lower, upper=upper,
193
m=m, n=n, quench=quench, dwell=dwell)
195
current_state, last_state, best_state = _state(), _state(), _state()
199
x0 = schedule.getstart_temp(best_state)
202
best_state.cost = 300e8
203
last_state.x = asarray(x0).copy()
204
fval = func(x0,*args)
206
last_state.cost = fval
207
if last_state.cost < best_state.cost:
208
best_state.cost = fval
209
best_state.x = asarray(x0).copy()
210
schedule.T = schedule.T0
211
fqueue = [100,300,500,700]
214
for n in range(dwell):
215
xnew = schedule.update_guess(x0)
216
fval = func(xnew,*args)
218
current_state.x = asarray(xnew).copy()
219
current_state.cost = fval
220
dE = current_state.cost - last_state.cost
221
if schedule.accept_test(dE):
223
last_state.x = current_state.x.copy()
224
last_state.cost = current_state.cost
225
if last_state.cost < best_state.cost:
226
best_state.x = last_state.x.copy()
227
best_state.cost = last_state.cost
228
schedule.update_temp()
230
# Stopping conditions
231
# 0) last saved values of f from each cooling step
232
# are all very similar (effectively cooled)
233
# 1) Tf is set and we are below it
234
# 2) maxeval is set and we are past it
235
# 3) maxiter is set and we are past it
236
# 4) maxaccept is set and we are past it
238
fqueue.append(squeeze(last_state.cost))
240
af = asarray(fqueue)*1.0
241
if all(abs((af-af[0])/af[0]) < feps):
243
if abs(af[-1]-best_state.cost) > feps*10:
245
print "Warning: Cooled to %f at %f but this is not" \
246
% (squeeze(last_state.cost), squeeze(last_state.x)) \
247
+ " the smallest point found."
249
if (Tf is not None) and (schedule.T < Tf):
252
if (maxeval is not None) and (schedule.feval > maxeval):
256
print "Warning: Maximum number of iterations exceeded."
259
if (maxaccept is not None) and (schedule.accepted > maxaccept):
264
return best_state.x, best_state.cost, schedule.T, \
265
schedule.feval, iter, schedule.accepted, retval
267
return best_state.x, retval
271
if __name__ == "__main__":
272
func = lambda x: cos(14.5*x-0.3) + (x+0.2)*x
273
#print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=1000,schedule=fast_sa())
274
#print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=1000,schedule=cauchy_sa())
275
#print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=1000,schedule=boltzmann_sa())
added
removed
Lib/optimize/anneal.py
