invariant.py @ 68a2c31

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

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

SasView

source: sasview/Invariant/invariant.py @ 68a2c31

Download in other formats: