Changeset 346bc88 in sasmodels for sasmodels/bumps_model.py
- Timestamp:
- Aug 31, 2015 10:24:45 PM (9 years ago)
- Branches:
- master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
- Children:
- 3e6aaad
- Parents:
- f224873
- git-author:
- Paul Kienzle <pkienzle@…> (08/31/15 22:12:39)
- git-committer:
- Paul Kienzle <pkienzle@…> (08/31/15 22:24:45)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
sasmodels/bumps_model.py
r0656250 r346bc88 1 1 """ 2 2 Wrap sasmodels for direct use by bumps. 3 4 :class:`Model` is a wrapper for the sasmodels kernel which defines a 5 bumps *Parameter* box for each kernel parameter. *Model* accepts keyword 6 arguments to set the initial value for each parameter. 7 8 :class:`Experiment` combines the *Model* function with a data file loaded by the 9 sasview data loader. *Experiment* takes a *cutoff* parameter controlling 10 how far the polydispersity integral extends. 11 12 A variety of helper functions are provided: 13 14 :func:`load_data` loads a sasview data file. 15 16 :func:`empty_data1D` creates an empty dataset, which is useful for plotting 17 a theory function before the data is measured. 18 19 :func:`empty_data2D` creates an empty 2D dataset. 20 21 :func:`set_beam_stop` masks the beam stop from the data. 22 23 :func:`set_half` selects the right or left half of the data, which can 24 be useful for shear measurements which have not been properly corrected 25 for path length and reflections. 26 27 :func:`set_top` cuts the top part off the data. 28 29 :func:`plot_data` plots the data file. 30 31 :func:`plot_theory` plots a calculated result from the model. 32 3 33 """ 4 34 5 35 import datetime 36 import warnings 6 37 7 38 import numpy as np 8 39 9 40 from . import sesans 41 from .resolution import Perfect1D, Pinhole1D, Slit1D 42 from .resolution2d import Pinhole2D 10 43 11 44 # CRUFT python 2.6 … … 20 53 21 54 55 # CRUFT: old style bumps wrapper which doesn't separate data and model 56 def BumpsModel(data, model, cutoff=1e-5, **kw): 57 warnings.warn("Use of BumpsModel is deprecated. Use bumps_model.Experiment instead.") 58 model = Model(model, **kw) 59 experiment = Experiment(data=data, model=model, cutoff=cutoff) 60 for k in model._parameter_names: 61 setattr(experiment, k, getattr(model, k)) 62 return experiment 63 64 65 22 66 def tic(): 23 67 """ … … 42 86 return data 43 87 44 45 def empty_data1D(q): 88 def plot_data(data, view='log'): 89 """ 90 Plot data loaded by the sasview loader. 91 """ 92 if hasattr(data, 'qx_data'): 93 _plot_2d_signal(data, data.data, view=view) 94 else: 95 # Note: kind of weird using the _plot_result1D to plot just the 96 # data, but it handles the masking and graph markup already, so 97 # do not repeat. 98 _plot_result1D(data, None, None, view) 99 100 def plot_theory(data, theory, view='log'): 101 if hasattr(data, 'qx_data'): 102 _plot_2d_signal(data, theory, view=view) 103 else: 104 _plot_result1D(data, theory, None, view, include_data=False) 105 106 107 def empty_data1D(q, resolution=0.05): 46 108 """ 47 109 Create empty 1D data using the given *q* as the x value. 48 110 49 Resolutions dq/q is5%.111 *resolution* dq/q defaults to 5%. 50 112 """ 51 113 … … 54 116 Iq = 100 * np.ones_like(q) 55 117 dIq = np.sqrt(Iq) 56 data = Data1D(q, Iq, dx= 0.05* q, dy=dIq)118 data = Data1D(q, Iq, dx=resolution * q, dy=dIq) 57 119 data.filename = "fake data" 58 120 data.qmin, data.qmax = q.min(), q.max() 121 data.mask = np.zeros(len(Iq), dtype='bool') 59 122 return data 60 123 61 124 62 def empty_data2D(qx, qy=None ):125 def empty_data2D(qx, qy=None, resolution=0.05): 63 126 """ 64 127 Create empty 2D data using the given mesh. … … 66 129 If *qy* is missing, create a square mesh with *qy=qx*. 67 130 68 Resolution dq/q is5%.131 *resolution* dq/q defaults to 5%. 69 132 """ 70 133 from sas.dataloader.data_info import Data2D, Detector … … 85 148 data.err_data = dIq 86 149 data.mask = mask 150 data.qmin = 1e-16 151 data.qmax = np.inf 87 152 88 153 # 5% dQ/Q resolution 89 data.dqx_data = 0.05 * Qx 90 data.dqy_data = 0.05 * Qy 154 if resolution != 0: 155 data.dqx_data = resolution * Qx 156 data.dqy_data = resolution * Qy 91 157 92 158 detector = Detector() … … 145 211 146 212 147 def plot_data(data, Iq, vmin=None, vmax=None, view='log'): 148 """ 149 Plot the target value for the data. This could be the data itself, 150 the theory calculation, or the residuals. 151 152 *scale* can be 'log' for log scale data, or 'linear'. 153 """ 154 from numpy.ma import masked_array 155 import matplotlib.pyplot as plt 156 if hasattr(data, 'qx_data'): 157 Iq = Iq + 0 158 valid = np.isfinite(Iq) 159 if view == 'log': 160 valid[valid] = (Iq[valid] > 0) 161 Iq[valid] = np.log10(Iq[valid]) 162 elif view == 'q4': 163 Iq[valid] = Iq*(data.qx_data[valid]**2+data.qy_data[valid]**2)**2 164 Iq[~valid | data.mask] = 0 165 #plottable = Iq 166 plottable = masked_array(Iq, ~valid | data.mask) 167 xmin, xmax = min(data.qx_data), max(data.qx_data) 168 ymin, ymax = min(data.qy_data), max(data.qy_data) 169 try: 170 if vmin is None: vmin = Iq[valid & ~data.mask].min() 171 if vmax is None: vmax = Iq[valid & ~data.mask].max() 172 except: 173 vmin, vmax = 0, 1 174 plt.imshow(plottable.reshape(128, 128), 175 interpolation='nearest', aspect=1, origin='upper', 176 extent=[xmin, xmax, ymin, ymax], vmin=vmin, vmax=vmax) 177 else: # 1D data 178 if view == 'linear' or view == 'q4': 179 #idx = np.isfinite(Iq) 180 scale = data.x**4 if view == 'q4' else 1.0 181 plt.plot(data.x, scale*Iq) #, '.') 182 else: 183 # Find the values that are finite and positive 184 idx = np.isfinite(Iq) 185 idx[idx] = (Iq[idx] > 0) 186 Iq[~idx] = np.nan 187 plt.loglog(data.x, Iq) 188 189 190 def _plot_result1D(data, theory, view): 213 def _plot_result1D(data, theory, resid, view, include_data=True): 191 214 """ 192 215 Plot the data and residuals for 1D data. … … 197 220 #data.y[data.y<0.05] = 0.5 198 221 mdata = masked_array(data.y, data.mask) 199 mdata[ np.isnan(mdata)] = masked222 mdata[~np.isfinite(mdata)] = masked 200 223 if view is 'log': 201 224 mdata[mdata <= 0] = masked 202 mtheory = masked_array(theory, mdata.mask)203 mresid = masked_array((theory - data.y) / data.dy, mdata.mask)204 225 205 226 scale = data.x**4 if view == 'q4' else 1.0 206 plt.subplot(121) 207 plt.errorbar(data.x, scale*mdata, yerr=data.dy) 208 plt.plot(data.x, scale*mtheory, '-', hold=True) 227 if resid is not None: 228 plt.subplot(121) 229 230 if include_data: 231 plt.errorbar(data.x, scale*mdata, yerr=data.dy, fmt='.') 232 if theory is not None: 233 mtheory = masked_array(theory, mdata.mask) 234 plt.plot(data.x, scale*mtheory, '-', hold=True) 235 plt.xscale(view) 209 236 plt.yscale('linear' if view == 'q4' else view) 210 plt.subplot(122) 211 plt.plot(data.x, mresid, 'x') 237 plt.xlabel('Q') 238 plt.ylabel('I(Q)') 239 if resid is not None: 240 mresid = masked_array(resid, mdata.mask) 241 plt.subplot(122) 242 plt.plot(data.x, mresid, 'x') 243 plt.ylabel('residuals') 244 plt.xlabel('Q') 245 plt.xscale(view) 212 246 213 247 # pylint: disable=unused-argument 214 def _plot_sesans(data, theory, view):248 def _plot_sesans(data, theory, resid, view): 215 249 import matplotlib.pyplot as plt 216 resid = (theory - data.y) / data.dy217 250 plt.subplot(121) 218 251 plt.errorbar(data.x, data.y, yerr=data.dy) … … 225 258 plt.ylabel('residuals (P/P0)') 226 259 227 def _plot_result2D(data, theory, view):260 def _plot_result2D(data, theory, resid, view): 228 261 """ 229 262 Plot the data and residuals for 2D data. 230 263 """ 231 264 import matplotlib.pyplot as plt 232 resid = (theory - data.data) / data.err_data 265 target = data.data[~data.mask] 266 if view == 'log': 267 vmin = min(target[target>0].min(), theory[theory>0].min()) 268 vmax = max(target.max(), theory.max()) 269 else: 270 vmin = min(target.min(), theory.min()) 271 vmax = max(target.max(), theory.max()) 272 #print vmin, vmax 233 273 plt.subplot(131) 234 plot_data(data, data.data, view=view) 274 _plot_2d_signal(data, target, view=view, vmin=vmin, vmax=vmax) 275 plt.title('data') 235 276 plt.colorbar() 236 277 plt.subplot(132) 237 plot_data(data, theory, view=view) 278 _plot_2d_signal(data, theory, view=view, vmin=vmin, vmax=vmax) 279 plt.title('theory') 238 280 plt.colorbar() 239 281 plt.subplot(133) 240 plot_data(data, resid, view='linear') 282 _plot_2d_signal(data, resid, view='linear') 283 plt.title('residuals') 241 284 plt.colorbar() 242 285 243 class BumpsModel(object): 244 """ 245 Return a bumps wrapper for a SAS model. 246 247 *data* is the data to be fitted. 248 249 *model* is the SAS model from :func:`core.load_model`. 250 251 *cutoff* is the integration cutoff, which avoids computing the 252 the SAS model where the polydispersity weight is low. 253 254 Model parameters can be initialized with additional keyword 255 arguments, or by assigning to model.parameter_name.value. 256 257 The resulting bumps model can be used directly in a FitProblem call. 258 """ 259 def __init__(self, data, model, cutoff=1e-5, **kw): 286 def _plot_2d_signal(data, signal, vmin=None, vmax=None, view='log'): 287 """ 288 Plot the target value for the data. This could be the data itself, 289 the theory calculation, or the residuals. 290 291 *scale* can be 'log' for log scale data, or 'linear'. 292 """ 293 import matplotlib.pyplot as plt 294 from numpy.ma import masked_array 295 296 image = np.zeros_like(data.qx_data) 297 image[~data.mask] = signal 298 valid = np.isfinite(image) 299 if view == 'log': 300 valid[valid] = (image[valid] > 0) 301 image[valid] = np.log10(image[valid]) 302 elif view == 'q4': 303 image[valid] *= (data.qx_data[valid]**2+data.qy_data[valid]**2)**2 304 image[~valid | data.mask] = 0 305 #plottable = Iq 306 plottable = masked_array(image, ~valid | data.mask) 307 xmin, xmax = min(data.qx_data), max(data.qx_data) 308 ymin, ymax = min(data.qy_data), max(data.qy_data) 309 # TODO: fix vmin, vmax so it is shared for theory/resid 310 vmin = vmax = None 311 try: 312 if vmin is None: vmin = image[valid & ~data.mask].min() 313 if vmax is None: vmax = image[valid & ~data.mask].max() 314 except: 315 vmin, vmax = 0, 1 316 #print vmin,vmax 317 plt.imshow(plottable.reshape(128, 128), 318 interpolation='nearest', aspect=1, origin='upper', 319 extent=[xmin, xmax, ymin, ymax], vmin=vmin, vmax=vmax) 320 321 322 class Model(object): 323 def __init__(self, kernel, **kw): 260 324 from bumps.names import Parameter 261 325 262 # remember inputs so we can inspect from outside 263 self.data = data 264 self.model = model 265 self.cutoff = cutoff 266 if hasattr(data, 'lam'): 267 self.data_type = 'sesans' 268 elif hasattr(data, 'qx_data'): 269 self.data_type = 'Iqxy' 270 else: 271 self.data_type = 'Iq' 272 273 partype = model.info['partype'] 274 275 # interpret data 276 if self.data_type == 'sesans': 277 q = sesans.make_q(data.sample.zacceptance, data.Rmax) 278 self.index = slice(None, None) 279 self.Iq = data.y 280 self.dIq = data.dy 281 self._theory = np.zeros_like(q) 282 q_vectors = [q] 283 elif self.data_type == 'Iqxy': 284 self.index = (data.mask == 0) & (~np.isnan(data.data)) 285 self.Iq = data.data[self.index] 286 self.dIq = data.err_data[self.index] 287 self._theory = np.zeros_like(data.data) 288 if not partype['orientation'] and not partype['magnetic']: 289 q_vectors = [np.sqrt(data.qx_data ** 2 + data.qy_data ** 2)] 290 else: 291 q_vectors = [data.qx_data, data.qy_data] 292 elif self.data_type == 'Iq': 293 self.index = (data.x >= data.qmin) & (data.x <= data.qmax) & ~np.isnan(data.y) 294 self.Iq = data.y[self.index] 295 self.dIq = data.dy[self.index] 296 self._theory = np.zeros_like(data.y) 297 q_vectors = [data.x] 298 else: 299 raise ValueError("Unknown data type") # never gets here 300 301 # Remember function inputs so we can delay loading the function and 302 # so we can save/restore state 303 self._fn_inputs = [v[self.index] for v in q_vectors] 304 self._fn = None 305 306 # define bumps parameters 326 self.kernel = kernel 327 partype = kernel.info['partype'] 328 307 329 pars = [] 308 for p in model.info['parameters']:330 for p in kernel.info['parameters']: 309 331 name, default, limits = p[0], p[2], p[3] 310 332 value = kw.pop(name, default) … … 313 335 for name in partype['pd-2d']: 314 336 for xpart, xdefault, xlimits in [ 315 316 317 318 337 ('_pd', 0, limits), 338 ('_pd_n', 35, (0, 1000)), 339 ('_pd_nsigma', 3, (0, 10)), 340 ('_pd_type', 'gaussian', None), 319 341 ]: 320 342 xname = name + xpart … … 328 350 raise TypeError("unexpected parameters: %s" 329 351 % (", ".join(sorted(kw.keys())))) 352 353 def parameters(self): 354 """ 355 Return a dictionary of parameters 356 """ 357 return dict((k, getattr(self, k)) for k in self._parameter_names) 358 359 class Experiment(object): 360 """ 361 Return a bumps wrapper for a SAS model. 362 363 *data* is the data to be fitted. 364 365 *model* is the SAS model from :func:`core.load_model`. 366 367 *cutoff* is the integration cutoff, which avoids computing the 368 the SAS model where the polydispersity weight is low. 369 370 Model parameters can be initialized with additional keyword 371 arguments, or by assigning to model.parameter_name.value. 372 373 The resulting bumps model can be used directly in a FitProblem call. 374 """ 375 def __init__(self, data, model, cutoff=1e-5): 376 377 # remember inputs so we can inspect from outside 378 self.data = data 379 self.model = model 380 self.cutoff = cutoff 381 if hasattr(data, 'lam'): 382 self.data_type = 'sesans' 383 elif hasattr(data, 'qx_data'): 384 self.data_type = 'Iqxy' 385 else: 386 self.data_type = 'Iq' 387 388 # interpret data 389 partype = model.kernel.info['partype'] 390 if self.data_type == 'sesans': 391 q = sesans.make_q(data.sample.zacceptance, data.Rmax) 392 self.index = slice(None, None) 393 self.Iq = data.y 394 self.dIq = data.dy 395 #self._theory = np.zeros_like(q) 396 q_vectors = [q] 397 elif self.data_type == 'Iqxy': 398 q = np.sqrt(data.qx_data**2 + data.qy_data**2) 399 qmin = getattr(data, 'qmin', 1e-16) 400 qmax = getattr(data, 'qmax', np.inf) 401 self.index = (~data.mask) & (~np.isnan(data.data)) \ 402 & (q >= qmin) & (q <= qmax) 403 self.Iq = data.data[self.index] 404 self.dIq = data.err_data[self.index] 405 self.resolution = Pinhole2D(data=data, index=self.index, 406 nsigma=3.0, accuracy='Low') 407 #self._theory = np.zeros_like(self.Iq) 408 if not partype['orientation'] and not partype['magnetic']: 409 raise ValueError("not 2D without orientation or magnetic parameters") 410 #qx,qy = self.resolution.q_calc 411 #q_vectors = [np.sqrt(qx**2 + qy**2)] 412 else: 413 q_vectors = self.resolution.q_calc 414 elif self.data_type == 'Iq': 415 self.index = (data.x >= data.qmin) & (data.x <= data.qmax) & ~np.isnan(data.y) 416 self.Iq = data.y[self.index] 417 self.dIq = data.dy[self.index] 418 if getattr(data, 'dx', None) is not None: 419 q, dq = data.x[self.index], data.dx[self.index] 420 if (dq>0).any(): 421 self.resolution = Pinhole1D(q, dq) 422 else: 423 self.resolution = Perfect1D(q) 424 elif (getattr(data, 'dxl', None) is not None and 425 getattr(data, 'dxw', None) is not None): 426 q = data.x[self.index] 427 width = data.dxh[self.index] # Note: dx 428 self.resolution = Slit1D(data.x[self.index], 429 width=data.dxh[self.index], 430 height=data.dxw[self.index]) 431 else: 432 self.resolution = Perfect1D(data.x[self.index]) 433 434 #self._theory = np.zeros_like(self.Iq) 435 q_vectors = [self.resolution.q_calc] 436 else: 437 raise ValueError("Unknown data type") # never gets here 438 439 # Remember function inputs so we can delay loading the function and 440 # so we can save/restore state 441 self._fn_inputs = [v for v in q_vectors] 442 self._fn = None 443 330 444 self.update() 331 445 … … 341 455 def parameters(self): 342 456 """ 343 344 """ 345 return dict((k, getattr(self, k)) for k in self._parameter_names)457 Return a dictionary of parameters 458 """ 459 return self.model.parameters() 346 460 347 461 def theory(self): 348 462 if 'theory' not in self._cache: 349 463 if self._fn is None: 350 q_input = self.model. make_input(self._fn_inputs)351 self._fn = self.model (q_input)352 353 fixed_pars = [getattr(self , p).value for p in self._fn.fixed_pars]464 q_input = self.model.kernel.make_input(self._fn_inputs) 465 self._fn = self.model.kernel(q_input) 466 467 fixed_pars = [getattr(self.model, p).value for p in self._fn.fixed_pars] 354 468 pd_pars = [self._get_weights(p) for p in self._fn.pd_pars] 355 469 #print fixed_pars,pd_pars 356 self._theory[self.index]= self._fn(fixed_pars, pd_pars, self.cutoff)470 Iq_calc = self._fn(fixed_pars, pd_pars, self.cutoff) 357 471 #self._theory[:] = self._fn.eval(pars, pd_pars) 358 472 if self.data_type == 'sesans': 359 473 result = sesans.hankel(self.data.x, self.data.lam * 1e-9, 360 474 self.data.sample.thickness / 10, 361 self._fn_inputs[0], self._theory)475 self._fn_inputs[0], Iq_calc) 362 476 self._cache['theory'] = result 363 477 else: 364 self._cache['theory'] = self._theory 478 Iq = self.resolution.apply(Iq_calc) 479 self._cache['theory'] = Iq 365 480 return self._cache['theory'] 366 481 367 482 def residuals(self): 368 483 #if np.any(self.err ==0): print "zeros in err" 369 return (self.theory() [self.index]- self.Iq) / self.dIq484 return (self.theory() - self.Iq) / self.dIq 370 485 371 486 def nllf(self): … … 381 496 Plot the data and residuals. 382 497 """ 383 data, theory = self.data, self.theory()498 data, theory, resid = self.data, self.theory(), self.residuals() 384 499 if self.data_type == 'Iq': 385 _plot_result1D(data, theory, view)500 _plot_result1D(data, theory, resid, view) 386 501 elif self.data_type == 'Iqxy': 387 _plot_result2D(data, theory, view)502 _plot_result2D(data, theory, resid, view) 388 503 elif self.data_type == 'sesans': 389 _plot_sesans(data, theory, view)504 _plot_sesans(data, theory, resid, view) 390 505 else: 391 506 raise ValueError("Unknown data type") 392 507 393 508 def simulate_data(self, noise=None): 394 print "noise", noise 395 if noise is None: 396 noise = self.dIq[self.index] 397 else: 398 noise = 0.01 * noise 399 self.dIq[self.index] = noise 400 y = self.theory() 401 y += y * np.random.randn(*y.shape) * noise 509 theory = self.theory() 510 if noise is not None: 511 self.dIq = theory*noise*0.01 512 dy = self.dIq 513 y = theory + np.random.randn(*dy.shape) * dy 514 self.Iq = y 402 515 if self.data_type == 'Iq': 516 self.data.dy[self.index] = dy 403 517 self.data.y[self.index] = y 404 518 elif self.data_type == 'Iqxy': … … 418 532 from . import weights 419 533 420 relative = self.model. info['partype']['pd-rel']421 limits = self.model. info['limits']534 relative = self.model.kernel.info['partype']['pd-rel'] 535 limits = self.model.kernel.info['limits'] 422 536 disperser, value, npts, width, nsigma = [ 423 getattr(self , par + ext)537 getattr(self.model, par + ext) 424 538 for ext in ('_pd_type', '', '_pd_n', '_pd', '_pd_nsigma')] 425 539 value, weight = weights.get_weights( … … 442 556 data = load_data('DEC07086.DAT') 443 557 set_beam_stop(data, 0.004) 444 plot_data(data , data.data)558 plot_data(data) 445 559 import matplotlib.pyplot as plt; plt.show() 446 560
Note: See TracChangeset
for help on using the changeset viewer.