qsmearing.py @ fc292d3

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload

Last change on this file since fc292d3 was 2ca00c4, checked in by Jae Cho <jhjcho@…>, 14 years ago
increase smearer npts for binning to get a reasonable smearing for USANS data
Property mode set to `100644`
File size: 20.8 KB

Line
1
2	#####################################################################
3	#This software was developed by the University of Tennessee as part of the
4	#Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
5	#project funded by the US National Science Foundation.
6	#See the license text in license.txt
7	#copyright 2008, University of Tennessee
8	######################################################################
9	import numpy
10	import math
11	import logging
12	import sys
13	import DataLoader.extensions.smearer as smearer
14	from DataLoader.smearing_2d import Smearer2D
15
16	def smear_selection(data1D, model = None):
17	"""
18	Creates the right type of smearer according
19	to the data.
20
21	The canSAS format has a rule that either
22	slit smearing data OR resolution smearing data
23	is available.
24
25	For the present purpose, we choose the one that
26	has none-zero data. If both slit and resolution
27	smearing arrays are filled with good data
28	(which should not happen), then we choose the
29	resolution smearing data.
30
31	:param data1D: Data1D object
32	:param model: sans.model instance
33	"""
34	# Sanity check. If we are not dealing with a SANS Data1D
35	# object, just return None
36	if data1D.__class__.__name__ not in ['Data1D', 'Theory1D']:
37	if data1D == None:
38	return None
39	elif data1D.dqx_data == None or data1D.dqy_data == None:
40	return None
41	return Smearer2D(data1D)
42
43	if not hasattr(data1D, "dx") and not hasattr(data1D, "dxl")\
44	and not hasattr(data1D, "dxw"):
45	return None
46
47	# Look for resolution smearing data
48	_found_resolution = False
49	if data1D.dx is not None and len(data1D.dx) == len(data1D.x):
50
51	# Check that we have non-zero data
52	if data1D.dx[0] > 0.0:
53	_found_resolution = True
54	#print "_found_resolution",_found_resolution
55	#print "data1D.dx[0]",data1D.dx[0],data1D.dxl[0]
56	# If we found resolution smearing data, return a QSmearer
57	if _found_resolution == True:
58	return QSmearer(data1D, model)
59
60	# Look for slit smearing data
61	_found_slit = False
62	if data1D.dxl is not None and len(data1D.dxl) == len(data1D.x) \
63	and data1D.dxw is not None and len(data1D.dxw) == len(data1D.x):
64
65	# Check that we have non-zero data
66	if data1D.dxl[0] > 0.0 or data1D.dxw[0] > 0.0:
67	_found_slit = True
68
69	# Sanity check: all data should be the same as a function of Q
70	for item in data1D.dxl:
71	if data1D.dxl[0] != item:
72	_found_resolution = False
73	break
74
75	for item in data1D.dxw:
76	if data1D.dxw[0] != item:
77	_found_resolution = False
78	break
79	# If we found slit smearing data, return a slit smearer
80	if _found_slit == True:
81	return SlitSmearer(data1D, model)
82	return None
83
84
85	class _BaseSmearer(object):
86
87	def __init__(self):
88	self.nbins = 0
89	self.nbins_low = 0
90	self.nbins_high = 0
91	self._weights = None
92	## Internal flag to keep track of C++ smearer initialization
93	self._init_complete = False
94	self._smearer = None
95	self.model = None
96
97	def __deepcopy__(self, memo={}):
98	"""
99	Return a valid copy of self.
100	Avoid copying the _smearer C object and force a matrix recompute
101	when the copy is used.
102	"""
103	result = _BaseSmearer()
104	result.nbins = self.nbins
105	return result
106
107	def _compute_matrix(self):
108	"""
109	"""
110	return NotImplemented
111
112	def get_bin_range(self, q_min=None, q_max=None):
113	"""
114
115	:param q_min: minimum q-value to smear
116	:param q_max: maximum q-value to smear
117
118	"""
119	# If this is the first time we call for smearing,
120	# initialize the C++ smearer object first
121	if not self._init_complete:
122	self._initialize_smearer()
123	if q_min == None:
124	q_min = self.min
125	if q_max == None:
126	q_max = self.max
127
128	_qmin_unsmeared, _qmax_unsmeared = self.get_unsmeared_range(q_min,
129	q_max)
130	_first_bin = None
131	_last_bin = None
132
133	#step = (self.max - self.min) / (self.nbins - 1.0)
134	# Find the first and last bin number in all extrapolated and real data
135	try:
136	for i in range(self.nbins):
137	q_i = smearer.get_q(self._smearer, i)
138	if (q_i >= _qmin_unsmeared) and (q_i <= _qmax_unsmeared):
139	# Identify first and last bin
140	if _first_bin is None:
141	_first_bin = i
142	else:
143	_last_bin = i
144	except:
145	msg = "_BaseSmearer.get_bin_range: "
146	msg += " error getting range\n %s" % sys.exc_value
147	raise RuntimeError, msg
148
149	# Find the first and last bin number only in the real data
150	_first_bin, _last_bin = self._get_unextrapolated_bin( \
151	_first_bin, _last_bin)
152
153	return _first_bin, _last_bin
154
155	def __call__(self, iq_in, first_bin = 0, last_bin = None):
156	"""
157	Perform smearing
158	"""
159	# If this is the first time we call for smearing,
160	# initialize the C++ smearer object first
161	if not self._init_complete:
162	self._initialize_smearer()
163
164	if last_bin is None or last_bin >= len(iq_in):
165	last_bin = len(iq_in) - 1
166	# Check that the first bin is positive
167	if first_bin < 0:
168	first_bin = 0
169
170	# With a model given, compute I for the extrapolated points and append
171	# to the iq_in
172	iq_in_temp = iq_in
173	if self.model != None:
174	temp_first, temp_last = self._get_extrapolated_bin( \
175	first_bin, last_bin)
176	if self.nbins_low > 1:
177	iq_in_low = self.model.evalDistribution( \
178	numpy.fabs(self.qvalues[0:self.nbins_low]))
179	iq_in_high = self.model.evalDistribution( \
180	self.qvalues[(len(self.qvalues) - \
181	self.nbins_high - 1):])
182	# Todo: find out who is sending iq[last_poin] = 0.
183	if iq_in[len(iq_in) - 1] == 0:
184	iq_in[len(iq_in) - 1] = iq_in_high[0]
185	# Append the extrapolated points to the data points
186	if self.nbins_low > 0:
187	iq_in_temp = numpy.append(iq_in_low, iq_in)
188	if self.nbins_high > 0:
189	iq_in_temp = numpy.append(iq_in_temp, iq_in_high[1:])
190	else:
191	temp_first = first_bin
192	temp_last = last_bin
193	#iq_in_temp = iq_in
194
195	# Sanity check
196	if len(iq_in_temp) != self.nbins:
197	msg = "Invalid I(q) vector: inconsistent array "
198	msg += " length %d != %s" % (len(iq_in_temp), str(self.nbins))
199	raise RuntimeError, msg
200
201	# Storage for smeared I(q)
202	iq_out = numpy.zeros(self.nbins)
203
204	smear_output = smearer.smear(self._smearer, iq_in_temp, iq_out,
205	#0, self.nbins - 1)
206	temp_first, temp_last)
207	#first_bin, last_bin)
208	if smear_output < 0:
209	msg = "_BaseSmearer: could not smear, code = %g" % smear_output
210	raise RuntimeError, msg
211
212	temp_first += self.nbins_low
213	temp_last = self.nbins - self.nbins_high
214	out = iq_out[temp_first: temp_last]
215
216	return out
217
218	def _initialize_smearer(self):
219	"""
220	"""
221	return NotImplemented
222
223
224	def _get_unextrapolated_bin(self, first_bin = 0, last_bin = 0):
225	"""
226	Get unextrapolated first bin and the last bin
227
228	: param first_bin: extrapolated first_bin
229	: param last_bin: extrapolated last_bin
230
231	: return fist_bin, last_bin: unextrapolated first and last bin
232	"""
233	# For first bin
234	if first_bin <= self.nbins_low:
235	first_bin = 0
236	else:
237	first_bin = first_bin - self.nbins_low
238	# For last bin
239	if last_bin >= (self.nbins - self.nbins_high):
240	last_bin = self.nbins - (self.nbins_high + self.nbins_low + 1)
241	elif last_bin >= self.nbins_low:
242	last_bin = last_bin - self.nbins_low
243	else:
244	last_bin = 0
245	return first_bin, last_bin
246
247	def _get_extrapolated_bin(self, first_bin = 0, last_bin = 0):
248	"""
249	Get extrapolated first bin and the last bin
250
251	: param first_bin: unextrapolated first_bin
252	: param last_bin: unextrapolated last_bin
253
254	: return first_bin, last_bin: extrapolated first and last bin
255	"""
256	# For the first bin
257	# In the case that needs low extrapolation data
258	first_bin = 0
259	# For last bin
260	if last_bin >= self.nbins - (self.nbins_high + self.nbins_low + 1):
261	# In the case that needs higher q extrapolation data
262	last_bin = self.nbins - 1
263	else:
264	# In the case that doesn't need higher q extrapolation data
265	last_bin += self.nbins_low
266
267	return first_bin, last_bin
268
269	class _SlitSmearer(_BaseSmearer):
270	"""
271	Slit smearing for I(q) array
272	"""
273
274	def __init__(self, nbins=None, width=None, height=None, min=None, max=None):
275	"""
276	Initialization
277
278	:param iq: I(q) array [cm-1]
279	:param width: slit width [A-1]
280	:param height: slit height [A-1]
281	:param min: Q_min [A-1]
282	:param max: Q_max [A-1]
283
284	"""
285	_BaseSmearer.__init__(self)
286	## Slit width in Q units
287	self.width = width
288	## Slit height in Q units
289	self.height = height
290	## Q_min (Min Q-value for I(q))
291	self.min = min
292	## Q_max (Max Q_value for I(q))
293	self.max = max
294	## Number of Q bins
295	self.nbins = nbins
296	## Number of points used in the smearing computation
297	self.npts = 3000
298	## Smearing matrix
299	self._weights = None
300	self.qvalues = None
301
302	def _initialize_smearer(self):
303	"""
304	Initialize the C++ smearer object.
305	This method HAS to be called before smearing
306	"""
307	#self._smearer = smearer.new_slit_smearer(self.width,
308	# self.height, self.min, self.max, self.nbins)
309	self._smearer = smearer.new_slit_smearer_with_q(self.width,
310	self.height, self.qvalues)
311	self._init_complete = True
312
313	def get_unsmeared_range(self, q_min, q_max):
314	"""
315	Determine the range needed in unsmeared-Q to cover
316	the smeared Q range
317	"""
318	# Range used for input to smearing
319	_qmin_unsmeared = q_min
320	_qmax_unsmeared = q_max
321	try:
322	_qmin_unsmeared = self.min
323	_qmax_unsmeared = self.max
324	except:
325	logging.error("_SlitSmearer.get_bin_range: %s" % sys.exc_value)
326	return _qmin_unsmeared, _qmax_unsmeared
327
328	class SlitSmearer(_SlitSmearer):
329	"""
330	Adaptor for slit smearing class and SANS data
331	"""
332	def __init__(self, data1D, model = None):
333	"""
334	Assumption: equally spaced bins of increasing q-values.
335
336	:param data1D: data used to set the smearing parameters
337	"""
338	# Initialization from parent class
339	super(SlitSmearer, self).__init__()
340
341	## Slit width
342	self.width = 0
343	self.nbins_low = 0
344	self.nbins_high = 0
345	self.model = model
346	if data1D.dxw is not None and len(data1D.dxw) == len(data1D.x):
347	self.width = data1D.dxw[0]
348	# Sanity check
349	for value in data1D.dxw:
350	if value != self.width:
351	msg = "Slit smearing parameters must "
352	msg += " be the same for all data"
353	raise RuntimeError, msg
354	## Slit height
355	self.height = 0
356	if data1D.dxl is not None and len(data1D.dxl) == len(data1D.x):
357	self.height = data1D.dxl[0]
358	# Sanity check
359	for value in data1D.dxl:
360	if value != self.height:
361	msg = "Slit smearing parameters must be"
362	msg += " the same for all data"
363	raise RuntimeError, msg
364	# If a model is given, get the q extrapolation
365	if self.model == None:
366	data1d_x = data1D.x
367	else:
368	# Take larger sigma
369	if self.height > self.width:
370	# The denominator (2.0) covers all the possible w^2 + h^2 range
371	sigma_in = data1D.dxl / 2.0
372	elif self.width > 0:
373	sigma_in = data1D.dxw / 2.0
374	else:
375	sigma_in = []
376
377	self.nbins_low, self.nbins_high, _, data1d_x = \
378	get_qextrapolate(sigma_in, data1D.x)
379
380	## Number of Q bins
381	self.nbins = len(data1d_x)
382	## Minimum Q
383	self.min = min(data1d_x)
384	## Maximum
385	self.max = max(data1d_x)
386	## Q-values
387	self.qvalues = data1d_x
388
389
390	class _QSmearer(_BaseSmearer):
391	"""
392	Perform Gaussian Q smearing
393	"""
394
395	def __init__(self, nbins=None, width=None, min=None, max=None):
396	"""
397	Initialization
398
399	:param nbins: number of Q bins
400	:param width: array standard deviation in Q [A-1]
401	:param min: Q_min [A-1]
402	:param max: Q_max [A-1]
403	"""
404	_BaseSmearer.__init__(self)
405	## Standard deviation in Q [A-1]
406	self.width = width
407	## Q_min (Min Q-value for I(q))
408	self.min = min
409	## Q_max (Max Q_value for I(q))
410	self.max = max
411	## Number of Q bins
412	self.nbins = nbins
413	## Smearing matrix
414	self._weights = None
415	self.qvalues = None
416
417	def _initialize_smearer(self):
418	"""
419	Initialize the C++ smearer object.
420	This method HAS to be called before smearing
421	"""
422	#self._smearer = smearer.new_q_smearer(numpy.asarray(self.width),
423	# self.min, self.max, self.nbins)
424	self._smearer = smearer.new_q_smearer_with_q(numpy.asarray(self.width),
425	self.qvalues)
426	self._init_complete = True
427
428	def get_unsmeared_range(self, q_min, q_max):
429	"""
430	Determine the range needed in unsmeared-Q to cover
431	the smeared Q range
432	Take 3 sigmas as the offset between smeared and unsmeared space
433	"""
434	# Range used for input to smearing
435	_qmin_unsmeared = q_min
436	_qmax_unsmeared = q_max
437	try:
438	offset = 3.0 * max(self.width)
439	_qmin_unsmeared = max([self.min, q_min - offset])
440	_qmax_unsmeared = min([self.max, q_max + offset])
441	except:
442	logging.error("_QSmearer.get_bin_range: %s" % sys.exc_value)
443	return _qmin_unsmeared, _qmax_unsmeared
444
445
446	class QSmearer(_QSmearer):
447	"""
448	Adaptor for Gaussian Q smearing class and SANS data
449	"""
450	def __init__(self, data1D, model = None):
451	"""
452	Assumption: equally spaced bins of increasing q-values.
453
454	:param data1D: data used to set the smearing parameters
455	"""
456	# Initialization from parent class
457	super(QSmearer, self).__init__()
458	data1d_x = []
459	self.nbins_low = 0
460	self.nbins_high = 0
461	self.model = model
462	## Resolution
463	#self.width = numpy.zeros(len(data1D.x))
464	if data1D.dx is not None and len(data1D.dx) == len(data1D.x):
465	self.width = data1D.dx
466
467	if self.model == None:
468	data1d_x = data1D.x
469	else:
470	self.nbins_low, self.nbins_high, self.width, data1d_x = \
471	get_qextrapolate(self.width, data1D.x)
472
473	## Number of Q bins
474	self.nbins = len(data1d_x)
475	## Minimum Q
476	self.min = min(data1d_x)
477	## Maximum
478	self.max = max(data1d_x)
479	## Q-values
480	self.qvalues = data1d_x
481
482
483	def get_qextrapolate(width, data_x):
484	"""
485	Make fake data_x points extrapolated outside of the data_x points
486
487	: param width: array of std of q resolution
488	: param Data1D.x: Data1D.x array
489
490	: return new_width, data_x_ext: extrapolated width array and x array
491
492	: assumption1: data_x is ordered from lower q to higher q
493	: assumption2: len(data) = len(width)
494	: assumption3: the distance between the data points is more compact
495	than the size of width
496	: Todo1: Make sure that the assumptions are correct for Data1D
497	: Todo2: This fixes the edge problem in Qsmearer but still needs to make
498	smearer interface
499	"""
500	# Length of the width
501	length = len(width)
502	width_low = math.fabs(width[0])
503	width_high = math.fabs(width[length -1])
504
505	# Compare width(dQ) to the data bin size and take smaller one as the bin
506	# size of the extrapolation; this will correct some weird behavior
507	# at the edge: This method was out (commented)
508	# because it becomes very expansive when
509	# bin size is very small comparing to the width.
510	# Now on, we will just give the bin size of the extrapolated points
511	# based on the width.
512	# Find bin sizes
513	#bin_size_low = math.fabs(data_x[1] - data_x[0])
514	#bin_size_high = math.fabs(data_x[length - 1] - data_x[length - 2])
515	# Let's set the bin size 1/3 of the width(sigma), it is good as long as
516	# the scattering is monotonous.
517	#if width_low < (bin_size_low):
518	bin_size_low = width_low / 10.0
519	#if width_high < (bin_size_high):
520	bin_size_high = width_high / 10.0
521
522	# Number of q points required below the 1st data point in order to extend
523	# them 3 times of the width (std)
524	nbins_low = math.ceil(3.0 * width_low / bin_size_low)
525	# Number of q points required above the last data point
526	nbins_high = math.ceil(3.0 * width_high / (bin_size_high))
527	# Make null q points
528	extra_low = numpy.zeros(nbins_low)
529	extra_high = numpy.zeros(nbins_high)
530	# Give extrapolated values
531	ind = 0
532	qvalue = data_x[0] - bin_size_low
533	#if qvalue > 0:
534	while(ind < nbins_low):
535	extra_low[nbins_low - (ind + 1)] = qvalue
536	qvalue -= bin_size_low
537	ind += 1
538	#if qvalue <= 0:
539	# break
540	# Redefine nbins_low
541	nbins_low = ind
542	# Reset ind for another extrapolation
543	ind = 0
544	qvalue = data_x[length -1] + bin_size_high
545	while(ind < nbins_high):
546	extra_high[ind] = qvalue
547	qvalue += bin_size_high
548	ind += 1
549	# Make a new qx array
550	if nbins_low > 0:
551	data_x_ext = numpy.append(extra_low, data_x)
552	else:
553	data_x_ext = data_x
554	data_x_ext = numpy.append(data_x_ext, extra_high)
555
556	# Redefine extra_low and high based on corrected nbins
557	# And note that it is not necessary for extra_width to be a non-zero
558	if nbins_low > 0:
559	extra_low = numpy.zeros(nbins_low)
560	extra_high = numpy.zeros(nbins_high)
561	# Make new width array
562	new_width = numpy.append(extra_low, width)
563	new_width = numpy.append(new_width, extra_high)
564
565	# nbins corrections due to the negative q value
566	nbins_low = nbins_low - len(data_x_ext[data_x_ext<0])
567	return nbins_low, nbins_high, \
568	new_width[data_x_ext>0], data_x_ext[data_x_ext>0]
569
570	if __name__ == '__main__':
571	x = 0.001 * numpy.arange(1, 11)
572	y = 12.0 - numpy.arange(1, 11)
573	print x
574	#for i in range(10): print i, 0.001 + i*0.008/9.0
575	#for i in range(100): print i, int(math.floor( (i/ (100/9.0)) ))
576	s = _SlitSmearer(nbins=10, width=0.0, height=0.005, min=0.001, max=0.010)
577	#s = _QSmearer(nbins=10, width=0.001, min=0.001, max=0.010)
578	s._compute_matrix()
579
580	sy = s(y)
581	print sy
582
583	if True:
584	for i in range(10):
585	print x[i], y[i], sy[i]
586	#print q, ' : ', s.weight(q), s._compute_iq(q)
587	#print q, ' : ', s(q), s._compute_iq(q)
588	#s._compute_iq(q)
589
590
591
592

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

SasView

source: sasview/DataLoader/qsmearing.py @ fc292d3

Download in other formats: