Changeset aafa962 in sasview


Ignore:
Timestamp:
Jan 16, 2010 12:44:34 PM (14 years ago)
Author:
Mathieu Doucet <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:
52070a1
Parents:
82703a1
Message:

invariant: minor code refactor

Location:
Invariant
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • Invariant/invariant.py

    r82703a1 raafa962  
    33    @author: Gervaise B. Alina/UTK 
    44""" 
    5  
     5#TODO: Need to decide if/how to use smearing for extrapolation 
    66import math  
    77import numpy 
    … …  
    2424        function given some x, y  
    2525    """ 
    26     def transform_values(self, a, b): 
     26    def linearize_data(self, x, y, dy=None): 
     27        """ 
     28            Linearize data so that a linear fit can be performed.  
     29            Filter out the data that can't be transformed. 
     30            @param x: array of q values 
     31            @param y: array of I(q) values 
     32            @param dy: array of dI(q) values 
     33        """ 
     34        return NotImplemented 
     35     
     36    def linearize_q_value(self, value): 
     37        """ 
     38            Transform the input q-value for linearization 
     39        """ 
     40        return NotImplemented 
     41 
     42    def extract_model_parameters(self, a, b): 
    2743        """ 
    2844            set private member 
    29         """ 
    30         return NotImplemented 
    31      
    32     def transform_value(self, value): 
    33         """ 
    34             @param value: float 
    35             return f(value) 
    3645        """ 
    3746        return NotImplemented 
    3847      
    39     def inverse_transform_value(self, value): 
    40         """ 
    41             reverse transform vaalue given a specific function 
    42             return f-1(value) 
    43             @param value : float 
     48    def evaluate_model(self, x): 
     49        """ 
     50            Returns an array f(x) values where f is the Transform function. 
    4451        """ 
    4552        return NotImplemented 
    4653     
    47     def get_extrapolated_data(self, data, q_start=Q_MINIMUM, q_end=Q_MAXIMUM): 
    48         """ 
    49             @return extrapolate data create from data 
    50         """ 
    51         return NotImplemented 
    52      
    53     def transform_data(self, x): 
    54         """ 
    55             @param x: vector of float 
    56             @param y : vector of float 
    57             return x, y transform given a specific function 
    58         """ 
    59         return NotImplemented 
    60      
    61     def inverse_transform_data(self, x, y): 
    62         """ 
    63             reverse transform x, y given a specific function 
    64         """ 
    65         return NotImplemented 
    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         return result_data 
    99    
    10054class Guinier(Transform): 
    10155    """ 
    … …  
    10862        self.radius = radius 
    10963     
    110     def transform_value(self, value): 
    111         """ 
    112             Square input value 
    113             @param value: of type float 
     64    def linearize_data(self, x, y, dy=None): 
     65        """ 
     66            Linearize data so that a linear fit can be performed.  
     67            Filter out the data that can't be transformed. 
     68            @param x: array of q values 
     69            @param y: array of I(q) values 
     70            @param dy: array of dI(q) values 
     71        """ 
     72        # Check that the vector lengths are equal 
     73        assert(len(x)==len(y)) 
     74        if dy is not None: 
     75            assert(len(x)==len(dy)) 
     76        else: 
     77            dy = numpy.array([math.sqrt(math.fabs(j)) for j in y]) 
     78         
     79        data_points = zip(x,y,dy) 
     80         
     81        # Transform the data 
     82        output_points = [(self.linearize_q_value(p[0]), 
     83                          math.log(p[1]), 
     84                          p[2]/p[1]) for p in data_points if p[0]>0 and p[1]>0] 
     85         
     86        x_out, y_out, dy_out = zip(*output_points) 
     87 
     88        return numpy.asarray(x_out), numpy.asarray(y_out), numpy.asarray(dy_out) 
     89     
     90    def linearize_q_value(self, value): 
     91        """ 
     92            Transform the input q-value for linearization 
     93            @param value: q-value 
     94            @return: q*q 
    11495        """ 
    11596        return value * value 
    11697     
    117     def transform_values(self, a, b): 
     98    def extract_model_parameters(self, a, b): 
    11899        """ 
    119100           assign new value to the scale and the radius 
    … …  
    125106        return a, b 
    126107         
    127     def transform_data(self, x): 
     108    def evaluate_model(self, x): 
    128109        """ 
    129110            return F(x)= scale* e-((radius*x)**2/3) 
    130111        """ 
    131112        return self._guinier(x) 
    132       
    133     def inverse_transform_data(self, x, y): 
    134         """ 
    135             this function finds x, y given equation1: y = ax + b 
    136             and equation2: y1 = scale* e-((radius*x1)**2/3). 
    137             where equation1, equation2 are equivalent 
    138             imply:  x = x1*x1 
    139                     y = ln(y1) 
    140                     b = math.exp(scale) 
    141                     a = -radius**2/3 
    142             @return  x, y 
    143         """ 
    144         result_x = numpy.array([i * i for i in x])    
    145         result_y = numpy.array([math.log(j) for j in y]) 
    146         return result_x, result_y  
    147          
     113              
    148114    def _guinier(self, x): 
    149115        """ 
    … …  
    173139        self.power = power 
    174140         
    175     def transform_values(self, a, b): 
     141    def linearize_data(self, x, y, dy=None): 
     142        """ 
     143            Linearize data so that a linear fit can be performed.  
     144            Filter out the data that can't be transformed. 
     145            @param x: array of q values 
     146            @param y: array of I(q) values 
     147            @param dy: array of dI(q) values 
     148        """ 
     149        # Check that the vector lengths are equal 
     150        assert(len(x)==len(y)) 
     151        if dy is not None: 
     152            assert(len(x)==len(dy)) 
     153        else: 
     154            dy = numpy.array([math.sqrt(math.fabs(j)) for j in y]) 
     155         
     156        data_points = zip(x,y,dy) 
     157         
     158        # Transform the data 
     159        output_points = [(self.linearize_q_value(p[0]), 
     160                          math.log(p[1]), 
     161                          p[2]/p[1]) for p in data_points if p[0]>0 and p[1]>0] 
     162         
     163        x_out, y_out, dy_out = zip(*output_points) 
     164 
     165        return numpy.asarray(x_out), numpy.asarray(y_out), numpy.asarray(dy_out) 
     166         
     167    def linearize_q_value(self, value): 
     168        """ 
     169            Transform the input q-value for linearization 
     170            @param value: q-value 
     171            @return: log(q) 
     172        """ 
     173        return math.log(value) 
     174     
     175    def extract_model_parameters(self, a, b): 
    176176        """ 
    177177            Assign new value to the scale and the power  
    … …  
    183183        return a, b 
    184184         
    185     def transform_value(self, value): 
    186         """ 
    187             compute log(value) 
    188         """ 
    189         return math.log(value) 
    190          
    191     def transform_data(self, x): 
     185    def evaluate_model(self, x): 
    192186        """ 
    193187            given a scale and a radius transform x, y using a power_law 
    … …  
    196190        return self._power_law(x) 
    197191        
    198     def inverse_transform_data(self, x, y): 
    199         """ 
    200             given a scale and a radius transform x, y using a inverse power_law 
    201             function 
    202         """ 
    203         result_x = numpy.array([math.log(i) for i in x])    
    204         result_y = numpy.array([math.log(j) for j in y]) 
    205          
    206         return result_x, result_y  
    207      
    208192    def _power_law(self, x): 
    209193        """ 
    … …  
    227211class Extrapolator: 
    228212    """ 
    229         compute f(x) 
     213        Extrapolate I(q) distribution using a given model 
    230214    """ 
    231215    def __init__(self, data): 
    … …  
    249233        self.set_fit_range() 
    250234       
    251     def set_fit_range(self ,qmin=None, qmax=None): 
     235    def set_fit_range(self, qmin=None, qmax=None): 
    252236        """ to set the fit range""" 
    253237        if qmin is not None: 
    … …  
    305289             
    306290            return a, b 
     291         
    307292 
    308293class InvariantCalculator(object): 
    … …  
    364349        return  new_data 
    365350      
    366     def _fit(self, function, qmin=Q_MINIMUM, qmax=Q_MAXIMUM, power=None): 
     351    def _fit(self, model, qmin=Q_MINIMUM, qmax=Q_MAXIMUM, power=None): 
    367352        """ 
    368353            fit data with function using  
    … …  
    379364                    for power_law will be the power value 
    380365        """ 
    381         transform = function 
    382         fit_x, fit_y = transform.inverse_transform_data(self._data.x, self._data.y) 
    383         fit_dy = transform.transform_error(y=self._data.y, dy=self._data.dy) 
    384         qmin = transform.transform_value(value=qmin) 
    385         qmax = transform.transform_value(value=qmax) 
    386          
     366        # Linearize the data in preparation for fitting 
     367        fit_x, fit_y, fit_dy = model.linearize_data(self._data.x, self._data.y, self._data.dy) 
     368       
     369        qmin = model.linearize_q_value(qmin) 
     370        qmax = model.linearize_q_value(qmax) 
     371         
     372        # Create a new Data1D object for processing 
    387373        fit_data = LoaderData1D(x=fit_x, y=fit_y, dy=fit_dy) 
    388374        fit_data.dxl = self._data.dxl 
    … …  
    392378        b, a = extrapolator.fit(power=power)  
    393379         
    394         return function.transform_values(a=a, b=b) 
     380        return model.extract_model_parameters(a=a, b=b) 
    395381     
    396382    def _get_qstar(self, data): 
    … …  
    633619        return self._get_qstar(data=data) 
    634620         
     621    def _get_extrapolated_data(self, model, npts=INTEGRATION_NSTEPS, 
     622                              q_start=Q_MINIMUM, q_end=Q_MAXIMUM): 
     623        """ 
     624            @return extrapolate data create from data 
     625        """ 
     626        #create new Data1D to compute the invariant 
     627        q = numpy.linspace(start=q_start, 
     628                               stop=q_end, 
     629                               num=npts, 
     630                               endpoint=True) 
     631        iq = model.evaluate_model(q) 
     632          
     633        # Determine whether we are extrapolating to high or low q-values 
     634        # If the data is slit smeared, get the index of the slit dimension array entry 
     635        # that we will use to smear the extrapolated data. 
     636        dxl = None 
     637        dxw = None 
     638 
     639        if self._data.is_slit_smeared(): 
     640            if q_start<min(self._data.x): 
     641                smear_index = 0 
     642            elif q_end>max(self._data.x): 
     643                smear_index = len(self._data.x)-1 
     644            else: 
     645                raise RuntimeError, "Extrapolation can only be evaluated for points outside the data Q range" 
     646             
     647            if self._data.dxl is not None : 
     648                dxl = numpy.ones(INTEGRATION_NSTEPS) 
     649                dxl = dxl * self._data.dxl[smear_index] 
     650                
     651            if self._data.dxw is not None : 
     652                dxw = numpy.ones(INTEGRATION_NSTEPS) 
     653                dxw = dxw * self._data.dxw[smear_index] 
     654              
     655        result_data = LoaderData1D(x=q, y=iq) 
     656        result_data.dxl = dxl 
     657        result_data.dxw = dxw 
     658        return result_data 
     659         
    635660    def get_extra_data_low(self): 
    636661        """ 
    … …  
    655680        """ 
    656681         
    657         # Data boundaries for fiiting 
     682        # Data boundaries for fitting 
    658683        qmin = self._data.x[0] 
    659684        qmax = self._data.x[self._low_extrapolation_npts - 1] 
     685         
    660686        # Extrapolate the low-Q data 
    661         #TODO: this fit fails. Fix it. 
    662         a, b = self._fit(function=self._low_extrapolation_function, 
    663                           qmin=qmin, 
    664                           qmax=qmax, 
    665                           power=self._low_extrapolation_power) 
     687        a, b = self._fit(model=self._low_extrapolation_function, 
     688                         qmin=qmin, 
     689                         qmax=qmax, 
     690                         power=self._low_extrapolation_power) 
    666691        #q_start point 
    667692        q_start = Q_MINIMUM 
    … …  
    669694            q_start = qmin/10 
    670695         
    671         data_min = self._low_extrapolation_function.get_extrapolated_data(data=self._data,  
    672                                             npts=INTEGRATION_NSTEPS, 
    673                               q_start=q_start, q_end=qmin, smear_indice=0) 
     696        data_min = self._get_extrapolated_data(model=self._low_extrapolation_function, 
     697                                               npts=INTEGRATION_NSTEPS, 
     698                                               q_start=q_start, q_end=qmin) 
    674699        return data_min 
    675700           
    … …  
    693718        qmax = self._data.x[x_len] 
    694719        q_end = Q_MAXIMUM 
    695         smear_indice = 0 
    696         if self._data.dxl is not None: 
    697             smear_indice = len(self._data.dxl) - 1 
    698720         
    699721        # fit the data with a model to get the appropriate parameters 
    700         a, b = self._fit(function=self._high_extrapolation_function, 
    701                                qmin=qmin, 
    702                                 qmax=qmax, 
    703                                 power=self._high_extrapolation_power) 
     722        a, b = self._fit(model=self._high_extrapolation_function, 
     723                         qmin=qmin, 
     724                         qmax=qmax, 
     725                         power=self._high_extrapolation_power) 
     726         
    704727        #create new Data1D to compute the invariant 
    705         data_max = self._high_extrapolation_function.get_extrapolated_data(data=self._data,  
    706                                             npts=INTEGRATION_NSTEPS, 
    707                                             q_start=qmax, q_end=q_end, 
    708                                             smear_indice=smear_indice) 
     728        data_max = self._get_extrapolated_data(model=self._high_extrapolation_function, 
     729                                               npts=INTEGRATION_NSTEPS, 
     730                                               q_start=qmax, q_end=q_end) 
    709731        return data_max 
    710732      

  • Invariant/test/PolySpheres.txt

    ref9ed58 raafa962  
     10 0 0 
    120.003   5.58459588740281        1.18158747955905 
    230.0031098987853131      5.58031083426704        0.590567038651574 

  • Invariant/test/utest_data_handling.py

    r6939bd4 raafa962  
    3131         
    3232        # Create invariant object. Background and scale left as defaults. 
    33         fit = invariant.FitFunctor(data=self.data) 
     33        fit = invariant.Extrapolator(data=self.data) 
    3434        a,b = fit.fit() 
    3535 
    … …  
    4848             
    4949        # Create invariant object. Background and scale left as defaults. 
    50         fit = invariant.FitFunctor(data=self.data) 
     50        fit = invariant.Extrapolator(data=self.data) 
    5151        a,b = fit.fit() 
    5252 
    … …  
    113113        """ 
    114114        self.scale = 1.5 
    115         self.rg = 70.0 
     115        self.rg = 30.0 
    116116        x = numpy.arange(0.0001, 0.1, 0.0001) 
    117117        y = numpy.asarray([self.scale * math.exp( -(q*self.rg)**2 / 3.0 ) for q in x]) 
    … …  
    129129         
    130130        self.assertEqual(inv._low_extrapolation_npts, 20) 
    131         self.assertEqual(inv._low_extrapolation_function.__name__, 'guinier') 
     131        self.assertEqual(inv._low_extrapolation_function.__class__, invariant.Guinier) 
    132132         
    133133        # Data boundaries for fiiting 
    … …  
    136136         
    137137        # Extrapolate the low-Q data 
    138         a, b = inv._fit(function=inv._low_extrapolation_function, 
     138        a, b = inv._fit(model=inv._low_extrapolation_function, 
    139139                          qmin=qmin, 
    140140                          qmax=qmax, 
    … …  
    172172         
    173173        self.assertEqual(inv._low_extrapolation_npts, 20) 
    174         self.assertEqual(inv._low_extrapolation_function.__name__, 'power_law') 
     174        self.assertEqual(inv._low_extrapolation_function.__class__, invariant.PowerLaw) 
    175175         
    176176        # Data boundaries for fitting 
    … …  
    179179         
    180180        # Extrapolate the low-Q data 
    181         a, b = inv._fit(function=inv._low_extrapolation_function, 
     181        a, b = inv._fit(model=inv._low_extrapolation_function, 
    182182                          qmin=qmin, 
    183183                          qmax=qmax, 
    … …  
    186186        self.assertAlmostEqual(self.scale, a, 6) 
    187187        self.assertAlmostEqual(self.m, b, 6) 
    188      
     188         
     189class TestLinearization(unittest.TestCase): 
     190     
     191    def test_guinier_incompatible_length(self): 
     192        g = invariant.Guinier() 
     193        self.assertRaises(AssertionError, g.linearize_data, [1],[1,2],None) 
     194        self.assertRaises(AssertionError, g.linearize_data, [1,2],[1,2],[1]) 
     195     
     196    def test_linearization(self): 
     197        """ 
     198            Check that the linearization process filters out points 
     199            that can't be transformed 
     200        """ 
     201        x = numpy.asarray(numpy.asarray([0,1,2,3])) 
     202        y = numpy.asarray(numpy.asarray([1,1,1,1])) 
     203        g = invariant.Guinier() 
     204        x_out, y_out, dy_out = g.linearize_data(x,y,None) 
     205        self.assertEqual(len(x_out), 3) 
     206        self.assertEqual(len(y_out), 3) 
     207        self.assertEqual(len(dy_out), 3) 
     208     

  • Invariant/test/utest_use_cases.py

    r46d50ca raafa962  
    33   
    44""" 
     5#TODO: there's no test for smeared extrapolation 
    56import unittest 
    67import numpy 
    … …  
    2425         
    2526        # Create invariant object. Background and scale left as defaults. 
    26         fit = invariant.FitFunctor(data=self.data) 
     27        fit = invariant.Extrapolator(data=self.data) 
    2728         
    2829        ##Without holding 
    … …  
    4041         
    4142        # Create invariant object. Background and scale left as defaults. 
    42         fit = invariant.FitFunctor(data=self.data) 
     43        fit = invariant.Extrapolator(data=self.data) 
    4344         
    4445        #With holding a = -power =4 
    … …  
    6263         
    6364        # Create invariant object. Background and scale left as defaults. 
    64         fit = invariant.FitFunctor(data=self.data) 
     65        fit = invariant.Extrapolator(data=self.data) 
    6566         
    6667        ##Without holding 
    … …  
    7879         
    7980        # Create invariant object. Background and scale left as defaults. 
    80         fit = invariant.FitFunctor(data=self.data) 
     81        fit = invariant.Extrapolator(data=self.data) 
    8182         
    8283        #With holding a = -power =4 
Note: See TracChangeset for help on using the changeset viewer.