Changeset 8cb0692 in sasview for src/sas/sascalc/data_util


Ignore:
Timestamp:
Mar 19, 2016 1:14:31 PM (9 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, ESS_GUI, ESS_GUI_Docs, ESS_GUI_batch_fitting, ESS_GUI_bumps_abstraction, ESS_GUI_iss1116, ESS_GUI_iss879, ESS_GUI_iss959, ESS_GUI_opencl, ESS_GUI_ordering, ESS_GUI_sync_sascalc, costrafo411, magnetic_scatt, release-4.1.1, release-4.1.2, release-4.2.2, release_4.0.1, ticket-1009, ticket-1094-headless, ticket-1242-2d-resolution, ticket-1243, ticket-1249, ticket885, unittest-saveload
Children:
ed0427a
Parents:
3c53d7f
Message:

use resolution functions from sasmodels

Location:
src/sas/sascalc/data_util
Files:
2 deleted
1 edited

Legend:

Unmodified
Added
Removed

  • src/sas/sascalc/data_util/qsmearing.py

    rfc18690 r8cb0692  
    1313import logging 
    1414import sys 
    15 import sas.models.sas_extension.smearer as smearer 
    16 from sas.sascalc.data_util.smearing_2d import Smearer2D 
    1715 
    18 def smear_selection(data1D, model = None): 
     16from sasmodels.resolution import Slit1D, Pinhole1D 
     17from sasmodels.resolution2d import Pinhole2D 
     18 
     19def smear_selection(data, model = None): 
    1920    """ 
    2021    Creates the right type of smearer according  
    … …  
    3132    resolution smearing data.  
    3233     
    33     :param data1D: Data1D object 
     34    :param data: Data1D object 
    3435    :param model: sas.model instance 
    3536    """ 
    3637    # Sanity check. If we are not dealing with a SAS Data1D 
    3738    # object, just return None 
    38     if  data1D.__class__.__name__ not in ['Data1D', 'Theory1D']: 
    39         if data1D == None: 
     39    if  data.__class__.__name__ not in ['Data1D', 'Theory1D']: 
     40        if data == None: 
    4041            return None 
    41         elif data1D.dqx_data == None or data1D.dqy_data == None: 
     42        elif data.dqx_data == None or data.dqy_data == None: 
    4243            return None 
    43         return Smearer2D(data1D) 
     44        return Pinhole2D(data) 
    4445     
    45     if  not hasattr(data1D, "dx") and not hasattr(data1D, "dxl")\ 
    46          and not hasattr(data1D, "dxw"): 
     46    if  not hasattr(data, "dx") and not hasattr(data, "dxl")\ 
     47         and not hasattr(data, "dxw"): 
    4748        return None 
    4849     
    4950    # Look for resolution smearing data 
    5051    _found_resolution = False 
    51     if data1D.dx is not None and len(data1D.dx) == len(data1D.x): 
     52    if data.dx is not None and len(data.dx) == len(data.x): 
    5253         
    5354        # Check that we have non-zero data 
    54         if data1D.dx[0] > 0.0: 
     55        if data.dx[0] > 0.0: 
    5556            _found_resolution = True 
    5657            #print "_found_resolution",_found_resolution 
    … …  
    5859    # If we found resolution smearing data, return a QSmearer 
    5960    if _found_resolution == True: 
    60         return QSmearer(data1D, model) 
    61         #return pinhole_smear(data1D, model) 
     61         return pinhole_smear(data, model) 
    6262 
    6363    # Look for slit smearing data 
    6464    _found_slit = False 
    65     if data1D.dxl is not None and len(data1D.dxl) == len(data1D.x) \ 
    66         and data1D.dxw is not None and len(data1D.dxw) == len(data1D.x): 
     65    if data.dxl is not None and len(data.dxl) == len(data.x) \ 
     66        and data.dxw is not None and len(data.dxw) == len(data.x): 
    6767         
    6868        # Check that we have non-zero data 
    69         if data1D.dxl[0] > 0.0 or data1D.dxw[0] > 0.0: 
     69        if data.dxl[0] > 0.0 or data.dxw[0] > 0.0: 
    7070            _found_slit = True 
    7171         
    7272        # Sanity check: all data should be the same as a function of Q 
    73         for item in data1D.dxl: 
    74             if data1D.dxl[0] != item: 
     73        for item in data.dxl: 
     74            if data.dxl[0] != item: 
    7575                _found_resolution = False 
    7676                break 
    7777             
    78         for item in data1D.dxw: 
    79             if data1D.dxw[0] != item: 
     78        for item in data.dxw: 
     79            if data.dxw[0] != item: 
    8080                _found_resolution = False 
    8181                break 
    8282    # If we found slit smearing data, return a slit smearer 
    8383    if _found_slit == True: 
    84         #return SlitSmearer(data1D, model) 
    85         return slit_smear(data1D, model) 
     84        return slit_smear(data, model) 
    8685    return None 
    87              
    88  
    89 class _BaseSmearer(object): 
    90     """ 
    91         Base class for smearers 
    92     """ 
    93     def __init__(self): 
    94         self.nbins = 0 
    95         self.nbins_low = 0 
    96         self.nbins_high = 0 
    97         self._weights = None 
    98         ## Internal flag to keep track of C++ smearer initialization 
    99         self._init_complete = False 
    100         self._smearer = None 
    101         self.model = None 
    102         self.min = None 
    103         self.max = None 
    104         self.qvalues = [] 
    105          
    106     def __deepcopy__(self, memo=None): 
    107         """ 
    108         Return a valid copy of self. 
    109         Avoid copying the _smearer C object and force a matrix recompute 
    110         when the copy is used.   
    111         """ 
    112         result = _BaseSmearer() 
    113         result.nbins = self.nbins 
    114         return result 
    115  
    116     def _compute_matrix(self): 
    117         """ 
    118             Place holder for matrix computation  
    119         """ 
    120         return NotImplemented 
    121  
    122     def get_unsmeared_range(self, q_min=None, q_max=None): 
    123         """ 
    124             Place holder for method returning unsmeared range 
    125         """ 
    126         return NotImplemented 
    127      
    128     def get_bin_range(self, q_min=None, q_max=None): 
    129         """ 
    130          
    131         :param q_min: minimum q-value to smear 
    132         :param q_max: maximum q-value to smear 
    133           
    134         """ 
    135         # If this is the first time we call for smearing, 
    136         # initialize the C++ smearer object first 
    137         if not self._init_complete: 
    138             self._initialize_smearer() 
    139         if q_min == None: 
    140             q_min = self.min 
    141         if q_max == None: 
    142             q_max = self.max 
    143  
    144         _qmin_unsmeared, _qmax_unsmeared = self.get_unsmeared_range(q_min, 
    145                                                                      q_max) 
    146         _first_bin = None 
    147         _last_bin  = None 
    148  
    149         #step = (self.max - self.min) / (self.nbins - 1.0) 
    150         # Find the first and last bin number in all extrapolated and real data 
    151         try: 
    152             for i in range(self.nbins): 
    153                 q_i = smearer.get_q(self._smearer, i) 
    154                 if (q_i >= _qmin_unsmeared) and (q_i <= _qmax_unsmeared): 
    155                     # Identify first and last bin 
    156                     if _first_bin is None: 
    157                         _first_bin = i 
    158                     else: 
    159                         _last_bin  = i 
    160         except: 
    161             msg = "_BaseSmearer.get_bin_range: " 
    162             msg += " error getting range\n  %s" % sys.exc_value 
    163             raise RuntimeError, msg 
    164     
    165         #  Find the first and last bin number only in the real data 
    166         _first_bin, _last_bin = self._get_unextrapolated_bin( \ 
    167                                         _first_bin, _last_bin) 
    168  
    169         return _first_bin, _last_bin 
    170  
    171     def __call__(self, iq_in, first_bin = 0, last_bin = None): 
    172         """ 
    173         Perform smearing 
    174         """ 
    175         # If this is the first time we call for smearing, 
    176         # initialize the C++ smearer object first 
    177         if not self._init_complete: 
    178             self._initialize_smearer() 
    179  
    180         if last_bin is None or last_bin >= len(iq_in): 
    181             last_bin = len(iq_in) - 1 
    182         # Check that the first bin is positive 
    183         if first_bin < 0: 
    184             first_bin = 0 
    185  
    186         # With a model given, compute I for the extrapolated points and append 
    187         # to the iq_in 
    188         iq_in_temp = iq_in 
    189         if self.model != None: 
    190             temp_first, temp_last = self._get_extrapolated_bin( \ 
    191                                                         first_bin, last_bin) 
    192             if self.nbins_low > 0: 
    193                 iq_in_low = self.model.evalDistribution( \ 
    194                                     numpy.fabs(self.qvalues[0:self.nbins_low])) 
    195             iq_in_high = self.model.evalDistribution( \ 
    196                                             self.qvalues[(len(self.qvalues) - \ 
    197                                             self.nbins_high - 1):]) 
    198             # Todo: find out who is sending iq[last_poin] = 0. 
    199             if iq_in[len(iq_in) - 1] == 0: 
    200                 iq_in[len(iq_in) - 1] = iq_in_high[0] 
    201             # Append the extrapolated points to the data points 
    202             if self.nbins_low > 0: 
    203                 iq_in_temp = numpy.append(iq_in_low, iq_in) 
    204             if self.nbins_high > 0: 
    205                 iq_in_temp = numpy.append(iq_in_temp, iq_in_high[1:]) 
    206         else: 
    207             temp_first = first_bin 
    208             temp_last = last_bin 
    209             #iq_in_temp = iq_in 
    210  
    211         # Sanity check 
    212         if len(iq_in_temp) != self.nbins: 
    213             msg = "Invalid I(q) vector: inconsistent array " 
    214             msg += " length %d != %s" % (len(iq_in_temp), str(self.nbins)) 
    215             raise RuntimeError, msg 
    216  
    217         # Storage for smeared I(q)    
    218         iq_out = numpy.zeros(self.nbins) 
    219  
    220         smear_output = smearer.smear(self._smearer, iq_in_temp, iq_out, 
    221                                       #0, self.nbins - 1) 
    222                                       temp_first, temp_last) 
    223                                       #first_bin, last_bin) 
    224         if smear_output < 0: 
    225             msg = "_BaseSmearer: could not smear, code = %g" % smear_output 
    226             raise RuntimeError, msg 
    227  
    228         temp_first = first_bin + self.nbins_low 
    229         temp_last = self.nbins - self.nbins_high 
    230         out = iq_out[temp_first: temp_last] 
    231  
    232         return out 
    233      
    234     def _initialize_smearer(self): 
    235         """ 
    236             Place holder for initializing data smearer 
    237         """ 
    238         return NotImplemented 
    239              
    240      
    241     def _get_unextrapolated_bin(self, first_bin = 0, last_bin = 0): 
    242         """ 
    243         Get unextrapolated first bin and the last bin 
    244          
    245         : param first_bin: extrapolated first_bin 
    246         : param last_bin: extrapolated last_bin 
    247          
    248         : return fist_bin, last_bin: unextrapolated first and last bin 
    249         """ 
    250         #  For first bin 
    251         if first_bin <= self.nbins_low: 
    252             first_bin = 0 
    253         else: 
    254             first_bin = first_bin - self.nbins_low 
    255         # For last bin 
    256         if last_bin >= (self.nbins - self.nbins_high): 
    257             last_bin  = self.nbins - (self.nbins_high + self.nbins_low + 1) 
    258         elif last_bin >= self.nbins_low: 
    259             last_bin = last_bin - self.nbins_low 
    260         else: 
    261             last_bin = 0 
    262         return first_bin, last_bin 
    263      
    264     def _get_extrapolated_bin(self, first_bin = 0, last_bin = 0): 
    265         """ 
    266         Get extrapolated first bin and the last bin 
    267          
    268         : param first_bin: unextrapolated first_bin 
    269         : param last_bin: unextrapolated last_bin 
    270          
    271         : return first_bin, last_bin: extrapolated first and last bin 
    272         """ 
    273         #  For the first bin 
    274         # In the case that needs low extrapolation data 
    275         first_bin = 0 
    276         # For last bin 
    277         if last_bin >= self.nbins - (self.nbins_high + self.nbins_low + 1): 
    278             # In the case that needs higher q extrapolation data  
    279             last_bin = self.nbins - 1 
    280         else: 
    281             # In the case that doesn't need higher q extrapolation data  
    282             last_bin += self.nbins_low 
    283  
    284         return first_bin, last_bin 
    285          
    286 class _SlitSmearer(_BaseSmearer): 
    287     """ 
    288     Slit smearing for I(q) array 
    289     """ 
    290      
    291     def __init__(self, nbins=None, width=None, height=None, min=None, max=None): 
    292         """ 
    293         Initialization 
    294              
    295         :param iq: I(q) array [cm-1] 
    296         :param width: slit width [A-1] 
    297         :param height: slit height [A-1] 
    298         :param min: Q_min [A-1] 
    299         :param max: Q_max [A-1] 
    300          
    301         """ 
    302         _BaseSmearer.__init__(self) 
    303         ## Slit width in Q units 
    304         self.width  = width 
    305         ## Slit height in Q units 
    306         self.height = height 
    307         ## Q_min (Min Q-value for I(q)) 
    308         self.min    = min 
    309         ## Q_max (Max Q_value for I(q)) 
    310         self.max    = max 
    311         ## Number of Q bins  
    312         self.nbins  = nbins 
    313         ## Number of points used in the smearing computation 
    314         self.npts   = 3000 
    315         ## Smearing matrix 
    316         self._weights = None 
    317         self.qvalues  = None 
    318          
    319     def _initialize_smearer(self): 
    320         """ 
    321         Initialize the C++ smearer object. 
    322         This method HAS to be called before smearing 
    323         """ 
    324         #self._smearer = smearer.new_slit_smearer(self.width, 
    325         # self.height, self.min, self.max, self.nbins) 
    326         self._smearer = smearer.new_slit_smearer_with_q(self.width,  
    327                                                     self.height, self.qvalues) 
    328         self._init_complete = True 
    329  
    330     def get_unsmeared_range(self, q_min, q_max): 
    331         """ 
    332         Determine the range needed in unsmeared-Q to cover 
    333         the smeared Q range 
    334         """ 
    335         # Range used for input to smearing 
    336         _qmin_unsmeared = q_min 
    337         _qmax_unsmeared = q_max  
    338         try: 
    339             _qmin_unsmeared = self.min 
    340             _qmax_unsmeared = self.max 
    341         except: 
    342             logging.error("_SlitSmearer.get_bin_range: %s" % sys.exc_value) 
    343         return _qmin_unsmeared, _qmax_unsmeared 
    344  
    345 class SlitSmearer(_SlitSmearer): 
    346     """ 
    347     Adaptor for slit smearing class and SAS data 
    348     """ 
    349     def __init__(self, data1D, model = None): 
    350         """ 
    351         Assumption: equally spaced bins of increasing q-values. 
    352          
    353         :param data1D: data used to set the smearing parameters 
    354         """ 
    355         # Initialization from parent class 
    356         super(SlitSmearer, self).__init__() 
    357          
    358         ## Slit width 
    359         self.width = 0 
    360         self.nbins_low = 0 
    361         self.nbins_high = 0 
    362         self.model = model 
    363         if data1D.dxw is not None and len(data1D.dxw) == len(data1D.x): 
    364             self.width = data1D.dxw[0] 
    365             # Sanity check 
    366             for value in data1D.dxw: 
    367                 if value != self.width: 
    368                     msg = "Slit smearing parameters must " 
    369                     msg += " be the same for all data" 
    370                     raise RuntimeError, msg 
    371         ## Slit height 
    372         self.height = 0 
    373         if data1D.dxl is not None and len(data1D.dxl) == len(data1D.x): 
    374             self.height = data1D.dxl[0] 
    375             # Sanity check 
    376             for value in data1D.dxl: 
    377                 if value != self.height: 
    378                     msg = "Slit smearing parameters must be" 
    379                     msg += " the same for all data" 
    380                     raise RuntimeError, msg 
    381         # If a model is given, get the q extrapolation 
    382         if self.model == None: 
    383             data1d_x = data1D.x 
    384         else: 
    385             # Take larger sigma 
    386             if self.height > self.width: 
    387                 # The denominator (2.0) covers all the possible w^2 + h^2 range 
    388                 sigma_in = data1D.dxl / 2.0 
    389             elif self.width > 0: 
    390                 sigma_in = data1D.dxw / 2.0 
    391             else: 
    392                 sigma_in = [] 
    393  
    394             self.nbins_low, self.nbins_high, _, data1d_x = \ 
    395                                 get_qextrapolate(sigma_in, data1D.x) 
    396  
    397         ## Number of Q bins 
    398         self.nbins = len(data1d_x) 
    399         ## Minimum Q  
    400         self.min = min(data1d_x) 
    401         ## Maximum 
    402         self.max = max(data1d_x) 
    403         ## Q-values 
    404         self.qvalues = data1d_x 
    405          
    406  
    407 class _QSmearer(_BaseSmearer): 
    408     """ 
    409     Perform Gaussian Q smearing 
    410     """ 
    411          
    412     def __init__(self, nbins=None, width=None, min=None, max=None): 
    413         """ 
    414         Initialization 
    415          
    416         :param nbins: number of Q bins 
    417         :param width: array standard deviation in Q [A-1] 
    418         :param min: Q_min [A-1] 
    419         :param max: Q_max [A-1] 
    420         """ 
    421         _BaseSmearer.__init__(self) 
    422         ## Standard deviation in Q [A-1] 
    423         self.width = width 
    424         ## Q_min (Min Q-value for I(q)) 
    425         self.min = min 
    426         ## Q_max (Max Q_value for I(q)) 
    427         self.max = max 
    428         ## Number of Q bins  
    429         self.nbins = nbins 
    430         ## Smearing matrix 
    431         self._weights = None 
    432         self.qvalues  = None 
    433          
    434     def _initialize_smearer(self): 
    435         """ 
    436         Initialize the C++ smearer object. 
    437         This method HAS to be called before smearing 
    438         """ 
    439         #self._smearer = smearer.new_q_smearer(numpy.asarray(self.width), 
    440         # self.min, self.max, self.nbins) 
    441         self._smearer = smearer.new_q_smearer_with_q(numpy.asarray(self.width), 
    442                                                       self.qvalues) 
    443         self._init_complete = True 
    444          
    445     def get_unsmeared_range(self, q_min, q_max): 
    446         """ 
    447         Determine the range needed in unsmeared-Q to cover 
    448         the smeared Q range 
    449         Take 3 sigmas as the offset between smeared and unsmeared space 
    450         """ 
    451         # Range used for input to smearing 
    452         _qmin_unsmeared = q_min 
    453         _qmax_unsmeared = q_max  
    454         try: 
    455             offset = 3.0 * max(self.width) 
    456             _qmin_unsmeared = self.min#max([self.min, q_min - offset]) 
    457             _qmax_unsmeared = self.max#min([self.max, q_max + offset]) 
    458         except: 
    459             logging.error("_QSmearer.get_bin_range: %s" % sys.exc_value) 
    460         return _qmin_unsmeared, _qmax_unsmeared 
    461          
    462      
    463 class QSmearer(_QSmearer): 
    464     """ 
    465     Adaptor for Gaussian Q smearing class and SAS data 
    466     """ 
    467     def __init__(self, data1D, model = None): 
    468         """ 
    469         Assumption: equally spaced bins of increasing q-values. 
    470          
    471         :param data1D: data used to set the smearing parameters 
    472         """ 
    473         # Initialization from parent class 
    474         super(QSmearer, self).__init__() 
    475         data1d_x = [] 
    476         self.nbins_low = 0 
    477         self.nbins_high = 0 
    478         self.model = model 
    479         ## Resolution 
    480         #self.width = numpy.zeros(len(data1D.x)) 
    481         if data1D.dx is not None and len(data1D.dx) == len(data1D.x): 
    482             self.width = data1D.dx 
    483          
    484         if self.model == None: 
    485             data1d_x = data1D.x 
    486         else: 
    487             self.nbins_low, self.nbins_high, self.width, data1d_x = \ 
    488                                 get_qextrapolate(self.width, data1D.x) 
    489  
    490         ## Number of Q bins 
    491         self.nbins = len(data1d_x) 
    492         ## Minimum Q  
    493         self.min = min(data1d_x) 
    494         ## Maximum 
    495         self.max = max(data1d_x) 
    496         ## Q-values 
    497         self.qvalues = data1d_x 
    498  
    499          
    500 def get_qextrapolate(width, data_x): 
    501     """ 
    502     Make fake data_x points extrapolated outside of the data_x points 
    503      
    504     :param width: array of std of q resolution 
    505     :param Data1D.x: Data1D.x array 
    506      
    507     :return new_width, data_x_ext: extrapolated width array and x array 
    508      
    509     :assumption1: data_x is ordered from lower q to higher q 
    510     :assumption2: len(data) = len(width) 
    511     :assumption3: the distance between the data points is more compact than the size of width  
    512     :Todo1: Make sure that the assumptions are correct for Data1D 
    513     :Todo2: This fixes the edge problem in Qsmearer but still needs to make smearer interface  
    514     """ 
    515     # Length of the width 
    516     length = len(width) 
    517     width_low = math.fabs(width[0])    
    518     width_high = math.fabs(width[length -1]) 
    519     nbins_low = 0.0  
    520     nbins_high = 0.0 
    521     # Compare width(dQ) to the data bin size and take smaller one as the bin  
    522     # size of the extrapolation; this will correct some weird behavior  
    523     # at the edge: This method was out (commented)  
    524     # because it becomes very expansive when 
    525     # bin size is very small comparing to the width. 
    526     # Now on, we will just give the bin size of the extrapolated points  
    527     # based on the width. 
    528     # Find bin sizes 
    529     #bin_size_low = math.fabs(data_x[1] - data_x[0]) 
    530     #bin_size_high = math.fabs(data_x[length - 1] - data_x[length - 2]) 
    531     # Let's set the bin size 1/3 of the width(sigma), it is good as long as 
    532     # the scattering is monotonous. 
    533     #if width_low < (bin_size_low): 
    534     bin_size_low = width_low / 10.0 
    535     #if width_high < (bin_size_high): 
    536     bin_size_high = width_high / 10.0 
    537          
    538     # Number of q points required below the 1st data point in order to extend 
    539     # them 3 times of the width (std) 
    540     if width_low > 0.0: 
    541         nbins_low = math.ceil(3.0 * width_low / bin_size_low) 
    542     # Number of q points required above the last data point 
    543     if width_high > 0.0: 
    544         nbins_high = math.ceil(3.0 * width_high / bin_size_high) 
    545     # Make null q points         
    546     extra_low = numpy.zeros(nbins_low) 
    547     extra_high = numpy.zeros(nbins_high) 
    548     # Give extrapolated values 
    549     ind = 0 
    550     qvalue = data_x[0] - bin_size_low 
    551     #if qvalue > 0: 
    552     while(ind < nbins_low): 
    553         extra_low[nbins_low - (ind + 1)] = qvalue 
    554         qvalue -= bin_size_low 
    555         ind += 1 
    556         #if qvalue <= 0: 
    557         #    break 
    558     # Redefine nbins_low 
    559     nbins_low = ind 
    560     # Reset ind for another extrapolation 
    561     ind = 0 
    562     qvalue = data_x[length -1] + bin_size_high 
    563     while(ind < nbins_high): 
    564         extra_high[ind] = qvalue 
    565         qvalue += bin_size_high 
    566         ind += 1 
    567     # Make a new qx array 
    568     if nbins_low > 0:   
    569         data_x_ext = numpy.append(extra_low, data_x) 
    570     else: 
    571         data_x_ext = data_x 
    572     data_x_ext = numpy.append(data_x_ext, extra_high) 
    573      
    574     # Redefine extra_low and high based on corrected nbins   
    575     # And note that it is not necessary for extra_width to be a non-zero  
    576     if nbins_low > 0:      
    577         extra_low = numpy.zeros(nbins_low) 
    578     extra_high = numpy.zeros(nbins_high)  
    579     # Make new width array 
    580     new_width = numpy.append(extra_low, width) 
    581     new_width = numpy.append(new_width, extra_high) 
    582      
    583     # nbins corrections due to the negative q value 
    584     nbins_low = nbins_low - len(data_x_ext[data_x_ext <= 0]) 
    585     return nbins_low, nbins_high, \ 
    586              new_width[data_x_ext > 0], data_x_ext[data_x_ext > 0] 
    58786 
    58887 
    589  
    590 from .resolution import Slit1D, Pinhole1D 
    59188class PySmear(object): 
    59289    """ 
Note: See TracChangeset for help on using the changeset viewer.