Changeset aafa962 in sasview
- Timestamp:
- Jan 16, 2010 12:44:34 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:
- 52070a1
- Parents:
- 82703a1
- Location:
- Invariant
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
Invariant/invariant.py
r82703a1 raafa962 3 3 @author: Gervaise B. Alina/UTK 4 4 """ 5 5 #TODO: Need to decide if/how to use smearing for extrapolation 6 6 import math 7 7 import numpy … … 24 24 function given some x, y 25 25 """ 26 def transform_values(self, a, b): 26 def linearize_data(self, x, y, dy=None): 27 """ 28 Linearize data so that a linear fit can be performed. 29 Filter out the data that can't be transformed. 30 @param x: array of q values 31 @param y: array of I(q) values 32 @param dy: array of dI(q) values 33 """ 34 return NotImplemented 35 36 def linearize_q_value(self, value): 37 """ 38 Transform the input q-value for linearization 39 """ 40 return NotImplemented 41 42 def extract_model_parameters(self, a, b): 27 43 """ 28 44 set private member 29 """30 return NotImplemented31 32 def transform_value(self, value):33 """34 @param value: float35 return f(value)36 45 """ 37 46 return NotImplemented 38 47 39 def inverse_transform_value(self, value): 40 """ 41 reverse transform vaalue given a specific function 42 return f-1(value) 43 @param value : float 48 def evaluate_model(self, x): 49 """ 50 Returns an array f(x) values where f is the Transform function. 44 51 """ 45 52 return NotImplemented 46 53 47 def get_extrapolated_data(self, data, q_start=Q_MINIMUM, q_end=Q_MAXIMUM):48 """49 @return extrapolate data create from data50 """51 return NotImplemented52 53 def transform_data(self, x):54 """55 @param x: vector of float56 @param y : vector of float57 return x, y transform given a specific function58 """59 return NotImplemented60 61 def inverse_transform_data(self, x, y):62 """63 reverse transform x, y given a specific function64 """65 return NotImplemented66 67 def transform_error(self, y , dy=None):68 if dy is None:69 dy = numpy.array([math.sqrt(j) for j in y])70 return dy / y71 72 def get_extrapolated_data(self, data, npts=INTEGRATION_NSTEPS,73 q_start=Q_MINIMUM, q_end=Q_MAXIMUM, smear_indice=0):74 """75 @return extrapolate data create from data76 """77 #create new Data1D to compute the invariant78 new_x = numpy.linspace(start=q_start,79 stop=q_end,80 num=npts,81 endpoint=True)82 new_y = self.transform_data(x=new_x)83 84 dxl = None85 dxw = None86 if data.dxl is not None :87 dxl = numpy.ones(INTEGRATION_NSTEPS)88 dxl = dxl * data.dxl[smear_indice]89 90 if data.dxw is not None :91 dxw = numpy.ones(INTEGRATION_NSTEPS)92 dxw = dxw * data.dxw[smear_indice]93 94 result_data = LoaderData1D(x=new_x, y=new_y)95 data.clone_without_data(clone=result_data)96 result_data.dxl = dxl97 result_data.dxw = dxw98 return result_data99 100 54 class Guinier(Transform): 101 55 """ … … 108 62 self.radius = radius 109 63 110 def transform_value(self, value): 111 """ 112 Square input value 113 @param value: of type float 64 def linearize_data(self, x, y, dy=None): 65 """ 66 Linearize data so that a linear fit can be performed. 67 Filter out the data that can't be transformed. 68 @param x: array of q values 69 @param y: array of I(q) values 70 @param dy: array of dI(q) values 71 """ 72 # Check that the vector lengths are equal 73 assert(len(x)==len(y)) 74 if dy is not None: 75 assert(len(x)==len(dy)) 76 else: 77 dy = numpy.array([math.sqrt(math.fabs(j)) for j in y]) 78 79 data_points = zip(x,y,dy) 80 81 # Transform the data 82 output_points = [(self.linearize_q_value(p[0]), 83 math.log(p[1]), 84 p[2]/p[1]) for p in data_points if p[0]>0 and p[1]>0] 85 86 x_out, y_out, dy_out = zip(*output_points) 87 88 return numpy.asarray(x_out), numpy.asarray(y_out), numpy.asarray(dy_out) 89 90 def linearize_q_value(self, value): 91 """ 92 Transform the input q-value for linearization 93 @param value: q-value 94 @return: q*q 114 95 """ 115 96 return value * value 116 97 117 def transform_values(self, a, b):98 def extract_model_parameters(self, a, b): 118 99 """ 119 100 assign new value to the scale and the radius … … 125 106 return a, b 126 107 127 def transform_data(self, x):108 def evaluate_model(self, x): 128 109 """ 129 110 return F(x)= scale* e-((radius*x)**2/3) 130 111 """ 131 112 return self._guinier(x) 132 133 def inverse_transform_data(self, x, y): 134 """ 135 this function finds x, y given equation1: y = ax + b 136 and equation2: y1 = scale* e-((radius*x1)**2/3). 137 where equation1, equation2 are equivalent 138 imply: x = x1*x1 139 y = ln(y1) 140 b = math.exp(scale) 141 a = -radius**2/3 142 @return x, y 143 """ 144 result_x = numpy.array([i * i for i in x]) 145 result_y = numpy.array([math.log(j) for j in y]) 146 return result_x, result_y 147 113 148 114 def _guinier(self, x): 149 115 """ … … 173 139 self.power = power 174 140 175 def transform_values(self, a, b): 141 def linearize_data(self, x, y, dy=None): 142 """ 143 Linearize data so that a linear fit can be performed. 144 Filter out the data that can't be transformed. 145 @param x: array of q values 146 @param y: array of I(q) values 147 @param dy: array of dI(q) values 148 """ 149 # Check that the vector lengths are equal 150 assert(len(x)==len(y)) 151 if dy is not None: 152 assert(len(x)==len(dy)) 153 else: 154 dy = numpy.array([math.sqrt(math.fabs(j)) for j in y]) 155 156 data_points = zip(x,y,dy) 157 158 # Transform the data 159 output_points = [(self.linearize_q_value(p[0]), 160 math.log(p[1]), 161 p[2]/p[1]) for p in data_points if p[0]>0 and p[1]>0] 162 163 x_out, y_out, dy_out = zip(*output_points) 164 165 return numpy.asarray(x_out), numpy.asarray(y_out), numpy.asarray(dy_out) 166 167 def linearize_q_value(self, value): 168 """ 169 Transform the input q-value for linearization 170 @param value: q-value 171 @return: log(q) 172 """ 173 return math.log(value) 174 175 def extract_model_parameters(self, a, b): 176 176 """ 177 177 Assign new value to the scale and the power … … 183 183 return a, b 184 184 185 def transform_value(self, value): 186 """ 187 compute log(value) 188 """ 189 return math.log(value) 190 191 def transform_data(self, x): 185 def evaluate_model(self, x): 192 186 """ 193 187 given a scale and a radius transform x, y using a power_law … … 196 190 return self._power_law(x) 197 191 198 def inverse_transform_data(self, x, y):199 """200 given a scale and a radius transform x, y using a inverse power_law201 function202 """203 result_x = numpy.array([math.log(i) for i in x])204 result_y = numpy.array([math.log(j) for j in y])205 206 return result_x, result_y207 208 192 def _power_law(self, x): 209 193 """ … … 227 211 class Extrapolator: 228 212 """ 229 compute f(x)213 Extrapolate I(q) distribution using a given model 230 214 """ 231 215 def __init__(self, data): … … 249 233 self.set_fit_range() 250 234 251 def set_fit_range(self ,qmin=None, qmax=None):235 def set_fit_range(self, qmin=None, qmax=None): 252 236 """ to set the fit range""" 253 237 if qmin is not None: … … 305 289 306 290 return a, b 291 307 292 308 293 class InvariantCalculator(object): … … 364 349 return new_data 365 350 366 def _fit(self, function, qmin=Q_MINIMUM, qmax=Q_MAXIMUM, power=None):351 def _fit(self, model, qmin=Q_MINIMUM, qmax=Q_MAXIMUM, power=None): 367 352 """ 368 353 fit data with function using … … 379 364 for power_law will be the power value 380 365 """ 381 transform = function 382 fit_x, fit_y = transform.inverse_transform_data(self._data.x, self._data.y) 383 fit_dy = transform.transform_error(y=self._data.y, dy=self._data.dy) 384 qmin = transform.transform_value(value=qmin) 385 qmax = transform.transform_value(value=qmax) 386 366 # Linearize the data in preparation for fitting 367 fit_x, fit_y, fit_dy = model.linearize_data(self._data.x, self._data.y, self._data.dy) 368 369 qmin = model.linearize_q_value(qmin) 370 qmax = model.linearize_q_value(qmax) 371 372 # Create a new Data1D object for processing 387 373 fit_data = LoaderData1D(x=fit_x, y=fit_y, dy=fit_dy) 388 374 fit_data.dxl = self._data.dxl … … 392 378 b, a = extrapolator.fit(power=power) 393 379 394 return function.transform_values(a=a, b=b)380 return model.extract_model_parameters(a=a, b=b) 395 381 396 382 def _get_qstar(self, data): … … 633 619 return self._get_qstar(data=data) 634 620 621 def _get_extrapolated_data(self, model, npts=INTEGRATION_NSTEPS, 622 q_start=Q_MINIMUM, q_end=Q_MAXIMUM): 623 """ 624 @return extrapolate data create from data 625 """ 626 #create new Data1D to compute the invariant 627 q = numpy.linspace(start=q_start, 628 stop=q_end, 629 num=npts, 630 endpoint=True) 631 iq = model.evaluate_model(q) 632 633 # Determine whether we are extrapolating to high or low q-values 634 # If the data is slit smeared, get the index of the slit dimension array entry 635 # that we will use to smear the extrapolated data. 636 dxl = None 637 dxw = None 638 639 if self._data.is_slit_smeared(): 640 if q_start<min(self._data.x): 641 smear_index = 0 642 elif q_end>max(self._data.x): 643 smear_index = len(self._data.x)-1 644 else: 645 raise RuntimeError, "Extrapolation can only be evaluated for points outside the data Q range" 646 647 if self._data.dxl is not None : 648 dxl = numpy.ones(INTEGRATION_NSTEPS) 649 dxl = dxl * self._data.dxl[smear_index] 650 651 if self._data.dxw is not None : 652 dxw = numpy.ones(INTEGRATION_NSTEPS) 653 dxw = dxw * self._data.dxw[smear_index] 654 655 result_data = LoaderData1D(x=q, y=iq) 656 result_data.dxl = dxl 657 result_data.dxw = dxw 658 return result_data 659 635 660 def get_extra_data_low(self): 636 661 """ … … 655 680 """ 656 681 657 # Data boundaries for fi iting682 # Data boundaries for fitting 658 683 qmin = self._data.x[0] 659 684 qmax = self._data.x[self._low_extrapolation_npts - 1] 685 660 686 # Extrapolate the low-Q data 661 #TODO: this fit fails. Fix it. 662 a, b = self._fit(function=self._low_extrapolation_function, 663 qmin=qmin, 664 qmax=qmax, 665 power=self._low_extrapolation_power) 687 a, b = self._fit(model=self._low_extrapolation_function, 688 qmin=qmin, 689 qmax=qmax, 690 power=self._low_extrapolation_power) 666 691 #q_start point 667 692 q_start = Q_MINIMUM … … 669 694 q_start = qmin/10 670 695 671 data_min = self._ low_extrapolation_function.get_extrapolated_data(data=self._data,672 npts=INTEGRATION_NSTEPS,673 q_start=q_start, q_end=qmin, smear_indice=0)696 data_min = self._get_extrapolated_data(model=self._low_extrapolation_function, 697 npts=INTEGRATION_NSTEPS, 698 q_start=q_start, q_end=qmin) 674 699 return data_min 675 700 … … 693 718 qmax = self._data.x[x_len] 694 719 q_end = Q_MAXIMUM 695 smear_indice = 0696 if self._data.dxl is not None:697 smear_indice = len(self._data.dxl) - 1698 720 699 721 # fit the data with a model to get the appropriate parameters 700 a, b = self._fit(function=self._high_extrapolation_function, 701 qmin=qmin, 702 qmax=qmax, 703 power=self._high_extrapolation_power) 722 a, b = self._fit(model=self._high_extrapolation_function, 723 qmin=qmin, 724 qmax=qmax, 725 power=self._high_extrapolation_power) 726 704 727 #create new Data1D to compute the invariant 705 data_max = self._high_extrapolation_function.get_extrapolated_data(data=self._data, 706 npts=INTEGRATION_NSTEPS, 707 q_start=qmax, q_end=q_end, 708 smear_indice=smear_indice) 728 data_max = self._get_extrapolated_data(model=self._high_extrapolation_function, 729 npts=INTEGRATION_NSTEPS, 730 q_start=qmax, q_end=q_end) 709 731 return data_max 710 732 -
Invariant/test/PolySpheres.txt
ref9ed58 raafa962 1 0 0 0 1 2 0.003 5.58459588740281 1.18158747955905 2 3 0.0031098987853131 5.58031083426704 0.590567038651574 -
Invariant/test/utest_data_handling.py
r6939bd4 raafa962 31 31 32 32 # Create invariant object. Background and scale left as defaults. 33 fit = invariant. FitFunctor(data=self.data)33 fit = invariant.Extrapolator(data=self.data) 34 34 a,b = fit.fit() 35 35 … … 48 48 49 49 # Create invariant object. Background and scale left as defaults. 50 fit = invariant. FitFunctor(data=self.data)50 fit = invariant.Extrapolator(data=self.data) 51 51 a,b = fit.fit() 52 52 … … 113 113 """ 114 114 self.scale = 1.5 115 self.rg = 70.0115 self.rg = 30.0 116 116 x = numpy.arange(0.0001, 0.1, 0.0001) 117 117 y = numpy.asarray([self.scale * math.exp( -(q*self.rg)**2 / 3.0 ) for q in x]) … … 129 129 130 130 self.assertEqual(inv._low_extrapolation_npts, 20) 131 self.assertEqual(inv._low_extrapolation_function.__ name__, 'guinier')131 self.assertEqual(inv._low_extrapolation_function.__class__, invariant.Guinier) 132 132 133 133 # Data boundaries for fiiting … … 136 136 137 137 # Extrapolate the low-Q data 138 a, b = inv._fit( function=inv._low_extrapolation_function,138 a, b = inv._fit(model=inv._low_extrapolation_function, 139 139 qmin=qmin, 140 140 qmax=qmax, … … 172 172 173 173 self.assertEqual(inv._low_extrapolation_npts, 20) 174 self.assertEqual(inv._low_extrapolation_function.__ name__, 'power_law')174 self.assertEqual(inv._low_extrapolation_function.__class__, invariant.PowerLaw) 175 175 176 176 # Data boundaries for fitting … … 179 179 180 180 # Extrapolate the low-Q data 181 a, b = inv._fit( function=inv._low_extrapolation_function,181 a, b = inv._fit(model=inv._low_extrapolation_function, 182 182 qmin=qmin, 183 183 qmax=qmax, … … 186 186 self.assertAlmostEqual(self.scale, a, 6) 187 187 self.assertAlmostEqual(self.m, b, 6) 188 188 189 class TestLinearization(unittest.TestCase): 190 191 def test_guinier_incompatible_length(self): 192 g = invariant.Guinier() 193 self.assertRaises(AssertionError, g.linearize_data, [1],[1,2],None) 194 self.assertRaises(AssertionError, g.linearize_data, [1,2],[1,2],[1]) 195 196 def test_linearization(self): 197 """ 198 Check that the linearization process filters out points 199 that can't be transformed 200 """ 201 x = numpy.asarray(numpy.asarray([0,1,2,3])) 202 y = numpy.asarray(numpy.asarray([1,1,1,1])) 203 g = invariant.Guinier() 204 x_out, y_out, dy_out = g.linearize_data(x,y,None) 205 self.assertEqual(len(x_out), 3) 206 self.assertEqual(len(y_out), 3) 207 self.assertEqual(len(dy_out), 3) 208 -
Invariant/test/utest_use_cases.py
r46d50ca raafa962 3 3 4 4 """ 5 #TODO: there's no test for smeared extrapolation 5 6 import unittest 6 7 import numpy … … 24 25 25 26 # Create invariant object. Background and scale left as defaults. 26 fit = invariant. FitFunctor(data=self.data)27 fit = invariant.Extrapolator(data=self.data) 27 28 28 29 ##Without holding … … 40 41 41 42 # Create invariant object. Background and scale left as defaults. 42 fit = invariant. FitFunctor(data=self.data)43 fit = invariant.Extrapolator(data=self.data) 43 44 44 45 #With holding a = -power =4 … … 62 63 63 64 # Create invariant object. Background and scale left as defaults. 64 fit = invariant. FitFunctor(data=self.data)65 fit = invariant.Extrapolator(data=self.data) 65 66 66 67 ##Without holding … … 78 79 79 80 # Create invariant object. Background and scale left as defaults. 80 fit = invariant. FitFunctor(data=self.data)81 fit = invariant.Extrapolator(data=self.data) 81 82 82 83 #With holding a = -power =4
Note: See TracChangeset
for help on using the changeset viewer.