source: trunk/python/asapfitter.py @ 975

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

Completed Ticket #7 - storing of fits.

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