source: trunk/python/asapfitter.py@ 958

Last change on this file since 958 was 943, checked in by mar637, 19 years ago

fixed a couple of variable mix ups

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