invariant.py @ 4cbb2f5

magnetic_scattrelease-4.2.2ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1249

Last change on this file since 4cbb2f5 was e090ba90, checked in by Paul Kienzle <pkienzle@…>, 6 years ago
remove errors and warnings from py37 tests of sascalc
Property mode set to `100644`
File size: 33.9 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 as np
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 = np.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 = np.asarray(x_out)
66	y_out = np.asarray(y_out)
67	dy_out = np.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 = np.array([self.dscale * math.exp(-((self.radius * q) ** 2 / 3)) \
161	for q in x])
162	p2 = np.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 np.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 = np.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 = np.array([self.dscale * math.pow(q, -self.power) for q in x])
235	p2 = np.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 np.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 = np.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 = np.zeros(len(self.data.x))
307
308	# Uncertainty
309	if type(self.data.dy) == np.ndarray and \
310	len(self.data.dy) == len(self.data.x) and \
311	np.all(self.data.dy > 0):
312	sigma = self.data.dy
313	else:
314	sigma = np.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 is not None:
332	sigma2 = linearized_data.dy * linearized_data.dy
333	a = -(power)
334	b = (np.sum(linearized_data.y / sigma2) \
335	- a * np.sum(linearized_data.x / sigma2)) / np.sum(1.0 / sigma2)
336
337
338	deltas = linearized_data.x * a + \
339	np.ones(len(linearized_data.x)) * b - linearized_data.y
340	residuals = np.sum(deltas * deltas / sigma2)
341
342	err = math.fabs(residuals) / np.sum(1.0 / sigma2)
343	return [a, b], [0, math.sqrt(err)]
344	else:
345	A = np.vstack([linearized_data.x / linearized_data.dy, 1.0 / linearized_data.dy]).T
346	p, residuals, _, _ = np.linalg.lstsq(A, linearized_data.y / linearized_data.dy,
347	rcond=None)
348
349	# Get the covariance matrix, defined as inv_cov = a_transposed * a
350	err = np.zeros(2)
351	try:
352	inv_cov = np.dot(A.transpose(), A)
353	cov = np.linalg.pinv(inv_cov)
354	err_matrix = math.fabs(residuals) * cov
355	err = [math.sqrt(err_matrix[0][0]), math.sqrt(err_matrix[1][1])]
356	except:
357	err = [-1.0, -1.0]
358
359	return p, err
360
361
362	class InvariantCalculator(object):
363	"""
364	Compute invariant if data is given.
365	Can provide volume fraction and surface area if the user provides
366	Porod constant and contrast values.
367
368	:precondition: the user must send a data of type DataLoader.Data1D
369	the user provide background and scale values.
370
371	:note: Some computations depends on each others.
372	"""
373	def __init__(self, data, background=0, scale=1):
374	"""
375	Initialize variables.
376
377	:param data: data must be of type DataLoader.Data1D
378	:param background: Background value. The data will be corrected
379	before processing
380	:param scale: Scaling factor for I(q). The data will be corrected
381	before processing
382	"""
383	# Background and scale should be private data member if the only way to
384	# change them are by instantiating a new object.
385	self._background = background
386	self._scale = scale
387	# slit height for smeared data
388	self._smeared = None
389	# The data should be private
390	self._data = self._get_data(data)
391	# get the dxl if the data is smeared: This is done only once on init.
392	if self._data.dxl is not None and self._data.dxl.all() > 0:
393	# assumes constant dxl
394	self._smeared = self._data.dxl[0]
395
396	# Since there are multiple variants of Q*, you should force the
397	# user to use the get method and keep Q* a private data member
398	self._qstar = None
399
400	# You should keep the error on Q* so you can reuse it without
401	# recomputing the whole thing.
402	self._qstar_err = 0
403
404	# Extrapolation parameters
405	self._low_extrapolation_npts = 4
406	self._low_extrapolation_function = Guinier()
407	self._low_extrapolation_power = None
408	self._low_extrapolation_power_fitted = None
409
410	self._high_extrapolation_npts = 4
411	self._high_extrapolation_function = PowerLaw()
412	self._high_extrapolation_power = None
413	self._high_extrapolation_power_fitted = None
414
415	# Extrapolation range
416	self._low_q_limit = Q_MINIMUM
417
418	def _get_data(self, data):
419	"""
420	:note: this function must be call before computing any type
421	of invariant
422
423	:return: new data = self._scale *data - self._background
424	"""
425	if not issubclass(data.__class__, LoaderData1D):
426	#Process only data that inherited from DataLoader.Data_info.Data1D
427	raise ValueError("Data must be of type DataLoader.Data1D")
428	#from copy import deepcopy
429	new_data = (self._scale * data) - self._background
430
431	# Check that the vector lengths are equal
432	assert len(new_data.x) == len(new_data.y)
433
434	# Verify that the errors are set correctly
435	if new_data.dy is None or len(new_data.x) != len(new_data.dy) or \
436	(min(new_data.dy) == 0 and max(new_data.dy) == 0):
437	new_data.dy = np.ones(len(new_data.x))
438	return new_data
439
440	def _fit(self, model, qmin=Q_MINIMUM, qmax=Q_MAXIMUM, power=None):
441	"""
442	fit data with function using
443	data = self._get_data()
444	fx = Functor(data , function)
445	y = data.y
446	slope, constant = linalg.lstsq(y,fx)
447
448	:param qmin: data first q value to consider during the fit
449	:param qmax: data last q value to consider during the fit
450	:param power : power value to consider for power-law
451	:param function: the function to use during the fit
452
453	:return a: the scale of the function
454	:return b: the other parameter of the function for guinier will be radius
455	for power_law will be the power value
456	"""
457	extrapolator = Extrapolator(data=self._data, model=model)
458	p, dp = extrapolator.fit(power=power, qmin=qmin, qmax=qmax)
459
460	return model.extract_model_parameters(constant=p[1], slope=p[0],
461	dconstant=dp[1], dslope=dp[0])
462
463	def _get_qstar(self, data):
464	"""
465	Compute invariant for pinhole data.
466	This invariant is given by: ::
467
468	q_star = x0*2 y0 dx0 +x12 y1 *dx1
469	+ ..+ xn*2 yn *dxn for non smeared data
470
471	q_star = dxl0 x0 y0 dx0 +dxl1 x1 y1 dx1
472	+ ..+ dlxn xn yn *dxn for smeared data
473
474	where n >= len(data.x)-1
475	dxl = slit height dQl
476	dxi = 1/2*(xi+1 - xi) + (xi - xi-1)
477	dx0 = (x1 - x0)/2
478	dxn = (xn - xn-1)/2
479
480	:param data: the data to use to compute invariant.
481
482	:return q_star: invariant value for pinhole data. q_star > 0
483	"""
484	if len(data.x) <= 1 or len(data.y) <= 1 or len(data.x) != len(data.y):
485	msg = "Length x and y must be equal"
486	msg += " and greater than 1; got x=%s, y=%s" % (len(data.x), len(data.y))
487	raise ValueError(msg)
488	else:
489	# Take care of smeared data
490	if self._smeared is None:
491	gx = data.x * data.x
492	# assumes that len(x) == len(dxl).
493	else:
494	gx = data.dxl * data.x
495
496	n = len(data.x) - 1
497	#compute the first delta q
498	dx0 = (data.x[1] - data.x[0]) / 2
499	#compute the last delta q
500	dxn = (data.x[n] - data.x[n - 1]) / 2
501	total = 0
502	total += gx[0] * data.y[0] * dx0
503	total += gx[n] * data.y[n] * dxn
504
505	if len(data.x) == 2:
506	return total
507	else:
508	#iterate between for element different
509	#from the first and the last
510	for i in range(1, n - 1):
511	dxi = (data.x[i + 1] - data.x[i - 1]) / 2
512	total += gx[i] * data.y[i] * dxi
513	return total
514
515	def _get_qstar_uncertainty(self, data):
516	"""
517	Compute invariant uncertainty with with pinhole data.
518	This uncertainty is given as follow: ::
519
520	dq_star = math.sqrt[(x0*2(dy0)dx0)*2 +
521	(x1*2 (dy1)dx1)2 + ..+ (xn2 (dyn)dxn)*2 ]
522	where n >= len(data.x)-1
523	dxi = 1/2*(xi+1 - xi) + (xi - xi-1)
524	dx0 = (x1 - x0)/2
525	dxn = (xn - xn-1)/2
526	dyn: error on dy
527
528	:param data:
529	:note: if data doesn't contain dy assume dy= math.sqrt(data.y)
530	"""
531	if len(data.x) <= 1 or len(data.y) <= 1 or \
532	len(data.x) != len(data.y) or \
533	(data.dy is not None and (len(data.dy) != len(data.y))):
534	msg = "Length of data.x and data.y must be equal"
535	msg += " and greater than 1; got x=%s, y=%s" % (len(data.x), 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	total = 0
556	total += (gx[0] * dy[0] * dx0) ** 2
557	total += (gx[n] * dy[n] * dxn) ** 2
558	if len(data.x) == 2:
559	return math.sqrt(total)
560	else:
561	#iterate between for element different
562	#from the first and the last
563	for i in range(1, n - 1):
564	dxi = (data.x[i + 1] - data.x[i - 1]) / 2
565	total += (gx[i] * dy[i] * dxi) ** 2
566	return math.sqrt(total)
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 = np.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 is not None:
583	result_data.dxl = self._smeared * np.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[int(self._low_extrapolation_npts - 1)]
614
615	# Extrapolate the low-Q data
616	p, _ = 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 = qmin * qmin * math.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[int(x_len - (self._high_extrapolation_npts - 1))]
653	qmax = self._data.x[int(x_len)]
654
655	# fit the data with a model to get the appropriate parameters
656	p, _ = 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, int(npts_in - 1))]
692
693	if q_start >= q_end:
694	return np.zeros(0), np.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 = int(self._high_extrapolation_npts)
718	_npts = len(self._data.x)
719	q_start = self._data.x[min(_npts, int(_npts - npts_in))]
720
721	if q_start >= q_end:
722	return np.zeros(0), np.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_low * dqs_low + dqs_hi * dqs_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/src/sas/sascalc/invariant/invariant.py @ 4cbb2f5

Download in other formats: