invariant.py @ 3c680c1

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

Last change on this file since 3c680c1 was dbfd307, checked in by Paul Kienzle <pkienzle@…>, 6 years ago
use np ≥ 1.14.0 default for lstsq rcond, but compatible with np < 1.14.0
Property mode set to `100644`
File size: 34.0 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	# CRUFT: numpy>=1.14.0 allows rcond=None for the following default
347	rcond = np.finfo(float).eps * max(A.shape)
348	p, residuals, _, _ = np.linalg.lstsq(A, linearized_data.y / linearized_data.dy,
349	rcond=rcond)
350
351	# Get the covariance matrix, defined as inv_cov = a_transposed * a
352	err = np.zeros(2)
353	try:
354	inv_cov = np.dot(A.transpose(), A)
355	cov = np.linalg.pinv(inv_cov)
356	err_matrix = math.fabs(residuals) * cov
357	err = [math.sqrt(err_matrix[0][0]), math.sqrt(err_matrix[1][1])]
358	except:
359	err = [-1.0, -1.0]
360
361	return p, err
362
363
364	class InvariantCalculator(object):
365	"""
366	Compute invariant if data is given.
367	Can provide volume fraction and surface area if the user provides
368	Porod constant and contrast values.
369
370	:precondition: the user must send a data of type DataLoader.Data1D
371	the user provide background and scale values.
372
373	:note: Some computations depends on each others.
374	"""
375	def __init__(self, data, background=0, scale=1):
376	"""
377	Initialize variables.
378
379	:param data: data must be of type DataLoader.Data1D
380	:param background: Background value. The data will be corrected
381	before processing
382	:param scale: Scaling factor for I(q). The data will be corrected
383	before processing
384	"""
385	# Background and scale should be private data member if the only way to
386	# change them are by instantiating a new object.
387	self._background = background
388	self._scale = scale
389	# slit height for smeared data
390	self._smeared = None
391	# The data should be private
392	self._data = self._get_data(data)
393	# get the dxl if the data is smeared: This is done only once on init.
394	if self._data.dxl is not None and self._data.dxl.all() > 0:
395	# assumes constant dxl
396	self._smeared = self._data.dxl[0]
397
398	# Since there are multiple variants of Q*, you should force the
399	# user to use the get method and keep Q* a private data member
400	self._qstar = None
401
402	# You should keep the error on Q* so you can reuse it without
403	# recomputing the whole thing.
404	self._qstar_err = 0
405
406	# Extrapolation parameters
407	self._low_extrapolation_npts = 4
408	self._low_extrapolation_function = Guinier()
409	self._low_extrapolation_power = None
410	self._low_extrapolation_power_fitted = None
411
412	self._high_extrapolation_npts = 4
413	self._high_extrapolation_function = PowerLaw()
414	self._high_extrapolation_power = None
415	self._high_extrapolation_power_fitted = None
416
417	# Extrapolation range
418	self._low_q_limit = Q_MINIMUM
419
420	def _get_data(self, data):
421	"""
422	:note: this function must be call before computing any type
423	of invariant
424
425	:return: new data = self._scale *data - self._background
426	"""
427	if not issubclass(data.__class__, LoaderData1D):
428	#Process only data that inherited from DataLoader.Data_info.Data1D
429	raise ValueError("Data must be of type DataLoader.Data1D")
430	#from copy import deepcopy
431	new_data = (self._scale * data) - self._background
432
433	# Check that the vector lengths are equal
434	assert len(new_data.x) == len(new_data.y)
435
436	# Verify that the errors are set correctly
437	if new_data.dy is None or len(new_data.x) != len(new_data.dy) or \
438	(min(new_data.dy) == 0 and max(new_data.dy) == 0):
439	new_data.dy = np.ones(len(new_data.x))
440	return new_data
441
442	def _fit(self, model, qmin=Q_MINIMUM, qmax=Q_MAXIMUM, power=None):
443	"""
444	fit data with function using
445	data = self._get_data()
446	fx = Functor(data , function)
447	y = data.y
448	slope, constant = linalg.lstsq(y,fx)
449
450	:param qmin: data first q value to consider during the fit
451	:param qmax: data last q value to consider during the fit
452	:param power : power value to consider for power-law
453	:param function: the function to use during the fit
454
455	:return a: the scale of the function
456	:return b: the other parameter of the function for guinier will be radius
457	for power_law will be the power value
458	"""
459	extrapolator = Extrapolator(data=self._data, model=model)
460	p, dp = extrapolator.fit(power=power, qmin=qmin, qmax=qmax)
461
462	return model.extract_model_parameters(constant=p[1], slope=p[0],
463	dconstant=dp[1], dslope=dp[0])
464
465	def _get_qstar(self, data):
466	"""
467	Compute invariant for pinhole data.
468	This invariant is given by: ::
469
470	q_star = x0*2 y0 dx0 +x12 y1 *dx1
471	+ ..+ xn*2 yn *dxn for non smeared data
472
473	q_star = dxl0 x0 y0 dx0 +dxl1 x1 y1 dx1
474	+ ..+ dlxn xn yn *dxn for smeared data
475
476	where n >= len(data.x)-1
477	dxl = slit height dQl
478	dxi = 1/2*(xi+1 - xi) + (xi - xi-1)
479	dx0 = (x1 - x0)/2
480	dxn = (xn - xn-1)/2
481
482	:param data: the data to use to compute invariant.
483
484	:return q_star: invariant value for pinhole data. q_star > 0
485	"""
486	if len(data.x) <= 1 or len(data.y) <= 1 or len(data.x) != len(data.y):
487	msg = "Length x and y must be equal"
488	msg += " and greater than 1; got x=%s, y=%s" % (len(data.x), len(data.y))
489	raise ValueError(msg)
490	else:
491	# Take care of smeared data
492	if self._smeared is None:
493	gx = data.x * data.x
494	# assumes that len(x) == len(dxl).
495	else:
496	gx = data.dxl * data.x
497
498	n = len(data.x) - 1
499	#compute the first delta q
500	dx0 = (data.x[1] - data.x[0]) / 2
501	#compute the last delta q
502	dxn = (data.x[n] - data.x[n - 1]) / 2
503	total = 0
504	total += gx[0] * data.y[0] * dx0
505	total += gx[n] * data.y[n] * dxn
506
507	if len(data.x) == 2:
508	return total
509	else:
510	#iterate between for element different
511	#from the first and the last
512	for i in range(1, n - 1):
513	dxi = (data.x[i + 1] - data.x[i - 1]) / 2
514	total += gx[i] * data.y[i] * dxi
515	return total
516
517	def _get_qstar_uncertainty(self, data):
518	"""
519	Compute invariant uncertainty with with pinhole data.
520	This uncertainty is given as follow: ::
521
522	dq_star = math.sqrt[(x0*2(dy0)dx0)*2 +
523	(x1*2 (dy1)dx1)2 + ..+ (xn2 (dyn)dxn)*2 ]
524	where n >= len(data.x)-1
525	dxi = 1/2*(xi+1 - xi) + (xi - xi-1)
526	dx0 = (x1 - x0)/2
527	dxn = (xn - xn-1)/2
528	dyn: error on dy
529
530	:param data:
531	:note: if data doesn't contain dy assume dy= math.sqrt(data.y)
532	"""
533	if len(data.x) <= 1 or len(data.y) <= 1 or \
534	len(data.x) != len(data.y) or \
535	(data.dy is not None and (len(data.dy) != len(data.y))):
536	msg = "Length of data.x and data.y must be equal"
537	msg += " and greater than 1; got x=%s, y=%s" % (len(data.x), len(data.y))
538	raise ValueError(msg)
539	else:
540	#Create error for data without dy error
541	if data.dy is None:
542	dy = math.sqrt(data.y)
543	else:
544	dy = data.dy
545	# Take care of smeared data
546	if self._smeared is None:
547	gx = data.x * data.x
548	# assumes that len(x) == len(dxl).
549	else:
550	gx = data.dxl * data.x
551
552	n = len(data.x) - 1
553	#compute the first delta
554	dx0 = (data.x[1] - data.x[0]) / 2
555	#compute the last delta
556	dxn = (data.x[n] - data.x[n - 1]) / 2
557	total = 0
558	total += (gx[0] * dy[0] * dx0) ** 2
559	total += (gx[n] * dy[n] * dxn) ** 2
560	if len(data.x) == 2:
561	return math.sqrt(total)
562	else:
563	#iterate between for element different
564	#from the first and the last
565	for i in range(1, n - 1):
566	dxi = (data.x[i + 1] - data.x[i - 1]) / 2
567	total += (gx[i] * dy[i] * dxi) ** 2
568	return math.sqrt(total)
569
570	def _get_extrapolated_data(self, model, npts=INTEGRATION_NSTEPS,
571	q_start=Q_MINIMUM, q_end=Q_MAXIMUM):
572	"""
573	:return: extrapolate data create from data
574	"""
575	#create new Data1D to compute the invariant
576	q = np.linspace(start=q_start,
577	stop=q_end,
578	num=npts,
579	endpoint=True)
580	iq = model.evaluate_model(q)
581	diq = model.evaluate_model_errors(q)
582
583	result_data = LoaderData1D(x=q, y=iq, dy=diq)
584	if self._smeared is not None:
585	result_data.dxl = self._smeared * np.ones(len(q))
586	return result_data
587
588	def get_data(self):
589	"""
590	:return: self._data
591	"""
592	return self._data
593
594	def get_extrapolation_power(self, range='high'):
595	"""
596	:return: the fitted power for power law function for a given
597	extrapolation range
598	"""
599	if range == 'low':
600	return self._low_extrapolation_power_fitted
601	return self._high_extrapolation_power_fitted
602
603	def get_qstar_low(self):
604	"""
605	Compute the invariant for extrapolated data at low q range.
606
607	Implementation:
608	data = self._get_extra_data_low()
609	return self._get_qstar()
610
611	:return q_star: the invariant for data extrapolated at low q.
612	"""
613	# Data boundaries for fitting
614	qmin = self._data.x[0]
615	qmax = self._data.x[int(self._low_extrapolation_npts - 1)]
616
617	# Extrapolate the low-Q data
618	p, _ = self._fit(model=self._low_extrapolation_function,
619	qmin=qmin,
620	qmax=qmax,
621	power=self._low_extrapolation_power)
622	self._low_extrapolation_power_fitted = p[0]
623
624	# Distribution starting point
625	self._low_q_limit = Q_MINIMUM
626	if Q_MINIMUM >= qmin:
627	self._low_q_limit = qmin / 10
628
629	data = self._get_extrapolated_data(\
630	model=self._low_extrapolation_function,
631	npts=INTEGRATION_NSTEPS,
632	q_start=self._low_q_limit, q_end=qmin)
633
634	# Systematic error
635	# If we have smearing, the shape of the I(q) distribution at low Q will
636	# may not be a Guinier or simple power law. The following is
637	# a conservative estimation for the systematic error.
638	err = qmin * qmin * math.fabs((qmin - self._low_q_limit) * \
639	(data.y[0] - data.y[INTEGRATION_NSTEPS - 1]))
640	return self._get_qstar(data), self._get_qstar_uncertainty(data) + err
641
642	def get_qstar_high(self):
643	"""
644	Compute the invariant for extrapolated data at high q range.
645
646	Implementation:
647	data = self._get_extra_data_high()
648	return self._get_qstar()
649
650	:return q_star: the invariant for data extrapolated at high q.
651	"""
652	# Data boundaries for fitting
653	x_len = len(self._data.x) - 1
654	qmin = self._data.x[int(x_len - (self._high_extrapolation_npts - 1))]
655	qmax = self._data.x[int(x_len)]
656
657	# fit the data with a model to get the appropriate parameters
658	p, _ = self._fit(model=self._high_extrapolation_function,
659	qmin=qmin,
660	qmax=qmax,
661	power=self._high_extrapolation_power)
662	self._high_extrapolation_power_fitted = p[0]
663
664	#create new Data1D to compute the invariant
665	data = self._get_extrapolated_data(\
666	model=self._high_extrapolation_function,
667	npts=INTEGRATION_NSTEPS,
668	q_start=qmax, q_end=Q_MAXIMUM)
669
670	return self._get_qstar(data), self._get_qstar_uncertainty(data)
671
672	def get_extra_data_low(self, npts_in=None, q_start=None, npts=20):
673	"""
674	Returns the extrapolated data used for the loew-Q invariant calculation.
675	By default, the distribution will cover the data points used for the
676	extrapolation. The number of overlap points is a parameter (npts_in).
677	By default, the maximum q-value of the distribution will be
678	the minimum q-value used when extrapolating for the purpose of the
679	invariant calculation.
680
681	:param npts_in: number of data points for which
682	the extrapolated data overlap
683	:param q_start: is the minimum value to uses for extrapolated data
684	:param npts: the number of points in the extrapolated distribution
685
686	"""
687	# Get extrapolation range
688	if q_start is None:
689	q_start = self._low_q_limit
690
691	if npts_in is None:
692	npts_in = self._low_extrapolation_npts
693	q_end = self._data.x[max(0, int(npts_in - 1))]
694
695	if q_start >= q_end:
696	return np.zeros(0), np.zeros(0)
697
698	return self._get_extrapolated_data(\
699	model=self._low_extrapolation_function,
700	npts=npts,
701	q_start=q_start, q_end=q_end)
702
703	def get_extra_data_high(self, npts_in=None, q_end=Q_MAXIMUM, npts=20):
704	"""
705	Returns the extrapolated data used for the high-Q invariant calculation.
706	By default, the distribution will cover the data points used for the
707	extrapolation. The number of overlap points is a parameter (npts_in).
708	By default, the maximum q-value of the distribution will be Q_MAXIMUM,
709	the maximum q-value used when extrapolating for the purpose of the
710	invariant calculation.
711
712	:param npts_in: number of data points for which the
713	extrapolated data overlap
714	:param q_end: is the maximum value to uses for extrapolated data
715	:param npts: the number of points in the extrapolated distribution
716	"""
717	# Get extrapolation range
718	if npts_in is None:
719	npts_in = int(self._high_extrapolation_npts)
720	_npts = len(self._data.x)
721	q_start = self._data.x[min(_npts, int(_npts - npts_in))]
722
723	if q_start >= q_end:
724	return np.zeros(0), np.zeros(0)
725
726	return self._get_extrapolated_data(\
727	model=self._high_extrapolation_function,
728	npts=npts,
729	q_start=q_start, q_end=q_end)
730
731	def set_extrapolation(self, range, npts=4, function=None, power=None):
732	"""
733	Set the extrapolation parameters for the high or low Q-range.
734	Note that this does not turn extrapolation on or off.
735
736	:param range: a keyword set the type of extrapolation . type string
737	:param npts: the numbers of q points of data to consider
738	for extrapolation
739	:param function: a keyword to select the function to use
740	for extrapolation.
741	of type string.
742	:param power: an power to apply power_low function
743
744	"""
745	range = range.lower()
746	if range not in ['high', 'low']:
747	raise ValueError("Extrapolation range should be 'high' or 'low'")
748	function = function.lower()
749	if function not in ['power_law', 'guinier']:
750	msg = "Extrapolation function should be 'guinier' or 'power_law'"
751	raise ValueError(msg)
752
753	if range == 'high':
754	if function != 'power_law':
755	msg = "Extrapolation only allows a power law at high Q"
756	raise ValueError(msg)
757	self._high_extrapolation_npts = npts
758	self._high_extrapolation_power = power
759	self._high_extrapolation_power_fitted = power
760	else:
761	if function == 'power_law':
762	self._low_extrapolation_function = PowerLaw()
763	else:
764	self._low_extrapolation_function = Guinier()
765	self._low_extrapolation_npts = npts
766	self._low_extrapolation_power = power
767	self._low_extrapolation_power_fitted = power
768
769	def get_qstar(self, extrapolation=None):
770	"""
771	Compute the invariant of the local copy of data.
772
773	:param extrapolation: string to apply optional extrapolation
774
775	:return q_star: invariant of the data within data's q range
776
777	:warning: When using setting data to Data1D ,
778	the user is responsible of
779	checking that the scale and the background are
780	properly apply to the data
781
782	"""
783	self._qstar = self._get_qstar(self._data)
784	self._qstar_err = self._get_qstar_uncertainty(self._data)
785
786	if extrapolation is None:
787	return self._qstar
788
789	# Compute invariant plus invariant of extrapolated data
790	extrapolation = extrapolation.lower()
791	if extrapolation == "low":
792	qs_low, dqs_low = self.get_qstar_low()
793	qs_hi, dqs_hi = 0, 0
794
795	elif extrapolation == "high":
796	qs_low, dqs_low = 0, 0
797	qs_hi, dqs_hi = self.get_qstar_high()
798
799	elif extrapolation == "both":
800	qs_low, dqs_low = self.get_qstar_low()
801	qs_hi, dqs_hi = self.get_qstar_high()
802
803	self._qstar += qs_low + qs_hi
804	self._qstar_err = math.sqrt(self._qstar_err * self._qstar_err \
805	+ dqs_low * dqs_low + dqs_hi * dqs_hi)
806
807	return self._qstar
808
809	def get_surface(self, contrast, porod_const, extrapolation=None):
810	"""
811	Compute the specific surface from the data.
812
813	Implementation::
814
815	V = self.get_volume_fraction(contrast, extrapolation)
816
817	Compute the surface given by:
818	surface = (2pi V(1- V)*porod_const)/ q_star
819
820	:param contrast: contrast value to compute the volume
821	:param porod_const: Porod constant to compute the surface
822	:param extrapolation: string to apply optional extrapolation
823
824	:return: specific surface
825	"""
826	# Compute the volume
827	volume = self.get_volume_fraction(contrast, extrapolation)
828	return 2 * math.pi * volume * (1 - volume) * \
829	float(porod_const) / self._qstar
830
831	def get_volume_fraction(self, contrast, extrapolation=None):
832	"""
833	Compute volume fraction is deduced as follow: ::
834
835	q_star = 2(picontrast)*2 volume( 1- volume)
836	for k = 10^(-8)q_star/(2(pi\|contrast\|)*2)
837	we get 2 values of volume:
838	with 1 - 4 * k >= 0
839	volume1 = (1- sqrt(1- 4*k))/2
840	volume2 = (1+ sqrt(1- 4*k))/2
841
842	q_star: the invariant value included extrapolation is applied
843	unit 1/A^(3)*1/cm
844	q_star = self.get_qstar()
845
846	the result returned will be 0 <= volume <= 1
847
848	:param contrast: contrast value provides by the user of type float.
849	contrast unit is 1/A^(2)= 10^(16)cm^(2)
850	:param extrapolation: string to apply optional extrapolation
851
852	:return: volume fraction
853
854	:note: volume fraction must have no unit
855	"""
856	if contrast <= 0:
857	raise ValueError("The contrast parameter must be greater than zero")
858
859	# Make sure Q star is up to date
860	self.get_qstar(extrapolation)
861
862	if self._qstar <= 0:
863	msg = "Invalid invariant: Invariant Q* must be greater than zero"
864	raise RuntimeError(msg)
865
866	# Compute intermediate constant
867	k = 1.e-8 * self._qstar / (2 * (math.pi * math.fabs(float(contrast))) ** 2)
868	# Check discriminant value
869	discrim = 1 - 4 * k
870
871	# Compute volume fraction
872	if discrim < 0:
873	msg = "Could not compute the volume fraction: negative discriminant"
874	raise RuntimeError(msg)
875	elif discrim == 0:
876	return 1 / 2
877	else:
878	volume1 = 0.5 * (1 - math.sqrt(discrim))
879	volume2 = 0.5 * (1 + math.sqrt(discrim))
880
881	if 0 <= volume1 and volume1 <= 1:
882	return volume1
883	elif 0 <= volume2 and volume2 <= 1:
884	return volume2
885	msg = "Could not compute the volume fraction: inconsistent results"
886	raise RuntimeError(msg)
887
888	def get_qstar_with_error(self, extrapolation=None):
889	"""
890	Compute the invariant uncertainty.
891	This uncertainty computation depends on whether or not the data is
892	smeared.
893
894	:param extrapolation: string to apply optional extrapolation
895
896	:return: invariant, the invariant uncertainty
897	"""
898	self.get_qstar(extrapolation)
899	return self._qstar, self._qstar_err
900
901	def get_volume_fraction_with_error(self, contrast, extrapolation=None):
902	"""
903	Compute uncertainty on volume value as well as the volume fraction
904	This uncertainty is given by the following equation: ::
905
906	dV = 0.5 * (4k dq_star) /(2* math.sqrt(1-k* q_star))
907
908	for k = 10^(-8)q_star/(2(pi\|contrast\|)*2)
909
910	q_star: the invariant value including extrapolated value if existing
911	dq_star: the invariant uncertainty
912	dV: the volume uncertainty
913
914	The uncertainty will be set to -1 if it can't be computed.
915
916	:param contrast: contrast value
917	:param extrapolation: string to apply optional extrapolation
918
919	:return: V, dV = volume fraction, error on volume fraction
920	"""
921	volume = self.get_volume_fraction(contrast, extrapolation)
922
923	# Compute error
924	k = 1.e-8 * self._qstar / (2 * (math.pi * math.fabs(float(contrast))) ** 2)
925	# Check value inside the sqrt function
926	value = 1 - k * self._qstar
927	if (value) <= 0:
928	uncertainty = -1
929	# Compute uncertainty
930	uncertainty = math.fabs((0.5 * 4 * k * \
931	self._qstar_err) / (2 * math.sqrt(1 - k * self._qstar)))
932
933	return volume, uncertainty
934
935	def get_surface_with_error(self, contrast, porod_const, extrapolation=None):
936	"""
937	Compute uncertainty of the surface value as well as the surface value.
938	The uncertainty is given as follow: ::
939
940	dS = porod_const 2pi[( dV -2VdV)/q_star
941	+ dq_star(v-v**2)
942
943	q_star: the invariant value
944	dq_star: the invariant uncertainty
945	V: the volume fraction value
946	dV: the volume uncertainty
947
948	:param contrast: contrast value
949	:param porod_const: porod constant value
950	:param extrapolation: string to apply optional extrapolation
951
952	:return S, dS: the surface, with its uncertainty
953	"""
954	# We get the volume fraction, with error
955	# get_volume_fraction_with_error calls get_volume_fraction
956	# get_volume_fraction calls get_qstar
957	# which computes Qstar and dQstar
958	v, dv = self.get_volume_fraction_with_error(contrast, extrapolation)
959
960	s = self.get_surface(contrast=contrast, porod_const=porod_const,
961	extrapolation=extrapolation)
962	ds = porod_const * 2 * math.pi * ((dv - 2 * v * dv) / self._qstar\
963	+ self._qstar_err * (v - v ** 2))
964
965	return s, ds

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

SasView

source: sasview/src/sas/sascalc/invariant/invariant.py @ 3c680c1

Download in other formats: