[346bc88] | 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 | """ |
---|
| 7 | from __future__ import division |
---|
| 8 | |
---|
| 9 | import numpy as np |
---|
| 10 | from numpy import pi, cos, sin, sqrt |
---|
| 11 | |
---|
| 12 | from .resolution import Resolution |
---|
| 13 | |
---|
| 14 | ## Singular point |
---|
| 15 | SIGMA_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 |
---|
| 18 | NSIGMA = 3.0 |
---|
| 19 | ## Defaults |
---|
[cd8dde1] | 20 | NR = {'xhigh':10, 'high':5, 'med':5, 'low':3} |
---|
| 21 | NPHI ={'xhigh':20, 'high':12, 'med':6, 'low':4} |
---|
[346bc88] | 22 | |
---|
| 23 | class 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 |
---|
[cd8dde1] | 43 | self.nr = NR[accuracy.lower()] |
---|
[346bc88] | 44 | ## number of bins in phi axis for over-sampling |
---|
[cd8dde1] | 45 | self.nphi = NPHI[accuracy.lower()] |
---|
[346bc88] | 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) |
---|
[3e6aaad] | 94 | bin_size = self.nsigma / nr |
---|
[346bc88] | 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 | """ |
---|
| 171 | if __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 |
---|
[9404dd3] | 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)) ))) |
---|
[346bc88] | 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. |
---|
[9404dd3] | 200 | print("Data length =", len(value)) |
---|
| 201 | print(" 2D linear function, I = 0 + 1*qy") |
---|
[346bc88] | 202 | text = " Gaussian weighted averaging on a 2D linear function will " |
---|
| 203 | text += "provides the results same as without the averaging." |
---|
[9404dd3] | 204 | print(text) |
---|
| 205 | print("qx_data", "qy_data", "I_nonsmear", "I_smeared") |
---|
[346bc88] | 206 | for ind in range(len(value)): |
---|
[9404dd3] | 207 | print(x[ind], y[ind], model.evalDistribution([x, y])[ind], value[ind]) |
---|
[346bc88] | 208 | |
---|
| 209 | |
---|
| 210 | if __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 |
---|
[9404dd3] | 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)) ))) |
---|
[346bc88] | 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. |
---|
[9404dd3] | 235 | print("Data length =",len(value), ", Data=",value) |
---|
[346bc88] | 236 | """ |
---|