source: trunk/python/asapfitter.py@ 1220

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

Merge from Release2.1.0b tag

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 19.8 KB
RevLine 
[113]1import _asap
[259]2from asap import rcParams
[723]3from asap import print_log
[1134]4from asap import NUM
[113]5
6class fitter:
7 """
8 The fitting class for ASAP.
9 """
[723]10
[113]11 def __init__(self):
12 """
13 Create a fitter object. No state is set.
14 """
15 self.fitter = _asap.fitter()
16 self.x = None
17 self.y = None
18 self.mask = None
19 self.fitfunc = None
[515]20 self.fitfuncs = None
[113]21 self.fitted = False
22 self.data = None
[515]23 self.components = 0
24 self._fittedrow = 0
[113]25 self._p = None
[515]26 self._selection = None
[113]27
28 def set_data(self, xdat, ydat, mask=None):
29 """
[158]30 Set the absissa and ordinate for the fit. Also set the mask
[113]31 indicationg valid points.
32 This can be used for data vectors retrieved from a scantable.
33 For scantable fitting use 'fitter.set_scan(scan, mask)'.
34 Parameters:
[158]35 xdat: the abcissa values
[113]36 ydat: the ordinate values
37 mask: an optional mask
[723]38
[113]39 """
40 self.fitted = False
41 self.x = xdat
42 self.y = ydat
43 if mask == None:
[1134]44 self.mask = NUM.ones(len(xdat))
[113]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:
[723]57 msg = "Please give a correct scan"
58 if rcParams['verbose']:
59 print msg
60 return
61 else:
62 raise TypeError(msg)
[113]63 self.fitted = False
64 self.data = thescan
[1075]65 self.mask = None
[113]66 if mask is None:
[1134]67 self.mask = NUM.ones(self.data.nchan())
[113]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 """
[723]82 #default poly order 0
[515]83 n=0
[113]84 if kwargs.has_key('poly'):
85 self.fitfunc = 'poly'
86 n = kwargs.get('poly')
[515]87 self.components = [n]
[113]88 elif kwargs.has_key('gauss'):
89 n = kwargs.get('gauss')
90 self.fitfunc = 'gauss'
[515]91 self.fitfuncs = [ 'gauss' for i in range(n) ]
92 self.components = [ 3 for i in range(n) ]
93 else:
[723]94 msg = "Invalid function type."
95 if rcParams['verbose']:
96 print msg
97 return
98 else:
99 raise TypeError(msg)
100
[113]101 self.fitter.setexpression(self.fitfunc,n)
102 return
[723]103
[1075]104 def fit(self, row=0, estimate=False):
[113]105 """
106 Execute the actual fitting process. All the state has to be set.
107 Parameters:
[1075]108 row: specify the row in the scantable
109 estimate: auto-compute an initial parameter set (default False)
110 This can be used to compute estimates even if fit was
111 called before.
[113]112 Example:
[515]113 s = scantable('myscan.asap')
114 s.set_cursor(thepol=1) # select second pol
[113]115 f = fitter()
116 f.set_scan(s)
117 f.set_function(poly=0)
[723]118 f.fit(row=0) # fit first row
[113]119 """
120 if ((self.x is None or self.y is None) and self.data is None) \
121 or self.fitfunc is None:
[723]122 msg = "Fitter not yet initialised. Please set data & fit function"
123 if rcParams['verbose']:
124 print msg
125 return
126 else:
127 raise RuntimeError(msg)
128
[113]129 else:
130 if self.data is not None:
[515]131 self.x = self.data._getabcissa(row)
132 self.y = self.data._getspectrum(row)
[723]133 from asap import asaplog
134 asaplog.push("Fitting:")
[943]135 i = row
136 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))
[1075]137 asaplog.push(out,False)
[515]138 self.fitter.setdata(self.x, self.y, self.mask)
[113]139 if self.fitfunc == 'gauss':
140 ps = self.fitter.getparameters()
[1075]141 if len(ps) == 0 or estimate:
[113]142 self.fitter.estimate()
[626]143 try:
[1075]144 converged = self.fitter.fit()
145 if not converged:
146 raise RuntimeError,"Fit didn't converge."
[626]147 except RuntimeError, msg:
[723]148 if rcParams['verbose']:
149 print msg
150 else:
151 raise
[515]152 self._fittedrow = row
[113]153 self.fitted = True
[723]154 print_log()
[113]155 return
156
[515]157 def store_fit(self):
[526]158 """
159 Store the fit parameters in the scantable.
160 """
[515]161 if self.fitted and self.data is not None:
162 pars = list(self.fitter.getparameters())
163 fixed = list(self.fitter.getfixedparameters())
[975]164 from asap.asapfit import asapfit
165 fit = asapfit()
166 fit.setparameters(pars)
167 fit.setfixedparameters(fixed)
168 fit.setfunctions(self.fitfuncs)
169 fit.setcomponents(self.components)
170 fit.setframeinfo(self.data._getcoordinfo())
171 self.data._addfit(fit,self._fittedrow)
[515]172
[1017]173 #def set_parameters(self, params, fixed=None, component=None):
174 def set_parameters(self,*args,**kwargs):
[526]175 """
176 Set the parameters to be fitted.
177 Parameters:
178 params: a vector of parameters
179 fixed: a vector of which parameters are to be held fixed
180 (default is none)
181 component: in case of multiple gaussians, the index of the
182 component
[1017]183 """
184 component = None
185 fixed = None
186 params = None
[1031]187
[1017]188 if len(args) and isinstance(args[0],dict):
189 kwargs = args[0]
190 if kwargs.has_key("fixed"): fixed = kwargs["fixed"]
191 if kwargs.has_key("params"): params = kwargs["params"]
192 if len(args) == 2 and isinstance(args[1], int):
193 component = args[1]
[515]194 if self.fitfunc is None:
[723]195 msg = "Please specify a fitting function first."
196 if rcParams['verbose']:
197 print msg
198 return
199 else:
200 raise RuntimeError(msg)
[515]201 if self.fitfunc == "gauss" and component is not None:
[1017]202 if not self.fitted and sum(self.fitter.getparameters()) == 0:
[1134]203 pars = list(NUM.zeros(len(self.components)*3))
204 fxd = list(NUM.zeros(len(pars)))
[515]205 else:
[723]206 pars = list(self.fitter.getparameters())
[515]207 fxd = list(self.fitter.getfixedparameters())
208 i = 3*component
209 pars[i:i+3] = params
210 fxd[i:i+3] = fixed
211 params = pars
[723]212 fixed = fxd
[113]213 self.fitter.setparameters(params)
214 if fixed is not None:
215 self.fitter.setfixedparameters(fixed)
[723]216 print_log()
[113]217 return
[515]218
[1217]219 def set_gauss_parameters(self, peak, centre, fwhm,
[1017]220 peakfixed=0, centerfixed=0,
[1217]221 fwhmfixed=0,
[515]222 component=0):
[113]223 """
[515]224 Set the Parameters of a 'Gaussian' component, set with set_function.
225 Parameters:
226 peak, centre, fhwm: The gaussian parameters
227 peakfixed,
228 centerfixed,
[1217]229 fwhmfixed: Optional parameters to indicate if
[515]230 the paramters should be held fixed during
231 the fitting process. The default is to keep
232 all parameters flexible.
[526]233 component: The number of the component (Default is the
234 component 0)
[515]235 """
236 if self.fitfunc != "gauss":
[723]237 msg = "Function only operates on Gaussian components."
238 if rcParams['verbose']:
239 print msg
240 return
241 else:
242 raise ValueError(msg)
[515]243 if 0 <= component < len(self.components):
[1217]244 d = {'params':[peak, centre, fwhm],
245 'fixed':[peakfixed, centerfixed, fwhmfixed]}
[1017]246 self.set_parameters(d, component)
[515]247 else:
[723]248 msg = "Please select a valid component."
249 if rcParams['verbose']:
250 print msg
251 return
252 else:
253 raise ValueError(msg)
254
[975]255 def get_area(self, component=None):
256 """
257 Return the area under the fitted gaussian component.
258 Parameters:
259 component: the gaussian component selection,
260 default (None) is the sum of all components
261 Note:
262 This will only work for gaussian fits.
263 """
264 if not self.fitted: return
265 if self.fitfunc == "gauss":
266 pars = list(self.fitter.getparameters())
267 from math import log,pi,sqrt
268 fac = sqrt(pi/log(16.0))
269 areas = []
270 for i in range(len(self.components)):
271 j = i*3
272 cpars = pars[j:j+3]
273 areas.append(fac * cpars[0] * cpars[2])
274 else:
275 return None
276 if component is not None:
277 return areas[component]
278 else:
279 return sum(areas)
280
[1075]281 def get_errors(self, component=None):
[515]282 """
[1075]283 Return the errors in the parameters.
284 Parameters:
285 component: get the errors for the specified component
286 only, default is all components
287 """
288 if not self.fitted:
289 msg = "Not yet fitted."
290 if rcParams['verbose']:
291 print msg
292 return
293 else:
294 raise RuntimeError(msg)
295 errs = list(self.fitter.geterrors())
296 cerrs = errs
297 if component is not None:
298 if self.fitfunc == "gauss":
299 i = 3*component
300 if i < len(errs):
301 cerrs = errs[i:i+3]
302 return cerrs
303
304 def get_parameters(self, component=None, errors=False):
305 """
[113]306 Return the fit paramters.
[526]307 Parameters:
308 component: get the parameters for the specified component
309 only, default is all components
[113]310 """
311 if not self.fitted:
[723]312 msg = "Not yet fitted."
313 if rcParams['verbose']:
314 print msg
315 return
316 else:
317 raise RuntimeError(msg)
[113]318 pars = list(self.fitter.getparameters())
319 fixed = list(self.fitter.getfixedparameters())
[1075]320 errs = list(self.fitter.geterrors())
[1039]321 area = []
[723]322 if component is not None:
[515]323 if self.fitfunc == "gauss":
324 i = 3*component
325 cpars = pars[i:i+3]
326 cfixed = fixed[i:i+3]
[1075]327 cerrs = errs[i:i+3]
[1039]328 a = self.get_area(component)
329 area = [a for i in range(3)]
[515]330 else:
331 cpars = pars
[723]332 cfixed = fixed
[1075]333 cerrs = errs
[515]334 else:
335 cpars = pars
336 cfixed = fixed
[1075]337 cerrs = errs
[1039]338 if self.fitfunc == "gauss":
339 for c in range(len(self.components)):
340 a = self.get_area(c)
341 area += [a for i in range(3)]
[1088]342 fpars = self._format_pars(cpars, cfixed, errors and cerrs, area)
[723]343 if rcParams['verbose']:
[515]344 print fpars
[1075]345 return {'params':cpars, 'fixed':cfixed, 'formatted': fpars,
346 'errors':cerrs}
[723]347
[1075]348 def _format_pars(self, pars, fixed, errors, area):
[113]349 out = ''
350 if self.fitfunc == 'poly':
351 c = 0
[515]352 for i in range(len(pars)):
353 fix = ""
354 if fixed[i]: fix = "(fixed)"
[1088]355 if errors :
356 out += ' p%d%s= %3.6f (%1.6f),' % (c,fix,pars[i], errors[i])
357 else:
358 out += ' p%d%s= %3.6f,' % (c,fix,pars[i])
[113]359 c+=1
[515]360 out = out[:-1] # remove trailing ','
[113]361 elif self.fitfunc == 'gauss':
362 i = 0
363 c = 0
[515]364 aunit = ''
365 ounit = ''
[113]366 if self.data:
[515]367 aunit = self.data.get_unit()
368 ounit = self.data.get_fluxunit()
[113]369 while i < len(pars):
[1039]370 if len(area):
371 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]372 else:
373 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]374 c+=1
375 i+=3
376 return out
[723]377
[113]378 def get_estimate(self):
379 """
[515]380 Return the parameter estimates (for non-linear functions).
[113]381 """
382 pars = self.fitter.getestimate()
[943]383 fixed = self.fitter.getfixedparameters()
[723]384 if rcParams['verbose']:
[1017]385 print self._format_pars(pars,fixed,None)
[113]386 return pars
387
388 def get_residual(self):
389 """
390 Return the residual of the fit.
391 """
392 if not self.fitted:
[723]393 msg = "Not yet fitted."
394 if rcParams['verbose']:
395 print msg
396 return
397 else:
398 raise RuntimeError(msg)
[113]399 return self.fitter.getresidual()
400
401 def get_chi2(self):
402 """
403 Return chi^2.
404 """
405 if not self.fitted:
[723]406 msg = "Not yet fitted."
407 if rcParams['verbose']:
408 print msg
409 return
410 else:
411 raise RuntimeError(msg)
[113]412 ch2 = self.fitter.getchi2()
[723]413 if rcParams['verbose']:
[113]414 print 'Chi^2 = %3.3f' % (ch2)
[723]415 return ch2
[113]416
417 def get_fit(self):
418 """
419 Return the fitted ordinate values.
420 """
421 if not self.fitted:
[723]422 msg = "Not yet fitted."
423 if rcParams['verbose']:
424 print msg
425 return
426 else:
427 raise RuntimeError(msg)
[113]428 return self.fitter.getfit()
429
430 def commit(self):
431 """
[526]432 Return a new scan where the fits have been commited (subtracted)
[113]433 """
434 if not self.fitted:
[723]435 msg = "Not yet fitted."
436 if rcParams['verbose']:
437 print msg
438 return
439 else:
440 raise RuntimeError(msg)
[975]441 from asap import scantable
442 if not isinstance(self.data, scantable):
[723]443 msg = "Not a scantable"
444 if rcParams['verbose']:
445 print msg
446 return
447 else:
448 raise TypeError(msg)
[113]449 scan = self.data.copy()
[259]450 scan._setspectrum(self.fitter.getresidual())
[723]451 print_log()
[1092]452 return scan
[113]453
[723]454 def plot(self, residual=False, components=None, plotparms=False, filename=None):
[113]455 """
456 Plot the last fit.
457 Parameters:
458 residual: an optional parameter indicating if the residual
459 should be plotted (default 'False')
[526]460 components: a list of components to plot, e.g [0,1],
461 -1 plots the total fit. Default is to only
462 plot the total fit.
463 plotparms: Inidicates if the parameter values should be present
464 on the plot
[113]465 """
466 if not self.fitted:
467 return
[723]468 if not self._p or self._p.is_dead:
469 if rcParams['plotter.gui']:
470 from asap.asaplotgui import asaplotgui as asaplot
471 else:
472 from asap.asaplot import asaplot
473 self._p = asaplot()
474 self._p.hold()
[113]475 self._p.clear()
[515]476 self._p.set_panels()
[652]477 self._p.palette(0)
[113]478 tlab = 'Spectrum'
[723]479 xlab = 'Abcissa'
[1017]480 ylab = 'Ordinate'
[1217]481 from matplotlib.numerix import ma,logical_not,array
482 m = NUM.ones(len(self.x))
483
[113]484 if self.data:
[515]485 tlab = self.data._getsourcename(self._fittedrow)
486 xlab = self.data._getabcissalabel(self._fittedrow)
487 m = self.data._getmask(self._fittedrow)
[626]488 ylab = self.data._get_ordinate_label()
[515]489
[1075]490 colours = ["#777777","#dddddd","red","orange","purple","green","magenta", "cyan"]
[652]491 self._p.palette(0,colours)
[515]492 self._p.set_line(label='Spectrum')
[1116]493 y = ma.masked_array(self.y,mask=logical_not(array(m,copy=False)))
[1088]494 self._p.plot(self.x, y)
[113]495 if residual:
[652]496 self._p.palette(1)
[515]497 self._p.set_line(label='Residual')
[1116]498 y = ma.masked_array(self.get_residual(),
499 mask=logical_not(array(m,copy=False)))
[1088]500 self._p.plot(self.x, y)
[652]501 self._p.palette(2)
[515]502 if components is not None:
503 cs = components
504 if isinstance(components,int): cs = [components]
[526]505 if plotparms:
[1031]506 self._p.text(0.15,0.15,str(self.get_parameters()['formatted']),size=8)
[515]507 n = len(self.components)
[652]508 self._p.palette(3)
[515]509 for c in cs:
510 if 0 <= c < n:
511 lab = self.fitfuncs[c]+str(c)
512 self._p.set_line(label=lab)
[1116]513 y = ma.masked_array(self.fitter.evaluate(c),
514 mask=logical_not(array(m,copy=False)))
[1088]515
516 self._p.plot(self.x, y)
[515]517 elif c == -1:
[652]518 self._p.palette(2)
[515]519 self._p.set_line(label="Total Fit")
[1116]520 y = ma.masked_array(self.fitter.getfit(),
521 mask=logical_not(array(m,copy=False)))
[1088]522 self._p.plot(self.x, y)
[515]523 else:
[652]524 self._p.palette(2)
[515]525 self._p.set_line(label='Fit')
[1116]526 y = ma.masked_array(self.fitter.getfit(),
527 mask=logical_not(array(m,copy=False)))
[1088]528 self._p.plot(self.x, y)
[723]529 xlim=[min(self.x),max(self.x)]
530 self._p.axes.set_xlim(xlim)
[113]531 self._p.set_axes('xlabel',xlab)
532 self._p.set_axes('ylabel',ylab)
533 self._p.set_axes('title',tlab)
534 self._p.release()
[723]535 if (not rcParams['plotter.gui']):
536 self._p.save(filename)
537 print_log()
[113]538
[1061]539 def auto_fit(self, insitu=None, plot=False):
[113]540 """
[515]541 Return a scan where the function is applied to all rows for
542 all Beams/IFs/Pols.
[723]543
[113]544 """
545 from asap import scantable
[515]546 if not isinstance(self.data, scantable) :
[723]547 msg = "Data is not a scantable"
548 if rcParams['verbose']:
549 print msg
550 return
551 else:
552 raise TypeError(msg)
[259]553 if insitu is None: insitu = rcParams['insitu']
554 if not insitu:
555 scan = self.data.copy()
556 else:
557 scan = self.data
[880]558 rows = xrange(scan.nrow())
[723]559 from asap import asaplog
[876]560 asaplog.push("Fitting:")
561 for r in rows:
[1031]562 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))
[880]563 asaplog.push(out, False)
[876]564 self.x = scan._getabcissa(r)
565 self.y = scan._getspectrum(r)
566 self.data = None
567 self.fit()
568 x = self.get_parameters()
[1061]569 if plot:
570 self.plot(residual=True)
571 x = raw_input("Accept fit ([y]/n): ")
572 if x.upper() == 'N':
573 continue
[880]574 scan._setspectrum(self.fitter.getresidual(), r)
[1061]575 if plot:
576 self._p.unmap()
577 self._p = None
[876]578 print_log()
579 return scan
[794]580
Note: See TracBrowser for help on using the repository browser.