import _asap class fitter: """ The fitting class for ASAP. """ def _verbose(self, *args): """ Set stdout output. """ if type(args[0]) is bool: self._vb = args[0] return elif len(args) == 0: return self._vb def __init__(self): """ Create a fitter object. No state is set. """ self.fitter = _asap.fitter() self.x = None self.y = None self.mask = None self.fitfunc = None self.fitted = False self.data = None self._p = None self._vb = True def set_data(self, xdat, ydat, mask=None): """ Set the abscissa and ordinate for the fit. Also set the mask indicationg valid points. This can be used for data vectors retrieved from a scantable. For scantable fitting use 'fitter.set_scan(scan, mask)'. Parameters: xdat: the abscissa values ydat: the ordinate values mask: an optional mask """ self.fitted = False self.x = xdat self.y = ydat if mask == None: from numarray import ones self.mask = ones(len(xdat)) else: self.mask = mask return def set_scan(self, thescan=None, mask=None): """ Set the 'data' (a scantable) of the fitter. Parameters: thescan: a scantable mask: a msk retireved from the scantable """ if not thescan: print "Please give a correct scan" self.fitted = False self.data = thescan if mask is None: from numarray import ones self.mask = ones(self.data.nchan()) else: self.mask = mask return def set_function(self, **kwargs): """ Set the function to be fit. Parameters: poly: use a polynomial of the order given gauss: fit the number of gaussian specified Example: fitter.set_function(gauss=2) # will fit two gaussians fitter.set_function(poly=3) # will fit a 3rd order polynomial """ #default poly order 0 self.fitfunc = 'poly' n=0 if kwargs.has_key('poly'): self.fitfunc = 'poly' n = kwargs.get('poly') elif kwargs.has_key('gauss'): n = kwargs.get('gauss') self.fitfunc = 'gauss' self.fitter.setexpression(self.fitfunc,n) return def fit(self): """ Execute the actual fitting process. All the state has to be set. Parameters: none Example: s= scantable('myscan.asap') f = fitter() f.set_scan(s) f.set_function(poly=0) f.fit() """ if ((self.x is None or self.y is None) and self.data is None) \ or self.fitfunc is None: print "Fitter not yet initialised. Please set data & fit function" return else: if self.data is not None: self.x = self.data.getabscissa() self.y = self.data.getspectrum() print "Fitting:" vb = self.data._verbose self.data._verbose(True) s = self.data.get_selection() self.data._verbose(vb) self.fitter.setdata(self.x,self.y,self.mask) if self.fitfunc == 'gauss': ps = self.fitter.getparameters() if len(ps) == 0: self.fitter.estimate() self.fitter.fit() self.fitted = True return def set_parameters(self, params, fixed=None): self.fitter.setparameters(params) if fixed is not None: self.fitter.setfixedparameters(fixed) return def get_parameters(self): """ Return the fit paramters. """ if not self.fitted: print "Not yet fitted." pars = list(self.fitter.getparameters()) fixed = list(self.fitter.getfixedparameters()) if self._vb: print self._format_pars(pars) return pars,fixed def _format_pars(self, pars): out = '' if self.fitfunc == 'poly': c = 0 for i in pars: out += ' p%d = %3.3f, ' % (c,i) c+=1 elif self.fitfunc == 'gauss': i = 0 c = 0 unit = '' if self.data: unit = self.data.get_unit() while i < len(pars): out += ' %d: peak = %3.3f , centre = %3.3f %s, FWHM = %3.3f %s \n' % (c,pars[i],pars[i+1],unit,pars[i+2],unit) c+=1 i+=3 return out def get_estimate(self): """ Return the paramter estimates (for non-linear functions). """ pars = self.fitter.getestimate() if self._vb: print self._format_pars(pars) return pars def get_residual(self): """ Return the residual of the fit. """ if not self.fitted: print "Not yet fitted." return self.fitter.getresidual() def get_chi2(self): """ Return chi^2. """ if not self.fitted: print "Not yet fitted." ch2 = self.fitter.getchi2() if self._vb: print 'Chi^2 = %3.3f' % (ch2) return ch2 def get_fit(self): """ Return the fitted ordinate values. """ if not self.fitted: print "Not yet fitted." return self.fitter.getfit() def commit(self): """ Return a new scan where teh fits have been commited. """ if not self.fitted: print "Not yet fitted." if self.data is not scantable: print "Only works with scantables" return scan = self.data.copy() scan.setspectrum(self.fitter.getresidual()) def plot(self, residual=False): """ Plot the last fit. Parameters: residual: an optional parameter indicating if the residual should be plotted (default 'False') """ if not self.fitted: return if not self._p: from asap.asaplot import ASAPlot self._p = ASAPlot() self._p.clear() tlab = 'Spectrum' xlab = 'Abscissa' if self.data: tlab = self.data._getsourcename(0) xlab = self.data.getabscissalabel(0) ylab = r'Flux' m = self.data.getmask(0) self._p.set_line(colour='blue',label='Spectrum') self._p.plot(self.x, self.y, m) if residual: self._p.set_line(colour='green',label='Residual') self._p.plot(self.x, self.get_residual(), m) self._p.set_line(colour='red',label='Fit') self._p.plot(self.x, self.get_fit(), m) self._p.set_axes('xlabel',xlab) self._p.set_axes('ylabel',ylab) self._p.set_axes('title',tlab) self._p.release() def auto_fit(self): """ Return a scan where the function is applied to all Beams/IFs/Pols. """ from asap import scantable if not isinstance(self.data,scantable) : print "Only works with scantables" return scan = self.data.copy() vb = scan._verbose scan._verbose(False) sel = scan.get_selection() for i in range(scan.nbeam()): scan.setbeam(i) for j in range(scan.nif()): scan.setif(j) for k in range(scan.npol()): scan.setpol(k) self.x = scan.getabscissa() self.y = scan.getspectrum() self.data = None self.fit() if self._vb: print "Fitting:" print 'Beam[%d], IF[%d], Pol[%d]' % (i,j,k) x = self.get_parameters() scan.setspectrum(self.fitter.getresidual()) scan.set_selection(sel[0],sel[1],sel[2]) scan._verbose(vb) return scan