Context Navigation

source: trunk/python/asapfitter.py @ 876

Visit:

Last change on this file since 876 was 876, checked in by mar637, 18 years ago
move to asap2 container
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 15.5 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: