Changeset cd2ced80 in sasview for DataLoader
- Timestamp:
- Jan 14, 2011 5:45:58 PM (14 years ago)
- 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
- Location:
- DataLoader
- Files:
-
- 3 added
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
DataLoader/extensions/smearer.cpp
ra74d650a rcd2ced80 133 133 // Check the length of the data 134 134 if (nbins<2) return; 135 int npts_h = height>0.0 ? npts : 1; 136 int npts_w = width>0.0 ? npts : 1; 137 138 // If both height and width are great than zero, 139 // modify the number of points in each direction so 140 // that the total number of points is still what 141 // the user would expect (downgrade resolution) 142 //if(npts_h>1 && npts_w>1){ 143 // npts_h = (int)ceil(sqrt((double)npts)); 144 // npts_w = npts_h; 145 //} 146 double shift_h, shift_w, hbin_size, wbin_size; 147 // Make sure height and width are all positive (FWMH/2) 148 // Assumption; height and width are all same for all q points 149 if(npts_h == 1){ 150 shift_h = 0.0; 151 } else { 152 shift_h = fabs(height); 153 } 154 if(npts_w == 1){ 155 shift_w = 0.0; 156 } else { 157 shift_w = fabs(width); 158 } 159 // size of the h bin and w bin 160 hbin_size = shift_h / nbins; 161 wbin_size = shift_w / nbins; 135 162 136 163 // Loop over all q-values 137 164 for(int i=0; i<nbins; i++) { 138 double q, q_min, q_max; 165 // Find Weights 166 // Find q where the resolution smearing calculation of I(q) occurs 167 double q, q_min, q_max, q_0; 139 168 get_bin_range(i, &q, &q_min, &q_max); 140 141 // For each q-value, compute the weight of each other q-bin 142 // in the I(q) array 143 int npts_h = height>0 ? npts : 1; 144 int npts_w = width>0 ? npts : 1; 145 146 // If both height and width are great than zero, 147 // modify the number of points in each direction so 148 // that the total number of points is still what 149 // the user would expect (downgrade resolution) 150 // Never down-grade npts_h. That will give incorrect slit smearing... 151 if(npts_h>1 && npts_w>1){ 152 npts_h = npts;//(int)ceil(sqrt((double)npts)); 153 // In general width is much smaller than height, so smaller npt_w 154 // should work fine. 155 // Todo: It is still very expansive in time. Think about better way. 156 npts_w = (int)ceil(npts_h / 100); 157 } 158 159 double shift_h, shift_w; 160 for(int k=0; k<npts_h; k++){ 161 if(npts_h==1){ 162 shift_h = 0; 163 } else { 164 shift_h = height/((double)npts_h-1.0) * (double)k; 165 } 166 for(int j=0; j<npts_w; j++){ 167 if(npts_w==1){ 168 shift_w = 0; 169 } else { 170 shift_w = width/((double)npts_w-1.0) * (double)j; 169 bool last_qpoint = true; 170 // Find q[0] value to normalize the weight later, 171 // otherwise, we will have a precision problem. 172 if (i == 0){ 173 q_0 = q; 174 } 175 // Loop over all qj-values 176 bool first_w = true; 177 for(int j=0; j<nbins; j++) { 178 double q_j, q_high, q_low; 179 // Calculate bin size of q_j 180 get_bin_range(j, &q_j, &q_low, &q_high); 181 // Check q_low that can not be negative. 182 if (q_low < 0.0){ 183 q_low = 0.0; 184 } 185 // default parameter values 186 (*weights)[i*nbins+j] = 0.0; 187 double shift_w = 0.0; 188 // Condition: zero slit smear. 189 if (npts_w == 1 && npts_h == 1){ 190 if(q_j == q) { 191 (*weights)[i*nbins+j] = 1.0; 171 192 } 172 double q_shifted = sqrt( ((q-shift_w)*(q-shift_w) + shift_h*shift_h) ); 173 174 // Find in which bin this shifted value belongs 175 int q_i=nbins; 176 if (even_binning) { 177 // This is kept for backward compatibility since the binning 178 // was originally defined differently for even bins. 179 q_i = (int)(floor( (q_shifted-qmin) /((qmax-qmin)/((double)nbins -1.0)) )); 180 } else { 181 for(int t=0; t<nbins; t++) { 182 double q_t, q_high, q_low; 183 get_bin_range(t, &q_t, &q_low, &q_high); 184 if(q_shifted>=q_low && q_shifted<q_high) { 185 q_i = t; 186 break; 193 } 194 //Condition:Smear weight integration for width >0 when the height (=0) does not present. 195 //Or height << width. 196 else if((npts_w!=1 && npts_h == 1)|| (npts_w!=1 && npts_h != 1 && width/height > 100.0)){ 197 shift_w = width; 198 //del_w = width/((double)npts_w-1.0); 199 double q_shifted_low = q - shift_w; 200 // High limit of the resolution range 201 double q_shifted_high = q + shift_w; 202 // Go through all the q_js for weighting those points 203 if(q_j >= q_shifted_low && q_j <= q_shifted_high) { 204 // The weighting factor comes, 205 // Give some weight (delq_bin) for the q_j within the resolution range 206 // Weight should be same for all qs except 207 // for the q bin size at j. 208 // Note that the division by q_0 is only due to the precision problem 209 // where q_high - q_low gets to very small. 210 // Later, it will be normalized again. 211 (*weights)[i*nbins+j] += (q_high - q_low)/q_0 ; 212 } 213 } 214 else{ 215 // Loop for width (;Height is analytical.) 216 // Condition: height >>> width, otherwise, below is not accurate enough. 217 // Smear weight numerical iteration for width >0 when the height (>0) presents. 218 // When width = 0, the numerical iteration will be skipped. 219 // The resolution calculation for the height is done by direct integration, 220 // assuming the I(q'=sqrt(q_j^2-(q+shift_w)^2)) is constant within a q' bin, [q_high, q_low]. 221 // In general, this weight numerical iteration for width >0 might be a rough approximation, 222 // but it must be good enough when height >>> width. 223 for(int k=(-npts_w + 1); k<npts_w; k++){ 224 if(npts_w!=1){ 225 shift_w = width/((double)npts_w-1.0)*(double)k; 226 } 227 // For each q-value, compute the weight of each other q-bin 228 // in the I(q) array 229 // Low limit of the resolution range 230 double q_shift = q + shift_w; 231 if (q_shift < 0.0){ 232 q_shift = 0.0; 233 } 234 double q_shifted_low = q_shift; 235 // High limit of the resolution range 236 double q_shifted_high = sqrt(q_shift * q_shift + shift_h * shift_h); 237 238 239 // Go through all the q_js for weighting those points 240 if(q_j >= q_shifted_low && q_j <= q_shifted_high) { 241 // The weighting factor comes, 242 // Give some weight (delq_bin) for the q_j within the resolution range 243 // Weight should be same for all qs except 244 // for the q bin size at j. 245 // Note that the division by q_0 is only due to the precision problem 246 // where q_high - q_low gets to very small. 247 // Later, it will be normalized again. 248 249 double q_shift_min = q - width; 250 251 double u = (q_j * q_j - (q_shift) * (q_shift)); 252 // The fabs below are not necessary but in case: the weight should never be imaginary. 253 // At the edge of each sub_width. weight += u(at q_high bin) - u(0), where u(0) = 0, 254 // and weighted by (2.0* npts_w -1.0)once for each q. 255 if (q == q_j) { 256 if (k==0) 257 (*weights)[i*nbins+j] += (sqrt(fabs((q_high)*(q_high)-q_shift * q_shift)))/q_0 * (2.0*double(npts_w)-1.0); 258 } 259 // For the rest of sub_width. weight += u(at q_high bin) - u(at q_low bin) 260 else if (u > 0.0){ 261 (*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 ; 187 262 } 188 263 } 189 264 } 190 191 // Skip the entries outside our I(q) range192 //TODO: be careful with edge effect193 if(q_i<nbins)194 (*weights)[i*nbins+q_i]++;195 265 } 196 266 } … … 277 347 278 348 for(int i=first_bin; i<=last_bin; i++){ 279 // Skip if weight is less than 1e-0 4(this value is much smaller than349 // Skip if weight is less than 1e-03(this value is much smaller than 280 350 // the weight at the 3*sigma distance 281 351 // Will speed up a little bit... 282 if ((*weights)[q_i*nbins+i] < 1.0e-00 4){352 if ((*weights)[q_i*nbins+i] < 1.0e-003){ 283 353 continue; 284 354 } … … 288 358 289 359 // Normalize counts 290 iq_out[q_i] = (counts>0.0) ? sum/counts : 0 ;360 iq_out[q_i] = (counts>0.0) ? sum/counts : 0.0; 291 361 } 292 362 } -
DataLoader/extensions/smearer.hh
r2ca00c4 rcd2ced80 61 61 protected: 62 62 // Number of points used in the smearing computation 63 static const int npts = 3000;63 static const int npts = 500; 64 64 65 65 public: -
DataLoader/qsmearing.py
r2ca00c4 rcd2ced80 157 157 Perform smearing 158 158 """ 159 import time 160 st = time.time() 159 161 # If this is the first time we call for smearing, 160 162 # initialize the C++ smearer object first … … 174 176 temp_first, temp_last = self._get_extrapolated_bin( \ 175 177 first_bin, last_bin) 176 if self.nbins_low > 1:178 if self.nbins_low > 0: 177 179 iq_in_low = self.model.evalDistribution( \ 178 180 numpy.fabs(self.qvalues[0:self.nbins_low])) … … 213 215 temp_last = self.nbins - self.nbins_high 214 216 out = iq_out[temp_first: temp_last] 215 217 print "time =", time.time() - st 216 218 return out 217 219 … … 549 551 # Make a new qx array 550 552 if nbins_low > 0: 551 553 data_x_ext = numpy.append(extra_low, data_x) 552 554 else: 553 555 data_x_ext = data_x 554 556 data_x_ext = numpy.append(data_x_ext, extra_high) 555 557 -
DataLoader/test/utest_smearing.py
rf867cd9 rcd2ced80 9 9 from DataLoader.data_info import Data1D, Data2D 10 10 from DataLoader.qsmearing import SlitSmearer, QSmearer, smear_selection 11 11 from sans.models.SphereModel import SphereModel 12 12 import os.path 13 13 … … 54 54 input = 12.0-numpy.arange(1,11) 55 55 output = s(input) 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. ] 56 # The following commented line was the correct output for even bins 57 # [see smearer.cpp for details] 58 #answer = [ 9.666, 9.056, 8.329, 7.494, 6.642, 5.721, 4.774, \ 59 # 3.824, 2.871, 2. ] 60 # The following answer was from numerical weighting algorithm. 61 #answer = [ 9.2302, 8.6806, 7.9533, 7.1673, 6.2889, 5.4, \ 62 # 4.5028, 3.5744, 2.6083, 2. ] 63 # For the new analytical algorithm, the small difference between 64 #these two could be from the first edge of the q bin size. 65 answer = [ 9.0618, 8.64018, 8.11868, 7.13916, 6.15285, 5.55556, \ 66 4.55842, 3.56061, 2.56235, 2. ] 59 67 for i in range(len(input)): 60 68 self.assertAlmostEqual(answer[i], output[i], 2) … … 89 97 5.00093415, 4.01898292, 3.15008701, 2.55214921] 90 98 for i in range(len(input)): 91 self.assertAlmostEqual(answer[i], output[i], 4) 92 93 99 self.assertAlmostEqual(answer[i], output[i], 2) 100 101 class smear_slit_h_w_tests(unittest.TestCase): 102 103 def setUp(self): 104 self.data = Loader().load("1000A_sphere_sm.xml") 105 self.model = SphereModel() 106 # The answer could be improved by developing better algorithm. 107 self.answer1 = Loader().load("slit_1000A_sphere_sm_w_0_0002.txt") 108 self.answer2 = Loader().load("slit_1000A_sphere_sm_h.txt") 109 self.answer3 = Loader().load("slit_1000A_sphere_sm_w_0_0001.txt") 110 # Get inputs 111 self.model.params['scale'] = 0.05 112 self.model.params['background'] = 0.01 113 self.model.params['radius'] = 10000.0 114 self.model.params['sldSolv'] = 6.3e-006 115 self.model.params['sldSph'] = 3.4e-006 116 117 def test_slit_h_w(self): 118 """ 119 Test identity slit smearing w/ h=0.117 w = 0.002 120 """ 121 # Set params and dQl 122 data = self.data 123 data.dxw = 0.0002 * numpy.ones(len(self.data.x)) 124 data.dxl = 0.117 * numpy.ones(len(self.data.x)) 125 # Create smearer for our data 126 s = SlitSmearer(data, self.model) 127 # Get smear 128 input = self.model.evalDistribution(data.x) 129 output = s(input) 130 # Get pre-calculated values 131 answer = self.answer1.y 132 # Now compare 133 for i in range(len(input)): 134 self.assertAlmostEqual(answer[i], output[i], 0) 135 136 def test_slit_h(self): 137 """ 138 Test identity slit smearing w/ h=0.117 w = 0.0 139 """ 140 # Set params and dQl 141 data = self.data 142 data.dxw = 0.0 * numpy.ones(len(self.data.x)) 143 data.dxl = 0.117 * numpy.ones(len(self.data.x)) 144 # Create smearer for our data 145 s = SlitSmearer(data, self.model) 146 # Get smear 147 input = self.model.evalDistribution(data.x) 148 output = s(input) 149 # Get pre-calculated values 150 answer = self.answer2.y 151 # Now compare 152 for i in range(len(input)): 153 self.assertAlmostEqual(answer[i], output[i], 0) 154 155 def test_slit_w(self): 156 """ 157 Test identity slit smearing w/ h=0.0 w = 0.001 158 """ 159 # Set params and dQl 160 data = self.data 161 data.dxw = 0.0001 * numpy.ones(len(self.data.x)) 162 data.dxl = 0.0 * numpy.ones(len(self.data.x)) 163 # Create smearer for our data 164 s = SlitSmearer(data, self.model) 165 # Get smear 166 input = self.model.evalDistribution(data.x) 167 output = s(input) 168 # Get pre-calculated values 169 answer = self.answer3.y 170 # Now compare 171 for i in range(len(input)): 172 if i <= 40: 173 self.assertAlmostEqual(answer[i], output[i], -3) 174 else: 175 self.assertAlmostEqual(answer[i], output[i], 0) 176 94 177 if __name__ == '__main__': 95 178 unittest.main()
Note: See TracChangeset
for help on using the changeset viewer.