Source code for sas.models.smearing_2d

"""
#This software was developed by the University of Tennessee as part of the
#Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
#project funded by the US National Science Foundation. 
#See the license text in license.txt
"""
import numpy
import math

## Singular point
SIGMA_ZERO = 1.0e-010
## Limit of how many sigmas to be covered for the Gaussian smearing
# default: 2.5 to cover 98.7% of Gaussian
LIMIT = 3.0
## Defaults
R_BIN = {'Xhigh':10, 'High':5, 'Med':5, 'Low':3}
PHI_BIN ={'Xhigh':20, 'High':12, 'Med':6, 'Low':4}   

[docs]class Smearer2D: """ Gaussian Q smearing class for SAS 2d data """ def __init__(self, data=None, model=None, index=None, limit=LIMIT, accuracy='Low', coords='polar', engine='c'): """ Assumption: equally spaced bins in dq_r, dq_phi space. :param data: 2d data used to set the smearing parameters :param model: model function :param index: 1d array with len(data) to define the range of the calculation: elements are given as True or False :param nr: number of bins in dq_r-axis :param nphi: number of bins in dq_phi-axis :param coord: coordinates [string], 'polar' or 'cartesian' :param engine: engine name [string]; 'c' or 'numpy' """ ## data self.data = data ## model self.model = model ## Accuracy: Higher stands for more sampling points in both directions ## of r and phi. self.accuracy = accuracy ## number of bins in r axis for over-sampling self.nr = R_BIN ## number of bins in phi axis for over-sampling self.nphi = PHI_BIN ## maximum nsigmas self.limit = limit self.index = index self.coords = coords self.smearer = True self._engine = engine self.qx_data = None self.qy_data = None self.q_data = None # dqx and dqy mean dq_parr and dq_perp self.dqx_data = None self.dqy_data = None self.phi_data = None
[docs] def get_data(self): """ Get qx_data, qy_data, dqx_data,dqy_data, and calculate phi_data=arctan(qx_data/qy_data) """ if self.data == None or self.data.__class__.__name__ == 'Data1D': return None if self.data.dqx_data == None or self.data.dqy_data == None: return None self.qx_data = self.data.qx_data[self.index] self.qy_data = self.data.qy_data[self.index] self.q_data = self.data.q_data[self.index] # Here dqx and dqy mean dq_parr and dq_perp self.dqx_data = self.data.dqx_data[self.index] self.dqy_data = self.data.dqy_data[self.index] self.phi_data = numpy.arctan(self.qx_data / self.qy_data) ## Remove singular points if exists self.dqx_data[self.dqx_data < SIGMA_ZERO] = SIGMA_ZERO self.dqy_data[self.dqy_data < SIGMA_ZERO] = SIGMA_ZERO return True
[docs] def set_accuracy(self, accuracy='Low'): """ Set accuracy. :param accuracy: string """ self.accuracy = accuracy
[docs] def set_smearer(self, smearer=True): """ Set whether or not smearer will be used :param smearer: smear object """ self.smearer = smearer
[docs] def set_data(self, data=None): """ Set data. :param data: DataLoader.Data_info type """ self.data = data
[docs] def set_model(self, model=None): """ Set model. :param model: sas.models instance """ self.model = model
[docs] def set_index(self, index=None): """ Set index. :param index: 1d arrays """ self.index = index
[docs] def get_value(self): """ Over sampling of r_nbins times phi_nbins, calculate Gaussian weights, then find smeared intensity """ valid = self.get_data() if valid == None: return valid # all zero values of dq if numpy.all(numpy.fabs(self.dqx_data <= 1.1e-10)) and \ numpy.all(numpy.fabs(self.dqy_data <= 1.1e-10)): self.smearer = False if self.smearer == False: return self.model.evalDistribution([self.qx_data, self.qy_data]) nr = self.nr[self.accuracy] nphi = self.nphi[self.accuracy] # Number of bins in the dqr direction (polar coordinate of dqx and dqy) bin_size = self.limit / nr # Total number of bins = # of bins # in dq_r-direction times # of bins in dq_phi-direction n_bins = nr * nphi # data length in the range of self.index len_data = len(self.qx_data) #len_datay = len(self.qy_data) if self._engine == 'c' and self.coords == 'polar': try: import sas.models.sas_extension.smearer2d_helper as smearer2dc smearc = smearer2dc.new_Smearer_helper(self.qx_data, self.qy_data, self.dqx_data, self.dqy_data, self.limit, nr, nphi, int(len_data)) weight_res = numpy.zeros(nr * nphi ) qx_res = numpy.zeros(nr * nphi * int(len_data)) qy_res = numpy.zeros(nr * nphi * int(len_data)) smearer2dc.smearer2d_helper(smearc, weight_res, qx_res, qy_res) except: raise else: # Mean values of dqr at each bins # starting from the half of bin size r = bin_size / 2.0 + numpy.arange(nr) * bin_size # mean values of qphi at each bines phi = numpy.arange(nphi) dphi = phi * 2.0 * math.pi / nphi dphi = dphi.repeat(nr) ## Transform to polar coordinate, # and set dphi at each data points ; 1d array dphi = dphi.repeat(len_data) q_phi = self.qy_data / self.qx_data # Starting angle is different between polar # and cartesian coordinates. #if self.coords != 'polar': # dphi += numpy.arctan( q_phi * self.dqx_data/ \ # self.dqy_data).repeat(n_bins).reshape(len_data,\ # n_bins).transpose().flatten() # The angle (phi) of the original q point q_phi = numpy.arctan(q_phi).repeat(n_bins).reshape(len_data,\ n_bins).transpose().flatten() ## Find Gaussian weight for each dq bins: The weight depends only # on r-direction (The integration may not need) weight_res = numpy.exp(-0.5 * ((r - bin_size / 2.0) * \ (r - bin_size / 2.0)))- \ numpy.exp(-0.5 * ((r + bin_size / 2.0 ) *\ (r + bin_size / 2.0))) # No needs of normalization here. #weight_res /= numpy.sum(weight_res) weight_res = weight_res.repeat(nphi).reshape(nr, nphi) weight_res = weight_res.transpose().flatten() ## Set dr for all dq bins for averaging dr = r.repeat(nphi).reshape(nr, nphi).transpose().flatten() ## Set dqr for all data points dqx = numpy.outer(dr, self.dqx_data).flatten() dqy = numpy.outer(dr, self.dqy_data).flatten() qx = self.qx_data.repeat(n_bins).reshape(len_data, \ n_bins).transpose().flatten() qy = self.qy_data.repeat(n_bins).reshape(len_data, \ n_bins).transpose().flatten() # The polar needs rotation by -q_phi if self.coords == 'polar': q_r = numpy.sqrt(qx * qx + qy * qy) qx_res = ((dqx*numpy.cos(dphi) + q_r) * numpy.cos(-q_phi) +\ dqy*numpy.sin(dphi) * numpy.sin(-q_phi)) qy_res = (-(dqx*numpy.cos(dphi) + q_r) * numpy.sin(-q_phi) +\ dqy*numpy.sin(dphi) * numpy.cos(-q_phi)) else: qx_res = qx + dqx*numpy.cos(dphi) qy_res = qy + dqy*numpy.sin(dphi) ## Evaluate all points val = self.model.evalDistribution([qx_res, qy_res]) ## Reshape into 2d array to use numpy weighted averaging value_res= val.reshape(n_bins, len(self.qx_data)) ## Averaging with Gaussian weighting: normalization included. value =numpy.average(value_res,axis=0, weights=weight_res) ## Return the smeared values in the range of self.index return value
""" if __name__ == '__main__': ## Test w/ 2D linear function x = 0.001*numpy.arange(1, 11) dx = numpy.ones(len(x))*0.0003 y = 0.001*numpy.arange(1, 11) dy = numpy.ones(len(x))*0.001 z = numpy.ones(10) dz = numpy.sqrt(z) from sas.dataloader import Data2D #for i in range(10): print i, 0.001 + i*0.008/9.0 #for i in range(100): print i, int(math.floor( (i/ (100/9.0)) )) out = Data2D() out.data = z out.qx_data = x out.qy_data = y out.dqx_data = dx out.dqy_data = dy out.q_data = numpy.sqrt(dx * dx + dy * dy) index = numpy.ones(len(x), dtype = bool) out.mask = index from sas.models.LineModel import LineModel model = LineModel() model.setParam("A", 0) smear = Smearer2D(out, model, index) #smear.set_accuracy('Xhigh') value = smear.get_value() ## All data are ones, so the smeared should also be ones. print "Data length =", len(value) print " 2D linear function, I = 0 + 1*qy" text = " Gaussian weighted averaging on a 2D linear function will " text += "provides the results same as without the averaging." print text print "qx_data", "qy_data", "I_nonsmear", "I_smeared" for ind in range(len(value)): print x[ind], y[ind], model.evalDistribution([x, y])[ind], value[ind] if __name__ == '__main__': ## Another Test w/ constant function x = 0.001*numpy.arange(1,11) dx = numpy.ones(len(x))*0.001 y = 0.001*numpy.arange(1,11) dy = numpy.ones(len(x))*0.001 z = numpy.ones(10) dz = numpy.sqrt(z) from DataLoader import Data2D #for i in range(10): print i, 0.001 + i*0.008/9.0 #for i in range(100): print i, int(math.floor( (i/ (100/9.0)) )) out = Data2D() out.data = z out.qx_data = x out.qy_data = y out.dqx_data = dx out.dqy_data = dy index = numpy.ones(len(x), dtype = bool) out.mask = index from sas.models.Constant import Constant model = Constant() value = Smearer2D(out,model,index).get_value() ## All data are ones, so the smeared values should also be ones. print "Data length =",len(value), ", Data=",value """