source: tags/asap1/python/asapfitter.py@ 1475

Last change on this file since 1475 was 794, checked in by mar637, 19 years ago

update from Release12

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