park.fitting.fitresult

1 import sys, math, time 2 import numpy 3 4 from park.util.formatnum import format_uncertainty as format 5

6 -class FitHandler(object):

7 """ 8 Abstract interface for fit thread handler. 9 10 The methods in this class are called by the optimizer as the fit 11 progresses. 12 13 Note that it is up to the optimizer to call the fit handler correctly, 14 reporting all status changes and maintaining the 'done' flag. 15 """ 16 done = False 17 """True when the fit job is complete""" 18 result = None 19 """The current best result of the fit""" 20

21 - def improvement(self):

22 """ 23 Called when a result is observed which is better than previous 24 results from the fit. 25 26 result is a FitResult object, with parameters, #calls and fitness. 27 """

28 - def error(self, msg):

29 """ 30 Model had an error; print traceback 31 """ 32 raise Exception(msg)

33 - def progress(self, current, expected):

34 """ 35 Called each cycle of the fit, reporting the current and the 36 expected amount of work. The meaning of these values is 37 optimizer dependent, but they can be converted into a percent 38 complete using (100*current)//expected. 39 40 Progress is updated each iteration of the fit, whatever that 41 means for the particular optimization algorithm. It is called 42 after any calls to improvement for the iteration so that the 43 update handler can control I/O bandwidth by suppressing 44 intermediate improvements until the fit is complete. 45 """

46 - def finalize(self):

47 """ 48 Fit is complete; best results are reported 49 """

50 - def abort(self):

51 """ 52 Fit was aborted. 53 """

54 55 # TODO: Checkpoint/restore. 56 # This should be driven by the computation, which knows when it is in a 57 # good position to save and restore state, rather than the executive, 58 # which at best could save and restore the state of a virtual machine 59 # or maybe ask the process nicely from a separate thread if it could 60 # please save itself.

61 -class FitController(object):

62 """ 63 Operations required to support functionality of GA refl. 64 """

65 - def __init__(self, pars):

66 self.parameterset = pars

67

68 - def set_bounds(self, p, range):

69 """ 70 Reset the bounds on parameter p. 71 72 Range can be:: 73 None: let p vary in (-inf,inf) 74 value: fix p at the value 75 (low,high): vary p in range (low,high) 76 """

77

78 - def scramble(self, p):

79 """ 80 Reset the fit population so that parameter p varies uniformly 81 within its range. 82 83 This is not a sensible option for all optimizers. 84 """

85

86 - def localopt(self, p=None, optimizer=None):

87 """ 88 Run a local optimizer at point p. 89 90 The result will be sent to the FitMonitor. The user may opt to 91 introduce this new individual back into the fit population. 92 93 If no p is given, then run from the best current value in the fit. 94 95 If no optimizer is given, then run a quasi-Newton optimizer 96 such as secant or dog-leg for functions without derivatives 97 or tnc for functions with derivatives. 98 99 Note: it is theoretically possible to split the minimization 100 into D0 and D1 portions, using different algorithms for each. 101 """

102

103 - def suggest(self, p):

104 """ 105 Introduce an individual to the fit population. 106 107 The suggested value may come from a local optimization from point 108 p or from the user playing with the parameters in the foreground 109 while the global optimization is running in the background. 110 """

111

112 - def lamarkian(self, p=None, optimizer=None):

113 """ 114 Lamarkian evolution. 115 116 Improve the fittest individual in the population using a local 117 optimizer and reintroduce it into the population. This strategy 118 is meaningful for a number of the stochastic optimizers including 119 simulated annealing since it allows the user to control the goal 120 directed behaviour. 121 122 This is a server-side combination of localopt and suggest, without 123 the need to cycle the results through FitMonitor. 124 """

125

126 - def attract(self, p, sigma=0.1):

127 """ 128 Search near p. 129 130 Sigma is the size of the attractor relative to the entire search space. 131 132 This is a user supplied heuristic telling the stochastic search 133 to prefer to look near point p when proposing new members of the 134 population. This will not affect the calculated objective value 135 136 The attractor translates into a scaled gaussian peak. 137 """

138

139 - def repel(self, p, I=1, sigma=0.1):

140 """ 141 Avoid p. 142 Avoid a particular region of the fit space when choosing new 143 individuals in the population. 144 145 Having identified a bad local optimum in the fit space, have 146 the stochastic optimizers avoid the this region using a penalty 147 function based 148 """

149 150

151 -class ConsoleUpdate(FitHandler):

152 """ 153 Print progress to the console. 154 """ 155 isbetter = False 156 """Record whether results improved since last update""" 157 progress_delta = 60 158 """Number of seconds between progress updates""" 159 improvement_delta = 5 160 """Number of seconds between improvement updates"""

161 - def __init__(self,quiet=False,progress_delta=60,improvement_delta=5):

162 """ 163 If quiet is true, only print out final summary, not progress and 164 improvements. 165 """ 166 self.progress_time = time.time() 167 self.progress_percent = 0 168 self.improvement_time = self.progress_time 169 self.isbetter = False 170 self.quiet = quiet 171 self.progress_delta = progress_delta 172 self.improvement_delta = improvement_delta

