source: sasview/park_integration/AbstractFitEngine.py @ b345c45

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 b345c45 was 027e8f2, checked in by Jae Cho <jhjcho@…>, 15 years ago

Fixed a 2D problem of not properly (calculating chi2 and fitting 2Ddata) with bounded Q range.

  • Property mode set to 100644
File size: 22.4 KB
Line 
1import logging, sys
2import park,numpy,math, copy
3from DataLoader.data_info import Data1D
4from DataLoader.data_info import Data2D
5class SansParameter(park.Parameter):
6    """
7        SANS model parameters for use in the PARK fitting service.
8        The parameter attribute value is redirected to the underlying
9        parameter value in the SANS model.
10    """
11    def __init__(self, name, model):
12        """
13            @param name: the name of the model parameter
14            @param model: the sans model to wrap as a park model
15        """
16        self._model, self._name = model,name
17        #set the value for the parameter of the given name
18        self.set(model.getParam(name))
19         
20    def _getvalue(self):
21        """
22            override the _getvalue of park parameter
23            @return value the parameter associates with self.name
24        """
25        return self._model.getParam(self.name)
26   
27    def _setvalue(self,value):
28        """
29            override the _setvalue pf park parameter
30            @param value: the value to set on a given parameter
31        """
32        self._model.setParam(self.name, value)
33       
34    value = property(_getvalue,_setvalue)
35   
36    def _getrange(self):
37        """
38            Override _getrange of park parameter
39            return the range of parameter
40        """
41        #if not  self.name in self._model.getDispParamList():
42        lo,hi = self._model.details[self.name][1:3]
43        if lo is None: lo = -numpy.inf
44        if hi is None: hi = numpy.inf
45        #else:
46            #lo,hi = self._model.details[self.name][1:]
47            #if lo is None: lo = -numpy.inf
48            #if hi is None: hi = numpy.inf
49        if lo >= hi:
50            raise ValueError,"wrong fit range for parameters"
51       
52        return lo,hi
53   
54    def _setrange(self,r):
55        """
56            override _setrange of park parameter
57            @param r: the value of the range to set
58        """
59        self._model.details[self.name][1:3] = r
60    range = property(_getrange,_setrange)
61   
62class Model(park.Model):
63    """
64        PARK wrapper for SANS models.
65    """
66    def __init__(self, sans_model, **kw):
67        """
68            @param sans_model: the sans model to wrap using park interface
69        """
70        park.Model.__init__(self, **kw)
71        self.model = sans_model
72        self.name = sans_model.name
73        #list of parameters names
74        self.sansp = sans_model.getParamList()
75        #list of park parameter
76        self.parkp = [SansParameter(p,sans_model) for p in self.sansp]
77        #list of parameterset
78        self.parameterset = park.ParameterSet(sans_model.name,pars=self.parkp)
79        self.pars=[]
80 
81 
82    def getParams(self,fitparams):
83        """
84            return a list of value of paramter to fit
85            @param fitparams: list of paramaters name to fit
86        """
87        list=[]
88        self.pars=[]
89        self.pars=fitparams
90        for item in fitparams:
91            for element in self.parkp:
92                 if element.name ==str(item):
93                     list.append(element.value)
94        return list
95   
96   
97    def setParams(self,paramlist, params):
98        """
99            Set value for parameters to fit
100            @param params: list of value for parameters to fit
101        """
102        try:
103            for i in range(len(self.parkp)):
104                for j in range(len(paramlist)):
105                    if self.parkp[i].name==paramlist[j]:
106                        self.parkp[i].value = params[j]
107                        self.model.setParam(self.parkp[i].name,params[j])
108        except:
109            raise
110 
111    def eval(self,x):
112        """
113            override eval method of park model.
114            @param x: the x value used to compute a function
115        """
116        try:
117                return self.model.evalDistribution(x)
118        except:
119                raise
120
121   
122class FitData1D(Data1D):
123    """
124        Wrapper class  for SANS data
125        FitData1D inherits from DataLoader.data_info.Data1D. Implements
126        a way to get residuals from data.
127    """
128    def __init__(self,x, y,dx= None, dy=None, smearer=None):
129        Data1D.__init__(self, x=numpy.array(x), y=numpy.array(y), dx=dx, dy=dy)
130        """
131            @param smearer: is an object of class QSmearer or SlitSmearer
132            that will smear the theory data (slit smearing or resolution
133            smearing) when set.
134           
135            The proper way to set the smearing object would be to
136            do the following:
137           
138            from DataLoader.qsmearing import smear_selection
139            smearer = smear_selection(some_data)
140            fitdata1d = FitData1D( x= [1,3,..,],
141                                    y= [3,4,..,8],
142                                    dx=None,
143                                    dy=[1,2...], smearer= smearer)
144           
145            Note that some_data _HAS_ to be of class DataLoader.data_info.Data1D
146           
147            Setting it back to None will turn smearing off.
148           
149        """
150       
151        self.smearer = smearer
152        if dy ==None or dy==[]:
153            self.dy= numpy.zeros(len(self.y)) 
154        else:
155            self.dy= numpy.asarray(dy)
156     
157        # For fitting purposes, replace zero errors by 1
158        #TODO: check validity for the rare case where only
159        # a few points have zero errors
160        self.dy[self.dy==0]=1
161       
162        ## Min Q-value
163        #Skip the Q=0 point, especially when y(q=0)=None at x[0].
164        if min (self.x) ==0.0 and self.x[0]==0 and not numpy.isfinite(self.y[0]):
165            self.qmin = min(self.x[self.x!=0])
166        else:                             
167            self.qmin= min (self.x)
168        ## Max Q-value
169        self.qmax = max (self.x)
170       
171        # Range used for input to smearing
172        self._qmin_unsmeared = self.qmin
173        self._qmax_unsmeared = self.qmax
174        # Identify the bin range for the unsmeared and smeared spaces
175        self.idx = (self.x>=self.qmin) & (self.x <= self.qmax)
176        self.idx_unsmeared = (self.x>=self._qmin_unsmeared) & (self.x <= self._qmax_unsmeared)
177 
178       
179       
180    def setFitRange(self,qmin=None,qmax=None):
181        """ to set the fit range"""
182        # Skip Q=0 point, (especially for y(q=0)=None at x[0]).
183        #ToDo: Fix this.
184        if qmin==0.0 and not numpy.isfinite(self.y[qmin]):
185            self.qmin = min(self.x[self.x!=0])
186        elif qmin!=None:                       
187            self.qmin = qmin           
188
189        if qmax !=None:
190            self.qmax = qmax
191           
192        # Determine the range needed in unsmeared-Q to cover
193        # the smeared Q range
194        self._qmin_unsmeared = self.qmin
195        self._qmax_unsmeared = self.qmax   
196       
197        self._first_unsmeared_bin = 0
198        self._last_unsmeared_bin  = len(self.x)-1
199       
200        if self.smearer!=None:
201            self._first_unsmeared_bin, self._last_unsmeared_bin = self.smearer.get_bin_range(self.qmin, self.qmax)
202            self._qmin_unsmeared = self.x[self._first_unsmeared_bin]
203            self._qmax_unsmeared = self.x[self._last_unsmeared_bin]
204           
205        # Identify the bin range for the unsmeared and smeared spaces
206        self.idx = (self.x>=self.qmin) & (self.x <= self.qmax)
207        self.idx_unsmeared = (self.x>=self._qmin_unsmeared) & (self.x <= self._qmax_unsmeared)
208 
209       
210    def getFitRange(self):
211        """
212            @return the range of data.x to fit
213        """
214        return self.qmin, self.qmax
215       
216    def residuals(self, fn):
217        """
218            Compute residuals.
219           
220            If self.smearer has been set, use if to smear
221            the data before computing chi squared.
222           
223            @param fn: function that return model value
224            @return residuals
225        """
226        # Compute theory data f(x)
227        fx= numpy.zeros(len(self.x))
228        fx[self.idx_unsmeared] = fn(self.x[self.idx_unsmeared])
229       
230        ## Smear theory data
231        if self.smearer is not None:
232            fx = self.smearer(fx, self._first_unsmeared_bin, self._last_unsmeared_bin)
233       
234        ## Sanity check
235        if numpy.size(self.dy)!= numpy.size(fx):
236            raise RuntimeError, "FitData1D: invalid error array %d <> %d" % (numpy.shape(self.dy),
237                                                                              numpy.size(fx))
238                                                                             
239        return (self.y[self.idx]-fx[self.idx])/self.dy[self.idx]
240     
241 
242       
243    def residuals_deriv(self, model, pars=[]):
244        """
245            @return residuals derivatives .
246            @note: in this case just return empty array
247        """
248        return []
249   
250   
251class FitData2D(Data2D):
252    """ Wrapper class  for SANS data """
253    def __init__(self,sans_data2d ,data=None, err_data=None):
254        Data2D.__init__(self, data= data, err_data= err_data)
255        """
256            Data can be initital with a data (sans plottable)
257            or with vectors.
258        """
259        self.x_bins_array = []
260        self.y_bins_array = []
261        self.res_err_image=[]
262        self.index_model=[]
263        self.qmin= None
264        self.qmax= None
265        self.set_data(sans_data2d )
266       
267       
268    def set_data(self, sans_data2d, qmin=None, qmax=None ):
269        """
270            Determine the correct x_bin and y_bin to fit
271        """
272        self.x_bins_array= numpy.reshape(sans_data2d.x_bins,
273                                         [1,len(sans_data2d.x_bins)])
274        self.y_bins_array = numpy.reshape(sans_data2d.y_bins,
275                                          [len(sans_data2d.y_bins),1])
276       
277        x_max = max(sans_data2d.xmin, sans_data2d.xmax)
278        y_max = max(sans_data2d.ymin, sans_data2d.ymax)
279       
280        ## fitting range
281        if qmin == None:
282            self.qmin = 1e-16
283        if qmax == None:
284            self.qmax = math.sqrt(x_max*x_max +y_max*y_max)
285        ## new error image for fitting purpose
286        if self.err_data== None or self.err_data ==[]:
287            self.res_err_data= numpy.zeros(len(self.y_bins),len(self.x_bins))
288        else:
289            self.res_err_data = copy.deepcopy(self.err_data)
290        self.res_err_data[self.res_err_data==0]=1
291       
292        self.radius= numpy.sqrt(self.x_bins_array**2 + self.y_bins_array**2)
293        self.index_model = (self.qmin <= self.radius)&(self.radius<= self.qmax)
294       
295       
296    def setFitRange(self,qmin=None,qmax=None):
297        """ to set the fit range"""
298        if qmin==0.0:
299            self.qmin = 1e-16
300        elif qmin!=None:                       
301            self.qmin = qmin           
302        if qmax!=None:
303            self.qmax= qmax
304           
305        self.radius= numpy.sqrt(self.x_bins_array**2 + self.y_bins_array**2)
306        self.index_model = (self.qmin <= self.radius)&(self.radius<= self.qmax)
307       
308    def getFitRange(self):
309        """
310            @return the range of data.x to fit
311        """
312        return self.qmin, self.qmax
313     
314    def residuals(self, fn): 
315       
316        res=self.index_model*(self.data - fn([self.x_bins_array,
317                             self.y_bins_array]))/self.res_err_data
318        return res.ravel() 
319       
320 
321    def residuals_deriv(self, model, pars=[]):
322        """
323            @return residuals derivatives .
324            @note: in this case just return empty array
325        """
326        return []
327   
328class FitAbort(Exception):
329    """
330        Exception raise to stop the fit
331    """
332    print"Creating fit abort Exception"
333
334
335class SansAssembly:
336    """
337         Sans Assembly class a class wrapper to be call in optimizer.leastsq method
338    """
339    def __init__(self,paramlist,Model=None , Data=None, curr_thread= None):
340        """
341            @param Model: the model wrapper fro sans -model
342            @param Data: the data wrapper for sans data
343        """
344        self.model = Model
345        self.data  = Data
346        self.paramlist=paramlist
347        self.curr_thread= curr_thread
348        self.res=[]
349        self.func_name="Functor"
350    def chisq(self, params):
351        """
352            Calculates chi^2
353            @param params: list of parameter values
354            @return: chi^2
355        """
356        sum = 0
357        for item in self.res:
358            sum += item*item
359        if len(self.res)==0:
360            return None
361        return sum/ len(self.res)
362   
363    def __call__(self,params):
364        """
365            Compute residuals
366            @param params: value of parameters to fit
367        """
368        #import thread
369        self.model.setParams(self.paramlist,params)
370        self.res= self.data.residuals(self.model.eval)
371        #if self.curr_thread != None :
372        #    try:
373        #        self.curr_thread.isquit()
374        #    except:
375        #        raise FitAbort,"stop leastsqr optimizer"   
376        return self.res
377   
378class FitEngine:
379    def __init__(self):
380        """
381            Base class for scipy and park fit engine
382        """
383        #List of parameter names to fit
384        self.paramList=[]
385        #Dictionnary of fitArrange element (fit problems)
386        self.fitArrangeDict={}
387       
388    def _concatenateData(self, listdata=[]):
389        """ 
390            _concatenateData method concatenates each fields of all data contains ins listdata.
391            @param listdata: list of data
392            @return Data: Data is wrapper class for sans plottable. it is created with all parameters
393             of data concatenanted
394            @raise: if listdata is empty  will return None
395            @raise: if data in listdata don't contain dy field ,will create an error
396            during fitting
397        """
398        #TODO: we have to refactor the way we handle data.
399        # We should move away from plottables and move towards the Data1D objects
400        # defined in DataLoader. Data1D allows data manipulations, which should be
401        # used to concatenate.
402        # In the meantime we should switch off the concatenation.
403        #if len(listdata)>1:
404        #    raise RuntimeError, "FitEngine._concatenateData: Multiple data files is not currently supported"
405        #return listdata[0]
406       
407        if listdata==[]:
408            raise ValueError, " data list missing"
409        else:
410            xtemp=[]
411            ytemp=[]
412            dytemp=[]
413            self.mini=None
414            self.maxi=None
415               
416            for item in listdata:
417                data=item.data
418                mini,maxi=data.getFitRange()
419                if self.mini==None and self.maxi==None:
420                    self.mini=mini
421                    self.maxi=maxi
422                else:
423                    if mini < self.mini:
424                        self.mini=mini
425                    if self.maxi < maxi:
426                        self.maxi=maxi
427                       
428                   
429                for i in range(len(data.x)):
430                    xtemp.append(data.x[i])
431                    ytemp.append(data.y[i])
432                    if data.dy is not None and len(data.dy)==len(data.y):   
433                        dytemp.append(data.dy[i])
434                    else:
435                        raise RuntimeError, "Fit._concatenateData: y-errors missing"
436            data= Data(x=xtemp,y=ytemp,dy=dytemp)
437            data.setFitRange(self.mini, self.maxi)
438            return data
439       
440       
441    def set_model(self,model,Uid,pars=[], constraints=[]):
442        """
443            set a model on a given uid in the fit engine.
444            @param model: sans.models type
445            @param Uid :is the key of the fitArrange dictionnary where model is saved as a value
446            @param pars: the list of parameters to fit
447            @param constraints: list of
448                tuple (name of parameter, value of parameters)
449                the value of parameter must be a string to constraint 2 different
450                parameters.
451                Example:
452                we want to fit 2 model M1 and M2 both have parameters A and B.
453                constraints can be:
454                 constraints = [(M1.A, M2.B+2), (M1.B= M2.A *5),...,]
455            @note : pars must contains only name of existing model's paramaters
456        """
457        if model == None:
458            raise ValueError, "AbstractFitEngine: Need to set model to fit"
459       
460        new_model= model
461        if not issubclass(model.__class__, Model):
462            new_model= Model(model)
463       
464        if len(constraints)>0:
465            for constraint in constraints:
466                name, value = constraint
467                try:
468                    new_model.parameterset[ str(name)].set( str(value) )
469                except:
470                    msg= "Fit Engine: Error occurs when setting the constraint"
471                    msg += " %s for parameter %s "%(value, name)
472                    raise ValueError, msg
473               
474        if len(pars) >0:
475            temp=[]
476            for item in pars:
477                if item in new_model.model.getParamList():
478                    temp.append(item)
479                    self.paramList.append(item)
480                else:
481                   
482                    msg = "wrong parameter %s used"%str(item)
483                    msg += "to set model %s. Choose"%str(new_model.model.name)
484                    msg += "parameter name within %s"%str(new_model.model.getParamList())
485                    raise ValueError,msg
486             
487            #A fitArrange is already created but contains dList only at Uid
488            if self.fitArrangeDict.has_key(Uid):
489                self.fitArrangeDict[Uid].set_model(new_model)
490                self.fitArrangeDict[Uid].pars= pars
491            else:
492            #no fitArrange object has been create with this Uid
493                fitproblem = FitArrange()
494                fitproblem.set_model(new_model)
495                fitproblem.pars= pars
496                self.fitArrangeDict[Uid] = fitproblem
497               
498        else:
499            raise ValueError, "park_integration:missing parameters"
500   
501    def set_data(self,data,Uid,smearer=None,qmin=None,qmax=None):
502        """ Receives plottable, creates a list of data to fit,set data
503            in a FitArrange object and adds that object in a dictionary
504            with key Uid.
505            @param data: data added
506            @param Uid: unique key corresponding to a fitArrange object with data
507        """
508        if data.__class__.__name__=='Data2D':
509            fitdata=FitData2D(sans_data2d=data, data=data.data, err_data= data.err_data)
510        else:
511            fitdata=FitData1D(x=data.x, y=data.y , dx= data.dx,dy=data.dy,smearer=smearer)
512       
513        fitdata.setFitRange(qmin=qmin,qmax=qmax)
514        #A fitArrange is already created but contains model only at Uid
515        if self.fitArrangeDict.has_key(Uid):
516            self.fitArrangeDict[Uid].add_data(fitdata)
517        else:
518        #no fitArrange object has been create with this Uid
519            fitproblem= FitArrange()
520            fitproblem.add_data(fitdata)
521            self.fitArrangeDict[Uid]=fitproblem   
522   
523    def get_model(self,Uid):
524        """
525            @param Uid: Uid is key in the dictionary containing the model to return
526            @return  a model at this uid or None if no FitArrange element was created
527            with this Uid
528        """
529        if self.fitArrangeDict.has_key(Uid):
530            return self.fitArrangeDict[Uid].get_model()
531        else:
532            return None
533   
534    def remove_Fit_Problem(self,Uid):
535        """remove   fitarrange in Uid"""
536        if self.fitArrangeDict.has_key(Uid):
537            del self.fitArrangeDict[Uid]
538           
539    def select_problem_for_fit(self,Uid,value):
540        """
541            select a couple of model and data at the Uid position in dictionary
542            and set in self.selected value to value
543            @param value: the value to allow fitting. can only have the value one or zero
544        """
545        if self.fitArrangeDict.has_key(Uid):
546             self.fitArrangeDict[Uid].set_to_fit( value)
547             
548             
549    def get_problem_to_fit(self,Uid):
550        """
551            return the self.selected value of the fit problem of Uid
552           @param Uid: the Uid of the problem
553        """
554        if self.fitArrangeDict.has_key(Uid):
555             self.fitArrangeDict[Uid].get_to_fit()
556   
557class FitArrange:
558    def __init__(self):
559        """
560            Class FitArrange contains a set of data for a given model
561            to perform the Fit.FitArrange must contain exactly one model
562            and at least one data for the fit to be performed.
563            model: the model selected by the user
564            Ldata: a list of data what the user wants to fit
565           
566        """
567        self.model = None
568        self.dList =[]
569        self.pars=[]
570        #self.selected  is zero when this fit problem is not schedule to fit
571        #self.selected is 1 when schedule to fit
572        self.selected = 0
573       
574    def set_model(self,model):
575        """
576            set_model save a copy of the model
577            @param model: the model being set
578        """
579        self.model = model
580       
581    def add_data(self,data):
582        """
583            add_data fill a self.dList with data to fit
584            @param data: Data to add in the list 
585        """
586        if not data in self.dList:
587            self.dList.append(data)
588           
589    def get_model(self):
590        """ @return: saved model """
591        return self.model   
592     
593    def get_data(self):
594        """ @return:  list of data dList"""
595        #return self.dList
596        return self.dList[0] 
597     
598    def remove_data(self,data):
599        """
600            Remove one element from the list
601            @param data: Data to remove from dList
602        """
603        if data in self.dList:
604            self.dList.remove(data)
605    def set_to_fit (self, value=0):
606        """
607           set self.selected to 0 or 1  for other values raise an exception
608           @param value: integer between 0 or 1
609        """
610        self.selected= value
611       
612    def get_to_fit(self):
613        """
614            @return self.selected value
615        """
616        return self.selected
Note: See TracBrowser for help on using the repository browser.