Changeset 5859862 in sasview


Ignore:
Timestamp:
Sep 20, 2009 3:41:40 PM (15 years ago)
Author:
Mathieu Doucet <doucetm@…>
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:
a65ffcb
Parents:
a1100c07
Message:

dataloader: fixed problem with uneven binning. Added utility method to get first and last bins in unsmeared q.

Location:
DataLoader
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • DataLoader/extensions/smearer.cpp

    ra3f8d58 r5859862  
    3333        // whether we need to compute one 
    3434        has_matrix = false; 
     35        even_binning = true; 
     36}; 
     37 
     38/** 
     39 * Constructor for BaseSmearer 
     40 * 
     41 * Used for uneven binning 
     42 * @param q: array of Q values 
     43 * @param nbins: number of Q bins 
     44 */ 
     45BaseSmearer :: BaseSmearer(double* q, int nbins) { 
     46        // Number of bins 
     47        this->nbins = nbins; 
     48        this->q_values = q; 
     49        // Flag to keep track of whether we have a smearing matrix or 
     50        // whether we need to compute one 
     51        has_matrix = false; 
     52        even_binning = false; 
    3553}; 
    3654 
    … …  
    5068}; 
    5169 
    52         /** 
    53          * Constructor for SlitSmearer 
    54          * 
    55          * @param width: array slit widths for each Q point, in Q units 
    56          * @param qmin: minimum Q value 
    57          * @param qmax: maximum Q value 
    58          * @param nbins: number of Q bins 
    59          */ 
     70/** 
     71 * Constructor for SlitSmearer 
     72 * 
     73 * @param width: slit width in Q units 
     74 * @param height: slit height in Q units 
     75 * @param q: array of Q values 
     76 * @param nbins: number of Q bins 
     77 */ 
     78SlitSmearer :: SlitSmearer(double width, double height, double* q, int nbins) : 
     79        BaseSmearer(q, nbins){ 
     80        this->height = height; 
     81        this->width = width; 
     82}; 
     83 
     84/** 
     85 * Constructor for QSmearer 
     86 * 
     87 * @param width: array slit widths for each Q point, in Q units 
     88 * @param qmin: minimum Q value 
     89 * @param qmax: maximum Q value 
     90 * @param nbins: number of Q bins 
     91 */ 
    6092QSmearer :: QSmearer(double* width, double qmin, double qmax, int nbins) : 
    6193        BaseSmearer(qmin, qmax, nbins){ 
    … …  
    6496 
    6597/** 
     98 * Constructor for QSmearer 
     99 * 
     100 * @param width: array slit widths for each Q point, in Q units 
     101 * @param q: array of Q values 
     102 * @param nbins: number of Q bins 
     103 */ 
     104QSmearer :: QSmearer(double* width, double* q, int nbins) : 
     105        BaseSmearer(q, nbins){ 
     106        this->width = width; 
     107}; 
     108 
     109/** 
    66110 * Compute the slit smearing matrix 
     111 * 
     112 * For even binning (q_min to q_max with nbins): 
     113 * 
     114 *   step = (q_max-q_min)/(nbins-1) 
     115 *   first bin goes from q_min to q_min+step 
     116 *   last bin goes from q_max to q_max+step 
     117 * 
     118 * For binning according to q array: 
     119 * 
     120 * Each q point represents a bin going from half the distance between it 
     121 * and the previous point to half the distance between it and the next point. 
     122 * 
     123 *    Example: bin i goes from (q_values[i-1]+q_values[i])/2 to (q_values[i]+q_values[i+1])/2 
     124 * 
     125 * The exceptions are the first and last bins, which are centered at the first and 
     126 * last q-values, respectively. The width of the first and last bins is the distance between 
     127 * their respective neighboring q-value. 
    67128 */ 
    68129void SlitSmearer :: compute_matrix(){ 
    69130 
    70131        weights = new vector<double>(nbins*nbins,0); 
     132 
     133        // Check the length of the data 
     134        if (nbins<2) return; 
    71135 
    72136        // Loop over all q-values 
    73137        for(int i=0; i<nbins; i++) { 
    74                 double q = qmin + (double)i*(qmax-qmin)/((double)nbins-1.0); 
     138                double q, q_min, q_max; 
     139                get_bin_range(i, &q, &q_min, &q_max); 
    75140 
    76141                // For each q-value, compute the weight of each other q-bin 
    … …  
    102167                                } 
    103168                                double q_shifted = sqrt( ((q-shift_w)*(q-shift_w) + shift_h*shift_h) ); 
    104                                 int q_i = (int)(floor( (q_shifted-qmin) /((qmax-qmin)/((double)nbins -1.0)) )); 
     169 
     170                                // Find in which bin this shifted value belongs 
     171                                int q_i=nbins; 
     172                                if (even_binning) { 
     173                                        // This is kept for backward compatibility since the binning 
     174                                        // was originally defined differently for even bins. 
     175                                        q_i = (int)(floor( (q_shifted-qmin) /((qmax-qmin)/((double)nbins -1.0)) )); 
     176                                } else { 
     177                                        for(int t=0; t<nbins; t++) { 
     178                                                double q_t, q_high, q_low; 
     179                                                get_bin_range(t, &q_t, &q_low, &q_high); 
     180                                                if(q_shifted>=q_low && q_shifted<q_high) { 
     181                                                        q_i = t; 
     182                                                        break; 
     183                                                } 
     184                                        } 
     185                                } 
    105186 
    106187                                // Skip the entries outside our I(q) range 
    … …  
    121202        // Loop over all q-values 
    122203        double step = (qmax-qmin)/((double)nbins-1.0); 
     204        double q, q_min, q_max; 
     205        double q_j, q_jmax, q_jmin; 
    123206        for(int i=0; i<nbins; i++) { 
    124                 double q = qmin + (double)i*step; 
    125                 double q_min = q - 0.5*step; 
    126                 double q_max = q + 0.5*step; 
     207                get_bin_range(i, &q, &q_min, &q_max); 
    127208 
    128209                for(int j=0; j<nbins; j++) { 
    129                         double q_j = qmin + (double)j*step; 
     210                        get_bin_range(j, &q_j, &q_jmin, &q_jmax); 
    130211 
    131212                        // Compute the fraction of the Gaussian contributing 
    … …  
    139220 
    140221/** 
     222 * Computes the Q range of a given bin of the Q distribution. 
     223 * The range is computed according the the data distribution that 
     224 * was given to the object at initialization. 
     225 * 
     226 * @param i: number of the bin in the distribution 
     227 * @param q: q-value of bin i 
     228 * @param q_min: lower bound of the bin 
     229 * @param q_max: higher bound of the bin 
     230 * 
     231 */ 
     232void BaseSmearer :: get_bin_range(int i, double* q, double* q_min, double* q_max) { 
     233        if (even_binning) { 
     234                double step = (qmax-qmin)/((double)nbins-1.0); 
     235                *q = qmin + (double)i*step; 
     236                *q_min = *q - 0.5*step; 
     237                *q_max = *q + 0.5*step; 
     238        } else { 
     239                *q = q_values[i]; 
     240                if (i==0) { 
     241                        double step = (q_values[1]-q_values[0])/2.0; 
     242                        *q_min = *q - step; 
     243                        *q_max = *q + step; 
     244                } else if (i==nbins-1) { 
     245                        double step = (q_values[i]-q_values[i-1])/2.0; 
     246                        *q_min = *q - step; 
     247                        *q_max = *q + step; 
     248                } else { 
     249                        *q_min = *q - (q_values[i]-q_values[i-1])/2.0; 
     250                        *q_max = *q + (q_values[i+1]-q_values[i])/2.0; 
     251                } 
     252        } 
     253} 
     254 
     255/** 
    141256 * Perform smearing by applying the smearing matrix to the input Q array 
    142257 */ 

  • DataLoader/extensions/smearer.hh

    ra3f8d58 r5859862  
    2424class BaseSmearer { 
    2525protected: 
     26        // Internal flag: true when the weights vector is filled 
    2627        bool has_matrix; 
    27     // Smearing matrix 
    28     vector<double>* weights; 
     28        // True when only qmin, qmax and nbins are given. Bins are equidistant 
     29        bool even_binning; 
     30        // Smearing matrix 
     31        vector<double>* weights; 
     32        // Q vector 
     33        double* q_values; 
    2934    // Q_min (Min Q-value for I(q)) 
    3035    double qmin; 
    … …  
    3742    // Constructor 
    3843    BaseSmearer(double qmin, double qmax, int nbins); 
     44    BaseSmearer(double* q, int nbins); 
    3945        // Smear function 
    4046        virtual void smear(double *, double *, int, int); 
    … …  
    4349        // Utility function to check the number of bins 
    4450        int get_nbins() { return nbins; } 
     51        // Get the q range of a particular bin 
     52        virtual void get_bin_range(int, double*, double*, double*); 
    4553}; 
    4654 
    … …  
    6371    // Constructor 
    6472        SlitSmearer(double width, double height, double qmin, double qmax, int nbins); 
     73        SlitSmearer(double width, double height, double* q, int nbins); 
    6574        // Compute the smearing matrix 
    6675        virtual void compute_matrix(); 
    … …  
    8089 
    8190    // Constructor 
    82     QSmearer(double* width,double qmin, double qmax, int nbins); 
     91    QSmearer(double* width, double qmin, double qmax, int nbins); 
     92    QSmearer(double* width, double* q, int nbins); 
    8393        // Compute the smearing matrix 
    8494        virtual void compute_matrix(); 

  • DataLoader/extensions/smearer_module.cpp

    ra3f8d58 r5859862  
    5454 
    5555/** 
     56 * Create a SlitSmearer as a python object by supplying a q array 
     57 */ 
     58PyObject * new_slit_smearer_with_q(PyObject *, PyObject *args) { 
     59        PyObject *qvalues_obj; 
     60        Py_ssize_t n_q; 
     61        double* qvalues; 
     62        double width; 
     63        double height; 
     64 
     65        if (!PyArg_ParseTuple(args, "ddO", &width, &height, &qvalues_obj)) return NULL; 
     66        OUTVECTOR(qvalues_obj, qvalues, n_q); 
     67        SlitSmearer* smearer = new SlitSmearer(width, height, qvalues, n_q); 
     68        return PyCObject_FromVoidPtr(smearer, del_slit_smearer); 
     69} 
     70 
     71/** 
    5672 * Delete a QSmearer object 
    5773 */ 
    … …  
    7692        OUTVECTOR(width_obj, width, n_w); 
    7793        QSmearer* smearer = new QSmearer(width, qmin, qmax, nbins); 
     94        return PyCObject_FromVoidPtr(smearer, del_q_smearer); 
     95} 
     96 
     97/** 
     98 * Create a QSmearer as a python object by supplying a q array 
     99 */ 
     100PyObject * new_q_smearer_with_q(PyObject *, PyObject *args) { 
     101        PyObject *width_obj; 
     102        Py_ssize_t n_w; 
     103        double* width; 
     104        PyObject *qvalues_obj; 
     105        Py_ssize_t n_q; 
     106        double* qvalues; 
     107 
     108        if (!PyArg_ParseTuple(args, "OO", &width_obj, &qvalues_obj)) return NULL; 
     109        OUTVECTOR(width_obj, width, n_w); 
     110        OUTVECTOR(qvalues_obj, qvalues, n_q); 
     111        QSmearer* smearer = new QSmearer(width, qvalues, n_q); 
    78112        return PyCObject_FromVoidPtr(smearer, del_q_smearer); 
    79113} 
    … …  
    114148} 
    115149 
     150/** 
     151 * Get the q-value of a given bin 
     152 */ 
     153PyObject * get_q(PyObject *, PyObject *args) { 
     154        PyObject *smear_obj; 
     155        int bin; 
     156 
     157        if (!PyArg_ParseTuple(args, "Oi", &smear_obj, &bin)) return NULL; 
     158 
     159        // Set the array pointers 
     160        void *temp = PyCObject_AsVoidPtr(smear_obj); 
     161        BaseSmearer* s = static_cast<BaseSmearer *>(temp); 
     162 
     163        double q, q_min, q_max; 
     164        s->get_bin_range(bin, &q, &q_min, &q_max); 
     165 
     166        return Py_BuildValue("d", q); 
     167} 
     168 
    116169 
    117170/** 
    … …  
    121174        {"new_slit_smearer", (PyCFunction)new_slit_smearer, METH_VARARGS, 
    122175                  "Create a new Slit Smearer object"}, 
     176        {"new_slit_smearer_with_q", (PyCFunction)new_slit_smearer_with_q, METH_VARARGS, 
     177                  "Create a new Slit Smearer object by supplying an array of Q values"}, 
    123178        {"new_q_smearer", (PyCFunction)new_q_smearer, METH_VARARGS, 
    124179                  "Create a new Q Smearer object"}, 
     180        {"new_q_smearer_with_q", (PyCFunction)new_q_smearer_with_q, METH_VARARGS, 
     181                  "Create a new Q Smearer object by supplying an array of Q values"}, 
    125182        {"smear",(PyCFunction)smear_input, METH_VARARGS, 
    126183                  "Smear the given input array"}, 
     184        {"get_q",(PyCFunction)get_q, METH_VARARGS, 
     185                  "Get the q-value of a given bin"}, 
    127186    {NULL} 
    128187}; 

  • DataLoader/qsmearing.py

    ra3f8d58 r5859862  
    1212import numpy 
    1313import math 
     14import logging, sys 
    1415import scipy.special 
    1516 
    … …  
    99100         
    100101    def _compute_matrix(self): return NotImplemented 
     102 
     103    def get_bin_range(self, q_min=None, q_max=None): 
     104        """ 
     105            @param q_min: minimum q-value to smear 
     106            @param q_max: maximum q-value to smear  
     107        """ 
     108        if q_min == None: 
     109            q_min = self.min 
     110         
     111        if q_max == None: 
     112            q_max = self.max 
     113         
     114        _qmin_unsmeared, _qmax_unsmeared = self.get_unsmeared_range(q_min, q_max) 
     115         
     116        _first_bin = None 
     117        _last_bin  = None 
     118 
     119        step = (self.max-self.min)/(self.nbins-1.0) 
     120        for i in range(self.nbins): 
     121            q_i = smearer.get_q(self._smearer, i) 
     122            if (q_i >= _qmin_unsmeared) and (q_i <= _qmax_unsmeared): 
     123                # Identify first and last bin 
     124                if _first_bin is None: 
     125                    _first_bin = i 
     126                else: 
     127                    _last_bin  = i 
     128                
     129        return _first_bin, _last_bin 
    101130 
    102131    def __call__(self, iq_in, first_bin=0, last_bin=None): 
    … …  
    161190            This method HAS to be called before smearing 
    162191        """ 
    163         self._smearer = smearer.new_slit_smearer(self.width, self.height, self.min, self.max, self.nbins) 
     192        #self._smearer = smearer.new_slit_smearer(self.width, self.height, self.min, self.max, self.nbins) 
     193        self._smearer = smearer.new_slit_smearer_with_q(self.width, self.height, self.qvalues) 
    164194        self._init_complete = True 
    165195 
     196    def get_unsmeared_range(self, q_min, q_max): 
     197        """ 
     198            Determine the range needed in unsmeared-Q to cover 
     199            the smeared Q range 
     200        """ 
     201        # Range used for input to smearing 
     202        _qmin_unsmeared = q_min 
     203        _qmax_unsmeared = q_max  
     204        try: 
     205            _qmin_unsmeared = self.min 
     206            _qmax_unsmeared = self.max 
     207        except: 
     208            logging.error("_SlitSmearer.get_bin_range: %s" % sys.exc_value) 
     209        return _qmin_unsmeared, _qmax_unsmeared 
    166210 
    167211class SlitSmearer(_SlitSmearer): 
    … …  
    199243        self.nbins = len(data1D.x) 
    200244        ## Minimum Q  
    201         self.min = data1D.x[0] 
     245        self.min = min(data1D.x) 
    202246        ## Maximum 
    203         self.max = data1D.x[len(data1D.x)-1]         
    204  
    205         #print "nbin,npts",self.nbins,self.npts 
     247        self.max = max(data1D.x) 
     248        ## Q-values 
     249        self.qvalues = data1D.x 
     250         
    206251 
    207252class _QSmearer(_BaseSmearer): 
    … …  
    236281            This method HAS to be called before smearing 
    237282        """ 
    238         self._smearer = smearer.new_q_smearer(numpy.asarray(self.width), self.min, self.max, self.nbins) 
     283        #self._smearer = smearer.new_q_smearer(numpy.asarray(self.width), self.min, self.max, self.nbins) 
     284        self._smearer = smearer.new_q_smearer_with_q(numpy.asarray(self.width), self.qvalues) 
    239285        self._init_complete = True 
     286         
     287    def get_unsmeared_range(self, q_min, q_max): 
     288        """ 
     289            Determine the range needed in unsmeared-Q to cover 
     290            the smeared Q range 
     291            Take 3 sigmas as the offset between smeared and unsmeared space 
     292        """ 
     293        # Range used for input to smearing 
     294        _qmin_unsmeared = q_min 
     295        _qmax_unsmeared = q_max  
     296        try: 
     297            offset = 3.0*max(self.width) 
     298            _qmin_unsmeared = max([self.min, q_min-offset]) 
     299            _qmax_unsmeared = min([self.max, q_max+offset]) 
     300        except: 
     301            logging.error("_QSmearer.get_bin_range: %s" % sys.exc_value) 
     302        return _qmin_unsmeared, _qmax_unsmeared 
     303         
     304         
    240305         
    241306class QSmearer(_QSmearer): 
    … …  
    260325        self.nbins = len(data1D.x) 
    261326        ## Minimum Q  
    262         self.min = data1D.x[0] 
     327        self.min = min(data1D.x) 
    263328        ## Maximum 
    264         self.max = data1D.x[len(data1D.x)-1]         
     329        self.max = max(data1D.x) 
     330        ## Q-values 
     331        self.qvalues = data1D.x 
    265332 
    266333 

  • DataLoader/test/utest_smearing.py

    ra3f8d58 r5859862  
    5454        input = 12.0-numpy.arange(1,11) 
    5555        output = s(input) 
    56         answer = [ 9.666,  9.056,  8.329,  7.494,  6.642,  5.721,  4.774,  3.824,  2.871, 2.   ] 
     56        # The following commented line was the correct output for even bins [see smearer.cpp for details]  
     57        #answer = [ 9.666,  9.056,  8.329,  7.494,  6.642,  5.721,  4.774,  3.824,  2.871, 2.   ] 
     58        answer = [ 9.2302,  8.6806,  7.9533,  7.1673,  6.2889,  5.4,     4.5028,  3.5744,  2.6083, 2.    ] 
    5759        for i in range(len(input)): 
    5860            self.assertAlmostEqual(answer[i], output[i], 2) 
Note: See TracChangeset for help on using the changeset viewer.