source: sasview/src/sas/fit/ParkFitting.py @ 386ffe1

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since 386ffe1 was 386ffe1, checked in by pkienzle, 9 years ago

remove scipy levenburg marquardt and park from ui

  • Property mode set to 100644
File size: 21.7 KB
Line 
1
2
3
4"""
5ParkFitting module contains SasParameter,Model,Data
6FitArrange, ParkFit,Parameter classes.All listed classes work together
7to perform a simple fit with park optimizer.
8"""
9
10_ = '''
11#import time
12import numpy
13import math
14from  numpy.linalg.linalg import LinAlgError
15#import park
16from park import fit
17from park import fitresult
18from  park.fitresult import FitParameter
19import park.simplex
20from park.assembly import Assembly
21from park.assembly import Part
22from park.fitmc import FitSimplex
23import park.fitmc
24from park.fit import Fitter
25from park.formatnum import format_uncertainty
26from sas.fit.AbstractFitEngine import FitEngine
27from sas.fit.AbstractFitEngine import FResult
28
29class SasParameter(park.Parameter):
30    """
31    SAS model parameters for use in the PARK fitting service.
32    The parameter attribute value is redirected to the underlying
33    parameter value in the SAS model.
34    """
35    def __init__(self, name, model, data):
36        """
37            :param name: the name of the model parameter
38            :param model: the sas model to wrap as a park model
39        """
40        park.Parameter.__init__(self, name)
41        #self._model, self._name = model, name
42        self.data = data
43        self.model = model
44        #set the value for the parameter of the given name
45        self.set(model.getParam(name))
46
47        # TODO: model is missing parameter ranges for dispersion parameters
48        if name not in model.details:
49            #print "setting details for",name
50            model.details[name] = ["", None, None]
51
52    def _getvalue(self):
53        """
54        override the _getvalue of park parameter
55
56        :return value the parameter associates with self.name
57
58        """
59        return self.model.getParam(self.name)
60
61    def _setvalue(self, value):
62        """
63        override the _setvalue pf park parameter
64
65        :param value: the value to set on a given parameter
66
67        """
68        self.model.setParam(self.name, value)
69
70    value = property(_getvalue, _setvalue)
71
72    def _getrange(self):
73        """
74        Override _getrange of park parameter
75        return the range of parameter
76        """
77        #if not  self.name in self._model.getDispParamList():
78        lo, hi = self.model.details[self.name][1:3]
79        if lo is None: lo = -numpy.inf
80        if hi is None: hi = numpy.inf
81        if lo > hi:
82            raise ValueError, "wrong fit range for parameters"
83
84        return lo, hi
85
86    def get_name(self):
87        """
88        """
89        return self._getname()
90
91    def _setrange(self, r):
92        """
93        override _setrange of park parameter
94
95        :param r: the value of the range to set
96
97        """
98        self.model.details[self.name][1:3] = r
99    range = property(_getrange, _setrange)
100
101
102class ParkModel(park.Model):
103    """
104    PARK wrapper for SAS models.
105    """
106    def __init__(self, sas_model, sas_data=None, **kw):
107        """
108        :param sas_model: the sas model to wrap using park interface
109
110        """
111        park.Model.__init__(self, **kw)
112        self.model = sas_model
113        self.name = sas_model.name
114        self.data = sas_data
115        #list of parameters names
116        self.sasp = sas_model.getParamList()
117        #list of park parameter
118        self.parkp = [SasParameter(p, sas_model, sas_data) for p in self.sasp]
119        #list of parameter set
120        self.parameterset = park.ParameterSet(sas_model.name, pars=self.parkp)
121        self.pars = []
122
123    def get_params(self, fitparams):
124        """
125        return a list of value of paramter to fit
126
127        :param fitparams: list of paramaters name to fit
128
129        """
130        list_params = []
131        self.pars = fitparams
132        for item in fitparams:
133            for element in self.parkp:
134                if element.name == str(item):
135                    list_params.append(element.value)
136        return list_params
137
138    def set_params(self, paramlist, params):
139        """
140        Set value for parameters to fit
141
142        :param params: list of value for parameters to fit
143
144        """
145        try:
146            for i in range(len(self.parkp)):
147                for j in range(len(paramlist)):
148                    if self.parkp[i].name == paramlist[j]:
149                        self.parkp[i].value = params[j]
150                        self.model.setParam(self.parkp[i].name, params[j])
151        except:
152            raise
153
154    def eval(self, x):
155        """
156            Override eval method of park model.
157
158            :param x: the x value used to compute a function
159        """
160        try:
161            return self.model.evalDistribution(x)
162        except:
163            raise
164
165    def eval_derivs(self, x, pars=[]):
166        """
167        Evaluate the model and derivatives wrt pars at x.
168
169        pars is a list of the names of the parameters for which derivatives
170        are desired.
171
172        This method needs to be specialized in the model to evaluate the
173        model function.  Alternatively, the model can implement is own
174        version of residuals which calculates the residuals directly
175        instead of calling eval.
176        """
177        return []
178
179
180class SasFitResult(fitresult.FitResult):
181    def __init__(self, *args, **kwrds):
182        fitresult.FitResult.__init__(self, *args, **kwrds)
183        self.theory = None
184        self.inputs = []
185       
186class SasFitSimplex(FitSimplex):
187    """
188    Local minimizer using Nelder-Mead simplex algorithm.
189
190    Simplex is robust and derivative free, though not very efficient.
191
192    This class wraps the bounds contrained Nelder-Mead simplex
193    implementation for `park.simplex.simplex`.
194    """
195    radius = 0.05
196    """Size of the initial simplex; this is a portion between 0 and 1"""
197    xtol = 1
198    #xtol = 1e-4
199    """Stop when simplex vertices are within xtol of each other"""
200    ftol = 5e-5
201    """Stop when vertex values are within ftol of each other"""
202    maxiter = None
203    """Maximum number of iterations before fit terminates"""
204    def __init__(self, ftol=5e-5):
205        self.ftol = ftol
206       
207    def fit(self, fitness, x0):
208        """Run the fit"""
209        self.cancel = False
210        pars = fitness.fit_parameters()
211        bounds = numpy.array([p.range for p in pars]).T
212        result = park.simplex.simplex(fitness, x0, bounds=bounds,
213                                 radius=self.radius, xtol=self.xtol,
214                                 ftol=self.ftol, maxiter=self.maxiter,
215                                 abort_test=self._iscancelled)
216        #print "calls:",result.calls
217        #print "simplex returned",result.x,result.fx
218        # Need to make our own copy of the fit results so that the
219        # values don't get stomped on by the next fit iteration.
220        fitpars = [SasFitParameter(pars[i].name,pars[i].range,v, pars[i].model, pars[i].data)
221                   for i,v in enumerate(result.x)]
222        res = SasFitResult(fitpars, result.calls, result.fx)
223        res.inputs = [(pars[i].model, pars[i].data) for i,v in enumerate(result.x)]
224        # Compute the parameter uncertainties from the jacobian
225        res.calc_cov(fitness)
226        return res
227     
228class SasFitter(Fitter):
229    """
230    """
231    def fit(self, fitness, handler):
232        """
233        Global optimizer.
234
235        This function should return immediately
236        """
237        # Determine initial value and bounds
238        pars = fitness.fit_parameters()
239        bounds = numpy.array([p.range for p in pars]).T
240        x0 = [p.value for p in pars]
241
242        # Initialize the monitor and results.
243        # Need to make our own copy of the fit results so that the
244        # values don't get stomped on by the next fit iteration.
245        handler.done = False
246        self.handler = handler
247        fitpars = [SasFitParameter(pars[i].name, pars[i].range, v,
248                                     pars[i].model, pars[i].data)
249                   for i,v in enumerate(x0)]
250        handler.result = fitresult.FitResult(fitpars, 0, numpy.NaN)
251
252        # Run the fit (fit should perform _progress and _improvement updates)
253        # This function may return before the fit is complete.
254        self._fit(fitness, x0, bounds)
255       
256class SasFitMC(SasFitter):
257    """
258    Monte Carlo optimizer.
259
260    This implements `park.fit.Fitter`.
261    """
262    localfit = SasFitSimplex()
263    start_points = 10
264    def __init__(self, localfit, start_points=10):
265        self.localfit = localfit
266        self.start_points = start_points
267       
268    def _fit(self, objective, x0, bounds):
269        """
270        Run a monte carlo fit.
271
272        This procedure maps a local optimizer across a set of initial points.
273        """
274        try:
275            park.fitmc.fitmc(objective, x0, bounds, self.localfit,
276                             self.start_points, self.handler)
277        except:
278            raise ValueError, "Fit did not converge.\n"
279       
280class SasPart(Part):
281    """
282    Part of a fitting assembly.  Part holds the model itself and
283    associated data.  The part can be initialized with a fitness
284    object or with a pair (model,data) for the default fitness function.
285
286    fitness (Fitness)
287        object implementing the `park.assembly.Fitness` interface.  In
288        particular, fitness should provide a parameterset attribute
289        containing a ParameterSet and a residuals method returning a vector
290        of residuals.
291    weight (dimensionless)
292        weight for the model.  See comments in assembly.py for details.
293    isfitted (boolean)
294        True if the model residuals should be included in the fit.
295        The model parameters may still be used in parameter
296        expressions, but there will be no comparison to the data.
297    residuals (vector)
298        Residuals for the model if they have been calculated, or None
299    degrees_of_freedom
300        Number of residuals minus number of fitted parameters.
301        Degrees of freedom for individual models does not make
302        sense in the presence of expressions combining models,
303        particularly in the case where a model has many parameters
304        but no data or many computed parameters.  The degrees of
305        freedom for the model is set to be at least one.
306    chisq
307        sum(residuals**2); use chisq/degrees_of_freedom to
308        get the reduced chisq value.
309
310        Get/set the weight on the given model.
311
312        assembly.weight(3) returns the weight on model 3 (0-origin)
313        assembly.weight(3,0.5) sets the weight on model 3 (0-origin)
314    """
315
316    def __init__(self, fitness, weight=1., isfitted=True):
317        Part.__init__(self, fitness=fitness, weight=weight,
318                       isfitted=isfitted)
319       
320        self.model, self.data = fitness[0], fitness[1]
321
322class SasFitParameter(FitParameter):
323    """
324    Fit result for an individual parameter.
325    """
326    def __init__(self, name, range, value, model, data):
327        FitParameter.__init__(self, name, range, value)
328        self.model = model
329        self.data = data
330       
331    def summarize(self):
332        """
333        Return parameter range string.
334
335        E.g.,  "       Gold .....|.... 5.2043 in [2,7]"
336        """
337        bar = ['.']*10
338        lo,hi = self.range
339        if numpy.isfinite(lo)and numpy.isfinite(hi):
340            portion = (self.value-lo)/(hi-lo)
341            if portion < 0: portion = 0.
342            elif portion >= 1: portion = 0.99999999
343            barpos = int(math.floor(portion*len(bar)))
344            bar[barpos] = '|'
345        bar = "".join(bar)
346        lostr = "[%g"%lo if numpy.isfinite(lo) else "(-inf"
347        histr = "%g]"%hi if numpy.isfinite(hi) else "inf)"
348        valstr = format_uncertainty(self.value, self.stderr)
349        model_name = str(None)
350        if self.model is not None:
351            model_name = self.model.name
352        data_name = str(None)
353        if self.data is not None:
354            data_name = self.data.name
355           
356        return "%25s %s %s in %s,%s, %s, %s"  % (self.name,bar,valstr,lostr,histr,
357                                                 model_name, data_name)
358    def __repr__(self):
359        #return "FitParameter('%s')"%self.name
360        return str(self.__class__)
361   
362class MyAssembly(Assembly):
363    def __init__(self, models, curr_thread=None):
364        """Build an assembly from a list of models."""
365        self.parts = []
366        for m in models:
367            self.parts.append(SasPart(m))
368        self.curr_thread = curr_thread
369        self.chisq = None
370        self._cancel = False
371        self.theory = None
372        self._reset()
373       
374    def fit_parameters(self):
375        """
376        Return an alphabetical list of the fitting parameters.
377
378        This function is called once at the beginning of a fit,
379        and serves as a convenient place to precalculate what
380        can be precalculated such as the set of fitting parameters
381        and the parameter expressions evaluator.
382        """
383        self.parameterset.setprefix()
384        self._fitparameters = self.parameterset.fitted
385        self._restraints = self.parameterset.restrained
386        pars = self.parameterset.flatten()
387        context = self.parameterset.gather_context()
388        self._fitexpression = park.expression.build_eval(pars,context)
389        #print "constraints",self._fitexpression.__doc__
390
391        self._fitparameters.sort(lambda a,b: cmp(a.path,b.path))
392        # Convert to fitparameter a object
393       
394        fitpars = [SasFitParameter(p.path,p.range,p.value, p.model, p.data)
395                   for p in self._fitparameters]
396        #print "fitpars", fitpars
397        return fitpars
398   
399    def extend_results_with_calculated_parameters(self, result):
400        """
401        Extend result from the fit with the calculated parameters.
402        """
403        calcpars = [SasFitParameter(p.path,p.range,p.value, p.model, p.data)
404                    for p in self.parameterset.computed]
405        result.parameters += calcpars
406        result.theory = self.theory
407
408    def eval(self):
409        """
410        Recalculate the theory functions, and from them, the
411        residuals and chisq.
412
413        :note: Call this after the parameters have been updated.
414        """
415        # Handle abort from a separate thread.
416        self._cancel = False
417        if self.curr_thread != None:
418            try:
419                self.curr_thread.isquit()
420            except:
421                self._cancel = True
422
423        # Evaluate the computed parameters
424        try:
425            self._fitexpression()
426        except NameError:
427            pass
428
429        # Check that the resulting parameters are in a feasible region.
430        if not self.isfeasible(): return numpy.inf
431
432        resid = []
433        k = len(self._fitparameters)
434        for m in self.parts:
435            # In order to support abort, need to be able to propagate an
436            # external abort signal from self.abort() into an abort signal
437            # for the particular model.  Can't see a way to do this which
438            # doesn't involve setting a state variable.
439            self._current_model = m
440            if self._cancel: return numpy.inf
441            if m.isfitted and m.weight != 0:
442                m.residuals, self.theory = m.fitness.residuals()
443                N = len(m.residuals)
444                m.degrees_of_freedom = N-k if N>k else 1
445                # dividing residuals by N in order to be consistent with Scipy
446                m.chisq = numpy.sum(m.residuals**2/N)
447                resid.append(m.weight*m.residuals)
448        self.residuals = numpy.hstack(resid)
449        N = len(self.residuals)
450        self.degrees_of_freedom = N-k if N>k else 1
451        self.chisq = numpy.sum(self.residuals**2)
452        return self.chisq/self.degrees_of_freedom
453   
454class ParkFit(FitEngine):
455    """
456    ParkFit performs the Fit.This class can be used as follow:
457    #Do the fit Park
458    create an engine: engine = ParkFit()
459    Use data must be of type plottable
460    Use a sas model
461   
462    Add data with a dictionnary of FitArrangeList where Uid is a key and data
463    is saved in FitArrange object.
464    engine.set_data(data,Uid)
465   
466    Set model parameter "M1"= model.name add {model.parameter.name:value}.
467   
468    ..note::
469
470       Set_param() if used must always preceded set_model() for the fit to be performed. ``engine.set_param( model,"M1", {'A':2,'B':4})``
471   
472    Add model with a dictionnary of FitArrangeList{} where Uid is a key
473    and model
474    is save in FitArrange object.
475    engine.set_model(model,Uid)
476   
477    engine.fit return chisqr,[model.parameter 1,2,..],[[err1....][..err2...]]
478    chisqr1, out1, cov1=engine.fit({model.parameter.name:value},qmin,qmax)
479   
480    ..note::
481
482        {model.parameter.name:value} is ignored in fit function since
483        the user should make sure to call set_param himself.
484       
485    """
486    def __init__(self):
487        """
488        Creates a dictionary (self.fitArrangeList={})of FitArrange elements
489        with Uid as keys
490        """
491        FitEngine.__init__(self)
492        self.fit_arrange_dict = {}
493        self.param_list = []
494       
495    def create_assembly(self, curr_thread, reset_flag=False):
496        """
497        Extract sasmodel and sasdata from
498        self.FitArrangelist ={Uid:FitArrange}
499        Create parkmodel and park data ,form a list couple of parkmodel
500        and parkdata
501        create an assembly self.problem=  park.Assembly([(parkmodel,parkdata)])
502        """
503        mylist = []
504        #listmodel = []
505        #i = 0
506        fitproblems = []
507        for fproblem in self.fit_arrange_dict.itervalues():
508            if fproblem.get_to_fit() == 1:
509                fitproblems.append(fproblem)
510        if len(fitproblems) == 0:
511            raise RuntimeError, "No Assembly scheduled for Park fitting."
512        for item in fitproblems:
513            model = item.get_model()
514            parkmodel = ParkModel(model.model, model.data)
515            parkmodel.pars = item.pars
516            if reset_flag:
517                # reset the initial value; useful for batch
518                for name in item.pars:
519                    ind = item.pars.index(name)
520                    parkmodel.model.setParam(name, item.vals[ind])
521
522            # set the constraints into the model
523            for p,v in item.constraints:
524                parkmodel.parameterset[str(p)].set(str(v))
525           
526            for p in parkmodel.parameterset:
527                ## does not allow status change for constraint parameters
528                if p.status != 'computed':
529                    if p.get_name() in item.pars:
530                        ## make parameters selected for
531                        #fit will be between boundaries
532                        p.set(p.range)         
533                    else:
534                        p.status = 'fixed'
535            data_list = item.get_data()
536            parkdata = data_list
537            fitness = (parkmodel, parkdata)
538            mylist.append(fitness)
539        self.problem = MyAssembly(models=mylist, curr_thread=curr_thread)
540       
541 
542    def fit(self, msg_q=None,
543            q=None, handler=None, curr_thread=None,
544            ftol=1.49012e-8, reset_flag=False):
545        """
546        Performs fit with park.fit module.It can  perform fit with one model
547        and a set of data, more than two fit of  one model and sets of data or
548        fit with more than two model associated with their set of data and
549        constraints
550       
551        :param pars: Dictionary of parameter names for the model and their
552            values.
553        :param qmin: The minimum value of data's range to be fit
554        :param qmax: The maximum value of data's range to be fit
555       
556        :note: all parameter are ignored most of the time.Are just there
557            to keep ScipyFit and ParkFit interface the same.
558           
559        :return: result.fitness Value of the goodness of fit metric
560        :return: result.pvec list of parameter with the best value
561            found during fitting
562        :return: result.cov Covariance matrix
563       
564        """
565        self.create_assembly(curr_thread=curr_thread, reset_flag=reset_flag)
566        localfit = SasFitSimplex()
567        localfit.ftol = ftol
568        localfit.xtol = 1e-6
569
570        # See `park.fitresult.FitHandler` for details.
571        fitter = SasFitMC(localfit=localfit, start_points=1)
572        if handler == None:
573            handler = fitresult.ConsoleUpdate(improvement_delta=0.1)
574       
575        result_list = []
576        try:
577            result = fit.fit(self.problem, fitter=fitter, handler=handler)
578            self.problem.extend_results_with_calculated_parameters(result)
579           
580        except LinAlgError:
581            raise ValueError, "SVD did not converge"
582
583        if result is None:
584            raise RuntimeError("park did not return a fit result")
585   
586        for m in self.problem.parts:
587            residuals, theory = m.fitness.residuals()
588            small_result = FResult(model=m.model, data=m.data.sas_data)
589            small_result.fitter_id = self.fitter_id
590            small_result.theory = theory
591            small_result.residuals = residuals
592            small_result.index = m.data.idx
593            small_result.fitness = result.fitness
594
595            # Extract the parameters that are part of this model; make sure
596            # they match the fitted parameters for this model, and place them
597            # in the same order as they occur in the model.
598            pars = {}
599            for p in result.parameters:
600                #if p.data.name == small_result.data.name and
601                if p.model.name == small_result.model.name:
602                    model_name, par_name = p.name.split('.', 1)
603                    pars[par_name] = (p.value, p.stderr)
604            #assert len(pars.keys()) == len(m.model.pars)
605            v,dv = zip(*[pars[p] for p in m.model.pars])
606            small_result.pvec = v
607            small_result.stderr = dv
608            small_result.param_list = m.model.pars
609
610            # normalize chisq by degrees of freedom
611            dof = len(small_result.residuals)-len(small_result.pvec)
612            small_result.fitness = numpy.sum(residuals**2)/dof
613
614            result_list.append(small_result)   
615        if q != None:
616            q.put(result_list)
617            return q
618        return result_list
619
620'''
Note: See TracBrowser for help on using the repository browser.