invariant.py @ c5607fa

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

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

SasView

source: sasview/Invariant/invariant.py @ c5607fa

Download in other formats: