Source Code for Module park.optim.fitmc

  1   
  2  # The job queue needs to be in a transaction/rollback protected database. 
  3  # If the server is rebooted, long running jobs should continue to work. 
  4  # 
  5  from __future__ import division 
  6  import numpy 
  7   
  8  from park.optim import pmap, fitutil 
  9  from park.fitting import fitresult 
 10   
 11  __all__ = ['fitmc', 'FitMCJob'] 
 12   
 13 -class LocalFit(object): 
 14      """ 
 15      Abstract interface for a local minimizer 
 16   
 17      See `park.fitmc.FitSimplex` for a concrete implementation. 
 18      """ 
 19 -    def fit(self, objective, x0): 
 20          """ 
 21          Minimize the value of a fitness function. 
 22   
 23          See `park.fitmc.Fitness` for the definition of the goodness of fit 
 24          object.  x0 is a vector containing the initial value for the fit. 
 25          """ 
 26 -    def abort(self): 
 27          """ 
 28          Cancel the fit.  This will be called from a separate execution 
 29          thread.  The fit should terminate as soon as possible after this 
 30          function has been called.  Ideally this would interrupt the 
 31          cur 
 32          """ 
 33   
 34 -class FitSimplex(LocalFit): 
 35      """ 
 36      Local minimizer using Nelder-Mead simplex algorithm. 
 37   
 38      Simplex is robust and derivative free, though not very efficient. 
 39   
 40      This class wraps the bounds contrained Nelder-Mead simplex 
 41      implementation for `park.simplex.simplex`. 
 42      """ 
 43      radius = 0.05 
 44      """Size of the initial simplex; this is a portion between 0 and 1""" 
 45      xtol = 1e-4 
 46      """Stop when simplex vertices are within xtol of each other""" 
 47      ftol = 1e-4 
 48      """Stop when vertex values are within ftol of each other""" 
 49      maxiter = None 
 50      """Maximum number of iterations before fit terminates""" 
 51 -    def fit(self, fitness, x0): 
 52          """Run the fit""" 
 53          import simplex 
 54   
 55          self.cancel = False 
 56          pars = fitness.fit_parameters() 
 57          bounds = numpy.array([p.range for p in pars]).T 
 58          result = simplex.simplex(fitness, x0, bounds=bounds, 
 59                                   radius=self.radius, xtol=self.xtol, 
 60                                   ftol=self.ftol, maxiter=self.maxiter, 
 61                                   abort_test=self._iscancelled) 
 62          #print "simplex returned",result.x,result.fx 
 63          # Need to make our own copy of the fit results so that the 
 64          # values don't get stomped on by the next fit iteration. 
 65          fitpars = [fitresult.FitParameter(pars[i].name,pars[i].range,v) 
 66                     for i,v in enumerate(result.x)] 
 67          res = fitresult.FitResult(fitpars, result.calls, result.fx) 
 68          # Compute the parameter uncertainties from the jacobian 
 69          res.calc_cov(fitness) 
 70          return res 
 71   
 72 -    def abort(self): 
 73          """Cancel the fit in progress from another thread of execution""" 
 74          # Effectively an atomic operation; rely on GIL to protect it. 
 75          self.cancel = True 
 76          # Abort the current function evaluation if possible. 
 77          if hasattr(fitness,'abort'): self.fitness.abort() 
 78   
 79 -    def _iscancelled(self): return self.cancel 
 80   
 81 -class FitCobyla(LocalFit): 
 82      """ 
 83      Local minimizer using Nelder-Mead simplex algorithm. 
 84   
 85      Simplex is robust and derivative free, though not very efficient. 
 86   
 87      This class wraps the bounds contrained Nelder-Mead simplex 
 88      implementation for `park.simplex.simplex`. 
 89      """ 
 90      radius = 0.05 
 91      """Size of the initial simplex; this is a portion between 0 and 1""" 
 92      xtol = 1e-4 
 93      """Stop when simplex vertices are within xtol of each other""" 
 94      ftol = 1e-4 
 95      """Stop when vertex values are within ftol of each other""" 
 96      maxiter = None 
 97      """Maximum number of iterations before fit terminates""" 
 98 -    def fit(self, fitness, x0): 
 99          """Run the fit""" 
100          from scipy.optimize import fmin_cobyla as cobyla 
101   
102          self.cancel = False 
103          pars = fitness.fit_parameters() 
104          bounds = numpy.array([p.range for p in pars]).T 
105          result = simplex.simplex(fitness, x0, bounds=bounds, 
106                                   radius=self.radius, xtol=self.xtol, 
107                                   ftol=self.ftol, maxiter=self.maxiter, 
108                                   abort_test=self._iscancelled) 
109          #print "simplex returned",result.x,result.fx 
110          # Need to make our own copy of the fit results so that the 
111          # values don't get stomped on by the next fit iteration. 
112          fitpars = [fitresult.FitParameter(pars[i].name,pars[i].range,v) 
113                     for i,v in enumerate(result.x)] 
114          res = fitresult.FitResult(fitpars, result.calls, result.fx) 
115          # Compute the parameter uncertainties from the jacobian 
116          res.calc_cov(fitness) 
117          return res 
118   
119 -    def abort(self): 
120          """Cancel the fit in progress from another thread of execution""" 
121          # Effectively an atomic operation; rely on GIL to protect it. 
122          self.cancel = True 
123          # Abort the current function evaluation if possible. 
124          if hasattr(fitness,'abort'): self.fitness.abort() 
125   
126 -    def _iscancelled(self): return self.cancel 
127   
128 -class MapMC(object): 
129      """ 
130      Evaluate a local fit at a particular start point. 
131   
132      This is the function to be mapped across a set of start points for the 
133      monte carlo map-reduce implementation. 
134   
135      See `park.pmap.Mapper` for required interface. 
136      """ 
137 -    def __init__(self, minimizer, fitness): 
138          self.minimizer, self.fitness = minimizer, fitness 
139 -    def __call__(self, x0): 
140          return self.minimizer.fit(self.fitness,x0) 
141 -    def abort(self): 
142          self.minimizer.abort() 
143   
144 -class CollectMC(object): 
145      """ 
146      Collect the results from multiple start points in a Monte Carlo fit engine. 
147   
148      See `park.pmap.Collector` for required interface. 
149      """ 
150 -    def __init__(self, number_expected, handler): 
151          self.number_expected = number_expected 
152          """Number of starting points to check with local optimizer""" 
153          self.iterations = 0 
154          self.best = None 
155          self.calls = 0 
156          self.handler = handler 
157          self.handler.done = False 
158 -    def __call__(self, result): 
159          # Keep track of the amount of work done 
160          self.iterations += 1 
161          self.calls += result.calls 
162          if self.best is None or result.fitness < self.best.fitness: 
163              self.best = result 
164              self.handler.result = result 
165              self.handler.improvement() 
166          # Progress should go after improvement in case the fit handler 
167          # wants to suppress intermediate improvements 
168          self.handler.progress(self.iterations, self.number_expected) 
169 -    def abort(self): 
170          #print "aborting" 
171          self.handler.done = True 
172          self.handler.abort() 
173 -    def finalize(self): 
174          #print "finalizing" 
175          self.handler.done = True 
176          self.handler.finalize() 
177 -    def error(self, msg): 
178          #print "error",msg 
179          self.handler.done = True 
180          self.handler.error(msg) 
181   
182 -def fitmc(fitness,localfit,n,handler): 
183      """ 
184      Run a monte carlo fit. 
185   
186      This procedure maps a local optimizer across a set of n initial points. 
187      The initial parameter value defined by the fitness parameters defines 
188      one initial point.  The remainder are randomly generated within the 
189      bounds of the problem. 
190   
191      localfit is the local optimizer to use.  It should be a bounded 
192      optimizer following the `park.fitmc.LocalFit` interface. 
193   
194      handler accepts updates to the current best set of fit parameters. 
195      See `park.fitresult.FitHandler` for details. 
196      """ 
197      P = fitutil.populate(n, fitness.fit_parameters(), radius=100, elitism=True) 
198      pmap.pmapreduce(MapMC(localfit,fitness), 
199                      CollectMC(n,handler), 
200                      P) 
201   
202 -class FitMC(object): 
203      """ 
204      Monte Carlo optimizer. 
205   
206      This implements `park.fit.Fitter`. 
207      """ 
208   
209 -    def __init__(self, 
210                   localfit=FitSimplex(), 
211                   start_points=10): 
212          """ 
213          Parameters are defined by `park.fitmc.fitmc`. 
214          """ 
215          self.localfit = localfit 
216          self.start_points = start_points 
217   
218 -    def fit(self, objective, handler): 
219          """ 
220          Run a monte carlo fit. 
221   
222          This procedure maps a local optimizer across a set of initial points. 
223          """ 
224          fitmc(objective, self.localfit, self.start_points, handler) 
225   
226 -def main(): 
227      """Multiple minima example""" 
228      import time, math 
229      class MultiMin(object): 
230          def fit_parameters(self): 
231              return [fitresult.FitParameter('x1',[-5,5],1)] 
232          def __call__(self, p): 
233              x=p[0] 
234              fx = x**2 + math.sin(2*math.pi*x+3*math.pi/2) 
235              return fx 
236      handler = fitresult.ConsoleUpdate() # Show updates on the console 
237      handler.progress_delta = 1          # Update user every second 
238      handler.improvement_delta = 0.1     # Show improvements almost immediately 
239      fitmc(MultiMin(), FitSimplex(), 5000, handler) 
240      while not handler.done: time.sleep(1) 
241   
242 -def assembly_example(): 
243      import park, time 
244      problem = park.assembly.example() 
245      handler = fitresult.ConsoleUpdate() # Show updates on the console 
246      handler.progress_delta = 1          # Update user every second 
247      handler.improvement_delta = 1       # Show improvements at the same rate 
248      #pmap.profile_mapper = True 
249      fitmc(problem, FitSimplex(), 500, handler) 
250      while not handler.done: time.sleep(1) 
251   
252  if __name__ == "__main__": 
253      #main() 
254      assembly_example() 
255