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

magnetic_scattrelease-4.2.2ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1249
Last change on this file since dbfd307 was dbfd307, checked in by Paul Kienzle <pkienzle@…>, 5 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
RevLine 
[959eb01]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"""
11This 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"""
18import math
19import numpy as np
20
21from sas.sascalc.dataloader.data_info import Data1D as LoaderData1D
22
23# The minimum q-value to be used when extrapolating
24Q_MINIMUM = 1e-5
25
26# The maximum q-value to be used when extrapolating
27Q_MAXIMUM = 10
28
29# Number of steps in the extrapolation
30INTEGRATION_NSTEPS = 1000
31
32class 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
106class 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-((radius*x)**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-((radius*x)**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
187class 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
264class 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
[7432acb]331        if power is not None:
[959eb01]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
[dbfd307]346            # CRUFT: numpy>=1.14.0 allows rcond=None for the following default
347            rcond = np.finfo(float).eps * max(A.shape)
[e090ba90]348            p, residuals, _, _ = np.linalg.lstsq(A, linearized_data.y / linearized_data.dy,
[dbfd307]349                                                 rcond=rcond)
[959eb01]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
364class 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.
[7432acb]394        if self._data.dxl is not None and self._data.dxl.all() > 0:
[959eb01]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
[574adc7]429            raise ValueError("Data must be of type DataLoader.Data1D")
[959eb01]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 +x1**2 *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))
[574adc7]489            raise ValueError(msg)
[959eb01]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
[574adc7]512                for i in range(1, n - 1):
[959eb01]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 + ..+ (xn**2 *(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))
[574adc7]538            raise ValueError(msg)
[959eb01]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
[574adc7]565                for i in range(1, n - 1):
[959eb01]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)
[7432acb]584        if self._smeared is not None:
[959eb01]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]
[b1f20d1]615        qmax = self._data.x[int(self._low_extrapolation_npts - 1)]
[959eb01]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
[b1f20d1]654        qmin = self._data.x[int(x_len - (self._high_extrapolation_npts - 1))]
655        qmax = self._data.x[int(x_len)]
[959eb01]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
[b1f20d1]693        q_end = self._data.x[max(0, int(npts_in - 1))]
[959eb01]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:
[b1f20d1]719            npts_in = int(self._high_extrapolation_npts)
[959eb01]720        _npts = len(self._data.x)
[b1f20d1]721        q_start = self._data.x[min(_npts, int(_npts - npts_in))]
[959eb01]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']:
[574adc7]747            raise ValueError("Extrapolation range should be 'high' or 'low'")
[959eb01]748        function = function.lower()
749        if function not in ['power_law', 'guinier']:
750            msg = "Extrapolation function should be 'guinier' or 'power_law'"
[574adc7]751            raise ValueError(msg)
[959eb01]752
753        if range == 'high':
754            if function != 'power_law':
755                msg = "Extrapolation only allows a power law at high Q"
[574adc7]756                raise ValueError(msg)
[959eb01]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 = (2*pi *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*(pi*contrast)**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:
[574adc7]857            raise ValueError("The contrast parameter must be greater than zero")
[959eb01]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"
[574adc7]864            raise RuntimeError(msg)
[959eb01]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"
[574adc7]874            raise RuntimeError(msg)
[959eb01]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"
[574adc7]886            raise RuntimeError(msg)
[959eb01]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 * (4*k* 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 *2*pi[( dV -2*V*dV)/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.