173

174 - def progress(self, k, n):

175 """ 176 Report on progress. 177 """ 178 if self.quiet: return 179 t = time.time() 180 p = int((100*k)//n) 181 182 # Show improvements if there are any 183 dt = t - self.improvement_time 184 if self.isbetter and dt > self.improvement_delta: 185 self.result.print_summary() 186 self.isbetter = False 187 self.improvement_time = t 188 189 # Update percent complete 190 dp = p-self.progress_percent 191 if dp < 1: return 192 dt = t - self.progress_time 193 if dt > self.progress_delta: 194 if 1 <= dp <= 2: 195 print "%d%% complete"%p 196 self.progress_percent = p 197 self.progress_time = t 198 elif 2 < dp <= 5: 199 if p//5 != self.progress_percent//5: 200 print "%d%% complete"%(5*(p//5)) 201 self.progress_percent = p 202 self.progress_time = t 203 else: 204 if p//10 != self.progress_percent//10: 205 print "%d%% complete"%(10*(p//10)) 206 self.progress_percent = p 207 self.progress_time = t

208

209 - def improvement(self):

210 """ 211 Called when a result is observed which is better than previous 212 results from the fit. 213 """ 214 self.isbetter = True

215

216 - def error(self, msg):

217 """ 218 Model had an error; print traceback 219 """ 220 if self.isbetter: 221 self.result.print_summary() 222 print msg

223

224 - def finalize(self):

225 if self.isbetter: 226 self.result.print_summary()

227

228 - def abort(self):

229 if self.isbetter: 230 self.result.print_summary()

231 232

233 -class FitParameter(object):

234 """ 235 Fit result for an individual parameter. 236 """

237 - def __init__(self, name, range, value):

238 self.name = name 239 self.range = range 240 self.value = value 241 self.stderr = None

242 - def summarize(self):

243 """ 244 Return parameter range string. 245 246 E.g., " Gold .....|.... 5.2043 in [2,7]" 247 """ 248 bar = ['.']*10 249 lo,hi = self.range 250 if numpy.isfinite(lo)and numpy.isfinite(hi): 251 portion = (self.value-lo)/(hi-lo) 252 if portion < 0: portion = 0. 253 elif portion >= 1: portion = 0.99999999 254 barpos = int(math.floor(portion*len(bar))) 255 bar[barpos] = '|' 256 bar = "".join(bar) 257 lostr = "[%g"%lo if numpy.isfinite(lo) else "(-inf" 258 histr = "%g]"%hi if numpy.isfinite(hi) else "inf)" 259 valstr = format(self.value, self.stderr) 260 return "%25s %s %s in %s,%s" % (self.name,bar,valstr,lostr,histr)

261 - def __repr__(self):

262 return "FitParameter('%s')"%self.name

263

264 -class FitResult(object):

265 """ 266 Container for reporting fit results. 267 """

268 - def __init__(self, parameters, calls, fitness):

269 self.parameters = parameters 270 """Fit parameter list, each with name, range and value attributes.""" 271 self.calls = calls 272 """Number of function calls""" 273 self.fitness = fitness 274 """Value of the goodness of fit metric""" 275 self.pvec = numpy.array([p.value for p in self.parameters]) 276 """Parameter vector""" 277 self.stderr = None 278 """Parameter uncertainties""" 279 self.cov = None 280 """Covariance matrix"""

281

282 - def calc_cov(self, fn):

283 """ 284 Return the covariance matrix inv(J'J) at point p. 285 """ 286 if hasattr(fn, 'jacobian'): 287 # Find cov of f at p 288 # cov(f,p) = inv(J'J) 289 # Use SVD 290 # J = U S V' 291 # J'J = (U S V')' (U S V') 292 # = V S' U' U S V' 293 # = V S S V' 294 # inv(J'J) = inv(V S S V') 295 # = inv(V') inv(S S) inv(V) 296 # = V inv (S S) V' 297 J = fn.jacobian(self.pvec) 298 u,s,vh = numpy.linalg.svd(J,0) 299 JTJinv = numpy.dot(vh.T.conj()/s**2,vh) 300 self.set_cov(JTJinv)

301

302 - def set_cov(self, cov):

303 """ 304 Return the covariance matrix inv(J'J) at point p. 305 """ 306 self.cov = cov 307 if cov is not None: 308 self.stderr = numpy.sqrt(numpy.diag(self.cov)) 309 # Set the uncertainties on the individual parameters 310 for k,p in enumerate(self.parameters): 311 p.stderr = self.stderr[k] 312 else: 313 self.stderr = None 314 # Reset the uncertainties on the individual parameters 315 for k,p in enumerate(self.parameters): 316 p.stderr = None

317

318 - def print_summary(self, fid=sys.stdout):

319 if self.parameters == None: return 320 for n,p in enumerate(self.parameters): 321 print >>fid, "P%-3d %s"%(n+1,p.summarize()) 322 print >>fid, "=== goodness of fit: %g"%(self.fitness)

323

Source Code for Module park.fitting.fitresult