invariant.py @ 383eeaa

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 383eeaa was 3bb37ef, checked in by Jae Cho <jhjcho@…>, 15 years ago
minor corrections of dx on edges of the Q* integration.
Property mode set to `100644`
File size: 31.4 KB

Line
1	"""
2	This module implements invariant and its related computations.
3	@author: Gervaise B. Alina/UTK
4	"""
5
6	import math
7	import numpy
8
9	from DataLoader.data_info import Data1D as LoaderData1D
10	from DataLoader.qsmearing import smear_selection
11
12
13	# The minimum q-value to be used when extrapolating
14	Q_MINIMUM = 1e-5
15
16	# The maximum q-value to be used when extrapolating
17	Q_MAXIMUM = 10
18
19	# Number of steps in the extrapolation
20	INTEGRATION_NSTEPS = 1000
21
22	def guinier(x, scale=1, radius=0.1):
23	"""
24	Compute a F(x) = scale* e-((radiusx)*2/3).
25	@param x: a vector of q values
26	@param scale: the scale value
27	@param radius: the guinier radius value
28	@return F(x)
29	"""
30	value = numpy.array([math.exp(-((radius * i)**2/3)) for i in x ])
31	return scale * value
32
33	def power_law(x, scale=1, power=4):
34	"""
35	F(x) = scale* (x)^(-power)
36	when power= 4. the model is porod
37	else power_law
38	The model has three parameters:
39	@param x: a vector of q values
40	@param power: power of the function
41	@param scale : scale factor value
42	@param F(x)
43	"""
44	value = numpy.array([ math.pow(i, -power) for i in x ])
45	return scale * value
46
47	class FitFunctor:
48	"""
49	compute f(x)
50	"""
51	def __init__(self, data):
52	"""
53	Determine a and b given a linear equation y = ax + b
54	@param Data: data containing x and y such as y = ax + b
55	"""
56	self.data = data
57	x_len = len(self.data.x) -1
58	#fitting range
59	self.qmin = self.data.x[0]
60	if self.qmin == 0:
61	self.qmin = Q_MINIMUM
62
63	self.qmax = self.data.x[x_len]
64	#Unsmeared q range
65	self._qmin_unsmeared = 0
66	self._qmax_unsmeared = self.data.x[x_len]
67
68	#bin for smear data
69	self._first_unsmeared_bin = 0
70	self._last_unsmeared_bin = x_len
71
72	# Identify the bin range for the unsmeared and smeared spaces
73	self.idx = (self.data.x >= self.qmin) & (self.data.x <= self.qmax)
74	self.idx_unsmeared = (self.data.x >= self._qmin_unsmeared) \
75	& (self.data.x <= self._qmax_unsmeared)
76
77	#get the smear object of data
78	self.smearer = smear_selection( self.data )
79
80
81	def set_fit_range(self ,qmin=None, qmax=None):
82	""" to set the fit range"""
83	if qmin is not None:
84	self.qmin = qmin
85	if qmax is not None:
86	self.qmax = qmax
87
88	# Determine the range needed in unsmeared-Q to cover
89	# the smeared Q range
90	self._qmin_unsmeared = self.qmin
91	self._qmax_unsmeared = self.qmax
92
93	self._first_unsmeared_bin = 0
94	self._last_unsmeared_bin = len(self.data.x)-1
95
96	if self.smearer!=None:
97	self._first_unsmeared_bin, self._last_unsmeared_bin = self.smearer.get_bin_range(self.qmin, self.qmax)
98	self._qmin_unsmeared = self.data.x[self._first_unsmeared_bin]
99	self._qmax_unsmeared = self.data.x[self._last_unsmeared_bin]
100
101	# Identify the bin range for the unsmeared and smeared spaces
102	self.idx = (self.data.x >= self.qmin) & (self.data.x <= self.qmax)
103	self.idx_unsmeared = (self.data.x >= self._qmin_unsmeared) \
104	& (self.data.x <= self._qmax_unsmeared)
105
106	def fit(self):
107	"""
108	Fit data for y = ax + b return a and b
109
110	"""
111	# Compute theory data f(x)
112	fx = numpy.zeros(len(self.data.x))
113	fx = self.data.y[self.idx_unsmeared ]
114	## Smear theory data
115	if self.smearer is not None:
116	fx = self.smearer(fx, self._first_unsmeared_bin,self._last_unsmeared_bin)
117
118	A = numpy.vstack([ self.data.x[self.idx],
119	numpy.ones(len(self.data.x[self.idx]))]).T
120
121	a, b = numpy.linalg.lstsq(A, fx)[0]
122	return a, b
123
124	class InvariantCalculator(object):
125	"""
126	Compute invariant if data is given.
127	Can provide volume fraction and surface area if the user provides
128	Porod constant and contrast values.
129	@precondition: the user must send a data of type DataLoader.Data1D
130	the user provide background and scale values.
131
132	@note: Some computations depends on each others.
133	"""
134
135
136	def __init__(self, data, background=0, scale=1 ):
137	"""
138	Initialize variables
139	@param data: data must be of type DataLoader.Data1D
140	@param contrast: contrast value of type float
141	@param pConst: Porod Constant of type float
142	"""
143	# Background and scale should be private data member if the only way to
144	# change them are by instantiating a new object.
145	self._background = background
146	self._scale = scale
147
148	# The data should be private
149	self._data = self._get_data(data)
150
151	# Since there are multiple variants of Q*, you should force the
152	# user to use the get method and keep Q* a private data member
153	self._qstar = None
154
155	# You should keep the error on Q* so you can reuse it without
156	# recomputing the whole thing.
157	self._qstar_err = 0
158
159	# Extrapolation parameters
160	self._low_extrapolation_npts = 4
161	self._low_extrapolation_function = guinier
162	self._low_extrapolation_power = 4
163
164	self._high_extrapolation_npts = 4
165	self._high_extrapolation_function = power_law
166	self._high_extrapolation_power = 4
167
168
169	def set_extrapolation(self, range, npts=4, function=None, power=4):
170	"""
171	Set the extrapolation parameters for the high or low Q-range.
172	Note that this does not turn extrapolation on or off.
173	@param range: a keyword set the type of extrapolation . type string
174	@param npts: the numbers of q points of data to consider for extrapolation
175	@param function: a keyword to select the function to use for extrapolation.
176	of type string.
177	@param power: an power to apply power_low function
178
179	"""
180	range = range.lower()
181	if range not in ['high', 'low']:
182	raise ValueError, "Extrapolation range should be 'high' or 'low'"
183	function = function.lower()
184	if function not in ['power_law', 'guinier']:
185	raise ValueError, "Extrapolation function should be 'guinier' or 'power_law'"
186
187	if range == 'high':
188	if function != 'power_law':
189	raise ValueError, "Extrapolation only allows a power law at high Q"
190	self._high_extrapolation_npts = npts
191	self._high_extrapolation_power = power
192	else:
193	if function == 'power_law':
194	self._low_extrapolation_function = power_law
195	else:
196	self._low_extrapolation_function = guinier
197	self._low_extrapolation_npts = npts
198	self._low_extrapolation_power = power
199
200	def _get_data(self, data):
201	"""
202	@note this function must be call before computing any type
203	of invariant
204	@return data= self._scale *data - self._background
205	"""
206	if not issubclass(data.__class__, LoaderData1D):
207	#Process only data that inherited from DataLoader.Data_info.Data1D
208	raise ValueError,"Data must be of type DataLoader.Data1D"
209	new_data = self._scale * data - self._background
210	new_data.dxl = data.dxl
211	new_data.dxw = data.dxw
212	return new_data
213
214	def _fit(self, function, qmin=Q_MINIMUM, qmax=Q_MAXIMUM):
215	"""
216	fit data with function using
217	data= self._get_data()
218	fx= Functor(data , function)
219	y = data.y
220	slope, constant = linalg.lstsq(y,fx)
221	@param qmin: data first q value to consider during the fit
222	@param qmax: data last q value to consider during the fit
223	@param function: the function to use during the fit
224	@return a: the scale of the function
225	@return b: the other parameter of the function for guinier will be radius
226	for power_law will be the power value
227	"""
228	if function.__name__ == "guinier":
229	fit_x = numpy.array([x * x for x in self._data.x])
230	qmin = qmin**2
231	qmax = qmax**2
232	fit_y = numpy.array([math.log(y) for y in self._data.y])
233	elif function.__name__ == "power_law":
234	fit_x = numpy.array([math.log(x) for x in self._data.x])
235	qmin = math.log(qmin)
236	qmax = math.log(qmax)
237	fit_y = numpy.array([math.log(y) for y in self._data.y])
238	else:
239	raise ValueError("Unknown function used to fit %s"%function.__name__)
240
241	fit_data = LoaderData1D(x=fit_x, y=fit_y)
242	fit_data.dxl = self._data.dxl
243	fit_data.dxw = self._data.dxw
244
245	functor = FitFunctor(data=fit_data)
246	functor.set_fit_range(qmin=qmin, qmax=qmax)
247	b, a = functor.fit()
248
249	if function.__name__ == "guinier":
250	# b is the radius value of the guinier function
251	b = math.sqrt(-3 * b)
252	if function.__name__ == "power_law":
253	b = -1 * b
254	# a is the scale of the guinier function
255	a = math.exp(a)
256	return a, b
257
258	def _get_qstar(self, data):
259	"""
260	Compute invariant for data
261	@param data: data to use to compute invariant.
262
263	"""
264	if data is None:
265	return 0
266	if data.is_slit_smeared():
267	return self._get_qstar_smear(data)
268	else:
269	return self._get_qstar_unsmear(data)
270
271
272	def get_qstar(self, extrapolation=None):
273	"""
274	Compute the invariant of the local copy of data.
275	Implementation:
276	if slit smear:
277	qstar_0 = self._get_qstar_smear()
278	else:
279	qstar_0 = self._get_qstar_unsmear()
280	if extrapolation is None:
281	return qstar_0
282	if extrapolation==low:
283	return qstar_0 + self._get_qstar_low()
284	elif extrapolation==high:
285	return qstar_0 + self._get_qstar_high()
286	elif extrapolation==both:
287	return qstar_0 + self._get_qstar_low() + self._get_qstar_high()
288
289	@param extrapolation: string to apply optional extrapolation
290	@return q_star: invariant of the data within data's q range
291
292	@warning: When using setting data to Data1D , the user is responsible of
293	checking that the scale and the background are properly apply to the data
294
295	@warning: if error occur self._get_qstar_low() or self._get_qstar_low()
296	their returned value will be ignored
297	"""
298	qstar_0 = self._get_qstar(data=self._data)
299
300	if extrapolation is None:
301	self._qstar = qstar_0
302	return self._qstar
303	# Compute invariant plus invaraint of extrapolated data
304	extrapolation = extrapolation.lower()
305	if extrapolation == "low":
306	self._qstar = qstar_0 + self._get_qstar_low()
307	return self._qstar
308	elif extrapolation == "high":
309	self._qstar = qstar_0 + self._get_qstar_high()
310	return self._qstar
311	elif extrapolation == "both":
312	self._qstar = qstar_0 + self._get_qstar_low() + self._get_qstar_high()
313	return self._qstar
314
315	def _get_qstar_unsmear(self, data):
316	"""
317	Compute invariant for pinhole data.
318	This invariant is given by:
319
320	q_star = x0*2 y0 dx0 +x12 y1 *dx1
321	+ ..+ xn*2 yn *dxn
322
323	where n >= len(data.x)-1
324	dxi = 1/2*(xi+1 - xi) + (xi - xi-1)
325	dx0 = (x1 - x0)/2
326	dxn = (xn - xn-1)/2
327	@param data: the data to use to compute invariant.
328	@return q_star: invariant value for pinhole data. q_star > 0
329	"""
330	if len(data.x) <= 1 or len(data.y) <= 1 or len(data.x)!= len(data.y):
331	msg = "Length x and y must be equal"
332	msg += " and greater than 1; got x=%s, y=%s"%(len(data.x), len(data.y))
333	raise ValueError, msg
334	else:
335	n = len(data.x)- 1
336	#compute the first delta q
337	dx0 = (data.x[1] - data.x[0])/2
338	#compute the last delta q
339	dxn = (data.x[n] - data.x[n-1])/2
340	sum = 0
341	sum += data.x[0] * data.x[0] * data.y[0] * dx0
342	sum += data.x[n] * data.x[n] * data.y[n] * dxn
343
344	if len(data.x) == 2:
345	return sum
346	else:
347	#iterate between for element different from the first and the last
348	for i in xrange(1, n-1):
349	dxi = (data.x[i+1] - data.x[i-1])/2
350	sum += data.x[i] * data.x[i] * data.y[i] * dxi
351	return sum
352
353	def _get_qstar_smear(self, data):
354	"""
355	Compute invariant for slit-smeared data.
356	This invariant is given by:
357	q_star = x0dxl y0dx0 + x1dxl y1 dx1
358	+ ..+ xndxl yn *dxn
359	where n >= len(data.x)-1
360	dxi = 1/2*(xi+1 - xi) + (xi - xi-1)
361	dx0 = (x1 - x0)/2
362	dxn = (xn - xn-1)/2
363	dxl: slit smear value
364
365	@return q_star: invariant value for slit smeared data.
366	"""
367	if not data.is_slit_smeared():
368	msg = "_get_qstar_smear need slit smear data "
369	msg += "Hint :dxl= %s , dxw= %s"%(str(data.dxl), str(data.dxw))
370	raise ValueError, msg
371
372	if len(data.x) <= 1 or len(data.y) <= 1 or len(data.x) != len(data.y)\
373	or len(data.x)!= len(data.dxl):
374	msg = "x, dxl, and y must be have the same length and greater than 1"
375	raise ValueError, msg
376	else:
377	n = len(data.x)-1
378	#compute the first delta
379	dx0 = (data.x[1] - data.x[0])/2
380	#compute the last delta
381	dxn = (data.x[n] - data.x[n-1])/2
382	sum = 0
383	sum += data.x[0] * data.dxl[0] * data.y[0] * dx0
384	sum += data.x[n] * data.dxl[n] * data.y[n] * dxn
385
386	if len(data.x)==2:
387	return sum
388	else:
389	#iterate between for element different from the first and the last
390	for i in xrange(1, n-1):
391	dxi = (data.x[i+1] - data.x[i-1])/2
392	sum += data.x[i] * data.dxl[i] * data.y[i] * dxi
393	return sum
394
395	def _get_qstar_uncertainty(self, data=None):
396	"""
397	Compute uncertainty of invariant value
398	Implementation:
399	if data is None:
400	data = self.data
401
402	if data.slit smear:
403	return self._get_qstar_smear_uncertainty(data)
404	else:
405	return self._get_qstar_unsmear_uncertainty(data)
406
407	@param: data use to compute the invariant which allow uncertainty
408	computation.
409	@return: uncertainty
410	"""
411	if data is None:
412	data = self.data
413
414	if data.is_slit_smeared():
415	return self._get_qstar_smear_uncertainty(data)
416	else:
417	return self._get_qstar_unsmear_uncertainty(data)
418
419	def _get_qstar_unsmear_uncertainty(self, data=None):
420	"""
421	Compute invariant uncertainty with with pinhole data.
422	This uncertainty is given as follow:
423	dq_star = math.sqrt[(x0*2(dy0)dx0)*2 +
424	(x1*2 (dy1)dx1)2 + ..+ (xn2 (dyn)dxn)*2 ]
425	where n >= len(data.x)-1
426	dxi = 1/2*(xi+1 - xi) + (xi - xi-1)
427	dx0 = (x1 - x0)/2
428	dxn = (xn - xn-1)/2
429	dyn: error on dy
430
431	@param data:
432	note: if data doesn't contain dy assume dy= math.sqrt(data.y)
433	"""
434	if data is None:
435	data = self.data
436
437	if len(data.x) <= 1 or len(data.y) <= 1 or \
438	len(self.data.x) != len(self.data.y):
439	msg = "Length of data.x and data.y must be equal"
440	msg += " and greater than 1; got x=%s, y=%s"%(len(data.x),
441	len(data.y))
442	raise ValueError, msg
443	else:
444	#Create error for data without dy error
445	if (data.dy is None) or (not data.dy):
446	dy = math.sqrt(y)
447	else:
448	dy = data.dy
449
450	n = len(data.x) - 1
451	#compute the first delta
452	dx0 = (data.x[1] - data.x[0])/2
453	#compute the last delta
454	dxn= (data.x[n] - data.x[n-1])/2
455	sum = 0
456	sum += (data.x[0] * data.x[0] * dy[0] * dx0)**2
457	sum += (data.x[n] * data.x[n] * dy[n] * dxn)**2
458	if len(data.x) == 2:
459	return math.sqrt(sum)
460	else:
461	#iterate between for element different from the first and the last
462	for i in xrange(1, n-1):
463	dxi = (data.x[i+1] - data.x[i-1])/2
464	sum += (data.x[i] * data.x[i] * dy[i] * dxi)**2
465	return math.sqrt(sum)
466
467	def _get_qstar_smear_uncertainty(self):
468	"""
469	Compute invariant uncertainty with slit smeared data.
470	This uncertainty is given as follow:
471	dq_star = x0dxl dy0 dx0 + x1dxl dy1 dx1
472	+ ..+ xndxl dyn *dxn
473	where n >= len(data.x)-1
474	dxi = 1/2*(xi+1 - xi) + (xi - xi-1)
475	dx0 = (x1 - x0)/2
476	dxn = (xn - xn-1)/2
477	dxl: slit smearing value
478	dyn : error on dy
479	@param data: data of type Data1D where the scale is applied
480	and the background is subtracted.
481
482	note: if data doesn't contain dy assume dy= math.sqrt(data.y)
483	"""
484	if data is None:
485	data = self._data
486
487	if not data.is_slit_smeared():
488	msg = "_get_qstar_smear_uncertainty need slit smear data "
489	msg += "Hint :dxl= %s , dxw= %s"%(str(data.dxl), str(data.dxw))
490	raise ValueError, msg
491
492	if len(data.x) <= 1 or len(data.y) <= 1 or len(data.x) != len(data.y)\
493	or len(data.x) != len(data.dxl):
494	msg = "x, dxl, and y must be have the same length and greater than 1"
495	raise ValueError, msg
496	else:
497	#Create error for data without dy error
498	if (data.dy is None) or (not data.dy):
499	dy = math.sqrt(y)
500	else:
501	dy = data.dy
502
503	n = len(data.x) - 1
504	#compute the first delta
505	dx0 = (data.x[1] - data.x[0])/2
506	#compute the last delta
507	dxn = (data.x[n] - data.x[n-1])/2
508	sum = 0
509	sum += (data.x[0] * data.dxl[0] * dy[0] * dx0)**2
510	sum += (data.x[n] * data.dxl[n] * dy[n] * dxn)**2
511
512	if len(data.x) == 2:
513	return math.sqrt(sum)
514	else:
515	#iterate between for element different from the first and the last
516	for i in xrange(1, n-1):
517	dxi = (data.x[i+1] - data.x[i-1])/2
518	sum += (data.x[i] * data.dxl[i] * dy[i] * dxi)**2
519	return math.sqrt(sum)
520
521	def get_surface(self,contrast, porod_const):
522	"""
523	Compute the surface of the data.
524
525	Implementation:
526	V= self.get_volume_fraction(contrast)
527
528	Compute the surface given by:
529	surface = (2pi V(1- V)*porod_const)/ q_star
530
531	@param contrast: contrast value to compute the volume
532	@param porod_const: Porod constant to compute the surface
533	@return: specific surface
534	"""
535	#Check whether we have Q star
536	if self._qstar is None:
537	self._qstar = self.get_star()
538	if self._qstar == 0:
539	raise RuntimeError("Cannot compute surface, invariant value is zero")
540	# Compute the volume
541	volume = self.get_volume_fraction(contrast)
542	return 2 * math.pi * volume (1 - volume) float(porod_const)/self._qstar
543
544
545	def get_volume_fraction(self, contrast):
546	"""
547	Compute volume fraction is deduced as follow:
548
549	q_star = 2(picontrast)*2 volume( 1- volume)
550	for k = 10^(-8)q_star/(2(pi\|contrast\|)*2)
551	we get 2 values of volume:
552	with 1 - 4 * k >= 0
553	volume1 = (1- sqrt(1- 4*k))/2
554	volume2 = (1+ sqrt(1- 4*k))/2
555
556	q_star: the invariant value included extrapolation is applied
557	unit 1/A^(3)*1/cm
558	q_star = self.get_qstar()
559
560	the result returned will be 0<= volume <= 1
561
562	@param contrast: contrast value provides by the user of type float.
563	contrast unit is 1/A^(2)= 10^(16)cm^(2)
564	@return: volume fraction
565	@note: volume fraction must have no unit
566	"""
567	if contrast < 0:
568	raise ValueError, "contrast must be greater than zero"
569
570	if self._qstar is None:
571	self._qstar = self.get_qstar()
572
573	if self._qstar < 0:
574	raise RuntimeError, "invariant must be greater than zero"
575
576	# Compute intermediate constant
577	k = 1.e-8 * self._qstar/(2 * (math.pi * math.fabs(float(contrast)))**2)
578	#Check discriminant value
579	discrim = 1 - 4 * k
580
581	# Compute volume fraction
582	if discrim < 0:
583	raise RuntimeError, "could not compute the volume fraction: negative discriminant"
584	elif discrim == 0:
585	return 1/2
586	else:
587	volume1 = 0.5 * (1 - math.sqrt(discrim))
588	volume2 = 0.5 * (1 + math.sqrt(discrim))
589
590	if 0 <= volume1 and volume1 <= 1:
591	return volume1
592	elif 0 <= volume2 and volume2 <= 1:
593	return volume2
594	raise RuntimeError, "could not compute the volume fraction: inconsistent results"
595
596	def _get_qstar_low(self):
597	"""
598	Compute the invariant for extrapolated data at low q range.
599
600	Implementation:
601	data = self.get_extra_data_low()
602	return self._get_qstar()
603
604	@return q_star: the invariant for data extrapolated at low q.
605	"""
606	data = self._get_extra_data_low()
607	return self._get_qstar(data=data)
608
609	def _get_qstar_high(self):
610	"""
611	Compute the invariant for extrapolated data at high q range.
612
613	Implementation:
614	data = self.get_extra_data_high()
615	return self._get_qstar()
616
617	@return q_star: the invariant for data extrapolated at high q.
618	"""
619	data = self._get_extra_data_high()
620	return self._get_qstar(data=data)
621
622	def _get_extra_data_low(self):
623	"""
624	This method creates a new data set from the invariant calculator.
625
626	It will use the extrapolation parameters kept as private data members.
627
628	self._low_extrapolation_npts is the number of data points to use in to fit.
629	self._low_extrapolation_function will be used as the fit function.
630
631
632
633	It takes npts first points of data, fits them with a given model
634	then uses the new parameters resulting from the fit to create a new data set.
635
636	The new data first point is Q_MINIMUM.
637
638	The last point of the new data is the first point of the original data.
639	the number of q points of this data is INTEGRATION_NSTEPS.
640
641	@return: a new data of type Data1D
642	"""
643	# Data boundaries for fiiting
644	qmin = self._data.x[0]
645	qmax = self._data.x[self._low_extrapolation_npts - 1]
646
647	try:
648	# fit the data with a model to get the appropriate parameters
649	a, b = self._fit(function=self._low_extrapolation_function,
650	qmin=qmin, qmax=qmax)
651	except:
652	return None
653
654	#q_start point
655	q_start = Q_MINIMUM
656	if Q_MINIMUM >= qmin:
657	q_start = qmin/10
658
659	#create new Data1D to compute the invariant
660	new_x = numpy.linspace(start=q_start,
661	stop=qmin,
662	num=INTEGRATION_NSTEPS,
663	endpoint=True)
664	new_y = self._low_extrapolation_function(x=new_x, scale=a, radius=b)
665	dxl = None
666	dxw = None
667	if self._data.dxl is not None:
668	dxl = numpy.ones(INTEGRATION_NSTEPS)
669	dxl = dxl * self._data.dxl[0]
670	if self._data.dxw is not None:
671	dxw = numpy.ones(INTEGRATION_NSTEPS)
672	dxw = dxw * self._data.dxw[0]
673
674	data_min = LoaderData1D(x=new_x, y=new_y)
675	data_min.dxl = dxl
676	data_min.dxw = dxw
677	self._data.clone_without_data( clone= data_min)
678
679	return data_min
680
681	def _get_extra_data_high(self):
682	"""
683	This method creates a new data from the invariant calculator.
684
685	It takes npts last points of data, fits them with a given model
686	(for this function only power_law will be use), then uses
687	the new parameters resulting from the fit to create a new data set.
688	The first point is the last point of data.
689	The last point of the new data is Q_MAXIMUM.
690	The number of q points of this data is INTEGRATION_NSTEPS.
691
692
693	@return: a new data of type Data1D
694	"""
695	# Data boundaries for fiiting
696	x_len = len(self._data.x) - 1
697	qmin = self._data.x[x_len - (self._high_extrapolation_npts - 1)]
698	qmax = self._data.x[x_len]
699
700	try:
701	# fit the data with a model to get the appropriate parameters
702	a, b = self._fit(function=self._high_extrapolation_function,
703	qmin=qmin, qmax=qmax)
704	except:
705	return None
706
707	#create new Data1D to compute the invariant
708	new_x = numpy.linspace(start=qmax,
709	stop=Q_MAXIMUM,
710	num=INTEGRATION_NSTEPS,
711	endpoint=True)
712
713	new_y = self._high_extrapolation_function(x=new_x, scale=a, power=b)
714
715	dxl = None
716	dxw = None
717	if self._data.dxl is not None:
718	dxl = numpy.ones(INTEGRATION_NSTEPS)
719	dxl = dxl * self._data.dxl[0]
720	if self._data.dxw is not None:
721	dxw = numpy.ones(INTEGRATION_NSTEPS)
722	dxw = dxw * self._data.dxw[0]
723
724	data_max = LoaderData1D(x=new_x, y=new_y)
725	data_max.dxl = dxl
726	data_max.dxw = dxw
727	self._data.clone_without_data(clone=data_max)
728
729	return data_max
730
731	def get_qstar_with_error(self, extrapolation=None):
732	"""
733	Compute the invariant uncertainty.
734	This uncertainty computation depends on whether or not the data is
735	smeared.
736	@return: invariant, the invariant uncertainty
737	return self._get_qstar(), self._get_qstar_smear_uncertainty()
738	"""
739	if self._qstar is None:
740	self._qstar = self.get_qstar(extrapolation=extrapolation)
741	if self._qstar_err is None:
742	self._qstar_err = self._get_qstar_smear_uncertainty()
743
744	return self._qstar, self._qstar_err
745
746	def get_volume_fraction_with_error(self, contrast):
747	"""
748	Compute uncertainty on volume value as well as the volume fraction
749	This uncertainty is given by the following equation:
750	dV = 0.5 * (4k dq_star) /(2* math.sqrt(1-k* q_star))
751
752	for k = 10^(-8)q_star/(2(pi\|contrast\|)*2)
753
754	q_star: the invariant value including extrapolated value if existing
755	dq_star: the invariant uncertainty
756	dV: the volume uncertainty
757	@param contrast: contrast value
758	@return: V, dV = self.get_volume_fraction_with_error(contrast), dV
759	"""
760	self._qstar, self._qstar_err = self.get_qstar_with_error()
761
762	volume = self.get_volume_fraction(contrast)
763	if self._qstar < 0:
764	raise ValueError, "invariant must be greater than zero"
765
766	k = 1.e-8 * self._qstar /(2 * (math.pi* math.fabs(float(contrast)))**2)
767	#check value inside the sqrt function
768	value = 1 - k * self._qstar
769	if (value) <= 0:
770	raise ValueError, "Cannot compute incertainty on volume"
771	# Compute uncertainty
772	uncertainty = (0.5 * 4 * k * self._qstar_err)/(2 * math.sqrt(1 - k * self._qstar))
773
774	return volume, math.fabs(uncertainty)
775
776	def get_surface_with_error(self, contrast, porod_const):
777	"""
778	Compute uncertainty of the surface value as well as thesurface value
779	this uncertainty is given as follow:
780
781	dS = porod_const 2pi[( dV -2VdV)/q_star
782	+ dq_star(v-v**2)
783
784	q_star: the invariant value including extrapolated value if existing
785	dq_star: the invariant uncertainty
786	V: the volume fraction value
787	dV: the volume uncertainty
788
789	@param contrast: contrast value
790	@param porod_const: porod constant value
791	@return S, dS: the surface, with its uncertainty
792	"""
793	v, dv = self.get_volume_fraction_with_error(contrast)
794	self._qstar, self._qstar_err = self.get_qstar_with_error()
795	if self._qstar <= 0:
796	raise ValueError, "invariant must be greater than zero"
797	ds = porod_const * 2 * math.pi * (( dv - 2 * v * dv)/ self._qstar\
798	+ self._qstar_err * ( v - v**2))
799	s = self.get_surface(contrast=contrast, porod_const=porod_const)
800	return s, ds

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

SasView

source: sasview/Invariant/invariant.py @ 383eeaa

Download in other formats: