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

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalcmagnetic_scattrelease-4.2.2ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since b1f20d1 was b1f20d1, checked in by krzywon, 5 years ago

Make invariant calculator numpy 1.12+ compatible. refs #888

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