source: sasview/src/sas/invariant/invariant.py @ 898a8b9

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 898a8b9 was 79492222, checked in by krzywon, 10 years ago

Changed the file and folder names to remove all SANS references.

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