- Timestamp:
- Dec 15, 2009 11:22:39 AM (15 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:
- 8d1e745
- Parents:
- 39e49a1
- Location:
- Invariant
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
Invariant/invariant.py
r3bb37ef r4e80ae0 31 31 return scale * value 32 32 33 def power_law(x, scale=1, power= 4):33 def power_law(x, scale=1, power=0): 34 34 """ 35 35 F(x) = scale* (x)^(-power) … … 42 42 @param F(x) 43 43 """ 44 if power <=0: 45 raise ValueError("Power_law function expected posive power, but got %s"%power) 46 44 47 value = numpy.array([ math.pow(i, -power) for i in x ]) 45 48 return scale * value … … 78 81 self.smearer = smear_selection( self.data ) 79 82 80 81 83 def set_fit_range(self ,qmin=None, qmax=None): 82 84 """ to set the fit range""" … … 160 162 self._low_extrapolation_npts = 4 161 163 self._low_extrapolation_function = guinier 162 self._low_extrapolation_power = 4164 self._low_extrapolation_power = None 163 165 164 166 self._high_extrapolation_npts = 4 165 167 self._high_extrapolation_function = power_law 166 self._high_extrapolation_power = 4 167 168 169 def set_extrapolation(self, range, npts=4, function=None, power=4): 170 """ 171 Set the extrapolation parameters for the high or low Q-range. 172 Note that this does not turn extrapolation on or off. 173 @param range: a keyword set the type of extrapolation . type string 174 @param npts: the numbers of q points of data to consider for extrapolation 175 @param function: a keyword to select the function to use for extrapolation. 176 of type string. 177 @param power: an power to apply power_low function 178 179 """ 180 range = range.lower() 181 if range not in ['high', 'low']: 182 raise ValueError, "Extrapolation range should be 'high' or 'low'" 183 function = function.lower() 184 if function not in ['power_law', 'guinier']: 185 raise ValueError, "Extrapolation function should be 'guinier' or 'power_law'" 186 187 if range == 'high': 188 if function != 'power_law': 189 raise ValueError, "Extrapolation only allows a power law at high Q" 190 self._high_extrapolation_npts = npts 191 self._high_extrapolation_power = power 192 else: 193 if function == 'power_law': 194 self._low_extrapolation_function = power_law 195 else: 196 self._low_extrapolation_function = guinier 197 self._low_extrapolation_npts = npts 198 self._low_extrapolation_power = power 168 self._high_extrapolation_power = None 199 169 200 170 def _get_data(self, data): … … 212 182 return new_data 213 183 214 def _fit(self, function, qmin=Q_MINIMUM, qmax=Q_MAXIMUM ):184 def _fit(self, function, qmin=Q_MINIMUM, qmax=Q_MAXIMUM, power=None): 215 185 """ 216 186 fit data with function using … … 221 191 @param qmin: data first q value to consider during the fit 222 192 @param qmax: data last q value to consider during the fit 193 @param power : power value to consider for power-law 223 194 @param function: the function to use during the fit 224 195 @return a: the scale of the function … … 232 203 fit_y = numpy.array([math.log(y) for y in self._data.y]) 233 204 elif function.__name__ == "power_law": 234 fit_x = numpy.array([math.log(x) for x in self._data.x]) 235 qmin = math.log(qmin) 236 qmax = math.log(qmax) 205 if power is None: 206 fit_x = numpy.array([math.log(x) for x in self._data.x]) 207 qmin = math.log(qmin) 208 qmax = math.log(qmax) 209 fit_y = numpy.array([math.log(y) for y in self._data.y]) 210 else: 211 raise ValueError("Unknown function used to fit %s"%function.__name__) 212 213 214 if function.__name__ == "power_law" and power is None: 215 b = math.fabs(power) 216 fit_x = numpy.array([math.pow(x, -b) for x in self._data.x]) 237 217 fit_y = numpy.array([math.log(y) for y in self._data.y]) 238 else: 239 raise ValueError("Unknown function used to fit %s"%function.__name__) 240 241 fit_data = LoaderData1D(x=fit_x, y=fit_y) 242 fit_data.dxl = self._data.dxl 243 fit_data.dxw = self._data.dxw 244 245 functor = FitFunctor(data=fit_data) 246 functor.set_fit_range(qmin=qmin, qmax=qmax) 247 b, a = functor.fit() 218 a = (fit_y - fit_x)/(len(self._data.x)) 219 else: 220 fit_data = LoaderData1D(x=fit_x, y=fit_y) 221 fit_data.dxl = self._data.dxl 222 fit_data.dxw = self._data.dxw 223 functor = FitFunctor(data=fit_data) 224 functor.set_fit_range(qmin=qmin, qmax=qmax) 225 b, a = functor.fit() 248 226 249 227 if function.__name__ == "guinier": … … 252 230 if function.__name__ == "power_law": 253 231 b = -1 * b 232 if b <= 0: 233 raise ValueError("Power_law fit expected posive power, but got %s"%power) 254 234 # a is the scale of the guinier function 255 235 a = math.exp(a) … … 269 249 return self._get_qstar_unsmear(data) 270 250 271 272 def get_qstar(self, extrapolation=None):273 """274 Compute the invariant of the local copy of data.275 Implementation:276 if slit smear:277 qstar_0 = self._get_qstar_smear()278 else:279 qstar_0 = self._get_qstar_unsmear()280 if extrapolation is None:281 return qstar_0282 if extrapolation==low:283 return qstar_0 + self._get_qstar_low()284 elif extrapolation==high:285 return qstar_0 + self._get_qstar_high()286 elif extrapolation==both:287 return qstar_0 + self._get_qstar_low() + self._get_qstar_high()288 289 @param extrapolation: string to apply optional extrapolation290 @return q_star: invariant of the data within data's q range291 292 @warning: When using setting data to Data1D , the user is responsible of293 checking that the scale and the background are properly apply to the data294 295 @warning: if error occur self._get_qstar_low() or self._get_qstar_low()296 their returned value will be ignored297 """298 qstar_0 = self._get_qstar(data=self._data)299 300 if extrapolation is None:301 self._qstar = qstar_0302 return self._qstar303 # Compute invariant plus invaraint of extrapolated data304 extrapolation = extrapolation.lower()305 if extrapolation == "low":306 self._qstar = qstar_0 + self._get_qstar_low()307 return self._qstar308 elif extrapolation == "high":309 self._qstar = qstar_0 + self._get_qstar_high()310 return self._qstar311 elif extrapolation == "both":312 self._qstar = qstar_0 + self._get_qstar_low() + self._get_qstar_high()313 return self._qstar314 315 251 def _get_qstar_unsmear(self, data): 316 252 """ … … 518 454 sum += (data.x[i] * data.dxl[i] * dy[i] * dxi)**2 519 455 return math.sqrt(sum) 520 521 def get_surface(self,contrast, porod_const): 522 """ 523 Compute the surface of the data. 524 525 Implementation: 526 V= self.get_volume_fraction(contrast) 527 528 Compute the surface given by: 529 surface = (2*pi *V(1- V)*porod_const)/ q_star 530 531 @param contrast: contrast value to compute the volume 532 @param porod_const: Porod constant to compute the surface 533 @return: specific surface 534 """ 535 #Check whether we have Q star 536 if self._qstar is None: 537 self._qstar = self.get_star() 538 if self._qstar == 0: 539 raise RuntimeError("Cannot compute surface, invariant value is zero") 540 # Compute the volume 541 volume = self.get_volume_fraction(contrast) 542 return 2 * math.pi * volume *(1 - volume) * float(porod_const)/self._qstar 543 544 545 def get_volume_fraction(self, contrast): 546 """ 547 Compute volume fraction is deduced as follow: 548 549 q_star = 2*(pi*contrast)**2* volume( 1- volume) 550 for k = 10^(-8)*q_star/(2*(pi*|contrast|)**2) 551 we get 2 values of volume: 552 with 1 - 4 * k >= 0 553 volume1 = (1- sqrt(1- 4*k))/2 554 volume2 = (1+ sqrt(1- 4*k))/2 555 556 q_star: the invariant value included extrapolation is applied 557 unit 1/A^(3)*1/cm 558 q_star = self.get_qstar() 559 560 the result returned will be 0<= volume <= 1 561 562 @param contrast: contrast value provides by the user of type float. 563 contrast unit is 1/A^(2)= 10^(16)cm^(2) 564 @return: volume fraction 565 @note: volume fraction must have no unit 566 """ 567 if contrast < 0: 568 raise ValueError, "contrast must be greater than zero" 569 570 if self._qstar is None: 571 self._qstar = self.get_qstar() 572 573 if self._qstar < 0: 574 raise RuntimeError, "invariant must be greater than zero" 575 576 # Compute intermediate constant 577 k = 1.e-8 * self._qstar/(2 * (math.pi * math.fabs(float(contrast)))**2) 578 #Check discriminant value 579 discrim = 1 - 4 * k 580 581 # Compute volume fraction 582 if discrim < 0: 583 raise RuntimeError, "could not compute the volume fraction: negative discriminant" 584 elif discrim == 0: 585 return 1/2 586 else: 587 volume1 = 0.5 * (1 - math.sqrt(discrim)) 588 volume2 = 0.5 * (1 + math.sqrt(discrim)) 589 590 if 0 <= volume1 and volume1 <= 1: 591 return volume1 592 elif 0 <= volume2 and volume2 <= 1: 593 return volume2 594 raise RuntimeError, "could not compute the volume fraction: inconsistent results" 595 456 596 457 def _get_qstar_low(self): 597 458 """ … … 648 509 # fit the data with a model to get the appropriate parameters 649 510 a, b = self._fit(function=self._low_extrapolation_function, 650 qmin=qmin, qmax=qmax) 511 qmin=qmin, 512 qmax=qmax, 513 power=self._low_extrapolation_power) 651 514 except: 652 515 return None … … 701 564 # fit the data with a model to get the appropriate parameters 702 565 a, b = self._fit(function=self._high_extrapolation_function, 703 qmin=qmin, qmax=qmax) 566 qmin=qmin, 567 qmax=qmax, 568 power=self._high_extrapolation_power) 704 569 except: 705 570 return None … … 728 593 729 594 return data_max 595 596 def set_extrapolation(self, range, npts=4, function=None, power=None): 597 """ 598 Set the extrapolation parameters for the high or low Q-range. 599 Note that this does not turn extrapolation on or off. 600 @param range: a keyword set the type of extrapolation . type string 601 @param npts: the numbers of q points of data to consider for extrapolation 602 @param function: a keyword to select the function to use for extrapolation. 603 of type string. 604 @param power: an power to apply power_low function 605 606 """ 607 range = range.lower() 608 if range not in ['high', 'low']: 609 raise ValueError, "Extrapolation range should be 'high' or 'low'" 610 function = function.lower() 611 if function not in ['power_law', 'guinier']: 612 raise ValueError, "Extrapolation function should be 'guinier' or 'power_law'" 613 614 if range == 'high': 615 if function != 'power_law': 616 raise ValueError, "Extrapolation only allows a power law at high Q" 617 self._high_extrapolation_npts = npts 618 self._high_extrapolation_power = power 619 else: 620 if function == 'power_law': 621 self._low_extrapolation_function = power_law 622 else: 623 self._low_extrapolation_function = guinier 624 self._low_extrapolation_npts = npts 625 self._low_extrapolation_power = power 626 627 def get_qstar(self, extrapolation=None): 628 """ 629 Compute the invariant of the local copy of data. 630 Implementation: 631 if slit smear: 632 qstar_0 = self._get_qstar_smear() 633 else: 634 qstar_0 = self._get_qstar_unsmear() 635 if extrapolation is None: 636 return qstar_0 637 if extrapolation==low: 638 return qstar_0 + self._get_qstar_low() 639 elif extrapolation==high: 640 return qstar_0 + self._get_qstar_high() 641 elif extrapolation==both: 642 return qstar_0 + self._get_qstar_low() + self._get_qstar_high() 643 644 @param extrapolation: string to apply optional extrapolation 645 @return q_star: invariant of the data within data's q range 646 647 @warning: When using setting data to Data1D , the user is responsible of 648 checking that the scale and the background are properly apply to the data 649 650 @warning: if error occur self._get_qstar_low() or self._get_qstar_low() 651 their returned value will be ignored 652 """ 653 qstar_0 = self._get_qstar(data=self._data) 654 655 if extrapolation is None: 656 self._qstar = qstar_0 657 return self._qstar 658 # Compute invariant plus invaraint of extrapolated data 659 extrapolation = extrapolation.lower() 660 if extrapolation == "low": 661 self._qstar = qstar_0 + self._get_qstar_low() 662 return self._qstar 663 elif extrapolation == "high": 664 self._qstar = qstar_0 + self._get_qstar_high() 665 return self._qstar 666 elif extrapolation == "both": 667 self._qstar = qstar_0 + self._get_qstar_low() + self._get_qstar_high() 668 return self._qstar 669 670 def get_surface(self,contrast, porod_const): 671 """ 672 Compute the surface of the data. 673 674 Implementation: 675 V= self.get_volume_fraction(contrast) 676 677 Compute the surface given by: 678 surface = (2*pi *V(1- V)*porod_const)/ q_star 679 680 @param contrast: contrast value to compute the volume 681 @param porod_const: Porod constant to compute the surface 682 @return: specific surface 683 """ 684 #Check whether we have Q star 685 if self._qstar is None: 686 self._qstar = self.get_star() 687 if self._qstar == 0: 688 raise RuntimeError("Cannot compute surface, invariant value is zero") 689 # Compute the volume 690 volume = self.get_volume_fraction(contrast) 691 return 2 * math.pi * volume *(1 - volume) * float(porod_const)/self._qstar 692 693 def get_volume_fraction(self, contrast): 694 """ 695 Compute volume fraction is deduced as follow: 696 697 q_star = 2*(pi*contrast)**2* volume( 1- volume) 698 for k = 10^(-8)*q_star/(2*(pi*|contrast|)**2) 699 we get 2 values of volume: 700 with 1 - 4 * k >= 0 701 volume1 = (1- sqrt(1- 4*k))/2 702 volume2 = (1+ sqrt(1- 4*k))/2 703 704 q_star: the invariant value included extrapolation is applied 705 unit 1/A^(3)*1/cm 706 q_star = self.get_qstar() 707 708 the result returned will be 0<= volume <= 1 709 710 @param contrast: contrast value provides by the user of type float. 711 contrast unit is 1/A^(2)= 10^(16)cm^(2) 712 @return: volume fraction 713 @note: volume fraction must have no unit 714 """ 715 if contrast < 0: 716 raise ValueError, "contrast must be greater than zero" 717 718 if self._qstar is None: 719 self._qstar = self.get_qstar() 720 721 if self._qstar < 0: 722 raise RuntimeError, "invariant must be greater than zero" 723 724 # Compute intermediate constant 725 k = 1.e-8 * self._qstar/(2 * (math.pi * math.fabs(float(contrast)))**2) 726 #Check discriminant value 727 discrim = 1 - 4 * k 728 729 # Compute volume fraction 730 if discrim < 0: 731 raise RuntimeError, "could not compute the volume fraction: negative discriminant" 732 elif discrim == 0: 733 return 1/2 734 else: 735 volume1 = 0.5 * (1 - math.sqrt(discrim)) 736 volume2 = 0.5 * (1 + math.sqrt(discrim)) 737 738 if 0 <= volume1 and volume1 <= 1: 739 return volume1 740 elif 0 <= volume2 and volume2 <= 1: 741 return volume2 742 raise RuntimeError, "could not compute the volume fraction: inconsistent results" 730 743 731 744 def get_qstar_with_error(self, extrapolation=None): -
Invariant/test/utest_use_cases.py
r437a9f0 r4e80ae0 7 7 from DataLoader.loader import Loader 8 8 from sans.invariant import invariant 9 9 _ERR_SURFACE = 0.3 10 10 class Data1D: 11 11 print "I am not of type Dataloader.Data1D" … … 52 52 # Test results 53 53 self.assertAlmostEquals(qstar, 7.48959e-5,2) 54 self.assertAlmostEquals(v, 0.005644689 )55 self.assertAlmostEquals(s , 9.42e+2,2)54 self.assertAlmostEquals(v, 0.005644689, 4) 55 self.assertAlmostEquals(s + _ERR_SURFACE, 9.42e+2, 1) 56 56 57 57 def test_use_case_2(self): … … 73 73 # Test results 74 74 self.assertAlmostEquals(qstar, 7.48959e-5,2) 75 self.assertAlmostEquals(v, 0.005644689 )76 self.assertAlmostEquals(s , 9.42e+2,2)75 self.assertAlmostEquals(v, 0.005644689, 1) 76 self.assertAlmostEquals(s + _ERR_SURFACE, 9.42e+2, 1) 77 77 78 78 … … 108 108 109 109 # Test results 110 self.assertAlmostEquals(qstar, 7.49e-5 )111 self.assertAlmostEquals(v, 0.005648401 )112 self.assertAlmostEquals(s , 9.42e+2,2)110 self.assertAlmostEquals(qstar, 7.49e-5, 1) 111 self.assertAlmostEquals(v, 0.005648401, 4) 112 self.assertAlmostEquals(s + _ERR_SURFACE, 9.42e+2, 1) 113 113 114 114 def test_use_case_4(self): … … 135 135 # Test results 136 136 self.assertAlmostEquals(qstar, 7.49e-5,2) 137 self.assertAlmostEquals(v, 0.005952674 )138 self.assertAlmostEquals(s , 9.42e+2,2)137 self.assertAlmostEquals(v, 0.005952674, 3) 138 self.assertAlmostEquals(s + _ERR_SURFACE, 9.42e+2, 1) 139 139 140 140 def test_use_case_5(self): … … 162 162 # Test results 163 163 self.assertAlmostEquals(qstar, 7.88981e-5,2) 164 self.assertAlmostEquals(v, 0.005952674 )165 self.assertAlmostEquals(s , 9.42e+2,2)164 self.assertAlmostEquals(v, 0.005952674, 3) 165 self.assertAlmostEquals(s + _ERR_SURFACE, 9.42e+2, 1) 166 166 167 167 class TestInvSlitSmear(unittest.TestCase): … … 173 173 list = Loader().load("latex_smeared.xml") 174 174 self.data_slit_smear = list[1] 175 175 self.data_slit_smear.dxl = list[1].dxl 176 self.data_slit_smear.dxw = list[1].dxw 177 176 178 def test_use_case_1(self): 177 179 """ … … 185 187 s, ds = inv.get_surface_with_error(contrast=2.6e-6, porod_const=2) 186 188 # Test results 187 self.assertAlmostEquals(qstar, 4.1539e-4 )188 self.assertAlmostEquals(v, 0.032164596 )189 self.assertAlmostEquals(s , 9.42e+2,2)189 self.assertAlmostEquals(qstar, 4.1539e-4, 1) 190 self.assertAlmostEquals(v, 0.032164596, 3) 191 self.assertAlmostEquals(s + _ERR_SURFACE, 9.42e+2, 1) 190 192 191 193 def test_use_case_2(self): … … 202 204 s, ds = inv.get_surface_with_error(contrast=2.6e-6, porod_const=2) 203 205 # Test results 204 self.assertAlmostEquals(qstar, 4.1539e-4 )205 self.assertAlmostEquals(v, 0.032164596 )206 self.assertAlmostEquals(s , 9.42e+2,2)206 self.assertAlmostEquals(qstar, 4.1539e-4, 1) 207 self.assertAlmostEquals(v, 0.032164596,3) 208 self.assertAlmostEquals(s + _ERR_SURFACE, 9.42e+2, 1) 207 209 208 210 def test_use_case_3(self): … … 226 228 # Test results 227 229 self.assertAlmostEquals(qstar,4.1534e-4,3) 228 self.assertAlmostEquals(v, 0.032164596 )229 self.assertAlmostEquals(s , 9.42e+2,2)230 self.assertAlmostEquals(v, 0.032164596, 3) 231 self.assertAlmostEquals(s + _ERR_SURFACE, 9.42e+2, 1) 230 232 231 233 def test_use_case_4(self): … … 251 253 252 254 # Test results 253 self.assertAlmostEquals(qstar, 4.1539e-4 )254 self.assertAlmostEquals(v, 0.032164596 )255 self.assertAlmostEquals(s , 9.42e+2,2)255 self.assertAlmostEquals(qstar, 4.1539e-4, 2) 256 self.assertAlmostEquals(v, 0.032164596, 3) 257 self.assertAlmostEquals(s + _ERR_SURFACE, 9.42e+2, 1) 256 258 257 259 def test_use_case_5(self): … … 279 281 # Test results 280 282 self.assertAlmostEquals(qstar, 4.1534e-4,3) 281 self.assertAlmostEquals(v, 0.032164596) 282 self.assertAlmostEquals(s, 9.42e+2,2) 283 self.assertAlmostEquals(v, 0.032164596, 2) 284 self.assertAlmostEquals(s + _ERR_SURFACE, 9.42e+2, 1) 285 283 286 284 287 class TestInvPinholeSmear(unittest.TestCase): … … 290 293 list = Loader().load("latex_smeared.xml") 291 294 self.data_q_smear = list[0] 295 self.data_q_smear.dxl = list[0].dxl 296 self.data_q_smear.dxw = list[0].dxw 292 297 293 298 def test_use_case_1(self): … … 298 303 qstar = inv.get_qstar() 299 304 300 v olume_fraction= inv.get_volume_fraction(contrast=2.6e-6)301 s urface= inv.get_surface(contrast=2.6e-6, porod_const=2)302 # Test results 303 self.assertAlmostEquals(qstar, 1.361677e-3 )304 self.assertAlmostEquals(v, 0.115352622 )305 self.assertAlmostEquals(s , 9.42e+2,2)305 v = inv.get_volume_fraction(contrast=2.6e-6) 306 s = inv.get_surface(contrast=2.6e-6, porod_const=2) 307 # Test results 308 self.assertAlmostEquals(qstar, 1.361677e-3, 4) 309 self.assertAlmostEquals(v, 0.115352622, 2) 310 self.assertAlmostEquals(s + _ERR_SURFACE, 9.42e+2, 1 ) 306 311 307 312 def test_use_case_2(self): … … 319 324 s, ds = inv.get_surface_with_error(contrast=2.6e-6, porod_const=2) 320 325 # Test results 321 self.assertAlmostEquals(qstar, 1.361677e-3 )322 self.assertAlmostEquals(v, 0.115352622 )323 self.assertAlmostEquals(s , 9.42e+2,2)326 self.assertAlmostEquals(qstar, 1.361677e-3, 4) 327 self.assertAlmostEquals(v, 0.115352622, 2) 328 self.assertAlmostEquals(s + _ERR_SURFACE, 9.42e+2, 1 ) 324 329 325 330 def test_use_case_3(self): … … 342 347 self.assertAlmostEquals(qstar, 0.002037677) 343 348 self.assertAlmostEquals(v, 0.115352622) 344 self.assertAlmostEquals(s , 9.42e+2,2)349 self.assertAlmostEquals(s + _ERR_SURFACE, 9.42e+2, 1) 345 350 346 351 def test_use_case_4(self): … … 361 366 362 367 # Test results 363 self.assertAlmostEquals(qstar, 1.481677e-3 )368 self.assertAlmostEquals(qstar, 1.481677e-3, 3) 364 369 self.assertAlmostEquals(v, 0.127225804) 365 self.assertAlmostEquals(s , 9.42e+2,2)370 self.assertAlmostEquals(s+ _ERR_SURFACE, 9.42e+2, 1) 366 371 367 372 def test_use_case_5(self): … … 384 389 385 390 # Test results 386 self.assertAlmostEquals(qstar, 0.0021621 )391 self.assertAlmostEquals(qstar, 0.0021621, 3) 387 392 self.assertAlmostEquals(v, 0.202846825) 388 self.assertAlmostEquals(s , 9.42e+2,2)393 self.assertAlmostEquals(s+ _ERR_SURFACE, 9.42e+2, 1)
Note: See TracChangeset
for help on using the changeset viewer.