Changeset 437a9f0 in sasview for Invariant


Ignore:
Timestamp:
Dec 1, 2009 4:21:19 PM (14 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:
7a108dd
Parents:
c5607fa
Message:

testing invariant

Location:
Invariant
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • Invariant/invariant.py

    ref9ed58 r437a9f0  
    88 
    99from DataLoader.data_info import Data1D as LoaderData1D 
    10 #from DataLoader.data_info import Data1D as LoaderData1D 
    1110from DataLoader.qsmearing import smear_selection 
    1211 
    … …  
    5049        compute f(x) 
    5150    """ 
    52     def __init__(self,data , function): 
    53         """ 
    54             @param function :the function used for computing residuals 
    55             @param Data: data used for computing residuals 
    56         """ 
    57         self.function = function 
     51    def __init__(self,data): 
     52        """ 
     53            Determine a and b given a linear equation y = ax + b 
     54            @param Data: data containing x and y  such as  y = ax + b  
     55        """ 
    5856        self.data  = data 
    5957        x_len = len(self.data.x) -1 
    … …  
    111109            
    112110        """ 
    113         fx = self.data.y 
     111        # Compute theory data f(x) 
     112        fx = numpy.zeros(len(self.data.x)) 
     113        fx = self.data.y[self.idx_unsmeared ] 
    114114        ## Smear theory data 
    115115        if self.smearer is not None: 
    116             fx = self.smearer(fx, self._first_unsmeared_bin, 
    117                                self._last_unsmeared_bin)                                                          
    118         A = numpy.vstack([ self.data.x, numpy.ones(len(self.data.x))]).T 
     116            fx = self.smearer(fx, self._first_unsmeared_bin,self._last_unsmeared_bin) 
     117     
     118        A = numpy.vstack([ self.data.x[self.idx], 
     119                           numpy.ones(len(self.data.x[self.idx]))]).T 
    119120 
    120121        a, b = numpy.linalg.lstsq(A, fx)[0] 
    … …  
    223224        """ 
    224225        if function.__name__ == "guinier": 
    225             fit_x = self._data.x * self._data.x 
     226            fit_x = numpy.array([x * x for x in self._data.x]) 
    226227            qmin = qmin**2 
    227228            qmax = qmax**2 
    228             fit_y = [math.log(y) for y in self._data.y] 
     229            fit_y = numpy.array([math.log(y) for y in self._data.y]) 
    229230        elif function.__name__ == "power_law": 
    230             fit_x = [math.log(x) for x in self._data.x] 
     231            fit_x = numpy.array([math.log(x) for x in self._data.x]) 
    231232            qmin = math.log(qmin) 
    232233            qmax = math.log(qmax) 
    233             fit_y = [math.log(y) for y in self._data.y] 
     234            fit_y = numpy.array([math.log(y) for y in self._data.y]) 
    234235        else: 
    235236            raise ValueError("Unknown function used to fit %s"%function.__name__) 
    … …  
    239240        fit_data.dxw = self._data.dxw    
    240241         
    241         functor = FitFunctor(data=fit_data, function= function) 
     242        functor = FitFunctor(data=fit_data) 
    242243        functor.set_fit_range(qmin=qmin, qmax=qmax) 
    243         a, b = functor.fit() 
     244        b, a = functor.fit() 
    244245         
    245246        if function.__name__ == "guinier": 
    246247            # b is the radius value of the guinier function 
    247             print "b",b 
    248248            b = math.sqrt(-3 * b) 
    249         
     249        if function.__name__ == "power_law": 
     250            b = -1 * b 
    250251        # a is the scale of the guinier function 
    251252        a = math.exp(a) 
    252         return a, math.fabs(b) 
     253        return a, b 
    253254     
    254255    def _get_qstar(self, data): 
    … …  
    446447            n = len(data.x) - 1 
    447448            #compute the first delta 
    448             dx0 = (data.x[1] - data.x[0])/2 
     449            dx0 = (data.x[1] + data.x[0])/2 
    449450            #compute the last delta 
    450451            dxn= data.x[n] - data.x[n-1] 
    … …  
    499500            n = len(data.x) - 1 
    500501            #compute the first delta 
    501             dx0 = (data.x[1] - data.x[0])/2 
     502            dx0 = (data.x[1] + data.x[0])/2 
    502503            #compute the last delta 
    503504            dxn = data.x[n] - data.x[n-1] 
    … …  
    570571            raise RuntimeError, "invariant must be greater than zero" 
    571572         
    572         print "self._qstar",self._qstar 
    573573        # Compute intermediate constant 
    574574        k =  1.e-8 * self._qstar/(2 * (math.pi * math.fabs(float(contrast)))**2) 
    575575        #Check discriminant value 
    576576        discrim = 1 - 4 * k 
    577         print "discrim",k,discrim 
     577         
    578578        # Compute volume fraction 
    579579        if discrim < 0: 
    580580            raise RuntimeError, "could not compute the volume fraction: negative discriminant" 
    581         elif discrim ==0: 
     581        elif discrim == 0: 
    582582            return 1/2 
    583583        else: 
    584             volume1 = 0.5 *(1 - math.sqrt(discrim)) 
    585             volume2 = 0.5 *(1 + math.sqrt(discrim)) 
     584            volume1 = 0.5 * (1 - math.sqrt(discrim)) 
     585            volume2 = 0.5 * (1 + math.sqrt(discrim)) 
    586586             
    587587            if 0 <= volume1 and volume1 <= 1: 
    … …  
    640640        # Data boundaries for fiiting 
    641641        qmin = self._data.x[0] 
    642         qmax = self._data.x[self._low_extrapolation_npts] 
     642        qmax = self._data.x[self._low_extrapolation_npts - 1] 
    643643         
    644644        try: 
    … …  
    648648        except: 
    649649            raise 
    650             return None 
     650            #return None 
    651651        
    652652        #create new Data1D to compute the invariant 
    … …  
    655655                               num=INTEGRATION_NSTEPS, 
    656656                               endpoint=True) 
    657         new_y = self._low_extrapolation_function(x=new_x, 
    658                                                           scale=a,  
    659                                                           radius=b) 
     657        new_y = self._low_extrapolation_function(x=new_x, scale=a, radius=b) 
    660658        dxl = None 
    661659        dxw = None 
    … …  
    664662            dxl = dxl * self._data.dxl[0] 
    665663        if self._data.dxw is not None: 
    666             dwl = numpy.ones(1, INTEGRATION_NSTEPS) 
    667             dwl = dwl * self._data.dwl[0] 
     664            dxw = numpy.ones(1, INTEGRATION_NSTEPS) 
     665            dxw = dxw * self._data.dxw[0] 
    668666             
    669667        data_min = LoaderData1D(x=new_x, y=new_y) 
    … …  
    690688        # Data boundaries for fiiting 
    691689        x_len = len(self._data.x) - 1 
    692         qmin = self._data.x[x_len - self._high_extrapolation_npts] 
     690        qmin = self._data.x[x_len - (self._high_extrapolation_npts - 1) ] 
    693691        qmax = self._data.x[x_len] 
    694692         
    … …  
    699697        except: 
    700698            raise 
    701             return None 
     699            #return None 
    702700        
    703701        #create new Data1D to compute the invariant 
    … …  
    706704                               num=INTEGRATION_NSTEPS, 
    707705                               endpoint=True) 
    708         new_y = self._high_extrapolation_function(x=new_x, 
    709                                                   scale=a,  
    710                                                   power=b) 
     706        
     707        new_y = self._high_extrapolation_function(x=new_x, scale=a, power=b) 
     708         
    711709        dxl = None 
    712710        dxw = None 
    … …  
    715713            dxl = dxl * self._data.dxl[0] 
    716714        if self._data.dxw is not None: 
    717             dwl = numpy.ones(1, INTEGRATION_NSTEPS) 
    718             dwl = dwl * self._data.dwl[0] 
     715            dxw = numpy.ones(1, INTEGRATION_NSTEPS) 
     716            dxw = dxw * self._data.dxw[0] 
    719717             
    720718        data_max = LoaderData1D(x=new_x, y=new_y) 
    … …  
    746744            dV = 0.5 * (4*k* dq_star) /(2* math.sqrt(1-k* q_star)) 
    747745                                  
    748             for k = 10^(8)*q_star/(2*(pi*|contrast|)**2) 
     746            for k = 10^(-8)*q_star/(2*(pi*|contrast|)**2) 
    749747             
    750748            q_star: the invariant value including extrapolated value if existing 
    … …  
    760758            raise ValueError, "invariant must be greater than zero" 
    761759         
    762         k =  1.e-8 * self._qstar /(2*(math.pi* math.fabs(float(contrast)))**2) 
     760        k =  1.e-8 * self._qstar /(2 * (math.pi* math.fabs(float(contrast)))**2) 
    763761        #check value inside the sqrt function 
    764762        value = 1 - k * self._qstar 
    … …  
    766764            raise ValueError, "Cannot compute incertainty on volume" 
    767765        # Compute uncertainty 
    768         uncertainty = (0.5 * 4 * k * self._qstar_err)/(2 * math.sqrt(1- k * self._qstar)) 
     766        uncertainty = (0.5 * 4 * k * self._qstar_err)/(2 * math.sqrt(1 - k * self._qstar)) 
    769767         
    770768        return volume, math.fabs(uncertainty) 

  • Invariant/test/utest_use_cases.py

    ref9ed58 r437a9f0  
    9393        #    parameters. For instance, you might not want to allow 'high' and 'guinier'. 
    9494        # The power parameter (not shown below) should default to 4. 
    95         inv.set_extrapolation(range='low', npts=20, function='guinier') 
     95        inv.set_extrapolation(range='low', npts=10, function='guinier') 
    9696         
    9797        # The version of the call without error 
    … …  
    108108         
    109109        # Test results 
    110         self.assertAlmostEquals(qstar, 7.49e-5,2) 
     110        self.assertAlmostEquals(qstar, 7.49e-5) 
    111111        self.assertAlmostEquals(v, 0.005648401) 
    112112        self.assertAlmostEquals(s, 9.42e+2,2) 
    … …  
    120120         
    121121        # Set the extrapolation parameters for the high-Q range 
    122         inv.set_extrapolation(range='high', npts=20, function='power_law', power=4) 
     122        inv.set_extrapolation(range='high', npts=10, function='power_law', power=4) 
    123123         
    124124        # The version of the call without error 
    … …  
    146146         
    147147        # Set the extrapolation parameters for the low- and high-Q ranges 
    148         inv.set_extrapolation(range='low', npts=5, function='guinier') 
    149         inv.set_extrapolation(range='high', npts=5, function='power_law', power=4) 
     148        inv.set_extrapolation(range='low', npts=10, function='guinier') 
     149        inv.set_extrapolation(range='high', npts=10, function='power_law', power=4) 
    150150         
    151151        # The version of the call without error 
    … …  
    161161         
    162162        # Test results 
    163         self.assertAlmostEquals(qstar, 7.90069e-5,2) 
     163        self.assertAlmostEquals(qstar, 7.88981e-5,2) 
    164164        self.assertAlmostEquals(v, 0.005952674) 
    165165        self.assertAlmostEquals(s, 9.42e+2,2) 
    … …  
    213213        inv = invariant.InvariantCalculator(data=self.data_slit_smear) 
    214214        # Set the extrapolation parameters for the low-Q range 
    215         inv.set_extrapolation(range='low', npts=20, function='guinier') 
     215        inv.set_extrapolation(range='low', npts=10, function='guinier') 
    216216        # The version of the call without error 
    217217        # At this point, we could still compute Q* without extrapolation by calling 
    … …  
    237237         
    238238        # Set the extrapolation parameters for the high-Q range 
    239         inv.set_extrapolation(range='high', npts=20, function='power_law', power=4) 
     239        inv.set_extrapolation(range='high', npts=10, function='power_law', power=4) 
    240240         
    241241        # The version of the call without error 
    … …  
    263263         
    264264        # Set the extrapolation parameters for the low- and high-Q ranges 
    265         inv.set_extrapolation(range='low', npts=20, function='guinier') 
    266         inv.set_extrapolation(range='high', npts=20, function='power_law', power=4) 
     265        inv.set_extrapolation(range='low', npts=10, function='guinier') 
     266        inv.set_extrapolation(range='high', npts=10, function='power_law', power=4) 
    267267         
    268268        # The version of the call without error 
    … …  
    351351        inv = invariant.InvariantCalculator(data=self.data_q_smear) 
    352352        # Set the extrapolation parameters for the high-Q range 
    353         inv.set_extrapolation(range='high', npts=20, function='power_law', power=4) 
     353        inv.set_extrapolation(range='high', npts=10, function='power_law', power=4) 
    354354        # The version of the call without error 
    355355        qstar = inv.get_qstar(extrapolation='high') 
    … …  
    372372        inv = invariant.InvariantCalculator(data=self.data_q_smear) 
    373373        # Set the extrapolation parameters for the low- and high-Q ranges 
    374         inv.set_extrapolation(range='low', npts=20, function='guinier') 
    375         inv.set_extrapolation(range='high', npts=20, function='power_law', power=4) 
     374        inv.set_extrapolation(range='low', npts=10, function='guinier') 
     375        inv.set_extrapolation(range='high', npts=10, function='power_law', power=4) 
    376376        # The version of the call without error 
    377377        # The function parameter defaults to None, then is picked to be 'power_law' for extrapolation='high' 
    … …  
    384384         
    385385        # Test results 
    386         self.assertAlmostEquals(qstar, 0.002157677) 
     386        self.assertAlmostEquals(qstar, 0.0021621) 
    387387        self.assertAlmostEquals(v, 0.202846825) 
    388388        self.assertAlmostEquals(s, 9.42e+2,2) 
Note: See TracChangeset for help on using the changeset viewer.