source: trunk/python/asapfitter.py @ 1073

Last change on this file since 1073 was 1061, checked in by mar637, 18 years ago

Added feature Ticket #42; The user can now plot the fits adn reject them.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 17.9 KB
RevLine 
[113]1import _asap
[259]2from asap import rcParams
[723]3from asap import print_log
[113]4
5class fitter:
6    """
7    The fitting class for ASAP.
8    """
[723]9
[113]10    def __init__(self):
11        """
12        Create a fitter object. No state is set.
13        """
14        self.fitter = _asap.fitter()
15        self.x = None
16        self.y = None
17        self.mask = None
18        self.fitfunc = None
[515]19        self.fitfuncs = None
[113]20        self.fitted = False
21        self.data = None
[515]22        self.components = 0
23        self._fittedrow = 0
[113]24        self._p = None
[515]25        self._selection = None
[113]26
27    def set_data(self, xdat, ydat, mask=None):
28        """
[158]29        Set the absissa and ordinate for the fit. Also set the mask
[113]30        indicationg valid points.
31        This can be used for data vectors retrieved from a scantable.
32        For scantable fitting use 'fitter.set_scan(scan, mask)'.
33        Parameters:
[158]34            xdat:    the abcissa values
[113]35            ydat:    the ordinate values
36            mask:    an optional mask
[723]37
[113]38        """
39        self.fitted = False
40        self.x = xdat
41        self.y = ydat
42        if mask == None:
43            from numarray import ones
44            self.mask = ones(len(xdat))
45        else:
46            self.mask = mask
47        return
48
49    def set_scan(self, thescan=None, mask=None):
50        """
51        Set the 'data' (a scantable) of the fitter.
52        Parameters:
53            thescan:     a scantable
54            mask:        a msk retireved from the scantable
55        """
56        if not thescan:
[723]57            msg = "Please give a correct scan"
58            if rcParams['verbose']:
59                print msg
60                return
61            else:
62                raise TypeError(msg)
[113]63        self.fitted = False
64        self.data = thescan
65        if mask is None:
66            from numarray import ones
67            self.mask = ones(self.data.nchan())
68        else:
69            self.mask = mask
70        return
71
72    def set_function(self, **kwargs):
73        """
74        Set the function to be fit.
75        Parameters:
76            poly:    use a polynomial of the order given
77            gauss:   fit the number of gaussian specified
78        Example:
79            fitter.set_function(gauss=2) # will fit two gaussians
80            fitter.set_function(poly=3)  # will fit a 3rd order polynomial
81        """
[723]82        #default poly order 0
[515]83        n=0
[113]84        if kwargs.has_key('poly'):
85            self.fitfunc = 'poly'
86            n = kwargs.get('poly')
[515]87            self.components = [n]
[113]88        elif kwargs.has_key('gauss'):
89            n = kwargs.get('gauss')
90            self.fitfunc = 'gauss'
[515]91            self.fitfuncs = [ 'gauss' for i in range(n) ]
92            self.components = [ 3 for i in range(n) ]
93        else:
[723]94            msg = "Invalid function type."
95            if rcParams['verbose']:
96                print msg
97                return
98            else:
99                raise TypeError(msg)
100
[113]101        self.fitter.setexpression(self.fitfunc,n)
102        return
[723]103
[515]104    def fit(self, row=0):
[113]105        """
106        Execute the actual fitting process. All the state has to be set.
107        Parameters:
[526]108            row:    specify the row in the scantable
[113]109        Example:
[515]110            s = scantable('myscan.asap')
111            s.set_cursor(thepol=1)        # select second pol
[113]112            f = fitter()
113            f.set_scan(s)
114            f.set_function(poly=0)
[723]115            f.fit(row=0)                  # fit first row
[113]116        """
117        if ((self.x is None or self.y is None) and self.data is None) \
118               or self.fitfunc is None:
[723]119            msg = "Fitter not yet initialised. Please set data & fit function"
120            if rcParams['verbose']:
121                print msg
122                return
123            else:
124                raise RuntimeError(msg)
125
[113]126        else:
127            if self.data is not None:
[515]128                self.x = self.data._getabcissa(row)
129                self.y = self.data._getspectrum(row)
[723]130                from asap import asaplog
131                asaplog.push("Fitting:")
[943]132                i = row
133                out = "Scan[%d] Beam[%d] IF[%d] Pol[%d] Cycle[%d]" % (self.data.getscan(i),self.data.getbeam(i),self.data.getif(i),self.data.getpol(i), self.data.getcycle(i))
[876]134                asaplog.push(out)
[515]135        self.fitter.setdata(self.x, self.y, self.mask)
[113]136        if self.fitfunc == 'gauss':
137            ps = self.fitter.getparameters()
138            if len(ps) == 0:
139                self.fitter.estimate()
[626]140        try:
141            self.fitter.fit()
142        except RuntimeError, msg:
[723]143            if rcParams['verbose']:
144                print msg
145            else:
146                raise
[515]147        self._fittedrow = row
[113]148        self.fitted = True
[723]149        print_log()
[113]150        return
151
[515]152    def store_fit(self):
[526]153        """
154        Store the fit parameters in the scantable.
155        """
[515]156        if self.fitted and self.data is not None:
157            pars = list(self.fitter.getparameters())
158            fixed = list(self.fitter.getfixedparameters())
[975]159            from asap.asapfit import asapfit
160            fit = asapfit()
161            fit.setparameters(pars)
162            fit.setfixedparameters(fixed)
163            fit.setfunctions(self.fitfuncs)
164            fit.setcomponents(self.components)
165            fit.setframeinfo(self.data._getcoordinfo())
166            self.data._addfit(fit,self._fittedrow)
[515]167
[1017]168    #def set_parameters(self, params, fixed=None, component=None):
169    def set_parameters(self,*args,**kwargs):
[526]170        """
171        Set the parameters to be fitted.
172        Parameters:
173              params:    a vector of parameters
174              fixed:     a vector of which parameters are to be held fixed
175                         (default is none)
176              component: in case of multiple gaussians, the index of the
177                         component
[1017]178        """
179        component = None
180        fixed = None
181        params = None
[1031]182
[1017]183        if len(args) and isinstance(args[0],dict):
184            kwargs = args[0]
185        if kwargs.has_key("fixed"): fixed = kwargs["fixed"]
186        if kwargs.has_key("params"): params = kwargs["params"]
187        if len(args) == 2 and isinstance(args[1], int):
188            component = args[1]
[515]189        if self.fitfunc is None:
[723]190            msg = "Please specify a fitting function first."
191            if rcParams['verbose']:
192                print msg
193                return
194            else:
195                raise RuntimeError(msg)
[515]196        if self.fitfunc == "gauss" and component is not None:
[1017]197            if not self.fitted and sum(self.fitter.getparameters()) == 0:
[515]198                from numarray import zeros
199                pars = list(zeros(len(self.components)*3))
200                fxd = list(zeros(len(pars)))
201            else:
[723]202                pars = list(self.fitter.getparameters())
[515]203                fxd = list(self.fitter.getfixedparameters())
204            i = 3*component
205            pars[i:i+3] = params
206            fxd[i:i+3] = fixed
207            params = pars
[723]208            fixed = fxd
[113]209        self.fitter.setparameters(params)
210        if fixed is not None:
211            self.fitter.setfixedparameters(fixed)
[723]212        print_log()
[113]213        return
[515]214
215    def set_gauss_parameters(self, peak, centre, fhwm,
[1017]216                             peakfixed=0, centerfixed=0,
217                             fhwmfixed=0,
[515]218                             component=0):
[113]219        """
[515]220        Set the Parameters of a 'Gaussian' component, set with set_function.
221        Parameters:
222            peak, centre, fhwm:  The gaussian parameters
223            peakfixed,
224            centerfixed,
225            fhwmfixed:           Optional parameters to indicate if
226                                 the paramters should be held fixed during
227                                 the fitting process. The default is to keep
228                                 all parameters flexible.
[526]229            component:           The number of the component (Default is the
230                                 component 0)
[515]231        """
232        if self.fitfunc != "gauss":
[723]233            msg = "Function only operates on Gaussian components."
234            if rcParams['verbose']:
235                print msg
236                return
237            else:
238                raise ValueError(msg)
[515]239        if 0 <= component < len(self.components):
[1017]240            d = {'params':[peak, centre, fhwm],
241                 'fixed':[peakfixed, centerfixed, fhwmfixed]}
242            self.set_parameters(d, component)
[515]243        else:
[723]244            msg = "Please select a valid  component."
245            if rcParams['verbose']:
246                print msg
247                return
248            else:
249                raise ValueError(msg)
250
[975]251    def get_area(self, component=None):
252        """
253        Return the area under the fitted gaussian component.
254        Parameters:
255              component:   the gaussian component selection,
256                           default (None) is the sum of all components
257        Note:
258              This will only work for gaussian fits.
259        """
260        if not self.fitted: return
261        if self.fitfunc == "gauss":
262            pars = list(self.fitter.getparameters())
263            from math import log,pi,sqrt
264            fac = sqrt(pi/log(16.0))
265            areas = []
266            for i in range(len(self.components)):
267                j = i*3
268                cpars = pars[j:j+3]
269                areas.append(fac * cpars[0] * cpars[2])
270        else:
271            return None
272        if component is not None:
273            return areas[component]
274        else:
275            return sum(areas)
276
[515]277    def get_parameters(self, component=None):
278        """
[113]279        Return the fit paramters.
[526]280        Parameters:
281             component:    get the parameters for the specified component
282                           only, default is all components
[113]283        """
284        if not self.fitted:
[723]285            msg = "Not yet fitted."
286            if rcParams['verbose']:
287                print msg
288                return
289            else:
290                raise RuntimeError(msg)
[113]291        pars = list(self.fitter.getparameters())
292        fixed = list(self.fitter.getfixedparameters())
[1039]293        area = []
[723]294        if component is not None:
[515]295            if self.fitfunc == "gauss":
296                i = 3*component
297                cpars = pars[i:i+3]
298                cfixed = fixed[i:i+3]
[1039]299                a = self.get_area(component)
300                area = [a for i in range(3)]
[515]301            else:
302                cpars = pars
[723]303                cfixed = fixed
[515]304        else:
305            cpars = pars
306            cfixed = fixed
[1039]307            if self.fitfunc == "gauss":
308                for c in range(len(self.components)):
309                  a = self.get_area(c)
310                  area += [a for i in range(3)]
311
312        fpars = self._format_pars(cpars, cfixed, area)
[723]313        if rcParams['verbose']:
[515]314            print fpars
[1031]315        return {'params':cpars, 'fixed':cfixed, 'formatted': fpars}
[723]316
[975]317    def _format_pars(self, pars, fixed, area):
[113]318        out = ''
319        if self.fitfunc == 'poly':
320            c = 0
[515]321            for i in range(len(pars)):
322                fix = ""
323                if fixed[i]: fix = "(fixed)"
[1061]324                out += '  p%d%s= %3.6f,' % (c,fix,pars[i])
[113]325                c+=1
[515]326            out = out[:-1]  # remove trailing ','
[113]327        elif self.fitfunc == 'gauss':
328            i = 0
329            c = 0
[515]330            aunit = ''
331            ounit = ''
[113]332            if self.data:
[515]333                aunit = self.data.get_unit()
334                ounit = self.data.get_fluxunit()
[113]335            while i < len(pars):
[1039]336                if len(area):
337                    out += '  %2d: peak = %3.3f %s , centre = %3.3f %s, FWHM = %3.3f %s\n      area = %3.3f %s %s\n' % (c,pars[i],ounit,pars[i+1],aunit,pars[i+2],aunit, area[i],ounit,aunit)
[1017]338                else:
339                    out += '  %2d: peak = %3.3f %s , centre = %3.3f %s, FWHM = %3.3f %s\n' % (c,pars[i],ounit,pars[i+1],aunit,pars[i+2],aunit,ounit,aunit)
[113]340                c+=1
341                i+=3
342        return out
[723]343
[113]344    def get_estimate(self):
345        """
[515]346        Return the parameter estimates (for non-linear functions).
[113]347        """
348        pars = self.fitter.getestimate()
[943]349        fixed = self.fitter.getfixedparameters()
[723]350        if rcParams['verbose']:
[1017]351            print self._format_pars(pars,fixed,None)
[113]352        return pars
353
354    def get_residual(self):
355        """
356        Return the residual of the fit.
357        """
358        if not self.fitted:
[723]359            msg = "Not yet fitted."
360            if rcParams['verbose']:
361                print msg
362                return
363            else:
364                raise RuntimeError(msg)
[113]365        return self.fitter.getresidual()
366
367    def get_chi2(self):
368        """
369        Return chi^2.
370        """
371        if not self.fitted:
[723]372            msg = "Not yet fitted."
373            if rcParams['verbose']:
374                print msg
375                return
376            else:
377                raise RuntimeError(msg)
[113]378        ch2 = self.fitter.getchi2()
[723]379        if rcParams['verbose']:
[113]380            print 'Chi^2 = %3.3f' % (ch2)
[723]381        return ch2
[113]382
383    def get_fit(self):
384        """
385        Return the fitted ordinate values.
386        """
387        if not self.fitted:
[723]388            msg = "Not yet fitted."
389            if rcParams['verbose']:
390                print msg
391                return
392            else:
393                raise RuntimeError(msg)
[113]394        return self.fitter.getfit()
395
396    def commit(self):
397        """
[526]398        Return a new scan where the fits have been commited (subtracted)
[113]399        """
400        if not self.fitted:
401            print "Not yet fitted."
[723]402            msg = "Not yet fitted."
403            if rcParams['verbose']:
404                print msg
405                return
406            else:
407                raise RuntimeError(msg)
[975]408        from asap import scantable
409        if not isinstance(self.data, scantable):
[723]410            msg = "Not a scantable"
411            if rcParams['verbose']:
412                print msg
413                return
414            else:
415                raise TypeError(msg)
[113]416        scan = self.data.copy()
[259]417        scan._setspectrum(self.fitter.getresidual())
[723]418        print_log()
[113]419
[723]420    def plot(self, residual=False, components=None, plotparms=False, filename=None):
[113]421        """
422        Plot the last fit.
423        Parameters:
424            residual:    an optional parameter indicating if the residual
425                         should be plotted (default 'False')
[526]426            components:  a list of components to plot, e.g [0,1],
427                         -1 plots the total fit. Default is to only
428                         plot the total fit.
429            plotparms:   Inidicates if the parameter values should be present
430                         on the plot
[113]431        """
432        if not self.fitted:
433            return
[723]434        if not self._p or self._p.is_dead:
435            if rcParams['plotter.gui']:
436                from asap.asaplotgui import asaplotgui as asaplot
437            else:
438                from asap.asaplot import asaplot
439            self._p = asaplot()
440        self._p.hold()
[113]441        self._p.clear()
[515]442        self._p.set_panels()
[652]443        self._p.palette(0)
[113]444        tlab = 'Spectrum'
[723]445        xlab = 'Abcissa'
[1017]446        ylab = 'Ordinate'
447        m = None
[113]448        if self.data:
[515]449            tlab = self.data._getsourcename(self._fittedrow)
450            xlab = self.data._getabcissalabel(self._fittedrow)
451            m = self.data._getmask(self._fittedrow)
[626]452            ylab = self.data._get_ordinate_label()
[515]453
[668]454        colours = ["#777777","#bbbbbb","red","orange","purple","green","magenta", "cyan"]
[652]455        self._p.palette(0,colours)
[515]456        self._p.set_line(label='Spectrum')
[113]457        self._p.plot(self.x, self.y, m)
458        if residual:
[652]459            self._p.palette(1)
[515]460            self._p.set_line(label='Residual')
[113]461            self._p.plot(self.x, self.get_residual(), m)
[652]462        self._p.palette(2)
[515]463        if components is not None:
464            cs = components
465            if isinstance(components,int): cs = [components]
[526]466            if plotparms:
[1031]467                self._p.text(0.15,0.15,str(self.get_parameters()['formatted']),size=8)
[515]468            n = len(self.components)
[652]469            self._p.palette(3)
[515]470            for c in cs:
471                if 0 <= c < n:
472                    lab = self.fitfuncs[c]+str(c)
473                    self._p.set_line(label=lab)
474                    self._p.plot(self.x, self.fitter.evaluate(c), m)
475                elif c == -1:
[652]476                    self._p.palette(2)
[515]477                    self._p.set_line(label="Total Fit")
[723]478                    self._p.plot(self.x, self.get_fit(), m)
[515]479        else:
[652]480            self._p.palette(2)
[515]481            self._p.set_line(label='Fit')
482            self._p.plot(self.x, self.get_fit(), m)
[723]483        xlim=[min(self.x),max(self.x)]
484        self._p.axes.set_xlim(xlim)
[113]485        self._p.set_axes('xlabel',xlab)
486        self._p.set_axes('ylabel',ylab)
487        self._p.set_axes('title',tlab)
488        self._p.release()
[723]489        if (not rcParams['plotter.gui']):
490            self._p.save(filename)
491        print_log()
[113]492
[1061]493    def auto_fit(self, insitu=None, plot=False):
[113]494        """
[515]495        Return a scan where the function is applied to all rows for
496        all Beams/IFs/Pols.
[723]497
[113]498        """
499        from asap import scantable
[515]500        if not isinstance(self.data, scantable) :
[723]501            msg = "Data is not a scantable"
502            if rcParams['verbose']:
503                print msg
504                return
505            else:
506                raise TypeError(msg)
[259]507        if insitu is None: insitu = rcParams['insitu']
508        if not insitu:
509            scan = self.data.copy()
510        else:
511            scan = self.data
[880]512        rows = xrange(scan.nrow())
[723]513        from asap import asaplog
[876]514        asaplog.push("Fitting:")
515        for r in rows:
[1031]516            out = " Scan[%d] Beam[%d] IF[%d] Pol[%d] Cycle[%d]" % (scan.getscan(r),scan.getbeam(r),scan.getif(r),scan.getpol(r), scan.getcycle(r))
[880]517            asaplog.push(out, False)
[876]518            self.x = scan._getabcissa(r)
519            self.y = scan._getspectrum(r)
520            self.data = None
521            self.fit()
522            x = self.get_parameters()
[1061]523            if plot:
524                self.plot(residual=True)
525                x = raw_input("Accept fit ([y]/n): ")
526                if x.upper() == 'N':
527                    continue
[880]528            scan._setspectrum(self.fitter.getresidual(), r)
[1061]529        if plot:
530            self._p.unmap()
531            self._p = None
[876]532        print_log()
533        return scan
[794]534
Note: See TracBrowser for help on using the repository browser.