source: tags/asap2beta/python/asapfitter.py@ 2803

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

Completed Ticket #7 - storing of fits.

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