invariant.py @ b11e127

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 b11e127 was 07eb396, checked in by Jae Cho <jhjcho@…>, 12 years ago
Invaraint freeze fix
Property mode set to `100644`
File size: 36.2 KB

Line
1	#####################################################################
2	#This software was developed by the University of Tennessee as part of the
3	#Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
4	#project funded by the US National Science Foundation.
5	#See the license text in license.txt
6	#copyright 2010, University of Tennessee
7	######################################################################
8
9	"""
10	This module implements invariant and its related computations.
11
12	:author: Gervaise B. Alina/UTK
13	:author: Mathieu Doucet/UTK
14	:author: Jae Cho/UTK
15
16	"""
17	import math
18	import numpy
19
20	from sans.dataloader.data_info import Data1D as LoaderData1D
21
22	# The minimum q-value to be used when extrapolating
23	Q_MINIMUM = 1e-5
24
25	# The maximum q-value to be used when extrapolating
26	Q_MAXIMUM = 10
27
28	# Number of steps in the extrapolation
29	INTEGRATION_NSTEPS = 1000
30
31	class Transform(object):
32	"""
33	Define interface that need to compute a function or an inverse
34	function given some x, y
35	"""
36
37	def linearize_data(self, data):
38	"""
39	Linearize data so that a linear fit can be performed.
40	Filter out the data that can't be transformed.
41
42	:param data: LoadData1D instance
43
44	"""
45	# Check that the vector lengths are equal
46	assert(len(data.x)==len(data.y))
47	if data.dy is not None:
48	assert(len(data.x)==len(data.dy))
49	dy = data.dy
50	else:
51	dy = numpy.ones(len(data.y))
52
53	# Transform the data
54	data_points = zip(data.x, data.y, dy)
55
56	output_points = [(self.linearize_q_value(p[0]),
57	math.log(p[1]),
58	p[2]/p[1]) for p in data_points if p[0]>0 and \
59	p[1]>0 and p[2]>0]
60
61	x_out, y_out, dy_out = zip(*output_points)
62
63	# Create Data1D object
64	x_out = numpy.asarray(x_out)
65	y_out = numpy.asarray(y_out)
66	dy_out = numpy.asarray(dy_out)
67	linear_data = LoaderData1D(x=x_out, y=y_out, dy=dy_out)
68
69	return linear_data
70
71	def get_allowed_bins(self, data):
72	"""
73	Goes through the data points and returns a list of boolean values
74	to indicate whether each points is allowed by the model or not.
75
76	:param data: Data1D object
77	"""
78	return [p[0]>0 and p[1]>0 and p[2]>0 for p in zip(data.x, data.y,
79	data.dy)]
80
81	def linearize_q_value(self, value):
82	"""
83	Transform the input q-value for linearization
84	"""
85	return NotImplemented
86
87	def extract_model_parameters(self, constant, slope, dconstant=0, dslope=0):
88	"""
89	set private member
90	"""
91	return NotImplemented
92
93	def evaluate_model(self, x):
94	"""
95	Returns an array f(x) values where f is the Transform function.
96	"""
97	return NotImplemented
98
99	def evaluate_model_errors(self, x):
100	"""
101	Returns an array of I(q) errors
102	"""
103	return NotImplemented
104
105	class Guinier(Transform):
106	"""
107	class of type Transform that performs operations related to guinier
108	function
109	"""
110	def __init__(self, scale=1, radius=60):
111	Transform.__init__(self)
112	self.scale = scale
113	self.radius = radius
114	## Uncertainty of scale parameter
115	self.dscale = 0
116	## Unvertainty of radius parameter
117	self.dradius = 0
118
119	def linearize_q_value(self, value):
120	"""
121	Transform the input q-value for linearization
122
123	:param value: q-value
124
125	:return: q*q
126	"""
127	return value * value
128
129	def extract_model_parameters(self, constant, slope, dconstant=0, dslope=0):
130	"""
131	assign new value to the scale and the radius
132	"""
133	self.scale = math.exp(constant)
134	if slope > 0:
135	slope = 0.0
136	self.radius = math.sqrt(-3 * slope)
137	# Errors
138	self.dscale = math.exp(constant)*dconstant
139	if slope == 0.0:
140	n_zero = -1.0e-24
141	self.dradius = -3.0/2.0/math.sqrt(-3 * n_zero)*dslope
142	else:
143	self.dradius = -3.0/2.0/math.sqrt(-3 * slope)*dslope
144
145	return [self.radius, self.scale], [self.dradius, self.dscale]
146
147	def evaluate_model(self, x):
148	"""
149	return F(x)= scale* e-((radiusx)*2/3)
150	"""
151	return self._guinier(x)
152
153	def evaluate_model_errors(self, x):
154	"""
155	Returns the error on I(q) for the given array of q-values
156
157	:param x: array of q-values
158	"""
159	p1 = numpy.array([self.dscale * math.exp(-((self.radius * q)**2/3)) \
160	for q in x])
161	p2 = numpy.array([self.scale * math.exp(-((self.radius * q)**2/3))\
162	* (-(q*2/3)) 2 * self.radius * self.dradius for q in x])
163	diq2 = p1p1 + p2p2
164	return numpy.array( [math.sqrt(err) for err in diq2] )
165
166	def _guinier(self, x):
167	"""
168	Retrive the guinier function after apply an inverse guinier function
169	to x
170	Compute a F(x) = scale* e-((radiusx)*2/3).
171
172	:param x: a vector of q values
173	:param scale: the scale value
174	:param radius: the guinier radius value
175
176	:return: F(x)
177	"""
178	# transform the radius of coming from the inverse guinier function to a
179	# a radius of a guinier function
180	if self.radius <= 0:
181	msg = "Rg expected positive value, but got %s"%self.radius
182	raise ValueError(msg)
183	value = numpy.array([math.exp(-((self.radius * i)**2/3)) for i in x ])
184	return self.scale * value
185
186	class PowerLaw(Transform):
187	"""
188	class of type transform that perform operation related to power_law
189	function
190	"""
191	def __init__(self, scale=1, power=4):
192	Transform.__init__(self)
193	self.scale = scale
194	self.power = power
195
196	def linearize_q_value(self, value):
197	"""
198	Transform the input q-value for linearization
199
200	:param value: q-value
201
202	:return: log(q)
203	"""
204	return math.log(value)
205
206	def extract_model_parameters(self, constant, slope, dconstant=0, dslope=0):
207	"""
208	Assign new value to the scale and the power
209	"""
210	self.power = -slope
211	self.scale = math.exp(constant)
212
213	# Errors
214	self.dscale = math.exp(constant)*dconstant
215	self.dpower = -dslope
216
217	return [self.power, self.scale], [self.dpower, self.dscale]
218
219	def evaluate_model(self, x):
220	"""
221	given a scale and a radius transform x, y using a power_law
222	function
223	"""
224	return self._power_law(x)
225
226	def evaluate_model_errors(self, x):
227	"""
228	Returns the error on I(q) for the given array of q-values
229	:param x: array of q-values
230	"""
231	p1 = numpy.array([self.dscale * math.pow(q, -self.power) for q in x])
232	p2 = numpy.array([self.scale * self.power * math.pow(q, -self.power-1)\
233	* self.dpower for q in x])
234	diq2 = p1p1 + p2p2
235	return numpy.array( [math.sqrt(err) for err in diq2] )
236
237	def _power_law(self, x):
238	"""
239	F(x) = scale* (x)^(-power)
240	when power= 4. the model is porod
241	else power_law
242	The model has three parameters: ::
243	1. x: a vector of q values
244	2. power: power of the function
245	3. scale : scale factor value
246
247	:param x: array
248	:return: F(x)
249	"""
250	if self.power <= 0:
251	msg = "Power_law function expected positive power,"
252	msg += " but got %s"%self.power
253	raise ValueError(msg)
254	if self.scale <= 0:
255	msg = "scale expected positive value, but got %s"%self.scale
256	raise ValueError(msg)
257
258	value = numpy.array([ math.pow(i, -self.power) for i in x ])
259	return self.scale * value
260
261	class Extrapolator:
262	"""
263	Extrapolate I(q) distribution using a given model
264	"""
265	def __init__(self, data, model=None):
266	"""
267	Determine a and b given a linear equation y = ax + b
268
269	If a model is given, it will be used to linearize the data before
270	the extrapolation is performed. If None,
271	a simple linear fit will be done.
272
273	:param data: data containing x and y such as y = ax + b
274	:param model: optional Transform object
275	"""
276	self.data = data
277	self.model = model
278
279	# Set qmin as the lowest non-zero value
280	self.qmin = Q_MINIMUM
281	for q_value in self.data.x:
282	if q_value > 0:
283	self.qmin = q_value
284	break
285	self.qmax = max(self.data.x)
286
287	def fit(self, power=None, qmin=None, qmax=None):
288	"""
289	Fit data for y = ax + b return a and b
290
291	:param power: a fixed, otherwise None
292	:param qmin: Minimum Q-value
293	:param qmax: Maximum Q-value
294	"""
295	if qmin is None:
296	qmin = self.qmin
297	if qmax is None:
298	qmax = self.qmax
299
300	# Identify the bin range for the fit
301	idx = (self.data.x >= qmin) & (self.data.x <= qmax)
302
303	fx = numpy.zeros(len(self.data.x))
304
305	# Uncertainty
306	if type(self.data.dy)==numpy.ndarray and \
307	len(self.data.dy)==len(self.data.x) and \
308	numpy.all(self.data.dy>0):
309	sigma = self.data.dy
310	else:
311	sigma = numpy.ones(len(self.data.x))
312
313	# Compute theory data f(x)
314	fx[idx] = self.data.y[idx]
315
316	# Linearize the data
317	if self.model is not None:
318	linearized_data = self.model.linearize_data(\
319	LoaderData1D(self.data.x[idx],
320	fx[idx],
321	dy = sigma[idx]))
322	else:
323	linearized_data = LoaderData1D(self.data.x[idx],
324	fx[idx],
325	dy = sigma[idx])
326
327	##power is given only for function = power_law
328	if power != None:
329	sigma2 = linearized_data.dy * linearized_data.dy
330	a = -(power)
331	b = (numpy.sum(linearized_data.y/sigma2) \
332	- a*numpy.sum(linearized_data.x/sigma2))/numpy.sum(1.0/sigma2)
333
334
335	deltas = linearized_data.x*a + \
336	numpy.ones(len(linearized_data.x))*b-linearized_data.y
337	residuals = numpy.sum(deltas*deltas/sigma2)
338
339	err = math.fabs(residuals) / numpy.sum(1.0/sigma2)
340	return [a, b], [0, math.sqrt(err)]
341	else:
342	A = numpy.vstack([ linearized_data.x/linearized_data.dy,
343	1.0/linearized_data.dy]).T
344	(p, residuals, rank, s) = numpy.linalg.lstsq(A,
345	linearized_data.y/linearized_data.dy)
346
347	# Get the covariance matrix, defined as inv_cov = a_transposed * a
348	err = numpy.zeros(2)
349	try:
350	inv_cov = numpy.dot(A.transpose(), A)
351	cov = numpy.linalg.pinv(inv_cov)
352	err_matrix = math.fabs(residuals) * cov
353	err = [math.sqrt(err_matrix[0][0]), math.sqrt(err_matrix[1][1])]
354	except:
355	err = [-1.0, -1.0]
356
357	return p, err
358
359
360	class InvariantCalculator(object):
361	"""
362	Compute invariant if data is given.
363	Can provide volume fraction and surface area if the user provides
364	Porod constant and contrast values.
365
366	:precondition: the user must send a data of type DataLoader.Data1D
367	the user provide background and scale values.
368
369	:note: Some computations depends on each others.
370	"""
371	def __init__(self, data, background=0, scale=1 ):
372	"""
373	Initialize variables.
374
375	:param data: data must be of type DataLoader.Data1D
376	:param background: Background value. The data will be corrected
377	before processing
378	:param scale: Scaling factor for I(q). The data will be corrected
379	before processing
380	"""
381	# Background and scale should be private data member if the only way to
382	# change them are by instantiating a new object.
383	self._background = background
384	self._scale = scale
385	# slit height for smeared data
386	self._smeared = None
387	# The data should be private
388	self._data = self._get_data(data)
389	# get the dxl if the data is smeared: This is done only once on init.
390	if self._data.dxl != None and self._data.dxl.all() >0:
391	# assumes constant dxl
392	self._smeared = self._data.dxl[0]
393
394	# Since there are multiple variants of Q*, you should force the
395	# user to use the get method and keep Q* a private data member
396	self._qstar = None
397
398	# You should keep the error on Q* so you can reuse it without
399	# recomputing the whole thing.
400	self._qstar_err = 0
401
402	# Extrapolation parameters
403	self._low_extrapolation_npts = 4
404	self._low_extrapolation_function = Guinier()
405	self._low_extrapolation_power = None
406	self._low_extrapolation_power_fitted = None
407
408	self._high_extrapolation_npts = 4
409	self._high_extrapolation_function = PowerLaw()
410	self._high_extrapolation_power = None
411	self._high_extrapolation_power_fitted = None
412
413	# Extrapolation range
414	self._low_q_limit = Q_MINIMUM
415
416	def _get_data(self, data):
417	"""
418	:note: this function must be call before computing any type
419	of invariant
420
421	:return: new data = self._scale *data - self._background
422	"""
423	if not issubclass(data.__class__, LoaderData1D):
424	#Process only data that inherited from DataLoader.Data_info.Data1D
425	raise ValueError,"Data must be of type DataLoader.Data1D"
426	#from copy import deepcopy
427	new_data = (self._scale * data) - self._background
428
429	# Check that the vector lengths are equal
430	assert(len(new_data.x)==len(new_data.y))
431
432	# Verify that the errors are set correctly
433	if new_data.dy is None or len(new_data.x) != len(new_data.dy) or \
434	(min(new_data.dy)==0 and max(new_data.dy)==0):
435	new_data.dy = numpy.ones(len(new_data.x))
436	return new_data
437
438	def _fit(self, model, qmin=Q_MINIMUM, qmax=Q_MAXIMUM, power=None):
439	"""
440	fit data with function using
441	data = self._get_data()
442	fx = Functor(data , function)
443	y = data.y
444	slope, constant = linalg.lstsq(y,fx)
445
446	:param qmin: data first q value to consider during the fit
447	:param qmax: data last q value to consider during the fit
448	:param power : power value to consider for power-law
449	:param function: the function to use during the fit
450
451	:return a: the scale of the function
452	:return b: the other parameter of the function for guinier will be radius
453	for power_law will be the power value
454	"""
455	extrapolator = Extrapolator(data=self._data, model=model)
456	p, dp = extrapolator.fit(power=power, qmin=qmin, qmax=qmax)
457
458	return model.extract_model_parameters(constant=p[1], slope=p[0],
459	dconstant=dp[1], dslope=dp[0])
460
461	def _get_qstar(self, data):
462	"""
463	Compute invariant for pinhole data.
464	This invariant is given by: ::
465
466	q_star = x0*2 y0 dx0 +x12 y1 *dx1
467	+ ..+ xn*2 yn *dxn for non smeared data
468
469	q_star = dxl0 x0 y0 dx0 +dxl1 x1 y1 dx1
470	+ ..+ dlxn xn yn *dxn for smeared data
471
472	where n >= len(data.x)-1
473	dxl = slit height dQl
474	dxi = 1/2*(xi+1 - xi) + (xi - xi-1)
475	dx0 = (x1 - x0)/2
476	dxn = (xn - xn-1)/2
477
478	:param data: the data to use to compute invariant.
479
480	:return q_star: invariant value for pinhole data. q_star > 0
481	"""
482	if len(data.x) <= 1 or len(data.y) <= 1 or len(data.x)!= len(data.y):
483	msg = "Length x and y must be equal"
484	msg += " and greater than 1; got x=%s, y=%s"%(len(data.x),
485	len(data.y))
486	raise ValueError, msg
487	else:
488	# Take care of smeared data
489	if self._smeared is None:
490	gx = data.x * data.x
491	# assumes that len(x) == len(dxl).
492	else:
493	gx = data.dxl * data.x
494
495	n = len(data.x)- 1
496	#compute the first delta q
497	dx0 = (data.x[1] - data.x[0])/2
498	#compute the last delta q
499	dxn = (data.x[n] - data.x[n-1])/2
500	sum = 0
501	sum += gx[0] * data.y[0] * dx0
502	sum += gx[n] * data.y[n] * dxn
503
504	if len(data.x) == 2:
505	return sum
506	else:
507	#iterate between for element different
508	#from the first and the last
509	for i in xrange(1, n-1):
510	dxi = (data.x[i+1] - data.x[i-1])/2
511	sum += gx[i] * data.y[i] * dxi
512	return sum
513
514	def _get_qstar_uncertainty(self, data):
515	"""
516	Compute invariant uncertainty with with pinhole data.
517	This uncertainty is given as follow: ::
518
519	dq_star = math.sqrt[(x0*2(dy0)dx0)*2 +
520	(x1*2 (dy1)dx1)2 + ..+ (xn2 (dyn)dxn)*2 ]
521	where n >= len(data.x)-1
522	dxi = 1/2*(xi+1 - xi) + (xi - xi-1)
523	dx0 = (x1 - x0)/2
524	dxn = (xn - xn-1)/2
525	dyn: error on dy
526
527	:param data:
528	:note: if data doesn't contain dy assume dy= math.sqrt(data.y)
529	"""
530	if len(data.x) <= 1 or len(data.y) <= 1 or \
531	len(data.x) != len(data.y) or \
532	(data.dy is not None and (len(data.dy) != len(data.y))):
533	msg = "Length of data.x and data.y must be equal"
534	msg += " and greater than 1; got x=%s, y=%s"%(len(data.x),
535	len(data.y))
536	raise ValueError, msg
537	else:
538	#Create error for data without dy error
539	if data.dy is None:
540	dy = math.sqrt(data.y)
541	else:
542	dy = data.dy
543	# Take care of smeared data
544	if self._smeared is None:
545	gx = data.x * data.x
546	# assumes that len(x) == len(dxl).
547	else:
548	gx = data.dxl * data.x
549
550	n = len(data.x) - 1
551	#compute the first delta
552	dx0 = (data.x[1] - data.x[0])/2
553	#compute the last delta
554	dxn= (data.x[n] - data.x[n-1])/2
555	sum = 0
556	sum += (gx[0] * dy[0] * dx0)**2
557	sum += (gx[n] * dy[n] * dxn)**2
558	if len(data.x) == 2:
559	return math.sqrt(sum)
560	else:
561	#iterate between for element different
562	#from the first and the last
563	for i in xrange(1, n-1):
564	dxi = (data.x[i+1] - data.x[i-1])/2
565	sum += (gx[i] * dy[i] * dxi)**2
566	return math.sqrt(sum)
567
568	def _get_extrapolated_data(self, model, npts=INTEGRATION_NSTEPS,
569	q_start=Q_MINIMUM, q_end=Q_MAXIMUM):
570	"""
571	:return: extrapolate data create from data
572	"""
573	#create new Data1D to compute the invariant
574	q = numpy.linspace(start=q_start,
575	stop=q_end,
576	num=npts,
577	endpoint=True)
578	iq = model.evaluate_model(q)
579	diq = model.evaluate_model_errors(q)
580
581	result_data = LoaderData1D(x=q, y=iq, dy=diq)
582	if self._smeared != None:
583	result_data.dxl = self._smeared * numpy.ones(len(q))
584	return result_data
585
586	def get_data(self):
587	"""
588	:return: self._data
589	"""
590	return self._data
591
592	def get_extrapolation_power(self, range='high'):
593	"""
594	:return: the fitted power for power law function for a given
595	extrapolation range
596	"""
597	if range == 'low':
598	return self._low_extrapolation_power_fitted
599	return self._high_extrapolation_power_fitted
600
601	def get_qstar_low(self):
602	"""
603	Compute the invariant for extrapolated data at low q range.
604
605	Implementation:
606	data = self._get_extra_data_low()
607	return self._get_qstar()
608
609	:return q_star: the invariant for data extrapolated at low q.
610	"""
611	# Data boundaries for fitting
612	qmin = self._data.x[0]
613	qmax = self._data.x[self._low_extrapolation_npts - 1]
614
615	# Extrapolate the low-Q data
616	p, dp = self._fit(model=self._low_extrapolation_function,
617	qmin=qmin,
618	qmax=qmax,
619	power=self._low_extrapolation_power)
620	self._low_extrapolation_power_fitted = p[0]
621
622	# Distribution starting point
623	self._low_q_limit = Q_MINIMUM
624	if Q_MINIMUM >= qmin:
625	self._low_q_limit = qmin/10
626
627	data = self._get_extrapolated_data(\
628	model=self._low_extrapolation_function,
629	npts=INTEGRATION_NSTEPS,
630	q_start=self._low_q_limit, q_end=qmin)
631
632	# Systematic error
633	# If we have smearing, the shape of the I(q) distribution at low Q will
634	# may not be a Guinier or simple power law. The following is
635	# a conservative estimation for the systematic error.
636	err = qminqminmath.fabs((qmin-self._low_q_limit)*\
637	(data.y[0] - data.y[INTEGRATION_NSTEPS-1]))
638	return self._get_qstar(data), self._get_qstar_uncertainty(data)+err
639
640	def get_qstar_high(self):
641	"""
642	Compute the invariant for extrapolated data at high q range.
643
644	Implementation:
645	data = self._get_extra_data_high()
646	return self._get_qstar()
647
648	:return q_star: the invariant for data extrapolated at high q.
649	"""
650	# Data boundaries for fitting
651	x_len = len(self._data.x) - 1
652	qmin = self._data.x[x_len - (self._high_extrapolation_npts - 1)]
653	qmax = self._data.x[x_len]
654
655	# fit the data with a model to get the appropriate parameters
656	p, dp = self._fit(model=self._high_extrapolation_function,
657	qmin=qmin,
658	qmax=qmax,
659	power=self._high_extrapolation_power)
660	self._high_extrapolation_power_fitted = p[0]
661
662	#create new Data1D to compute the invariant
663	data = self._get_extrapolated_data(\
664	model=self._high_extrapolation_function,
665	npts=INTEGRATION_NSTEPS,
666	q_start=qmax, q_end=Q_MAXIMUM)
667
668	return self._get_qstar(data), self._get_qstar_uncertainty(data)
669
670	def get_extra_data_low(self, npts_in=None, q_start=None, npts=20):
671	"""
672	Returns the extrapolated data used for the loew-Q invariant calculation.
673	By default, the distribution will cover the data points used for the
674	extrapolation. The number of overlap points is a parameter (npts_in).
675	By default, the maximum q-value of the distribution will be
676	the minimum q-value used when extrapolating for the purpose of the
677	invariant calculation.
678
679	:param npts_in: number of data points for which
680	the extrapolated data overlap
681	:param q_start: is the minimum value to uses for extrapolated data
682	:param npts: the number of points in the extrapolated distribution
683
684	"""
685	# Get extrapolation range
686	if q_start is None:
687	q_start = self._low_q_limit
688
689	if npts_in is None:
690	npts_in = self._low_extrapolation_npts
691	q_end = self._data.x[max(0, npts_in-1)]
692
693	if q_start >= q_end:
694	return numpy.zeros(0), numpy.zeros(0)
695
696	return self._get_extrapolated_data(\
697	model=self._low_extrapolation_function,
698	npts=npts,
699	q_start=q_start, q_end=q_end)
700
701	def get_extra_data_high(self, npts_in=None, q_end=Q_MAXIMUM, npts=20):
702	"""
703	Returns the extrapolated data used for the high-Q invariant calculation.
704	By default, the distribution will cover the data points used for the
705	extrapolation. The number of overlap points is a parameter (npts_in).
706	By default, the maximum q-value of the distribution will be Q_MAXIMUM,
707	the maximum q-value used when extrapolating for the purpose of the
708	invariant calculation.
709
710	:param npts_in: number of data points for which the
711	extrapolated data overlap
712	:param q_end: is the maximum value to uses for extrapolated data
713	:param npts: the number of points in the extrapolated distribution
714	"""
715	# Get extrapolation range
716	if npts_in is None:
717	npts_in = self._high_extrapolation_npts
718	_npts = len(self._data.x)
719	q_start = self._data.x[min(_npts, _npts-npts_in)]
720
721	if q_start >= q_end:
722	return numpy.zeros(0), numpy.zeros(0)
723
724	return self._get_extrapolated_data(\
725	model=self._high_extrapolation_function,
726	npts=npts,
727	q_start=q_start, q_end=q_end)
728
729	def set_extrapolation(self, range, npts=4, function=None, power=None):
730	"""
731	Set the extrapolation parameters for the high or low Q-range.
732	Note that this does not turn extrapolation on or off.
733
734	:param range: a keyword set the type of extrapolation . type string
735	:param npts: the numbers of q points of data to consider
736	for extrapolation
737	:param function: a keyword to select the function to use
738	for extrapolation.
739	of type string.
740	:param power: an power to apply power_low function
741
742	"""
743	range = range.lower()
744	if range not in ['high', 'low']:
745	raise ValueError, "Extrapolation range should be 'high' or 'low'"
746	function = function.lower()
747	if function not in ['power_law', 'guinier']:
748	msg = "Extrapolation function should be 'guinier' or 'power_law'"
749	raise ValueError, msg
750
751	if range == 'high':
752	if function != 'power_law':
753	msg = "Extrapolation only allows a power law at high Q"
754	raise ValueError, msg
755	self._high_extrapolation_npts = npts
756	self._high_extrapolation_power = power
757	self._high_extrapolation_power_fitted = power
758	else:
759	if function == 'power_law':
760	self._low_extrapolation_function = PowerLaw()
761	else:
762	self._low_extrapolation_function = Guinier()
763	self._low_extrapolation_npts = npts
764	self._low_extrapolation_power = power
765	self._low_extrapolation_power_fitted = power
766
767	def get_qstar(self, extrapolation=None):
768	"""
769	Compute the invariant of the local copy of data.
770
771	:param extrapolation: string to apply optional extrapolation
772
773	:return q_star: invariant of the data within data's q range
774
775	:warning: When using setting data to Data1D ,
776	the user is responsible of
777	checking that the scale and the background are
778	properly apply to the data
779
780	"""
781	self._qstar = self._get_qstar(self._data)
782	self._qstar_err = self._get_qstar_uncertainty(self._data)
783
784	if extrapolation is None:
785	return self._qstar
786
787	# Compute invariant plus invariant of extrapolated data
788	extrapolation = extrapolation.lower()
789	if extrapolation == "low":
790	qs_low, dqs_low = self.get_qstar_low()
791	qs_hi, dqs_hi = 0, 0
792
793	elif extrapolation == "high":
794	qs_low, dqs_low = 0, 0
795	qs_hi, dqs_hi = self.get_qstar_high()
796
797	elif extrapolation == "both":
798	qs_low, dqs_low = self.get_qstar_low()
799	qs_hi, dqs_hi = self.get_qstar_high()
800
801	self._qstar += qs_low + qs_hi
802	self._qstar_err = math.sqrt(self._qstar_err*self._qstar_err \
803	+ dqs_lowdqs_low + dqs_hidqs_hi)
804
805	return self._qstar
806
807	def get_surface(self, contrast, porod_const, extrapolation=None):
808	"""
809	Compute the specific surface from the data.
810
811	Implementation::
812
813	V = self.get_volume_fraction(contrast, extrapolation)
814
815	Compute the surface given by:
816	surface = (2pi V(1- V)*porod_const)/ q_star
817
818	:param contrast: contrast value to compute the volume
819	:param porod_const: Porod constant to compute the surface
820	:param extrapolation: string to apply optional extrapolation
821
822	:return: specific surface
823	"""
824	# Compute the volume
825	volume = self.get_volume_fraction(contrast, extrapolation)
826	return 2 * math.pi * volume (1 - volume) \
827	float(porod_const)/self._qstar
828
829	def get_volume_fraction(self, contrast, extrapolation=None):
830	"""
831	Compute volume fraction is deduced as follow: ::
832
833	q_star = 2(picontrast)*2 volume( 1- volume)
834	for k = 10^(-8)q_star/(2(pi\|contrast\|)*2)
835	we get 2 values of volume:
836	with 1 - 4 * k >= 0
837	volume1 = (1- sqrt(1- 4*k))/2
838	volume2 = (1+ sqrt(1- 4*k))/2
839
840	q_star: the invariant value included extrapolation is applied
841	unit 1/A^(3)*1/cm
842	q_star = self.get_qstar()
843
844	the result returned will be 0 <= volume <= 1
845
846	:param contrast: contrast value provides by the user of type float.
847	contrast unit is 1/A^(2)= 10^(16)cm^(2)
848	:param extrapolation: string to apply optional extrapolation
849
850	:return: volume fraction
851
852	:note: volume fraction must have no unit
853	"""
854	if contrast <= 0:
855	raise ValueError, "The contrast parameter must be greater than zero"
856
857	# Make sure Q star is up to date
858	self.get_qstar(extrapolation)
859
860	if self._qstar <= 0:
861	msg = "Invalid invariant: Invariant Q* must be greater than zero"
862	raise RuntimeError, msg
863
864	# Compute intermediate constant
865	k = 1.e-8 * self._qstar/(2 * (math.pi * math.fabs(float(contrast)))**2)
866	# Check discriminant value
867	discrim = 1 - 4 * k
868
869	# Compute volume fraction
870	if discrim < 0:
871	msg = "Could not compute the volume fraction: negative discriminant"
872	raise RuntimeError, msg
873	elif discrim == 0:
874	return 1/2
875	else:
876	volume1 = 0.5 * (1 - math.sqrt(discrim))
877	volume2 = 0.5 * (1 + math.sqrt(discrim))
878
879	if 0 <= volume1 and volume1 <= 1:
880	return volume1
881	elif 0 <= volume2 and volume2 <= 1:
882	return volume2
883	msg = "Could not compute the volume fraction: inconsistent results"
884	raise RuntimeError, msg
885
886	def get_qstar_with_error(self, extrapolation=None):
887	"""
888	Compute the invariant uncertainty.
889	This uncertainty computation depends on whether or not the data is
890	smeared.
891
892	:param extrapolation: string to apply optional extrapolation
893
894	:return: invariant, the invariant uncertainty
895	"""
896	self.get_qstar(extrapolation)
897	return self._qstar, self._qstar_err
898
899	def get_volume_fraction_with_error(self, contrast, extrapolation=None):
900	"""
901	Compute uncertainty on volume value as well as the volume fraction
902	This uncertainty is given by the following equation: ::
903
904	dV = 0.5 * (4k dq_star) /(2* math.sqrt(1-k* q_star))
905
906	for k = 10^(-8)q_star/(2(pi\|contrast\|)*2)
907
908	q_star: the invariant value including extrapolated value if existing
909	dq_star: the invariant uncertainty
910	dV: the volume uncertainty
911
912	The uncertainty will be set to -1 if it can't be computed.
913
914	:param contrast: contrast value
915	:param extrapolation: string to apply optional extrapolation
916
917	:return: V, dV = volume fraction, error on volume fraction
918	"""
919	volume = self.get_volume_fraction(contrast, extrapolation)
920
921	# Compute error
922	k = 1.e-8 * self._qstar /(2 * (math.pi* math.fabs(float(contrast)))**2)
923	# Check value inside the sqrt function
924	value = 1 - k * self._qstar
925	if (value) <= 0:
926	uncertainty = -1
927	# Compute uncertainty
928	uncertainty = math.fabs((0.5 * 4 * k * \
929	self._qstar_err)/(2 * math.sqrt(1 - k * self._qstar)))
930
931	return volume, uncertainty
932
933	def get_surface_with_error(self, contrast, porod_const, extrapolation=None):
934	"""
935	Compute uncertainty of the surface value as well as the surface value.
936	The uncertainty is given as follow: ::
937
938	dS = porod_const 2pi[( dV -2VdV)/q_star
939	+ dq_star(v-v**2)
940
941	q_star: the invariant value
942	dq_star: the invariant uncertainty
943	V: the volume fraction value
944	dV: the volume uncertainty
945
946	:param contrast: contrast value
947	:param porod_const: porod constant value
948	:param extrapolation: string to apply optional extrapolation
949
950	:return S, dS: the surface, with its uncertainty
951	"""
952	# We get the volume fraction, with error
953	# get_volume_fraction_with_error calls get_volume_fraction
954	# get_volume_fraction calls get_qstar
955	# which computes Qstar and dQstar
956	v, dv = self.get_volume_fraction_with_error(contrast, extrapolation)
957
958	s = self.get_surface(contrast=contrast, porod_const=porod_const,
959	extrapolation=extrapolation)
960	ds = porod_const * 2 * math.pi * (( dv - 2 * v * dv)/ self._qstar\
961	+ self._qstar_err * ( v - v**2))
962
963	return s, ds

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

SasView

source: sasview/sansinvariant/src/sans/invariant/invariant.py @ b11e127

Download in other formats: