qsmearing.py @ 7241d56

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 7241d56 was 7241d56, checked in by Jae Cho <jhjcho@…>, 13 years ago
improved edge behavior of Qsmearer
Property mode set to `100644`
File size: 19.7 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)
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	iq_in_low = self.model.evalDistribution( \
177	numpy.fabs(self.qvalues[0:self.nbins_low]))
178	iq_in_high = self.model.evalDistribution( \
179	self.qvalues[(len(self.qvalues) - \
180	self.nbins_high - 1):])
181	# Todo: find out who is sending iq[last_poin] = 0.
182	if iq_in[len(iq_in) - 1] == 0:
183	iq_in[len(iq_in) - 1] = iq_in_high[0]
184	# Append the extrapolated points to the data points
185	if self.nbins_low > 0:
186	iq_in_temp = numpy.append(iq_in_low, iq_in)
187	if self.nbins_high > 0:
188	iq_in_temp = numpy.append(iq_in_temp, iq_in_high[1:])
189	else:
190	temp_first = first_bin
191	temp_last = last_bin
192	iq_in_temp = iq_in
193
194	# Sanity check
195	if len(iq_in_temp) != self.nbins:
196	msg = "Invalid I(q) vector: inconsistent array "
197	msg += " length %d != %s" % (len(iq_in_temp), str(self.nbins))
198	raise RuntimeError, msg
199
200	# Storage for smeared I(q)
201	iq_out = numpy.zeros(self.nbins)
202
203	smear_output = smearer.smear(self._smearer, iq_in_temp, iq_out,
204	#0, self.nbins - 1)
205	temp_first, temp_last)
206	#first_bin, last_bin)
207	if smear_output < 0:
208	msg = "_BaseSmearer: could not smear, code = %g" % smear_output
209	raise RuntimeError, msg
210
211	temp_first += self.nbins_low
212	temp_last = self.nbins - self.nbins_high
213
214	return iq_out[temp_first: temp_last]
215
216	def _initialize_smearer(self):
217	"""
218	"""
219	return NotImplemented
220
221	def set_model(self, model):
222	"""
223	Set model
224	"""
225	if model != None:
226	self.model = model
227
228
229	def _get_unextrapolated_bin(self, first_bin = 0, last_bin = 0):
230	"""
231	Get unextrapolated first bin and the last bin
232
233	: param first_bin: extrapolated first_bin
234	: param last_bin: extrapolated last_bin
235
236	: return fist_bin, last_bin: unextrapolated first and last bin
237	"""
238	# For first bin
239	if first_bin <= self.nbins_low:
240	first_bin = 0
241	else:
242	first_bin = first_bin - self.nbins_low
243	# For last bin
244	if last_bin >= (self.nbins - self.nbins_high):
245	last_bin = self.nbins - (self.nbins_high + self.nbins_low + 1)
246	elif last_bin >= self.nbins_low:
247	last_bin = last_bin - self.nbins_low
248	else:
249	last_bin = 0
250	return first_bin, last_bin
251
252	def _get_extrapolated_bin(self, first_bin = 0, last_bin = 0):
253	"""
254	Get extrapolated first bin and the last bin
255
256	: param first_bin: unextrapolated first_bin
257	: param last_bin: unextrapolated last_bin
258
259	: return first_bin, last_bin: extrapolated first and last bin
260	"""
261	# For the first bin
262	# In the case that needs low extrapolation data
263	first_bin = 0
264	# For last bin
265	if last_bin >= self.nbins - (self.nbins_high + self.nbins_low + 1):
266	# In the case that needs higher q extrapolation data
267	last_bin = self.nbins - 1
268	else:
269	# In the case that doesn't need higher q extrapolation data
270	last_bin += self.nbins_low
271
272	return first_bin, last_bin
273
274	class _SlitSmearer(_BaseSmearer):
275	"""
276	Slit smearing for I(q) array
277	"""
278
279	def __init__(self, nbins=None, width=None, height=None, min=None, max=None):
280	"""
281	Initialization
282
283	:param iq: I(q) array [cm-1]
284	:param width: slit width [A-1]
285	:param height: slit height [A-1]
286	:param min: Q_min [A-1]
287	:param max: Q_max [A-1]
288
289	"""
290	_BaseSmearer.__init__(self)
291	## Slit width in Q units
292	self.width = width
293	## Slit height in Q units
294	self.height = height
295	## Q_min (Min Q-value for I(q))
296	self.min = min
297	## Q_max (Max Q_value for I(q))
298	self.max = max
299	## Number of Q bins
300	self.nbins = nbins
301	## Number of points used in the smearing computation
302	self.npts = 1000
303	## Smearing matrix
304	self._weights = None
305	self.qvalues = None
306
307	def _initialize_smearer(self):
308	"""
309	Initialize the C++ smearer object.
310	This method HAS to be called before smearing
311	"""
312	#self._smearer = smearer.new_slit_smearer(self.width,
313	# self.height, self.min, self.max, self.nbins)
314	self._smearer = smearer.new_slit_smearer_with_q(self.width,
315	self.height, self.qvalues)
316	self._init_complete = True
317
318	def get_unsmeared_range(self, q_min, q_max):
319	"""
320	Determine the range needed in unsmeared-Q to cover
321	the smeared Q range
322	"""
323	# Range used for input to smearing
324	_qmin_unsmeared = q_min
325	_qmax_unsmeared = q_max
326	try:
327	_qmin_unsmeared = self.min
328	_qmax_unsmeared = self.max
329	except:
330	logging.error("_SlitSmearer.get_bin_range: %s" % sys.exc_value)
331	return _qmin_unsmeared, _qmax_unsmeared
332
333	class SlitSmearer(_SlitSmearer):
334	"""
335	Adaptor for slit smearing class and SANS data
336	"""
337	def __init__(self, data1D):
338	"""
339	Assumption: equally spaced bins of increasing q-values.
340
341	:param data1D: data used to set the smearing parameters
342	"""
343	# Initialization from parent class
344	super(SlitSmearer, self).__init__()
345
346	## Slit width
347	self.width = 0
348	self.nbins_low = 0
349	self.nbins_high = 0
350	if data1D.dxw is not None and len(data1D.dxw) == len(data1D.x):
351	self.width = data1D.dxw[0]
352	# Sanity check
353	for value in data1D.dxw:
354	if value != self.width:
355	msg = "Slit smearing parameters must "
356	msg += " be the same for all data"
357	raise RuntimeError, msg
358	## Slit height
359	self.height = 0
360	if data1D.dxl is not None and len(data1D.dxl) == len(data1D.x):
361	self.height = data1D.dxl[0]
362	# Sanity check
363	for value in data1D.dxl:
364	if value != self.height:
365	msg = "Slit smearing parameters must be"
366	msg += " the same for all data"
367	raise RuntimeError, msg
368
369	## Number of Q bins
370	self.nbins = len(data1D.x)
371	## Minimum Q
372	self.min = min(data1D.x)
373	## Maximum
374	self.max = max(data1D.x)
375	## Q-values
376	self.qvalues = data1D.x
377
378
379	class _QSmearer(_BaseSmearer):
380	"""
381	Perform Gaussian Q smearing
382	"""
383
384	def __init__(self, nbins=None, width=None, min=None, max=None):
385	"""
386	Initialization
387
388	:param nbins: number of Q bins
389	:param width: array standard deviation in Q [A-1]
390	:param min: Q_min [A-1]
391	:param max: Q_max [A-1]
392	"""
393	_BaseSmearer.__init__(self)
394	## Standard deviation in Q [A-1]
395	self.width = width
396	## Q_min (Min Q-value for I(q))
397	self.min = min
398	## Q_max (Max Q_value for I(q))
399	self.max = max
400	## Number of Q bins
401	self.nbins = nbins
402	## Smearing matrix
403	self._weights = None
404	self.qvalues = None
405
406	def _initialize_smearer(self):
407	"""
408	Initialize the C++ smearer object.
409	This method HAS to be called before smearing
410	"""
411	#self._smearer = smearer.new_q_smearer(numpy.asarray(self.width),
412	# self.min, self.max, self.nbins)
413	self._smearer = smearer.new_q_smearer_with_q(numpy.asarray(self.width),
414	self.qvalues)
415	self._init_complete = True
416
417	def get_unsmeared_range(self, q_min, q_max):
418	"""
419	Determine the range needed in unsmeared-Q to cover
420	the smeared Q range
421	Take 3 sigmas as the offset between smeared and unsmeared space
422	"""
423	# Range used for input to smearing
424	_qmin_unsmeared = q_min
425	_qmax_unsmeared = q_max
426	try:
427	offset = 3.0 * max(self.width)
428	_qmin_unsmeared = max([self.min, q_min - offset])
429	_qmax_unsmeared = min([self.max, q_max + offset])
430	except:
431	logging.error("_QSmearer.get_bin_range: %s" % sys.exc_value)
432	return _qmin_unsmeared, _qmax_unsmeared
433
434
435	class QSmearer(_QSmearer):
436	"""
437	Adaptor for Gaussian Q smearing class and SANS data
438	"""
439	def __init__(self, data1D, model = None):
440	"""
441	Assumption: equally spaced bins of increasing q-values.
442
443	:param data1D: data used to set the smearing parameters
444	"""
445	# Initialization from parent class
446	super(QSmearer, self).__init__()
447	data1d_x = []
448	self.nbins_low = 0
449	self.nbins_high = 0
450	self.model = model
451	## Resolution
452	#self.width = numpy.zeros(len(data1D.x))
453	if data1D.dx is not None and len(data1D.dx) == len(data1D.x):
454	self.width = data1D.dx
455
456	if self.model == None:
457	data1d_x = data1D.x
458	else:
459	self.nbins_low, self.nbins_high, self.width, data1d_x = \
460	get_qextrapolate(self.width, data1D.x)
461
462	## Number of Q bins
463	self.nbins = len(data1d_x)
464	## Minimum Q
465	self.min = min(data1d_x)
466	## Maximum
467	self.max = max(data1d_x)
468	## Q-values
469	self.qvalues = data1d_x
470
471
472	def get_qextrapolate(width, data_x):
473	"""
474	Make fake data_x points extrapolated outside of the data_x points
475
476	: param width: array of std of q resolution
477	: param Data1D.x: Data1D.x array
478
479	: return new_width, data_x_ext: extrapolated width array and x array
480
481	: assumption1: data_x is ordered from lower q to higher q
482	: assumption2: len(data) = len(width)
483	: assumption3: the distance between the data points is more compact
484	than the size of width
485	: Todo1: Make sure that the assumptions are correct for Data1D
486	: Todo2: This fixes the edge problem in Qsmearer but still needs to make
487	smearer interface
488	"""
489	# Length of the width
490	length = len(width)
491	width_low = math.fabs(width[0])
492	width_high = math.fabs(width[length -1])
493	# Find bin sizes
494	bin_size_low = math.fabs(data_x[1] - data_x[0])
495	bin_size_high = math.fabs(data_x[length - 1] - data_x[length - 2])
496	# Compare width(dQ) to the data bin size and take smaller one as the bin
497	# size of the extrapolation; this will correct some weird behavior
498	# at the edge
499	if width_low < (bin_size_low):
500	bin_size_low = width_low / 2.0
501	if width_high < (bin_size_high):
502	bin_size_high = width_high / 2.0
503
504	# Number of q points required below the 1st data point in order to extend
505	# them 3 times of the width (std)
506	nbins_low = math.ceil(3 * width_low / bin_size_low)
507	# Number of q points required above the last data point
508	nbins_high = math.ceil(3 * width_high / (bin_size_high))
509	# Make null q points
510	extra_low = numpy.zeros(nbins_low)
511	extra_high = numpy.zeros(nbins_high)
512	# Give extrapolated values
513	ind = 0
514	qvalue = data_x[0] - bin_size_low
515	while(ind < nbins_low):
516	extra_low[nbins_low - (ind + 1)] = qvalue
517	qvalue -= bin_size_low
518	ind += 1
519	# Remove the points <= 0
520	#extra_low = extra_low[extra_low > 0]
521	#nbins_low = len(extra_low)
522	# Reset ind for another extrapolation
523	ind = 0
524	qvalue = data_x[length -1] + bin_size_high
525	while(ind < nbins_high):
526	extra_high[ind] = qvalue
527	qvalue += bin_size_high
528	ind += 1
529	# Make a new qx array
530	data_x_ext = numpy.append(extra_low, data_x)
531	data_x_ext = numpy.append(data_x_ext, extra_high)
532
533	# Redefine extra_low and high based on corrected nbins
534	# And note that it is not necessary for extra_width to be a non-zero
535	extra_low = numpy.zeros(nbins_low)
536	extra_high = numpy.zeros(nbins_high)
537	# Make new width array
538	new_width = numpy.append(extra_low, width)
539	new_width = numpy.append(new_width, extra_high)
540
541	return nbins_low, nbins_high, new_width, data_x_ext
542
543	if __name__ == '__main__':
544	x = 0.001 * numpy.arange(1, 11)
545	y = 12.0 - numpy.arange(1, 11)
546	print x
547	#for i in range(10): print i, 0.001 + i*0.008/9.0
548	#for i in range(100): print i, int(math.floor( (i/ (100/9.0)) ))
549	s = _SlitSmearer(nbins=10, width=0.0, height=0.005, min=0.001, max=0.010)
550	#s = _QSmearer(nbins=10, width=0.001, min=0.001, max=0.010)
551	s._compute_matrix()
552
553	sy = s(y)
554	print sy
555
556	if True:
557	for i in range(10):
558	print x[i], y[i], sy[i]
559	#print q, ' : ', s.weight(q), s._compute_iq(q)
560	#print q, ' : ', s(q), s._compute_iq(q)
561	#s._compute_iq(q)
562
563
564
565

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

SasView

source: sasview/DataLoader/qsmearing.py @ 7241d56

Download in other formats: