Changeset cd2ced80 in sasview for DataLoader

Timestamp:

Jan 14, 2011 5:45:58 PM (13 years ago)

Author:

Jae Cho <jhjcho@…>

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:

8fb8b0c

Parents:

67c7e89

Message:

new algorithm for slit smearing and its test

Location:

DataLoader

Files:

• extensions/smearer.cpp (modified) (3 diffs)
• extensions/smearer.hh (modified) (1 diff)
• qsmearing.py (modified) (4 diffs)
• test/slit_1000A_sphere_sm_h.txt (added)
• test/slit_1000A_sphere_sm_w_0_0001.txt (added)
• test/slit_1000A_sphere_sm_w_0_0002.txt (added)
• test/utest_smearing.py (modified) (3 diffs)

DataLoader/extensions/smearer.cpp

@@ -133 +133 @@
         // Check the length of the data
         if (nbins<2) return;
+        int npts_h = height>0.0 ? npts : 1;
+        int npts_w = width>0.0 ? npts : 1;
+
+        // If both height and width are great than zero,
+        // modify the number of points in each direction so
+        // that the total number of points is still what
+        // the user would expect (downgrade resolution)
+        //if(npts_h>1 && npts_w>1){
+        //      npts_h = (int)ceil(sqrt((double)npts));
+        //      npts_w = npts_h;
+        //}
+        double shift_h, shift_w, hbin_size, wbin_size;
+        // Make sure height and width are all positive (FWMH/2)
+        // Assumption; height and width are all same for all q points
+        if(npts_h == 1){
+                shift_h = 0.0;
+        } else {
+                shift_h = fabs(height);
+        }
+        if(npts_w == 1){
+                shift_w = 0.0;
+        } else {
+                shift_w = fabs(width);
+        }
+        // size of the h bin and w bin
+        hbin_size = shift_h / nbins;
+        wbin_size = shift_w / nbins;
 
         // Loop over all q-values
         for(int i=0; i<nbins; i++) {
-                double q, q_min, q_max;
+                // Find Weights
+                // Find q where the resolution smearing calculation of I(q) occurs
+                double q, q_min, q_max, q_0;
                 get_bin_range(i, &q, &q_min, &q_max);
-
-                // For each q-value, compute the weight of each other q-bin
-                // in the I(q) array
-                int npts_h = height>0 ? npts : 1;
-                int npts_w = width>0 ? npts : 1;
-
-                // If both height and width are great than zero,
-                // modify the number of points in each direction so
-                // that the total number of points is still what
-                // the user would expect (downgrade resolution)
-                // Never down-grade npts_h. That will give incorrect slit smearing...
-                if(npts_h>1 && npts_w>1){
-                        npts_h = npts;//(int)ceil(sqrt((double)npts));
-                        // In general width is much smaller than height, so smaller npt_w
-                        // should work fine.
-                        // Todo: It is still very expansive in time. Think about better way.
-                        npts_w = (int)ceil(npts_h / 100);
-                }
-
-                double shift_h, shift_w;
-                for(int k=0; k<npts_h; k++){
-                        if(npts_h==1){
-                                shift_h = 0;
-                        } else {
-                                shift_h = height/((double)npts_h-1.0) * (double)k;
-                        }
-                        for(int j=0; j<npts_w; j++){
-                                if(npts_w==1){
-                                        shift_w = 0;
-                                } else {
-                                        shift_w = width/((double)npts_w-1.0) * (double)j;
+                bool last_qpoint = true;
+                // Find q[0] value to normalize the weight later,
+                //  otherwise, we will have a precision problem.
+                if (i == 0){
+                        q_0 = q;
+                }
+                // Loop over all qj-values
+                bool first_w = true;
+                for(int j=0; j<nbins; j++) {
+                        double q_j, q_high, q_low;
+                        // Calculate bin size of q_j
+                        get_bin_range(j, &q_j, &q_low, &q_high);
+                        // Check q_low that can not be negative.
+                        if (q_low < 0.0){
+                                q_low = 0.0;
+                        }
+                        // default parameter values
+                        (*weights)[i*nbins+j] = 0.0;
+                        double shift_w = 0.0;
+                        // Condition: zero slit smear.
+                        if (npts_w == 1 && npts_h == 1){
+                                if(q_j == q) {
+                                        (*weights)[i*nbins+j] = 1.0;
                                 }
-                                double q_shifted = sqrt( ((q-shift_w)*(q-shift_w) + shift_h*shift_h) );
-
-                                // 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;
+                        }
+                        //Condition:Smear weight integration for width >0 when the height (=0) does not present.
+                        //Or height << width.
+                        else if((npts_w!=1 && npts_h == 1)|| (npts_w!=1 && npts_h != 1 && width/height > 100.0)){
+                                shift_w = width;
+                                //del_w = width/((double)npts_w-1.0);
+                                double q_shifted_low = q - shift_w;
+                                // High limit of the resolution range
+                                double q_shifted_high = q + shift_w;
+                                // Go through all the q_js for weighting those points
+                                if(q_j >= q_shifted_low && q_j <= q_shifted_high) {
+                                        // The weighting factor comes,
+                                        // Give some weight (delq_bin) for the q_j within the resolution range
+                                        // Weight should be same for all qs except
+                                        // for the q bin size at j.
+                                        // Note that the division by q_0 is only due to the precision problem
+                                        // where q_high - q_low gets to very small.
+                                        // Later, it will be normalized again.
+                                        (*weights)[i*nbins+j] += (q_high - q_low)/q_0 ;
+                                }
+                        }
+                        else{
+                                // Loop for width (;Height is analytical.)
+                                // Condition: height >>> width, otherwise, below is not accurate enough.
+                                // Smear weight numerical iteration for width >0 when the height (>0) presents.
+                                // When width = 0, the numerical iteration will be skipped.
+                                // The resolution calculation for the height is done by direct integration,
+                                // assuming the I(q'=sqrt(q_j^2-(q+shift_w)^2)) is constant within a q' bin, [q_high, q_low].
+                                // In general, this weight numerical iteration for width >0 might be a rough approximation,
+                                // but it must be good enough when height >>> width.
+                                for(int k=(-npts_w + 1); k<npts_w; k++){
+                                        if(npts_w!=1){
+                                                shift_w = width/((double)npts_w-1.0)*(double)k;
+                                        }
+                                        // For each q-value, compute the weight of each other q-bin
+                                        // in the I(q) array
+                                        // Low limit of the resolution range
+                                        double q_shift = q + shift_w;
+                                        if (q_shift < 0.0){
+                                                q_shift = 0.0;
+                                        }
+                                        double q_shifted_low = q_shift;
+                                        // High limit of the resolution range
+                                        double q_shifted_high = sqrt(q_shift * q_shift + shift_h * shift_h);
+
+
+                                        // Go through all the q_js for weighting those points
+                                        if(q_j >= q_shifted_low && q_j <= q_shifted_high) {
+                                                // The weighting factor comes,
+                                                // Give some weight (delq_bin) for the q_j within the resolution range
+                                                // Weight should be same for all qs except
+                                                // for the q bin size at j.
+                                                // Note that the division by q_0 is only due to the precision problem
+                                                // where q_high - q_low gets to very small.
+                                                // Later, it will be normalized again.
+
+                                                double q_shift_min = q - width;
+
+                                                double u = (q_j * q_j - (q_shift) * (q_shift));
+                                                // The fabs below are not necessary but in case: the weight should never be imaginary.
+                                                // At the edge of each sub_width. weight += u(at q_high bin) - u(0), where u(0) = 0,
+                                                // and weighted by (2.0* npts_w -1.0)once for each q.
+                                                if (q == q_j) {
+                                                        if (k==0)
+                                                                (*weights)[i*nbins+j] += (sqrt(fabs((q_high)*(q_high)-q_shift * q_shift)))/q_0 * (2.0*double(npts_w)-1.0);
+                                                }
+                                                // For the rest of sub_width. weight += u(at q_high bin) - u(at q_low bin)
+                                                else if (u > 0.0){
+                                                        (*weights)[i*nbins+j] += (sqrt(fabs((q_high)*(q_high)- q_shift * q_shift))-sqrt(fabs((q_low)*(q_low)- q_shift * q_shift)))/q_0 ;
                                                 }
                                         }
                                 }
-
-                                // Skip the entries outside our I(q) range
-                                //TODO: be careful with edge effect
-                                if(q_i<nbins)
-                                        (*weights)[i*nbins+q_i]++;
                         }
                 }
@@ -277 +347 @@
 
                 for(int i=first_bin; i<=last_bin; i++){
-                        // Skip if weight is less than 1e-04(this value is much smaller than
+                        // Skip if weight is less than 1e-03(this value is much smaller than
                         // the weight at the 3*sigma distance
                         // Will speed up a little bit...
-                        if ((*weights)[q_i*nbins+i] < 1.0e-004){
+                        if ((*weights)[q_i*nbins+i] < 1.0e-003){
                                 continue;
                         }
@@ -288 +358 @@
 
                 // Normalize counts
-                iq_out[q_i] = (counts>0.0) ? sum/counts : 0;
+                iq_out[q_i] = (counts>0.0) ? sum/counts : 0.0;
         }
 }

DataLoader/extensions/smearer.hh

@@ -61 +61 @@
 protected:
     // Number of points used in the smearing computation
-    static const int npts   = 3000;
+    static const int npts   = 500;
 
 public:

DataLoader/qsmearing.py

@@ -157 +157 @@
         Perform smearing
         """
+        import time
+        st = time.time()
         # If this is the first time we call for smearing,
         # initialize the C++ smearer object first
@@ -174 +176 @@
             temp_first, temp_last = self._get_extrapolated_bin( \
                                                         first_bin, last_bin)
-            if self.nbins_low > 1:
+            if self.nbins_low > 0:
                 iq_in_low = self.model.evalDistribution( \
                                     numpy.fabs(self.qvalues[0:self.nbins_low]))
@@ -213 +215 @@
         temp_last = self.nbins - self.nbins_high
         out = iq_out[temp_first: temp_last]
-
+        print "time =", time.time() - st
         return out
     
@@ -549 +551 @@
     # Make a new qx array
     if nbins_low > 0:  
-        data_x_ext = numpy.append(extra_low, data_x)
+        data_x_ext = numpy.append(extra_low, data_x)
     else:
-        data_x_ext = data_x
+        data_x_ext = data_x
     data_x_ext = numpy.append(data_x_ext, extra_high)
     

DataLoader/test/utest_smearing.py

@@ -9 +9 @@
 from DataLoader.data_info import Data1D, Data2D
 from DataLoader.qsmearing import SlitSmearer, QSmearer, smear_selection
- 
+from sans.models.SphereModel import SphereModel
 import os.path
 
@@ -54 +54 @@
         input = 12.0-numpy.arange(1,11)
         output = s(input)
-        # The following commented line was the correct output for even bins [see smearer.cpp for details] 
-        #answer = [ 9.666,  9.056,  8.329,  7.494,  6.642,  5.721,  4.774,  3.824,  2.871, 2.   ]
-        answer = [ 9.2302,  8.6806,  7.9533,  7.1673,  6.2889,  5.4,     4.5028,  3.5744,  2.6083, 2.    ]
+        # The following commented line was the correct output for even bins 
+        # [see smearer.cpp for details] 
+        #answer = [ 9.666,  9.056,  8.329,  7.494,  6.642,  5.721,  4.774, \
+        #  3.824,  2.871, 2.   ]
+        # The following answer was from numerical weighting algorithm.
+        #answer = [ 9.2302,  8.6806,  7.9533,  7.1673,  6.2889,  5.4,   \
+        #  4.5028,  3.5744,  2.6083, 2.    ]
+        # For the new analytical algorithm, the small difference between 
+        #these two could be from the first edge of the q bin size. 
+        answer = [ 9.0618,  8.64018,  8.11868,  7.13916,  6.15285,  5.55556,  \
+                     4.55842,  3.56061,  2.56235, 2.    ]
         for i in range(len(input)):
             self.assertAlmostEqual(answer[i], output[i], 2)
@@ -89 +97 @@
                   5.00093415,   4.01898292,   3.15008701,   2.55214921]
         for i in range(len(input)):
-            self.assertAlmostEqual(answer[i], output[i], 4)
-      
-
+            self.assertAlmostEqual(answer[i], output[i], 2)
+            
+class smear_slit_h_w_tests(unittest.TestCase):
+    
+    def setUp(self):
+        self.data = Loader().load("1000A_sphere_sm.xml")
+        self.model = SphereModel()
+        # The answer could be improved by developing better algorithm.
+        self.answer1 = Loader().load("slit_1000A_sphere_sm_w_0_0002.txt")
+        self.answer2 = Loader().load("slit_1000A_sphere_sm_h.txt")
+        self.answer3 = Loader().load("slit_1000A_sphere_sm_w_0_0001.txt")
+        # Get inputs
+        self.model.params['scale'] = 0.05
+        self.model.params['background'] = 0.01
+        self.model.params['radius'] = 10000.0
+        self.model.params['sldSolv'] = 6.3e-006
+        self.model.params['sldSph'] = 3.4e-006
+        
+    def test_slit_h_w(self):
+        """
+            Test identity slit smearing w/ h=0.117 w = 0.002
+        """
+        # Set params and dQl
+        data = self.data
+        data.dxw = 0.0002 * numpy.ones(len(self.data.x))
+        data.dxl = 0.117 * numpy.ones(len(self.data.x))
+        # Create smearer for our data
+        s = SlitSmearer(data, self.model)
+        # Get smear
+        input  = self.model.evalDistribution(data.x)
+        output = s(input)
+        # Get pre-calculated values
+        answer = self.answer1.y
+        # Now compare
+        for i in range(len(input)):
+            self.assertAlmostEqual(answer[i], output[i], 0)
+     
+    def test_slit_h(self):
+        """
+            Test identity slit smearing w/ h=0.117 w = 0.0
+        """
+        # Set params and dQl
+        data = self.data
+        data.dxw = 0.0 * numpy.ones(len(self.data.x))
+        data.dxl = 0.117 * numpy.ones(len(self.data.x))
+        # Create smearer for our data
+        s = SlitSmearer(data, self.model)
+        # Get smear
+        input  = self.model.evalDistribution(data.x)
+        output = s(input)
+        # Get pre-calculated values
+        answer = self.answer2.y
+        # Now compare
+        for i in range(len(input)):
+            self.assertAlmostEqual(answer[i], output[i], 0)
+     
+    def test_slit_w(self):
+        """
+            Test identity slit smearing w/ h=0.0 w = 0.001
+        """
+        # Set params and dQl
+        data = self.data
+        data.dxw = 0.0001 * numpy.ones(len(self.data.x))
+        data.dxl = 0.0 * numpy.ones(len(self.data.x))
+        # Create smearer for our data
+        s = SlitSmearer(data, self.model)
+        # Get smear
+        input  = self.model.evalDistribution(data.x)
+        output = s(input)
+        # Get pre-calculated values
+        answer = self.answer3.y
+        # Now compare
+        for i in range(len(input)):
+             if i <= 40:
+                 self.assertAlmostEqual(answer[i], output[i], -3)    
+             else:
+                self.assertAlmostEqual(answer[i], output[i], 0)       
+                
 if __name__ == '__main__':
     unittest.main()