Changeset 59a41066 in sasview for Invariant


Ignore:
Timestamp:
Jan 14, 2010 2:31:24 PM (15 years ago)
Author:
Gervaise Alina <gervyh@…>
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:
82703a1
Parents:
3082632
Message:

refactoring [incomplete] not all test passing

File:
1 edited

Legend:

Unmodified
Added
Removed

  • Invariant/invariant.py

    r3082632 r59a41066  
    2020INTEGRATION_NSTEPS = 1000 
    2121 
    22 def guinier(x, scale=1, radius=60): 
     22class Transform(object): 
    2323    """ 
    24         Compute a F(x) = scale* e-((radius*x)**2/3). 
    25         @param x: a vector of q values  
    26         @param scale: the scale value 
    27         @param radius: the guinier radius value 
    28         @return F(x) 
    29     """    
    30     if radius <= 0: 
    31         raise ValueError("Rg expected positive value, but got %s"%radius)   
    32     value = numpy.array([math.exp(-((radius * i)**2/3)) for i in x ])  
    33     scale = numpy.sqrt(scale*scale)       
    34     if scale == 0: 
    35         raise ValueError("scale expected positive value, but got %s"%scale)  
    36     return scale * value 
    37  
    38 def power_law(x, scale=1, power=4): 
     24        Define interface that need to compute a function or an inverse 
     25        function given some x, y  
    3926    """ 
    40         F(x) = scale* (x)^(-power) 
    41             when power= 4. the model is porod  
    42             else power_law 
    43         The model has three parameters:  
    44         @param x: a vector of q values 
    45         @param power: power of the function 
    46         @param scale : scale factor value 
    47         @param F(x) 
    48     """   
    49     if power <=0: 
    50         raise ValueError("Power_law function expected positive power, but got %s"%power) 
    51      
    52     value = numpy.array([ math.pow(i, -power) for i in x ])   
    53     scale = numpy.sqrt(scale*scale) 
    54     if scale == 0: 
    55         raise ValueError("scale expected positive value, but got %s"%scale)  
    56     return scale * value 
    57  
    58 class FitFunctor: 
     27    def set_value(self, a, b): 
     28        """ 
     29            set private member 
     30        """ 
     31    def transform_value(self, value): 
     32        """ 
     33            Can use one of the function transform_x2 or transform_logx 
     34        """ 
     35    def transform_x2(self, value): 
     36        """ 
     37            @param value: float 
     38        """ 
     39        return value **2 
     40     
     41    def transform_to_logx(self, value): 
     42        """ 
     43            @param value: float 
     44        """ 
     45        return math.log(value) 
     46         
     47    def inverse_transform_value(self, value): 
     48        """ 
     49            reverse transform vaalue given a specific function 
     50            return f-1(value) 
     51            @param value : float 
     52        """ 
     53    def get_extrapolated_data(self, data, q_start=Q_MINIMUM, q_end=Q_MAXIMUM): 
     54        """ 
     55            @return extrapolate data create from data 
     56        """ 
     57    def transform_data(self, x): 
     58        """ 
     59            @param x: vector of float 
     60            @param y : vector of float 
     61            return x, y transform given a specific function 
     62        """ 
     63    def inverse_transform_data(self, x, y): 
     64        """ 
     65            reverse transform x, y given a specific function 
     66        """ 
     67    def transform_error(self, y , dy=None): 
     68        if dy is None: 
     69            dy = numpy.array([math.sqrt(j) for j in y]) 
     70        return dy / y 
     71     
     72    def get_extrapolated_data(self, data, npts=INTEGRATION_NSTEPS, 
     73                              q_start=Q_MINIMUM, q_end=Q_MAXIMUM, smear_indice=0): 
     74        """ 
     75            @return extrapolate data create from data 
     76        """ 
     77        #create new Data1D to compute the invariant 
     78        new_x = numpy.linspace(start=q_start, 
     79                               stop=q_end, 
     80                               num=npts, 
     81                               endpoint=True) 
     82        new_y = self.transform_data(x=new_x) 
     83          
     84        dxl = None 
     85        dxw = None 
     86        if data.dxl is not None : 
     87            dxl = numpy.ones(INTEGRATION_NSTEPS) 
     88            dxl = dxl * data.dxl[smear_indice] 
     89            
     90        if data.dxw is not None : 
     91            dxw = numpy.ones(INTEGRATION_NSTEPS) 
     92            dxw = dxw * data.dxw[smear_indice] 
     93          
     94        result_data = LoaderData1D(x=new_x, y=new_y) 
     95        data.clone_without_data(clone=result_data) 
     96        result_data.dxl = dxl 
     97        result_data.dxw = dxw 
     98         
     99        
     100        return result_data 
     101   
     102class Guinier(Transform): 
     103    """ 
     104        class of type Transform that performs operations related to guinier  
     105        function 
     106    """ 
     107    def __init__(self, scale=1, radius=60): 
     108        Transform.__init__(self) 
     109        self.scale = scale 
     110        self.radius = radius 
     111     
     112    def transform_value(self, value): 
     113        return self.transform_x2(value) 
     114     
     115    def set_value(self, a, b): 
     116        # b is the radius value of the guinier function 
     117        if b >=0 : 
     118                raise ValueError("Guinier fit was not converged") 
     119        else: 
     120                b = math.sqrt(-3 * b) 
     121        a = math.exp(a) 
     122     
     123        self.scale = a 
     124        self.radius = b 
     125         
     126    def transform_data(self, x): 
     127        """ 
     128            given a scale and a radius transform x, y using a guinier 
     129            function 
     130        """ 
     131        return self._guinier(x) 
     132      
     133    def inverse_transform_data(self, x, y): 
     134        """ 
     135            given a scale and a radius transform x, y using a inverse guinier 
     136            function 
     137        """ 
     138        result_x = numpy.array([i * i for i in x])    
     139        result_y = numpy.array([math.log(j) for j in y]) 
     140         
     141        return result_x, result_y  
     142         
     143    def _guinier(self, x): 
     144        """ 
     145            Retrive the guinier function after apply an inverse guinier function 
     146            to x 
     147            Compute a F(x) = scale* e-((radius*x)**2/3). 
     148            @param x: a vector of q values  
     149            @param scale: the scale value 
     150            @param radius: the guinier radius value 
     151            @return F(x) 
     152        """    
     153        # transform the radius of coming from the inverse guinier function to a  
     154        # a radius of a guinier function 
     155        if self.radius <= 0: 
     156            raise ValueError("Rg expected positive value, but got %s"%self.radius)   
     157 
     158        value = numpy.array([math.exp(-((self.radius * i)**2/3)) for i in x ])  
     159        
     160        return self.scale * value 
     161 
     162class PowerLaw(Transform): 
     163    """ 
     164        class of type transform that perform operation related to power_law  
     165        function 
     166    """ 
     167    def __init__(self, scale=1, power=4): 
     168        Transform.__init__(self) 
     169        self.scale = scale 
     170        self.power = power 
     171         
     172    def set_value(self, a, b): 
     173        b = -1 * b 
     174        if b <= 0: 
     175                raise ValueError("Power_law fit expected posive power, but got %s"%power) 
     176        a = math.exp(a) 
     177        self.scale = a 
     178        self.power = b 
     179         
     180    def transform_data(self, x, a=None, b=None): 
     181        """ 
     182            given a scale and a radius transform x, y using a power_law 
     183            function 
     184        """ 
     185        if  a is not None: 
     186            self.scale = a 
     187        if  b is not None: 
     188            self.power = b 
     189             
     190        return self._power_law(x) 
     191        
     192    def inverse_transform_data(self, x, y): 
     193        """ 
     194            given a scale and a radius transform x, y using a inverse power_law 
     195            function 
     196        """ 
     197        result_x = numpy.array([i * i for i in x])    
     198        result_y = numpy.array([math.log(j) for j in y]) 
     199         
     200        return result_x, result_y  
     201     
     202    def _power_law(self, x): 
     203        """ 
     204            F(x) = scale* (x)^(-power) 
     205                when power= 4. the model is porod  
     206                else power_law 
     207            The model has three parameters:  
     208            @param x: a vector of q values 
     209            @param power: power of the function 
     210            @param scale : scale factor value 
     211            @param F(x) 
     212        """ 
     213        if self.power <= 0: 
     214            raise ValueError("Power_law function expected positive power, but got %s"%power) 
     215        if self.scale <= 0: 
     216            raise ValueError("scale expected positive value, but got %s"%self.scale)  
     217        
     218        value = numpy.array([ math.pow(i, -self.power) for i in x ])   
     219        return self.scale * value 
     220 
     221class Extrapolator: 
    59222    """ 
    60223        compute f(x) 
    … …  
    66229        """ 
    67230        self.data  = data 
    68          
     231        
    69232        # Set qmin as the lowest non-zero value 
    70233        self.qmin = Q_MINIMUM 
    71234        for q_value in self.data.x: 
    72             if q_value>0:  
     235            if q_value > 0:  
    73236                self.qmin = q_value 
    74237                break 
    … …  
    76239        
    77240        #get the smear object of data 
    78         self.smearer = smear_selection( self.data ) 
     241        self.smearer = smear_selection(self.data) 
    79242        # Set the q-range information to allow smearing 
    80243        self.set_fit_range() 
    … …  
    91254        self._qmin_unsmeared = self.qmin 
    92255        self._qmax_unsmeared = self.qmax     
     256         
    93257        if self.smearer is not None: 
    94258            self._first_unsmeared_bin, self._last_unsmeared_bin = self.smearer.get_bin_range(self.qmin, self.qmax) 
    … …  
    107271        """ 
    108272        fx = numpy.zeros(len(self.data.x)) 
    109         # Uncertainty 
     273 
     274         # Uncertainty 
    110275        if type(self.data.dy)==numpy.ndarray and len(self.data.dy)==len(self.data.x): 
    111276            sigma = self.data.dy 
    112277        else: 
    113278            sigma = numpy.ones(len(self.data.x)) 
    114  
     279             
    115280        # Compute theory data f(x) 
    116281        fx[self.idx_unsmeared] = self.data.y[self.idx_unsmeared] 
    … …  
    120285         
    121286        fx[self.idx_unsmeared] = fx[self.idx_unsmeared]/sigma[self.idx_unsmeared] 
    122  
     287         
    123288        ##power is given only for function = power_law     
    124289        if power != None: 
    … …  
    171336        # Extrapolation parameters 
    172337        self._low_extrapolation_npts = 4 
    173         self._low_extrapolation_function = guinier 
     338        self._low_extrapolation_function = Guinier() 
    174339        self._low_extrapolation_power = None 
    175340    
    176341        self._high_extrapolation_npts = 4 
    177         self._high_extrapolation_function = power_law 
     342        self._high_extrapolation_function = PowerLaw() 
    178343        self._high_extrapolation_power = None 
    179344         
    … …  
    192357        # Copy data that is not copied by the operations 
    193358        #TODO: fix this in DataLoader 
    194         new_data.dxl = data.dxl 
    195         new_data.dxw = data.dxw         
     359        #new_data.dxl = data.dxl 
     360        #new_data.dxw = data.dxw         
    196361 
    197362        return new_data 
    … …  
    212377                    for power_law will be the power value 
    213378        """ 
    214         #fit_x = numpy.array([math.log(x) for x in self._data.x]) 
    215         if function.__name__ == "guinier":    
    216             fit_x = numpy.array([x * x for x in self._data.x])      
    217             qmin = qmin**2 
    218             qmax = qmax**2 
    219             fit_y = numpy.array([math.log(y) for y in self._data.y]) 
    220             fit_dy = numpy.array([y for y in self._data.y]) 
    221             fit_dy = numpy.array([dy for dy in self._data.dy])/fit_dy 
    222  
    223         elif function.__name__ == "power_law": 
    224             fit_x = numpy.array([math.log(x) for x in self._data.x]) 
    225             qmin = math.log(qmin) 
    226             qmax = math.log(qmax) 
    227             fit_y = numpy.array([math.log(y) for y in self._data.y]) 
    228             fit_dy = numpy.array([y for y in self._data.y]) 
    229             fit_dy = numpy.array([dy for dy in self._data.dy])/fit_dy 
    230  
    231         else: 
    232             raise ValueError("Unknown function used to fit %s"%function.__name__) 
    233         
    234          
    235         #else: 
     379        transform = function 
     380        fit_x, fit_y = transform.inverse_transform_data(self._data.x, self._data.y) 
     381        fit_dy = transform.transform_error(y=self._data.y, dy=self._data.dy) 
     382        qmin = transform.transform_value(value=qmin) 
     383        qmax = transform.transform_value(value=qmax) 
     384         
    236385        fit_data = LoaderData1D(x=fit_x, y=fit_y, dy=fit_dy) 
    237386        fit_data.dxl = self._data.dxl 
    238387        fit_data.dxw = self._data.dxw    
    239         functor = FitFunctor(data=fit_data) 
     388        functor = Extrapolator(data=fit_data) 
    240389        functor.set_fit_range(qmin=qmin, qmax=qmax) 
    241390        b, a = functor.fit(power=power)          
    242        
    243                    
    244         if function.__name__ == "guinier": 
    245             # b is the radius value of the guinier function 
    246             if b>=0: 
    247                 raise ValueError("Guinier fit was not converged") 
    248             else: 
    249                 b = math.sqrt(-3 * b) 
    250  
    251  
    252         if function.__name__ == "power_law": 
    253             b = -1 * b 
    254             if b <= 0: 
    255                 raise ValueError("Power_law fit expected posive power, but got %s"%power) 
    256         # a is the scale of the guinier function 
    257         a = math.exp(a) 
    258  
    259         return a, b 
     391        return b, a 
    260392     
    261393    def _get_qstar(self, data): 
    … …  
    332464              or len(data.x)!= len(data.dxl): 
    333465            msg = "x, dxl, and y must be have the same length and greater than 1" 
     466            msg +=" len x=%s, y=%s, dxl=%s"%(len(data.x),len(data.y),len(data.dxl)) 
    334467            raise ValueError, msg 
    335468        else: 
    … …  
    438571            note: if data doesn't contain dy assume dy= math.sqrt(data.y) 
    439572        """ 
    440         if data is None: 
    441             data = self._data 
     573        #if data is None: 
     574        #    data = self._data 
    442575             
    443576        if not data.is_slit_smeared(): 
    … …  
    529662        # Extrapolate the low-Q data 
    530663        #TODO: this fit fails. Fix it. 
    531         a, b = self._fit(function=self._low_extrapolation_function, 
     664        b, a = self._fit(function=self._low_extrapolation_function, 
    532665                          qmin=qmin, 
    533666                          qmax=qmax, 
    … …  
    539672            q_start = qmin/10 
    540673             
    541         #create new Data1D to compute the invariant 
    542         new_x = numpy.linspace(start=q_start, 
    543                                stop=qmin, 
    544                                num=INTEGRATION_NSTEPS, 
    545                                endpoint=True) 
    546         new_y = self._low_extrapolation_function(x=new_x, scale=a, radius=b) 
    547         dxl = None 
    548         dxw = None 
    549         if self._data.dxl is not None: 
    550             dxl = numpy.ones(INTEGRATION_NSTEPS) 
    551             dxl = dxl * self._data.dxl[0] 
    552         if self._data.dxw is not None: 
    553             dxw = numpy.ones(INTEGRATION_NSTEPS) 
    554             dxw = dxw * self._data.dxw[0] 
    555  
    556         data_min = LoaderData1D(x=new_x, y=new_y) 
    557         data_min.dxl = dxl 
    558         data_min.dxw = dxw 
    559         self._data.clone_without_data( clone= data_min) 
    560  
     674        self._low_extrapolation_function.set_value(a= a, b=b) 
     675        data_min = self._low_extrapolation_function.get_extrapolated_data(data=self._data,  
     676                                            npts=INTEGRATION_NSTEPS, 
     677                              q_start=q_start, q_end=qmin, smear_indice=0) 
     678        
    561679        return data_min 
    562680           
    … …  
    577695        # Data boundaries for fiiting 
    578696        x_len = len(self._data.x) - 1 
    579         qmin = self._data.x[x_len - (self._high_extrapolation_npts - 1)] 
     697        q_start = self._data.x[x_len - (self._high_extrapolation_npts - 1)] 
    580698        qmax = self._data.x[x_len] 
    581          
    582         try: 
    583             # fit the data with a model to get the appropriate parameters 
    584             a, b = self._fit(function=self._high_extrapolation_function, 
    585                                    qmin=qmin, 
    586                                     qmax=qmax, 
    587                                     power=self._high_extrapolation_power) 
    588         except: 
    589             return None 
    590         
     699        smear_indice = x_len 
     700         
     701        # fit the data with a model to get the appropriate parameters 
     702        b, a = self._fit(function=self._high_extrapolation_function, 
     703                               qmin=q_start, 
     704                                qmax=qmax, 
     705                                power=self._high_extrapolation_power) 
     706   
    591707        #create new Data1D to compute the invariant 
    592         new_x = numpy.linspace(start=qmax, 
    593                                stop=Q_MAXIMUM, 
    594                                num=INTEGRATION_NSTEPS, 
    595                                endpoint=True) 
    596         
    597         new_y = self._high_extrapolation_function(x=new_x, scale=a, power=b) 
    598          
    599         dxl = None 
    600         dxw = None 
    601         if self._data.dxl is not None: 
    602             dxl = numpy.ones(INTEGRATION_NSTEPS) 
    603             dxl = dxl * self._data.dxl[0] 
    604         if self._data.dxw is not None: 
    605             dxw = numpy.ones(INTEGRATION_NSTEPS) 
    606             dxw = dxw * self._data.dxw[0] 
    607              
    608         data_max = LoaderData1D(x=new_x, y=new_y) 
    609         data_max.dxl = dxl 
    610         data_max.dxw = dxw 
    611         self._data.clone_without_data(clone=data_max) 
    612      
     708        self._high_extrapolation_function.set_value(a=a, b=b) 
     709        data_max = self._high_extrapolation_function.get_extrapolated_data(data=self._data,  
     710                                            npts=INTEGRATION_NSTEPS, 
     711                                            q_start=q_start, q_end=qmax, 
     712                                            smear_indice=smear_indice) 
    613713        return data_max 
    614714      
    … …  
    638738        else: 
    639739            if function == 'power_law': 
    640                 self._low_extrapolation_function = power_law 
     740                self._low_extrapolation_function = PowerLaw() 
    641741            else: 
    642                 self._low_extrapolation_function = guinier 
     742                self._low_extrapolation_function = Guinier() 
    643743            self._low_extrapolation_npts  = npts 
    644744            self._low_extrapolation_power = power 
Note: See TracChangeset for help on using the changeset viewer.