qsmearing.py @ 16d7079

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

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

SasView

source: sasview/DataLoader/qsmearing.py @ 16d7079

Download in other formats: