|
casa
$Rev:20696$
|
|
casa
$Rev:20696$
|
00001 import _asap 00002 from asap.parameters import rcParams 00003 from asap.logging import asaplog, asaplog_post_dec 00004 from asap.utils import _n_bools, mask_and 00005 from numpy import ndarray 00006 00007 class fitter: 00008 """ 00009 The fitting class for ASAP. 00010 """ 00011 def __init__(self): 00012 """ 00013 Create a fitter object. No state is set. 00014 """ 00015 self.fitter = _asap.fitter() 00016 self.x = None 00017 self.y = None 00018 self.mask = None 00019 self.fitfunc = None 00020 self.fitfuncs = None 00021 self.fitted = False 00022 self.data = None 00023 self.components = 0 00024 self._fittedrow = 0 00025 self._p = None 00026 self._selection = None 00027 self.uselinear = False 00028 self._constraints = [] 00029 00030 def set_data(self, xdat, ydat, mask=None): 00031 """ 00032 Set the absissa and ordinate for the fit. Also set the mask 00033 indicating valid points. 00034 This can be used for data vectors retrieved from a scantable. 00035 For scantable fitting use 'fitter.set_scan(scan, mask)'. 00036 Parameters: 00037 xdat: the abcissa values 00038 ydat: the ordinate values 00039 mask: an optional mask 00040 00041 """ 00042 self.fitted = False 00043 self.x = xdat 00044 self.y = ydat 00045 if mask == None: 00046 self.mask = _n_bools(len(xdat), True) 00047 else: 00048 self.mask = mask 00049 return 00050 00051 @asaplog_post_dec 00052 def set_scan(self, thescan=None, mask=None): 00053 """ 00054 Set the 'data' (a scantable) of the fitter. 00055 Parameters: 00056 thescan: a scantable 00057 mask: a msk retrieved from the scantable 00058 """ 00059 if not thescan: 00060 msg = "Please give a correct scan" 00061 raise TypeError(msg) 00062 self.fitted = False 00063 self.data = thescan 00064 self.mask = None 00065 if mask is None: 00066 self.mask = _n_bools(self.data.nchan(), True) 00067 else: 00068 self.mask = mask 00069 return 00070 00071 @asaplog_post_dec 00072 def set_function(self, **kwargs): 00073 """ 00074 Set the function to be fit. 00075 Parameters: 00076 poly: use a polynomial of the order given with nonlinear 00077 least squares fit 00078 lpoly: use polynomial of the order given with linear least 00079 squares fit 00080 gauss: fit the number of gaussian specified 00081 lorentz: fit the number of lorentzian specified 00082 sinusoid: fit the number of sinusoid specified 00083 Example: 00084 fitter.set_function(poly=3) # will fit a 3rd order polynomial 00085 # via nonlinear method 00086 fitter.set_function(lpoly=3) # will fit a 3rd order polynomial 00087 # via linear method 00088 fitter.set_function(gauss=2) # will fit two gaussians 00089 fitter.set_function(lorentz=2) # will fit two lorentzians 00090 fitter.set_function(sinusoid=3) # will fit three sinusoids 00091 """ 00092 #default poly order 0 00093 n=0 00094 if kwargs.has_key('poly'): 00095 self.fitfunc = 'poly' 00096 self.fitfuncs = ['poly'] 00097 n = kwargs.get('poly') 00098 self.components = [n+1] 00099 self.uselinear = False 00100 elif kwargs.has_key('lpoly'): 00101 self.fitfunc = 'poly' 00102 self.fitfuncs = ['lpoly'] 00103 n = kwargs.get('lpoly') 00104 self.components = [n+1] 00105 self.uselinear = True 00106 elif kwargs.has_key('gauss'): 00107 n = kwargs.get('gauss') 00108 self.fitfunc = 'gauss' 00109 self.fitfuncs = [ 'gauss' for i in range(n) ] 00110 self.components = [ 3 for i in range(n) ] 00111 self.uselinear = False 00112 elif kwargs.has_key('lorentz'): 00113 n = kwargs.get('lorentz') 00114 self.fitfunc = 'lorentz' 00115 self.fitfuncs = [ 'lorentz' for i in range(n) ] 00116 self.components = [ 3 for i in range(n) ] 00117 self.uselinear = False 00118 elif kwargs.has_key('sinusoid'): 00119 n = kwargs.get('sinusoid') 00120 self.fitfunc = 'sinusoid' 00121 self.fitfuncs = [ 'sinusoid' for i in range(n) ] 00122 self.components = [ 3 for i in range(n) ] 00123 self.uselinear = False 00124 elif kwargs.has_key('expression'): 00125 self.uselinear = False 00126 raise RuntimeError("Not yet implemented") 00127 else: 00128 msg = "Invalid function type." 00129 raise TypeError(msg) 00130 00131 self.fitter.setexpression(self.fitfunc,n) 00132 self._constraints = [] 00133 self.fitted = False 00134 return 00135 00136 @asaplog_post_dec 00137 def fit(self, row=0, estimate=False): 00138 """ 00139 Execute the actual fitting process. All the state has to be set. 00140 Parameters: 00141 row: specify the row in the scantable 00142 estimate: auto-compute an initial parameter set (default False) 00143 This can be used to compute estimates even if fit was 00144 called before. 00145 Example: 00146 s = scantable('myscan.asap') 00147 s.set_cursor(thepol=1) # select second pol 00148 f = fitter() 00149 f.set_scan(s) 00150 f.set_function(poly=0) 00151 f.fit(row=0) # fit first row 00152 """ 00153 if ((self.x is None or self.y is None) and self.data is None) \ 00154 or self.fitfunc is None: 00155 msg = "Fitter not yet initialised. Please set data & fit function" 00156 raise RuntimeError(msg) 00157 00158 if self.data is not None: 00159 self.x = self.data._getabcissa(row) 00160 self.y = self.data._getspectrum(row) 00161 #self.mask = mask_and(self.mask, self.data._getmask(row)) 00162 if len(self.x) == len(self.mask): 00163 self.mask = mask_and(self.mask, self.data._getmask(row)) 00164 else: 00165 asaplog.push('lengths of data and mask are not the same. ' 00166 'preset mask will be ignored') 00167 asaplog.post('WARN','asapfit.fit') 00168 self.mask=self.data._getmask(row) 00169 asaplog.push("Fitting:") 00170 i = row 00171 out = "Scan[%d] Beam[%d] IF[%d] Pol[%d] Cycle[%d]" % ( 00172 self.data.getscan(i), 00173 self.data.getbeam(i), 00174 self.data.getif(i), 00175 self.data.getpol(i), 00176 self.data.getcycle(i)) 00177 00178 asaplog.push(out, False) 00179 00180 self.fitter.setdata(self.x, self.y, self.mask) 00181 if self.fitfunc == 'gauss' or self.fitfunc == 'lorentz': 00182 ps = self.fitter.getparameters() 00183 if len(ps) == 0 or estimate: 00184 self.fitter.estimate() 00185 fxdpar = list(self.fitter.getfixedparameters()) 00186 if len(fxdpar) and fxdpar.count(0) == 0: 00187 raise RuntimeError,"No point fitting, if all parameters are fixed." 00188 if self._constraints: 00189 for c in self._constraints: 00190 self.fitter.addconstraint(c[0]+[c[-1]]) 00191 if self.uselinear: 00192 converged = self.fitter.lfit() 00193 else: 00194 converged = self.fitter.fit() 00195 if not converged: 00196 raise RuntimeError,"Fit didn't converge." 00197 self._fittedrow = row 00198 self.fitted = True 00199 return 00200 00201 def store_fit(self, filename=None): 00202 """ 00203 Save the fit parameters. 00204 Parameters: 00205 filename: if specified save as an ASCII file, if None (default) 00206 store it in the scnatable 00207 """ 00208 if self.fitted and self.data is not None: 00209 pars = list(self.fitter.getparameters()) 00210 fixed = list(self.fitter.getfixedparameters()) 00211 from asap.asapfit import asapfit 00212 fit = asapfit() 00213 fit.setparameters(pars) 00214 fit.setfixedparameters(fixed) 00215 fit.setfunctions(self.fitfuncs) 00216 fit.setcomponents(self.components) 00217 fit.setframeinfo(self.data._getcoordinfo()) 00218 if filename is not None: 00219 import os 00220 filename = os.path.expandvars(os.path.expanduser(filename)) 00221 if os.path.exists(filename): 00222 raise IOError("File '%s' exists." % filename) 00223 fit.save(filename) 00224 else: 00225 self.data._addfit(fit,self._fittedrow) 00226 00227 @asaplog_post_dec 00228 def set_parameters(self,*args,**kwargs): 00229 """ 00230 Set the parameters to be fitted. 00231 Parameters: 00232 params: a vector of parameters 00233 fixed: a vector of which parameters are to be held fixed 00234 (default is none) 00235 component: in case of multiple gaussians/lorentzians/sinusoidals, 00236 the index of the target component 00237 """ 00238 component = None 00239 fixed = None 00240 params = None 00241 00242 if len(args) and isinstance(args[0],dict): 00243 kwargs = args[0] 00244 if kwargs.has_key("fixed"): fixed = kwargs["fixed"] 00245 if kwargs.has_key("params"): params = kwargs["params"] 00246 if len(args) == 2 and isinstance(args[1], int): 00247 component = args[1] 00248 if self.fitfunc is None: 00249 msg = "Please specify a fitting function first." 00250 raise RuntimeError(msg) 00251 if (self.fitfunc == "gauss" or self.fitfunc == "lorentz" 00252 or self.fitfunc == "sinusoid") and component is not None: 00253 if not self.fitted and sum(self.fitter.getparameters()) == 0: 00254 pars = _n_bools(len(self.components)*3, False) 00255 fxd = _n_bools(len(pars), False) 00256 else: 00257 pars = list(self.fitter.getparameters()) 00258 fxd = list(self.fitter.getfixedparameters()) 00259 i = 3*component 00260 pars[i:i+3] = params 00261 fxd[i:i+3] = fixed 00262 params = pars 00263 fixed = fxd 00264 self.fitter.setparameters(params) 00265 if fixed is not None: 00266 self.fitter.setfixedparameters(fixed) 00267 return 00268 00269 @asaplog_post_dec 00270 def set_gauss_parameters(self, peak, centre, fwhm, 00271 peakfixed=0, centrefixed=0, 00272 fwhmfixed=0, 00273 component=0): 00274 """ 00275 Set the Parameters of a 'Gaussian' component, set with set_function. 00276 Parameters: 00277 peak, centre, fwhm: The gaussian parameters 00278 peakfixed, 00279 centrefixed, 00280 fwhmfixed: Optional parameters to indicate if 00281 the paramters should be held fixed during 00282 the fitting process. The default is to keep 00283 all parameters flexible. 00284 component: The number of the component (Default is the 00285 component 0) 00286 """ 00287 if self.fitfunc != "gauss": 00288 msg = "Function only operates on Gaussian components." 00289 raise ValueError(msg) 00290 if 0 <= component < len(self.components): 00291 d = {'params':[peak, centre, fwhm], 00292 'fixed':[peakfixed, centrefixed, fwhmfixed]} 00293 self.set_parameters(d, component) 00294 else: 00295 msg = "Please select a valid component." 00296 raise ValueError(msg) 00297 00298 @asaplog_post_dec 00299 def set_lorentz_parameters(self, peak, centre, fwhm, 00300 peakfixed=0, centrefixed=0, 00301 fwhmfixed=0, 00302 component=0): 00303 """ 00304 Set the Parameters of a 'Lorentzian' component, set with set_function. 00305 Parameters: 00306 peak, centre, fwhm: The lorentzian parameters 00307 peakfixed, 00308 centrefixed, 00309 fwhmfixed: Optional parameters to indicate if 00310 the paramters should be held fixed during 00311 the fitting process. The default is to keep 00312 all parameters flexible. 00313 component: The number of the component (Default is the 00314 component 0) 00315 """ 00316 if self.fitfunc != "lorentz": 00317 msg = "Function only operates on Lorentzian components." 00318 raise ValueError(msg) 00319 if 0 <= component < len(self.components): 00320 d = {'params':[peak, centre, fwhm], 00321 'fixed':[peakfixed, centrefixed, fwhmfixed]} 00322 self.set_parameters(d, component) 00323 else: 00324 msg = "Please select a valid component." 00325 raise ValueError(msg) 00326 00327 @asaplog_post_dec 00328 def set_sinusoid_parameters(self, ampl, period, x0, 00329 amplfixed=0, periodfixed=0, 00330 x0fixed=0, 00331 component=0): 00332 """ 00333 Set the Parameters of a 'Sinusoidal' component, set with set_function. 00334 Parameters: 00335 ampl, period, x0: The sinusoidal parameters 00336 amplfixed, 00337 periodfixed, 00338 x0fixed: Optional parameters to indicate if 00339 the paramters should be held fixed during 00340 the fitting process. The default is to keep 00341 all parameters flexible. 00342 component: The number of the component (Default is the 00343 component 0) 00344 """ 00345 if self.fitfunc != "sinusoid": 00346 msg = "Function only operates on Sinusoidal components." 00347 raise ValueError(msg) 00348 if 0 <= component < len(self.components): 00349 d = {'params':[ampl, period, x0], 00350 'fixed': [amplfixed, periodfixed, x0fixed]} 00351 self.set_parameters(d, component) 00352 else: 00353 msg = "Please select a valid component." 00354 raise ValueError(msg) 00355 00356 00357 def add_constraint(self, xpar, y): 00358 """Add parameter constraints to the fit. This is done by setting up 00359 linear equations for the related parameters. 00360 00361 For example a two component gaussian fit where the amplitudes are 00362 constraint by amp1 = 2*amp2 and paramaters for the two components 00363 in the order [amp1,peakv1,sigma1,amp2,peakv2,sigma2] 00364 needs a constraint 00365 00366 add_constraint([1, 0, 0, -2, 0, 0], 0) 00367 00368 as the linear equation is 00369 00370 1*amp1 + 0*peakv1 + 0*sigma1 -2*amp2 + 0*peakv2 + 0*sigma2 = 0 00371 00372 and similarly for a velocity difference of v2-v1=17 00373 00374 add_constraint([0.,-1.,0.,0.,1.,0.], 17.) 00375 00376 """ 00377 self._constraints.append((xpar, y)) 00378 00379 00380 def get_area(self, component=None): 00381 """ 00382 Return the area under the fitted gaussian/lorentzian component. 00383 Parameters: 00384 component: the gaussian/lorentzian component selection, 00385 default (None) is the sum of all components 00386 Note: 00387 This will only work for gaussian/lorentzian fits. 00388 """ 00389 if not self.fitted: return 00390 if self.fitfunc == "gauss" or self.fitfunc == "lorentz": 00391 pars = list(self.fitter.getparameters()) 00392 from math import log,pi,sqrt 00393 if self.fitfunc == "gauss": 00394 fac = sqrt(pi/log(16.0)) 00395 elif self.fitfunc == "lorentz": 00396 fac = pi/2.0 00397 areas = [] 00398 for i in range(len(self.components)): 00399 j = i*3 00400 cpars = pars[j:j+3] 00401 areas.append(fac * cpars[0] * cpars[2]) 00402 else: 00403 return None 00404 if component is not None: 00405 return areas[component] 00406 else: 00407 return sum(areas) 00408 00409 @asaplog_post_dec 00410 def get_errors(self, component=None): 00411 """ 00412 Return the errors in the parameters. 00413 Parameters: 00414 component: get the errors for the specified component 00415 only, default is all components 00416 """ 00417 if not self.fitted: 00418 msg = "Not yet fitted." 00419 raise RuntimeError(msg) 00420 errs = list(self.fitter.geterrors()) 00421 cerrs = errs 00422 if component is not None: 00423 if self.fitfunc == "gauss" or self.fitfunc == "lorentz" \ 00424 or self.fitfunc == "sinusoid": 00425 i = 3*component 00426 if i < len(errs): 00427 cerrs = errs[i:i+3] 00428 return cerrs 00429 00430 00431 @asaplog_post_dec 00432 def get_parameters(self, component=None, errors=False): 00433 """ 00434 Return the fit paramters. 00435 Parameters: 00436 component: get the parameters for the specified component 00437 only, default is all components 00438 """ 00439 if not self.fitted: 00440 msg = "Not yet fitted." 00441 raise RuntimeError(msg) 00442 pars = list(self.fitter.getparameters()) 00443 fixed = list(self.fitter.getfixedparameters()) 00444 errs = list(self.fitter.geterrors()) 00445 area = [] 00446 if component is not None: 00447 if self.fitfunc == "poly" or self.fitfunc == "lpoly": 00448 cpars = pars 00449 cfixed = fixed 00450 cerrs = errs 00451 else: 00452 i = 3*component 00453 cpars = pars[i:i+3] 00454 cfixed = fixed[i:i+3] 00455 cerrs = errs[i:i+3] 00456 if self.fitfunc == "gauss" or self.fitfunc == "lorentz": 00457 a = self.get_area(component) 00458 area = [a for i in range(3)] 00459 else: 00460 cpars = pars 00461 cfixed = fixed 00462 cerrs = errs 00463 if self.fitfunc == "gauss" or self.fitfunc == "lorentz": 00464 for c in range(len(self.components)): 00465 a = self.get_area(c) 00466 area += [a for i in range(3)] 00467 fpars = self._format_pars(cpars, cfixed, errors and cerrs, area) 00468 asaplog.push(fpars) 00469 return {'params':cpars, 'fixed':cfixed, 'formatted': fpars, 00470 'errors':cerrs} 00471 00472 def _format_pars(self, pars, fixed, errors, area): 00473 out = '' 00474 if self.fitfunc == "poly" or self.fitfunc == "lpoly": 00475 c = 0 00476 for i in range(len(pars)): 00477 fix = "" 00478 if len(fixed) and fixed[i]: fix = "(fixed)" 00479 out += " p%d%s= %3.6f" % (c, fix, pars[i]) 00480 if errors : out += " (%1.6f)" % errors[i] 00481 out += "," 00482 c+=1 00483 out = out[:-1] # remove trailing ',' 00484 else: 00485 i = 0 00486 c = 0 00487 if self.fitfunc == "gauss" or self.fitfunc == "lorentz": 00488 pnam = ["peak", "centre", "FWHM"] 00489 elif self.fitfunc == "sinusoid": 00490 pnam = ["amplitude", "period", "x0"] 00491 aunit = "" 00492 ounit = "" 00493 if self.data: 00494 aunit = self.data.get_unit() 00495 ounit = self.data.get_fluxunit() 00496 while i < len(pars): 00497 fix0 = fix1 = fix2 = "" 00498 if i < len(fixed)-2: 00499 if fixed[i]: fix0 = "(fixed)" 00500 if fixed[i+1]: fix1 = "(fixed)" 00501 if fixed[i+2]: fix2 = "(fixed)" 00502 out += " %2d: " % c 00503 out += "%s%s = %3.3f %s, " % (pnam[0], fix0, pars[i], ounit) 00504 out += "%s%s = %3.3f %s, " % (pnam[1], fix1, pars[i+1], aunit) 00505 out += "%s%s = %3.3f %s\n" % (pnam[2], fix2, pars[i+2], aunit) 00506 if len(area): out += " area = %3.3f %s %s\n" % (area[i], 00507 ounit, 00508 aunit) 00509 c+=1 00510 i+=3 00511 return out 00512 00513 00514 @asaplog_post_dec 00515 def get_estimate(self): 00516 """ 00517 Return the parameter estimates (for non-linear functions). 00518 """ 00519 pars = self.fitter.getestimate() 00520 fixed = self.fitter.getfixedparameters() 00521 asaplog.push(self._format_pars(pars,fixed,None,None)) 00522 return pars 00523 00524 @asaplog_post_dec 00525 def get_residual(self): 00526 """ 00527 Return the residual of the fit. 00528 """ 00529 if not self.fitted: 00530 msg = "Not yet fitted." 00531 raise RuntimeError(msg) 00532 return self.fitter.getresidual() 00533 00534 @asaplog_post_dec 00535 def get_chi2(self): 00536 """ 00537 Return chi^2. 00538 """ 00539 if not self.fitted: 00540 msg = "Not yet fitted." 00541 raise RuntimeError(msg) 00542 ch2 = self.fitter.getchi2() 00543 asaplog.push( 'Chi^2 = %3.3f' % (ch2) ) 00544 return ch2 00545 00546 @asaplog_post_dec 00547 def get_fit(self): 00548 """ 00549 Return the fitted ordinate values. 00550 """ 00551 if not self.fitted: 00552 msg = "Not yet fitted." 00553 raise RuntimeError(msg) 00554 return self.fitter.getfit() 00555 00556 @asaplog_post_dec 00557 def commit(self): 00558 """ 00559 Return a new scan where the fits have been commited (subtracted) 00560 """ 00561 if not self.fitted: 00562 msg = "Not yet fitted." 00563 raise RuntimeError(msg) 00564 from asap import scantable 00565 if not isinstance(self.data, scantable): 00566 msg = "Not a scantable" 00567 raise TypeError(msg) 00568 scan = self.data.copy() 00569 scan._setspectrum(self.fitter.getresidual()) 00570 return scan 00571 00572 @asaplog_post_dec 00573 def plot(self, residual=False, components=None, plotparms=False, 00574 filename=None): 00575 """ 00576 Plot the last fit. 00577 Parameters: 00578 residual: an optional parameter indicating if the residual 00579 should be plotted (default 'False') 00580 components: a list of components to plot, e.g [0,1], 00581 -1 plots the total fit. Default is to only 00582 plot the total fit. 00583 plotparms: Inidicates if the parameter values should be present 00584 on the plot 00585 """ 00586 from matplotlib import rc as rcp 00587 if not self.fitted: 00588 return 00589 #if not self._p or self._p.is_dead: 00590 if not (self._p and self._p._alive()): 00591 from asap.asapplotter import new_asaplot 00592 del self._p 00593 self._p = new_asaplot(rcParams['plotter.gui']) 00594 self._p.hold() 00595 self._p.clear() 00596 rcp('lines', linewidth=1) 00597 self._p.set_panels() 00598 self._p.palette(0) 00599 tlab = 'Spectrum' 00600 xlab = 'Abcissa' 00601 ylab = 'Ordinate' 00602 from numpy import ma,logical_not,logical_and,array 00603 m = self.mask 00604 if self.data: 00605 tlab = self.data._getsourcename(self._fittedrow) 00606 xlab = self.data._getabcissalabel(self._fittedrow) 00607 if self.data._getflagrow(self._fittedrow): 00608 m = [False] 00609 else: 00610 m = logical_and(self.mask, 00611 array(self.data._getmask(self._fittedrow), 00612 copy=False)) 00613 00614 ylab = self.data._get_ordinate_label() 00615 00616 colours = ["#777777","#dddddd","red","orange","purple","green", 00617 "magenta", "cyan"] 00618 nomask=True 00619 for i in range(len(m)): 00620 nomask = nomask and m[i] 00621 if len(m) == 1: 00622 m = m[0] 00623 invm = (not m) 00624 else: 00625 invm = logical_not(m) 00626 label0='Masked Region' 00627 label1='Spectrum' 00628 if ( nomask ): 00629 label0=label1 00630 else: 00631 y = ma.masked_array( self.y, mask = m ) 00632 self._p.palette(1,colours) 00633 self._p.set_line( label = label1 ) 00634 self._p.plot( self.x, y ) 00635 self._p.palette(0,colours) 00636 self._p.set_line(label=label0) 00637 y = ma.masked_array(self.y,mask=invm) 00638 self._p.plot(self.x, y) 00639 if residual: 00640 self._p.palette(7) 00641 self._p.set_line(label='Residual') 00642 y = ma.masked_array(self.get_residual(), 00643 mask=invm) 00644 self._p.plot(self.x, y) 00645 self._p.palette(2) 00646 if components is not None: 00647 cs = components 00648 if isinstance(components,int): cs = [components] 00649 if plotparms: 00650 self._p.text(0.15,0.15, 00651 str(self.get_parameters()['formatted']),size=8) 00652 n = len(self.components) 00653 self._p.palette(3) 00654 for c in cs: 00655 if 0 <= c < n: 00656 lab = self.fitfuncs[c]+str(c) 00657 self._p.set_line(label=lab) 00658 y = ma.masked_array(self.fitter.evaluate(c), mask=invm) 00659 00660 self._p.plot(self.x, y) 00661 elif c == -1: 00662 self._p.palette(2) 00663 self._p.set_line(label="Total Fit") 00664 y = ma.masked_array(self.fitter.getfit(), 00665 mask=invm) 00666 self._p.plot(self.x, y) 00667 else: 00668 self._p.palette(2) 00669 self._p.set_line(label='Fit') 00670 y = ma.masked_array(self.fitter.getfit(),mask=invm) 00671 self._p.plot(self.x, y) 00672 xlim=[min(self.x),max(self.x)] 00673 self._p.axes.set_xlim(xlim) 00674 self._p.set_axes('xlabel',xlab) 00675 self._p.set_axes('ylabel',ylab) 00676 self._p.set_axes('title',tlab) 00677 self._p.release() 00678 if (not rcParams['plotter.gui']): 00679 self._p.save(filename) 00680 00681 @asaplog_post_dec 00682 def auto_fit(self, insitu=None, plot=False): 00683 """ 00684 Return a scan where the function is applied to all rows for 00685 all Beams/IFs/Pols. 00686 00687 """ 00688 from asap import scantable 00689 if not isinstance(self.data, scantable) : 00690 msg = "Data is not a scantable" 00691 raise TypeError(msg) 00692 if insitu is None: insitu = rcParams['insitu'] 00693 if not insitu: 00694 scan = self.data.copy() 00695 else: 00696 scan = self.data 00697 rows = xrange(scan.nrow()) 00698 # Save parameters of baseline fits as a class attribute. 00699 # NOTICE: This does not reflect changes in scantable! 00700 if len(rows) > 0: self.blpars=[] 00701 asaplog.push("Fitting:") 00702 for r in rows: 00703 out = " Scan[%d] Beam[%d] IF[%d] Pol[%d] Cycle[%d]" % ( 00704 scan.getscan(r), 00705 scan.getbeam(r), 00706 scan.getif(r), 00707 scan.getpol(r), 00708 scan.getcycle(r) 00709 ) 00710 asaplog.push(out, False) 00711 self.x = scan._getabcissa(r) 00712 self.y = scan._getspectrum(r) 00713 #self.mask = mask_and(self.mask, scan._getmask(r)) 00714 if len(self.x) == len(self.mask): 00715 self.mask = mask_and(self.mask, self.data._getmask(row)) 00716 else: 00717 asaplog.push('lengths of data and mask are not the same. ' 00718 'preset mask will be ignored') 00719 asaplog.post('WARN','asapfit.fit') 00720 self.mask=self.data._getmask(row) 00721 self.data = None 00722 self.fit() 00723 x = self.get_parameters() 00724 fpar = self.get_parameters() 00725 if plot: 00726 self.plot(residual=True) 00727 x = raw_input("Accept fit ([y]/n): ") 00728 if x.upper() == 'N': 00729 self.blpars.append(None) 00730 continue 00731 scan._setspectrum(self.fitter.getresidual(), r) 00732 self.blpars.append(fpar) 00733 if plot: 00734 self._p.quit() 00735 del self._p 00736 self._p = None 00737 return scan
1.8.0
