source: branches/alma/python/asapfitter.py@ 2804

Last change on this file since 2804 was 1782, checked in by Kana Sugimoto, 14 years ago

New Development: Yes

JIRA Issue: No (merge)

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s): single dish package

Description:

Merged some unimplemented developments from the trunk.
[Details]

  1. Restored PKSmsg.cc, PKSmsg.h to alma branch from trunk
  2. Updated pks_matchs.cc and src/SConscript to ones in trunk
  3. uncommented "@print_log_dec"
  4. implemented bug fixes to add(), sub(), mul(), and div() in scantable.py
  5. implemented get_region() and get_point() in asaplotbase.py
  6. activated interactive keyword for text(), arrow(), annotate(), axvline(), axhline(), axvspan(), and axhspan() in asapplotter.py


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