Source Code for Module park.optim.simplex

  1  #__docformat__ = "restructuredtext en" 
  2  # ******NOTICE*************** 
  3  # from optimize.py module by Travis E. Oliphant 
  4  # 
  5  # You may copy and use this module as you see fit with no 
  6  # guarantee implied provided you keep this notice in all copies. 
  7  # *****END NOTICE************ 
  8  # 
  9  # Modified by Paul Kienzle to support bounded minimization 
 10  """ 
 11  Downhill simplex optimizer. 
 12  """ 
 13   
 14  __all__ = ['simplex'] 
 15   
 16  __docformat__ = "restructuredtext en" 
 17   
 18  import numpy 
 19  __version__="0.7" 
 20   
 21 -def wrap_function(function, bounds): 
 22      ncalls = [0] 
 23      if bounds is not None: 
 24          lo, hi = [numpy.asarray(v) for v in bounds] 
 25          def function_wrapper(x): 
 26              ncalls[0] += 1 
 27              if numpy.any((x<lo)|(x>hi)): 
 28                  return numpy.inf 
 29              else: 
 30                  #function(x) 
 31                  return function(x) 
 32      else: 
 33          def function_wrapper(x): 
 34              ncalls[0] += 1 
 35              return function(x) 
 36      return ncalls, function_wrapper 
 37   
 38 -class Result: 
 39      """ 
 40      Results from the fit. 
 41   
 42      x : ndarray 
 43          Best parameter set 
 44      fx : float 
 45          Best value 
 46      iters : int 
 47          Number of iterations 
 48      calls : int 
 49          Number of function calls 
 50      success : boolean 
 51          True if the fit completed successful, false if terminated early 
 52          because of too many iterations. 
 53      """ 
 54 -    def __init__(self, x, fx, iters, calls, status): 
 55          self.x,self.fx,self.iters,self.calls=x,fx,iters,calls 
 56          self.status = status 
 57 -    def __str__(self): 
 58          return "Minimum %g at %s after %d calls"%(self.fx,self.x,self.calls) 
 59   
 60 -def dont_abort(): return False 
 61   
 62 -def simplex(f, x0=None, bounds=None, radius=0.05, 
 63              xtol=1e-4, ftol=1e-4, maxiter=None, 
 64              update_handler=None, abort_test=dont_abort): 
 65      """ 
 66      Minimize a function using Nelder-Mead downhill simplex algorithm. 
 67   
 68      This optimizer is also known as Amoeba (from Numerical Recipes) and 
 69      the Nealder-Mead simplex algorithm.  This is not the simplex algorithm 
 70      for solving constrained linear systems. 
 71   
 72      Downhill simplex is a robust derivative free algorithm for finding 
 73      minima.  It proceeds by choosing a set of points (the simplex) forming 
 74      an n-dimensional triangle, and transforming that triangle so that the 
 75      worst vertex is improved, either by stretching, shrinking or reflecting 
 76      it about the center of the triangle.  This algorithm is not known for 
 77      its speed, but for its simplicity and robustness, and is a good algorithm 
 78      to start your problem with. 
 79   
 80      *Parameters*: 
 81   
 82          f : callable f(x,*args) 
 83              The objective function to be minimized. 
 84          x0 : ndarray 
 85              Initial guess. 
 86          bounds : (ndarray,ndarray) or None 
 87              Bounds on the parameter values for the function. 
 88          radius: float 
 89              Size of the initial simplex.  For bounded parameters (those 
 90              which have finite lower and upper bounds), radius is clipped 
 91              to a value in (0,0.5] representing the portion of the 
 92              range to use as the size of the initial simplex. 
 93   
 94      *Returns*: Result (`park.simplex.Result`) 
 95   
 96          x : ndarray 
 97              Parameter that minimizes function. 
 98          fx : float 
 99              Value of function at minimum: ``fopt = func(xopt)``. 
100          iters : int 
101              Number of iterations performed. 
102          calls : int 
103              Number of function calls made. 
104          success : boolean 
105              True if fit completed successfully. 
106   
107      *Other Parameters*: 
108   
109          xtol : float 
110              Relative error in xopt acceptable for convergence. 
111          ftol : number 
112              Relative error in func(xopt) acceptable for convergence. 
113          maxiter : int=200*N 
114              Maximum number of iterations to perform.  Defaults 
115          update_handler : callable 
116              Called after each iteration, as callback(xk,fxk), where xk 
117              is the current parameter vector and fxk is the function value. 
118              Returns True if the fit should continue. 
119   
120   
121      *Notes* 
122   
123          Uses a Nelder-Mead simplex algorithm to find the minimum of 
124          function of one or more variables. 
125   
126      """ 
127      fcalls, func = wrap_function(f, bounds) 
128      x0 = numpy.asfarray(x0).flatten() 
129      #print "x0",x0 
130      N = len(x0) 
131      rank = len(x0.shape) 
132      if not -1 < rank < 2: 
133          raise ValueError, "Initial guess must be a scalar or rank-1 sequence." 
134   
135      if maxiter is None: 
136          maxiter = N * 200 
137   
138      rho = 1; chi = 2; psi = 0.5; sigma = 0.5; 
139   
140      if rank == 0: 
141          sim = numpy.zeros((N+1,), dtype=x0.dtype) 
142      else: 
143          sim = numpy.zeros((N+1,N), dtype=x0.dtype) 
144      fsim = numpy.zeros((N+1,), float) 
145      sim[0] = x0 
146      fsim[0] = func(x0) 
147   
148      # Metropolitan simplex: simplex has vertices at x0 and at 
149      # x0 + j*radius for each unit vector j.  Radius is a percentage 
150      # change from the initial value, or just the radius if the initial 
151      # value is 0.  For bounded problems, the radius is a percentage of 
152      # the bounded range in dimension j. 
153      val = x0*(1+radius) 
154      val[val == 0] = radius 
155      if bounds is not None: 
156          radius = numpy.clip(radius,0,0.5) 
157          lo,hi = [numpy.asarray(v) for v in bounds] 
158   
159          # Keep the initial simplex inside the bounds 
160          x0[x0<lo] = lo 
161          x0[x0>hi] = hi 
162          bounded = ~numpy.isinf(lo) & ~numpy.isinf(hi) 
163          val[bounded] = x0[bounded] + (hi[bounded]-lo[bounded])*radius 
164          val[val<lo] = lo 
165          val[val>hi] = hi 
166   
167          # If the initial point was at or beyond an upper bound, then bounds 
168          # projection will put x0 and x0+j*radius at the same point.  We 
169          # need to detect these collisions and reverse the radius step 
170          # direction when such collisions occur.  The only time the collision 
171          # can occur at the lower bound is when upper and lower bound are 
172          # identical.  In that case, we are already done. 
173          collision = val==x0 
174          if numpy.any(collision): 
175              reverse = x0*(1-radius) 
176              reverse[reverse==0] = -radius 
177              reverse[bounded] = x0[bounded] - (hi[bounded]-lo[bounded])*radius 
178              val[collision] = reverse[collision] 
179   
180          # Make tolerance relative for bounded parameters 
181          tol = numpy.ones(x0.shape)*xtol 
182          tol[bounded] = (hi[bounded]-lo[bounded])*xtol 
183          xtol = tol 
184   
185      # Compute values at the simplex vertices 
186      for k in range(0,N): 
187          y = x0+0 
188          y[k] = val[k] 
189          sim[k+1] = y 
190          fsim[k+1] = func(y) 
191   
192      #print sim 
193      ind = numpy.argsort(fsim) 
194      fsim = numpy.take(fsim,ind,0) 
195      # sort so sim[0,:] has the lowest function value 
196      sim = numpy.take(sim,ind,0) 
197      #print sim 
198   
199      iterations = 1 
200      while iterations < maxiter: 
201          if numpy.all(abs(sim[1:]-sim[0]) <= xtol) \ 
202              and max(abs(fsim[0]-fsim[1:])) <= ftol: 
203              #print abs(sim[1:]-sim[0]) 
204              break 
205   
206          xbar = numpy.sum(sim[:-1],0) / N 
207          xr = (1+rho)*xbar - rho*sim[-1] 
208          #print "xbar" ,xbar,rho,sim[-1],N 
209          #break 
210          fxr = func(xr) 
211          doshrink = 0 
212   
213          if fxr < fsim[0]: 
214              xe = (1+rho*chi)*xbar - rho*chi*sim[-1] 
215              fxe = func(xe) 
216   
217              if fxe < fxr: 
218                  sim[-1] = xe 
219                  fsim[-1] = fxe 
220              else: 
221                  sim[-1] = xr 
222                  fsim[-1] = fxr 
223          else: # fsim[0] <= fxr 
224              if fxr < fsim[-2]: 
225                  sim[-1] = xr 
226                  fsim[-1] = fxr 
227              else: # fxr >= fsim[-2] 
228                  # Perform contraction 
229                  if fxr < fsim[-1]: 
230                      xc = (1+psi*rho)*xbar - psi*rho*sim[-1] 
231                      fxc = func(xc) 
232   
233                      if fxc <= fxr: 
234                          sim[-1] = xc 
235                          fsim[-1] = fxc 
236                      else: 
237                          doshrink=1 
238                  else: 
239                      # Perform an inside contraction 
240                      xcc = (1-psi)*xbar + psi*sim[-1] 
241                      fxcc = func(xcc) 
242   
243                      if fxcc < fsim[-1]: 
244                          sim[-1] = xcc 
245                          fsim[-1] = fxcc 
246                      else: 
247                          doshrink = 1 
248   
249                  if doshrink: 
250                      for j in xrange(1,N+1): 
251                          sim[j] = sim[0] + sigma*(sim[j] - sim[0]) 
252                          fsim[j] = func(sim[j]) 
253   
254          ind = numpy.argsort(fsim) 
255          sim = numpy.take(sim,ind,0) 
256          fsim = numpy.take(fsim,ind,0) 
257          if update_handler is not None: 
258              update_handler(sim[0],fsim[0]) 
259          iterations += 1 
260          if abort_test(): break 
261   
262      status = 0 if iterations < maxiter else 1 
263      res = Result(sim[0],fsim[0],iterations,fcalls[0], status) 
264      return res 
265   
266 -def main(): 
267      import time 
268      def rosen(x):  # The Rosenbrock function 
269          return numpy.sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0,axis=0) 
270   
271   
272      x0 = [0.8,1.2,0.7] 
273      print "Nelder-Mead Simplex" 
274      print "===================" 
275      start = time.time() 
276      x = simplex(rosen,x0) 
277      print x 
278      print "Time:",time.time() - start 
279   
280      x0 = [0]*3 
281      print "Nelder-Mead Simplex" 
282      print "===================" 
283      print "starting at zero" 
284      start = time.time() 
285      x = simplex(rosen,x0) 
286      print x 
287      print "Time:",time.time() - start 
288   
289      x0 = [0.8,1.2,0.7] 
290      lo,hi = [0]*3, [1]*3 
291      print "Bounded Nelder-Mead Simplex" 
292      print "===========================" 
293      start = time.time() 
294      x = simplex(rosen,x0,bounds=(lo,hi)) 
295      print x 
296      print "Time:",time.time() - start 
297   
298   
299      x0 = [0.8,1.2,0.7] 
300      lo,hi = [0.999]*3, [1.001]*3 
301      print "Bounded Nelder-Mead Simplex" 
302      print "===========================" 
303      print "tight bounds" 
304      print "simplex is smaller than 1e-7 in every dimension, but you can't" 
305      print "see this without uncommenting the print statement simplex function" 
306      start = time.time() 
307      x = simplex(rosen,x0,bounds=(lo,hi),xtol=1e-4) 
308      print x 
309      print "Time:",time.time() - start 
310   
311   
312      x0 = [0]*3 
313      hi,lo = [-0.999]*3, [-1.001]*3 
314      print "Bounded Nelder-Mead Simplex" 
315      print "===========================" 
316      print "tight bounds, x0=0 outside bounds from above" 
317      start = time.time() 
318      x = simplex(lambda x:rosen(-x),x0,bounds=(lo,hi),xtol=1e-4) 
319      print x 
320      print "Time:",time.time() - start 
321   
322   
323      x0 = [0.8,1.2,0.7] 
324      lo,hi = [-numpy.inf]*3, [numpy.inf]*3 
325      print "Bounded Nelder-Mead Simplex" 
326      print "===========================" 
327      print "infinite bounds" 
328      start = time.time() 
329      x = simplex(rosen,x0,bounds=(lo,hi),xtol=1e-4) 
330      print x 
331      print "Time:",time.time() - start 
332   
333  if __name__ == "__main__": 
334      main() 
335