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

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 cc694d0 was ba7dceb, checked in by Gervaise Alina <gervyh@…>, 13 years ago

prepare fitting engine for multiprocessing

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