- Timestamp:
- Dec 17, 2009 12:19:56 PM (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:
- 472b11c
- Parents:
- 8d1e745
- Location:
- Invariant
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
Invariant/invariant.py
r4e80ae0 r9b6497bb 6 6 import math 7 7 import numpy 8 import scipy 8 9 9 10 from DataLoader.data_info import Data1D as LoaderData1D … … 20 21 INTEGRATION_NSTEPS = 1000 21 22 22 def guinier(x, scale=1, radius= 0.1):23 def guinier(x, scale=1, radius=60): 23 24 """ 24 25 Compute a F(x) = scale* e-((radius*x)**2/3). … … 28 29 @return F(x) 29 30 """ 31 if radius <= 0: 32 raise ValueError("Rg expected positive value, but got %s"%radius) 30 33 value = numpy.array([math.exp(-((radius * i)**2/3)) for i in x ]) 34 scale = numpy.sqrt(scale*scale) 35 if scale == 0: 36 raise ValueError("scale expected positive value, but got %s"%scale) 31 37 return scale * value 32 38 33 def power_law(x, scale=1, power= 0):39 def power_law(x, scale=1, power=4): 34 40 """ 35 41 F(x) = scale* (x)^(-power) … … 43 49 """ 44 50 if power <=0: 45 raise ValueError("Power_law function expected posi ve power, but got %s"%power)51 raise ValueError("Power_law function expected positive power, but got %s"%power) 46 52 47 53 value = numpy.array([ math.pow(i, -power) for i in x ]) 54 scale = numpy.sqrt(scale*scale) 55 if scale == 0: 56 raise ValueError("scale expected positive value, but got %s"%scale) 48 57 return scale * value 49 58 … … 111 120 112 121 """ 122 123 fx = numpy.zeros(len(self.data.x)) 124 sigma = numpy.zeros(len(self.data.x)) 125 126 #define dy^2 127 sigma = self.data.dy[self.idx_unsmeared ] 128 sigma2 = sigma*sigma 129 113 130 # Compute theory data f(x) 114 fx = numpy.zeros(len(self.data.x)) 115 fx = self.data.y[self.idx_unsmeared ] 131 fx = self.data.y[self.idx_unsmeared ]/sigma2 116 132 ## Smear theory data 117 133 if self.smearer is not None: 118 134 fx = self.smearer(fx, self._first_unsmeared_bin,self._last_unsmeared_bin) 119 135 120 A = numpy.vstack([ self.data.x[self.idx] ,121 numpy.ones(len(self.data.x[self.idx])) ]).T136 A = numpy.vstack([ self.data.x[self.idx]/sigma2, 137 numpy.ones(len(self.data.x[self.idx]))/sigma2]).T 122 138 123 139 a, b = numpy.linalg.lstsq(A, fx)[0] … … 178 194 raise ValueError,"Data must be of type DataLoader.Data1D" 179 195 new_data = self._scale * data - self._background 196 new_data.dy = data.dy 180 197 new_data.dxl = data.dxl 181 198 new_data.dxw = data.dxw … … 199 216 if function.__name__ == "guinier": 200 217 fit_x = numpy.array([x * x for x in self._data.x]) 218 201 219 qmin = qmin**2 202 220 qmax = qmax**2 203 221 fit_y = numpy.array([math.log(y) for y in self._data.y]) 222 fit_dy = numpy.array([y for y in self._data.y]) 223 fit_dy = numpy.array([dy for dy in self._data.dy])/fit_dy 224 204 225 elif function.__name__ == "power_law": 205 226 if power is None: 206 227 fit_x = numpy.array([math.log(x) for x in self._data.x]) 228 207 229 qmin = math.log(qmin) 208 230 qmax = math.log(qmax) 209 fit_y = numpy.array([math.log(y) for y in self._data.y]) 231 fit_y = numpy.array([math.log(y) for y in self._data.y]) 232 fit_dy = numpy.array([y for y in self._data.y]) 233 fit_dy = numpy.array([dy for dy in self._data.dy])/fit_dy 234 210 235 else: 211 236 raise ValueError("Unknown function used to fit %s"%function.__name__) 212 213 214 if function.__name__ == "power_law" and power is None: 237 238 if function.__name__ == "power_law" and power != None: 215 239 b = math.fabs(power) 216 fit_x = numpy.array([math.pow(x, -b) for x in self._data.x])217 240 fit_y = numpy.array([math.log(y) for y in self._data.y]) 218 a = (fit_y - fit_x)/(len(self._data.x)) 219 else: 220 fit_data = LoaderData1D(x=fit_x, y=fit_y) 241 fit_dy = numpy.array([y for y in self._data.y]) 242 fit_dy = numpy.array([dy for dy in self._data.dy])/fit_dy 243 sigma2 = fit_dy*fit_dy 244 a = scipy.sum(fit_y/sigma2) - scipy.sum(fit_x/sigma2*b)/scipy.sum(sigma2) 245 else: 246 fit_data = LoaderData1D(x=fit_x, y=fit_y, dy=fit_dy) 221 247 fit_data.dxl = self._data.dxl 222 248 fit_data.dxw = self._data.dxw … … 225 251 b, a = functor.fit() 226 252 253 227 254 if function.__name__ == "guinier": 228 255 # b is the radius value of the guinier function 229 b = math.sqrt(-3 * b) 256 if b>=0: 257 raise ValueError("Guinier fit was not converged") 258 else: 259 b = math.sqrt(-3 * b) 260 261 230 262 if function.__name__ == "power_law": 231 b = -1 * b 263 if power == None: 264 b = -1 * b 232 265 if b <= 0: 233 266 raise ValueError("Power_law fit expected posive power, but got %s"%power) 234 267 # a is the scale of the guinier function 235 268 a = math.exp(a) 269 236 270 return a, b 237 271 … … 502 536 @return: a new data of type Data1D 503 537 """ 538 504 539 # Data boundaries for fiiting 505 540 qmin = self._data.x[0] … … 514 549 except: 515 550 return None 516 551 517 552 #q_start point 518 553 q_start = Q_MINIMUM … … 534 569 dxw = numpy.ones(INTEGRATION_NSTEPS) 535 570 dxw = dxw * self._data.dxw[0] 536 571 537 572 data_min = LoaderData1D(x=new_x, y=new_y) 538 573 data_min.dxl = dxl 539 574 data_min.dxw = dxw 540 575 self._data.clone_without_data( clone= data_min) 541 576 542 577 return data_min 543 578 -
Invariant/test/utest_use_cases.py
r4e80ae0 r9b6497bb 345 345 346 346 # Test results 347 self.assertAlmostEquals(qstar, 0.00 2037677)347 self.assertAlmostEquals(qstar, 0.00138756,2) 348 348 self.assertAlmostEquals(v, 0.115352622) 349 349 self.assertAlmostEquals(s + _ERR_SURFACE, 9.42e+2, 1) … … 389 389 390 390 # Test results 391 self.assertAlmostEquals(qstar, 0.00 21621, 3)391 self.assertAlmostEquals(qstar, 0.00460319, 3) 392 392 self.assertAlmostEquals(v, 0.202846825) 393 393 self.assertAlmostEquals(s+ _ERR_SURFACE, 9.42e+2, 1)
Note: See TracChangeset
for help on using the changeset viewer.