source: sasmodels/sasmodels/resolution2d.py @ 7ae2b7f

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

still more linting; ignore numpy types

  • Property mode set to 100644
File size: 9.1 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  # type: ignore
10from numpy import pi, cos, sin, sqrt  # type: ignore
11
12from . import resolution
13from .resolution import Resolution
14
15## Singular point
16SIGMA_ZERO = 1.0e-010
17## Limit of how many sigmas to be covered for the Gaussian smearing
18# default: 2.5 to cover 98.7% of Gaussian
19NSIGMA = 3.0
20## Defaults
21NR = {'xhigh':10, 'high':5, 'med':5, 'low':3}
22NPHI = {'xhigh':20, 'high':12, 'med':6, 'low':4}
23
24## Defaults
25N_SLIT_PERP = {'xhigh':1000, 'high':500, 'med':200, 'low':50}
26N_SLIT_PERP_DOC = ", ".join("%s=%d"%(name,value) for value,name in
27                            sorted((2*v+1,k) for k,v in N_SLIT_PERP.items()))
28
29class Pinhole2D(Resolution):
30    """
31    Gaussian Q smearing class for SAS 2d data
32    """
33
34    def __init__(self, data=None, index=None,
35                 nsigma=NSIGMA, accuracy='Low', coords='polar'):
36        """
37        Assumption: equally spaced bins in dq_r, dq_phi space.
38
39        :param data: 2d data used to set the smearing parameters
40        :param index: 1d array with len(data) to define the range
41         of the calculation: elements are given as True or False
42        :param nr: number of bins in dq_r-axis
43        :param nphi: number of bins in dq_phi-axis
44        :param coord: coordinates [string], 'polar' or 'cartesian'
45        """
46        ## Accuracy: Higher stands for more sampling points in both directions
47        ## of r and phi.
48        ## number of bins in r axis for over-sampling
49        self.nr = NR[accuracy.lower()]
50        ## number of bins in phi axis for over-sampling
51        self.nphi = NPHI[accuracy.lower()]
52        ## maximum nsigmas
53        self.nsigma = nsigma
54        self.coords = coords
55        self._init_data(data, index)
56
57    def _init_data(self, data, index):
58        """
59        Get qx_data, qy_data, dqx_data,dqy_data,
60        and calculate phi_data=arctan(qx_data/qy_data)
61        """
62        # TODO: maybe don't need to hold copy of qx,qy,dqx,dqy,data,index
63        # just need q_calc and weights
64        self.data = data
65        self.index = index
66
67        self.qx_data = data.qx_data[index]
68        self.qy_data = data.qy_data[index]
69        self.q_data = data.q_data[index]
70
71        dqx = getattr(data, 'dqx_data', None)
72        dqy = getattr(data, 'dqy_data', None)
73        if dqx is not None and dqy is not None:
74            # Here dqx and dqy mean dq_parr and dq_perp
75            self.dqx_data = dqx[index]
76            self.dqy_data = dqy[index]
77            ## Remove singular points if exists
78            self.dqx_data[self.dqx_data < SIGMA_ZERO] = SIGMA_ZERO
79            self.dqy_data[self.dqy_data < SIGMA_ZERO] = SIGMA_ZERO
80            qx_calc, qy_calc, weights = self._calc_res()
81            self.q_calc = [qx_calc, qy_calc]
82            self.q_calc_weights = weights
83        else:
84            # No resolution information
85            self.dqx_data = self.dqy_data = None
86            self.q_calc = [self.qx_data, self.qy_data]
87            self.q_calc_weights = None
88
89        #self.phi_data = np.arctan(self.qx_data / self.qy_data)
90
91    def _calc_res(self):
92        """
93        Over sampling of r_nbins times phi_nbins, calculate Gaussian weights,
94        then find smeared intensity
95        """
96        nr, nphi = self.nr, self.nphi
97        # Total number of bins = # of bins
98        nbins = nr * nphi
99        # Number of bins in the dqr direction (polar coordinate of dqx and dqy)
100        bin_size = self.nsigma / nr
101        # in dq_r-direction times # of bins in dq_phi-direction
102        # data length in the range of self.index
103        nq = len(self.qx_data)
104
105        # Mean values of dqr at each bins
106        # starting from the half of bin size
107        r = bin_size / 2.0 + np.arange(nr) * bin_size
108        # mean values of qphi at each bines
109        phi = np.arange(nphi)
110        dphi = phi * 2.0 * pi / nphi
111        dphi = dphi.repeat(nr)
112
113        ## Transform to polar coordinate,
114        #  and set dphi at each data points ; 1d array
115        dphi = dphi.repeat(nq)
116        q_phi = self.qy_data / self.qx_data
117
118        # Starting angle is different between polar
119        #  and cartesian coordinates.
120        #if self.coords != 'polar':
121        #    dphi += np.arctan( q_phi * self.dqx_data/ \
122        #                  self.dqy_data).repeat(nbins).reshape(nq,\
123        #                                nbins).transpose().flatten()
124
125        # The angle (phi) of the original q point
126        q_phi = np.arctan(q_phi).repeat(nbins)\
127            .reshape(nq, nbins).transpose().flatten()
128        ## Find Gaussian weight for each dq bins: The weight depends only
129        #  on r-direction (The integration may not need)
130        weight_res = (np.exp(-0.5 * (r - bin_size / 2.0)**2)  -
131                      np.exp(-0.5 * (r + bin_size / 2.0)**2))
132        # No needs of normalization here.
133        #weight_res /= np.sum(weight_res)
134        weight_res = weight_res.repeat(nphi).reshape(nr, nphi)
135        weight_res = weight_res.transpose().flatten()
136
137        ## Set dr for all dq bins for averaging
138        dr = r.repeat(nphi).reshape(nr, nphi).transpose().flatten()
139        ## Set dqr for all data points
140        dqx = np.outer(dr, self.dqx_data).flatten()
141        dqy = np.outer(dr, self.dqy_data).flatten()
142
143        qx = self.qx_data.repeat(nbins)\
144            .reshape(nq, nbins).transpose().flatten()
145        qy = self.qy_data.repeat(nbins)\
146            .reshape(nq, nbins).transpose().flatten()
147
148        # The polar needs rotation by -q_phi
149        if self.coords == 'polar':
150            q_r = sqrt(qx**2 + qy**2)
151            qx_res = ((dqx*cos(dphi) + q_r) * cos(-q_phi)
152                      + dqy*sin(dphi) * sin(-q_phi))
153            qy_res = (-(dqx*cos(dphi) + q_r) * sin(-q_phi)
154                      + dqy*sin(dphi) * cos(-q_phi))
155        else:
156            qx_res = qx + dqx*cos(dphi)
157            qy_res = qy + dqy*sin(dphi)
158
159
160        return qx_res, qy_res, weight_res
161
162    def apply(self, theory):
163        if self.q_calc_weights is not None:
164            # TODO: interpolate rather than recomputing all the different qx,qy
165            # Resolution needs to be applied
166            nq, nbins = len(self.qx_data), self.nr * self.nphi
167            ## Reshape into 2d array to use np weighted averaging
168            theory = np.reshape(theory, (nbins, nq))
169            ## Averaging with Gaussian weighting: normalization included.
170            value = np.average(theory, axis=0, weights=self.q_calc_weights)
171            ## Return the smeared values in the range of self.index
172            return value
173        else:
174            return theory
175
176
177class Slit2D(Resolution):
178    """
179    Slit aperture with resolution function on an oriented sample.
180
181    *q* points at which the data is measured.
182
183    *qx_width* slit width in qx
184
185    *qy_width* slit height in qy; current implementation requires a fixed
186    qy_width for all q points.
187
188    *q_calc* is the list of q points to calculate, or None if this
189    should be estimated from the *q* and *qx_width*.
190
191    *accuracy* determines the number of *qy* points to compute for each *q*.
192    The values are stored in sasmodels.resolution2d.N_SLIT_PERP.  The default
193    values are: %s
194    """
195    __doc__ = __doc__%N_SLIT_PERP_DOC
196    def __init__(self, q, qx_width, qy_width=0., q_calc=None, accuracy='low'):
197        # Remember what q and width was used even though we won't need them
198        # after the weight matrix is constructed
199        self.q, self.qx_width, self.qy_width = q, qx_width, qy_width
200
201        # Allow independent resolution on each qx point even though it is not
202        # needed in practice.  Set qy_width to the maximum qy width.
203        if np.isscalar(qx_width):
204            qx_width = np.ones(len(q))*qx_width
205        else:
206            qx_width = np.asarray(qx_width)
207        if not np.isscalar(qy_width):
208            qy_width = np.max(qy_width)
209
210        # Build grid of qx, qy points
211        if q_calc is not None:
212            qx_calc = np.sort(q_calc)
213        else:
214            qx_calc = resolution.pinhole_extend_q(q, qx_width, nsigma=3)
215        qy_min, qy_max = np.log10(np.min(q)), np.log10(qy_width)
216        qy_calc = np.logspace(qy_min, qy_max, N_SLIT_PERP[accuracy])
217        qy_calc = np.hstack((-qy_calc[::-1], 0, qy_calc))
218        self.q_calc = [v.flatten() for v in np.meshgrid(qx_calc, qy_calc)]
219        self.qx_calc, self.qy_calc = qx_calc, qy_calc
220        self.nx, self.ny = len(qx_calc), len(qy_calc)
221        self.dy = 2*qy_width/self.ny
222
223        # Build weight matrix for resolution integration
224        if np.any(qx_width > 0):
225            self.weights = resolution.pinhole_resolution(qx_calc, q,
226                    np.maximum(qx_width, resolution.MINIMUM_RESOLUTION))
227        elif len(qx_calc)==len(q) and np.all(qx_calc == q):
228            self.weights = None
229        else:
230            raise ValueError("Slit2D fails with q_calc != q")
231
232    def apply(self, theory):
233        Iq = np.trapz(theory.reshape(self.ny, self.nx), axis=0, x=self.qy_calc)
234        if self.weights is not None:
235            Iq = resolution.apply_resolution_matrix(self.weights, Iq)
236        return Iq
237
Note: See TracBrowser for help on using the repository browser.