Changeset f8aa738 in sasview for src/sas/sascalc/data_util
- Timestamp:
- Mar 21, 2016 5:17:48 AM (9 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:
- 80b1df3
- Parents:
- 240a2d2
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
src/sas/sascalc/data_util/qsmearing.py
r240a2d2 rf8aa738 13 13 import logging 14 14 import sys 15 import sas.models.sas_extension.smearer as smearer16 from sas.sascalc.data_util.smearing_2d import Smearer2D17 15 18 def smear_selection(data1D, model = None): 16 from sasmodels.resolution import Slit1D, Pinhole1D 17 from sasmodels.resolution2d import Pinhole2D 18 19 def smear_selection(data, model = None): 19 20 """ 20 Creates the right type of smearer according 21 Creates the right type of smearer according 21 22 to the data. 22 23 23 The canSAS format has a rule that either 24 24 slit smearing data OR resolution smearing data 25 is available. 26 25 is available. 26 27 27 For the present purpose, we choose the one that 28 28 has none-zero data. If both slit and resolution 29 smearing arrays are filled with good data 29 smearing arrays are filled with good data 30 30 (which should not happen), then we choose the 31 resolution smearing data. 32 33 :param data 1D: Data1D object31 resolution smearing data. 32 33 :param data: Data1D object 34 34 :param model: sas.model instance 35 35 """ 36 36 # Sanity check. If we are not dealing with a SAS Data1D 37 37 # object, just return None 38 if data 1D.__class__.__name__ not in ['Data1D', 'Theory1D']:39 if data 1D== None:38 if data.__class__.__name__ not in ['Data1D', 'Theory1D']: 39 if data == None: 40 40 return None 41 elif data 1D.dqx_data == None or data1D.dqy_data == None:41 elif data.dqx_data == None or data.dqy_data == None: 42 42 return None 43 return Smearer2D(data1D)44 45 if not hasattr(data 1D, "dx") and not hasattr(data1D, "dxl")\46 and not hasattr(data 1D, "dxw"):43 return Pinhole2D(data) 44 45 if not hasattr(data, "dx") and not hasattr(data, "dxl")\ 46 and not hasattr(data, "dxw"): 47 47 return None 48 48 49 49 # Look for resolution smearing data 50 50 _found_resolution = False 51 if data 1D.dx is not None and len(data1D.dx) == len(data1D.x):52 51 if data.dx is not None and len(data.dx) == len(data.x): 52 53 53 # Check that we have non-zero data 54 if data 1D.dx[0] > 0.0:54 if data.dx[0] > 0.0: 55 55 _found_resolution = True 56 56 #print "_found_resolution",_found_resolution … … 58 58 # If we found resolution smearing data, return a QSmearer 59 59 if _found_resolution == True: 60 return QSmearer(data1D, model) 61 #return pinhole_smear(data1D, model) 60 return pinhole_smear(data, model) 62 61 63 62 # Look for slit smearing data 64 63 _found_slit = False 65 if data 1D.dxl is not None and len(data1D.dxl) == len(data1D.x) \66 and data 1D.dxw is not None and len(data1D.dxw) == len(data1D.x):67 64 if data.dxl is not None and len(data.dxl) == len(data.x) \ 65 and data.dxw is not None and len(data.dxw) == len(data.x): 66 68 67 # Check that we have non-zero data 69 if data 1D.dxl[0] > 0.0 or data1D.dxw[0] > 0.0:68 if data.dxl[0] > 0.0 or data.dxw[0] > 0.0: 70 69 _found_slit = True 71 70 72 71 # Sanity check: all data should be the same as a function of Q 73 for item in data 1D.dxl:74 if data 1D.dxl[0] != item:72 for item in data.dxl: 73 if data.dxl[0] != item: 75 74 _found_resolution = False 76 75 break 77 78 for item in data 1D.dxw:79 if data 1D.dxw[0] != item:76 77 for item in data.dxw: 78 if data.dxw[0] != item: 80 79 _found_resolution = False 81 80 break 82 81 # If we found slit smearing data, return a slit smearer 83 82 if _found_slit == True: 84 #return SlitSmearer(data1D, model) 85 return slit_smear(data1D, model) 83 return slit_smear(data, model) 86 84 return None 87 88 89 class _BaseSmearer(object):90 """91 Base class for smearers92 """93 def __init__(self):94 self.nbins = 095 self.nbins_low = 096 self.nbins_high = 097 self._weights = None98 ## Internal flag to keep track of C++ smearer initialization99 self._init_complete = False100 self._smearer = None101 self.model = None102 self.min = None103 self.max = None104 self.qvalues = []105 106 def __deepcopy__(self, memo=None):107 """108 Return a valid copy of self.109 Avoid copying the _smearer C object and force a matrix recompute110 when the copy is used.111 """112 result = _BaseSmearer()113 result.nbins = self.nbins114 return result115 116 def _compute_matrix(self):117 """118 Place holder for matrix computation119 """120 return NotImplemented121 122 def get_unsmeared_range(self, q_min=None, q_max=None):123 """124 Place holder for method returning unsmeared range125 """126 return NotImplemented127 128 def get_bin_range(self, q_min=None, q_max=None):129 """130 131 :param q_min: minimum q-value to smear132 :param q_max: maximum q-value to smear133 134 """135 # If this is the first time we call for smearing,136 # initialize the C++ smearer object first137 if not self._init_complete:138 self._initialize_smearer()139 if q_min == None:140 q_min = self.min141 if q_max == None:142 q_max = self.max143 144 _qmin_unsmeared, _qmax_unsmeared = self.get_unsmeared_range(q_min,145 q_max)146 _first_bin = None147 _last_bin = None148 149 #step = (self.max - self.min) / (self.nbins - 1.0)150 # Find the first and last bin number in all extrapolated and real data151 try:152 for i in range(self.nbins):153 q_i = smearer.get_q(self._smearer, i)154 if (q_i >= _qmin_unsmeared) and (q_i <= _qmax_unsmeared):155 # Identify first and last bin156 if _first_bin is None:157 _first_bin = i158 else:159 _last_bin = i160 except:161 msg = "_BaseSmearer.get_bin_range: "162 msg += " error getting range\n %s" % sys.exc_value163 raise RuntimeError, msg164 165 # Find the first and last bin number only in the real data166 _first_bin, _last_bin = self._get_unextrapolated_bin( \167 _first_bin, _last_bin)168 169 return _first_bin, _last_bin170 171 def __call__(self, iq_in, first_bin = 0, last_bin = None):172 """173 Perform smearing174 """175 # If this is the first time we call for smearing,176 # initialize the C++ smearer object first177 if not self._init_complete:178 self._initialize_smearer()179 180 if last_bin is None or last_bin >= len(iq_in):181 last_bin = len(iq_in) - 1182 # Check that the first bin is positive183 if first_bin < 0:184 first_bin = 0185 186 # With a model given, compute I for the extrapolated points and append187 # to the iq_in188 iq_in_temp = iq_in189 if self.model != None:190 temp_first, temp_last = self._get_extrapolated_bin( \191 first_bin, last_bin)192 if self.nbins_low > 0:193 iq_in_low = self.model.evalDistribution( \194 numpy.fabs(self.qvalues[0:self.nbins_low]))195 iq_in_high = self.model.evalDistribution( \196 self.qvalues[(len(self.qvalues) - \197 self.nbins_high - 1):])198 # Todo: find out who is sending iq[last_poin] = 0.199 if iq_in[len(iq_in) - 1] == 0:200 iq_in[len(iq_in) - 1] = iq_in_high[0]201 # Append the extrapolated points to the data points202 if self.nbins_low > 0:203 iq_in_temp = numpy.append(iq_in_low, iq_in)204 if self.nbins_high > 0:205 iq_in_temp = numpy.append(iq_in_temp, iq_in_high[1:])206 else:207 temp_first = first_bin208 temp_last = last_bin209 #iq_in_temp = iq_in210 211 # Sanity check212 if len(iq_in_temp) != self.nbins:213 msg = "Invalid I(q) vector: inconsistent array "214 msg += " length %d != %s" % (len(iq_in_temp), str(self.nbins))215 raise RuntimeError, msg216 217 # Storage for smeared I(q)218 iq_out = numpy.zeros(self.nbins)219 220 smear_output = smearer.smear(self._smearer, iq_in_temp, iq_out,221 #0, self.nbins - 1)222 temp_first, temp_last)223 #first_bin, last_bin)224 if smear_output < 0:225 msg = "_BaseSmearer: could not smear, code = %g" % smear_output226 raise RuntimeError, msg227 228 temp_first = first_bin + self.nbins_low229 temp_last = self.nbins - self.nbins_high230 out = iq_out[temp_first: temp_last]231 232 return out233 234 def _initialize_smearer(self):235 """236 Place holder for initializing data smearer237 """238 return NotImplemented239 240 241 def _get_unextrapolated_bin(self, first_bin = 0, last_bin = 0):242 """243 Get unextrapolated first bin and the last bin244 245 : param first_bin: extrapolated first_bin246 : param last_bin: extrapolated last_bin247 248 : return fist_bin, last_bin: unextrapolated first and last bin249 """250 # For first bin251 if first_bin <= self.nbins_low:252 first_bin = 0253 else:254 first_bin = first_bin - self.nbins_low255 # For last bin256 if last_bin >= (self.nbins - self.nbins_high):257 last_bin = self.nbins - (self.nbins_high + self.nbins_low + 1)258 elif last_bin >= self.nbins_low:259 last_bin = last_bin - self.nbins_low260 else:261 last_bin = 0262 return first_bin, last_bin263 264 def _get_extrapolated_bin(self, first_bin = 0, last_bin = 0):265 """266 Get extrapolated first bin and the last bin267 268 : param first_bin: unextrapolated first_bin269 : param last_bin: unextrapolated last_bin270 271 : return first_bin, last_bin: extrapolated first and last bin272 """273 # For the first bin274 # In the case that needs low extrapolation data275 first_bin = 0276 # For last bin277 if last_bin >= self.nbins - (self.nbins_high + self.nbins_low + 1):278 # In the case that needs higher q extrapolation data279 last_bin = self.nbins - 1280 else:281 # In the case that doesn't need higher q extrapolation data282 last_bin += self.nbins_low283 284 return first_bin, last_bin285 286 class _SlitSmearer(_BaseSmearer):287 """288 Slit smearing for I(q) array289 """290 291 def __init__(self, nbins=None, width=None, height=None, min=None, max=None):292 """293 Initialization294 295 :param iq: I(q) array [cm-1]296 :param width: slit width [A-1]297 :param height: slit height [A-1]298 :param min: Q_min [A-1]299 :param max: Q_max [A-1]300 301 """302 _BaseSmearer.__init__(self)303 ## Slit width in Q units304 self.width = width305 ## Slit height in Q units306 self.height = height307 ## Q_min (Min Q-value for I(q))308 self.min = min309 ## Q_max (Max Q_value for I(q))310 self.max = max311 ## Number of Q bins312 self.nbins = nbins313 ## Number of points used in the smearing computation314 self.npts = 3000315 ## Smearing matrix316 self._weights = None317 self.qvalues = None318 319 def _initialize_smearer(self):320 """321 Initialize the C++ smearer object.322 This method HAS to be called before smearing323 """324 #self._smearer = smearer.new_slit_smearer(self.width,325 # self.height, self.min, self.max, self.nbins)326 self._smearer = smearer.new_slit_smearer_with_q(self.width,327 self.height, self.qvalues)328 self._init_complete = True329 330 def get_unsmeared_range(self, q_min, q_max):331 """332 Determine the range needed in unsmeared-Q to cover333 the smeared Q range334 """335 # Range used for input to smearing336 _qmin_unsmeared = q_min337 _qmax_unsmeared = q_max338 try:339 _qmin_unsmeared = self.min340 _qmax_unsmeared = self.max341 except:342 logging.error("_SlitSmearer.get_bin_range: %s" % sys.exc_value)343 return _qmin_unsmeared, _qmax_unsmeared344 345 class SlitSmearer(_SlitSmearer):346 """347 Adaptor for slit smearing class and SAS data348 """349 def __init__(self, data1D, model = None):350 """351 Assumption: equally spaced bins of increasing q-values.352 353 :param data1D: data used to set the smearing parameters354 """355 # Initialization from parent class356 super(SlitSmearer, self).__init__()357 358 ## Slit width359 self.width = 0360 self.nbins_low = 0361 self.nbins_high = 0362 self.model = model363 if data1D.dxw is not None and len(data1D.dxw) == len(data1D.x):364 self.width = data1D.dxw[0]365 # Sanity check366 for value in data1D.dxw:367 if value != self.width:368 msg = "Slit smearing parameters must "369 msg += " be the same for all data"370 raise RuntimeError, msg371 ## Slit height372 self.height = 0373 if data1D.dxl is not None and len(data1D.dxl) == len(data1D.x):374 self.height = data1D.dxl[0]375 # Sanity check376 for value in data1D.dxl:377 if value != self.height:378 msg = "Slit smearing parameters must be"379 msg += " the same for all data"380 raise RuntimeError, msg381 # If a model is given, get the q extrapolation382 if self.model == None:383 data1d_x = data1D.x384 else:385 # Take larger sigma386 if self.height > self.width:387 # The denominator (2.0) covers all the possible w^2 + h^2 range388 sigma_in = data1D.dxl / 2.0389 elif self.width > 0:390 sigma_in = data1D.dxw / 2.0391 else:392 sigma_in = []393 394 self.nbins_low, self.nbins_high, _, data1d_x = \395 get_qextrapolate(sigma_in, data1D.x)396 397 ## Number of Q bins398 self.nbins = len(data1d_x)399 ## Minimum Q400 self.min = min(data1d_x)401 ## Maximum402 self.max = max(data1d_x)403 ## Q-values404 self.qvalues = data1d_x405 406 407 class _QSmearer(_BaseSmearer):408 """409 Perform Gaussian Q smearing410 """411 412 def __init__(self, nbins=None, width=None, min=None, max=None):413 """414 Initialization415 416 :param nbins: number of Q bins417 :param width: array standard deviation in Q [A-1]418 :param min: Q_min [A-1]419 :param max: Q_max [A-1]420 """421 _BaseSmearer.__init__(self)422 ## Standard deviation in Q [A-1]423 self.width = width424 ## Q_min (Min Q-value for I(q))425 self.min = min426 ## Q_max (Max Q_value for I(q))427 self.max = max428 ## Number of Q bins429 self.nbins = nbins430 ## Smearing matrix431 self._weights = None432 self.qvalues = None433 434 def _initialize_smearer(self):435 """436 Initialize the C++ smearer object.437 This method HAS to be called before smearing438 """439 #self._smearer = smearer.new_q_smearer(numpy.asarray(self.width),440 # self.min, self.max, self.nbins)441 self._smearer = smearer.new_q_smearer_with_q(numpy.asarray(self.width),442 self.qvalues)443 self._init_complete = True444 445 def get_unsmeared_range(self, q_min, q_max):446 """447 Determine the range needed in unsmeared-Q to cover448 the smeared Q range449 Take 3 sigmas as the offset between smeared and unsmeared space450 """451 # Range used for input to smearing452 _qmin_unsmeared = q_min453 _qmax_unsmeared = q_max454 try:455 offset = 3.0 * max(self.width)456 _qmin_unsmeared = self.min#max([self.min, q_min - offset])457 _qmax_unsmeared = self.max#min([self.max, q_max + offset])458 except:459 logging.error("_QSmearer.get_bin_range: %s" % sys.exc_value)460 return _qmin_unsmeared, _qmax_unsmeared461 462 463 class QSmearer(_QSmearer):464 """465 Adaptor for Gaussian Q smearing class and SAS data466 """467 def __init__(self, data1D, model = None):468 """469 Assumption: equally spaced bins of increasing q-values.470 471 :param data1D: data used to set the smearing parameters472 """473 # Initialization from parent class474 super(QSmearer, self).__init__()475 data1d_x = []476 self.nbins_low = 0477 self.nbins_high = 0478 self.model = model479 ## Resolution480 #self.width = numpy.zeros(len(data1D.x))481 if data1D.dx is not None and len(data1D.dx) == len(data1D.x):482 self.width = data1D.dx483 484 if self.model == None:485 data1d_x = data1D.x486 else:487 self.nbins_low, self.nbins_high, self.width, data1d_x = \488 get_qextrapolate(self.width, data1D.x)489 490 ## Number of Q bins491 self.nbins = len(data1d_x)492 ## Minimum Q493 self.min = min(data1d_x)494 ## Maximum495 self.max = max(data1d_x)496 ## Q-values497 self.qvalues = data1d_x498 499 500 def get_qextrapolate(width, data_x):501 """502 Make fake data_x points extrapolated outside of the data_x points503 504 :param width: array of std of q resolution505 :param Data1D.x: Data1D.x array506 507 :return new_width, data_x_ext: extrapolated width array and x array508 509 :assumption1: data_x is ordered from lower q to higher q510 :assumption2: len(data) = len(width)511 :assumption3: the distance between the data points is more compact than the size of width512 :Todo1: Make sure that the assumptions are correct for Data1D513 :Todo2: This fixes the edge problem in Qsmearer but still needs to make smearer interface514 """515 # Length of the width516 length = len(width)517 width_low = math.fabs(width[0])518 width_high = math.fabs(width[length -1])519 nbins_low = 0.0520 nbins_high = 0.0521 # Compare width(dQ) to the data bin size and take smaller one as the bin522 # size of the extrapolation; this will correct some weird behavior523 # at the edge: This method was out (commented)524 # because it becomes very expansive when525 # bin size is very small comparing to the width.526 # Now on, we will just give the bin size of the extrapolated points527 # based on the width.528 # Find bin sizes529 #bin_size_low = math.fabs(data_x[1] - data_x[0])530 #bin_size_high = math.fabs(data_x[length - 1] - data_x[length - 2])531 # Let's set the bin size 1/3 of the width(sigma), it is good as long as532 # the scattering is monotonous.533 #if width_low < (bin_size_low):534 bin_size_low = width_low / 10.0535 #if width_high < (bin_size_high):536 bin_size_high = width_high / 10.0537 538 # Number of q points required below the 1st data point in order to extend539 # them 3 times of the width (std)540 if width_low > 0.0:541 nbins_low = math.ceil(3.0 * width_low / bin_size_low)542 # Number of q points required above the last data point543 if width_high > 0.0:544 nbins_high = math.ceil(3.0 * width_high / bin_size_high)545 # Make null q points546 extra_low = numpy.zeros(nbins_low)547 extra_high = numpy.zeros(nbins_high)548 # Give extrapolated values549 ind = 0550 qvalue = data_x[0] - bin_size_low551 #if qvalue > 0:552 while(ind < nbins_low):553 extra_low[nbins_low - (ind + 1)] = qvalue554 qvalue -= bin_size_low555 ind += 1556 #if qvalue <= 0:557 # break558 # Redefine nbins_low559 nbins_low = ind560 # Reset ind for another extrapolation561 ind = 0562 qvalue = data_x[length -1] + bin_size_high563 while(ind < nbins_high):564 extra_high[ind] = qvalue565 qvalue += bin_size_high566 ind += 1567 # Make a new qx array568 if nbins_low > 0:569 data_x_ext = numpy.append(extra_low, data_x)570 else:571 data_x_ext = data_x572 data_x_ext = numpy.append(data_x_ext, extra_high)573 574 # Redefine extra_low and high based on corrected nbins575 # And note that it is not necessary for extra_width to be a non-zero576 if nbins_low > 0:577 extra_low = numpy.zeros(nbins_low)578 extra_high = numpy.zeros(nbins_high)579 # Make new width array580 new_width = numpy.append(extra_low, width)581 new_width = numpy.append(new_width, extra_high)582 583 # nbins corrections due to the negative q value584 nbins_low = nbins_low - len(data_x_ext[data_x_ext <= 0])585 return nbins_low, nbins_high, \586 new_width[data_x_ext > 0], data_x_ext[data_x_ext > 0]587 85 588 86 589 590 from .resolution import Slit1D, Pinhole1D591 87 class PySmear(object): 592 88 """ … … 601 97 """ 602 98 Apply the resolution function to the data. 603 604 99 Note that this is called with iq_in matching data.x, but with 605 100 iq_in[first_bin:last_bin] set to theory values for these bins, 606 101 and the remainder left undefined. The first_bin, last_bin values 607 102 should be those returned from get_bin_range. 608 609 103 The returned value is of the same length as iq_in, with the range 610 104 first_bin:last_bin set to the resolution smeared values. … … 626 120 """ 627 121 For a given q_min, q_max, find the corresponding indices in the data. 628 629 122 Returns first, last. 630 631 123 Note that these are indexes into q from the data, not the q_calc 632 124 needed by the resolution function. Note also that these are the
Note: See TracChangeset
for help on using the changeset viewer.