source: sasview/sansmodels/src/sans/models/smearing_2d.py @ 8bac371

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since 8bac371 was 87615a48, checked in by Jae Cho <jhjcho@…>, 13 years ago

added c-extension for smearer2d_helper(but locked it for now)

  • Property mode set to 100644
File size: 11.3 KB
Line 
1#####################################################################
2#This software was developed by the University of Tennessee as part of the
3#Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
4#project funded by the US National Science Foundation.
5#See the license text in license.txt
6#copyright 2008, University of Tennessee
7######################################################################
8
9## TODO: Need test,and check Gaussian averaging
10import numpy
11import math
12
13## Singular point
14SIGMA_ZERO = 1.0e-010
15## Limit of how many sigmas to be covered for the Gaussian smearing
16# default: 2.5 to cover 98.7% of Gaussian
17LIMIT = 3.0
18## Defaults
19R_BIN = {'Xhigh':10, 'High':5,'Med':5,'Low':3}
20PHI_BIN ={'Xhigh':20,'High':12,'Med':6,'Low':4}   
21
22class Smearer2D:
23    """
24    Gaussian Q smearing class for SANS 2d data
25    """
26     
27    def __init__(self, data=None, model=None, index=None, 
28                 limit=LIMIT, accuracy='Low', coords='polar', engine='python'):
29        """
30        Assumption: equally spaced bins in dq_r, dq_phi space.
31       
32        :param data: 2d data used to set the smearing parameters
33        :param model: model function
34        :param index: 1d array with len(data) to define the range of the calculation: elements are given as True or False
35        :param nr: number of bins in dq_r-axis
36        :param nphi: number of bins in dq_phi-axis
37        :param coord: coordinates [string], 'polar' or 'cartesian'
38        :param engine: engine name [string]; 'c' or 'numpy'
39        """
40        ## data
41        self.data = data
42        ## model
43        self.model = model
44        ## Accuracy: Higher stands for more sampling points in both directions of r and phi.
45        self.accuracy = accuracy
46        ## number of bins in r axis for over-sampling
47        self.nr = R_BIN
48        ## number of bins in phi axis for over-sampling
49        self.nphi = PHI_BIN
50        ## maximum nsigmas
51        self.limit = limit
52        self.index = index
53        self.coords = coords
54        self.smearer = True
55        self._engine = engine
56       
57       
58    def get_data(self):   
59        """
60        get qx_data, qy_data, dqx_data,dqy_data,and calculate phi_data=arctan(qx_data/qy_data)
61        """
62        if self.data == None or self.data.__class__.__name__ == 'Data1D':
63            return None
64        if self.data.dqx_data == None or self.data.dqy_data == None:
65            return None
66        self.qx_data = self.data.qx_data[self.index]
67        self.qy_data = self.data.qy_data[self.index]
68        self.q_data = self.data.q_data[self.index]
69        # Here dqx and dqy mean dq_parr and dq_perp
70        self.dqx_data = self.data.dqx_data[self.index]
71        self.dqy_data = self.data.dqy_data[self.index]
72        self.phi_data = numpy.arctan(self.qx_data/self.qy_data)
73        ## Remove singular points if exists
74        self.dqx_data[self.dqx_data<SIGMA_ZERO]=SIGMA_ZERO
75        self.dqy_data[self.dqy_data<SIGMA_ZERO]=SIGMA_ZERO
76        return True
77   
78    def set_accuracy(self, accuracy='Low'):   
79        """
80        Set accuracy.
81       
82        :param accuracy:  string
83        """
84        self.accuracy = accuracy
85
86    def set_smearer(self, smearer=True):
87        """
88        Set whether or not smearer will be used
89       
90        :param smearer: smear object
91       
92        """
93        self.smearer = smearer
94       
95    def set_data(self, data=None):   
96        """
97        Set data.
98       
99        :param data: DataLoader.Data_info type
100        """
101        self.data = data
102 
103           
104    def set_model(self, model=None):   
105        """
106        Set model.
107       
108        :param model: sans.models instance
109        """
110        self.model = model 
111           
112    def set_index(self, index=None):   
113        """
114        Set index.
115       
116        :param index: 1d arrays
117        """
118        self.index = index       
119   
120    def get_value(self):
121        """
122        Over sampling of r_nbins times phi_nbins, calculate Gaussian weights,
123        then find smeared intensity
124        """   
125        valid = self.get_data()
126        if valid == None:
127            return valid
128        # all zero values of dq
129        if numpy.all(numpy.fabs(self.dqx_data <= 1.1e-10)) and \
130                        numpy.all(numpy.fabs(self.dqy_data <= 1.1e-10)):
131            self.smearer = False
132 
133        if self.smearer == False:
134            return self.model.evalDistribution([self.qx_data, self.qy_data]) 
135
136        nr = self.nr[self.accuracy]
137        nphi = self.nphi[self.accuracy]
138        # Number of bins in the dqr direction (polar coordinate of dqx and dqy)
139        bin_size = self.limit / nr
140        # Total number of bins = # of bins
141        # in dq_r-direction times # of bins in dq_phi-direction
142        n_bins = nr * nphi
143        # data length in the range of self.index
144        len_data = len(self.qx_data)
145        len_datay = len(self.qy_data)
146
147        if self._engine == 'c' and self.coords == 'polar':
148            try:
149                import sans_extension.smearer2d_helper as smearer2dc
150                smearc = smearer2dc.new_Smearer_helper(self.qx_data, self.qy_data,
151                                              self.dqx_data, self.dqy_data,
152                                              self.limit, nr, nphi, int(len_data))
153                weight_res = numpy.zeros(nr * nphi )
154                qx_res = numpy.zeros(nr * nphi * int(len_data))
155                qy_res = numpy.zeros(nr * nphi * int(len_data))
156                smearer2dc.smearer2d_helper(smearc,weight_res, qx_res, qy_res)
157            except:
158                raise 
159        else:
160            # Mean values of dqr at each bins ,starting from the half of bin size
161            r = bin_size / 2.0 + numpy.arange(nr) * bin_size
162            # mean values of qphi at each bines
163            phi = numpy.arange(nphi)
164            dphi = phi * 2.0 * math.pi / nphi
165            dphi = dphi.repeat(nr)
166   
167            ## Transform to polar coordinate,
168            #  and set dphi at each data points ; 1d array
169            dphi = dphi.repeat(len_data)
170            q_phi = self.qy_data / self.qx_data
171           
172            # Starting angle is different between polar and cartesian coordinates.
173            #if self.coords != 'polar':
174            #    dphi += numpy.arctan( q_phi * self.dqx_data/ \
175            #                     self.dqy_data).repeat(n_bins).reshape(len_data,\
176            #                                        n_bins).transpose().flatten()
177   
178            # The angle (phi) of the original q point
179            q_phi = numpy.arctan(q_phi).repeat(n_bins).reshape(len_data,\
180                                                    n_bins).transpose().flatten()
181            ## Find Gaussian weight for each dq bins: The weight depends only
182            #  on r-direction (The integration may not need)
183            weight_res = numpy.exp(-0.5 * ((r - bin_size / 2.0) * \
184                                        (r - bin_size / 2.0)))- \
185                                        numpy.exp(-0.5 * ((r + bin_size / 2.0 ) *\
186                                        (r + bin_size / 2.0)))
187            # No needs of normalization here.
188            #weight_res /= numpy.sum(weight_res)
189            weight_res = weight_res.repeat(nphi).reshape(nr, nphi)
190   
191            weight_res = weight_res.transpose().flatten()
192           
193            ## Set dr for all dq bins for averaging
194            dr = r.repeat(nphi).reshape(nr,nphi).transpose().flatten()
195            ## Set dqr for all data points
196            dqx = numpy.outer(dr,self.dqx_data).flatten()
197            dqy = numpy.outer(dr,self.dqy_data).flatten()
198   
199            qx = self.qx_data.repeat(n_bins).reshape(len_data,\
200                                                 n_bins).transpose().flatten()
201            qy = self.qy_data.repeat(n_bins).reshape(len_data,\
202                                                 n_bins).transpose().flatten()
203   
204            # The polar needs rotation by -q_phi
205            if self.coords == 'polar':
206                q_r = numpy.sqrt(qx * qx + qy * qy)
207                qx_res = ((dqx*numpy.cos(dphi) + q_r) * numpy.cos(-q_phi) +\
208                               dqy*numpy.sin(dphi) * numpy.sin(-q_phi))
209                qy_res = (-(dqx*numpy.cos(dphi) + q_r) * numpy.sin(-q_phi) +\
210                               dqy*numpy.sin(dphi) * numpy.cos(-q_phi))
211            else:
212                qx_res = qx +  dqx*numpy.cos(dphi)
213                qy_res = qy +  dqy*numpy.sin(dphi)
214           
215        ## Evaluate all points
216        val = self.model.evalDistribution([qx_res, qy_res]) 
217        ## Reshape into 2d array to use numpy weighted averaging
218        value_res= val.reshape(n_bins,len(self.qx_data))
219        ## Averaging with Gaussian weighting: normalization included.
220        value =numpy.average(value_res,axis=0, weights=weight_res)
221        ## Return the smeared values in the range of self.index
222        return value
223   
224if __name__ == '__main__':
225    ## Test w/ 2D linear function
226    x = 0.001*numpy.arange(1,11)
227    dx = numpy.ones(len(x))*0.0003
228    y = 0.001*numpy.arange(1,11)
229    dy = numpy.ones(len(x))*0.001
230    z = numpy.ones(10)
231    dz = numpy.sqrt(z)
232   
233    from DataLoader import Data2D
234    #for i in range(10): print i, 0.001 + i*0.008/9.0
235    #for i in range(100): print i, int(math.floor( (i/ (100/9.0)) ))
236    out = Data2D()
237    out.data = z
238    out.qx_data = x
239    out.qy_data = y
240    out.dqx_data = dx
241    out.dqy_data = dy
242    out.q_data = numpy.sqrt(dx * dx + dy * dy)
243    index = numpy.ones(len(x), dtype = bool)
244    out.mask = index
245    from sans.models.LineModel import LineModel
246    model = LineModel()
247    model.setParam("A", 0)
248
249    smear = Smearer2D(out,model,index)
250    #smear.set_accuracy('Xhigh')
251    value = smear.get_value()
252    ## All data are ones, so the smeared should also be ones.
253    print "Data length =",len(value)
254    print " 2D linear function, I = 0 + 1*qy"
255    text = " Gaussian weighted averaging on a 2D linear function will "
256    text += "provides the results same as without the averaging."
257    print text
258    print "qx_data", "qy_data", "I_nonsmear", "I_smeared"
259    for ind in range(len(value)):
260        print x[ind],y[ind],model.evalDistribution([x,y])[ind], value[ind]
261       
262"""
263for i in range(len(qx_res)/(128*128)):
264    k = i * 128*128 +64
265
266    print qx_res[k]-qqx[k], qy_res[k]-qqy[k]
267print qqx[64],qqy[64]
268""" 
269"""   
270if __name__ == '__main__':
271    ## Another Test w/ constant function
272    x = 0.001*numpy.arange(1,11)
273    dx = numpy.ones(len(x))*0.001
274    y = 0.001*numpy.arange(1,11)
275    dy = numpy.ones(len(x))*0.001
276    z = numpy.ones(10)
277    dz = numpy.sqrt(z)
278   
279    from DataLoader import Data2D
280    #for i in range(10): print i, 0.001 + i*0.008/9.0
281    #for i in range(100): print i, int(math.floor( (i/ (100/9.0)) ))
282    out = Data2D()
283    out.data = z
284    out.qx_data = x
285    out.qy_data = y
286    out.dqx_data = dx
287    out.dqy_data = dy
288    index = numpy.ones(len(x), dtype = bool)
289    out.mask = index
290    from sans.models.Constant import Constant
291    model = Constant()
292
293    value = Smearer2D(out,model,index).get_value()
294    ## All data are ones, so the smeared values should also be ones.
295    print "Data length =",len(value), ", Data=",value
296"""   
Note: See TracBrowser for help on using the repository browser.