source: sasview/Invariant/invariant.py @ 935ae18

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since 935ae18 was c4f52e3, checked in by Mathieu Doucet <doucetm@…>, 15 years ago

invariant: minor clean up

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