invariant.py @ d7ee5866

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalcmagnetic_scattrelease-4.1.1release-4.1.2release-4.2.2ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload

Last change on this file since d7ee5866 was b699768, checked in by Piotr Rozyczko <piotr.rozyczko@…>, 9 years ago
Initial commit of the refactored SasCalc? module.
Property mode set to `100644`
File size: 34.8 KB

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

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

SasView

source: sasview/src/sas/sascalc/invariant/invariant.py @ d7ee5866

Download in other formats: