Context Navigation

source: trunk/python/asapfitter.py@ 1046

Visit:

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

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

Download in other formats: