source: sasview/DataLoader/qsmearing.py @ 5385a9e

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 5385a9e was 1b02da1d, checked in by Jae Cho <jhjcho@…>, 14 years ago

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