source: trunk/python/asapfitter.py @ 1826

Last change on this file since 1826 was 1826, checked in by Malte Marquarding, 14 years ago

Tidy up of imports (now imported from asap.). Also fixed some whitespace/tab issues
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 26.1 KB
Line 
1import _asap
2from asap.parameters import rcParams
3from asap.logging import asaplog, print_log, print_log_dec
4from asap.utils import _n_bools, mask_and
5
6
7class fitter:
8    """
9    The fitting class for ASAP.
10    """
11
12    def __init__(self):
13        """
14        Create a fitter object. No state is set.
15        """
16        self.fitter = _asap.fitter()
17        self.x = None
18        self.y = None
19        self.mask = None
20        self.fitfunc = None
21        self.fitfuncs = None
22        self.fitted = False
23        self.data = None
24        self.components = 0
25        self._fittedrow = 0
26        self._p = None
27        self._selection = None
28        self.uselinear = False
29
30    def set_data(self, xdat, ydat, mask=None):
31        """
32        Set the absissa and ordinate for the fit. Also set the mask
33        indicationg valid points.
34        This can be used for data vectors retrieved from a scantable.
35        For scantable fitting use 'fitter.set_scan(scan, mask)'.
36        Parameters:
37            xdat:    the abcissa values
38            ydat:    the ordinate values
39            mask:    an optional mask
40
41        """
42        self.fitted = False
43        self.x = xdat
44        self.y = ydat
45        if mask == None:
46            self.mask = _n_bools(len(xdat), True)
47        else:
48            self.mask = mask
49        return
50
51    def set_scan(self, thescan=None, mask=None):
52        """
53        Set the 'data' (a scantable) of the fitter.
54        Parameters:
55            thescan:     a scantable
56            mask:        a msk retrieved from the scantable
57        """
58        if not thescan:
59            msg = "Please give a correct scan"
60            if rcParams['verbose']:
61                #print msg
62                asaplog.push(msg)
63                print_log('ERROR')
64                return
65            else:
66                raise TypeError(msg)
67        self.fitted = False
68        self.data = thescan
69        self.mask = None
70        if mask is None:
71            self.mask = _n_bools(self.data.nchan(), True)
72        else:
73            self.mask = mask
74        return
75
76    def set_function(self, **kwargs):
77        """
78        Set the function to be fit.
79        Parameters:
80            poly:    use a polynomial of the order given with nonlinear least squares fit
81            lpoly:   use polynomial of the order given with linear least squares fit
82            gauss:   fit the number of gaussian specified
83            lorentz: fit the number of lorentzian specified
84        Example:
85            fitter.set_function(poly=3)  # will fit a 3rd order polynomial via nonlinear method
86            fitter.set_function(lpoly=3)  # will fit a 3rd order polynomial via linear method
87            fitter.set_function(gauss=2) # will fit two gaussians
88            fitter.set_function(lorentz=2) # will fit two lorentzians
89        """
90        #default poly order 0
91        n=0
92        if kwargs.has_key('poly'):
93            self.fitfunc = 'poly'
94            n = kwargs.get('poly')
95            self.components = [n]
96            self.uselinear = False
97        elif kwargs.has_key('lpoly'):
98            self.fitfunc = 'poly'
99            n = kwargs.get('lpoly')
100            self.components = [n]
101            self.uselinear = True
102        elif kwargs.has_key('gauss'):
103            n = kwargs.get('gauss')
104            self.fitfunc = 'gauss'
105            self.fitfuncs = [ 'gauss' for i in range(n) ]
106            self.components = [ 3 for i in range(n) ]
107            self.uselinear = False
108        elif kwargs.has_key('lorentz'):
109            n = kwargs.get('lorentz')
110            self.fitfunc = 'lorentz'
111            self.fitfuncs = [ 'lorentz' for i in range(n) ]
112            self.components = [ 3 for i in range(n) ]
113            self.uselinear = False
114        else:
115            msg = "Invalid function type."
116            if rcParams['verbose']:
117                #print msg
118                asaplog.push(msg)
119                print_log('ERROR')
120                return
121            else:
122                raise TypeError(msg)
123
124        self.fitter.setexpression(self.fitfunc,n)
125        self.fitted = False
126        return
127
128    @print_log_dec
129    def fit(self, row=0, estimate=False):
130        """
131        Execute the actual fitting process. All the state has to be set.
132        Parameters:
133            row:        specify the row in the scantable
134            estimate:   auto-compute an initial parameter set (default False)
135                        This can be used to compute estimates even if fit was
136                        called before.
137        Example:
138            s = scantable('myscan.asap')
139            s.set_cursor(thepol=1)        # select second pol
140            f = fitter()
141            f.set_scan(s)
142            f.set_function(poly=0)
143            f.fit(row=0)                  # fit first row
144        """
145        if ((self.x is None or self.y is None) and self.data is None) \
146               or self.fitfunc is None:
147            msg = "Fitter not yet initialised. Please set data & fit function"
148            if rcParams['verbose']:
149                #print msg
150                asaplog.push(msg)
151                print_log('ERROR')
152                return
153            else:
154                raise RuntimeError(msg)
155
156        else:
157            if self.data is not None:
158                self.x = self.data._getabcissa(row)
159                self.y = self.data._getspectrum(row)
160                self.mask = mask_and(self.mask, self.data._getmask(row))
161                asaplog.push("Fitting:")
162                i = row
163                out = "Scan[%d] Beam[%d] IF[%d] Pol[%d] Cycle[%d]" % (self.data.getscan(i),
164                                                                      self.data.getbeam(i),
165                                                                      self.data.getif(i),
166                                                                      self.data.getpol(i),
167                                                                      self.data.getcycle(i))
168                asaplog.push(out,False)
169        self.fitter.setdata(self.x, self.y, self.mask)
170        if self.fitfunc == 'gauss' or self.fitfunc == 'lorentz':
171            ps = self.fitter.getparameters()
172            if len(ps) == 0 or estimate:
173                self.fitter.estimate()
174        try:
175            fxdpar = list(self.fitter.getfixedparameters())
176            if len(fxdpar) and fxdpar.count(0) == 0:
177                 raise RuntimeError,"No point fitting, if all parameters are fixed."
178            if self.uselinear:
179                converged = self.fitter.lfit()
180            else:
181                converged = self.fitter.fit()
182            if not converged:
183                raise RuntimeError,"Fit didn't converge."
184        except RuntimeError, msg:
185            if rcParams['verbose']:
186                #print msg
187                print_log()
188                asaplog.push(str(msg))
189                print_log('ERROR')
190            else:
191                raise
192        self._fittedrow = row
193        self.fitted = True
194        print_log()
195        return
196
197    def store_fit(self, filename=None):
198        """
199        Save the fit parameters.
200        Parameters:
201            filename:    if specified save as an ASCII file, if None (default)
202                         store it in the scnatable
203        """
204        if self.fitted and self.data is not None:
205            pars = list(self.fitter.getparameters())
206            fixed = list(self.fitter.getfixedparameters())
207            from asap.asapfit import asapfit
208            fit = asapfit()
209            fit.setparameters(pars)
210            fit.setfixedparameters(fixed)
211            fit.setfunctions(self.fitfuncs)
212            fit.setcomponents(self.components)
213            fit.setframeinfo(self.data._getcoordinfo())
214            if filename is not None:
215                import os
216                filename = os.path.expandvars(os.path.expanduser(filename))
217                if os.path.exists(filename):
218                    raise IOError("File '%s' exists." % filename)
219                fit.save(filename)
220            else:
221                self.data._addfit(fit,self._fittedrow)
222
223    @print_log_dec
224    def set_parameters(self,*args,**kwargs):
225        """
226        Set the parameters to be fitted.
227        Parameters:
228              params:    a vector of parameters
229              fixed:     a vector of which parameters are to be held fixed
230                         (default is none)
231              component: in case of multiple gaussians, the index of the
232                         component
233        """
234        component = None
235        fixed = None
236        params = None
237
238        if len(args) and isinstance(args[0],dict):
239            kwargs = args[0]
240        if kwargs.has_key("fixed"): fixed = kwargs["fixed"]
241        if kwargs.has_key("params"): params = kwargs["params"]
242        if len(args) == 2 and isinstance(args[1], int):
243            component = args[1]
244        if self.fitfunc is None:
245            msg = "Please specify a fitting function first."
246            if rcParams['verbose']:
247                #print msg
248                asaplog.push(msg)
249                print_log('ERROR')
250                return
251            else:
252                raise RuntimeError(msg)
253        if (self.fitfunc == "gauss" or self.fitfunc == 'lorentz') and component is not None:
254            if not self.fitted and sum(self.fitter.getparameters()) == 0:
255                pars = _n_bools(len(self.components)*3, False)
256                fxd = _n_bools(len(pars), False)
257            else:
258                pars = list(self.fitter.getparameters())
259                fxd = list(self.fitter.getfixedparameters())
260            i = 3*component
261            pars[i:i+3] = params
262            fxd[i:i+3] = fixed
263            params = pars
264            fixed = fxd
265        self.fitter.setparameters(params)
266        if fixed is not None:
267            self.fitter.setfixedparameters(fixed)
268        print_log()
269        return
270
271    def set_gauss_parameters(self, peak, centre, fwhm,
272                             peakfixed=0, centrefixed=0,
273                             fwhmfixed=0,
274                             component=0):
275        """
276        Set the Parameters of a 'Gaussian' component, set with set_function.
277        Parameters:
278            peak, centre, fwhm:  The gaussian parameters
279            peakfixed,
280            centrefixed,
281            fwhmfixed:           Optional parameters to indicate if
282                                 the paramters should be held fixed during
283                                 the fitting process. The default is to keep
284                                 all parameters flexible.
285            component:           The number of the component (Default is the
286                                 component 0)
287        """
288        if self.fitfunc != "gauss":
289            msg = "Function only operates on Gaussian components."
290            if rcParams['verbose']:
291                #print msg
292                asaplog.push(msg)
293                print_log('ERROR')
294                return
295            else:
296                raise ValueError(msg)
297        if 0 <= component < len(self.components):
298            d = {'params':[peak, centre, fwhm],
299                 'fixed':[peakfixed, centrefixed, fwhmfixed]}
300            self.set_parameters(d, component)
301        else:
302            msg = "Please select a valid  component."
303            if rcParams['verbose']:
304                #print msg
305                asaplog.push(msg)
306                print_log('ERROR')
307                return
308            else:
309                raise ValueError(msg)
310
311    def set_lorentz_parameters(self, peak, centre, fwhm,
312                             peakfixed=0, centrefixed=0,
313                             fwhmfixed=0,
314                             component=0):
315        """
316        Set the Parameters of a 'Lorentzian' component, set with set_function.
317        Parameters:
318            peak, centre, fwhm:  The gaussian parameters
319            peakfixed,
320            centrefixed,
321            fwhmfixed:           Optional parameters to indicate if
322                                 the paramters should be held fixed during
323                                 the fitting process. The default is to keep
324                                 all parameters flexible.
325            component:           The number of the component (Default is the
326                                 component 0)
327        """
328        if self.fitfunc != "lorentz":
329            msg = "Function only operates on Lorentzian components."
330            if rcParams['verbose']:
331                #print msg
332                asaplog.push(msg)
333                print_log('ERROR')
334                return
335            else:
336                raise ValueError(msg)
337        if 0 <= component < len(self.components):
338            d = {'params':[peak, centre, fwhm],
339                 'fixed':[peakfixed, centrefixed, fwhmfixed]}
340            self.set_parameters(d, component)
341        else:
342            msg = "Please select a valid  component."
343            if rcParams['verbose']:
344                #print msg
345                asaplog.push(msg)
346                print_log('ERROR')
347                return
348            else:
349                raise ValueError(msg)
350
351    def get_area(self, component=None):
352        """
353        Return the area under the fitted gaussian/lorentzian component.
354        Parameters:
355              component:   the gaussian/lorentzian component selection,
356                           default (None) is the sum of all components
357        Note:
358              This will only work for gaussian/lorentzian fits.
359        """
360        if not self.fitted: return
361        if self.fitfunc == "gauss" or self.fitfunc == "lorentz":
362            pars = list(self.fitter.getparameters())
363            from math import log,pi,sqrt
364            if self.fitfunc == "gauss":
365                fac = sqrt(pi/log(16.0))
366            elif self.fitfunc == "lorentz":
367                fac = pi/2.0
368            areas = []
369            for i in range(len(self.components)):
370                j = i*3
371                cpars = pars[j:j+3]
372                areas.append(fac * cpars[0] * cpars[2])
373        else:
374            return None
375        if component is not None:
376            return areas[component]
377        else:
378            return sum(areas)
379
380    def get_errors(self, component=None):
381        """
382        Return the errors in the parameters.
383        Parameters:
384            component:    get the errors for the specified component
385                          only, default is all components
386        """
387        if not self.fitted:
388            msg = "Not yet fitted."
389            if rcParams['verbose']:
390                #print msg
391                asaplog.push(msg)
392                print_log('ERROR')
393                return
394            else:
395                raise RuntimeError(msg)
396        errs = list(self.fitter.geterrors())
397        cerrs = errs
398        if component is not None:
399            if self.fitfunc == "gauss" or self.fitfunc == "lorentz":
400                i = 3*component
401                if i < len(errs):
402                    cerrs = errs[i:i+3]
403        return cerrs
404
405    def get_parameters(self, component=None, errors=False):
406        """
407        Return the fit paramters.
408        Parameters:
409             component:    get the parameters for the specified component
410                           only, default is all components
411        """
412        if not self.fitted:
413            msg = "Not yet fitted."
414            if rcParams['verbose']:
415                #print msg
416                asaplog.push(msg)
417                print_log('ERROR')
418                return
419            else:
420                raise RuntimeError(msg)
421        pars = list(self.fitter.getparameters())
422        fixed = list(self.fitter.getfixedparameters())
423        errs = list(self.fitter.geterrors())
424        area = []
425        if component is not None:
426            if self.fitfunc == "gauss" or self.fitfunc == "lorentz":
427                i = 3*component
428                cpars = pars[i:i+3]
429                cfixed = fixed[i:i+3]
430                cerrs = errs[i:i+3]
431                a = self.get_area(component)
432                area = [a for i in range(3)]
433            else:
434                cpars = pars
435                cfixed = fixed
436                cerrs = errs
437        else:
438            cpars = pars
439            cfixed = fixed
440            cerrs = errs
441            if self.fitfunc == "gauss" or self.fitfunc == "lorentz":
442                for c in range(len(self.components)):
443                  a = self.get_area(c)
444                  area += [a for i in range(3)]
445        fpars = self._format_pars(cpars, cfixed, errors and cerrs, area)
446        if rcParams['verbose']:
447            #print fpars
448            asaplog.push(fpars)
449            print_log()
450        return {'params':cpars, 'fixed':cfixed, 'formatted': fpars,
451                'errors':cerrs}
452
453    def _format_pars(self, pars, fixed, errors, area):
454        out = ''
455        if self.fitfunc == 'poly':
456            c = 0
457            for i in range(len(pars)):
458                fix = ""
459                if len(fixed) and fixed[i]: fix = "(fixed)"
460                if errors :
461                    out += '  p%d%s= %3.6f (%1.6f),' % (c,fix,pars[i], errors[i])
462                else:
463                    out += '  p%d%s= %3.6f,' % (c,fix,pars[i])
464                c+=1
465            out = out[:-1]  # remove trailing ','
466        elif self.fitfunc == 'gauss' or self.fitfunc == 'lorentz':
467            i = 0
468            c = 0
469            aunit = ''
470            ounit = ''
471            if self.data:
472                aunit = self.data.get_unit()
473                ounit = self.data.get_fluxunit()
474            while i < len(pars):
475                if len(area):
476                    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)
477                else:
478                    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)
479                c+=1
480                i+=3
481        return out
482
483    def get_estimate(self):
484        """
485        Return the parameter estimates (for non-linear functions).
486        """
487        pars = self.fitter.getestimate()
488        fixed = self.fitter.getfixedparameters()
489        if rcParams['verbose']:
490            #print self._format_pars(pars,fixed,None)
491            asaplog.push(self._format_pars(pars,fixed,None))
492            print_log()
493        return pars
494
495    def get_residual(self):
496        """
497        Return the residual of the fit.
498        """
499        if not self.fitted:
500            msg = "Not yet fitted."
501            if rcParams['verbose']:
502                #print msg
503                asaplog.push(msg)
504                print_log('ERROR')
505                return
506            else:
507                raise RuntimeError(msg)
508        return self.fitter.getresidual()
509
510    def get_chi2(self):
511        """
512        Return chi^2.
513        """
514        if not self.fitted:
515            msg = "Not yet fitted."
516            if rcParams['verbose']:
517                #print msg
518                asaplog.push(msg)
519                print_log('ERROR')
520                return
521            else:
522                raise RuntimeError(msg)
523        ch2 = self.fitter.getchi2()
524        if rcParams['verbose']:
525            #print 'Chi^2 = %3.3f' % (ch2)
526            asaplog.push( 'Chi^2 = %3.3f' % (ch2) )
527            print_log()
528        return ch2
529
530    def get_fit(self):
531        """
532        Return the fitted ordinate values.
533        """
534        if not self.fitted:
535            msg = "Not yet fitted."
536            if rcParams['verbose']:
537                #print msg
538                asaplog.push(msg)
539                print_log('ERROR')
540                return
541            else:
542                raise RuntimeError(msg)
543        return self.fitter.getfit()
544
545    @print_log_dec
546    def commit(self):
547        """
548        Return a new scan where the fits have been commited (subtracted)
549        """
550        if not self.fitted:
551            msg = "Not yet fitted."
552            if rcParams['verbose']:
553                #print msg
554                asaplog.push(msg)
555                print_log('ERROR')
556                return
557            else:
558                raise RuntimeError(msg)
559        from asap import scantable
560        if not isinstance(self.data, scantable):
561            msg = "Not a scantable"
562            if rcParams['verbose']:
563                #print msg
564                asaplog.push(msg)
565                print_log('ERROR')
566                return
567            else:
568                raise TypeError(msg)
569        scan = self.data.copy()
570        scan._setspectrum(self.fitter.getresidual())
571        print_log()
572        return scan
573
574    @print_log_dec
575    def plot(self, residual=False, components=None, plotparms=False,
576             filename=None):
577        """
578        Plot the last fit.
579        Parameters:
580            residual:    an optional parameter indicating if the residual
581                         should be plotted (default 'False')
582            components:  a list of components to plot, e.g [0,1],
583                         -1 plots the total fit. Default is to only
584                         plot the total fit.
585            plotparms:   Inidicates if the parameter values should be present
586                         on the plot
587        """
588        if not self.fitted:
589            return
590        if not self._p or self._p.is_dead:
591            if rcParams['plotter.gui']:
592                from asap.asaplotgui import asaplotgui as asaplot
593            else:
594                from asap.asaplot import asaplot
595            self._p = asaplot()
596        self._p.hold()
597        self._p.clear()
598        self._p.set_panels()
599        self._p.palette(0)
600        tlab = 'Spectrum'
601        xlab = 'Abcissa'
602        ylab = 'Ordinate'
603        from numpy import ma,logical_not,logical_and,array
604        m = self.mask
605        if self.data:
606            tlab = self.data._getsourcename(self._fittedrow)
607            xlab = self.data._getabcissalabel(self._fittedrow)
608            m =  logical_and(self.mask,
609                             array(self.data._getmask(self._fittedrow),
610                                   copy=False))
611
612            ylab = self.data._get_ordinate_label()
613
614        colours = ["#777777","#dddddd","red","orange","purple","green","magenta", "cyan"]
615        nomask=True
616        for i in range(len(m)):
617            nomask = nomask and m[i]
618        label0='Masked Region'
619        label1='Spectrum'
620        if ( nomask ):
621            label0=label1
622        else:
623            y = ma.masked_array( self.y, mask = m )
624            self._p.palette(1,colours)
625            self._p.set_line( label = label1 )
626            self._p.plot( self.x, y )
627        self._p.palette(0,colours)
628        self._p.set_line(label=label0)
629        y = ma.masked_array(self.y,mask=logical_not(m))
630        self._p.plot(self.x, y)
631        if residual:
632            self._p.palette(7)
633            self._p.set_line(label='Residual')
634            y = ma.masked_array(self.get_residual(),
635                                  mask=logical_not(m))
636            self._p.plot(self.x, y)
637        self._p.palette(2)
638        if components is not None:
639            cs = components
640            if isinstance(components,int): cs = [components]
641            if plotparms:
642                self._p.text(0.15,0.15,str(self.get_parameters()['formatted']),size=8)
643            n = len(self.components)
644            self._p.palette(3)
645            for c in cs:
646                if 0 <= c < n:
647                    lab = self.fitfuncs[c]+str(c)
648                    self._p.set_line(label=lab)
649                    y = ma.masked_array(self.fitter.evaluate(c),
650                                          mask=logical_not(m))
651
652                    self._p.plot(self.x, y)
653                elif c == -1:
654                    self._p.palette(2)
655                    self._p.set_line(label="Total Fit")
656                    y = ma.masked_array(self.fitter.getfit(),
657                                          mask=logical_not(m))
658                    self._p.plot(self.x, y)
659        else:
660            self._p.palette(2)
661            self._p.set_line(label='Fit')
662            y = ma.masked_array(self.fitter.getfit(),
663                                  mask=logical_not(m))
664            self._p.plot(self.x, y)
665        xlim=[min(self.x),max(self.x)]
666        self._p.axes.set_xlim(xlim)
667        self._p.set_axes('xlabel',xlab)
668        self._p.set_axes('ylabel',ylab)
669        self._p.set_axes('title',tlab)
670        self._p.release()
671        if (not rcParams['plotter.gui']):
672            self._p.save(filename)
673        print_log()
674
675    @print_log_dec
676    def auto_fit(self, insitu=None, plot=False):
677        """
678        Return a scan where the function is applied to all rows for
679        all Beams/IFs/Pols.
680
681        """
682        from asap import scantable
683        if not isinstance(self.data, scantable) :
684            msg = "Data is not a scantable"
685            if rcParams['verbose']:
686                #print msg
687                asaplog.push(msg)
688                print_log('ERROR')
689                return
690            else:
691                raise TypeError(msg)
692        if insitu is None: insitu = rcParams['insitu']
693        if not insitu:
694            scan = self.data.copy()
695        else:
696            scan = self.data
697        rows = xrange(scan.nrow())
698        # Save parameters of baseline fits as a class attribute.
699        # NOTICE: This does not reflect changes in scantable!
700        if len(rows) > 0: self.blpars=[]
701        asaplog.push("Fitting:")
702        for r in rows:
703            out = " Scan[%d] Beam[%d] IF[%d] Pol[%d] Cycle[%d]" % (scan.getscan(r),
704                                                                   scan.getbeam(r),
705                                                                   scan.getif(r),
706                                                                   scan.getpol(r),
707                                                                   scan.getcycle(r))
708            asaplog.push(out, False)
709            self.x = scan._getabcissa(r)
710            self.y = scan._getspectrum(r)
711            self.mask = mask_and(self.mask, scan._getmask(r))
712            self.data = None
713            self.fit()
714            x = self.get_parameters()
715            fpar = self.get_parameters()
716            if plot:
717                self.plot(residual=True)
718                x = raw_input("Accept fit ([y]/n): ")
719                if x.upper() == 'N':
720                    self.blpars.append(None)
721                    continue
722            scan._setspectrum(self.fitter.getresidual(), r)
723            self.blpars.append(fpar)
724        if plot:
725            self._p.unmap()
726            self._p = None
727        print_log()
728        return scan
Note: See TracBrowser for help on using the repository browser.