source: sasview/park_integration/src/sans/fit/ParkFitting.py @ 31469d50

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 31469d50 was 1cff677, checked in by Gervaise Alina <gervyh@…>, 13 years ago

working on get thing data and model from result of fit

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