Changeset 9b6497bb in sasview for Invariant


Ignore:
Timestamp:
Dec 17, 2009 10:19:56 AM (14 years ago)
Author:
Jae Cho <jhjcho@…>
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:
472b11c
Parents:
8d1e745
Message:

math corrections

Location:
Invariant
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • Invariant/invariant.py

    r4e80ae0 r9b6497bb  
    66import math  
    77import numpy 
     8import scipy 
    89 
    910from DataLoader.data_info import Data1D as LoaderData1D 
    … …  
    2021INTEGRATION_NSTEPS = 1000 
    2122 
    22 def guinier(x, scale=1, radius=0.1): 
     23def guinier(x, scale=1, radius=60): 
    2324    """ 
    2425        Compute a F(x) = scale* e-((radius*x)**2/3). 
    … …  
    2829        @return F(x) 
    2930    """    
     31    if radius <= 0: 
     32        raise ValueError("Rg expected positive value, but got %s"%radius)   
    3033    value = numpy.array([math.exp(-((radius * i)**2/3)) for i in x ])  
     34    scale = numpy.sqrt(scale*scale)       
     35    if scale == 0: 
     36        raise ValueError("scale expected positive value, but got %s"%scale)  
    3137    return scale * value 
    3238 
    33 def power_law(x, scale=1, power=0): 
     39def power_law(x, scale=1, power=4): 
    3440    """ 
    3541        F(x) = scale* (x)^(-power) 
    … …  
    4349    """   
    4450    if power <=0: 
    45         raise ValueError("Power_law function expected posive power, but got %s"%power) 
     51        raise ValueError("Power_law function expected positive power, but got %s"%power) 
    4652     
    4753    value = numpy.array([ math.pow(i, -power) for i in x ])   
     54    scale = numpy.sqrt(scale*scale) 
     55    if scale == 0: 
     56        raise ValueError("scale expected positive value, but got %s"%scale)  
    4857    return scale * value 
    4958 
    … …  
    111120            
    112121        """ 
     122         
     123        fx = numpy.zeros(len(self.data.x)) 
     124        sigma = numpy.zeros(len(self.data.x)) 
     125 
     126        #define dy^2 
     127        sigma = self.data.dy[self.idx_unsmeared ] 
     128        sigma2 = sigma*sigma 
     129 
    113130        # Compute theory data f(x) 
    114         fx = numpy.zeros(len(self.data.x)) 
    115         fx = self.data.y[self.idx_unsmeared ] 
     131        fx = self.data.y[self.idx_unsmeared ]/sigma2 
    116132        ## Smear theory data 
    117133        if self.smearer is not None: 
    118134            fx = self.smearer(fx, self._first_unsmeared_bin,self._last_unsmeared_bin) 
    119135     
    120         A = numpy.vstack([ self.data.x[self.idx], 
    121                            numpy.ones(len(self.data.x[self.idx]))]).T 
     136        A = numpy.vstack([ self.data.x[self.idx]/sigma2, 
     137                           numpy.ones(len(self.data.x[self.idx]))/sigma2]).T 
    122138 
    123139        a, b = numpy.linalg.lstsq(A, fx)[0] 
    … …  
    178194            raise ValueError,"Data must be of type DataLoader.Data1D" 
    179195        new_data = self._scale * data - self._background 
     196        new_data.dy = data.dy 
    180197        new_data.dxl = data.dxl 
    181198        new_data.dxw = data.dxw 
    … …  
    199216        if function.__name__ == "guinier": 
    200217            fit_x = numpy.array([x * x for x in self._data.x]) 
     218             
    201219            qmin = qmin**2 
    202220            qmax = qmax**2 
    203221            fit_y = numpy.array([math.log(y) for y in self._data.y]) 
     222            fit_dy = numpy.array([y for y in self._data.y]) 
     223            fit_dy = numpy.array([dy for dy in self._data.dy])/fit_dy 
     224 
    204225        elif function.__name__ == "power_law": 
    205226            if power is None: 
    206227                fit_x = numpy.array([math.log(x) for x in self._data.x]) 
     228 
    207229                qmin = math.log(qmin) 
    208230                qmax = math.log(qmax) 
    209                 fit_y = numpy.array([math.log(y) for y in self._data.y]) 
     231            fit_y = numpy.array([math.log(y) for y in self._data.y]) 
     232            fit_dy = numpy.array([y for y in self._data.y]) 
     233            fit_dy = numpy.array([dy for dy in self._data.dy])/fit_dy 
     234 
    210235        else: 
    211236            raise ValueError("Unknown function used to fit %s"%function.__name__) 
    212              
    213          
    214         if function.__name__ == "power_law" and  power is None: 
     237         
     238        if function.__name__ == "power_law" and  power != None: 
    215239            b = math.fabs(power) 
    216             fit_x = numpy.array([math.pow(x, -b) for x in self._data.x]) 
    217240            fit_y = numpy.array([math.log(y) for y in self._data.y]) 
    218             a = (fit_y - fit_x)/(len(self._data.x)) 
    219         else: 
    220             fit_data = LoaderData1D(x=fit_x, y=fit_y) 
     241            fit_dy = numpy.array([y for y in self._data.y]) 
     242            fit_dy = numpy.array([dy for dy in self._data.dy])/fit_dy 
     243            sigma2 = fit_dy*fit_dy 
     244            a = scipy.sum(fit_y/sigma2) - scipy.sum(fit_x/sigma2*b)/scipy.sum(sigma2) 
     245        else: 
     246            fit_data = LoaderData1D(x=fit_x, y=fit_y, dy=fit_dy) 
    221247            fit_data.dxl = self._data.dxl 
    222248            fit_data.dxw = self._data.dxw    
    … …  
    225251            b, a = functor.fit() 
    226252         
     253                   
    227254        if function.__name__ == "guinier": 
    228255            # b is the radius value of the guinier function 
    229             b = math.sqrt(-3 * b) 
     256            if b>=0: 
     257                raise ValueError("Guinier fit was not converged") 
     258            else: 
     259                b = math.sqrt(-3 * b) 
     260 
     261 
    230262        if function.__name__ == "power_law": 
    231             b = -1 * b 
     263            if power == None: 
     264                b = -1 * b 
    232265            if b <= 0: 
    233266                raise ValueError("Power_law fit expected posive power, but got %s"%power) 
    234267        # a is the scale of the guinier function 
    235268        a = math.exp(a) 
     269         
    236270        return a, b 
    237271     
    … …  
    502536            @return: a new data of type Data1D 
    503537        """ 
     538         
    504539        # Data boundaries for fiiting 
    505540        qmin = self._data.x[0] 
    … …  
    514549        except: 
    515550            return None 
    516         
     551         
    517552        #q_start point 
    518553        q_start = Q_MINIMUM 
    … …  
    534569            dxw = numpy.ones(INTEGRATION_NSTEPS) 
    535570            dxw = dxw * self._data.dxw[0] 
    536              
     571 
    537572        data_min = LoaderData1D(x=new_x, y=new_y) 
    538573        data_min.dxl = dxl 
    539574        data_min.dxw = dxw 
    540575        self._data.clone_without_data( clone= data_min) 
    541          
     576 
    542577        return data_min 
    543578           

  • Invariant/test/utest_use_cases.py

    r4e80ae0 r9b6497bb  
    345345         
    346346        # Test results 
    347         self.assertAlmostEquals(qstar, 0.002037677) 
     347        self.assertAlmostEquals(qstar, 0.00138756,2) 
    348348        self.assertAlmostEquals(v, 0.115352622) 
    349349        self.assertAlmostEquals(s + _ERR_SURFACE, 9.42e+2, 1) 
    … …  
    389389         
    390390        # Test results 
    391         self.assertAlmostEquals(qstar, 0.0021621, 3) 
     391        self.assertAlmostEquals(qstar, 0.00460319, 3) 
    392392        self.assertAlmostEquals(v, 0.202846825) 
    393393        self.assertAlmostEquals(s+ _ERR_SURFACE, 9.42e+2, 1) 
Note: See TracChangeset for help on using the changeset viewer.