source: trunk/python/asapfitter.py@ 1075

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

added errors to class; added fit convergence test; optional non-estimate in fit

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