source: trunk/python/asapfitter.py@ 250

Last change on this file since 250 was 190, checked in by mar637, 20 years ago

Bug fix: Handling destruction of window via self._is_dead

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 8.5 KB
Line 
1import _asap
2
3class fitter:
4 """
5 The fitting class for ASAP.
6 """
7 def _verbose(self, *args):
8 """
9 Set stdout output.
10 """
11 if type(args[0]) is bool:
12 self._vb = args[0]
13 return
14 elif len(args) == 0:
15 return self._vb
16
17 def __init__(self):
18 """
19 Create a fitter object. No state is set.
20 """
21 self.fitter = _asap.fitter()
22 self.x = None
23 self.y = None
24 self.mask = None
25 self.fitfunc = None
26 self.fitted = False
27 self.data = None
28 self._p = None
29 self._vb = True
30
31 def set_data(self, xdat, ydat, mask=None):
32 """
33 Set the absissa and ordinate for the fit. Also set the mask
34 indicationg valid points.
35 This can be used for data vectors retrieved from a scantable.
36 For scantable fitting use 'fitter.set_scan(scan, mask)'.
37 Parameters:
38 xdat: the abcissa values
39 ydat: the ordinate values
40 mask: an optional mask
41
42 """
43 self.fitted = False
44 self.x = xdat
45 self.y = ydat
46 if mask == None:
47 from numarray import ones
48 self.mask = ones(len(xdat))
49 else:
50 self.mask = mask
51 return
52
53 def set_scan(self, thescan=None, mask=None):
54 """
55 Set the 'data' (a scantable) of the fitter.
56 Parameters:
57 thescan: a scantable
58 mask: a msk retireved from the scantable
59 """
60 if not thescan:
61 print "Please give a correct scan"
62 self.fitted = False
63 self.data = thescan
64 if mask is None:
65 from numarray import ones
66 self.mask = ones(self.data.nchan())
67 else:
68 self.mask = mask
69 return
70
71 def set_function(self, **kwargs):
72 """
73 Set the function to be fit.
74 Parameters:
75 poly: use a polynomial of the order given
76 gauss: fit the number of gaussian specified
77 Example:
78 fitter.set_function(gauss=2) # will fit two gaussians
79 fitter.set_function(poly=3) # will fit a 3rd order polynomial
80 """
81 #default poly order 0
82 self.fitfunc = 'poly'
83 n=0
84 if kwargs.has_key('poly'):
85 self.fitfunc = 'poly'
86 n = kwargs.get('poly')
87 elif kwargs.has_key('gauss'):
88 n = kwargs.get('gauss')
89 self.fitfunc = 'gauss'
90
91 self.fitter.setexpression(self.fitfunc,n)
92 return
93
94 def fit(self):
95 """
96 Execute the actual fitting process. All the state has to be set.
97 Parameters:
98 none
99 Example:
100 s= scantable('myscan.asap')
101 f = fitter()
102 f.set_scan(s)
103 f.set_function(poly=0)
104 f.fit()
105 """
106 if ((self.x is None or self.y is None) and self.data is None) \
107 or self.fitfunc is None:
108 print "Fitter not yet initialised. Please set data & fit function"
109 return
110 else:
111 if self.data is not None:
112 self.x = self.data.getabcissa()
113 self.y = self.data.getspectrum()
114 print "Fitting:"
115 vb = self.data._verbose
116 self.data._verbose(True)
117 s = self.data.get_selection()
118 self.data._verbose(vb)
119
120 self.fitter.setdata(self.x,self.y,self.mask)
121 if self.fitfunc == 'gauss':
122 ps = self.fitter.getparameters()
123 if len(ps) == 0:
124 self.fitter.estimate()
125 self.fitter.fit()
126 self.fitted = True
127 return
128
129 def set_parameters(self, params, fixed=None):
130 self.fitter.setparameters(params)
131 if fixed is not None:
132 self.fitter.setfixedparameters(fixed)
133 return
134
135 def get_parameters(self):
136 """
137 Return the fit paramters.
138
139 """
140 if not self.fitted:
141 print "Not yet fitted."
142 pars = list(self.fitter.getparameters())
143 fixed = list(self.fitter.getfixedparameters())
144 if self._vb:
145 print self._format_pars(pars)
146 return pars,fixed
147
148 def _format_pars(self, pars):
149 out = ''
150 if self.fitfunc == 'poly':
151 c = 0
152 for i in pars:
153 out += ' p%d = %3.3f, ' % (c,i)
154 c+=1
155 elif self.fitfunc == 'gauss':
156 i = 0
157 c = 0
158 unit = ''
159 if self.data:
160 unit = self.data.get_unit()
161 while i < len(pars):
162 out += ' %d: peak = %3.3f , centre = %3.3f %s, FWHM = %3.3f %s \n' % (c,pars[i],pars[i+1],unit,pars[i+2],unit)
163 c+=1
164 i+=3
165 return out
166
167 def get_estimate(self):
168 """
169 Return the paramter estimates (for non-linear functions).
170 """
171 pars = self.fitter.getestimate()
172 if self._vb:
173 print self._format_pars(pars)
174 return pars
175
176
177 def get_residual(self):
178 """
179 Return the residual of the fit.
180 """
181 if not self.fitted:
182 print "Not yet fitted."
183 return self.fitter.getresidual()
184
185 def get_chi2(self):
186 """
187 Return chi^2.
188 """
189
190 if not self.fitted:
191 print "Not yet fitted."
192 ch2 = self.fitter.getchi2()
193 if self._vb:
194 print 'Chi^2 = %3.3f' % (ch2)
195 return ch2
196
197 def get_fit(self):
198 """
199 Return the fitted ordinate values.
200 """
201 if not self.fitted:
202 print "Not yet fitted."
203 return self.fitter.getfit()
204
205 def commit(self):
206 """
207 Return a new scan where teh fits have been commited.
208 """
209 if not self.fitted:
210 print "Not yet fitted."
211 if self.data is not scantable:
212 print "Only works with scantables"
213 return
214 scan = self.data.copy()
215 scan.setspectrum(self.fitter.getresidual())
216
217 def plot(self, residual=False):
218 """
219 Plot the last fit.
220 Parameters:
221 residual: an optional parameter indicating if the residual
222 should be plotted (default 'False')
223 """
224 if not self.fitted:
225 return
226 if not self._p:
227 from asap.asaplot import ASAPlot
228 self._p = ASAPlot()
229 if self._.is_dead:
230 from asap.asaplot import ASAPlot
231 self._p = ASAPlot()
232 self._p.clear()
233 tlab = 'Spectrum'
234 xlab = 'Abcissa'
235 if self.data:
236 tlab = self.data._getsourcename(0)
237 xlab = self.data.getabcissalabel(0)
238 ylab = r'Flux'
239 m = self.data.getmask(0)
240 self._p.set_line(colour='blue',label='Spectrum')
241 self._p.plot(self.x, self.y, m)
242 if residual:
243 self._p.set_line(colour='green',label='Residual')
244 self._p.plot(self.x, self.get_residual(), m)
245 self._p.set_line(colour='red',label='Fit')
246 self._p.plot(self.x, self.get_fit(), m)
247
248 self._p.set_axes('xlabel',xlab)
249 self._p.set_axes('ylabel',ylab)
250 self._p.set_axes('title',tlab)
251 self._p.release()
252
253
254 def auto_fit(self):
255 """
256 Return a scan where the function is applied to all rows for all Beams/IFs/Pols.
257
258 """
259 from asap import scantable
260 if not isinstance(self.data,scantable) :
261 print "Only works with scantables"
262 return
263 scan = self.data.copy()
264 vb = scan._verbose
265 scan._verbose(False)
266 sel = scan.get_selection()
267 rows = range(scan.nrow())
268 for i in range(scan.nbeam()):
269 scan.setbeam(i)
270 for j in range(scan.nif()):
271 scan.setif(j)
272 for k in range(scan.npol()):
273 scan.setpol(k)
274 if self._vb:
275 print "Fitting:"
276 print 'Beam[%d], IF[%d], Pol[%d]' % (i,j,k)
277 for iRow in rows:
278 self.x = scan.getabcissa(iRow)
279 self.y = scan.getspectrum(iRow)
280 self.data = None
281 self.fit()
282 x = self.get_parameters()
283 scan.setspectrum(self.fitter.getresidual(),iRow)
284 scan.set_selection(sel[0],sel[1],sel[2])
285 scan._verbose(vb)
286 return scan
Note: See TracBrowser for help on using the repository browser.