source: trunk/python/asapfitter.py@ 917

Last change on this file since 917 was 880, checked in by mar637, 19 years ago

added linefinder

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