[aa36f96] | 1 | |
---|
| 2 | |
---|
| 3 | |
---|
[792db7d5] | 4 | """ |
---|
[aa36f96] | 5 | ParkFitting module contains SansParameter,Model,Data |
---|
[b2f25dc5] | 6 | FitArrange, ParkFit,Parameter classes.All listed classes work together |
---|
| 7 | to perform a simple fit with park optimizer. |
---|
[792db7d5] | 8 | """ |
---|
[89f3b66] | 9 | #import time |
---|
[dfd3577] | 10 | import numpy |
---|
[c7e8898] | 11 | import math |
---|
[d59cb8c] | 12 | from numpy.linalg.linalg import LinAlgError |
---|
[b2f25dc5] | 13 | #import park |
---|
| 14 | from park import fit |
---|
| 15 | from park import fitresult |
---|
[1cff677] | 16 | from park.fitresult import FitParameter |
---|
| 17 | import park.simplex |
---|
[b2f25dc5] | 18 | from park.assembly import Assembly |
---|
[1cff677] | 19 | from park.assembly import Part |
---|
[b2f25dc5] | 20 | from park.fitmc import FitSimplex |
---|
[dfd3577] | 21 | import park.fitmc |
---|
[b2f25dc5] | 22 | from park.fitmc import FitMC |
---|
[1cff677] | 23 | from park.fit import Fitter |
---|
| 24 | from park.formatnum import format_uncertainty |
---|
[61cb28d] | 25 | #from Loader import Load |
---|
[b2f25dc5] | 26 | from sans.fit.AbstractFitEngine import FitEngine |
---|
[1cff677] | 27 | |
---|
[634ca14] | 28 | class SansFitResult(fitresult.FitResult): |
---|
| 29 | def __init__(self, *args, **kwrds): |
---|
| 30 | fitresult.FitResult.__init__(self, *args, **kwrds) |
---|
[425e49ca] | 31 | self.theory = None |
---|
[634ca14] | 32 | self.inputs = [] |
---|
| 33 | |
---|
[1cff677] | 34 | class SansFitSimplex(FitSimplex): |
---|
| 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 = 1 |
---|
| 46 | #xtol = 1e-4 |
---|
| 47 | """Stop when simplex vertices are within xtol of each other""" |
---|
[425e49ca] | 48 | ftol = 5e-5 |
---|
[1cff677] | 49 | """Stop when vertex values are within ftol of each other""" |
---|
| 50 | maxiter = None |
---|
| 51 | """Maximum number of iterations before fit terminates""" |
---|
[d59cb8c] | 52 | def __init__(self, ftol=5e-5): |
---|
| 53 | self.ftol = ftol |
---|
| 54 | |
---|
[1cff677] | 55 | def fit(self, fitness, x0): |
---|
| 56 | """Run the fit""" |
---|
| 57 | self.cancel = False |
---|
| 58 | pars = fitness.fit_parameters() |
---|
| 59 | bounds = numpy.array([p.range for p in pars]).T |
---|
| 60 | result = park.simplex.simplex(fitness, x0, bounds=bounds, |
---|
| 61 | radius=self.radius, xtol=self.xtol, |
---|
| 62 | ftol=self.ftol, maxiter=self.maxiter, |
---|
| 63 | abort_test=self._iscancelled) |
---|
| 64 | #print "calls:",result.calls |
---|
| 65 | #print "simplex returned",result.x,result.fx |
---|
| 66 | # Need to make our own copy of the fit results so that the |
---|
| 67 | # values don't get stomped on by the next fit iteration. |
---|
| 68 | fitpars = [SansFitParameter(pars[i].name,pars[i].range,v, pars[i].model, pars[i].data) |
---|
| 69 | for i,v in enumerate(result.x)] |
---|
[634ca14] | 70 | res = SansFitResult(fitpars, result.calls, result.fx) |
---|
| 71 | res.inputs = [(pars[i].model, pars[i].data) for i,v in enumerate(result.x)] |
---|
[1cff677] | 72 | # Compute the parameter uncertainties from the jacobian |
---|
| 73 | res.calc_cov(fitness) |
---|
| 74 | return res |
---|
| 75 | |
---|
| 76 | class SansFitter(Fitter): |
---|
| 77 | """ |
---|
| 78 | """ |
---|
| 79 | def fit(self, fitness, handler): |
---|
| 80 | """ |
---|
| 81 | Global optimizer. |
---|
| 82 | |
---|
| 83 | This function should return immediately |
---|
| 84 | """ |
---|
| 85 | # Determine initial value and bounds |
---|
| 86 | pars = fitness.fit_parameters() |
---|
| 87 | bounds = numpy.array([p.range for p in pars]).T |
---|
| 88 | x0 = [p.value for p in pars] |
---|
| 89 | |
---|
| 90 | # Initialize the monitor and results. |
---|
| 91 | # Need to make our own copy of the fit results so that the |
---|
| 92 | # values don't get stomped on by the next fit iteration. |
---|
| 93 | handler.done = False |
---|
| 94 | self.handler = handler |
---|
| 95 | fitpars = [SansFitParameter(pars[i].name, pars[i].range, v, |
---|
| 96 | pars[i].model, pars[i].data) |
---|
| 97 | for i,v in enumerate(x0)] |
---|
| 98 | handler.result = fitresult.FitResult(fitpars, 0, numpy.NaN) |
---|
[d4b0687] | 99 | |
---|
[1cff677] | 100 | # Run the fit (fit should perform _progress and _improvement updates) |
---|
| 101 | # This function may return before the fit is complete. |
---|
| 102 | self._fit(fitness, x0, bounds) |
---|
| 103 | |
---|
| 104 | class SansFitMC(SansFitter): |
---|
| 105 | """ |
---|
| 106 | Monte Carlo optimizer. |
---|
| 107 | |
---|
| 108 | This implements `park.fit.Fitter`. |
---|
| 109 | """ |
---|
| 110 | localfit = SansFitSimplex() |
---|
| 111 | start_points = 10 |
---|
[d59cb8c] | 112 | def __init__(self, localfit, start_points=10): |
---|
| 113 | self.localfit = localfit |
---|
| 114 | self.start_points = start_points |
---|
| 115 | |
---|
[1cff677] | 116 | def _fit(self, objective, x0, bounds): |
---|
| 117 | """ |
---|
| 118 | Run a monte carlo fit. |
---|
| 119 | |
---|
| 120 | This procedure maps a local optimizer across a set of initial points. |
---|
| 121 | """ |
---|
| 122 | park.fitmc.fitmc(objective, x0, bounds, self.localfit, |
---|
| 123 | self.start_points, self.handler) |
---|
| 124 | |
---|
| 125 | |
---|
| 126 | class SansPart(Part): |
---|
| 127 | """ |
---|
| 128 | Part of a fitting assembly. Part holds the model itself and |
---|
| 129 | associated data. The part can be initialized with a fitness |
---|
| 130 | object or with a pair (model,data) for the default fitness function. |
---|
[fadea71] | 131 | |
---|
[1cff677] | 132 | fitness (Fitness) |
---|
| 133 | object implementing the `park.assembly.Fitness` interface. In |
---|
| 134 | particular, fitness should provide a parameterset attribute |
---|
| 135 | containing a ParameterSet and a residuals method returning a vector |
---|
| 136 | of residuals. |
---|
| 137 | weight (dimensionless) |
---|
| 138 | weight for the model. See comments in assembly.py for details. |
---|
| 139 | isfitted (boolean) |
---|
| 140 | True if the model residuals should be included in the fit. |
---|
| 141 | The model parameters may still be used in parameter |
---|
| 142 | expressions, but there will be no comparison to the data. |
---|
| 143 | residuals (vector) |
---|
| 144 | Residuals for the model if they have been calculated, or None |
---|
| 145 | degrees_of_freedom |
---|
| 146 | Number of residuals minus number of fitted parameters. |
---|
| 147 | Degrees of freedom for individual models does not make |
---|
| 148 | sense in the presence of expressions combining models, |
---|
| 149 | particularly in the case where a model has many parameters |
---|
| 150 | but no data or many computed parameters. The degrees of |
---|
| 151 | freedom for the model is set to be at least one. |
---|
| 152 | chisq |
---|
| 153 | sum(residuals**2); use chisq/degrees_of_freedom to |
---|
| 154 | get the reduced chisq value. |
---|
[dfd3577] | 155 | |
---|
[1cff677] | 156 | Get/set the weight on the given model. |
---|
| 157 | |
---|
| 158 | assembly.weight(3) returns the weight on model 3 (0-origin) |
---|
| 159 | assembly.weight(3,0.5) sets the weight on model 3 (0-origin) |
---|
| 160 | """ |
---|
| 161 | |
---|
| 162 | def __init__(self, fitness, weight=1., isfitted=True): |
---|
| 163 | Part.__init__(self, fitness=fitness, weight=weight, |
---|
| 164 | isfitted=isfitted) |
---|
| 165 | |
---|
| 166 | self.model, self.data = fitness[0], fitness[1] |
---|
| 167 | |
---|
| 168 | class SansFitParameter(FitParameter): |
---|
| 169 | """ |
---|
| 170 | Fit result for an individual parameter. |
---|
| 171 | """ |
---|
| 172 | def __init__(self, name, range, value, model, data): |
---|
| 173 | FitParameter.__init__(self, name, range, value) |
---|
| 174 | self.model = model |
---|
| 175 | self.data = data |
---|
| 176 | |
---|
| 177 | def summarize(self): |
---|
| 178 | """ |
---|
| 179 | Return parameter range string. |
---|
| 180 | |
---|
| 181 | E.g., " Gold .....|.... 5.2043 in [2,7]" |
---|
| 182 | """ |
---|
| 183 | bar = ['.']*10 |
---|
| 184 | lo,hi = self.range |
---|
| 185 | if numpy.isfinite(lo)and numpy.isfinite(hi): |
---|
| 186 | portion = (self.value-lo)/(hi-lo) |
---|
| 187 | if portion < 0: portion = 0. |
---|
| 188 | elif portion >= 1: portion = 0.99999999 |
---|
| 189 | barpos = int(math.floor(portion*len(bar))) |
---|
| 190 | bar[barpos] = '|' |
---|
| 191 | bar = "".join(bar) |
---|
| 192 | lostr = "[%g"%lo if numpy.isfinite(lo) else "(-inf" |
---|
| 193 | histr = "%g]"%hi if numpy.isfinite(hi) else "inf)" |
---|
| 194 | valstr = format_uncertainty(self.value, self.stderr) |
---|
| 195 | model_name = str(None) |
---|
| 196 | if self.model is not None: |
---|
| 197 | model_name = self.model.name |
---|
| 198 | data_name = str(None) |
---|
| 199 | if self.data is not None: |
---|
| 200 | data_name = self.data.name |
---|
| 201 | |
---|
| 202 | return "%25s %s %s in %s,%s, %s, %s" % (self.name,bar,valstr,lostr,histr, |
---|
| 203 | model_name, data_name) |
---|
| 204 | def __repr__(self): |
---|
| 205 | #return "FitParameter('%s')"%self.name |
---|
| 206 | return str(self.__class__) |
---|
| 207 | |
---|
[dfd3577] | 208 | class MyAssembly(Assembly): |
---|
| 209 | def __init__(self, models, curr_thread=None): |
---|
| 210 | Assembly.__init__(self, models) |
---|
| 211 | self.curr_thread = curr_thread |
---|
[bd5ac39] | 212 | self.chisq = None |
---|
| 213 | self._cancel = False |
---|
[e03a14e] | 214 | self.theory = None |
---|
[dfd3577] | 215 | |
---|
[1cff677] | 216 | def fit_parameters(self): |
---|
| 217 | """ |
---|
| 218 | Return an alphabetical list of the fitting parameters. |
---|
| 219 | |
---|
| 220 | This function is called once at the beginning of a fit, |
---|
| 221 | and serves as a convenient place to precalculate what |
---|
| 222 | can be precalculated such as the set of fitting parameters |
---|
| 223 | and the parameter expressions evaluator. |
---|
| 224 | """ |
---|
| 225 | self.parameterset.setprefix() |
---|
| 226 | self._fitparameters = self.parameterset.fitted |
---|
| 227 | self._restraints = self.parameterset.restrained |
---|
| 228 | pars = self.parameterset.flatten() |
---|
| 229 | context = self.parameterset.gather_context() |
---|
| 230 | self._fitexpression = park.expression.build_eval(pars,context) |
---|
| 231 | #print "constraints",self._fitexpression.__doc__ |
---|
| 232 | |
---|
[cef2da8] | 233 | self._fitparameters.sort(lambda a,b: cmp(a.path,b.path)) |
---|
[1cff677] | 234 | # Convert to fitparameter a object |
---|
| 235 | |
---|
| 236 | fitpars = [SansFitParameter(p.path,p.range,p.value, p.model, p.data) |
---|
| 237 | for p in self._fitparameters] |
---|
| 238 | #print "fitpars", fitpars |
---|
| 239 | return fitpars |
---|
| 240 | |
---|
| 241 | def all_results(self, result): |
---|
| 242 | """ |
---|
| 243 | Extend result from the fit with the calculated parameters. |
---|
| 244 | """ |
---|
| 245 | calcpars = [SansFitParameter(p.path,p.range,p.value, p.model, p.data) |
---|
| 246 | for p in self.parameterset.computed] |
---|
| 247 | result.parameters += calcpars |
---|
[e03a14e] | 248 | result.theory = self.theory |
---|
[1cff677] | 249 | |
---|
[dfd3577] | 250 | def eval(self): |
---|
| 251 | """ |
---|
| 252 | Recalculate the theory functions, and from them, the |
---|
| 253 | residuals and chisq. |
---|
| 254 | |
---|
| 255 | :note: Call this after the parameters have been updated. |
---|
| 256 | """ |
---|
| 257 | # Handle abort from a separate thread. |
---|
| 258 | self._cancel = False |
---|
| 259 | if self.curr_thread != None: |
---|
| 260 | try: |
---|
| 261 | self.curr_thread.isquit() |
---|
| 262 | except: |
---|
| 263 | self._cancel = True |
---|
| 264 | |
---|
| 265 | # Evaluate the computed parameters |
---|
[bd5ac39] | 266 | try: |
---|
| 267 | self._fitexpression() |
---|
| 268 | except NameError: |
---|
| 269 | pass |
---|
[dfd3577] | 270 | |
---|
| 271 | # Check that the resulting parameters are in a feasible region. |
---|
| 272 | if not self.isfeasible(): return numpy.inf |
---|
| 273 | |
---|
| 274 | resid = [] |
---|
| 275 | k = len(self._fitparameters) |
---|
| 276 | for m in self.parts: |
---|
| 277 | # In order to support abort, need to be able to propagate an |
---|
| 278 | # external abort signal from self.abort() into an abort signal |
---|
| 279 | # for the particular model. Can't see a way to do this which |
---|
| 280 | # doesn't involve setting a state variable. |
---|
| 281 | self._current_model = m |
---|
| 282 | if self._cancel: return numpy.inf |
---|
| 283 | if m.isfitted and m.weight != 0: |
---|
[e03a14e] | 284 | m.residuals, self.theory = m.fitness.residuals() |
---|
[dfd3577] | 285 | N = len(m.residuals) |
---|
| 286 | m.degrees_of_freedom = N-k if N>k else 1 |
---|
[1fa4db39] | 287 | # dividing residuals by N in order to be consistent with Scipy |
---|
[c7e8898] | 288 | m.chisq = numpy.sum(m.residuals**2/N) |
---|
| 289 | resid.append(m.weight*m.residuals/math.sqrt(N)) |
---|
[dfd3577] | 290 | self.residuals = numpy.hstack(resid) |
---|
| 291 | N = len(self.residuals) |
---|
| 292 | self.degrees_of_freedom = N-k if N>k else 1 |
---|
| 293 | self.chisq = numpy.sum(self.residuals**2) |
---|
| 294 | return self.chisq |
---|
| 295 | |
---|
[4c718654] | 296 | class ParkFit(FitEngine): |
---|
[7705306] | 297 | """ |
---|
[aa36f96] | 298 | ParkFit performs the Fit.This class can be used as follow: |
---|
| 299 | #Do the fit Park |
---|
| 300 | create an engine: engine = ParkFit() |
---|
| 301 | Use data must be of type plottable |
---|
| 302 | Use a sans model |
---|
| 303 | |
---|
| 304 | Add data with a dictionnary of FitArrangeList where Uid is a key and data |
---|
| 305 | is saved in FitArrange object. |
---|
| 306 | engine.set_data(data,Uid) |
---|
| 307 | |
---|
| 308 | Set model parameter "M1"= model.name add {model.parameter.name:value}. |
---|
| 309 | |
---|
| 310 | :note: Set_param() if used must always preceded set_model() |
---|
| 311 | for the fit to be performed. |
---|
| 312 | engine.set_param( model,"M1", {'A':2,'B':4}) |
---|
| 313 | |
---|
[b2f25dc5] | 314 | Add model with a dictionnary of FitArrangeList{} where Uid is a key |
---|
| 315 | and model |
---|
[aa36f96] | 316 | is save in FitArrange object. |
---|
| 317 | engine.set_model(model,Uid) |
---|
| 318 | |
---|
| 319 | engine.fit return chisqr,[model.parameter 1,2,..],[[err1....][..err2...]] |
---|
| 320 | chisqr1, out1, cov1=engine.fit({model.parameter.name:value},qmin,qmax) |
---|
| 321 | |
---|
| 322 | :note: {model.parameter.name:value} is ignored in fit function since |
---|
[792db7d5] | 323 | the user should make sure to call set_param himself. |
---|
[aa36f96] | 324 | |
---|
[7705306] | 325 | """ |
---|
[916a15f] | 326 | def __init__(self): |
---|
[792db7d5] | 327 | """ |
---|
[aa36f96] | 328 | Creates a dictionary (self.fitArrangeList={})of FitArrange elements |
---|
| 329 | with Uid as keys |
---|
[792db7d5] | 330 | """ |
---|
[b2f25dc5] | 331 | FitEngine.__init__(self) |
---|
| 332 | self.fit_arrange_dict = {} |
---|
| 333 | self.param_list = [] |
---|
[37d9521] | 334 | |
---|
[dfd3577] | 335 | def create_assembly(self, curr_thread): |
---|
[7705306] | 336 | """ |
---|
[b2f25dc5] | 337 | Extract sansmodel and sansdata from |
---|
| 338 | self.FitArrangelist ={Uid:FitArrange} |
---|
| 339 | Create parkmodel and park data ,form a list couple of parkmodel |
---|
| 340 | and parkdata |
---|
[792db7d5] | 341 | create an assembly self.problem= park.Assembly([(parkmodel,parkdata)]) |
---|
[7705306] | 342 | """ |
---|
[89f3b66] | 343 | mylist = [] |
---|
[b2f25dc5] | 344 | #listmodel = [] |
---|
| 345 | #i = 0 |
---|
[89f3b66] | 346 | fitproblems = [] |
---|
[b2f25dc5] | 347 | for fproblem in self.fit_arrange_dict.itervalues(): |
---|
[89f3b66] | 348 | if fproblem.get_to_fit() == 1: |
---|
[a9e04aa] | 349 | fitproblems.append(fproblem) |
---|
[89f3b66] | 350 | if len(fitproblems) == 0: |
---|
[a9e04aa] | 351 | raise RuntimeError, "No Assembly scheduled for Park fitting." |
---|
| 352 | return |
---|
| 353 | for item in fitproblems: |
---|
| 354 | parkmodel = item.get_model() |
---|
[9e85792] | 355 | for p in parkmodel.parameterset: |
---|
[fe886ee] | 356 | ## does not allow status change for constraint parameters |
---|
[89f3b66] | 357 | if p.status != 'computed': |
---|
[b2f25dc5] | 358 | if p.get_name()in item.pars: |
---|
| 359 | ## make parameters selected for |
---|
| 360 | #fit will be between boundaries |
---|
[89f3b66] | 361 | p.set(p.range) |
---|
[aed7c57] | 362 | else: |
---|
[89f3b66] | 363 | p.status = 'fixed' |
---|
[b2f25dc5] | 364 | data_list = item.get_data() |
---|
| 365 | parkdata = data_list |
---|
[89f3b66] | 366 | fitness = (parkmodel, parkdata) |
---|
[ca6d914] | 367 | mylist.append(fitness) |
---|
[dfd3577] | 368 | self.problem = MyAssembly(models=mylist, curr_thread=curr_thread) |
---|
[126a761] | 369 | |
---|
[dfd3577] | 370 | |
---|
[93de635d] | 371 | def fit(self, q=None, handler=None, curr_thread=None, ftol=1.49012e-8): |
---|
[7705306] | 372 | """ |
---|
[aa36f96] | 373 | Performs fit with park.fit module.It can perform fit with one model |
---|
| 374 | and a set of data, more than two fit of one model and sets of data or |
---|
[b2f25dc5] | 375 | fit with more than two model associated with their set of data and |
---|
| 376 | constraints |
---|
[aa36f96] | 377 | |
---|
[b2f25dc5] | 378 | :param pars: Dictionary of parameter names for the model and their |
---|
| 379 | values. |
---|
[aa36f96] | 380 | :param qmin: The minimum value of data's range to be fit |
---|
| 381 | :param qmax: The maximum value of data's range to be fit |
---|
| 382 | |
---|
[b2f25dc5] | 383 | :note: all parameter are ignored most of the time.Are just there |
---|
| 384 | to keep ScipyFit and ParkFit interface the same. |
---|
[aa36f96] | 385 | |
---|
| 386 | :return: result.fitness Value of the goodness of fit metric |
---|
[b2f25dc5] | 387 | :return: result.pvec list of parameter with the best value |
---|
| 388 | found during fitting |
---|
[aa36f96] | 389 | :return: result.cov Covariance matrix |
---|
| 390 | |
---|
[7705306] | 391 | """ |
---|
[dfd3577] | 392 | self.create_assembly(curr_thread=curr_thread) |
---|
[1cff677] | 393 | localfit = SansFitSimplex() |
---|
[93de635d] | 394 | localfit.ftol = ftol |
---|
[681f0dc] | 395 | |
---|
[916a15f] | 396 | # See `park.fitresult.FitHandler` for details. |
---|
[1cff677] | 397 | fitter = SansFitMC(localfit=localfit, start_points=1) |
---|
[681f0dc] | 398 | if handler == None: |
---|
[89f3b66] | 399 | handler = fitresult.ConsoleUpdate(improvement_delta=0.1) |
---|
[d59cb8c] | 400 | try: |
---|
| 401 | result = fit.fit(self.problem, fitter=fitter, handler=handler) |
---|
| 402 | self.problem.all_results(result) |
---|
| 403 | except LinAlgError: |
---|
| 404 | raise ValueError, "SVD did not converge" |
---|
| 405 | |
---|
[89f3b66] | 406 | if result != None: |
---|
| 407 | if q != None: |
---|
[fd6b789] | 408 | q.put(result) |
---|
| 409 | return q |
---|
[48882d1] | 410 | return result |
---|
[ee5b04c] | 411 | else: |
---|
| 412 | raise ValueError, "SVD did not converge" |
---|
[aa36f96] | 413 | |
---|