source: sasview/src/sas/models/smearing_2d.py @ 827484cf

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 827484cf was fd5ac0d, checked in by krzywon, 10 years ago

I have completed the removal of all SANS references.
I will build, run, and run all unit tests before pushing.

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