source: sasmodels/sasmodels/resolution2d.py @ cf404cb

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since cf404cb was 9404dd3, checked in by Paul Kienzle <pkienzle@…>, 8 years ago

python 3.x support

  • Property mode set to 100644
File size: 8.5 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"""
7from __future__ import division
8
9import numpy as np
10from numpy import pi, cos, sin, sqrt
11
12from .resolution import Resolution
13
14## Singular point
15SIGMA_ZERO = 1.0e-010
16## Limit of how many sigmas to be covered for the Gaussian smearing
17# default: 2.5 to cover 98.7% of Gaussian
18NSIGMA = 3.0
19## Defaults
20NR = {'xhigh':10, 'high':5, 'med':5, 'low':3}
21NPHI ={'xhigh':20, 'high':12, 'med':6, 'low':4}
22
23class Pinhole2D(Resolution):
24    """
25    Gaussian Q smearing class for SAS 2d data
26    """
27     
28    def __init__(self, data=None, index=None,
29                 nsigma=NSIGMA, accuracy='Low', coords='polar'):
30        """
31        Assumption: equally spaced bins in dq_r, dq_phi space.
32       
33        :param data: 2d data used to set the smearing parameters
34        :param index: 1d array with len(data) to define the range
35         of the calculation: elements are given as True or False
36        :param nr: number of bins in dq_r-axis
37        :param nphi: number of bins in dq_phi-axis
38        :param coord: coordinates [string], 'polar' or 'cartesian'
39        """
40        ## Accuracy: Higher stands for more sampling points in both directions
41        ## of r and phi.
42        ## number of bins in r axis for over-sampling
43        self.nr = NR[accuracy.lower()]
44        ## number of bins in phi axis for over-sampling
45        self.nphi = NPHI[accuracy.lower()]
46        ## maximum nsigmas
47        self.nsigma= nsigma
48        self.coords = coords
49        self._init_data(data, index)
50
51    def _init_data(self, data, index):
52        """
53        Get qx_data, qy_data, dqx_data,dqy_data,
54        and calculate phi_data=arctan(qx_data/qy_data)
55        """
56        # TODO: maybe don't need to hold copy of qx,qy,dqx,dqy,data,index
57        # just need q_calc and weights
58        self.data = data
59        self.index = index
60
61        self.qx_data = data.qx_data[index]
62        self.qy_data = data.qy_data[index]
63        self.q_data = data.q_data[index]
64
65        dqx = getattr(data, 'dqx_data', None)
66        dqy = getattr(data, 'dqy_data', None)
67        if dqx is not None and dqy is not None:
68            # Here dqx and dqy mean dq_parr and dq_perp
69            self.dqx_data = dqx[index]
70            self.dqy_data = dqy[index]
71            ## Remove singular points if exists
72            self.dqx_data[self.dqx_data < SIGMA_ZERO] = SIGMA_ZERO
73            self.dqy_data[self.dqy_data < SIGMA_ZERO] = SIGMA_ZERO
74            qx_calc, qy_calc, weights = self._calc_res()
75            self.q_calc = [qx_calc, qy_calc]
76            self.q_calc_weights = weights
77        else:
78            # No resolution information
79            self.dqx_data = self.dqy_data = None
80            self.q_calc = [self.qx_data, self.qy_data]
81            self.q_calc_weights = None
82
83        #self.phi_data = np.arctan(self.qx_data / self.qy_data)
84
85    def _calc_res(self):
86        """
87        Over sampling of r_nbins times phi_nbins, calculate Gaussian weights,
88        then find smeared intensity
89        """   
90        nr, nphi = self.nr, self.nphi
91        # Total number of bins = # of bins
92        nbins = nr * nphi
93        # Number of bins in the dqr direction (polar coordinate of dqx and dqy)
94        bin_size = self.nsigma / nr
95        # in dq_r-direction times # of bins in dq_phi-direction
96        # data length in the range of self.index
97        nq = len(self.qx_data)
98
99        # Mean values of dqr at each bins
100        # starting from the half of bin size
101        r = bin_size / 2.0 + np.arange(nr) * bin_size
102        # mean values of qphi at each bines
103        phi = np.arange(nphi)
104        dphi = phi * 2.0 * pi / nphi
105        dphi = dphi.repeat(nr)
106
107        ## Transform to polar coordinate,
108        #  and set dphi at each data points ; 1d array
109        dphi = dphi.repeat(nq)
110        q_phi = self.qy_data / self.qx_data
111
112        # Starting angle is different between polar
113        #  and cartesian coordinates.
114        #if self.coords != 'polar':
115        #    dphi += np.arctan( q_phi * self.dqx_data/ \
116        #                  self.dqy_data).repeat(nbins).reshape(nq,\
117        #                                nbins).transpose().flatten()
118
119        # The angle (phi) of the original q point
120        q_phi = np.arctan(q_phi).repeat(nbins)\
121            .reshape(nq, nbins).transpose().flatten()
122        ## Find Gaussian weight for each dq bins: The weight depends only
123        #  on r-direction (The integration may not need)
124        weight_res = (np.exp(-0.5 * (r - bin_size / 2.0)**2)  -
125                      np.exp(-0.5 * (r + bin_size / 2.0)**2))
126        # No needs of normalization here.
127        #weight_res /= np.sum(weight_res)
128        weight_res = weight_res.repeat(nphi).reshape(nr, nphi)
129        weight_res = weight_res.transpose().flatten()
130
131        ## Set dr for all dq bins for averaging
132        dr = r.repeat(nphi).reshape(nr, nphi).transpose().flatten()
133        ## Set dqr for all data points
134        dqx = np.outer(dr, self.dqx_data).flatten()
135        dqy = np.outer(dr, self.dqy_data).flatten()
136
137        qx = self.qx_data.repeat(nbins)\
138            .reshape(nq, nbins).transpose().flatten()
139        qy = self.qy_data.repeat(nbins)\
140            .reshape(nq, nbins).transpose().flatten()
141
142        # The polar needs rotation by -q_phi
143        if self.coords == 'polar':
144            q_r = sqrt(qx**2 + qy**2)
145            qx_res = ( (dqx*cos(dphi) + q_r) * cos(-q_phi) +
146                           dqy*sin(dphi) * sin(-q_phi))
147            qy_res = (-(dqx*cos(dphi) + q_r) * sin(-q_phi) +
148                           dqy*sin(dphi) * cos(-q_phi))
149        else:
150            qx_res = qx +  dqx*cos(dphi)
151            qy_res = qy +  dqy*sin(dphi)
152
153
154        return qx_res, qy_res, weight_res
155
156    def apply(self, theory):
157        if self.q_calc_weights is not None:
158            # TODO: interpolate rather than recomputing all the different qx,qy
159            # Resolution needs to be applied
160            nq, nbins = len(self.qx_data), self.nr * self.nphi
161            ## Reshape into 2d array to use np weighted averaging
162            theory = np.reshape(theory, (nbins, nq))
163            ## Averaging with Gaussian weighting: normalization included.
164            value =np.average(theory, axis=0, weights=self.q_calc_weights)
165            ## Return the smeared values in the range of self.index
166            return value
167        else:
168            return theory
169
170"""
171if __name__ == '__main__':
172    ## Test w/ 2D linear function
173    x = 0.001*np.arange(1, 11)
174    dx = np.ones(len(x))*0.0003
175    y = 0.001*np.arange(1, 11)
176    dy = np.ones(len(x))*0.001
177    z = np.ones(10)
178    dz = sqrt(z)
179
180    from sas.dataloader import Data2D
181    #for i in range(10): print(i, 0.001 + i*0.008/9.0)
182    #for i in range(100): print(i, int(math.floor( (i/ (100/9.0)) )))
183    out = Data2D()
184    out.data = z
185    out.qx_data = x
186    out.qy_data = y
187    out.dqx_data = dx
188    out.dqy_data = dy
189    out.q_data = sqrt(dx * dx + dy * dy)
190    index = np.ones(len(x), dtype = bool)
191    out.mask = index
192    from sas.models.LineModel import LineModel
193    model = LineModel()
194    model.setParam("A", 0)
195
196    smear = Smearer2D(out, model, index)
197    #smear.set_accuracy('Xhigh')
198    value = smear.get_value()
199    ## All data are ones, so the smeared should also be ones.
200    print("Data length =", len(value))
201    print(" 2D linear function, I = 0 + 1*qy")
202    text = " Gaussian weighted averaging on a 2D linear function will "
203    text += "provides the results same as without the averaging."
204    print(text)
205    print("qx_data", "qy_data", "I_nonsmear", "I_smeared")
206    for ind in range(len(value)):
207        print(x[ind], y[ind], model.evalDistribution([x, y])[ind], value[ind])
208
209
210if __name__ == '__main__':
211    ## Another Test w/ constant function
212    x = 0.001*np.arange(1,11)
213    dx = np.ones(len(x))*0.001
214    y = 0.001*np.arange(1,11)
215    dy = np.ones(len(x))*0.001
216    z = np.ones(10)
217    dz = sqrt(z)
218
219    from DataLoader import Data2D
220    #for i in range(10): print(i, 0.001 + i*0.008/9.0)
221    #for i in range(100): print(i, int(math.floor( (i/ (100/9.0)) )))
222    out = Data2D()
223    out.data = z
224    out.qx_data = x
225    out.qy_data = y
226    out.dqx_data = dx
227    out.dqy_data = dy
228    index = np.ones(len(x), dtype = bool)
229    out.mask = index
230    from sas.models.Constant import Constant
231    model = Constant()
232
233    value = Smearer2D(out,model,index).get_value()
234    ## All data are ones, so the smeared values should also be ones.
235    print("Data length =",len(value), ", Data=",value)
236"""
Note: See TracBrowser for help on using the repository browser.