Changeset d4bf55e9 in sasview for sansmodels


Ignore:
Timestamp:
Mar 10, 2011 11:41:51 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:
ae83ad3
Parents:
510e7ad
Message:

Now 2d smearing uses polar symmetry (not rectangular)

Location:
sansmodels/src/sans/models
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • sansmodels/src/sans/models/LineModel.py

    r988130c6 rd4bf55e9  
    66from sans.models.BaseComponent import BaseComponent 
    77import math 
     8import numpy 
    89 
    910 
    … …  
    7071        else: 
    7172            return self._line(x) 
     73         
     74    def evalDistribution(self, qdist): 
     75        """ 
     76        Evaluate a distribution of q-values. 
     77         
     78        * For 1D, a numpy array is expected as input: 
     79         
     80            evalDistribution(q) 
     81             
     82        where q is a numpy array. 
     83         
     84         
     85        * For 2D, a list of numpy arrays are expected: [qx_prime,qy_prime], 
     86          where 1D arrays, 
     87                 
     88        :param qdist: ndarray of scalar q-values or list [qx,qy]  
     89                    where qx,qy are 1D ndarrays  
     90         
     91        """ 
     92        if qdist.__class__.__name__ == 'list': 
     93            # Check whether we have a list of ndarrays [qx,qy] 
     94            if len(qdist)!=2 or \ 
     95                qdist[0].__class__.__name__ != 'ndarray' or \ 
     96                qdist[1].__class__.__name__ != 'ndarray': 
     97                    raise RuntimeError, "evalDistribution expects a list of 2 ndarrays" 
     98                 
     99            # Extract qx and qy for code clarity 
     100            qx = qdist[0] 
     101            qy = qdist[1] 
     102            #For 2D, Z = A + B * Y,  
     103            # so that it keeps its linearity in y-direction. 
     104            # calculate q_r component for 2D isotropic 
     105            q =  qy 
     106            # vectorize the model function runXY 
     107            v_model = numpy.vectorize(self.runXY,otypes=[float]) 
     108            # calculate the scattering 
     109            iq_array = v_model(q) 
     110 
     111            return iq_array 
     112                 
     113        elif qdist.__class__.__name__ == 'ndarray': 
     114                # We have a simple 1D distribution of q-values 
     115                v_model = numpy.vectorize(self.runXY,otypes=[float]) 
     116                iq_array = v_model(qdist) 
     117 
     118                return iq_array 
     119             
     120        else: 
     121            mesg = "evalDistribution is expecting an ndarray of scalar q-values" 
     122            mesg += " or a list [qx,qy] where qx,qy are 2D ndarrays." 
     123            raise RuntimeError, mesg 
     124         
     125     
    72126    
    73127 

  • sansmodels/src/sans/models/smearing_2d.py

    r04eb1a4 rd4bf55e9  
    1414## Limit of how many sigmas to be covered for the Gaussian smearing 
    1515# default: 2.5 to cover 98.7% of Gaussian 
    16 LIMIT = 2.5 
     16LIMIT = 3.0 
    1717## Defaults 
    1818R_BIN = {'Xhigh':10.0, 'High':5.0,'Med':5.0,'Low':3.0} 
    … …  
    2525      
    2626    def __init__(self, data=None, model=None, index=None,  
    27                  limit=LIMIT, accuracy='Low'): 
     27                 limit=LIMIT, accuracy='Low', coords='polar'): 
    2828        """ 
    2929        Assumption: equally spaced bins in dq_r, dq_phi space. 
    … …  
    3434        :param nr: number of bins in dq_r-axis 
    3535        :param nphi: number of bins in dq_phi-axis 
     36        :param coord: coordinates [string], 'polar' or 'cartesian' 
    3637        """ 
    3738        ## data 
    … …  
    4849        self.limit = limit 
    4950        self.index = index 
     51        self.coords = coords 
    5052        self.smearer = True 
    5153         
    … …  
    6163        self.qx_data = self.data.qx_data[self.index] 
    6264        self.qy_data = self.data.qy_data[self.index] 
     65        self.q_data = self.data.q_data[self.index] 
     66        # Here dqx and dqy mean dq_parr and dq_perp 
    6367        self.dqx_data = self.data.dqx_data[self.index] 
    6468        self.dqy_data = self.data.dqy_data[self.index] 
    … …  
    113117    def get_value(self): 
    114118        """ 
    115         Over sampling of r_nbins times phi_nbins, calculate Gaussian weights, then find smeared intensity 
     119        Over sampling of r_nbins times phi_nbins, calculate Gaussian weights,  
     120        then find smeared intensity 
    116121        For the default values, this is equivalent (but by using numpy array  
    117122        the speed optimized by a factor of ten)to the following: :: 
    118123        
    119124        
    120             Remove the singular points if exists 
    121             self.dqx_data[self.dqx_data==0]=SIGMA_ZERO 
    122             self.dqy_data[self.dqy_data==0]=SIGMA_ZERO 
    123              
    124             for phi in range(0,4): 
    125                 for r in range(0,5): 
    126                     n = (phi)*5+(r) 
    127                     r = r+0.25 
    128                     dphi = phi*2.0*math.pi/4.0 + numpy.arctan(self.qy_data[index_model]/self.dqy_data[index_model]/self.qx_data[index_model]*/self.dqx_data[index_model]) 
    129                     dq = r*sqrt( self.dqx_data[index_model]*self.dqx_data[index_model] \ 
    130                         + self.dqy_data[index_model]*self.dqy_data[index_model] ) 
    131                     #integrant of exp(-0.5*r*r) r dr at each bins : The integration may not need. 
    132                     weight_res[n] = e^{(-0.5*((r-0.25)*(r-0.25)))}- e^{(-0.5*((r-0.25)*(r-0.25)))} 
    133                     #if phi != 0 and r != 0: 
    134                     qx_res = numpy.append(qx_res,self.qx_data[index_model]+ dq * cos(dphi)) 
    135                     qy_res = numpy.append(qy_res,self.qy_data[index_model]+ dq * sin(dphi)) 
     125        Remove the singular points if exists 
     126        self.dqx_data[self.dqx_data==0]=SIGMA_ZERO 
     127        self.dqy_data[self.dqy_data==0]=SIGMA_ZERO 
     128         
     129        for phi in range(0,4): 
     130            for r in range(0,5): 
     131                n = (phi)*5+(r) 
     132                r = r+0.25 
     133                dphi = phi*2.0*math.pi/4.0 + numpy.arctan( \ 
     134                        self.qy_data[index_model]/self.dqy_data[index_model]/ \ 
     135                        self.qx_data[index_model]*/self.dqx_data[index_model]) 
     136                dq = r*sqrt( self.dqx_data[index_model]*\ 
     137                        self.dqx_data[index_model] \ 
     138                    + self.dqy_data[index_model]*self.dqy_data[index_model] ) 
     139                #integrant of exp(-0.5*r*r) r dr at each bins :  
     140                The integration may not need. 
     141                weight_res[n] = e^{(-0.5*((r-0.25)*(r-0.25)))}- \ 
     142                                e^{(-0.5*((r-0.25)*(r-0.25)))} 
     143                #if phi != 0 and r != 0: 
     144                qx_res = numpy.append(qx_res,self.qx_data[index_model]+ \ 
     145                                                            dq * cos(dphi)) 
     146                qy_res = numpy.append(qy_res,self.qy_data[index_model]+ \ 
     147                                                            dq * sin(dphi)) 
    136148                     
    137149        Then compute I(qx_res,qy_res) and do weighted averaging.  
    … …  
    145157                        numpy.all(numpy.fabs(self.dqy_data <= 1.1e-10)): 
    146158            self.smearer = False 
     159   
    147160        if self.smearer == False: 
    148             return self.model.evalDistribution([self.qx_data,self.qy_data])  
     161            return self.model.evalDistribution([self.qx_data, self.qy_data])  
    149162 
    150163        nr = self.nr[self.accuracy] 
    … …  
    156169 
    157170        # Number of bins in the dqr direction (polar coordinate of dqx and dqy) 
    158         bin_size = self.limit/nr 
    159         # Total number of bins = # of bins in dq_r-direction times # of bins in dq_phi-direction 
     171        bin_size = self.limit / nr 
     172        # Total number of bins = # of bins  
     173        # in dq_r-direction times # of bins in dq_phi-direction 
    160174        n_bins = nr * nphi 
    161175        # Mean values of dqr at each bins ,starting from the half of bin size 
    162         r = bin_size/2.0+numpy.arange(nr)*bin_size 
     176        r = bin_size / 2.0 + numpy.arange(nr) * bin_size 
    163177        # mean values of qphi at each bines 
    164178        phi = numpy.arange(nphi) 
    165         dphi = phi*2.0*math.pi/nphi 
     179        dphi = phi * 2.0 * math.pi / nphi 
    166180        dphi = dphi.repeat(nr) 
    167         ## Transform to polar coordinate, and set dphi at each data points ; 1d array 
    168         dphi = dphi.repeat(len_data)+numpy.arctan(self.qy_data*self.dqx_data/self.qx_data/self.dqy_data).repeat(n_bins).reshape(len_data,n_bins).transpose().flatten() 
    169         ## Find Gaussian weight for each dq bins: The weight depends only on r-direction (The integration may not need) 
    170         weight_res = numpy.exp(-0.5*((r-bin_size/2.0)*(r-bin_size/2.0)))-numpy.exp(-0.5*((r+bin_size/2.0)*(r+bin_size/2.0))) 
    171         weight_res = weight_res.repeat(nphi).reshape(nr,nphi).transpose().flatten() 
     181 
     182        ## Transform to polar coordinate,  
     183        #  and set dphi at each data points ; 1d array 
     184        dphi = dphi.repeat(len_data) 
     185        q_phi = self.qy_data / self.qx_data 
     186         
     187        # Starting angle is different between polar and cartesian coordinates. 
     188        #if self.coords != 'polar': 
     189        #    dphi += numpy.arctan( q_phi * self.dqx_data/ \ 
     190        #                     self.dqy_data).repeat(n_bins).reshape(len_data,\ 
     191        #                                        n_bins).transpose().flatten() 
     192 
     193        # The angle (phi) of the original q point 
     194        q_phi = numpy.arctan(q_phi).repeat(n_bins).reshape(len_data,\ 
     195                                                n_bins).transpose().flatten() 
     196        ## Find Gaussian weight for each dq bins: The weight depends only  
     197        #  on r-direction (The integration may not need) 
     198        weight_res = numpy.exp(-0.5 * ((r - bin_size / 2.0) * \ 
     199                                    (r - bin_size / 2.0)))- \ 
     200                                    numpy.exp(-0.5 * ((r + bin_size / 2.0 ) *\ 
     201                                    (r + bin_size / 2.0))) 
     202        weight_res /= numpy.sum(weight_res) 
     203        weight_res = weight_res.repeat(nphi).reshape(nr, nphi) 
     204 
     205        weight_res = weight_res.transpose().flatten() 
    172206         
    173207        ## Set dr for all dq bins for averaging 
    … …  
    176210        dqx = numpy.outer(dr,self.dqx_data).flatten() 
    177211        dqy = numpy.outer(dr,self.dqy_data).flatten() 
    178         qx = self.qx_data.repeat(n_bins).reshape(len_data,n_bins).transpose().flatten() 
    179         qy = self.qy_data.repeat(n_bins).reshape(len_data,n_bins).transpose().flatten() 
    180  
    181          
    182         ## Over-sampled qx_data and qy_data. 
    183         qx_res = qx+ dqx*numpy.cos(dphi) 
    184         qy_res = qy+ dqy*numpy.sin(dphi) 
    185          
     212 
     213        qx = self.qx_data.repeat(n_bins).reshape(len_data,\ 
     214                                             n_bins).transpose().flatten() 
     215        qy = self.qy_data.repeat(n_bins).reshape(len_data,\ 
     216                                             n_bins).transpose().flatten() 
     217 
     218        # The polar needs rotation by -q_phi 
     219        if self.coords == 'polar': 
     220            qx_res = qx + (dqx*numpy.cos(dphi) * numpy.cos(-q_phi) +\ 
     221                           dqy*numpy.sin(dphi) * numpy.sin(-q_phi)) 
     222            qy_res = qy + (-dqx*numpy.cos(dphi) * numpy.sin(-q_phi) +\ 
     223                           dqy*numpy.sin(dphi) * numpy.cos(-q_phi)) 
     224        else: 
     225            qx_res = qx +  dqx*numpy.cos(dphi) 
     226            qy_res = qx +  dqy*numpy.sin(dphi) 
     227 
    186228        ## Evaluate all points 
    187         val = self.model.evalDistribution([qx_res,qy_res])  
     229        val = self.model.evalDistribution([qx_res, qy_res])  
    188230 
    189231        ## Reshape into 2d array to use numpy weighted averaging 
    … …  
    191233 
    192234        ## Averaging with Gaussian weighting: normalization included. 
    193         value =numpy.average(value_res,axis=0,weights=weight_res) 
     235        value =numpy.average(value_res,axis=0, weights=weight_res) 
    194236 
    195237        ## Return the smeared values in the range of self.index 
    … …  
    199241    ## Test w/ 2D linear function 
    200242    x = 0.001*numpy.arange(1,11) 
    201     dx = numpy.ones(len(x))*0.001 
     243    dx = numpy.ones(len(x))*0.0003 
    202244    y = 0.001*numpy.arange(1,11) 
    203245    dy = numpy.ones(len(x))*0.001 
    … …  
    214256    out.dqx_data = dx 
    215257    out.dqy_data = dy 
     258    out.q_data = numpy.sqrt(dx * dx + dy * dy) 
    216259    index = numpy.ones(len(x), dtype = bool) 
    217260    out.mask = index 
    … …  
    225268    ## All data are ones, so the smeared should also be ones. 
    226269    print "Data length =",len(value) 
    227     print " 2D linear function, I = 0 + 1*qx*qy" 
    228     print " Gaussian weighted averaging on a 2D linear function will provides the results same as without the averaging." 
     270    print " 2D linear function, I = 0 + 1*qy" 
     271    text = " Gaussian weighted averaging on a 2D linear function will " 
     272    text += "provides the results same as without the averaging." 
     273    print text 
    229274    print "qx_data", "qy_data", "I_nonsmear", "I_smeared" 
    230275    for ind in range(len(value)): 
    231276        print x[ind],y[ind],model.evalDistribution([x,y])[ind], value[ind] 
    232    
     277         
     278""" 
     279for i in range(len(qx_res)/(128*128)): 
     280    k = i * 128*128 +64 
     281 
     282    print qx_res[k]-qqx[k], qy_res[k]-qqy[k] 
     283print qqx[64],qqy[64] 
     284"""  
    233285"""     
    234286if __name__ == '__main__': 
Note: See TracChangeset for help on using the changeset viewer.