source: sasview/park_integration/src/sans/fit/ParkFitting.py @ 91a43d0

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 91a43d0 was 1fa4db39, checked in by Jae Cho <jhjcho@…>, 13 years ago

added comments

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