Changeset c76bdc3 in sasview


Ignore:
Timestamp:
Mar 4, 2015 12:30:43 PM (9 years ago)
Author:
Doucet, Mathieu <doucetm@…>
Branches:
master, ESS_GUI, ESS_GUI_Docs, ESS_GUI_batch_fitting, ESS_GUI_bumps_abstraction, ESS_GUI_iss1116, ESS_GUI_iss879, ESS_GUI_iss959, ESS_GUI_opencl, ESS_GUI_ordering, ESS_GUI_sync_sascalc, costrafo411, magnetic_scatt, release-4.1.1, release-4.1.2, release-4.2.2, release_4.0.1, ticket-1009, ticket-1094-headless, ticket-1242-2d-resolution, ticket-1243, ticket-1249, ticket885, unittest-saveload
Children:
ab13d33
Parents:
5f8fc78
Message:

pylint fixes

Location:
src/sas
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • src/sas/invariant/invariant.py

    r79492222 rc76bdc3  
     1# pylint: disable=invalid-name 
    12##################################################################### 
    23#This software was developed by the University of Tennessee as part of the 
    34#Distributed Data Analysis of Neutron Scattering Experiments (DANSE) 
    4 #project funded by the US National Science Foundation.  
     5#project funded by the US National Science Foundation. 
    56#See the license text in license.txt 
    67#copyright 2010, University of Tennessee 
     
    1516 
    1617""" 
    17 import math  
     18import math 
    1819import numpy 
    1920 
     
    2122 
    2223# The minimum q-value to be used when extrapolating 
    23 Q_MINIMUM  = 1e-5 
     24Q_MINIMUM = 1e-5 
    2425 
    2526# The maximum q-value to be used when extrapolating 
    26 Q_MAXIMUM  = 10 
     27Q_MAXIMUM = 10 
    2728 
    2829# Number of steps in the extrapolation 
     
    3233    """ 
    3334    Define interface that need to compute a function or an inverse 
    34     function given some x, y  
     35    function given some x, y 
    3536    """ 
    36      
     37 
    3738    def linearize_data(self, data): 
    3839        """ 
    39         Linearize data so that a linear fit can be performed.  
     40        Linearize data so that a linear fit can be performed. 
    4041        Filter out the data that can't be transformed. 
    41          
     42 
    4243        :param data: LoadData1D instance 
    43          
     44 
    4445        """ 
    4546        # Check that the vector lengths are equal 
    46         assert(len(data.x)==len(data.y)) 
     47        assert len(data.x) == len(data.y) 
    4748        if data.dy is not None: 
    48             assert(len(data.x)==len(data.dy)) 
     49            assert len(data.x) == len(data.dy) 
    4950            dy = data.dy 
    5051        else: 
    5152            dy = numpy.ones(len(data.y)) 
    52              
     53 
    5354        # Transform the data 
    5455        data_points = zip(data.x, data.y, dy) 
     
    5657        output_points = [(self.linearize_q_value(p[0]), 
    5758                          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          
     59                          p[2] / p[1]) for p in data_points if p[0] > 0 and \ 
     60                          p[1] > 0 and p[2] > 0] 
     61 
    6162        x_out, y_out, dy_out = zip(*output_points) 
    62         
     63 
    6364        # Create Data1D object 
    6465        x_out = numpy.asarray(x_out) 
     
    6667        dy_out = numpy.asarray(dy_out) 
    6768        linear_data = LoaderData1D(x=x_out, y=y_out, dy=dy_out) 
    68          
     69 
    6970        return linear_data 
    70      
     71 
    7172    def get_allowed_bins(self, data): 
    7273        """ 
    7374        Goes through the data points and returns a list of boolean values 
    7475        to indicate whether each points is allowed by the model or not. 
    75          
     76 
    7677        :param data: Data1D object 
    7778        """ 
    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          
     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 
    8182    def linearize_q_value(self, value): 
    8283        """ 
     
    9091        """ 
    9192        return NotImplemented 
    92       
     93 
    9394    def evaluate_model(self, x): 
    9495        """ 
     
    9697        """ 
    9798        return NotImplemented 
    98      
     99 
    99100    def evaluate_model_errors(self, x): 
    100101        """ 
     
    102103        """ 
    103104        return NotImplemented 
    104      
     105 
    105106class Guinier(Transform): 
    106107    """ 
    107     class of type Transform that performs operations related to guinier  
     108    class of type Transform that performs operations related to guinier 
    108109    function 
    109110    """ 
     
    113114        self.radius = radius 
    114115        ## Uncertainty of scale parameter 
    115         self.dscale  = 0 
     116        self.dscale = 0 
    116117        ## Unvertainty of radius parameter 
    117118        self.dradius = 0 
    118          
     119 
    119120    def linearize_q_value(self, value): 
    120121        """ 
    121122        Transform the input q-value for linearization 
    122          
     123 
    123124        :param value: q-value 
    124          
     125 
    125126        :return: q*q 
    126127        """ 
    127128        return value * value 
    128      
     129 
    129130    def extract_model_parameters(self, constant, slope, dconstant=0, dslope=0): 
    130131        """ 
     
    136137        self.radius = math.sqrt(-3 * slope) 
    137138        # Errors 
    138         self.dscale = math.exp(constant)*dconstant 
     139        self.dscale = math.exp(constant) * dconstant 
    139140        if slope == 0.0: 
    140141            n_zero = -1.0e-24 
    141             self.dradius = -3.0/2.0/math.sqrt(-3 * n_zero)*dslope 
     142            self.dradius = -3.0 / 2.0 / math.sqrt(-3 * n_zero) * dslope 
    142143        else: 
    143             self.dradius = -3.0/2.0/math.sqrt(-3 * slope)*dslope 
    144          
     144            self.dradius = -3.0 / 2.0 / math.sqrt(-3 * slope) * dslope 
     145 
    145146        return [self.radius, self.scale], [self.dradius, self.dscale] 
    146          
     147 
    147148    def evaluate_model(self, x): 
    148149        """ 
     
    150151        """ 
    151152        return self._guinier(x) 
    152               
     153 
    153154    def evaluate_model_errors(self, x): 
    154155        """ 
    155156        Returns the error on I(q) for the given array of q-values 
    156          
     157 
    157158        :param x: array of q-values 
    158159        """ 
    159         p1 = numpy.array([self.dscale * math.exp(-((self.radius * q)**2/3)) \ 
     160        p1 = numpy.array([self.dscale * math.exp(-((self.radius * q) ** 2 / 3)) \ 
    160161                          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               
     162        p2 = numpy.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 numpy.array([math.sqrt(err) for err in diq2]) 
     166 
    166167    def _guinier(self, x): 
    167168        """ 
     
    169170        to x 
    170171        Compute a F(x) = scale* e-((radius*x)**2/3). 
    171          
    172         :param x: a vector of q values  
     172 
     173        :param x: a vector of q values 
    173174        :param scale: the scale value 
    174175        :param radius: the guinier radius value 
    175          
     176 
    176177        :return: F(x) 
    177         """    
    178         # transform the radius of coming from the inverse guinier function to a  
     178        """ 
     179        # transform the radius of coming from the inverse guinier function to a 
    179180        # a radius of a guinier function 
    180181        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 ])  
     182            msg = "Rg expected positive value, but got %s" % self.radius 
     183            raise ValueError(msg) 
     184        value = numpy.array([math.exp(-((self.radius * i) ** 2 / 3)) for i in x]) 
    184185        return self.scale * value 
    185186 
    186187class PowerLaw(Transform): 
    187188    """ 
    188     class of type transform that perform operation related to power_law  
     189    class of type transform that perform operation related to power_law 
    189190    function 
    190191    """ 
     
    193194        self.scale = scale 
    194195        self.power = power 
    195     
     196        self.dscale = 0.0 
     197        self.dpower = 0.0 
     198 
    196199    def linearize_q_value(self, value): 
    197200        """ 
    198201        Transform the input q-value for linearization 
    199          
     202 
    200203        :param value: q-value 
    201          
     204 
    202205        :return: log(q) 
    203206        """ 
    204207        return math.log(value) 
    205      
     208 
    206209    def extract_model_parameters(self, constant, slope, dconstant=0, dslope=0): 
    207210        """ 
    208         Assign new value to the scale and the power  
     211        Assign new value to the scale and the power 
    209212        """ 
    210213        self.power = -slope 
    211214        self.scale = math.exp(constant) 
    212          
     215 
    213216        # Errors 
    214         self.dscale = math.exp(constant)*dconstant 
     217        self.dscale = math.exp(constant) * dconstant 
    215218        self.dpower = -dslope 
    216          
     219 
    217220        return [self.power, self.scale], [self.dpower, self.dscale] 
    218          
     221 
    219222    def evaluate_model(self, x): 
    220223        """ 
     
    223226        """ 
    224227        return self._power_law(x) 
    225      
     228 
    226229    def evaluate_model_errors(self, x): 
    227230        """ 
     
    230233        """ 
    231234        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)\ 
     235        p2 = numpy.array([self.scale * self.power * math.pow(q, -self.power - 1)\ 
    233236                           * 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        diq2 = p1 * p1 + p2 * p2 
     238        return numpy.array([math.sqrt(err) for err in diq2]) 
     239 
    237240    def _power_law(self, x): 
    238241        """ 
    239242        F(x) = scale* (x)^(-power) 
    240             when power= 4. the model is porod  
     243            when power= 4. the model is porod 
    241244            else power_law 
    242245        The model has three parameters: :: 
     
    244247            2. power: power of the function 
    245248            3. scale : scale factor value 
    246          
     249 
    247250        :param x: array 
    248251        :return: F(x) 
     
    250253        if self.power <= 0: 
    251254            msg = "Power_law function expected positive power," 
    252             msg += " but got %s"%self.power 
     255            msg += " but got %s" % self.power 
    253256            raise ValueError(msg) 
    254257        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 ])   
     258            msg = "scale expected positive value, but got %s" % self.scale 
     259            raise ValueError(msg) 
     260 
     261        value = numpy.array([math.pow(i, -self.power) for i in x]) 
    259262        return self.scale * value 
    260263 
    261 class Extrapolator: 
     264class Extrapolator(object): 
    262265    """ 
    263266    Extrapolate I(q) distribution using a given model 
     
    266269        """ 
    267270        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  
     271 
     272        If a model is given, it will be used to linearize the data before 
    270273        the extrapolation is performed. If None, 
    271274        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 
     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 
    277280        self.model = model 
    278         
     281 
    279282        # Set qmin as the lowest non-zero value 
    280283        self.qmin = Q_MINIMUM 
    281284        for q_value in self.data.x: 
    282             if q_value > 0:  
     285            if q_value > 0: 
    283286                self.qmin = q_value 
    284287                break 
    285288        self.qmax = max(self.data.x) 
    286               
     289 
    287290    def fit(self, power=None, qmin=None, qmax=None): 
    288291        """ 
    289292        Fit data for y = ax + b  return a and b 
    290         
     293 
    291294        :param power: a fixed, otherwise None 
    292295        :param qmin: Minimum Q-value 
     
    297300        if qmax is None: 
    298301            qmax = self.qmax 
    299              
     302 
    300303        # Identify the bin range for the fit 
    301304        idx = (self.data.x >= qmin) & (self.data.x <= qmax) 
    302          
     305 
    303306        fx = numpy.zeros(len(self.data.x)) 
    304307 
    305308        # 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        if type(self.data.dy) == numpy.ndarray and \ 
     310            len(self.data.dy) == len(self.data.x) and \ 
     311            numpy.all(self.data.dy > 0): 
    309312            sigma = self.data.dy 
    310313        else: 
     
    313316        # Compute theory data f(x) 
    314317        fx[idx] = self.data.y[idx] 
    315          
     318 
    316319        # Linearize the data 
    317320        if self.model is not None: 
    318321            linearized_data = self.model.linearize_data(\ 
    319322                                            LoaderData1D(self.data.x[idx], 
    320                                                                 fx[idx], 
    321                                                             dy = sigma[idx])) 
     323                                                         fx[idx], 
     324                                                         dy=sigma[idx])) 
    322325        else: 
    323326            linearized_data = LoaderData1D(self.data.x[idx], 
    324327                                           fx[idx], 
    325                                            dy = sigma[idx]) 
    326          
    327         ##power is given only for function = power_law     
     328                                           dy=sigma[idx]) 
     329 
     330        ##power is given only for function = power_law 
    328331        if power != None: 
    329332            sigma2 = linearized_data.dy * linearized_data.dy 
    330333            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) 
     334            b = (numpy.sum(linearized_data.y / sigma2) \ 
     335                 - a * numpy.sum(linearized_data.x / sigma2)) / numpy.sum(1.0 / sigma2) 
     336 
     337 
     338            deltas = linearized_data.x * a + \ 
     339                    numpy.ones(len(linearized_data.x)) * b - linearized_data.y 
     340            residuals = numpy.sum(deltas * deltas / sigma2) 
     341 
     342            err = math.fabs(residuals) / numpy.sum(1.0 / sigma2) 
    340343            return [a, b], [0, math.sqrt(err)] 
    341344        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              
     345            A = numpy.vstack([linearized_data.x / linearized_data.dy, 1.0 / linearized_data.dy]).T 
     346            (p, residuals, _, _) = numpy.linalg.lstsq(A, linearized_data.y / linearized_data.dy) 
     347 
    347348            # Get the covariance matrix, defined as inv_cov = a_transposed * a 
    348349            err = numpy.zeros(2) 
     
    354355            except: 
    355356                err = [-1.0, -1.0] 
    356                  
     357 
    357358            return p, err 
    358          
     359 
    359360 
    360361class InvariantCalculator(object): 
     
    363364    Can provide volume fraction and surface area if the user provides 
    364365    Porod constant  and contrast values. 
    365      
     366 
    366367    :precondition:  the user must send a data of type DataLoader.Data1D 
    367368                    the user provide background and scale values. 
    368                      
    369     :note: Some computations depends on each others.  
     369 
     370    :note: Some computations depends on each others. 
    370371    """ 
    371     def __init__(self, data, background=0, scale=1 ): 
     372    def __init__(self, data, background=0, scale=1): 
    372373        """ 
    373374        Initialize variables. 
    374          
     375 
    375376        :param data: data must be of type DataLoader.Data1D 
    376         :param background: Background value. The data will be corrected  
     377        :param background: Background value. The data will be corrected 
    377378            before processing 
    378379        :param scale: Scaling factor for I(q). The data will be corrected 
     
    388389        self._data = self._get_data(data) 
    389390        # 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        if self._data.dxl != None and self._data.dxl.all() > 0: 
    391392            # assumes constant dxl 
    392393            self._smeared = self._data.dxl[0] 
    393        
     394 
    394395        # Since there are multiple variants of Q*, you should force the 
    395396        # user to use the get method and keep Q* a private data member 
    396397        self._qstar = None 
    397          
     398 
    398399        # You should keep the error on Q* so you can reuse it without 
    399400        # recomputing the whole thing. 
    400401        self._qstar_err = 0 
    401          
     402 
    402403        # Extrapolation parameters 
    403404        self._low_extrapolation_npts = 4 
     
    405406        self._low_extrapolation_power = None 
    406407        self._low_extrapolation_power_fitted = None 
    407     
     408 
    408409        self._high_extrapolation_npts = 4 
    409410        self._high_extrapolation_function = PowerLaw() 
    410411        self._high_extrapolation_power = None 
    411412        self._high_extrapolation_power_fitted = None 
    412          
     413 
    413414        # Extrapolation range 
    414415        self._low_q_limit = Q_MINIMUM 
    415          
     416 
    416417    def _get_data(self, data): 
    417418        """ 
    418419        :note: this function must be call before computing any type 
    419420         of invariant 
    420           
     421 
    421422        :return: new data = self._scale *data - self._background 
    422423        """ 
    423424        if not issubclass(data.__class__, LoaderData1D): 
    424425            #Process only data that inherited from DataLoader.Data_info.Data1D 
    425             raise ValueError,"Data must be of type DataLoader.Data1D" 
     426            raise ValueError, "Data must be of type DataLoader.Data1D" 
    426427        #from copy import deepcopy 
    427         new_data = (self._scale * data) - self._background    
    428          
     428        new_data = (self._scale * data) - self._background 
     429 
    429430        # Check that the vector lengths are equal 
    430         assert(len(new_data.x)==len(new_data.y)) 
    431          
     431        assert len(new_data.x) == len(new_data.y) 
     432 
    432433        # Verify that the errors are set correctly 
    433434        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))   
     435            (min(new_data.dy) == 0 and max(new_data.dy) == 0): 
     436            new_data.dy = numpy.ones(len(new_data.x)) 
    436437        return  new_data 
    437       
     438 
    438439    def _fit(self, model, qmin=Q_MINIMUM, qmax=Q_MAXIMUM, power=None): 
    439440        """ 
    440         fit data with function using  
     441        fit data with function using 
    441442        data = self._get_data() 
    442443        fx = Functor(data , function) 
    443444        y = data.y 
    444445        slope, constant = linalg.lstsq(y,fx) 
    445          
     446 
    446447        :param qmin: data first q value to consider during the fit 
    447448        :param qmax: data last q value to consider during the fit 
    448         :param power : power value to consider for power-law  
     449        :param power : power value to consider for power-law 
    449450        :param function: the function to use during the fit 
    450          
     451 
    451452        :return a: the scale of the function 
    452453        :return b: the other parameter of the function for guinier will be radius 
     
    454455        """ 
    455456        extrapolator = Extrapolator(data=self._data, model=model) 
    456         p, dp = extrapolator.fit(power=power, qmin=qmin, qmax=qmax)  
    457          
     457        p, dp = extrapolator.fit(power=power, qmin=qmin, qmax=qmax) 
     458 
    458459        return model.extract_model_parameters(constant=p[1], slope=p[0], 
    459460                                              dconstant=dp[1], dslope=dp[0]) 
    460      
     461 
    461462    def _get_qstar(self, data): 
    462463        """ 
    463464        Compute invariant for pinhole data. 
    464465        This invariant is given by: :: 
    465      
    466             q_star = x0**2 *y0 *dx0 +x1**2 *y1 *dx1  
     466 
     467            q_star = x0**2 *y0 *dx0 +x1**2 *y1 *dx1 
    467468                        + ..+ xn**2 *yn *dxn    for non smeared data 
    468                          
    469             q_star = dxl0 *x0 *y0 *dx0 +dxl1 *x1 *y1 *dx1  
     469 
     470            q_star = dxl0 *x0 *y0 *dx0 +dxl1 *x1 *y1 *dx1 
    470471                        + ..+ dlxn *xn *yn *dxn    for smeared data 
    471                          
     472 
    472473            where n >= len(data.x)-1 
    473474            dxl = slit height dQl 
     
    475476            dx0 = (x1 - x0)/2 
    476477            dxn = (xn - xn-1)/2 
    477              
     478 
    478479        :param data: the data to use to compute invariant. 
    479          
     480 
    480481        :return q_star: invariant value for pinhole data. q_star > 0 
    481482        """ 
    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)) 
     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)) 
    486486            raise ValueError, msg 
    487487        else: 
     
    490490                gx = data.x * data.x 
    491491            # assumes that len(x) == len(dxl). 
    492             else:                
    493                 gx = data.dxl * data.x  
    494                  
    495             n = len(data.x)- 1 
     492            else: 
     493                gx = data.dxl * data.x 
     494 
     495            n = len(data.x) - 1 
    496496            #compute the first delta q 
    497             dx0 = (data.x[1] - data.x[0])/2 
     497            dx0 = (data.x[1] - data.x[0]) / 2 
    498498            #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              
     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 
    504504            if len(data.x) == 2: 
    505                 return sum 
     505                return total 
    506506            else: 
    507507                #iterate between for element different 
    508508                #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              
     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 
    514514    def _get_qstar_uncertainty(self, data): 
    515515        """ 
    516516        Compute invariant uncertainty with with pinhole data. 
    517517        This uncertainty is given as follow: :: 
    518          
     518 
    519519           dq_star = math.sqrt[(x0**2*(dy0)*dx0)**2 + 
    520520                (x1**2 *(dy1)*dx1)**2 + ..+ (xn**2 *(dyn)*dxn)**2 ] 
     
    524524        dxn = (xn - xn-1)/2 
    525525        dyn: error on dy 
    526         
     526 
    527527        :param data: 
    528528        :note: if data doesn't contain dy assume dy= math.sqrt(data.y) 
    529         """           
     529        """ 
    530530        if len(data.x) <= 1 or len(data.y) <= 1 or \ 
    531531            len(data.x) != len(data.y) or \ 
    532532            (data.dy is not None and (len(data.dy) != len(data.y))): 
    533533            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)) 
     534            msg += " and greater than 1; got x=%s, y=%s" % (len(data.x), len(data.y)) 
    536535            raise ValueError, msg 
    537536        else: 
    538537            #Create error for data without dy error 
    539538            if data.dy is None: 
    540                 dy = math.sqrt(data.y)  
     539                dy = math.sqrt(data.y) 
    541540            else: 
    542541                dy = data.dy 
     
    547546            else: 
    548547                gx = data.dxl * data.x 
    549    
     548 
    550549            n = len(data.x) - 1 
    551550            #compute the first delta 
    552             dx0 = (data.x[1] - data.x[0])/2 
     551            dx0 = (data.x[1] - data.x[0]) / 2 
    553552            #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 
     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 
    558557            if len(data.x) == 2: 
    559                 return math.sqrt(sum) 
     558                return math.sqrt(total) 
    560559            else: 
    561560                #iterate between for element different 
    562561                #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          
     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 
    568567    def _get_extrapolated_data(self, model, npts=INTEGRATION_NSTEPS, 
    569                               q_start=Q_MINIMUM, q_end=Q_MAXIMUM): 
     568                               q_start=Q_MINIMUM, q_end=Q_MAXIMUM): 
    570569        """ 
    571570        :return: extrapolate data create from data 
     
    578577        iq = model.evaluate_model(q) 
    579578        diq = model.evaluate_model_errors(q) 
    580           
     579 
    581580        result_data = LoaderData1D(x=q, y=iq, dy=diq) 
    582581        if self._smeared != None: 
    583582            result_data.dxl = self._smeared * numpy.ones(len(q)) 
    584583        return result_data 
    585      
     584 
    586585    def get_data(self): 
    587586        """ 
     
    589588        """ 
    590589        return self._data 
    591      
     590 
    592591    def get_extrapolation_power(self, range='high'): 
    593592        """ 
     
    598597            return self._low_extrapolation_power_fitted 
    599598        return self._high_extrapolation_power_fitted 
    600      
     599 
    601600    def get_qstar_low(self): 
    602601        """ 
    603602        Compute the invariant for extrapolated data at low q range. 
    604          
     603 
    605604        Implementation: 
    606605            data = self._get_extra_data_low() 
    607606            return self._get_qstar() 
    608              
     607 
    609608        :return q_star: the invariant for data extrapolated at low q. 
    610609        """ 
     
    612611        qmin = self._data.x[0] 
    613612        qmax = self._data.x[self._low_extrapolation_npts - 1] 
    614          
     613 
    615614        # 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) 
     615        p, _ = self._fit(model=self._low_extrapolation_function, 
     616                         qmin=qmin, 
     617                         qmax=qmax, 
     618                         power=self._low_extrapolation_power) 
    620619        self._low_extrapolation_power_fitted = p[0] 
    621          
     620 
    622621        # Distribution starting point 
    623622        self._low_q_limit = Q_MINIMUM 
    624623        if Q_MINIMUM >= qmin: 
    625             self._low_q_limit = qmin/10 
    626          
     624            self._low_q_limit = qmin / 10 
     625 
    627626        data = self._get_extrapolated_data(\ 
    628627                                    model=self._low_extrapolation_function, 
    629                                             npts=INTEGRATION_NSTEPS, 
    630                                         q_start=self._low_q_limit, q_end=qmin) 
    631          
     628                                    npts=INTEGRATION_NSTEPS, 
     629                                    q_start=self._low_q_limit, q_end=qmin) 
     630 
    632631        # Systematic error 
    633632        # 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  
     633        # may not be a Guinier or simple power law. The following is 
    635634        # 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          
     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 
    640639    def get_qstar_high(self): 
    641640        """ 
    642641        Compute the invariant for extrapolated data at high q range. 
    643          
     642 
    644643        Implementation: 
    645644            data = self._get_extra_data_high() 
    646645            return self._get_qstar() 
    647              
     646 
    648647        :return q_star: the invariant for data extrapolated at high q. 
    649648        """ 
     
    652651        qmin = self._data.x[x_len - (self._high_extrapolation_npts - 1)] 
    653652        qmax = self._data.x[x_len] 
    654          
     653 
    655654        # 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) 
     655        p, _ = self._fit(model=self._high_extrapolation_function, 
     656                         qmin=qmin, 
     657                         qmax=qmax, 
     658                         power=self._high_extrapolation_power) 
    660659        self._high_extrapolation_power_fitted = p[0] 
    661          
     660 
    662661        #create new Data1D to compute the invariant 
    663662        data = self._get_extrapolated_data(\ 
    664663                                    model=self._high_extrapolation_function, 
    665                                            npts=INTEGRATION_NSTEPS, 
    666                                            q_start=qmax, q_end=Q_MAXIMUM)         
    667          
     664                                    npts=INTEGRATION_NSTEPS, 
     665                                    q_start=qmax, q_end=Q_MAXIMUM) 
     666 
    668667        return self._get_qstar(data), self._get_qstar_uncertainty(data) 
    669      
     668 
    670669    def get_extra_data_low(self, npts_in=None, q_start=None, npts=20): 
    671670        """ 
    672671        Returns the extrapolated data used for the loew-Q invariant calculation. 
    673         By default, the distribution will cover the data points used for the  
     672        By default, the distribution will cover the data points used for the 
    674673        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          
     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 
    679678        :param npts_in: number of data points for which 
    680679            the extrapolated data overlap 
    681680        :param q_start: is the minimum value to uses for extrapolated data 
    682681        :param npts: the number of points in the extrapolated distribution 
    683             
     682 
    684683        """ 
    685684        # Get extrapolation range 
    686685        if q_start is None: 
    687686            q_start = self._low_q_limit 
    688              
     687 
    689688        if npts_in is None: 
    690689            npts_in = self._low_extrapolation_npts 
    691         q_end = self._data.x[max(0, npts_in-1)] 
    692          
     690        q_end = self._data.x[max(0, npts_in - 1)] 
     691 
    693692        if q_start >= q_end: 
    694693            return numpy.zeros(0), numpy.zeros(0) 
     
    696695        return self._get_extrapolated_data(\ 
    697696                                    model=self._low_extrapolation_function, 
    698                                            npts=npts, 
    699                                            q_start=q_start, q_end=q_end) 
    700            
     697                                    npts=npts, 
     698                                    q_start=q_start, q_end=q_end) 
     699 
    701700    def get_extra_data_high(self, npts_in=None, q_end=Q_MAXIMUM, npts=20): 
    702701        """ 
    703702        Returns the extrapolated data used for the high-Q invariant calculation. 
    704         By default, the distribution will cover the data points used for the  
     703        By default, the distribution will cover the data points used for the 
    705704        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          
     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 
    710709        :param npts_in: number of data points for which the 
    711710            extrapolated data overlap 
     
    717716            npts_in = self._high_extrapolation_npts 
    718717        _npts = len(self._data.x) 
    719         q_start = self._data.x[min(_npts, _npts-npts_in)] 
    720          
     718        q_start = self._data.x[min(_npts, _npts - npts_in)] 
     719 
    721720        if q_start >= q_end: 
    722721            return numpy.zeros(0), numpy.zeros(0) 
    723          
     722 
    724723        return self._get_extrapolated_data(\ 
    725724                                model=self._high_extrapolation_function, 
    726                                            npts=npts, 
    727                                            q_start=q_start, q_end=q_end) 
    728       
     725                                npts=npts, 
     726                                q_start=q_start, q_end=q_end) 
     727 
    729728    def set_extrapolation(self, range, npts=4, function=None, power=None): 
    730729        """ 
    731730        Set the extrapolation parameters for the high or low Q-range. 
    732731        Note that this does not turn extrapolation on or off. 
    733          
     732 
    734733        :param range: a keyword set the type of extrapolation . type string 
    735734        :param npts: the numbers of q points of data to consider 
     
    739738            of type string. 
    740739        :param power: an power to apply power_low function 
    741                  
     740 
    742741        """ 
    743742        range = range.lower() 
     
    748747            msg = "Extrapolation function should be 'guinier' or 'power_law'" 
    749748            raise ValueError, msg 
    750          
     749 
    751750        if range == 'high': 
    752751            if function != 'power_law': 
    753752                msg = "Extrapolation only allows a power law at high Q" 
    754753                raise ValueError, msg 
    755             self._high_extrapolation_npts  = npts 
     754            self._high_extrapolation_npts = npts 
    756755            self._high_extrapolation_power = power 
    757756            self._high_extrapolation_power_fitted = power 
     
    761760            else: 
    762761                self._low_extrapolation_function = Guinier() 
    763             self._low_extrapolation_npts  = npts 
     762            self._low_extrapolation_npts = npts 
    764763            self._low_extrapolation_power = power 
    765764            self._low_extrapolation_power_fitted = power 
    766          
     765 
    767766    def get_qstar(self, extrapolation=None): 
    768767        """ 
    769768        Compute the invariant of the local copy of data. 
    770         
    771         :param extrapolation: string to apply optional extrapolation  
    772             
     769 
     770        :param extrapolation: string to apply optional extrapolation 
     771 
    773772        :return q_star: invariant of the data within data's q range 
    774          
     773 
    775774        :warning: When using setting data to Data1D , 
    776775            the user is responsible of 
    777776            checking that the scale and the background are 
    778777            properly apply to the data 
    779          
     778 
    780779        """ 
    781780        self._qstar = self._get_qstar(self._data) 
    782781        self._qstar_err = self._get_qstar_uncertainty(self._data) 
    783          
     782 
    784783        if extrapolation is None: 
    785784            return self._qstar 
    786          
     785 
    787786        # Compute invariant plus invariant of extrapolated data 
    788         extrapolation = extrapolation.lower()     
     787        extrapolation = extrapolation.lower() 
    789788        if extrapolation == "low": 
    790789            qs_low, dqs_low = self.get_qstar_low() 
    791             qs_hi, dqs_hi   = 0, 0 
    792              
     790            qs_hi, dqs_hi = 0, 0 
     791 
    793792        elif extrapolation == "high": 
    794793            qs_low, dqs_low = 0, 0 
    795             qs_hi, dqs_hi   = self.get_qstar_high() 
    796              
     794            qs_hi, dqs_hi = self.get_qstar_high() 
     795 
    797796        elif extrapolation == "both": 
    798797            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          
     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 
    805804        return self._qstar 
    806         
     805 
    807806    def get_surface(self, contrast, porod_const, extrapolation=None): 
    808807        """ 
    809808        Compute the specific surface from the data. 
    810          
     809 
    811810        Implementation:: 
    812          
     811 
    813812          V =  self.get_volume_fraction(contrast, extrapolation) 
    814      
     813 
    815814          Compute the surface given by: 
    816815            surface = (2*pi *V(1- V)*porod_const)/ q_star 
    817             
     816 
    818817        :param contrast: contrast value to compute the volume 
    819         :param porod_const: Porod constant to compute the surface  
     818        :param porod_const: Porod constant to compute the surface 
    820819        :param extrapolation: string to apply optional extrapolation 
    821          
    822         :return: specific surface  
     820 
     821        :return: specific surface 
    823822        """ 
    824823        # Compute the volume 
    825824        volume = self.get_volume_fraction(contrast, extrapolation) 
    826         return 2 * math.pi * volume *(1 - volume) * \ 
    827             float(porod_const)/self._qstar 
    828          
     825        return 2 * math.pi * volume * (1 - volume) * \ 
     826            float(porod_const) / self._qstar 
     827 
    829828    def get_volume_fraction(self, contrast, extrapolation=None): 
    830829        """ 
    831830        Compute volume fraction is deduced as follow: :: 
    832          
     831 
    833832            q_star = 2*(pi*contrast)**2* volume( 1- volume) 
    834833            for k = 10^(-8)*q_star/(2*(pi*|contrast|)**2) 
     
    837836                 volume1 = (1- sqrt(1- 4*k))/2 
    838837                 volume2 = (1+ sqrt(1- 4*k))/2 
    839             
     838 
    840839            q_star: the invariant value included extrapolation is applied 
    841840                         unit  1/A^(3)*1/cm 
    842841                    q_star = self.get_qstar() 
    843                      
     842 
    844843            the result returned will be 0 <= volume <= 1 
    845          
     844 
    846845        :param contrast: contrast value provides by the user of type float. 
    847846                 contrast unit is 1/A^(2)= 10^(16)cm^(2) 
    848847        :param extrapolation: string to apply optional extrapolation 
    849          
     848 
    850849        :return: volume fraction 
    851          
     850 
    852851        :note: volume fraction must have no unit 
    853852        """ 
    854853        if contrast <= 0: 
    855             raise ValueError, "The contrast parameter must be greater than zero"   
    856          
     854            raise ValueError, "The contrast parameter must be greater than zero" 
     855 
    857856        # Make sure Q star is up to date 
    858857        self.get_qstar(extrapolation) 
    859          
     858 
    860859        if self._qstar <= 0: 
    861860            msg = "Invalid invariant: Invariant Q* must be greater than zero" 
    862861            raise RuntimeError, msg 
    863          
     862 
    864863        # Compute intermediate constant 
    865         k =  1.e-8 * self._qstar/(2 * (math.pi * math.fabs(float(contrast)))**2) 
     864        k = 1.e-8 * self._qstar / (2 * (math.pi * math.fabs(float(contrast))) ** 2) 
    866865        # Check discriminant value 
    867866        discrim = 1 - 4 * k 
    868          
     867 
    869868        # Compute volume fraction 
    870869        if discrim < 0: 
     
    872871            raise RuntimeError, msg 
    873872        elif discrim == 0: 
    874             return 1/2 
     873            return 1 / 2 
    875874        else: 
    876875            volume1 = 0.5 * (1 - math.sqrt(discrim)) 
    877876            volume2 = 0.5 * (1 + math.sqrt(discrim)) 
    878              
     877 
    879878            if 0 <= volume1 and volume1 <= 1: 
    880879                return volume1 
    881             elif 0 <= volume2 and volume2 <= 1:  
    882                 return volume2  
     880            elif 0 <= volume2 and volume2 <= 1: 
     881                return volume2 
    883882            msg = "Could not compute the volume fraction: inconsistent results" 
    884883            raise RuntimeError, msg 
    885      
     884 
    886885    def get_qstar_with_error(self, extrapolation=None): 
    887886        """ 
     
    889888        This uncertainty computation depends on whether or not the data is 
    890889        smeared. 
    891          
     890 
    892891        :param extrapolation: string to apply optional extrapolation 
    893          
     892 
    894893        :return: invariant, the invariant uncertainty 
    895         """     
     894        """ 
    896895        self.get_qstar(extrapolation) 
    897896        return self._qstar, self._qstar_err 
    898      
     897 
    899898    def get_volume_fraction_with_error(self, contrast, extrapolation=None): 
    900899        """ 
    901900        Compute uncertainty on volume value as well as the volume fraction 
    902901        This uncertainty is given by the following equation: :: 
    903          
     902 
    904903            dV = 0.5 * (4*k* dq_star) /(2* math.sqrt(1-k* q_star)) 
    905                                   
     904 
    906905            for k = 10^(-8)*q_star/(2*(pi*|contrast|)**2) 
    907              
     906 
    908907            q_star: the invariant value including extrapolated value if existing 
    909908            dq_star: the invariant uncertainty 
    910909            dV: the volume uncertainty 
    911          
     910 
    912911        The uncertainty will be set to -1 if it can't be computed. 
    913          
    914         :param contrast: contrast value  
     912 
     913        :param contrast: contrast value 
    915914        :param extrapolation: string to apply optional extrapolation 
    916          
     915 
    917916        :return: V, dV = volume fraction, error on volume fraction 
    918917        """ 
    919918        volume = self.get_volume_fraction(contrast, extrapolation) 
    920          
     919 
    921920        # Compute error 
    922         k =  1.e-8 * self._qstar /(2 * (math.pi* math.fabs(float(contrast)))**2) 
     921        k = 1.e-8 * self._qstar / (2 * (math.pi * math.fabs(float(contrast))) ** 2) 
    923922        # Check value inside the sqrt function 
    924923        value = 1 - k * self._qstar 
     
    927926        # Compute uncertainty 
    928927        uncertainty = math.fabs((0.5 * 4 * k * \ 
    929                         self._qstar_err)/(2 * math.sqrt(1 - k * self._qstar))) 
    930          
     928                        self._qstar_err) / (2 * math.sqrt(1 - k * self._qstar))) 
     929 
    931930        return volume, uncertainty 
    932      
     931 
    933932    def get_surface_with_error(self, contrast, porod_const, extrapolation=None): 
    934933        """ 
    935934        Compute uncertainty of the surface value as well as the surface value. 
    936935        The uncertainty is given as follow: :: 
    937          
     936 
    938937            dS = porod_const *2*pi[( dV -2*V*dV)/q_star 
    939938                 + dq_star(v-v**2) 
    940                   
     939 
    941940            q_star: the invariant value 
    942941            dq_star: the invariant uncertainty 
    943942            V: the volume fraction value 
    944943            dV: the volume uncertainty 
    945          
     944 
    946945        :param contrast: contrast value 
    947         :param porod_const: porod constant value  
     946        :param porod_const: porod constant value 
    948947        :param extrapolation: string to apply optional extrapolation 
    949          
     948 
    950949        :return S, dS: the surface, with its uncertainty 
    951950        """ 
     
    956955        v, dv = self.get_volume_fraction_with_error(contrast, extrapolation) 
    957956 
    958         s = self.get_surface(contrast=contrast, porod_const=porod_const,  
     957        s = self.get_surface(contrast=contrast, porod_const=porod_const, 
    959958                             extrapolation=extrapolation) 
    960         ds = porod_const * 2 * math.pi * (( dv - 2 * v * dv)/ self._qstar\ 
    961                  + self._qstar_err * ( v - v**2)) 
     959        ds = porod_const * 2 * math.pi * ((dv - 2 * v * dv) / self._qstar\ 
     960                 + self._qstar_err * (v - v ** 2)) 
    962961 
    963962        return s, ds 
  • src/sas/pr/invertor.py

    r5f8fc78 rc76bdc3  
    110110    slit_height = None 
    111111    slit_width = None 
     112    ## Error 
     113    err = 0 
     114    ## Data to invert 
     115    x = None 
     116    y = None 
    112117 
    113118    def __init__(self): 
    114119        Cinvertor.__init__(self) 
    115         self.err = None 
    116         self.d_max = None 
     120        self.d_max = 200.0 
    117121        self.q_min = self.set_qmin(-1.0) 
    118122        self.q_max = self.set_qmax(-1.0) 
    119         self.x = None 
    120         self.y = None 
    121123        self.has_bck = False 
    122124 
Note: See TracChangeset for help on using the changeset viewer.