invariant.py @ c451be9

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 c451be9 was 4d270706, checked in by Gervaise Alina <gervyh@…>, 15 years ago
change on method name and fit
Property mode set to `100644`
File size: 34.0 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 sans.guiframe.dataFitting 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	fit_x = numpy.array([x * x for x in self._data.x])
236	qmin = qmin**2
237	qmax = qmax**2
238	fit_y = numpy.array([math.log(y) for y in self._data.y])
239	fit_dy = numpy.array([y for y in self._data.y])
240	fit_dy = numpy.array([dy for dy in self._data.dy])/fit_dy
241
242	elif function.__name__ == "power_law":
243	fit_x = numpy.array([math.log(x) for x in self._data.x])
244	qmin = math.log(qmin)
245	qmax = math.log(qmax)
246	fit_y = numpy.array([math.log(y) for y in self._data.y])
247	fit_dy = numpy.array([y for y in self._data.y])
248	fit_dy = numpy.array([dy for dy in self._data.dy])/fit_dy
249
250	else:
251	raise ValueError("Unknown function used to fit %s"%function.__name__)
252
253
254	#else:
255	fit_data = LoaderData1D(x=fit_x, y=fit_y, dy=fit_dy)
256	fit_data.dxl = self._data.dxl
257	fit_data.dxw = self._data.dxw
258	functor = FitFunctor(data=fit_data)
259	functor.set_fit_range(qmin=qmin, qmax=qmax)
260	b, a = functor.fit(power=power)
261
262
263	if function.__name__ == "guinier":
264	# b is the radius value of the guinier function
265	if b>=0:
266	raise ValueError("Guinier fit was not converged")
267	else:
268	b = math.sqrt(-3 * b)
269
270
271	if function.__name__ == "power_law":
272	b = -1 * b
273	if b <= 0:
274	raise ValueError("Power_law fit expected posive power, but got %s"%power)
275	# a is the scale of the guinier function
276	a = math.exp(a)
277
278	return a, b
279
280	def _get_qstar(self, data):
281	"""
282	Compute invariant for data
283	@param data: data to use to compute invariant.
284
285	"""
286	if data is None:
287	return 0
288	if data.is_slit_smeared():
289	return self._get_qstar_smear(data)
290	else:
291	return self._get_qstar_unsmear(data)
292
293	def _get_qstar_unsmear(self, data):
294	"""
295	Compute invariant for pinhole data.
296	This invariant is given by:
297
298	q_star = x0*2 y0 dx0 +x12 y1 *dx1
299	+ ..+ xn*2 yn *dxn
300
301	where n >= len(data.x)-1
302	dxi = 1/2*(xi+1 - xi) + (xi - xi-1)
303	dx0 = (x1 - x0)/2
304	dxn = (xn - xn-1)/2
305	@param data: the data to use to compute invariant.
306	@return q_star: invariant value for pinhole data. q_star > 0
307	"""
308	if len(data.x) <= 1 or len(data.y) <= 1 or len(data.x)!= len(data.y):
309	msg = "Length x and y must be equal"
310	msg += " and greater than 1; got x=%s, y=%s"%(len(data.x), len(data.y))
311	raise ValueError, msg
312	else:
313	n = len(data.x)- 1
314	#compute the first delta q
315	dx0 = (data.x[1] - data.x[0])/2
316	#compute the last delta q
317	dxn = (data.x[n] - data.x[n-1])/2
318	sum = 0
319	sum += data.x[0] * data.x[0] * data.y[0] * dx0
320	sum += data.x[n] * data.x[n] * data.y[n] * dxn
321
322	if len(data.x) == 2:
323	return sum
324	else:
325	#iterate between for element different from the first and the last
326	for i in xrange(1, n-1):
327	dxi = (data.x[i+1] - data.x[i-1])/2
328	sum += data.x[i] * data.x[i] * data.y[i] * dxi
329	return sum
330
331	def _get_qstar_smear(self, data):
332	"""
333	Compute invariant for slit-smeared data.
334	This invariant is given by:
335	q_star = x0dxl y0dx0 + x1dxl y1 dx1
336	+ ..+ xndxl yn *dxn
337	where n >= len(data.x)-1
338	dxi = 1/2*(xi+1 - xi) + (xi - xi-1)
339	dx0 = (x1 - x0)/2
340	dxn = (xn - xn-1)/2
341	dxl: slit smear value
342
343	@return q_star: invariant value for slit smeared data.
344	"""
345	if not data.is_slit_smeared():
346	msg = "_get_qstar_smear need slit smear data "
347	msg += "Hint :dxl= %s , dxw= %s"%(str(data.dxl), str(data.dxw))
348	raise ValueError, msg
349
350	if len(data.x) <= 1 or len(data.y) <= 1 or len(data.x) != len(data.y)\
351	or len(data.x)!= len(data.dxl):
352	msg = "x, dxl, and y must be have the same length and greater than 1"
353	raise ValueError, msg
354	else:
355	n = len(data.x)-1
356	#compute the first delta
357	dx0 = (data.x[1] - data.x[0])/2
358	#compute the last delta
359	dxn = (data.x[n] - data.x[n-1])/2
360	sum = 0
361	sum += data.x[0] * data.dxl[0] * data.y[0] * dx0
362	sum += data.x[n] * data.dxl[n] * data.y[n] * dxn
363
364	if len(data.x)==2:
365	return sum
366	else:
367	#iterate between for element different from the first and the last
368	for i in xrange(1, n-1):
369	dxi = (data.x[i+1] - data.x[i-1])/2
370	sum += data.x[i] * data.dxl[i] * data.y[i] * dxi
371	return sum
372
373	def _get_qstar_uncertainty(self, data=None):
374	"""
375	Compute uncertainty of invariant value
376	Implementation:
377	if data is None:
378	data = self.data
379
380	if data.slit smear:
381	return self._get_qstar_smear_uncertainty(data)
382	else:
383	return self._get_qstar_unsmear_uncertainty(data)
384
385	@param: data use to compute the invariant which allow uncertainty
386	computation.
387	@return: uncertainty
388	"""
389	if data is None:
390	data = self.data
391
392	if data.is_slit_smeared():
393	return self._get_qstar_smear_uncertainty(data)
394	else:
395	return self._get_qstar_unsmear_uncertainty(data)
396
397	def _get_qstar_unsmear_uncertainty(self, data=None):
398	"""
399	Compute invariant uncertainty with with pinhole data.
400	This uncertainty is given as follow:
401	dq_star = math.sqrt[(x0*2(dy0)dx0)*2 +
402	(x1*2 (dy1)dx1)2 + ..+ (xn2 (dyn)dxn)*2 ]
403	where n >= len(data.x)-1
404	dxi = 1/2*(xi+1 - xi) + (xi - xi-1)
405	dx0 = (x1 - x0)/2
406	dxn = (xn - xn-1)/2
407	dyn: error on dy
408
409	@param data:
410	note: if data doesn't contain dy assume dy= math.sqrt(data.y)
411	"""
412	if data is None:
413	data = self.data
414
415	if len(data.x) <= 1 or len(data.y) <= 1 or \
416	len(self.data.x) != len(self.data.y):
417	msg = "Length of data.x and data.y must be equal"
418	msg += " and greater than 1; got x=%s, y=%s"%(len(data.x),
419	len(data.y))
420	raise ValueError, msg
421	else:
422	#Create error for data without dy error
423	if (data.dy is None) or (not data.dy):
424	dy = math.sqrt(y)
425	else:
426	dy = data.dy
427
428	n = len(data.x) - 1
429	#compute the first delta
430	dx0 = (data.x[1] - data.x[0])/2
431	#compute the last delta
432	dxn= (data.x[n] - data.x[n-1])/2
433	sum = 0
434	sum += (data.x[0] * data.x[0] * dy[0] * dx0)**2
435	sum += (data.x[n] * data.x[n] * dy[n] * dxn)**2
436	if len(data.x) == 2:
437	return math.sqrt(sum)
438	else:
439	#iterate between for element different from the first and the last
440	for i in xrange(1, n-1):
441	dxi = (data.x[i+1] - data.x[i-1])/2
442	sum += (data.x[i] * data.x[i] * dy[i] * dxi)**2
443	return math.sqrt(sum)
444
445	def _get_qstar_smear_uncertainty(self):
446	"""
447	Compute invariant uncertainty with slit smeared data.
448	This uncertainty is given as follow:
449	dq_star = x0dxl dy0 dx0 + x1dxl dy1 dx1
450	+ ..+ xndxl dyn *dxn
451	where n >= len(data.x)-1
452	dxi = 1/2*(xi+1 - xi) + (xi - xi-1)
453	dx0 = (x1 - x0)/2
454	dxn = (xn - xn-1)/2
455	dxl: slit smearing value
456	dyn : error on dy
457	@param data: data of type Data1D where the scale is applied
458	and the background is subtracted.
459
460	note: if data doesn't contain dy assume dy= math.sqrt(data.y)
461	"""
462	if data is None:
463	data = self._data
464
465	if not data.is_slit_smeared():
466	msg = "_get_qstar_smear_uncertainty need slit smear data "
467	msg += "Hint :dxl= %s , dxw= %s"%(str(data.dxl), str(data.dxw))
468	raise ValueError, msg
469
470	if len(data.x) <= 1 or len(data.y) <= 1 or len(data.x) != len(data.y)\
471	or len(data.x) != len(data.dxl):
472	msg = "x, dxl, and y must be have the same length and greater than 1"
473	raise ValueError, msg
474	else:
475	#Create error for data without dy error
476	if (data.dy is None) or (not data.dy):
477	dy = math.sqrt(y)
478	else:
479	dy = data.dy
480
481	n = len(data.x) - 1
482	#compute the first delta
483	dx0 = (data.x[1] - data.x[0])/2
484	#compute the last delta
485	dxn = (data.x[n] - data.x[n-1])/2
486	sum = 0
487	sum += (data.x[0] * data.dxl[0] * dy[0] * dx0)**2
488	sum += (data.x[n] * data.dxl[n] * dy[n] * dxn)**2
489
490	if len(data.x) == 2:
491	return math.sqrt(sum)
492	else:
493	#iterate between for element different from the first and the last
494	for i in xrange(1, n-1):
495	dxi = (data.x[i+1] - data.x[i-1])/2
496	sum += (data.x[i] * data.dxl[i] * dy[i] * dxi)**2
497	return math.sqrt(sum)
498
499	def get_qstar_low(self):
500	"""
501	Compute the invariant for extrapolated data at low q range.
502
503	Implementation:
504	data = self.get_extra_data_low()
505	return self._get_qstar()
506
507	@return q_star: the invariant for data extrapolated at low q.
508	"""
509	data = self.get_extra_data_low()
510	return self._get_qstar(data=data)
511
512	def get_qstar_high(self):
513	"""
514	Compute the invariant for extrapolated data at high q range.
515
516	Implementation:
517	data = self.get_extra_data_high()
518	return self._get_qstar()
519
520	@return q_star: the invariant for data extrapolated at high q.
521	"""
522	data = self.get_extra_data_high()
523	return self._get_qstar(data=data)
524
525	def get_extra_data_low(self):
526	"""
527	This method creates a new data set from the invariant calculator.
528
529	It will use the extrapolation parameters kept as private data members.
530
531	self._low_extrapolation_npts is the number of data points to use in to fit.
532	self._low_extrapolation_function will be used as the fit function.
533
534
535
536	It takes npts first points of data, fits them with a given model
537	then uses the new parameters resulting from the fit to create a new data set.
538
539	The new data first point is Q_MINIMUM.
540
541	The last point of the new data is the first point of the original data.
542	the number of q points of this data is INTEGRATION_NSTEPS.
543
544	@return: a new data of type Data1D
545	"""
546
547	# Data boundaries for fiiting
548	qmin = self._data.x[0]
549	qmax = self._data.x[self._low_extrapolation_npts - 1]
550	try:
551	# fit the data with a model to get the appropriate parameters
552	a, b = self._fit(function=self._low_extrapolation_function,
553	qmin=qmin,
554	qmax=qmax,
555	power=self._low_extrapolation_power)
556	except:
557	#raise
558	return None
559
560	#q_start point
561	q_start = Q_MINIMUM
562	if Q_MINIMUM >= qmin:
563	q_start = qmin/10
564
565	#create new Data1D to compute the invariant
566	new_x = numpy.linspace(start=q_start,
567	stop=qmin,
568	num=INTEGRATION_NSTEPS,
569	endpoint=True)
570	new_y = self._low_extrapolation_function(x=new_x, scale=a, radius=b)
571	dxl = None
572	dxw = None
573	if self._data.dxl is not None:
574	dxl = numpy.ones(INTEGRATION_NSTEPS)
575	dxl = dxl * self._data.dxl[0]
576	if self._data.dxw is not None:
577	dxw = numpy.ones(INTEGRATION_NSTEPS)
578	dxw = dxw * self._data.dxw[0]
579
580	data_min = LoaderData1D(x=new_x, y=new_y)
581	data_min.dxl = dxl
582	data_min.dxw = dxw
583	self._data.clone_without_data( clone= data_min)
584
585	return data_min
586
587	def get_extra_data_high(self):
588	"""
589	This method creates a new data from the invariant calculator.
590
591	It takes npts last points of data, fits them with a given model
592	(for this function only power_law will be use), then uses
593	the new parameters resulting from the fit to create a new data set.
594	The first point is the last point of data.
595	The last point of the new data is Q_MAXIMUM.
596	The number of q points of this data is INTEGRATION_NSTEPS.
597
598
599	@return: a new data of type Data1D
600	"""
601	# Data boundaries for fiiting
602	x_len = len(self._data.x) - 1
603	qmin = self._data.x[x_len - (self._high_extrapolation_npts - 1)]
604	qmax = self._data.x[x_len]
605
606	try:
607	# fit the data with a model to get the appropriate parameters
608	a, b = self._fit(function=self._high_extrapolation_function,
609	qmin=qmin,
610	qmax=qmax,
611	power=self._high_extrapolation_power)
612	except:
613	return None
614
615	#create new Data1D to compute the invariant
616	new_x = numpy.linspace(start=qmax,
617	stop=Q_MAXIMUM,
618	num=INTEGRATION_NSTEPS,
619	endpoint=True)
620
621	new_y = self._high_extrapolation_function(x=new_x, scale=a, power=b)
622
623	dxl = None
624	dxw = None
625	if self._data.dxl is not None:
626	dxl = numpy.ones(INTEGRATION_NSTEPS)
627	dxl = dxl * self._data.dxl[0]
628	if self._data.dxw is not None:
629	dxw = numpy.ones(INTEGRATION_NSTEPS)
630	dxw = dxw * self._data.dxw[0]
631
632	data_max = LoaderData1D(x=new_x, y=new_y)
633	data_max.dxl = dxl
634	data_max.dxw = dxw
635	self._data.clone_without_data(clone=data_max)
636
637	return data_max
638
639	def set_extrapolation(self, range, npts=4, function=None, power=None):
640	"""
641	Set the extrapolation parameters for the high or low Q-range.
642	Note that this does not turn extrapolation on or off.
643	@param range: a keyword set the type of extrapolation . type string
644	@param npts: the numbers of q points of data to consider for extrapolation
645	@param function: a keyword to select the function to use for extrapolation.
646	of type string.
647	@param power: an power to apply power_low function
648
649	"""
650	range = range.lower()
651	if range not in ['high', 'low']:
652	raise ValueError, "Extrapolation range should be 'high' or 'low'"
653	function = function.lower()
654	if function not in ['power_law', 'guinier']:
655	raise ValueError, "Extrapolation function should be 'guinier' or 'power_law'"
656
657	if range == 'high':
658	if function != 'power_law':
659	raise ValueError, "Extrapolation only allows a power law at high Q"
660	self._high_extrapolation_npts = npts
661	self._high_extrapolation_power = power
662	else:
663	if function == 'power_law':
664	self._low_extrapolation_function = power_law
665	else:
666	self._low_extrapolation_function = guinier
667	self._low_extrapolation_npts = npts
668	self._low_extrapolation_power = power
669
670	def get_qstar(self, extrapolation=None):
671	"""
672	Compute the invariant of the local copy of data.
673	Implementation:
674	if slit smear:
675	qstar_0 = self._get_qstar_smear()
676	else:
677	qstar_0 = self._get_qstar_unsmear()
678	if extrapolation is None:
679	return qstar_0
680	if extrapolation==low:
681	return qstar_0 + self.get_qstar_low()
682	elif extrapolation==high:
683	return qstar_0 + self.get_qstar_high()
684	elif extrapolation==both:
685	return qstar_0 + self.get_qstar_low() + self.get_qstar_high()
686
687	@param extrapolation: string to apply optional extrapolation
688	@return q_star: invariant of the data within data's q range
689
690	@warning: When using setting data to Data1D , the user is responsible of
691	checking that the scale and the background are properly apply to the data
692
693	@warning: if error occur self.get_qstar_low() or self.get_qstar_low()
694	their returned value will be ignored
695	"""
696	qstar_0 = self._get_qstar(data=self._data)
697
698	if extrapolation is None:
699	self._qstar = qstar_0
700	return self._qstar
701	# Compute invariant plus invaraint of extrapolated data
702	extrapolation = extrapolation.lower()
703	if extrapolation == "low":
704	self._qstar = qstar_0 + self.get_qstar_low()
705	return self._qstar
706	elif extrapolation == "high":
707	self._qstar = qstar_0 + self.get_qstar_high()
708	return self._qstar
709	elif extrapolation == "both":
710	self._qstar = qstar_0 + self.get_qstar_low() + self.get_qstar_high()
711	return self._qstar
712
713	def get_surface(self,contrast, porod_const):
714	"""
715	Compute the surface of the data.
716
717	Implementation:
718	V= self.get_volume_fraction(contrast)
719
720	Compute the surface given by:
721	surface = (2pi V(1- V)*porod_const)/ q_star
722
723	@param contrast: contrast value to compute the volume
724	@param porod_const: Porod constant to compute the surface
725	@return: specific surface
726	"""
727	if contrast is None or porod_const is None:
728	return
729	#Check whether we have Q star
730	if self._qstar is None:
731	self._qstar = self.get_star()
732	if self._qstar == 0:
733	raise RuntimeError("Cannot compute surface, invariant value is zero")
734	# Compute the volume
735	volume = self.get_volume_fraction(contrast)
736	return 2 * math.pi * volume (1 - volume) float(porod_const)/self._qstar
737
738	def get_volume_fraction(self, contrast):
739	"""
740	Compute volume fraction is deduced as follow:
741
742	q_star = 2(picontrast)*2 volume( 1- volume)
743	for k = 10^(-8)q_star/(2(pi\|contrast\|)*2)
744	we get 2 values of volume:
745	with 1 - 4 * k >= 0
746	volume1 = (1- sqrt(1- 4*k))/2
747	volume2 = (1+ sqrt(1- 4*k))/2
748
749	q_star: the invariant value included extrapolation is applied
750	unit 1/A^(3)*1/cm
751	q_star = self.get_qstar()
752
753	the result returned will be 0<= volume <= 1
754
755	@param contrast: contrast value provides by the user of type float.
756	contrast unit is 1/A^(2)= 10^(16)cm^(2)
757	@return: volume fraction
758	@note: volume fraction must have no unit
759	"""
760	if contrast is None:
761	return
762	if contrast < 0:
763	raise ValueError, "contrast must be greater than zero"
764
765	if self._qstar is None:
766	self._qstar = self.get_qstar()
767
768	if self._qstar < 0:
769	raise RuntimeError, "invariant must be greater than zero"
770
771	# Compute intermediate constant
772	k = 1.e-8 * self._qstar/(2 * (math.pi * math.fabs(float(contrast)))**2)
773	#Check discriminant value
774	discrim = 1 - 4 * k
775
776	# Compute volume fraction
777	if discrim < 0:
778	raise RuntimeError, "could not compute the volume fraction: negative discriminant"
779	elif discrim == 0:
780	return 1/2
781	else:
782	volume1 = 0.5 * (1 - math.sqrt(discrim))
783	volume2 = 0.5 * (1 + math.sqrt(discrim))
784
785	if 0 <= volume1 and volume1 <= 1:
786	return volume1
787	elif 0 <= volume2 and volume2 <= 1:
788	return volume2
789	raise RuntimeError, "could not compute the volume fraction: inconsistent results"
790
791	def get_qstar_with_error(self, extrapolation=None):
792	"""
793	Compute the invariant uncertainty.
794	This uncertainty computation depends on whether or not the data is
795	smeared.
796	@return: invariant, the invariant uncertainty
797	return self._get_qstar(), self._get_qstar_smear_uncertainty()
798	"""
799	if self._qstar is None:
800	self._qstar = self.get_qstar(extrapolation=extrapolation)
801	if self._qstar_err is None:
802	self._qstar_err = self._get_qstar_smear_uncertainty()
803
804	return self._qstar, self._qstar_err
805
806	def get_volume_fraction_with_error(self, contrast):
807	"""
808	Compute uncertainty on volume value as well as the volume fraction
809	This uncertainty is given by the following equation:
810	dV = 0.5 * (4k dq_star) /(2* math.sqrt(1-k* q_star))
811
812	for k = 10^(-8)q_star/(2(pi\|contrast\|)*2)
813
814	q_star: the invariant value including extrapolated value if existing
815	dq_star: the invariant uncertainty
816	dV: the volume uncertainty
817	@param contrast: contrast value
818	@return: V, dV = self.get_volume_fraction_with_error(contrast), dV
819	"""
820	if contrast is None:
821	return
822	self._qstar, self._qstar_err = self.get_qstar_with_error()
823
824	volume = self.get_volume_fraction(contrast)
825	if self._qstar < 0:
826	raise ValueError, "invariant must be greater than zero"
827
828	k = 1.e-8 * self._qstar /(2 * (math.pi* math.fabs(float(contrast)))**2)
829	#check value inside the sqrt function
830	value = 1 - k * self._qstar
831	if (value) <= 0:
832	raise ValueError, "Cannot compute incertainty on volume"
833	# Compute uncertainty
834	uncertainty = (0.5 * 4 * k * self._qstar_err)/(2 * math.sqrt(1 - k * self._qstar))
835
836	return volume, math.fabs(uncertainty)
837
838	def get_surface_with_error(self, contrast, porod_const):
839	"""
840	Compute uncertainty of the surface value as well as thesurface value
841	this uncertainty is given as follow:
842
843	dS = porod_const 2pi[( dV -2VdV)/q_star
844	+ dq_star(v-v**2)
845
846	q_star: the invariant value including extrapolated value if existing
847	dq_star: the invariant uncertainty
848	V: the volume fraction value
849	dV: the volume uncertainty
850
851	@param contrast: contrast value
852	@param porod_const: porod constant value
853	@return S, dS: the surface, with its uncertainty
854	"""
855	if contrast is None or porod_const is None:
856	return
857	v, dv = self.get_volume_fraction_with_error(contrast)
858	self._qstar, self._qstar_err = self.get_qstar_with_error()
859	if self._qstar <= 0:
860	raise ValueError, "invariant must be greater than zero"
861	ds = porod_const * 2 * math.pi * (( dv - 2 * v * dv)/ self._qstar\
862	+ self._qstar_err * ( v - v**2))
863	s = self.get_surface(contrast=contrast, porod_const=porod_const)
864	return s, ds

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

SasView

source: sasview/Invariant/invariant.py @ c451be9

Download in other formats: