Changeset 723
- Timestamp:
- 11/18/05 14:43:38 (19 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/python/asapfitter.py
r668 r723 1 1 import _asap 2 2 from asap import rcParams 3 from asap import print_log 3 4 4 5 class fitter: … … 6 7 The fitting class for ASAP. 7 8 """ 8 def _verbose(self, *args): 9 """ 10 Set stdout output. 11 """ 12 if type(args[0]) is bool: 13 self._vb = args[0] 14 return 15 elif len(args) == 0: 16 return self._vb 17 9 18 10 def __init__(self): 19 11 """ … … 31 23 self._fittedrow = 0 32 24 self._p = None 33 self._vb = True34 25 self._selection = None 35 26 … … 44 35 ydat: the ordinate values 45 36 mask: an optional mask 46 37 47 38 """ 48 39 self.fitted = False … … 64 55 """ 65 56 if not thescan: 66 print "Please give a correct scan" 57 msg = "Please give a correct scan" 58 if rcParams['verbose']: 59 print msg 60 return 61 else: 62 raise TypeError(msg) 67 63 self.fitted = False 68 64 self.data = thescan … … 84 80 fitter.set_function(poly=3) # will fit a 3rd order polynomial 85 81 """ 86 #default poly order 0 82 #default poly order 0 87 83 n=0 88 84 if kwargs.has_key('poly'): … … 96 92 self.components = [ 3 for i in range(n) ] 97 93 else: 98 print "Invalid function type." 99 return 94 msg = "Invalid function type." 95 if rcParams['verbose']: 96 print msg 97 return 98 else: 99 raise TypeError(msg) 100 100 101 self.fitter.setexpression(self.fitfunc,n) 101 102 return 102 103 103 104 def fit(self, row=0): 104 105 """ … … 112 113 f.set_scan(s) 113 114 f.set_function(poly=0) 114 f.fit(row=0) # fit first row 115 f.fit(row=0) # fit first row 115 116 """ 116 117 if ((self.x is None or self.y is None) and self.data is None) \ 117 118 or self.fitfunc is None: 118 print "Fitter not yet initialised. Please set data & fit function" 119 return 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 120 126 else: 121 127 if self.data is not None: 122 128 self.x = self.data._getabcissa(row) 123 129 self.y = self.data._getspectrum(row) 124 print "Fitting:" 125 vb = self.data._vb 126 self.data._vb = True 130 from asap import asaplog 131 asaplog.push("Fitting:") 127 132 self.selection = self.data.get_cursor() 128 self.data._vb = vb129 133 self.fitter.setdata(self.x, self.y, self.mask) 130 134 if self.fitfunc == 'gauss': … … 135 139 self.fitter.fit() 136 140 except RuntimeError, msg: 137 print msg 141 if rcParams['verbose']: 142 print msg 143 else: 144 raise 138 145 self._fittedrow = row 139 146 self.fitted = True 147 print_log() 140 148 return 141 149 … … 161 169 """ 162 170 if self.fitfunc is None: 163 print "Please specify a fitting function first." 164 return 171 msg = "Please specify a fitting function first." 172 if rcParams['verbose']: 173 print msg 174 return 175 else: 176 raise RuntimeError(msg) 165 177 if self.fitfunc == "gauss" and component is not None: 166 178 if not self.fitted: … … 169 181 fxd = list(zeros(len(pars))) 170 182 else: 171 pars = list(self.fitter.getparameters()) 183 pars = list(self.fitter.getparameters()) 172 184 fxd = list(self.fitter.getfixedparameters()) 173 185 i = 3*component … … 175 187 fxd[i:i+3] = fixed 176 188 params = pars 177 fixed = fxd 189 fixed = fxd 178 190 self.fitter.setparameters(params) 179 191 if fixed is not None: 180 192 self.fitter.setfixedparameters(fixed) 193 print_log() 181 194 return 182 195 … … 199 212 """ 200 213 if self.fitfunc != "gauss": 201 print "Function only operates on Gaussian components." 202 return 214 msg = "Function only operates on Gaussian components." 215 if rcParams['verbose']: 216 print msg 217 return 218 else: 219 raise ValueError(msg) 203 220 if 0 <= component < len(self.components): 204 221 self.set_parameters([peak, centre, fhwm], … … 206 223 component) 207 224 else: 208 print "Please select a valid component." 209 return 210 225 msg = "Please select a valid component." 226 if rcParams['verbose']: 227 print msg 228 return 229 else: 230 raise ValueError(msg) 231 211 232 def get_parameters(self, component=None): 212 233 """ … … 217 238 """ 218 239 if not self.fitted: 219 print "Not yet fitted." 240 msg = "Not yet fitted." 241 if rcParams['verbose']: 242 print msg 243 return 244 else: 245 raise RuntimeError(msg) 220 246 pars = list(self.fitter.getparameters()) 221 247 fixed = list(self.fitter.getfixedparameters()) 222 if component is not None: 248 if component is not None: 223 249 if self.fitfunc == "gauss": 224 250 i = 3*component … … 227 253 else: 228 254 cpars = pars 229 cfixed = fixed 255 cfixed = fixed 230 256 else: 231 257 cpars = pars 232 258 cfixed = fixed 233 259 fpars = self._format_pars(cpars, cfixed) 234 if self._vb:260 if rcParams['verbose']: 235 261 print fpars 236 262 return cpars, cfixed, fpars 237 263 238 264 def _format_pars(self, pars, fixed): 239 265 out = '' … … 259 285 i+=3 260 286 return out 261 287 262 288 def get_estimate(self): 263 289 """ … … 265 291 """ 266 292 pars = self.fitter.getestimate() 267 if self._vb:293 if rcParams['verbose']: 268 294 print self._format_pars(pars) 269 295 return pars 270 271 296 272 297 def get_residual(self): … … 275 300 """ 276 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: 277 344 print "Not yet fitted." 278 return self.fitter.getresidual() 279 280 def get_chi2(self): 281 """ 282 Return chi^2. 283 """ 284 if not self.fitted: 285 print "Not yet fitted." 286 ch2 = self.fitter.getchi2() 287 if self._vb: 288 print 'Chi^2 = %3.3f' % (ch2) 289 return ch2 290 291 def get_fit(self): 292 """ 293 Return the fitted ordinate values. 294 """ 295 if not self.fitted: 296 print "Not yet fitted." 297 return self.fitter.getfit() 298 299 def commit(self): 300 """ 301 Return a new scan where the fits have been commited (subtracted) 302 """ 303 if not self.fitted: 304 print "Not yet fitted." 345 msg = "Not yet fitted." 346 if rcParams['verbose']: 347 print msg 348 return 349 else: 350 raise RuntimeError(msg) 305 351 if self.data is not scantable: 306 print "Only works with scantables" 307 return 352 msg = "Not a scantable" 353 if rcParams['verbose']: 354 print msg 355 return 356 else: 357 raise TypeError(msg) 308 358 scan = self.data.copy() 309 359 scan._setspectrum(self.fitter.getresidual()) 310 311 def plot(self, residual=False, components=None, plotparms=False): 360 print_log() 361 362 def plot(self, residual=False, components=None, plotparms=False, filename=None): 312 363 """ 313 364 Plot the last fit. … … 323 374 if not self.fitted: 324 375 return 325 if not self._p: 326 from asap.asaplot import ASAPlot 327 self._p = ASAPlot() 328 if self._p.is_dead: 329 from asap.asaplot import ASAPlot 330 self._p = ASAPlot() 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() 331 383 self._p.clear() 332 384 self._p.set_panels() 333 385 self._p.palette(0) 334 386 tlab = 'Spectrum' 335 xlab = 'Abcissa' 387 xlab = 'Abcissa' 336 388 m = () 337 389 if self.data: … … 365 417 self._p.palette(2) 366 418 self._p.set_line(label="Total Fit") 367 self._p.plot(self.x, self.get_fit(), m) 419 self._p.plot(self.x, self.get_fit(), m) 368 420 else: 369 421 self._p.palette(2) 370 422 self._p.set_line(label='Fit') 371 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) 372 426 self._p.set_axes('xlabel',xlab) 373 427 self._p.set_axes('ylabel',ylab) 374 428 self._p.set_axes('title',tlab) 375 429 self._p.release() 430 if (not rcParams['plotter.gui']): 431 self._p.save(filename) 432 print_log() 376 433 377 434 def auto_fit(self, insitu=None): … … 379 436 Return a scan where the function is applied to all rows for 380 437 all Beams/IFs/Pols. 381 438 382 439 """ 383 440 from asap import scantable 384 441 if not isinstance(self.data, scantable) : 385 print "Only works with scantables" 386 return 442 msg = "Data is not a scantable" 443 if rcParams['verbose']: 444 print msg 445 return 446 else: 447 raise TypeError(msg) 387 448 if insitu is None: insitu = rcParams['insitu'] 388 449 if not insitu: … … 390 451 else: 391 452 scan = self.data 392 vb = scan._vb393 scan._vb = False394 453 sel = scan.get_cursor() 395 454 rows = range(scan.nrow()) 455 from asap import asaplog 396 456 for i in range(scan.nbeam()): 397 457 scan.setbeam(i) … … 400 460 for k in range(scan.npol()): 401 461 scan.setpol(k) 402 if self._vb:403 print "Fitting:"404 print 'Beam[%d], IF[%d], Pol[%d]' % (i,j,k)462 asaplog.push("Fitting:") 463 out = 'Beam[%d], IF[%d], Pol[%d]' % (i,j,k) 464 asaplog.push(out) 405 465 for iRow in rows: 406 466 self.x = scan._getabcissa(iRow) 407 467 self.y = scan._getspectrum(iRow) 408 468 self.data = None 409 self.fit() 469 self.fit() 410 470 x = self.get_parameters() 411 471 scan._setspectrum(self.fitter.getresidual(),iRow) 412 472 scan.set_cursor(sel[0],sel[1],sel[2]) 413 scan._vb = vb473 print_log() 414 474 if not insitu: 415 475 return scan
Note:
See TracChangeset
for help on using the changeset viewer.