invariant.py @ 6ef7ac5a

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

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

SasView

source: sasview/Invariant/invariant.py @ 6ef7ac5a

Download in other formats: