source: sasview/sansmodels/src/sans/models/qsmearing.py @ f40f743

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since f40f743 was ef8d42f, checked in by Mathieu Doucet <doucetm@…>, 13 years ago

Fixing code style problems and bugs

  • Property mode set to 100644
File size: 20.7 KB
Line 
1"""
2    Handle Q smearing
3"""
4#####################################################################
5#This software was developed by the University of Tennessee as part of the
6#Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
7#project funded by the US National Science Foundation.
8#See the license text in license.txt
9#copyright 2008, University of Tennessee
10######################################################################
11import numpy
12import math
13import logging
14import sys
15import sans.models.sans_extension.smearer as smearer 
16from sans.models.smearing_2d import Smearer2D
17
18def smear_selection(data1D, model = None):
19    """
20    Creates the right type of smearer according
21    to the data.
22
23    The canSAS format has a rule that either
24    slit smearing data OR resolution smearing data
25    is available.
26   
27    For the present purpose, we choose the one that
28    has none-zero data. If both slit and resolution
29    smearing arrays are filled with good data
30    (which should not happen), then we choose the
31    resolution smearing data.
32   
33    :param data1D: Data1D object
34    :param model: sans.model instance
35    """
36    # Sanity check. If we are not dealing with a SANS Data1D
37    # object, just return None
38    if  data1D.__class__.__name__ not in ['Data1D', 'Theory1D']:
39        if data1D == None:
40            return None
41        elif data1D.dqx_data == None or data1D.dqy_data == None:
42            return None
43        return Smearer2D(data1D)
44   
45    if  not hasattr(data1D, "dx") and not hasattr(data1D, "dxl")\
46         and not hasattr(data1D, "dxw"):
47        return None
48   
49    # Look for resolution smearing data
50    _found_resolution = False
51    if data1D.dx is not None and len(data1D.dx) == len(data1D.x):
52       
53        # Check that we have non-zero data
54        if data1D.dx[0] > 0.0:
55            _found_resolution = True
56            #print "_found_resolution",_found_resolution
57            #print "data1D.dx[0]",data1D.dx[0],data1D.dxl[0]
58    # If we found resolution smearing data, return a QSmearer
59    if _found_resolution == True:
60        return QSmearer(data1D, model)
61
62    # Look for slit smearing data
63    _found_slit = False
64    if data1D.dxl is not None and len(data1D.dxl) == len(data1D.x) \
65        and data1D.dxw is not None and len(data1D.dxw) == len(data1D.x):
66       
67        # Check that we have non-zero data
68        if data1D.dxl[0] > 0.0 or data1D.dxw[0] > 0.0:
69            _found_slit = True
70       
71        # Sanity check: all data should be the same as a function of Q
72        for item in data1D.dxl:
73            if data1D.dxl[0] != item:
74                _found_resolution = False
75                break
76           
77        for item in data1D.dxw:
78            if data1D.dxw[0] != item:
79                _found_resolution = False
80                break
81    # If we found slit smearing data, return a slit smearer
82    if _found_slit == True:
83        return SlitSmearer(data1D, model)
84    return None
85           
86
87class _BaseSmearer(object):
88    """
89        Base class for smearers
90    """
91    def __init__(self):
92        self.nbins = 0
93        self.nbins_low = 0
94        self.nbins_high = 0
95        self._weights = None
96        ## Internal flag to keep track of C++ smearer initialization
97        self._init_complete = False
98        self._smearer = None
99        self.model = None
100        self.min = None
101        self.max = None
102        self.qvalues = []
103       
104    def __deepcopy__(self, memo=None):
105        """
106        Return a valid copy of self.
107        Avoid copying the _smearer C object and force a matrix recompute
108        when the copy is used. 
109        """
110        result = _BaseSmearer()
111        result.nbins = self.nbins
112        return result
113
114    def _compute_matrix(self):
115        """
116            Place holder for matrix computation
117        """
118        return NotImplemented
119
120    def get_unsmeared_range(self, q_min=None, q_max=None):
121        """
122            Place holder for method returning unsmeared range
123        """
124        return NotImplemented
125   
126    def get_bin_range(self, q_min=None, q_max=None):
127        """
128       
129        :param q_min: minimum q-value to smear
130        :param q_max: maximum q-value to smear
131         
132        """
133        # If this is the first time we call for smearing,
134        # initialize the C++ smearer object first
135        if not self._init_complete:
136            self._initialize_smearer()
137        if q_min == None:
138            q_min = self.min
139        if q_max == None:
140            q_max = self.max
141
142        _qmin_unsmeared, _qmax_unsmeared = self.get_unsmeared_range(q_min,
143                                                                     q_max)
144        _first_bin = None
145        _last_bin  = None
146
147        #step = (self.max - self.min) / (self.nbins - 1.0)
148        # Find the first and last bin number in all extrapolated and real data
149        try:
150            for i in range(self.nbins):
151                q_i = smearer.get_q(self._smearer, i)
152                if (q_i >= _qmin_unsmeared) and (q_i <= _qmax_unsmeared):
153                    # Identify first and last bin
154                    if _first_bin is None:
155                        _first_bin = i
156                    else:
157                        _last_bin  = i
158        except:
159            msg = "_BaseSmearer.get_bin_range: "
160            msg += " error getting range\n  %s" % sys.exc_value
161            raise RuntimeError, msg
162   
163        #  Find the first and last bin number only in the real data
164        _first_bin, _last_bin = self._get_unextrapolated_bin( \
165                                        _first_bin, _last_bin)
166
167        return _first_bin, _last_bin
168
169    def __call__(self, iq_in, first_bin = 0, last_bin = None):
170        """
171        Perform smearing
172        """
173        # If this is the first time we call for smearing,
174        # initialize the C++ smearer object first
175        if not self._init_complete:
176            self._initialize_smearer()
177
178        if last_bin is None or last_bin >= len(iq_in):
179            last_bin = len(iq_in) - 1
180        # Check that the first bin is positive
181        if first_bin < 0:
182            first_bin = 0
183
184        # With a model given, compute I for the extrapolated points and append
185        # to the iq_in
186        iq_in_temp = iq_in
187        if self.model != None:
188            temp_first, temp_last = self._get_extrapolated_bin( \
189                                                        first_bin, last_bin)
190            if self.nbins_low > 0:
191                iq_in_low = self.model.evalDistribution( \
192                                    numpy.fabs(self.qvalues[0:self.nbins_low]))
193            iq_in_high = self.model.evalDistribution( \
194                                            self.qvalues[(len(self.qvalues) - \
195                                            self.nbins_high - 1):])
196            # Todo: find out who is sending iq[last_poin] = 0.
197            if iq_in[len(iq_in) - 1] == 0:
198                iq_in[len(iq_in) - 1] = iq_in_high[0]
199            # Append the extrapolated points to the data points
200            if self.nbins_low > 0:                             
201                iq_in_temp = numpy.append(iq_in_low, iq_in)
202            if self.nbins_high > 0:
203                iq_in_temp = numpy.append(iq_in_temp, iq_in_high[1:])
204        else:
205            temp_first = first_bin
206            temp_last = last_bin
207            #iq_in_temp = iq_in
208
209        # Sanity check
210        if len(iq_in_temp) != self.nbins:
211            msg = "Invalid I(q) vector: inconsistent array "
212            msg += " length %d != %s" % (len(iq_in_temp), str(self.nbins))
213            raise RuntimeError, msg
214
215        # Storage for smeared I(q)   
216        iq_out = numpy.zeros(self.nbins)
217
218        smear_output = smearer.smear(self._smearer, iq_in_temp, iq_out,
219                                      #0, self.nbins - 1)
220                                      temp_first, temp_last)
221                                      #first_bin, last_bin)
222        if smear_output < 0:
223            msg = "_BaseSmearer: could not smear, code = %g" % smear_output
224            raise RuntimeError, msg
225
226        temp_first = first_bin + self.nbins_low
227        temp_last = self.nbins - self.nbins_high
228        out = iq_out[temp_first: temp_last]
229
230        return out
231   
232    def _initialize_smearer(self):
233        """
234            Place holder for initializing data smearer
235        """
236        return NotImplemented
237           
238   
239    def _get_unextrapolated_bin(self, first_bin = 0, last_bin = 0):
240        """
241        Get unextrapolated first bin and the last bin
242       
243        : param first_bin: extrapolated first_bin
244        : param last_bin: extrapolated last_bin
245       
246        : return fist_bin, last_bin: unextrapolated first and last bin
247        """
248        #  For first bin
249        if first_bin <= self.nbins_low:
250            first_bin = 0
251        else:
252            first_bin = first_bin - self.nbins_low
253        # For last bin
254        if last_bin >= (self.nbins - self.nbins_high):
255            last_bin  = self.nbins - (self.nbins_high + self.nbins_low + 1)
256        elif last_bin >= self.nbins_low:
257            last_bin = last_bin - self.nbins_low
258        else:
259            last_bin = 0
260        return first_bin, last_bin
261   
262    def _get_extrapolated_bin(self, first_bin = 0, last_bin = 0):
263        """
264        Get extrapolated first bin and the last bin
265       
266        : param first_bin: unextrapolated first_bin
267        : param last_bin: unextrapolated last_bin
268       
269        : return first_bin, last_bin: extrapolated first and last bin
270        """
271        #  For the first bin
272        # In the case that needs low extrapolation data
273        first_bin = 0
274        # For last bin
275        if last_bin >= self.nbins - (self.nbins_high + self.nbins_low + 1):
276            # In the case that needs higher q extrapolation data
277            last_bin = self.nbins - 1
278        else:
279            # In the case that doesn't need higher q extrapolation data
280            last_bin += self.nbins_low
281
282        return first_bin, last_bin
283       
284class _SlitSmearer(_BaseSmearer):
285    """
286    Slit smearing for I(q) array
287    """
288   
289    def __init__(self, nbins=None, width=None, height=None, min=None, max=None):
290        """
291        Initialization
292           
293        :param iq: I(q) array [cm-1]
294        :param width: slit width [A-1]
295        :param height: slit height [A-1]
296        :param min: Q_min [A-1]
297        :param max: Q_max [A-1]
298       
299        """
300        _BaseSmearer.__init__(self)
301        ## Slit width in Q units
302        self.width  = width
303        ## Slit height in Q units
304        self.height = height
305        ## Q_min (Min Q-value for I(q))
306        self.min    = min
307        ## Q_max (Max Q_value for I(q))
308        self.max    = max
309        ## Number of Q bins
310        self.nbins  = nbins
311        ## Number of points used in the smearing computation
312        self.npts   = 3000
313        ## Smearing matrix
314        self._weights = None
315        self.qvalues  = None
316       
317    def _initialize_smearer(self):
318        """
319        Initialize the C++ smearer object.
320        This method HAS to be called before smearing
321        """
322        #self._smearer = smearer.new_slit_smearer(self.width,
323        # self.height, self.min, self.max, self.nbins)
324        self._smearer = smearer.new_slit_smearer_with_q(self.width, 
325                                                    self.height, self.qvalues)
326        self._init_complete = True
327
328    def get_unsmeared_range(self, q_min, q_max):
329        """
330        Determine the range needed in unsmeared-Q to cover
331        the smeared Q range
332        """
333        # Range used for input to smearing
334        _qmin_unsmeared = q_min
335        _qmax_unsmeared = q_max
336        try:
337            _qmin_unsmeared = self.min
338            _qmax_unsmeared = self.max
339        except:
340            logging.error("_SlitSmearer.get_bin_range: %s" % sys.exc_value)
341        return _qmin_unsmeared, _qmax_unsmeared
342
343class SlitSmearer(_SlitSmearer):
344    """
345    Adaptor for slit smearing class and SANS data
346    """
347    def __init__(self, data1D, model = None):
348        """
349        Assumption: equally spaced bins of increasing q-values.
350       
351        :param data1D: data used to set the smearing parameters
352        """
353        # Initialization from parent class
354        super(SlitSmearer, self).__init__()
355       
356        ## Slit width
357        self.width = 0
358        self.nbins_low = 0
359        self.nbins_high = 0
360        self.model = model
361        if data1D.dxw is not None and len(data1D.dxw) == len(data1D.x):
362            self.width = data1D.dxw[0]
363            # Sanity check
364            for value in data1D.dxw:
365                if value != self.width:
366                    msg = "Slit smearing parameters must "
367                    msg += " be the same for all data"
368                    raise RuntimeError, msg
369        ## Slit height
370        self.height = 0
371        if data1D.dxl is not None and len(data1D.dxl) == len(data1D.x):
372            self.height = data1D.dxl[0]
373            # Sanity check
374            for value in data1D.dxl:
375                if value != self.height:
376                    msg = "Slit smearing parameters must be"
377                    msg += " the same for all data"
378                    raise RuntimeError, msg
379        # If a model is given, get the q extrapolation
380        if self.model == None:
381            data1d_x = data1D.x
382        else:
383            # Take larger sigma
384            if self.height > self.width:
385                # The denominator (2.0) covers all the possible w^2 + h^2 range
386                sigma_in = data1D.dxl / 2.0
387            elif self.width > 0:
388                sigma_in = data1D.dxw / 2.0
389            else:
390                sigma_in = []
391
392            self.nbins_low, self.nbins_high, _, data1d_x = \
393                                get_qextrapolate(sigma_in, data1D.x)
394
395        ## Number of Q bins
396        self.nbins = len(data1d_x)
397        ## Minimum Q
398        self.min = min(data1d_x)
399        ## Maximum
400        self.max = max(data1d_x)
401        ## Q-values
402        self.qvalues = data1d_x
403       
404
405class _QSmearer(_BaseSmearer):
406    """
407    Perform Gaussian Q smearing
408    """
409       
410    def __init__(self, nbins=None, width=None, min=None, max=None):
411        """
412        Initialization
413       
414        :param nbins: number of Q bins
415        :param width: array standard deviation in Q [A-1]
416        :param min: Q_min [A-1]
417        :param max: Q_max [A-1]
418        """
419        _BaseSmearer.__init__(self)
420        ## Standard deviation in Q [A-1]
421        self.width = width
422        ## Q_min (Min Q-value for I(q))
423        self.min = min
424        ## Q_max (Max Q_value for I(q))
425        self.max = max
426        ## Number of Q bins
427        self.nbins = nbins
428        ## Smearing matrix
429        self._weights = None
430        self.qvalues  = None
431       
432    def _initialize_smearer(self):
433        """
434        Initialize the C++ smearer object.
435        This method HAS to be called before smearing
436        """
437        #self._smearer = smearer.new_q_smearer(numpy.asarray(self.width),
438        # self.min, self.max, self.nbins)
439        self._smearer = smearer.new_q_smearer_with_q(numpy.asarray(self.width),
440                                                      self.qvalues)
441        self._init_complete = True
442       
443    def get_unsmeared_range(self, q_min, q_max):
444        """
445        Determine the range needed in unsmeared-Q to cover
446        the smeared Q range
447        Take 3 sigmas as the offset between smeared and unsmeared space
448        """
449        # Range used for input to smearing
450        _qmin_unsmeared = q_min
451        _qmax_unsmeared = q_max
452        try:
453            offset = 3.0 * max(self.width)
454            _qmin_unsmeared = self.min#max([self.min, q_min - offset])
455            _qmax_unsmeared = self.max#min([self.max, q_max + offset])
456        except:
457            logging.error("_QSmearer.get_bin_range: %s" % sys.exc_value)
458        return _qmin_unsmeared, _qmax_unsmeared
459       
460   
461class QSmearer(_QSmearer):
462    """
463    Adaptor for Gaussian Q smearing class and SANS data
464    """
465    def __init__(self, data1D, model = None):
466        """
467        Assumption: equally spaced bins of increasing q-values.
468       
469        :param data1D: data used to set the smearing parameters
470        """
471        # Initialization from parent class
472        super(QSmearer, self).__init__()
473        data1d_x = []
474        self.nbins_low = 0
475        self.nbins_high = 0
476        self.model = model
477        ## Resolution
478        #self.width = numpy.zeros(len(data1D.x))
479        if data1D.dx is not None and len(data1D.dx) == len(data1D.x):
480            self.width = data1D.dx
481       
482        if self.model == None:
483            data1d_x = data1D.x
484        else:
485            self.nbins_low, self.nbins_high, self.width, data1d_x = \
486                                get_qextrapolate(self.width, data1D.x)
487
488        ## Number of Q bins
489        self.nbins = len(data1d_x)
490        ## Minimum Q
491        self.min = min(data1d_x)
492        ## Maximum
493        self.max = max(data1d_x)
494        ## Q-values
495        self.qvalues = data1d_x
496
497       
498def get_qextrapolate(width, data_x):
499    """
500    Make fake data_x points extrapolated outside of the data_x points
501   
502    : param width: array of std of q resolution
503    : param Data1D.x: Data1D.x array
504   
505    : return new_width, data_x_ext: extrapolated width array and x array
506   
507    : assumption1: data_x is ordered from lower q to higher q
508    : assumption2: len(data) = len(width)
509    : assumption3: the distance between the data points is more compact
510            than the size of width
511    : Todo1: Make sure that the assumptions are correct for Data1D
512    : Todo2: This fixes the edge problem in Qsmearer but still needs to make
513            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]
Note: See TracBrowser for help on using the repository browser.