Changeset 5859862 in sasview for DataLoader/extensions

Timestamp:

Sep 20, 2009 5: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/extensions

Files:

• smearer.cpp (modified) (6 diffs)
• smearer.hh (modified) (5 diffs)
• smearer_module.cpp (modified) (4 diffs)

DataLoader/extensions/smearer.cpp

@@ -33 +33 @@
         // whether we need to compute one
         has_matrix = false;
+        even_binning = true;
+};
+
+/**
+ * Constructor for BaseSmearer
+ *
+ * Used for uneven binning
+ * @param q: array of Q values
+ * @param nbins: number of Q bins
+ */
+BaseSmearer :: BaseSmearer(double* q, int nbins) {
+        // Number of bins
+        this->nbins = nbins;
+        this->q_values = q;
+        // Flag to keep track of whether we have a smearing matrix or
+        // whether we need to compute one
+        has_matrix = false;
+        even_binning = false;
 };
 
@@ -50 +68 @@
 };
 
-        /**
-         * Constructor for SlitSmearer
-         *
-         * @param width: array slit widths for each Q point, in Q units
-         * @param qmin: minimum Q value
-         * @param qmax: maximum Q value
-         * @param nbins: number of Q bins
-         */
+/**
+ * Constructor for SlitSmearer
+ *
+ * @param width: slit width in Q units
+ * @param height: slit height in Q units
+ * @param q: array of Q values
+ * @param nbins: number of Q bins
+ */
+SlitSmearer :: SlitSmearer(double width, double height, double* q, int nbins) :
+        BaseSmearer(q, nbins){
+        this->height = height;
+        this->width = width;
+};
+
+/**
+ * Constructor for QSmearer
+ *
+ * @param width: array slit widths for each Q point, in Q units
+ * @param qmin: minimum Q value
+ * @param qmax: maximum Q value
+ * @param nbins: number of Q bins
+ */
 QSmearer :: QSmearer(double* width, double qmin, double qmax, int nbins) :
         BaseSmearer(qmin, qmax, nbins){
@@ -64 +96 @@
 
 /**
+ * Constructor for QSmearer
+ *
+ * @param width: array slit widths for each Q point, in Q units
+ * @param q: array of Q values
+ * @param nbins: number of Q bins
+ */
+QSmearer :: QSmearer(double* width, double* q, int nbins) :
+        BaseSmearer(q, nbins){
+        this->width = width;
+};
+
+/**
  * Compute the slit smearing matrix
+ *
+ * For even binning (q_min to q_max with nbins):
+ *
+ *   step = (q_max-q_min)/(nbins-1)
+ *   first bin goes from q_min to q_min+step
+ *   last bin goes from q_max to q_max+step
+ *
+ * For binning according to q array:
+ *
+ * Each q point represents a bin going from half the distance between it
+ * and the previous point to half the distance between it and the next point.
+ *
+ *    Example: bin i goes from (q_values[i-1]+q_values[i])/2 to (q_values[i]+q_values[i+1])/2
+ *
+ * The exceptions are the first and last bins, which are centered at the first and
+ * last q-values, respectively. The width of the first and last bins is the distance between
+ * their respective neighboring q-value.
  */
 void SlitSmearer :: compute_matrix(){
 
         weights = new vector<double>(nbins*nbins,0);
+
+        // Check the length of the data
+        if (nbins<2) return;
 
         // Loop over all q-values
         for(int i=0; i<nbins; i++) {
-                double q = qmin + (double)i*(qmax-qmin)/((double)nbins-1.0);
+                double q, q_min, q_max;
+                get_bin_range(i, &q, &q_min, &q_max);
 
                 // For each q-value, compute the weight of each other q-bin
@@ -102 +167 @@
                                 }
                                 double q_shifted = sqrt( ((q-shift_w)*(q-shift_w) + shift_h*shift_h) );
-                                int q_i = (int)(floor( (q_shifted-qmin) /((qmax-qmin)/((double)nbins -1.0)) ));
+
+                                // Find in which bin this shifted value belongs
+                                int q_i=nbins;
+                                if (even_binning) {
+                                        // This is kept for backward compatibility since the binning
+                                        // was originally defined differently for even bins.
+                                        q_i = (int)(floor( (q_shifted-qmin) /((qmax-qmin)/((double)nbins -1.0)) ));
+                                } else {
+                                        for(int t=0; t<nbins; t++) {
+                                                double q_t, q_high, q_low;
+                                                get_bin_range(t, &q_t, &q_low, &q_high);
+                                                if(q_shifted>=q_low && q_shifted<q_high) {
+                                                        q_i = t;
+                                                        break;
+                                                }
+                                        }
+                                }
 
                                 // Skip the entries outside our I(q) range
@@ -121 +202 @@
         // Loop over all q-values
         double step = (qmax-qmin)/((double)nbins-1.0);
+        double q, q_min, q_max;
+        double q_j, q_jmax, q_jmin;
         for(int i=0; i<nbins; i++) {
-                double q = qmin + (double)i*step;
-                double q_min = q - 0.5*step;
-                double q_max = q + 0.5*step;
+                get_bin_range(i, &q, &q_min, &q_max);
 
                 for(int j=0; j<nbins; j++) {
-                        double q_j = qmin + (double)j*step;
+                        get_bin_range(j, &q_j, &q_jmin, &q_jmax);
 
                         // Compute the fraction of the Gaussian contributing
@@ -139 +220 @@
 
 /**
+ * Computes the Q range of a given bin of the Q distribution.
+ * The range is computed according the the data distribution that
+ * was given to the object at initialization.
+ *
+ * @param i: number of the bin in the distribution
+ * @param q: q-value of bin i
+ * @param q_min: lower bound of the bin
+ * @param q_max: higher bound of the bin
+ *
+ */
+void BaseSmearer :: get_bin_range(int i, double* q, double* q_min, double* q_max) {
+        if (even_binning) {
+                double step = (qmax-qmin)/((double)nbins-1.0);
+                *q = qmin + (double)i*step;
+                *q_min = *q - 0.5*step;
+                *q_max = *q + 0.5*step;
+        } else {
+                *q = q_values[i];
+                if (i==0) {
+                        double step = (q_values[1]-q_values[0])/2.0;
+                        *q_min = *q - step;
+                        *q_max = *q + step;
+                } else if (i==nbins-1) {
+                        double step = (q_values[i]-q_values[i-1])/2.0;
+                        *q_min = *q - step;
+                        *q_max = *q + step;
+                } else {
+                        *q_min = *q - (q_values[i]-q_values[i-1])/2.0;
+                        *q_max = *q + (q_values[i+1]-q_values[i])/2.0;
+                }
+        }
+}
+
+/**
  * Perform smearing by applying the smearing matrix to the input Q array
  */

DataLoader/extensions/smearer.hh

@@ -24 +24 @@
 class BaseSmearer {
 protected:
+        // Internal flag: true when the weights vector is filled
         bool has_matrix;
-    // Smearing matrix
-    vector<double>* weights;
+        // True when only qmin, qmax and nbins are given. Bins are equidistant
+        bool even_binning;
+        // Smearing matrix
+        vector<double>* weights;
+        // Q vector
+        double* q_values;
     // Q_min (Min Q-value for I(q))
     double qmin;
@@ -37 +42 @@
     // Constructor
     BaseSmearer(double qmin, double qmax, int nbins);
+    BaseSmearer(double* q, int nbins);
         // Smear function
         virtual void smear(double *, double *, int, int);
@@ -43 +49 @@
         // Utility function to check the number of bins
         int get_nbins() { return nbins; }
+        // Get the q range of a particular bin
+        virtual void get_bin_range(int, double*, double*, double*);
 };
 
@@ -63 +71 @@
     // Constructor
         SlitSmearer(double width, double height, double qmin, double qmax, int nbins);
+        SlitSmearer(double width, double height, double* q, int nbins);
         // Compute the smearing matrix
         virtual void compute_matrix();
@@ -80 +89 @@
 
     // Constructor
-    QSmearer(double* width,double qmin, double qmax, int nbins);
+    QSmearer(double* width, double qmin, double qmax, int nbins);
+    QSmearer(double* width, double* q, int nbins);
         // Compute the smearing matrix
         virtual void compute_matrix();

DataLoader/extensions/smearer_module.cpp

@@ -54 +54 @@
 
 /**
+ * Create a SlitSmearer as a python object by supplying a q array
+ */
+PyObject * new_slit_smearer_with_q(PyObject *, PyObject *args) {
+        PyObject *qvalues_obj;
+        Py_ssize_t n_q;
+        double* qvalues;
+        double width;
+        double height;
+
+        if (!PyArg_ParseTuple(args, "ddO", &width, &height, &qvalues_obj)) return NULL;
+        OUTVECTOR(qvalues_obj, qvalues, n_q);
+        SlitSmearer* smearer = new SlitSmearer(width, height, qvalues, n_q);
+        return PyCObject_FromVoidPtr(smearer, del_slit_smearer);
+}
+
+/**
  * Delete a QSmearer object
  */
@@ -76 +92 @@
         OUTVECTOR(width_obj, width, n_w);
         QSmearer* smearer = new QSmearer(width, qmin, qmax, nbins);
+        return PyCObject_FromVoidPtr(smearer, del_q_smearer);
+}
+
+/**
+ * Create a QSmearer as a python object by supplying a q array
+ */
+PyObject * new_q_smearer_with_q(PyObject *, PyObject *args) {
+        PyObject *width_obj;
+        Py_ssize_t n_w;
+        double* width;
+        PyObject *qvalues_obj;
+        Py_ssize_t n_q;
+        double* qvalues;
+
+        if (!PyArg_ParseTuple(args, "OO", &width_obj, &qvalues_obj)) return NULL;
+        OUTVECTOR(width_obj, width, n_w);
+        OUTVECTOR(qvalues_obj, qvalues, n_q);
+        QSmearer* smearer = new QSmearer(width, qvalues, n_q);
         return PyCObject_FromVoidPtr(smearer, del_q_smearer);
 }
@@ -114 +148 @@
 }
 
+/**
+ * Get the q-value of a given bin
+ */
+PyObject * get_q(PyObject *, PyObject *args) {
+        PyObject *smear_obj;
+        int bin;
+
+        if (!PyArg_ParseTuple(args, "Oi", &smear_obj, &bin)) return NULL;
+
+        // Set the array pointers
+        void *temp = PyCObject_AsVoidPtr(smear_obj);
+        BaseSmearer* s = static_cast<BaseSmearer *>(temp);
+
+        double q, q_min, q_max;
+        s->get_bin_range(bin, &q, &q_min, &q_max);
+
+        return Py_BuildValue("d", q);
+}
+
 
 /**
@@ -121 +174 @@
         {"new_slit_smearer", (PyCFunction)new_slit_smearer, METH_VARARGS,
                   "Create a new Slit Smearer object"},
+        {"new_slit_smearer_with_q", (PyCFunction)new_slit_smearer_with_q, METH_VARARGS,
+                  "Create a new Slit Smearer object by supplying an array of Q values"},
         {"new_q_smearer", (PyCFunction)new_q_smearer, METH_VARARGS,
                   "Create a new Q Smearer object"},
+        {"new_q_smearer_with_q", (PyCFunction)new_q_smearer_with_q, METH_VARARGS,
+                  "Create a new Q Smearer object by supplying an array of Q values"},
         {"smear",(PyCFunction)smear_input, METH_VARARGS,
                   "Smear the given input array"},
+        {"get_q",(PyCFunction)get_q, METH_VARARGS,
+                  "Get the q-value of a given bin"},
     {NULL}
 };